Bayesian Skyline Plots disagree with range size changes based on Species Distribution Models for Holarctic birds

During the Quaternary, large climate oscillations impacted the distribution and demography of species globally. Two approaches have played a major role in reconstructing changes through time: Bayesian Skyline Plots (BSPs), which reconstruct population fluctuations based on genetic data, and Species Distribution Models (SDMs), which allow us to back‐cast the range occupied by a species based on its climatic preferences. In this paper, we contrast these two approaches by applying them to a large data set of 102 Holarctic bird species, for which both mitochondrial DNA sequences and distribution maps are available, to reconstruct their dynamics since the Last Glacial Maximum (LGM). Most species experienced an increase in effective population size (Ne, as estimated by BSPs) as well as an increase in geographical range (as reconstructed by SDMs) since the LGM; however, we found no correlation between the magnitude of changes in Ne and range size. The only clear signal we could detect was a later and greater increase in Ne for wetland birds compared to species that live in other habitats, a probable consequence of a delayed and more extensive increase in the extent of this habitat type after the LGM. The lack of correlation between SDM and BSP reconstructions could not be reconciled even when range shifts were considered. We suggest that this pattern might be linked to changes in population densities, which can be independent of range changes, and caution that interpreting either SDMs or BSPs independently is problematic and potentially misleading.

In the absence of direct evidence of species' prehistoric response to alterations in climatic conditions, sophisticated statistical methods have been developed to infer demographic history from more indirect sources. Different types of data, and consequently different approaches, are generally favoured by different fields. One of the most widely used genetic methods for inferring demographic history is the "skyline plot," a family of graphical, nonparametric methods introduced by Pybus et al. (2000). Grounded in the principles of Kingman's coalescent theory (Kingman, 1982), the "skylineframework" uses DNA sequences to reconstruct a gene tree. The rate of coalescent events within it can be used to infer how the population size changed over time: periods of low coalescent rates imply a large population while a high density of coalescences implies a small population. Although skyline plots have been used to reconstruct demographic histories for extant and extinct species (Stiller et al., 2010), and across taxa that include vertebrates (Lu et al., 2012;Vignaud et al., 2014), invertebrates (Sanchez et al., 2016;Villalta et al., 2018) and even bacteria (Segawa et al., 2018), comparative multispecies studies are only now emerging (Burbrink et al., 2016).
Skyline plots have been used extensively to infer the response of species during the Last Glacial Maximum (LGM), and are often paired with climatic reconstructions to infer the changes in available habitat for a given species (Calderón et al., 2016;Foote et al., 2013;Lorenzen et al., 2012). One popular approach for reconstructing possible changes in available habitat for a species through time is bioclimatic Species Distribution Models (SDMs) (Elith & Leathwick, 2009). Modelling algorithms combine data on occurrences with environmental information to describe a species' current distribution and estimate how changes in environmental variables may have influenced this range through time. The underlying logic in linking these approaches is that, assuming limited population structure and appropriate sampling, a skyline plot could provide an indication of changes in total population size, which, in turn, may affect the extent of the range occupied by a species. However, the association between effective population sizes (N e ), as reconstructed by skyline plots, and species ranges is generally assumed, not tested.
There are many reasons why N e might not be a good proxy for species ranges. Population structure is known to be a confounding effect, and the recommendation to counter it is to pool samples from multiple locations (Heller et al., 2013). However, even with this sampling scheme, there might be a mismatch in the two quantities if mean population density, and thus N e , was affected by climate change differently to total range extent. For example, an increase in mean population density, and thus population size, might occur without a change in range, if the quality of habitat and its carrying capacity increased without a change in its extent (Figure 1a; Figures S1, S2) (Fordham et al., 2012). Given the positive relationship generally observed between range extent and mean local population density (Connor et al., 2000), N e would also be expected to increase by a greater proportion than range extent under climatic amelioration.
Another possibility is that, without substantial gene flow, skyline plots may reconstruct the population dynamics of the sampled locations rather than the whole species (Miller et al., 2018).
Another scenario that might lead to a disconnect between local N e and range size arises during range shifts: sampled locations, suitable for a species today, might have been only marginally suitable in the past; that is, what is now the core area occupied by a species (where it is abundant, and sampling is more likely) might be inhabited by populations that in the past were at low densities because the area was only marginally suitable. Skylines from such populations would reveal a strong increase in N e reflecting the local amelioration of conditions, irrespective of broader range changes ( Figure 1b; Figures S1, S3). A similarly confounded signal will be found in the more extreme scenario where, as a result of a sizeable range shift, the sampled populations inhabit areas which were completely unsuitable in the past, and thus have undergone a founder event. Such populations would be characterized by an increase in N e as they recovered from the local bottleneck associated with the founder event ( Figure 1c; Figures S1, S4).

