Receding ice drove parallel expansions in Southern Ocean penguins

Significance We analyze population genomic datasets across 3 penguin genera to test for demographic shifts driven by historical climate events. Numerous species inhabiting coastlines affected by heavy sea ice during the Last Glacial Maximum show genomic signatures of near-simultaneous population expansions associated with postglacial warming, contrasting with stable or declining demographic histories inferred for species occupying consistently ice free habitats. Shallow population genomic structure detected within species distributed across the vast Southern Ocean likely provides further evidence for recent demographic shifts and recent genetic exchange among populations. Our results demonstrate dramatic, ecosystem-wide responses to climate change and highlight the potential for future biological shifts in the Southern Ocean as global warming continues. Climate shifts are key drivers of ecosystem change. Despite the critical importance of Antarctica and the Southern Ocean for global climate, the extent of climate-driven ecological change in this region remains controversial. In particular, the biological effects of changing sea ice conditions are poorly understood. We hypothesize that rapid postglacial reductions in sea ice drove biological shifts across multiple widespread Southern Ocean species. We test for demographic shifts driven by climate events over recent millennia by analyzing population genomic datasets spanning 3 penguin genera (Eudyptes, Pygoscelis, and Aptenodytes). Demographic analyses for multiple species (macaroni/royal, eastern rockhopper, Adélie, gentoo, king, and emperor) currently inhabiting southern coastlines affected by heavy sea ice conditions during the Last Glacial Maximum (LGM) yielded genetic signatures of near-simultaneous population expansions associated with postglacial warming. Populations of the ice-adapted emperor penguin are inferred to have expanded slightly earlier than those of species requiring ice-free terrain. These concerted high-latitude expansion events contrast with relatively stable or declining demographic histories inferred for 4 penguin species (northern rockhopper, western rockhopper, Fiordland crested, and Snares crested) that apparently persisted throughout the LGM in ice-free habitats. Limited genetic structure detected in all ice-affected species across the vast Southern Ocean may reflect both rapid postglacial colonization of subantarctic and Antarctic shores, in addition to recent genetic exchange among populations. Together, these analyses highlight dramatic, ecosystem-wide responses to past Southern Ocean climate change and suggest potential for further shifts as warming continues.

