Strong and lasting impacts of past global warming on baleen whales and their prey

Abstract Global warming is affecting the population dynamics and trophic interactions across a wide range of ecosystems and habitats. Translating these real‐time effects into their long‐term consequences remains a challenge. The rapid and extreme warming period that occurred after the Last Glacial Maximum (LGM) during the Pleistocene–Holocene transition (7–12 thousand years ago) provides an opportunity to gain insights into the long‐term responses of natural populations to periods with global warming. The effects of this post‐LGM warming period have been assessed in many terrestrial taxa, whereas insights into the impacts of rapid global warming on marine taxa remain limited, especially for megafauna. In order to understand how large‐scale climate fluctuations during the post‐LGM affected baleen whales and their prey, we conducted an extensive, large‐scale analysis of the long‐term effects of the post‐LGM warming on abundance and inter‐ocean connectivity in eight baleen whale and seven prey (fish and invertebrates) species across the Southern and the North Atlantic Ocean; two ocean basins that differ in key oceanographic features. The analysis was based upon 7032 mitochondrial DNA sequences as well as genome‐wide DNA sequence variation in 100 individuals. The estimated temporal changes in genetic diversity during the last 30,000 years indicated that most baleen whale populations underwent post‐LGM expansions in both ocean basins. The increase in baleen whale abundance during the Holocene was associated with simultaneous changes in their prey and climate. Highly correlated, synchronized and exponential increases in abundance in both baleen whales and their prey in the Southern Ocean were indicative of a dramatic increase in ocean productivity. In contrast, the demographic fluctuations observed in baleen whales and their prey in the North Atlantic Ocean were subtle, varying across taxa and time. Perhaps most important was the observation that the ocean‐wide expansions and decreases in abundance that were initiated by the post‐LGM global warming, continued for millennia after global temperatures stabilized, reflecting persistent, long‐lasting impacts of global warming on marine fauna.