F I G U R E 1
The top of each panel represents two schematic maps of the species range at the LGM (right) and present (left). Colour density within each virtual range shows population density. Blue crosses represent the geographically stable genetic sampling location. (a) Increasing population size without range change recovers an increasing BSP profile, as seen in the lower panel. (b) Core area today was only marginally suitable in the past, range size remains the same but local amelioration leads to a local increase in N e . (c) Modern sampled populations inhabit areas outside the past species range. BSP recovers a steep increase in N e associated with a founder event. See Material S1 for simulations supporting these expectations The extent to which changes in population density and bottlenecks related to range shifts can override the signals linked to changes in range size and the overall metapopulation size is unknown. We collected from GenBank a comprehensive data set of mitochondrial DNA (mtDNA) sequences from many species of Holarctic birds, reconstructing their population dynamics using Bayesian Skyline Plots (BSPs). A simple prediction, based on the relative changes of habitat types as reconstructed from the pollen record (Allen et al., 2010), is that species associated with closed habitats (e.g., forests) should have increased since the LGM, whilst species from open and semiclosed habitats (grasslands and steppes) should show a decrease. However, species have more complex niche requirements than a simple association with a broad habitat type, and a more realistic prediction is that N e should change in line with changes in extent of their geographical range. We therefore modelled changes in potential range size between the LGM and the present, using palaeoclimate reconstructions and SDMs, and investigated the relationship between changes in N e and range size. Reconstructing individual ranges in the past is challenging because species may depend upon habitats whose extent is not easily reconstructed (e.g., wetlands); therefore, we also compared demographic profiles of species grouped by major biomes to investigate whether there was any consistent pattern in their response to climatic amelioration in the Holocene.

| Raw genetic data
We collated summary information on all available avian NADH dehydrogenase subunit 2 (ND2) and cytochrome b (cytb) sequences in GenBank, two of the most frequently uploaded avian mitochondrial genes, and screened for Holarctic species using the list of Voous (1977). Only species with more than 10 accessions for either gene were retained. When species had sufficient data for both ND2 and cytb, sequences for each gene were extracted and handled as distinct data sets.

| Alignment
Sequence data for each species/gene combination come from multiple studies, so comparable regions were found by aligning sequences in mega (version 7.0; (Kumar et al., 2016)) using clustalw (Thompson et al., 1994). Sequences for each taxon were trimmed to the longest common section between all samples. If inclusion of a single sequence required the loss of >200 bp from more than 50% of the other sequences, that sample was excluded as were all positions containing insertions, deletions or sequencing ambiguities. When studies uploaded one copy of each haplotype, associated frequency data were used to generate the appropriate number of copies in our database. If frequencies were unpublished the data were excluded.
Sample sizes varied from 11 to 453 sequences per species, with lengths from 236 to 1,137 bp.

| Median joining networks (MJNs)
For each species/gene combination, we built an MJN in popart (Leigh & Bryant, 2015). If the MJN contained long branches (single branch with 30 or more substitutions, indicative of profound population substructure) the species sampling location was reviewed. If clear geographical separation or grouping was found, the data were divided as appropriate and treated as discrete data sets. Single samples with >30 mutations on a branch were considered extreme outliers and dropped from alignments.

| Mutation rate
Recent work (Nabholz et al., 2016) proposes that body mass can be used to accurately calculate taxon-specific per year substitution rates and provides a correction factor for rate variation according to body mass as well as major mtDNA loci. We created data set-specific molecular evolution rates using average body mass data for each of our species (Dunning, 2007)
A gene-specific correction factor (0.052 for ND2, −0.069 for cytb) was then applied (see Table S3). "Calibration set 4" was used as it includes younger species splits that should lead to estimates more appropriate for the within-species dynamics investigated here. Due to the uncertainty surrounding mutation rates, analyses based on "Calibration set 2" were also run (data not shown).