Climate shifts are key drivers of ecosystem change. Despite the critical importance of Antarctica and the Southern Ocean for global climate, the extent of climate-driven ecological change in this region remains controversial. In particular, the biological effects of changing sea ice conditions are poorly understood. We hypothesize that rapid postglacial reductions in sea ice drove biological shifts across multiple widespread Southern Ocean species. We test for demographic shifts driven by climate events over recent millennia by analyzing population genomic datasets spanning 3 penguin genera (Eudyptes, Pygoscelis, and Aptenodytes). Demographic analyses for multiple species (macaroni/royal, eastern rockhopper, Adélie, gentoo, king, and emperor) currently inhabiting southern coastlines affected by heavy sea ice conditions during the Last Glacial Maximum (LGM) yielded genetic signatures of near-simultaneous population expansions associated with postglacial warming. Populations of the ice-adapted emperor penguin are inferred to have expanded slightly earlier than those of species requiring ice-free terrain. These concerted high-latitude expansion events contrast with relatively stable or declining demographic histories inferred for 4 penguin species (northern rockhopper, western rockhopper, Fiordland crested, and Snares crested) that apparently persisted throughout the LGM in ice-free habitats. Limited genetic structure detected in all ice-affected species across the vast Southern Ocean may reflect both rapid postglacial colonization of subantarctic and Antarctic shores, in addition to recent genetic exchange among populations. Together, these analyses highlight dramatic, ecosystem-wide responses to past Southern Ocean climate change and suggest potential for further shifts as warming continues.
Sphenisciformes | climate change | Last Glacial Maximum | refugia | genomics C limate change is substantially impacting the abundance and distribution of wildlife, with many species' ranges shifting poleward as a result of climate warming (1). Similar shifts occurred after the Last Glacial Maximum (LGM; 18,000 to 25,000 y ago) (2,3), as temperate refugial populations of many species expanded into high latitudes. While such range shifts may be readily achieved on continents [where terrestrial habitats are essentially continuous (4)], the challenges are more pronounced for isolated or fragmented populations that rely on long-distance dispersal (5,6). For instance, many high-latitude coastal and terrestrial ecosystems of Significance We analyze population genomic datasets across 3 penguin genera to test for demographic shifts driven by historical climate events. Numerous species inhabiting coastlines affected by heavy sea ice during the Last Glacial Maximum show genomic signatures of near-simultaneous population expansions associated with postglacial warming, contrasting with stable or declining demographic histories inferred for species occupying consistently ice free habitats. Shallow population genomic structure detected within species distributed across the vast Southern Ocean likely provides further evidence for recent demographic shifts and recent genetic exchange among populations. Our results demonstrate dramatic, ecosystem-wide responses to climate change and highlight the potential for future biological shifts in the Southern Ocean as global warming continues.  the Southern Hemisphere are isolated by vast ocean gaps (Fig. 1). Southern Ocean circumpolar fronts (including the Subtropical Front and the Antarctic Polar Front) may present additional physical and thermal barriers to southward range expansion of isolated southern coastal populations (7,8).
Understanding past shifts in species distributions is crucial for forecasting responses to contemporary and future climate change. Currently, there is considerable uncertainty surrounding the extent to which high-latitude wildlife populations might have persisted in the Southern Ocean throughout the LGM versus the extent of post-LGM expansion (6,9,10). Recent genetic data, however, hint at major ecosystem-wide change following reductions in southern winter sea ice (9,11,12). Importantly, past expansions can be reconstructed via genetic analysis of modern populations (2,13). While several studies of Southern Ocean species have detected limited population genetic structure, consistent with recent demographic shifts and/or gene flow (11,12,(14)(15)(16)(17)(18), a comprehensive genome-wide assessment of Southern Ocean wildlife is lacking. Moreover, as responses to climate change can potentially vary among species (12,19,20), distinguishing between concerted (multispecies) versus idiosyncratic (single species) shifts may be crucial to forecasting responses to future climate change (21).
Penguins (Sphenisciformes) are iconic marine birds that inhabit all major southern landmasses, with their greatest species diversity in Antarctica and the subantarctic ( Fig. 1 and SI Appendix, Fig. S1). Although most penguins are natally philopatric (22), some can disperse vast distances traversing major Southern Ocean fronts (23,24) and represent important components of both coastal and marine ecosystems (25). Here we analyze several thousand single nucleotide polymorphisms (SNPs) across 11 Antarctic, subantarctic, and temperate penguin species to test for concerted responses to climate change. We detect genomic signatures of population expansion in multiple species currently distributed largely within the LGM sea ice zone, consistent with concerted recolonization of Antarctic and subantarctic coasts during post-LGM warming. In contrast, demographic histories inferred for 4 temperate penguin species are relatively stable or declining. Our results suggest consistent population dynamics across a species-rich high-latitude assemblage in response to postglacial ice reduction and demonstrate the potential for rapid change to Southern Ocean ecosystems under future warming.

Results
Demographic reconstructions of effective population sizes (N e ) for 11 penguin species using CubSFS (26), SNAPP (SNP and AFLP package for phylogenetic analysis) (27), Tajima's D (28), and Multidice (29) were based on 3,000 to 13,000 SNPs per species (SI Appendix, Tables S1-S3). Macaroni and royal (Eudyptes chrysolophus chrysolophus/E. c. schlegeli) penguins were considered a single species based on structure/F ST (fixation index) analyses (see also ref. 18), whereas Snares crested (E. robustus) and the northern rockhopper (E. moseleyi) penguins were excluded from some analyses due to their small sample sizes ( Fig. 1 and SI Appendix, Tables S4-S5). These analyses revealed comparable postglacial N e expansions for 6 3A). Notably, these 7 species all predominantly occur south of the LGM sea ice limit ( The expansion time frames inferred for most southern lineages (20,000 to 15,000 y ago) correspond to a period of rapid post-LGM warming (31) ( Fig. 2A and SI Appendix, Table S1). These reconstructions suggest populations of the ice-adapted emperor penguin expanded earlier than those of most other southern penguin lineages which require ice-free terrain (see also refs. 15 and 32). The magnitude of inferred postglacial N e expansions is, on average, a 2.7-fold increase (ranging from 1.19-to 4.4-fold increase) ( Fig. 2A and SI Appendix, Table S1). We detected some variation in the outcomes of different demographic analyses for particular species, perhaps a reflection of varying sensitivity of different model-based approaches and/or biological signal. For example, the CubSFS analysis contrasted with other approaches in suggesting chinstrap penguin populations expanded prior to the LGM and declined following the LGM. Overall, however, there is broad support for "stable/declining" demographic trajectories for species inhabiting LGM ice-free regions versus predominantly "expanding" trajectories for LGM ice-affected species (Fig. 3A).
We used Multidice to test for synchronous versus asynchronous expansions across the 7 expanded species identified based on our demographic analyses (Fig. 3A). To this end, we modeled a single expansion event within the last 50,000 y in which up to 7 species coexpanded. The synchronous expansion event was inferred to have occurred 20,779 to 24,804 y ago, depending on the summary statistics chosen (Fig. 3B and SI Appendix, Table S3). While only 2 or 3 of these southern species were inferred to have expanded simultaneously (SI Appendix, Table S3), minor differences between inferred expansion timings ( Fig. 2A) likely hindered the ability for Multidice to detect a single expansion event corresponding to all expanding species.
Tests for intraspecific genomic divergence across the ranges of individual species (including previous analyses of Pygoscelis and Aptenodytes species; see refs. 11, 12, 33, and 34) consistently revealed shallow genetic structure within species ( Fig. 1 and SI  Appendix, Figs. S4-S7 and Tables S6 and S7). In all cases apart from gentoo penguins, we found that panmixia (K [number of genetic clusters] = 1) was supported but that using location priors found evidence for additional fine-scale structure, as previously reported (11,12,14,33,34). Such patterns are consistent with post-LGM demographic and biogeographic expansions (for southern LGM sea ice species) and recent genetic exchange among populations. Specifically, F ST , principal coordinates analyses (PCoA), Structure, discriminant analysis of principal components (DAPC), SNAPP, and phylogenetic analyses for Eudyptes, Pygoscelis, and Aptenodytes all revealed relatively shallow within-species genomic structure among southern populations ( Fig. 1 and SI Appendix, Figs. S4-S7 and Tables S6 and S7) (14,15,34,35). In contrast to the recent genetic exchange inferred within most species and between macaroni and royal penguins (17,18,35), these analyses detected little or no admixture among species (SI Appendix, Figs. S4-S7).

Discussion
Our study detected broadly consistent genome-wide signatures of post-LGM expansion across penguin species that currently breed south of the LGM sea ice zone (Fig. 3A). By contrast, 4 species currently breeding north of the LGM sea ice zone exhibited genetic signatures of relatively stable or declining demographies (Fig. 3A). Although estimates of precise LGM breeding ranges for penguins remain elusive (but see ref. 36), our findings are consistent with the hypothesis of (6) that during the LGM, many Southern Ocean species retreated to ice-free refugia (e.g., Gough, Amsterdam, and Falklands islands; southern South America; and New Zealand's southern islands) ( Fig. 2B; see refs. 6 and 9). Indeed, several recent studies have suggested that post-LGM reductions in sea ice were accompanied by rapid recolonization of high-latitude shores (8,9,12) (Fig. 2C). Recent demographic studies of penguins (Adélie, emperor, and king) (15,16,32) and the southern elephant seal (37), for example, inferred rapid postglacial recolonization events. By contrast, recent snow petrel analyses provide only limited evidence for such postglacial shifts (10). The choice of mutation rate, and possibly time dependency issues, might play some part in these apparently conflicting patterns among taxa. Some contrasting responses among species may also stem from interspecific ecological differences (e.g., variation in feeding ecology, philopatry, habitat preferences). Shifting oceanographic and coastal environmental features associated with postglacial warming may also have impacted local species.
While most LGM coasts are now inundated (Fig. 2B), some potential LGM refugia may be suggested on the basis of current distributions (e.g., the eastern rockhopper penguin likely expanded south from the Auckland, Campbell, and Antipodes islands; Fig.  2C). Previous studies have concluded that the Southern Ocean's circumpolar fronts can represent important barriers to dispersal for many marine species (7,8), including penguins (14,38). However, several penguin species can clearly traverse such boundaries (23,24), and this exceptional dispersal ability may help to explain their apparently rapid biogeographic shifts in response to changing climate (see also ref. 37).
While CubSFS suggested the chinstrap penguin may have declined following the LGM, Tajima's D and SNAPP supported population expansion for this species, comparable to results for other southern species (Fig. 3A). This anomaly may perhaps reflect issues with the mutation rate and/or generation time used or may indicate an idiosyncratic ecological response for this southern species (e.g., variation in feeding ecology, philopatry, habitat preferences, sensitivity to oceanographic fronts). Based on evidence from combined demographic analyses (Fig. 3A), the suggestion that chinstrap penguins have declined since the LGM should be treated with some caution.
A consistent finding of our study is the lack of major genomewide differentiation across the ranges of most penguin species, including several species showing circumpolar near homogeneity LGM breeding ranges in both B and C are uncertain. Maps adapted from ref. 6. Copyright (2012) with permission from Elsevier. As the Snares crested penguin was included in other demographic analyses (Fig. 3), the species is shown in both B and C. Colored symbols (squares and circles) are consistent with Figs. 1 and 3. (15, 16) (i.e., K = 1; F ST < 0.02; Fig. 1 and SI Appendix, Table S6). These relatively shallow F ST values contrast with more substantial structure and evidence for multiple Southern Ocean refugia in white-chinned petrels [K = 3; F ST > 0. 10 (39)]. While biallelic markers such as the SNPs analyzed here are theoretically capable of yielding F ST as high as 1 (i.e., fixed differences at all loci), we note that the upper range of this parameter can be limited by allele frequency distribution (40), and thus, these values should be treated with some caution. While the use of location priors at higher values of K reveals additional, fine-scale population differentiation ( Fig. 1 and SI Appendix, Fig. S4) (see also refs. 14, 33, and 34), such structure can potentially evolve rapidly (e.g., ref. 41). Interestingly, the relatively shallow differentiation observed within and among some colonies [e.g., emperor (14,34)] may also provide additional evidence of recent or ongoing gene flow and admixture, sometimes over vast distances (Fig. 1). Subtle population differentiation detectable with location priors might reflect the influence of contemporary oceanographic fronts and/or changes in local sea ice conditions, as previously suggested by (11,14,16,17), and may have considerable relevance over ecological time frames (e.g., conservation management, studies of migration).
Understanding how biota responded to past climate change is essential for predicting species distributions and population sizes under future climate projections and for developing appropriate conservation management strategies (11,42). As global temperatures continue to increase, midlatitude biota will continue to shift toward the poles (8) or, alternatively, may face extinction (6,8). Many penguin populations are currently declining or are predicted to decline as warming continues (43)(44)(45). Some of the northernmost colonies of Adélie and emperor penguins have already disappeared (43,46), and in the case of emperor penguins, these changes have been linked directly to reductions in sea ice (47). By contrast, populations of gentoo penguin are apparently expanding their ranges southward as the climate warms (48). Our study broadly demonstrates the demographic sensitivity of Southern Ocean wildlife to the effects of past climate change (49), highlighting the potential for future shifts under anthropogenic climate change. DArT-Seq Library Preparation and Filtering. DNA was extracted from 428 Eudyptes penguin samples spanning 6 species (Fig. 1 and SI Appendix, Fig. S1 and Table S4; macaroni/royal penguins were considered a single species; see refs. 17, 18 and 35) using a modified Qiagen DNeasy Blood and Tissue kit. Library preparation and SNP discovery were performed on the 282 highestquality DNA extracts using Diversity Arrays Technology Pty. Ltd. (DArT-seq) in Canberra, Australia (50). Each sample was processed following (51) and was sequenced across 3 lanes on an Illumina Hi-Seq 2500. Sequences were processed using in-house proprietary DArT analytical pipelines. We used DartR version 1.1.6 (52) in R version 3.5.1 (R Core Team, 2018) (53) to filter the DArT-seq data for 10 separate Eudyptes datasets [based on previous systematic discussions (17,35) and SI Appendix, Table S5]. For these Eudyptes datasets, we filtered on reproducibility (t = 1) and filtered out monomorphic loci, loci with call rates <0.95%, all individuals with call rates <0.90%, all loci with trimmed sequence tags, and all loci that departed from Hardy-Weinberg equilibrium in any colony (P = 0.05 following Bonferroni correction). We also obtained filtered RAD-seq datasets from an additional 5 penguin species generated and examined by refs. 14, 33, and 34, comprising Adélie (n = 87), gentoo (n = 36), chinstrap (n = 44), king (n = 64), and emperor (n = 110) penguins (SI Appendix, Table S8). See SI Appendix for details.

Materials and Methods
Phylogenomic Analysis and Population Structure. To clarify the evolutionary relationships among our Eudyptes samples sequenced in this study, we created a maximum likelihood phylogeny using RAxML-HPC version 8.2.1 (54) (SI Appendix, Fig. S5). We undertook similar population structure analyses for Eudyptes as previously implemented for Pygoscelis and Aptenodytes in ref. 14 as follows: We calculated population summary statistics, including the number of private alleles, observed and expected heterozygosity, the inbreeding coefficient, and global and pairwise F ST (Fig. 1 and SI Appendix, Tables S5-S7). Genetic clusters were visualized using 3 methods: PCoA using adegenet (55) (SI Appendix, Fig. S4), the Evanno method (56) in Structure version 2.3.4 (57) to estimate the most likely K ( Fig. 1 and SI Appendix, Fig. S4), and DAPC using adegenet (SI Appendix, Fig. S4). We used the SNAPP tree set analyzer in BEAST version 2.4.7 (27,58) to investigate gene flow between closely related Eudyptes species (SI Appendix, Fig. S6) based on our results and systematic discussions of refs. 17 and 35. While SNAPP analyses were previously generated for the emperor, king, and gentoo datasets (14,33,34), we also undertook SNAPP analyses for the chinstrap and Adélie penguin datasets obtained from (14) (SI Appendix, Fig.  S7). See SI Appendix for details.
Testing for Demographic Expansions. We reconstructed population histories for 11 Eudyptes, Pygoscelis, and Aptenodytes species over the last 1,000,000 y by estimating the time and magnitude of demographic changes using 4 different approaches (northern rockhopper and Snares crested penguins were excluded from some analyses due to low sample size). Specifically, we reconstructed the demographic histories using CubSFS ( Fig. 2A and SI Appendix, Figs. S2 and S3 and Table S1), obtained Tajima's D (SI Appendix, Table S2), identified the change in theta values as inferred by our previous SNAPP analyses (SI Appendix, Figs. S8 and S9), and tested for synchronous expansion using Multidice (Fig. 3B and SI Appendix, Table S3). As the Eudyptes and Pygoscelis/Aptenodytes datasets were obtained using different pipelines (DArT-seq versus RAD-seq), we applied further stringent filtering to ensure consistency between the datasets (SI Appendix). While previous studies reported shallow population genetic structure within most Pygoscelis and Aptenodytes species (11, 14-16, 33, 34) (Fig. 1), given the relatively shallow F ST values involved and reported K = 1 (14,33,34), we consider this fine-scale structure likely to have evolved recently and broadly consistent with a scenario of high gene flow suitable for combined demographic analysis [with the exception of the gentoo penguin (14)]. To account for the deeper genetic structure observed in gentoo populations [e.g., 4 distinct lineages (14)] (Fig. 1), we limited most subsequent analyses of this species to one lineage ( Fig. 1)   the folded site frequency spectrum down to increase the number of segregating sites using EasySFS. We then adjusted the number of monomorphic sites in our site frequency spectrum to reflect the total number of monomorphic loci within each species following down projection. For all analyses, we assumed a generation time of 8 to 14 y (from ref. 42). For all models, we used a mutation rate of 2.6 × 10 −7 per locus per generation (12,15).