Adding ecology into phylogeography: ecological niche models and phylogeography in tandem reveals the demographic history of the subalpine warbler complex

ABSTRACT

Capsule: This study documents evidence of interglacial refugia during the Last Interglacial for birds in the Mediterranean region, and emphasizes the importance of the Last Interglacial on the geographic distribution and genetic structure of Mediterranean species.

Aims: We focused on the historical biogeography of the subalpine warbler complex: Subalpine Warbler Sylvia cantillans and Moltoni’s Warbler Sylvia subalpina; we tested if this Mediterranean bird complex shared a similar demographic fate as the present-day widespread species in the temperate zones of Europe, through the late Quaternary glacial-interglacial cycles.

Methods: An ecological niche model was developed to predict the geographic distribution of the subalpine warblers under the past (the Last Interglacial and the Last Glacial Maximum) and the present bioclimatic conditions. Additionally, Bayesian Skyline Plot analysis was used to assess effective population size changes over the history of the subalpine warbler complex.

Results: During the Last Glacial Maximum, the subalpine warblers almost reached their current distribution in the Mediterranean region; yet, unlike the widespread temperate bird species, they survived the Last Interglacial in allopatric refugia in the Mediterranean region.

Conclusion: A unique biogeographic pattern was revealed, indicating the importance of the Last Interglacial on current distributional patterns and demographic histories of common bird species in the Mediterranean region. This study suggests that Mediterranean biogeography is far more complex than previously assumed, and so deserves further study and more attention.

Understanding the impact that climate fluctuations have had on species’ distributions and their demographic histories is a key element in biogeography. Numerous studies have shown that the Quaternary glacial-interglacial cycles had an immense effect on dispersal patterns and demographic structures of species, particularly those of vertebrates (Hewitt 1999). Birds in the Palearctic is one such well-studied group, especially in terms of their distributional and demographic histories during the last 130,000 years (Dai et al. 2011, Zhao et al. 2012, Hung et al. 2012, Pellegrino et al. 2014, Perktaş et al. 2015, Kamp et al. 2019). Of all major historical biogeographic events that crafted species’ demographic histories, special emphasis should be given to the last ice age, which reached its maximum 26,500 years before present (Clark et al. 2009). European bird species, that are widespread and well-adapted to present temperate climates, survived in small restricted and climatically favourable areas referred to as refugia during this time period. Some southern parts of Europe that have been postulated as glacial refugia include Iberia, Italy, the Balkans, and Anatolia, from where these present-day species might have expanded their ranges northwards (Cooper et al. 1995, Hewitt 1996, Perktaş et al. 2011, Pons et al. 2011, Perktaş & Quintero 2013). However, no detailed study discusses if the biogeographic history of birds currently restricted to these Mediterranean presumptive refugial sites is at all concordant with that of birds that are widespread across Europe today.

The subalpine warbler complex is a species complex comprising Moltoni’s Warbler Sylvia subalpina and three subspecies of Subalpine Warbler Sylvia cantillans (del Hoyo 2006, McInerny et al. 2018), with an intriguingly restricted distributional pattern, providing an ideal case to explore the aforementioned biogeographic question (Shirihai et al. 2001). This species complex has been well studied in terms of vocalization, phylogeography, and taxonomy (Brambilla et al. 2008a, 2008b, 2010, Svensson 2013). Throughout its distributional area, a broad spectrum of Mediterranean habitat types, such as garrigue, maquis, and oak woodlands, is key for the survival of at least 75% of the breeding populations (Tucker & Evans 1997). According to Svensson (2013) and Dickinson & Christidis (2014), morphological and molecular characters support a three-way split of the complex: the Western Subalpine Warbler Sylvia inornata, the Moltoni’s Warbler S. subalpina, and the Eastern Subalpine Warbler S. cantillans (Figure 1). Following McInerny et al. (2018), here we treat Eastern Subalpine Warbler Sylvia cantillans albistriata and Western Subalpine Warbler Sylvia cantillans cantillans as distinct subspecies of Subalpine Warbler, Moltoni’s Warbler as a distinct species and the three taxa together as the ‘subalpine warbler complex’.

Figure 1. Distribution range of the subalpine warbler species complex. Sampling points and hypothetical species limits were re-drawn based on Brambilla et al. (2008b) and Svensson (2013).