| BSP analysis
There now exist a range of related skyline plot methods (Ho & Shapiro, 2011). Here, we focus on BSP (Drummond et al., 2005) because its relative simplicity and inherent robustness make it popular in the field and particularly suitable for the quality of data sets in our study. For each data set, BSP analyses were implemented in beast2 (Bouckaert et al., 2014) using a strict clock with taxon-specific body mass/gene mutation rates, running 300 million steps sampled every 30,000 steps, with the first 10% discarded as burn-in. "bModelTest" was used to select the most appropriate site model and parameters for individual analyses (Bouckaert, 2015) whilst "bGroupSizes" and Sampling strategy can have important effects on BSP profiles (Chikhi et al., 2010;Heller et al., 2013;Städler et al., 2009).
Geolocation information is only available for a small proportion of samples; for a subset of five species with large sample sizes and good geolocation coverage, we confirmed that using all samples provided representative profiles compared to more structured subsampling (Materials S2, Figure S5A-E). Furthermore, we also confirmed (Materials S3, Figure S6) that, for these species, BSPs recovered profiles in line with those from Extended Bayesian Skyline Plots (EBSPs) (Heled & Drummond, 2008), which are more powerful but also require larger data sets (EBSPs were designed to take advantage of multilocus data).

| Inclusion criteria
Where data were available for both ND2 and cytb, BSP profiles were compared along with summary statistics on each data set. When profiles agreed, the best supported data set was retained (e.g., largest sample size, longest sequences) to represent that species' history. If profiles were discordant but there was a clear disparity in the data sets' quality, we kept the profile from the stronger data set.
If data quality did not vary but profiles were discordant, we conservatively rejected both profiles. For inclusion in further analysis, profiles needed to have a history deeper than 5000 years ago (ka) but shallower than 1 million years ago (deeper histories probably resulted from problematic fits). Also, because of the sample sizes and sequence lengths available for many species, we had to limit our BSP analysis to detect only one major size change (i.e., 3 change points). For this reason, if the largest change occurred before 60 ka, we excluded the species as we would have been blind to any smaller changes that occurred in the period of interest.

| Habitat classification
The species considered are associated with a wide range of habitat types, especially in terms of the predominant vegetation. Given the need for a small number of categories, no classification system can be perfect. We used the expert ornithological opinion of one

| Timing of expansion
Identifying a population's point of expansion from a BSP profile proved difficult given the wide range of shapes present: some population sizes changed little or gradually, while others showed sharp and/or multiple points of inflection. We attempted to develop algorithms that could automatically capture the inflection point but found them to be unreliable, and opted for visual scoring by E.F.M.
Scoring was done prior to the running of SDMs, and thus blind to the expected changes from the range reconstructions.

| Size change
To compare the magnitude of estimated population expansions, the relative population size change between 21 ka and the present day was calculated. N e at 21 ka was interpolated for each BSP profile and, where profiles were shorter than 21 ka, the population was assumed to be the size at the start of profile. Given the uncertainty of dating with molecular methods, and as there is a known risk of false signals in BSP profiles at very recent times (Heller et al., 2013), we also ran our analysis using size estimates from 60 ka/the earli-

| Phylogenetic correction
All analyses included a phylogenetic correction based on the most complete molecular phylogeny of extant birds (www.birdt ree. org, Jetz et al., 2012). Jetz et al. (2012) present two phylogenetic backbones, respectively Ericson and Hackett ("backbone-E" and "backbone-H"), considered equiprobable. For each backbone, we generated 1000 trees, randomly resolving polytomies in each, and repeated all analyses for each tree with a "corPagel" phylogenetic correction from the "ape" package in r. We also provide the results based ordinary least squares (i.e., without phylogenetic correction).

| Species Distribution Models
Our pipeline is provided as an R script in the Supporting Information.

| Climate reconstructions
To identify areas suitable for species through time, we used a 0.5° resolution data set for 19 bioclimatic variables (Net Primary productivity, Leaf Area Index and all the BioClim variables (Hijmans et al., 2005) excluding BIO2 and BIO3) that are available for the last 21,000 years in 1000-year time steps (Beyer et al., 2020). This data set was constructed from a combination of HadCM3 climate simulations of the last 120,000 years, high-resolution HadAM3H simulations of the last 21,000 years and empirical present-day data. The data were downscaled and bias-corrected using the Delta Method (Beyer et al., 2019).