| INTRODUC TI ON
biological change (Tornqvist & Hijma, 2012). Increases in global temperatures by an average of 15°C during the Pleistocene-Holocene transition led to large-scale physical and environmental changes, which in turn led to extensive species redistributions and changes in abundance in a wide range of terrestrial taxa (Brüniche-Olsen et al., 2021;Hewitt, 2000;Lorenzen et al., 2011;Lyons et al., 2010).
Although the effects of the Pleistocene-Holocene transition on terrestrial ecosystems have been well documented, the responses of marine megafauna, such as baleen whales, to these changes remain poorly understood.
With a few exceptions, baleen whales are globally distributed megafauna that feed on invertebrates and fish. Most baleen whale populations undertake extensive, seasonal migrations between low-latitude winter breeding grounds and high-latitude summer feeding areas (Lockyer & Brown, 1981). The trophic position and ocean-wide range of baleen whales suggests that they are subject to environmental and ecological changes across entire ocean basins.
Consequently, it is reasonable to hypothesize that the ocean-wide ecological changes caused by past periods of global warming are mirrored in the long-term demographic history of baleen whales.
Analyses of intra-specific genetic variation can be utilized to gain insights into the responses of species to past climate oscillations and how past climate oscillations structured contemporary ecosystems (for an overview see Colella et al., 2020). Herein, DNA sequences from more than 7000 specimens were employed to infer the demographic histories during the last 30 ky in eight baleen whale species and seven prey species (i.e., fish, krill, and copepods); a period when Earth and its oceans underwent substantial global warming. In order to disentangle the intertwined effects of baleen whale feeding ecology and oceanic context, this study assessed the long-term demographic responses in both baleen whales and their prey across two ocean basins with contrasting oceanographic characteristics (Seibold et al., 2018); the Southern Ocean and the North Atlantic Ocean.
The Southern and the North Atlantic Oceans differ in largescale oceanographic features ( Figure 1). The Southern Ocean is a large ocean basin dominated by the wide and persistent Antarctic Circumpolar Current (Figure 1a; Tynan, 1998) and a comparatively stable pelagic food web that is largely centered around a single species, the Antarctic krill (Euphausia superba; Hopkins, 1985). In contrast, the much smaller North Atlantic Ocean is influenced by multiple, interacting, warm and cold ocean currents as well as continental run-off and cyclic climate oscillations (Figure 1b;O'Hare et al., 2005;Rossby, 1996). The pelagic food web in the North Atlantic Ocean is more complex and dynamic in terms of baleen whale prey compared to the Southern Ocean (Kortsch et al., 2019). In the Southern Ocean, most baleen whale species feed primarily on krill, whereas the same baleen whale species have a more diverse diet in the North Atlantic Ocean (Bluhm & Gradinger, 2008;Kawamura, 1980;Víkingsson et al., 2014). The two oceans also differ in key physical features, such as the extent of seasonal sea ice, as well as the velocity of sea ice reduction during the global warming after the Last Glacial Maximum (LGM), 19-26 kya (Figure 1c-f; Clark et al., 2009;Spindler, 1990).
To the best of our knowledge, no previous study has assessed the correlation in the long-term demographic responses to climate change across multiple predators and their prey on such large spatial scales.

| Taxon selection
This study focused on eight baleen whale species as well as seven fish and invertebrate species, representing baleen whale prey ( Figure 2; Table 1  and the capelin (Mallotus villosus). These species are known baleen whale prey or occupy a similar ecological niche as known baleen whale prey species thus were used as proxies for baleen whale prey (Table 1). Antarctic krill is found in the Southern Ocean, whereas the Northern krill, the copepods Ce. typicus and Ca. helgolandicus, as well as the fishes, are all Northern Hemisphere species. The copepod, P. abdominalis, has a wider distribution, including the Northern and Southern Hemisphere (Figure 2g-l).

| Sample collection
The specimens were collected in the North and South Atlantic Oceans, the southwestern Indian Ocean, and the Southern Ocean.
The latter three areas are collectively referred to as the Southern Ocean (Table 1).
Skin samples from free-ranging baleen whales were collected using remote biopsy sampling techniques (Palsbøll et al., 1991).
Additionally, some samples were obtained during necropsies from beached carcasses, fisheries bycatches as well as local subsistence or commercial whaling operations prior to the moratorium on commercial whaling in 1986. Tissue samples were preserved in 5 m NaCl with 20% dimethyl sulfoxide (Amos & Hoelzel, 1991) and stored at −20 or −80°C. Total-cell DNA was extracted from tissue samples using either standard phenol and chloroform extraction processes (Sambrook et al., 1989) or DNeasy™ columns (Qiagen, Inc.) following the manufacturer's instructions.
Multi-locus microsatellite genotypes (data not included) were employed to identify and remove duplicate samples collected from the same individual and to identify closely related individuals sampled in a non-independent manner (e.g., mother and calf). Closely related individuals sampled at random, that is, during different sightings, were retained in the data set (Waples & Anderson, 2017).

| Mitochondrial DNA sequence data
Mitochondrial DNA (mtDNA) sequence data were obtained from different regions of the mitochondrial genome (Table 1) (CYTB), and 16S ribosomal DNA (16S rDNA). mtDNA sequence data were either generated during this study or obtained from previously published studies (Table 1). The vast majority of published baleen whale mtDNA sequences cover the 5′-end of mitochondrial CR. Therefore, using these sequences facilitates access to large data sets, which comprise the most variable part of the mitochondrial genome. Accordingly, new baleen whale mtDNA sequences generated for this study targeted the same region in the manner described below. The Northern krill data set was generated using a single-strand conformation polymorphism (SSCP) protocol to detect different DNA sequences (i.e., haplotypes).
Only different electromorhps were sequenced implying that some sequence variation may have gone undetected.

| Data processing and sequence alignment
Mitochondrial DNA sequences were aligned using the CLUSTALW algorithm (Thompson et al., 1994) as implemented in MEGA (ver. 6.0, Tamura et al., 2013) with default parameter settings. All DNA sequences of a particular locus and species were trimmed to the same length (Table 1).

| Diversity estimation
Descriptive genetic diversity indices were estimated as implemented in DNASP v.6.12.03 (Rozas et al., 2017) for each species and ocean basin. The indices were; the number of segregating sites (S), number of haplotypes (h), haplotype (H) and nucleotide (π) diversity (Nei, 1987), as well as Tajima's D (Tajima, 1989) and Fu's F statistic (Fu, 1997). Nucleotide sites subject to insertions or deletions and missing data were excluded from the analyses.

| Estimation of changes in genetic diversity and migration rates from single-locus DNA sequences
Temporal changes in regional genetic diversity (Θ) and immigration rates (M) were employed as proxies for inferences of changes in abundance and connectivity, respectively. The genetic diversity in a population is determined by the composite parameter Θ, the product of the mutation rate, and the effective population size (Watterson, 1975), which in turn, is a function of the census population size or abundance. M denotes the probability (scaled with the mutation rate) per generation that an individual is an immigrant. The temporal changes in Θ and M were estimated using the approach implemented in MIGRATE-N (ver. 3.6.6, Beerli & Felsenstein, 1999, 2001. The specific analyses parameter values are tabulated in Table S1. The maximum sample size was set to 250 DNA sequences. For data sets with more than 250 DNA sequences. For larger data sets, 250 DNA sequences were selected at random (without replacement). Random sub-sampling was also employed to reduce sample sizes to the smallest number in order to avoid biases due to confidence interval [CI] of the final estimates was approximated to

assuming a normal distribution of m).
Possible effects of intra-oceanic, population genetic structure were assessed by comparing the outcome of estimates based on pooled and spatially distinct samples. Population samples were analyzed separately when there was a discernible effect of spatial population genetic structure on the final estimate (e.g., Northern krill).
Following Sasaki et al. (2005), the two nominal right whale species were treated as different oceanic populations of a single species, due to their low degree of genetic divergence among oceanic populations of southern and North Atlantic right whales (Rosenbaum et al., 2000), which were similar to inter-oceanic divergences estimated among con-specific populations in other baleen whale species (e.g., Jackson et al., 2014).
Converting the time estimates obtained with MIGRATE-N into calendar years necessitated estimates of generational mutation rates. A range of previously reported mutation rates (in the target or closely related species) was explored; the findings are described in the Supplementary Material "Notes on mutation rates," Tables S2 and S3, and Figures S2 and S3. For the mitochondrial genome, the same generational mutation rate was assumed for all species. However, the final annual mutation rate differs among species because the generation times are species specific. In the case of the CR, a generational mutation rate at 1.125 × 10 −6 per site was employed. This rate was within the range of previously reported estimates ranging from 2 × 10 −7 to 2 × 10 −5 per site per generation (Supplementary Material "Notes on mutation rates," Figure S2). In the case of DNA sequences obtained from other mitochondrial genes (i.e., other than the CR), a generational mutation rate at 3.4 × 10 −7 per site was employed. This value was within the range previously reported for the coding regions in the mitochondrial genome or the entire mitochondrial genome, 2 × 10 −8 -2 × 10 −4 per site (Supplementary Material "Notes on mutation rates," Figure S2).
The consistency among estimates obtained from DNA sequences collected from different regions of the mitochondrial genome was assessed by comparing the estimate of temporal changes in Θ inferred from different mitochondrial genes in the same species (Tables S1 and S4; Figure S3). Conspecific temporal changes in Θ, estimated from mtDNA sequences, were also compared to similar estimates obtained from genome-wide, nuclear data produced by next generation sequencing in three selected baleen whale species (North Atlantic common minke whale, North Atlantic fin whale and southern right whale), as described below (see Section 2.4 below and Table S4; Figure S3).

| Data processing
In the case of the quaddRAD library, PCR clones were removed using the clone_filter script implemented in the software suite STACKS (ver. 1.47, Catchen et al., 2013). Low quality reads were removed from the raw FASTQ files with process_radtags (STACKS, ver. 1.47) using default settings. The output from process_radtags was concatenated in paired-and single-end files for each species. The remaining reads were aligned to a reference genome with BOWTIE2 (ver. 2.2.8, Langmead & Salzberg, 2012) as an "end-to-end" alignment, employing the setting very_sensitive (i.e., D 20, R 3, N 0, L 20, and i S,1,0.50). The maximum fragment size for concordant pairedend alignments was set at 600 and discordant alignments discarded.
Reads from the common minke whale and fin whale were aligned against the common minke whale genome (Yim et al., 2014) and against the bowhead whale genome (Keane et al., 2015) in case of the southern right whale.
The folded site frequency spectrum (SFS, Fisher, 1930;Nielsen et al., 2012) Korneliussen et al., 2014). SNPs with a genotype quality below 20 and a mapping quality below 10 were discarded. Only SNPs genotyped in a minimum of 80% of the individuals were included in the estimation. For each species, two SFSs were estimated, using a minimum read coverage at 10 and 2, respectively. The two data sets differed in terms of the total number of SNPs (Table 1 and  Due to the uncertainties surrounding determining mutation rates, we focused our interpretations on Θ rather than on the effective population size (N e ). The latter is related to Θ and the mutation rate as follows: N e = Θ 4 , where denotes the generational mutation rate.

| Temperature data
Surface air temperature (SAT) estimates for the Southern Hemisphere were inferred from deuterium measurements collected from the Antarctic EPICA Dome C Ice Core and obtained from Jouzel et al. (2007).
For the Northern Hemisphere, continental atmospheric temperatures between 40 and 80°N (calibrated with oxygen isotope records from 57 sediment cores) were obtained from Bintanja et al. (2005).

| Maps of ocean circulation, sea ice reconstructions and species ranges
Maps were generated with ARCGIS® (ver.

| Correlations among baleen whales, prey, and climate
Pearson's correlation coefficients were estimated using R (ver.

| RE SULTS
In total, 4761 mtDNA sequences from eight baleen whale species and 2271 mtDNA sequences from seven prey species (fish and invertebrates) were analyzed to infer temporal changes in Θ and M. In addition, between 14,304 and 62,579 genome-wide SNPs were analyzed in a total of 100 individuals from three baleen whale species in order to assess if the temporal changes in Θ estimated from the mtDNA sequences likely reflected the genome-wide genetic diversity ( Table 1).
The mean genetic diversity was higher among the Southern In contrast, the exponential increases in Θ observed in the fin, common minke, and bowhead whale were delayed until 6-8 kya ( Figure 4d).
Temporal changes in Θ in key prey species, such as krill, copepods, capelin, and herring were estimated as well in order to assess if the post-LGM changes in Θ observed among the baleen whales were correlated with prey availability (Figure 2g-l, Table 1).
In the Southern Ocean, large exponential and synchronous increases in Θ were estimated in Antarctic krill and P. abdominalis, starting around 15-23 kya (Figure 4b). These increases were congruent with the increases observed among the baleen whales in the Southern Ocean. As was the case for the baleen whales, the F I G U R E 3 Haplotype and nucleotide diversity (π) from baleen whales of the Southern and North Atlantic oceans. Numbers represent mean estimates based upon mitochondrial control region DNA sequences. Sample sizes are listed in parentheses with the common species names temporal changes in Θ varied considerably among the prey species in the North Atlantic Ocean. An increase in Θ immediately after the LGM was observed in all invertebrates (i.e., Northern krill, and three copepods, Ca. helgolandicus, Ce. typicus, and P. abdominalis), whereas the estimated increases in Θ in fishes (capelin and herring) were observed 6-8 kya (Figure 4e). In the Southern Ocean, significant and strong correlations (r = .88-1.00, p < .0005) were observed between the krill foraging baleen whales (i.e., the Antarctic minke, southern right, humpback, blue and fin whale) and the two prey species (i.e., Antarctic krill and P. abdominalis, Figure 5a). Ocean (Figure 6b).
The temporal changes in Θ inferred from the genome-wide SNP data in three selected baleen whale species were consistent with and gave support to those inferred from the mtDNA sequence variation (Figure 7), that is, baleen whale population expansions during the Pleistocene-Holocene transition (see Supplementary Material "Notes on genome-wide SNP genotype analyses"). Increasing the read depth of the genome from ×2 to ×10 reduced the number of SNPs by ~50% (Table 1), which in turn muted the degree of the estimated expansions, likely an effect of removing rare variants. The estimated temporal changes in Θ were consistent with the changes estimated from mtDNA sequences in the common minke and the southern right whale at both read depths. While in the fin whale, the median estimate of the temporal changes in Θ did not indicate an exponential expansion at a read depth at ×10, the 95% confidence bands at both read depths were nearly identical (Figure 7; Figure S4).
The estimates of M (i.e., inter-ocean basin connectivity) were subject to substantial uncertainty ( Figure S6). Overall, the interoceanic intra-specific connectivity in baleen whales increased

| Impacts of past global warming
The Pleistocene-Holocene transition after the LGM represented a period when global temperatures rose rapidly, which, in turn, led to drastic reductions in polar sea ice extent and substantial sea level rises (Abrantes, 2000;Clark et al., 2009). The net result of these oceanographic changes was an overall expansion of suitable marine habitat and concurrent, large-scale changes in ocean circulation that resulted in increased primary productivity (Tsandev et al., 2008). These changes were not instantaneous but rather they were set in motion, creating a long-term re-configuration of the marine Ocean is dominated by a single, major and persistent current; the Antarctic Circumpolar Current (Tynan, 1998). In contrast, the North Atlantic Ocean is influenced by several smaller, less stable, interacting cold and warm water currents as well as periodic climate oscillations (O'Hare et al., 2005;Rossby, 1996). The Southern Ocean pelagic food web is also based mostly on Antarctic krill (Hopkins, 1985), whereas the North Atlantic Ocean pelagic food web basis comprises a diverse zooplankton community subject to fluctuations in abundances, in part due to the abovementioned periodic climate oscillations (Pershing et al., 2005). The temporal trend in the rate of increase or decrease observed among North Atlantic baleen whales and their prey species appeared to change substantially 6-10 kya; either transitioning from an increase to a decrease in abundance or from a stable to an exponentially increasing abundance. The so-called "8.2 kya event" (Alley et al., 1997) took place during this epoch. The 8.2 kya event was caused by a massive discharge of glacial melt water that was released into the western North Atlantic Ocean from proglacial lakes, which led to a precipitous drop in global ocean temperatures (Alley et al., 1997;Barber et al., 1999). Sediment core diatom records suggest that the 8.2 kya event resulted in a major shift in phytoplankton composition consistent with a reduction in primary productivity, in particular in the western North Atlantic Ocean (Caissie et al., 2010;Harland et al., 2016). The decline in primary productivity could have led to a reduction of prey and deteriorating environmental conditions, possibly leading to declines in some baleen whale populations. However, apart from the copepod P. abdominalis, the few North Atlantic prey species included in this study did not appear to have declined in abundance after the 8.2 kya event.
Samples from the herring originated mostly from the Baltic Sea. Due to potential founder effects, the demographic history of Baltic Sea populations might differ from those of the North Atlantic.
Additional sampling across space and additional species at lower trophic levels may yield further insights. The level of inter-specific competition among the baleen whales or competition with other marine mammals may also have been affected by the concomitant environmental changes due to the 8.2 kya event in turn leading to changes in abundance.
The temporal differences in the rate of change in abundance observed between most Southern Ocean baleen whale and the common minke whale most likely stem from ecological differences.
Common minke whales are distributed at comparatively lower latitudes and feed on myctophid fishes (Kato & Fujise, 2000). In contrast, the other Southern Ocean baleen whale species in this study, occupy higher latitudes, closer to the sea ice where they feed primarily on krill (Kawamura, 1980).
The strong correlation observed in the changes in abundance between baleen whale species and their prey, for example, between the krill-eating baleen whales and Antarctic krill, suggest that baleen whales responded in synchrony to ocean-wide ecological changes, such as prey availability (Ims & Andreassen, 2000;Korpimӓki et al., 2005;Seyboth et al., 2016), geological conditions (Hansen et al., 2013), or a combination of the two. Similar post-LGM increases have been observed in krill-eating Antarctic penguin species (Frugone et al., 2018;Trucchi et al., 2014;Younger, Clucas, 2015;Younger, Emmerson, 2015), suggesting that the apparent, large increase in primary productivity was sufficient to drive expansions in abundances across all krill predators.

F I G U R E 7
Temporal changes in Θ estimated from mtDNA control region sequences and genome-wide SNP genotypes of two levels of minimum read depth. Estimated demographic history employing mtDNA control region and genome-wide SNP genotypes (at a minimum coverage at ×2 and ×10, respectively) in North Atlantic common minke whale (a), North Atlantic fin whale (b) and southern right whale (c). Time in thousands of years ago (kya) is along the horizontal axis, the estimates of Θ are along the vertical axis. Red and blue shading denotes the Holocene and Pleistocene, respectively. The darkest blue indicates the LGM Numerous studies, either at local or regional scales, have proposed that retreating polar sea ice extent and upwelling led to an increase in the abundance of key zooplankton species, such as krill, copepods, and other meso-zooplankton species (Loeb et al., 1997;Moore, 2016;Stenseth et al., 2002). However, to the best of our knowledge, this study is the first to LGM increases in abundance at the lower trophic levels (e.g., krill and copepods) appeared to precede the increases in abundance at the higher trophic levels (e.g., fishes and baleen whales), although the precision of the exact temporal placement within and among species makes it difficult to draw affirmative conclusions. This initial increase in prey abundance was observed in both ocean basins, suggesting a bottom-up (White, 1978) enrichment of the oceans during the initial warming phase during the Pleistocene-Holocene transition ( Figure 8). This is consistent with previous paleo-oceanographic modeling (Radi & de Vernal, 2008;Tsandev et al., 2008) and a shift in phytoplankton composition from perennial pelagic to seasonal sea-ice-associated species during the Pleistocene-Holocene transition. The latter species are indicative of high levels of primary productivity (Caissie et al., 2010;Harland et al., 2016).

| Critical evaluation of methods
We assessed the impact of employing different mitochondrial genes and compared the results from mitochondrial genes with genomewide SNPs in a subset of species. The results from the different genes ( Figure S3), or between the nuclear and mitochondrial genome (Figure 7), were congruent in terms of the long-term demographic trends. This outcome suggests that the inferences drawn from single-locus mitochondrial genes were robust, which is unsurprising given the long time and large spatial scale of the estimations.
A relevant, but unavoidable caveat to this study, as well as similar assessments based on genetic inference methods, is the inherent uncertainty of the applied mutation rates, which is consequential in terms of inferring the timing of changes in Θ and M. Although the choice of mutation rates has implications for placing the detected changes at a specific annual time, the demographic trends are unaffected to the choice of mutation rate. In this study, the chosen generational mutation rates were within the range of published generational mutation rate estimates in vertebrate and invertebrate species, such as baleen whales, fish, and crustaceans (Supplementary Material "Notes on mutation rates," Tables S2 and S3). However, the possibility that the annual estimates are biased cannot be excluded in this or any similar genetic assessments. Nevertheless, the overall consistency observed in this study among the annual estimates based upon different mitochondrial regions with those estimates obtained from genome-wide SNPs (employing different mutation rates) was reassuring.
Many and different processes affect demographic change, such as fluctuations in habitat and resource availability (e.g., Foote et al., 2013), prey preference (e.g., Fleming et al., 2016;Víkingsson et al., 2014), as well as inter-and intra-specific competition (e.g., Moore, 2016;Moore et al., 2019). Disentangling the relative contributions of all the processes that may contribute to the observed temporal and spatial trends in demographic changes is non-trivial, and will require substantial efforts. The strength of this study lies with the overall sampling design, which includes extensive horizontal (i.e., same trophic level) and vertical (i.e., prey and predators) ecological sampling. The inclusion of multiple species can be viewed as "pseudo-replication": in effect each species, at a given trophic level, represents a single response to the underlying, global processes, such as global warming. The overall consistency in the observed demographic trends across multiple taxa and trophic levels suggests that our analysis captured fundamental drivers of change.

| Implication for future global warming
The rapid rise in global temperatures during the Pleistocene-Holocene transition plateaued ~10 kya, yet most vertebrate and invertebrate taxa included in this study continued to increase in abundance in both ocean basins (until ~1 kya, the most recent time point included in this analysis). In other words, the Pleistocene-Holocene transition appeared set into motion long-term oceanographic and ecological transitions that continued to change abundance and connectivity among baleen whales and their prey during an additional ~10 ky. This observation raises the possibility that the current global warming has already set processes in motion that will result in longterm and wide-ranging rearrangements in marine ecosystems that will continue for millennia after temperatures stabilize. The stabilization of global temperatures during the Holocene may also be a subsequent contributing factor, although exactly how constancy per se would facilitate change (here in abundance and migration rates) is not readily evident. Global warming impacts a wide variety of known and unknown processes; abiotic (e.g., the 8.2 kya event) and biotic, which in turn possibly interact as well. The analysis presented here provides a broad picture, but is insufficient to discern among the many more detailed plausible interactions and processes.
A number of recent studies have attempted to predict the effects of current global warming on marine mammals from contemporary, short-term field observations, and known aspects of each species' ecology (Laidre et al., 2008;Moore & Huntington, 2008;Tulloch et al., 2019). Some baleen whale populations, such as humpback, fin, and blue whales, appear to arrive earlier on the summer foraging grounds (Ramp et al., 2015) and at increasingly higher latitudes, potentially increasing competition with obligate polar species, such as the bowhead whale (Moore, 2016;Moore et al., 2019). Similarly, changes in the distribution and migratory routes linked to sea ice concentrations (Druckenmiller et al., 2018), plasticity in diet preferences associated with changes at lower tropic levels (Fortune et al., 2020), and behavioral adaptations due to increasing killer whale predation (Matthews et al., 2020) have been observed in bowhead whales. Although the findings reported here suggested that most baleen whale species may benefit from the global warming, the results also suggest that the effects of increasing temperatures on baleen whale abundance and migration are complex, with potential wide ranging and long-lasting impacts. In addition, the rate of current ongoing global warming is higher than the post-LGM warming and will likely reaching higher temperatures (Bova et al., 2021;Foster et al., 2017). Regional oceanographic conditions, including temperature, F I G U R E 8 Bottom-up control model of the demographic responses of baleen whales during the Pleistocene-Holocene transition. Red and blue shading denotes the Holocene and Pleistocene, respectively. The darkest blue shading defines the LGM annual sea ice dynamics, and prey base, appear to have been the main driver of the long-term responses seen in baleen whales during past epochs characterized by global warming. However, the rapid and extreme sea ice loss, predicted in the coming decades due to current global warming, including unprecedented ice-free summers (which did not apply to the Pleistocene-Holocene transition) leaves considerable uncertainty as to the future of all whales.

ACK N OWLED G M ENTS
The baleen whale illustrations were reproduced with permission from Frédérique Lucas. Prey illustrations were drawn by Ligia Arreola. Hans J. Skaug and Conor Ryan kindly provided sample materials. Steve Beissinger, Robin Waples, Richard Hudson, and two anonymous reviewers provided helpful suggestions.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.