Earlier studies on the subalpine warbler complex analysed mitochondrial DNA (mtDNA, the cyt-b gene) along much of its distributional range (Brambilla et al. 2008b, 2010). These studies, for the most part, solved taxonomic problems within the complex, yet did not delve into a detailed review of the demographic history of the complex. Integrating distributional analyses (e.g. ecological niche modelling) with phylogeographic studies can offer a better understanding of the complex demographic patterns of species (Gür 2013, Perktaş et al. 2015, 2017). As new questions arise concerning the ecological and evolutionary mechanisms responsible for shaping the demographic history of species, examining the distributional shifts and ecological niche differentiations is an effective approach to explore plausible answers. In this paper, we aim to develop detailed distributional projections from ecological niche models and to integrate these results with demographic analyses based on mtDNA data of the species complex. Therefore, this study can be thought of as an extension of the work of Brambilla et al. (2008a, 2008b, 2010), with distributional projections (i.e. ecological niche modelling). With our results, we discuss whether Mediterranean bird species have a similar demographic history as present-day widespread species in the temperate zones of Europe.

Methods

Ecological niche modelling

Input data

We collected occurrence data for our species complex from eBird (www.ebird.org) data portal for the last 57 years. Considering the migratory nature of this complex throughout its distribution range, we only collected data for the breeding season (Cramp 1992, Shirihai et al. 2001), i.e. the time period ranging from early May to the end of June, since the records of nesting and hatching behaviour are concentrated within this time period (Shirihai et al. 2001, Thévenot et al. 2003). We collected 2910 occurrence points for subalpine warblers, including that of the Western Subalpine Warbler from Iberia and southern France, and the Eastern Subalpine Warbler from Anatolia, Aegean Islands, and the Balkans. We did not use any e-Bird-based occurrence points of the Moltoni’s Warbler and the southern and central Italian population of the Eastern Subalpine Warbler. Even though the species limits based on genetic data were quite obvious, the vague geographic limits in occurrence data permitted us only to use the published occurrence points (Brambilla et al. 2008b) for the Moltoni’s Warbler (n = 11) and the Italian population of the Eastern Subalpine Warbler (n = 10).