| Species data
Species occurrences were downloaded from the GBIF database (https://www.gbif.org) without any filtering (links in Table S4).
Occurrences were then filtered based on coordinate accuracy (maximum error: 10 km), keeping only observations within the breeding and resident geographical ranges from Birdlife (BirdLife International, 2018). Points were regridded based on the palaeoclimatic reconstructions (0.5° x 0.5°) and, to work on presence/absence not frequency, only one presence per grid cell was kept.
For each species, this cleaned data set was used to select a subset of the bioclimatic variables available in the palaeoclimatic reconstructions (Beyer et al., 2020). We used a standardized approach in which we only removed highly correlated variables (r > .7), as these can lead to model instability (Guisan et al., 2017). A correla- in North America.
To reduce the risk of geographical bias from opportunistic observations affecting the SDMs, the data set was thinned using the r

| Model fitting
SDM was performed using the r package biomod2 (Thuiller et al., 2019) for all species with more than 10 occurrences after filtering and thinning (Stockwell & Peterson, 2002). The thinned data set was used as presences, the landmass of Eurasia or North America as background, and the same number of pseudo-absences as presences were randomly drawn five times from land outside the BirdLife resident and breeding masks, creating five independent data sets for further analysis. We found that drawing pseudo-absences in this way, from outside the masks, was the most effective strategy for retrieving SDMs consistent with species' modern-day ranges. By confirming that the estimated distributions recovered for the modern day were in accord with the BirdLife range predictions, we were confident that the modelled niche being projected into the past was as accurate as possible ( Figure S7).
Following Bagchi et al. (2013), models were run independently for each of the five pseudoabsence data sets using four different algorithms: generalized linear models (GLMs), generalized boosting method (GBM), generalized additive models (GAMs) and random forest. Model evaluation was performed by spatial cross-validation (Roberts et al., 2017), namely splitting the data sets (presences and all five runs of pseudoabsences) based on latitudinal bands in America and longitudinal bands in Eurasia ( Figure S8) with the r package BlockCV (Valavi et al., 2019), and using 4/5 of the splits to calibrate the model and the remaining 1/5 to evaluate it. As a data split containing only absences cannot be used for evaluation, and given the great variety of distribution of the species analysed, we maximized the probability of having at least some presences in all data splits by creating 15 spatial blocks encompassing the whole region of interest, either North America or Europe. Each block was given an ID, numbered sequentially 1-5 repeating, and the 15 blocks were assembled into these five working data splits grouped by the assigned ID numbers.
The models were run five times (once for each pseudoabsence run) for each of the four mentioned algorithms, using in turn four of the five defined data splits to calibrate and one to evaluate based on the TSS (threshold =0.7).
A full ensemble, combining all pseudo-absences sets and algorithms (Araújo & New, 2007), was built using only models with TSS > 0.7 av-

| Range change and overlap
In order to calculate the range of each species, such continuous projections were transformed into binary (either 0-absence more likely, or 1-presence more likely) using 0.5 as a threshold. The binary projections were then used to estimate the climatically suitable area (km 2 ) for each species now and in the LGM. To do so, we reprojected the rasters to the Eckert IV equal-area pseudocylindrical projection setting the grid size to 50 × 50 km and multiplied the number of cells occupied in each period, and their overlap, by 2500 km 2 (cell area).

| Range size comparison between sample species and all Holarctic species
We used BirdLife data on the Extent of Occurrence (km 2 ) to define range size of all Holarctic bird species, those with a mean range latitude of above 20°N. We then compared this full data set to a subset of range sizes for the species included in our study using a Wilcoxon rank sum test.

| Summary of available BSPs
A scan of GenBank yielded a preliminary data set of 208 species having a minimum of 10 individuals sequenced for either the ND2 or cytb genes. Data sets were then discarded because of: insufficient haplotypes for demographic reconstruction (<5), insufficient sequence length (<200 bp), sequences across studies not from comparable gene sections, unpublished haplotype frequencies, inappropriate sampling strategies used by the original study (e.g., nonrandom) or extensive population substructure (full details in Materials and Methods). Application of these criteria left 167 data sets for BSP analyses. They were analysed with beast2 using a BSP and, with one exception (king eider, Somateria spectabilis), converged successfully.
In beast2, we resized the "bGroupSizes" parameter, potentially constraining the level of detail recoverable in the profiles but providing comparable estimates between variable-quality data sets, analysed using the same settings.
Data for both ND2 and cytb were available for 28 species. For 18 species, expansion times and profiles across them were consistent, and a single profile was selected to illustrate the population history.
For 10 species, the two genes showed discordant demographic histories but, in seven cases, one data set was of appreciably lower quality (e.g., fewer samples, shorter sequences) and was removed. Three species were rejected because the two genes gave discordant profiles despite appearing to be of comparable quality (details in Table S1).
Further profiles had to be excluded as they were too deep (limit ≥ 1,000,000 years before the present, n = 9) or too short (limit ≤ 5,000 years before the present, n = 4) to be informative for There was no indication that species associated with particular habitats were more/less likely to be excluded (p = .77, Fisher's exact test, excluding "Other" category as it had too few species for testing). Skyline profiles encompass a wide range of shapes, variously exhibiting single sharp inflection points, gradual changes in size and multiple points of change (see Figure S9). No significant differences were found in the proportion of ND2/cytb genes in each habitat type (p = .78, Fisher's exact test, excluding "Other" category), nor in the proportions of species from the Palearctic, Nearctic or Holarctic in each habitat type (p = .10, Fisher's exact test, excluding "Other" category).

| Direction and magnitude of demographic change
Only

| Direction and magnitude of change in potential geographical range
To investigate the plausibility of changes in the extent of climatically suitable area contributing to the overwhelming majority of profiles showing an increase in N e , we created individual SDMs for each of our species. We had credible SDMs for 96 species (Table S2); five species had to be excluded as there were insufficient observation points left for analysis after data thinning, and one species was rejected as its SDMs led to a present-day projection much larger than the observed range. For the valid SDMs, we generated ranges for the LGM (21 ka) and present day and quantified the changes in their size.
As was the case for N e trajectories across the BSPs, many species

| Timing of change
We next explored the relationship between the timing of the dominant population size change (since 60 ka) and habitat type (excluding associated with the LGM. However, using the rate from "Calibration set 2" (which includes older nodes than set 4) would recover older expansion dates (data not shown); as such dates would correspond to periods of high ice coverage, they seem less likely.
For an alternative view of when the expansions occurred, we used an MSI. The MSI depicts normalized changes in size averaged across species within each habitat type for each time point (Figure 3b). Despite exploiting a different aspect of the BSP profile shape, mean change at each time point rather than a single mean date of maximum change, MSI profiles reveal a pattern that is strongly supportive of the previous result, where wetland species expand appreciably later than species in the other three habitats.
We also confirmed there was no correlation between body mass and timing of population size change (p = .171, Figure S12).

| DISCUSS ION
We gathered a large collection of mtDNA data sets from many bird species to look for habitat-associated trends in population size through time. Although variable data quality may lead to uncertainties about the magnitude of any particular change in population size we detect, the direction of change is relatively robust (Grant, 2015;Grant et al., 2012). Of 102 species, only four species show an overall decrease in effective population size (N e ). Changes during the last deglaciation in the modelled extent of the geographical range also indicated increases for most species, though their proportion was lower than for N e . However, we could find no association across species between the direction or magnitude of change in N e and range reconstructions.
Species with large ranges at present are the probable "winners" in their response to climatic change and are more likely to show an increase in range and population size from the LGM. Widespread species might also be more likely to have been sampled, and indeed, we found that, based on BirdLife breeding range data, the species studied here have significantly larger modern-day ranges than Holarctic species as a whole (Wilcoxon rank sum p < .001, Figure   S13). Most species sampled also showed an increase in range size based on SDMs. However, the proportion of expanding species according to SDMs was much lower than that observed for BSPs, thus failing to fully explain the ubiquity of expanding BSPs, and there was no statistical association between changes since the LGM as reconstructed by BSPs and SDMs.
Colonization bottlenecks during range shifts can, in principle, lead to an increase BSP irrespective of the overall change in range size, as long as migration is low enough (i.e., if the BSP captures the local dynamics in a given population/small geographical region rather than the whole range). However, if this mechanism were important, we would expect species exhibiting an increase in BSP despite a range contraction to be associated with large shifts. In fact, this is not the case for the majority of species with increasing BSP profiles and decreasing or stable SDM profiles, ( Figure S14). Therefore, colonization bottlenecks do not seem to explain the ubiquity of expanding BSPs in our data set.
Increases in migration can also lead to an increasing BSP without any change in census population sizes. The potential role of migration in producing counterintuitive N e estimates when assuming panmictic populations for a whole species has received much attention in the context of interpreting cross-coalescence (e.g., Multiple Sequentially Markovian Coalescent (MSMC)) profiles (Mazet et al., 2015). However, it is difficult to envisage a scenario where migration would increase significantly in the face of a range contraction; the effect of migration is more likely to be seen during an expansion, when previously isolated fragments are reconnected. Thus, it seems unlikely that migration can explain our results.
A final, more likely but difficult-to-test explanation is that population densities have increased since the LGM. Thus, even for species that have experienced a range contraction, there might have been changes in local population dynamics such that the average density is higher at present. This decoupling of range extent and density makes interpretation of the SDMs and BSPs very challenging. To resolve them, there is a need to use species abundance models that explicitly predict population densities rather than presence/absence (Howard et al., 2014;Johnston et al., 2015;Potts & Elith, 2006 with other species, availability of appropriate prey types) which might determine the distribution of a species (Araújo & Guisan, 2006). If those neglected factors are important enough and are not adequately predicted by (i.e., correlated with) a climatic variable, the SDMs will fail to properly reconstruct the past range. Furthermore, SDMs assume that the niche of a species has not changed through time, but some species might have adapted to different conditions over the last 20,000 years (Townsend, 2011;Veloz et al., 2012).
Thus, we caution that some of the mismatch between BSPs and SDMs might also be attributed to shortcomings of the latter.
The timings we find for when population expansions occurred agree broadly with those of changes in climate after the LGM. Dates were based on mitochondrial mutation rates calibrated for body size and based on a calibration set that included relatively young species splits (Nabholz et al., 2016), and thus likely to give faster mutation rates (i.e., less affected by selection) that were appropriate for within-species analysis (Ho, 2007;Howell et al., 2003;Penny, 2005).
Using a calibration set that included older splits (Nabholz et al., 2016) would lead to much older (well before the last glaciation), and thus less realistic changes in N e . However, we strongly caution that mutation rates calibrated by bird body size, whilst the best available option for comparative analysis, are likely to be noisy, and individual species estimates should not be overly interpreted. Ideally, we would need taxon-specific mutation rates (Hope et al., 2014) which are not available for the number of species investigated here. Having said this, the fact that species from the same habitat tend to yield broadly similar profiles gives us confidence that the relative timings are robust, even if the absolute values have room for improvement.
We found that changes in N e for wetland-associated species have occurred more recently, and have been significantly larger, than those from terrestrial habitats. Although the significance is marginal when based on point estimates for the date of most rapidly changing size, the finding is supported by MSI analysis, which also reveals a similar pattern. Compared to many terrestrial habitats, wetlands tend to be less stable: factors such as local water levels, meltwater from retreating ice sheets, and soil erosion all play into wetland development and could have delayed their stabilization after the LGM.
Although reconstruction of wetland environments and the modelling of wetland recovery is difficult (Fan & Miguez-Macho, 2011;Kaplan, 2002;Lafleur, 2008;Valdes et al., 2005), analysis of pollen across Eurasia shows that species associated with wetlands such as sphagnum moss and alder trees exhibit much later expansions compared with terrestrial species (Allen et al., 1999;Giesecke et al., 2017). Indeed, the expansion of alder relative to other trees (Allen et al., 1999) matches closely the later expansion we find for wetland vs. nonwetland birds.
BSP analysis is powerful but depends on a number of assumptions that are rarely met in real data, most notably the use of a random sample of individuals drawn from a panmictic population (Chikhi et al., 2010;Heller et al., 2013;Pannell, 2003). Consequently, most profiles should be seen as approximations that are easy to over-interpret (Grant, 2015;Miller et al., 2018). However, increasing numbers of public domain data sets allow studies based on characteristics averaged across multiple profiles constructed from species or populations that share habitats or traits. This approach has its own challenges and the data need stringent filtering. Here, we constructed and inspected network diagrams for each data set, identifying species with genetic outliers and evidence of strong population substructure to be divided or excluded. It is noteworthy that the use of public domain data sets can additionally restrict the level of control a study has on sampling scheme, a key tool for limiting the impact of violating the panmixia assumption (Chikhi et al., 2010;Heller et al., 2013). We explored the impact of sampling scheme on five of the data sets from this study and found the shape and trend of the median line to be generally robust to dramatic alterations in sampling scheme ( Figure S5A-E). Although we recognize that structural biases undoubtedly exist in the data uploaded to GenBank, this gave us confidence that structural bias was not driving the demographic patterns recovered. However, the large number of species investigated allowed us to see a clear pattern of population expansion in almost all species following the LGM, irrespective of their range dynamics inferred from SDMs, and a tendency for the expansion to occur later in wetland species. The near-ubiquitous signal of expansion suggests a possible decoupling of range size and local densities, implying a need for carefully interpretation of BSP and SDMs to describe species-wide responses.

ACK N OWLED G EM ENTS
We thank Professor Brian Huntley and a group of palaeoecologists gathered to honour him upon his retirement for pointing out the similarity of our finding of a later increase in effective population size of wetland birds and a similar delay in the expansion of some wetland plants during the last deglaciation.

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