To get the best-homogenized and unbiased distribution of occurrences avoiding spatial autocorrelation, we first eliminated the duplicates and artefacts by manually checking all data points. Then we used a 50 km distance-filter (in SDM Toolbox [Brown 2014] in ArcGIS ArcGIS version 10.5.1. [http://www.arcgis.com]) to further homogenize the occurrence dataset. The process yielded a tally of 59 records from Anatolia, the Aegean Islands, and the Balkans for the Eastern Subalpine Warbler, and 110 records from Iberia and France for the Western Subalpine Warbler. These rarefied records were then subjected to the ecological niche modelling process.

We used the WorldClim climatic data archive (Hijmans et al. 2005) to obtain the Climate data (1950–2000) at a spatial resolution of 2.5 arc min (4.6 km at the equator). We obtained parallel data sets for the Last Glacial Maximum (LGM; approximately 21 kilo years before the present), including model outputs for the CCSM4, MPI-ESM-P, and MIROC-ESM simulations, and for the Last Interglacial [LIG; approximately 120–140 kilo years before the present, see Otto-Bliesner et al. (2006) for details]. These climate data sets offer three temporal pictures of climate over the Pleistocene and recent interglacial, glacial, and present-day conditions. We masked these bioclimatic data to include only −19° to 50°E and 20° to 50°N.

Ecological niche models

We combined occurrence data for the subalpine warbler complex and restricted the first model output to the species’ historically accessible areas by means of dispersal (i.e. the M area in the BAM diagram framework in Soberon & Peterson (2005) and Barve et al. (2011)). We hypothesized the accessible area for the lineages of subalpine warblers based on the genetic limits, taxonomic evaluation (Svensson 2013), and main habitat types in the Mediterranean region (e.g. Mediterranean forest and shrubland; Tucker & Evans 1997). To characterize M areas for the subalpine warbler complex, we created a minimum convex polygon around the occurrence points, resulting in a 200 km buffer zone.

Reducing climate data based on species ecology (Cramp 1992, Shirihai et al. 2001) is often preferred to reduce the candidate predictor in cases of limited available occurrence points for any given species (Elith & Leathwick 2009). In our case, we just had a handful of occurrence points for the Moltoni’s Warbler. Thus, we reduced the climatic data and used only the following six climatic variables: Bio1 (Annual Mean Temperature), Bio4 (Temperature Seasonality), Bio5 (Max Temperature of Warmest Month), Bio12 (Annual Precipitation), Bio14 (Precipitation of Driest Month), and Bio15 (Precipitation Seasonality). These six bioclimatic variables characterize well the distribution area of the subalpine warbler complex during the breeding season (also see Tellería et al. 2016 for a similar methodological approach). The distribution range of the subalpine warbler complex is almost equivalent to that of the Mediterranean basin biodiversity hotspot. Thus, we obtained the six bioclimatic variables for the analysis using the approximate limit of the Western Palearctic region. We then performed a model calibration, which is a crucial step before the final projection (Ingenloff et al. 2017), using the maximum entropy machine-learning algorithm in the software Maxent version 3.3.3k (Phillips et al., 2006), developing the models for the M areas of the species complex with different feature types combinations (linear, quadratic, product, hinge, and threshold), and regularization multipliers (0.1, 0.2, 0.5, 1, 2, 5, 10). We compared the models using the Akaike Information Criterion corrected (AICc) for small sample sizes using the software ENMTools version 1.4.4 (Warren et al. 2010). The model with the lowest AICc score was selected as the most accurate one. To test the model significance for each species, we used 50% of test data which was based on random subsamples, and a partial receiver operating characteristic (pROC) approach via the PARTIALROC function in the R package ‘ENMG_ADGETS’ (Barve & Barve 2013). Finally, we ran all models with 10 replications for (1) the present day, (2) the LGM (CCSM4, MIROC-ESM and MPI-ESM-P), and (3) LIG across the whole Mediterranean region, including the range for the Mediterranean basin biodiversity hotspot. At the last step, the 10-percentile training presence threshold approach was used to convert model outputs to binary predictions (Radosavljevic & Anderson, 2014, also see Ülker et al. 2018).

Historical demography

The Bayesian skyline plot

We examined the demographic events throughout the history of the subalpine warbler complex based on the cyt-b gene of mtDNA obtained from Brambilla et al. (2008b). Using these data, we performed a Bayesian skyline plot analysis (Heled & Drummond 2012) in BEAST version 1.7.5 (Drummond et al. 2012). This analysis uses coalescent approaches to estimate effective population size changes through time. In order to compare the results from this analysis with those of the niche modelling ones (see below), we combined mtDNA sequences (Brambilla et al. 2008b) from all populations in these assessments (Perktas et al. 2015).

Best-fit substitution models for the skyline plots were identified for the sequences in MEGA version X (Kumar et al. 2018). To date demographic events, we used the 2% per million years mutation rate for the cyt-b for (Drovetski 2003, Brito 2005, Pereira & Baker 2006, Weir & Schluter 2008). We used a generation time of one year for the subalpine warblers, based on Cramp (1992). Multiple independent skyline plot runs were performed using the following parameters: linear models 100 × 10⁶ steps, parameters sampled every 10,000 steps, and burn-in 10%. We used a strict molecular clock model under a uniform distribution of priors. Effective sample size values of the parameters were over 200 for each run, indicating that 100 million steps were sufficient to assess population fluctuations over the history of the subalpine warbler complex.

Isolation-by-distance

Genetic distance among populations of subalpine warblers was estimated using Fst statistic using DNAsp version 6.12.01 (Rozas et al. 2017). Then, Fst was used to generate a genetic distance matrix to evaluate genetic isolation-by-distance pattern. A matrix of genetic distances between 12 populations was estimated based on the formula, Fst/(1-Fst). A matrix of geographic distance (km) was estimated by the geographic distance matrix generator (Ersts, P.J. [Internet] Geographic Distance Matrix Generator, version 1.2.3. American Museum of Natural History, Center for Biodiversity and Conservation. Available from http://biodiversityinformatics.amnh.org/open_source/gdmg). A Mantel test with 10,000 random permutations was performed between the genetic and the logarithmic geographic matrices (Slatkin 1993, Rousset 1997).

Results

Ecological niche modelling

The model calibration result based on AICc values showed that the best model included a regularization multiplier of 2 for the subalpine warbler complex and five feature types: linear, quadratic, product, hinge and threshold. On the basis of the model calibration results, the ecological niche modelling results provided a high area under the curve value for the training data (AUC = 0.751, sd = 0.005), indicating that the model had a very good predictive ability. Further analysis, partial ROC statistics showed that the model provided predictions of the geographical range of the subalpine warbler complex that were significantly better than random expectations (P < .05). Response curves of bioclimatic variables showed that the subalpine warblers occurred in an environment characterized by particularly temperature seasonality and annual precipitation. Two variables (BIO4 and BIO12) substantially contributed to the model (>65%).

Under the present bioclimatic conditions, the model prediction was compatible with the known distribution range of the subalpine warbler complex (Figure 2); however, ample overprediction in the present model showed suitable areas for the species in parts of Middle and Northern Europe. Hence, this prediction suggests that the subalpine warbler complex is ‘almost’ at equilibrium with the climate (‘almost’ at equilibrium due to overprediction in the present prediction). The models were unstable under certain past bioclimatic conditions, in particular, during the transitional periods between the Last Interglacial and the Last Glacial Maximum (Figure 2). Of particular interest is that our study revealed that the distributional pattern of the subalpine warbler complex in the Mediterranean region during the Last Interglacial was extremely restricted and limited to allopatric refugia (Figure 2).

Figure 2. Ecological niche model showing the geographical distribution of the subalpine warbler complex under the Present (Pre., 1950–2000) and the reconstructed past (the Last Glacial Maximum, 22 KYBP; and the Last Interglacial, 130 KYBP) bioclimatic conditions.

Historical demography

The GTR + G+I model was identified as the best fit for the Bayesian skyline plots (corrected AICc = 5222.498). The Bayesian skyline plot indicated a distinct pattern of population contraction followed by population expansion, over the late Quaternary. Our analyses unanimously indicated that both population contraction and expansion happened before the Last Glacial Maximum (Figure 3).

Figure 3. The Bayesian skyline plot for the subalpine warbler complex. Based on 2% mutation rate of the cyt-b gene of mtDNA, the Last Glacial Maximum line is located on the time axis.

In addition to the skyline plot analysis, genetic distances among populations showed a positive correlation to geographic distance, which signalled evidence for isolation-by-distance (r = 0.338, P < .05). Additionally, some population comparisons exhibited high F_ST values, perhaps suggesting high genetic differentiation in separate geographies in the Mediterranean region.

Discussion

In this study, we integrated ecological niche modelling with phylogeographic analyses in an effort to describe the effects of late Quaternary climate fluctuations on the demographic history of the subalpine warbler species complex. By employing ecological niche modelling, we were able to assess climate-driven distributional shifts of the species complex (under the species-climate equilibrium through time; Nogués-Bravo 2009, Perktaş et al. 2015, Perktaş & Gür 2015). Phylogeographic techniques allowed us to assess the geographic isolation of closely related lineages, using the geographic variation of genetic diversity, in order to explore the dynamics of the species’ range (e.g. expansion–contraction; Avise 2000). Phylogeographic assessment based on mtDNA data generally reveals the recent history (e.g. the Last Glacial Maximum, the Last Interglacial) of a species (Hewitt 1999). As a consequence, results from mtDNA phylogeography are often compatible with ecological niche modelling outcomes (Gür 2013, Perktaş et al. 2015). From this viewpoint, we established an integrative perspective to test the demographic history of the subalpine warbler complex, to explore the parallels between the demographic history of a Mediterranean bird and a present-day widespread species complex in the temperate zones of Europe.

The assessment on the species’ demographic history suggested that populations of the subalpine warbler complex survived in allopatric interglacial refugia along the Mediterranean region. Ecological niche models further revealed that the species complex nearly reached its present distribution range before the Last Glacial Maximum, such interglacial refugia patterns have not been described for the Mediterranean region previously, rendering this a significant result. Earlier studies in the Western Palearctic region discussed the effects of the Last Glacial Maximum on the formation of three allopatric refugia in the Mediterranean region, without any mention of the Last Interglacial and its plausible influences on the process (Brito 2005, Pons et al. 2011, Perktaş et al. 2011). The only other evidence for such interglacial refugia were reported from Anatolia, for the Kruper’s Nuthatch Sitta krueperi (Perktaş et al. 2015) and the Anatolian Ground Squirrel Spermophilus xanthoprymnus (Gür 2013).

Even the well-known ‘refugia-within-refugia’ model failed to describe or discuss the effect of the Last Interglacial in the Mediterranean region (Gómez & Lunt 2007). The model focused on the areas with high genetic diversity (i.e. high allelic richness) in each separate refugium in the Mediterranean region (Iberia, Italy, and the Balkans) and discussed the causes of variation of genetic diversity under the effect of the Last Glacial Maximum in the context of an ‘expansion-contraction’ model (Provan & Bennett 2008) in the Western Palearctic Region. However, the vast majority of phylogeographic studies from the region (Canestrelli et al. 2007, 2008, 2012) focused on different organisms that had limited distributions in these well-known refugia, and the results were mostly coherent with a refugia-within-refugia model. Yet, outcomes of these studies (on the demographic history) based on mtDNA showed a substantial expansion pattern, which started before the Last Glacial Maximum, similarly as in the case of this study.

Our ecological niche modelling prediction for the present distribution of the subalpine warbler complex was largely concordant with its known distribution. The only artefact in the model was the overprediction of the taxa in the middle and further north in Europe. However, some recent records showed that subalpine warblers could breed further north of the Mediterranean Region due to current climate change effects (see e-Bird records). Therefore, this prediction indicates the species-climate equilibrium and stability of ecological niches of the subalpine warblers based on bioclimatic data (Nogues-Bravo 2009), which could add confidence to the model. Model output showed that the most prominent bioclimatic variables are temperature seasonality (BIO4) and annual precipitation (BIO12). Areas with relatively high precipitation and temperature seasonality in the Mediterranean region were more suitable areas for this species complex. High precipitation could probably affect the productivity and, therefore, could affect food availability for migrating bird species. Hence, relatively dry but productive areas could be accepted as a suitable habitat for subalpine warblers. Our modelling results do not confirm glacial range contraction and interglacial range expansion (i.e. expansion–contraction model, Provan & Bennett 2008) for subalpine warblers. However, they do confirm the suggestion of Stewart et al. (2010) that continental species are adapted to drier climates with greater seasonality, and these species have been in refugia during the interglacial periods. Therefore, this study adds to the evidence that interglacial range contraction might be proposed as a general biogeographic pattern not only for high altitude species (Provan & Bennett 2008, Stewart et al. 2010) but also for Mediterranean species, such as the subalpine warblers.

The subalpine warbler complex is phylogeographically structured into four mtDNA lineages that diverged most probably before and at the beginning of Pleistocene, approximately between 2.5 and 1.7 million years ago based on 2% sequence divergence (Brambilla et al. 2008b). These authors suggested a constant demographic history for the eastern lineage, without a detailed discussion for the estimate of population expansion for the three clades, or whether the expansion event was concordant with the demographic history of widespread bird species in temperate regions in Europe (Brito 2005, Pons et al. 2011, Perktaş et al. 2011, Perktaş & Quintero 2013). In addition to this, in four clades, almost all haplotypes were closely related, yet each clade was geographically isolated (Avise’s phylogeographic category III, Avise 2000). This result indicates that contemporary gene flow among clades has been low enough, promoting genetic divergence among the populations of the subalpine warblers. Taken all together, the Bayesian Skyline Plot and an isolation-by-distance pattern might be consistent with both ecological niche modelling and observed mtDNA phylogeographic pattern (Brambilla et al. 2008b). Thus, the most reasonable explanation for the subalpine warbler complex demographic history might be the long-term isolation within three different interglacial refugia in the Mediterranean Region, followed by a recent expansion from these refugia.

Petit et al. (2005) suggested that long-term persistence of isolated populations have been a common phenomenon in the Mediterranean region. Thus, unique patterns of geographic variation of genetic diversity of different organisms (such as the subalpine warbler complex) have not been unexpected in this region. Based on IPCC (2001), climate change is expected to be the most prominent factor for extinction risk in the Mediterranean region. This prediction concurs with the findings in this study. Past climate fluctuations and corresponding range dynamics of subalpine warblers left their marks in its genetic population structure, and a comprehensive ecological assessment on this pattern helped unravel their robust history in the late Quaternary.

The subalpine warbler complex mostly use scrubland and shrubland vegetation types but not forests in the Mediterranean region (Shirihai et al. 2001, Brambilla et al. 2006). After the Last Interglacial, this region experienced dryness during the interglacial-glacial transition phase, as a result, during the Last Glacial Maximum, forest cover significantly declined throughout the Mediterranean region, according to pollen data (Fletcher & Goñi 2008, Combourieu-Nebout et al. 2009). The decline in forest biomass during the interglacial-glacial transition phases and more pronouncedly in the glacial periods is a well-known consequence of the effect of glacial-interglacial cycles on ecosystems (Birks & Birks 2004). With the onset of the Holocene, although sclerophyllous oak forests expanded in the region, the Mediterranean shrublands remained as one of the important vegetation types (Pons & Reille 1988, Fletcher & Goñi 2008, Combourieu-Nebout et al. 2009). The historical vegetation changes since the Last Interglacial in the Mediterranean region suggest that the subalpine warblers had more suitable habitat during the Last Glacial Maximum and the Holocene. All the conclusions from previous studies on the vegetation history in the Mediterranean region were concordant with our assessments on the historical demography of the subalpine warbler complex. However, detailed phylogeographic assessments (e.g. multi-loci assessment, genomic perspective, etc.) on the subalpine warblers will add to our understanding of the historical biogeographic patterns for birds in the Mediterranean region.

Acknowledgements

The authors are grateful to University of Kansas, ENM group members to their help during the preparation process of this manuscript. The authors extent their special thanks to Banu Ş. Önder, an anonymous reviewer, and the editor for their valuable comments to earlier version of this manuscript.

ORCID

Utku Perktaş http://orcid.org/0000-0002-5988-7289