Abstract
Global warming is increasingly exacerbating biodiversity loss. Populations locally adapted to spatially heterogeneous environments may respond differentially to climate change, but this intraspecific variation has only recently been considered when modelling vulnerability under climate change. Here, we incorporate intraspecific variation in genomic offset and ecological niche modelling to estimate climate change-driven vulnerability in two bird species in the Sino-Himalayan Mountains. We found that the cold-tolerant populations show higher genomic offset but risk less challenge for niche suitability decline under future climate than the warm-tolerant populations. Based on a genome-niche index estimated by combining genomic offset and niche suitability change, we identified the populations with the least genome-niche interruption as potential donors for evolutionary rescue, i.e., the populations tolerant to climate change. We evaluated potential rescue routes via a landscape genetic analysis. Overall, we demonstrate that the integration of genomic offset, niche suitability modelling, and landscape connectivity can improve climate change-driven vulnerability assessments and facilitate effective conservation management.
Similar content being viewed by others
Introduction
Anthropogenic climate change is one of the primary drivers of environmental change and global biodiversity loss1,2. Rapid climate changes can lead to shifts in species ranges, population decline and even extinction if the changes exceed the physiological tolerance of organisms3,4. As these effects become increasingly profound in light of global warming, how organisms can respond to environmental changes is becoming a focus of basic research5,6. Modelling changes in the distribution range, suitable climatic conditions and vegetation types of species under different climate scenarios has provided considerable insights into the impacts of climate change on biodiversity (e.g. refs. 6,7,8,9,10). However, as these models rely solely on abiotic and biotic environmental changes, ecologic genomics, which investigates how the genomic variants of a species vary along current environmental gradients and how much genetic change has been required to keep up with climate change-driven environmental changes, has only recently been integrated into modelling species responses to climate change11,12,13.
These ecological genomic studies have shown different genotype-climate associations between populations and suggested that local populations may respond differentially to climate change11,12,13. Nevertheless, intraspecific variation has rarely been incorporated into ecological niche models that often assume a uniform climate response between populations11,14,15,16. There is an urgent need to incorporate intraspecific genomic variation in modelling habitat suitability in the context of climate change, as such information is necessary for understanding fine-scale estimates of climate change-driven vulnerability (e.g. ref. 17).
Mountainous areas harbour exceptional biodiversity and endemism but are highly vulnerable to climate change18. This is because the complex topography within a rather small geographical area leads to dramatic ecological stratification in mountain regions. Species are often confined to spatially heterogeneous environments and locally adapted to diversified climate conditions18,19,20. These population-specific adaptations likely drive different responses under climate change because the populations likely track their own optimal environmental conditions20. Despite having high potential, this intraspecific variation has not been considered in vulnerability estimates driven by climate change in mountainous species. Herein, we integrate ecological genomics and niche modelling to investigate the population-specific responses to future climate change in two mountainous birds in the Sino-Himalayan Mountains. Together with the mountains of southwestern China (hereafter the Southwest Mountains) and the East Himalayan hotspots21, these areas harbour unique temperate biodiversity of high elevation but are now under serious threat from the rapidly changing climate22.
Our focal species, the Green-backed tit (Parus monticolus) and Elliotâs laughingthrush (Trochalopteron elliotii), are two mountainous birds that are mainly found in the Sino-Himalayan Mountains. The two species are distributed across different elevation ranges; P. monticolus is a mid-elevational bird commonly found between 1500 and 2500âm above sea level (a.s.l.), whereas T. elliotii lives at higher elevations between 2000 and 4500âm a.s.l.23. Previous phylogeographic studies have revealed that these birds are often restricted to ecologically and topographically heterogeneous areas, suggesting that they have adapted locally to disparate ecological zones24,25,26. As such, these species offer an excellent study system to examine the intraspecific variation in genotype-climate associations and how this would influence climate change-driven vulnerability.
Here, we combine ecological genomics and niche modelling to evaluate the population-specific responses of the two species to future climate change. We consider both the genomic offset (i.e., a measure of the mismatch in genotype-climate association between current and potential future climates11) and niche suitability change (i.e., a measure of the difference in niche suitability between current and future ecological niches) for potential evolutionary rescue; that is, only populations with minor genome-niche interruption can serve as a desired store of speciesâ survival for climate change (i.e., being a donor). We then use a landscape genetic approach to evaluate potential rescue costs (a framework is demonstrated in Fig. 1). Overall, as genomic offset and niche suitability change show different population-specific vulnerability estimates under future climate change, we emphasise the importance and strength of combining ecological genomics, niche modelling and landscape connectivity to guide effective mitigation efforts on increasingly threatened biodiversity.
Results
Species distribution
The distribution ranges of T. elliotii and P. monticolus cover the Sino-Himalayan Mountains and the mid-elevational mountains in Central China (Fig. 2a). Across their ranges, the average elevation declines from west to east and from north to south. The temperature and precipitation decrease from southeast to northwest and from south to north. The current distributions of T. elliotii and P. monticolus partially overlap but differ in certain ecological zones and elevation ranges (Fig. 2b). They both inhabit in the southwest mountainous zone, the western mountainous plateau zone, and parts of the Loess Plateau zone and the east meadow zone. However, T. elliotii also occurs in the southern Tibetan zone, whereas P. monticolus also lives in the eastern Himalayan zone. Given the considerable environmental and climatic variation across their distribution ranges, the two species provide an excellent opportunity to study intraspecific responses to environmental change following future climate change.
De novo genome assembly and annotation of T. elliotii
We generated a reference genome for T. elliotii to facilitate mapping of the resequencing data. In total, 157 gigabases (Gb) data were sequenced for a T. elliotii individual. After cleaning and quality control of the raw data, all qualified data were assembled into a genome of 1.12âGb with an N50 contig length of 95âKb and N50 scaffold length of 2.702âMb. Using a homologue-based approach, we annotated 17,585 protein-coding genes. Using Benchmarking Universal Single-Copy Orthologs (BUSCO, aves_odb927) as the reference gene set, we estimated that the assembly contains 90% complete single-copy BUSCOs, 2% complete duplicated BUSCOs, and 4.5% fragmented BUSCOs. We used this genome for the subsequent ecological genomic analyses.
Intraspecific variation in genotype-climate association
We integrated population genomics and climatic data to identify climate-associated genomic variation. We generated genome-wide resequencing data from 55 and 58 individuals across the distribution ranges of T. elliotii and P. monticolus, respectively (Supplementary Data 1). It is noteworthy that the strategy to obtain deep-sequenced whole-genome datasets (at average coverage of 19.14à and 15.48à for T. elliotii and P. monticolus, respectively, Supplementary Data 2) allowed us to capture maximal genomic variation, in contrast to most previous studies that are based on reduced-representation genomic datasets (see ref. 13). Using a de novo genome of T. elliotii generated in this study and an assembly of the great tit (Parus major) genome28, we identified 10.3 and 3.9 million single-nucleotide polymorphisms (SNPs) in T. elliotii and P. monticolus, respectively.
We used gradientForest29, a machine-learning regression tree-based approach, to first identify the climatic variables that are most closely associated with the genetic variation in the two species and then transform multidimensional climatic variables into multidimensional genetic space (i.e., turnover in allele frequencies)11. We found that climatic variables related to seasonal changes in temperature and precipitation were strongly associated with the genomic variants in both species. Of the 19 climatic variables tested (Supplementary Table 1), the top five uncorrelated variables identified were BIO3 (isothermality, i.e., mean diurnal range divided by the temperature annual range), BIO18 (precipitation of the warmest quarter), BIO9 (mean temperature of the driest quarter), BIO19 (precipitation of the coldest quarter) and BIO5 (max temperature of the warmest month) for P. monticolus, while the top five climatic variables for T. elliotii were BIO2 (mean diurnal range, i.e., the mean of the monthly differences between the maximum and minimum temperatures), BIO10 (mean temperature of the warmest quarter), BIO7 (temperature annual range), BIO19 and BIO4 (temperature seasonality) (Fig. 3a, b and Supplementary Table 2).
Out of 50,000 randomly extracted SNPs, we identified 5446 and 7294 SNPs that were significantly associated with the top climatic variables (R2â>â0, see Methods) for T. elliotii and P. monticolus, respectively. Using an aggregate community-level turnover function across these SNPs, we visualised the genotype-climate turnover surface across the distribution ranges of the two species using the first three principal components (PCs) of the gradientForest outputs. We found that the genotype-climate associations vary from the western parts (i.e., the southern Tibetan zone and eastern Himalayan zones) to the eastern (i.e., the western mountainous plateau zone) and southern parts of the distribution ranges (i.e., the southwest mountainous zone) (Fig. 3a, b). These results suggest that the two species show intraspecific variation in their genotype-climate associations, which are likely subject to local adaptation to heterogeneous climatic conditions.
Genomic offset to future climate change
Based on the genotype-climate associations across the distribution ranges of the two species, we predicted which part of their ranges might be most vulnerable to future climate change. We used gradientForest29 to calculate genomic offset as the Euclidean distance between current and projected future genomic compositions. Genomic offset is thus a measure of how much genetic change is needed to adjust to the new climate conditions, and the populations with the greatest genomic offset are those that have to adjust the most. To account for variability between climate models, we predicted the genomic offset for each species using four different climate models, MPI-ESM-LR30, CCSM431, MICRO-532 and CNRM-CM5-233 (Supplementary Table 3), under the moderate scenario (representative concentration pathway, RCP 4.5âW/m2) and the worst scenario (RCP 8.5âW/m2) of greenhouse gas emission trajectories in 2050 and 5070, respectively. Genomic offset estimates revealed similar but magnitude-dependent spatial patterns under the four RCP emission scenarios considered at the 2050 and 2070 horizons (Supplementary Fig. 1). These results consistently showed that the populations in the western part (T. elliotii in the southern Tibetan zone and P. monticolus in the eastern Himalayan zone) exhibited greater genomic offset than the populations inhabiting the eastern (the Loess Plateau zone, eastern meadow zone and western mountainous plateau zone) and the southern parts of the distribution ranges (the southwest mountainous zone, two-tailed Wilcoxon rank-sum test, false discovery rate (FDR)-adjusted Pâ<â0.001, Fig. 3c, d). Altogether, these results suggest that intraspecific variation in genotype-climate associations of the two species likely drives different genomic offsets in response to future climate change.
Intraspecific variation in genotype-climate association drives different degrees of genomic offsets
To investigate how the intraspecific variation in the genotype-climate association has driven the different degrees of genomic offsets in response to climate change, we identified genomic variants that are significantly associated with each of the top climatic variables (as identified by the gradientForest analysis) using three complementary approaches: latent factor mixed model (LFMM34), redundancy analysis (RDA) and distance-based redundancy analysis (dbRDA)35,36. A total of 72 and 798 SNPs were identified by all three methods for T. elliotii and P. monticolus, respectively (Fig. 4a, b). These SNPs are widely distributed across the genomes, with 25 and 204 SNPs located in the coding sequence and promoter regions (5k upstream and downstream of genes) across 23 and 147 genes for T. elliotii and P. monticolus, respectively (Supplementary Table 4). We then carried out functional enrichment analyses and found that these genes are functionally enriched in catalytic and metabolic processes (Supplementary Fig. 2), with some previously documented as being important in climate adaptation. In particular, CRB1, the only climate-associated gene identified in both T. elliotii and P. monticolus, has also been identified in other cold-tolerant vertebrates, i.e., the Adélie penguin (Pygoscelis adeliae)37 and the Emperor penguin (Aptenodytes forsteri)38. The relationships between allele frequency variation in the SNP of CRB1 and the mean diurnal range (BIO2, T. elliotii), as well as the isothermality (BIO3, P. monticolus), suggest a potential role of this gene in adaptation to extreme temperatures (Fig. 4c, d). In addition, we also detected some genes that are related to heat and cold tolerance, i.e., VPS5339, BRS340 and FASN41, as well as adaptation to hypoxic conditions, i.e., MYH742, SLC19A143 and ARHGAP3944.
We then carried out population genetic structure analyses based on these SNPs using Admixture v1.345 to define the climate-tolerant groups. We found Kâ=â3 to have the smallest cross-validation error in both species and thus to be the optimal number of clusters explaining the variation among individuals (Supplementary Fig. 3). Based on the proportion of the respective cluster within each individual (>0.6), we classified 33% of T. elliotii individuals as tolerant to cold-dry conditions (shrubland and coniferous forests in the southern Tibetan zone and southwest mountainous zone), 15% as tolerant to warm-humid conditions (subtropical broadleaf forests in the southwest mountainous zone), and 16% to warm-dry conditions (temperate broadleaf forests in the western mountainous plateau zone). In P. monticolus, 10% of the individuals were classified as cold-dry tolerant (coniferous forests in an eastern Himalayan zone), 40% of the individuals were classified as warm-humid tolerant (subtropical broadleaf forests in the southwest mountainous zone), and 48% of the individuals were classified as warm-dry tolerant (temperate broadleaf forest in the western mountainous plateau zone, Fig. 4e, f).
We then compared genomic offset values among the three groups using outputs calculated from gradientForest analysis (i.e., R2 positive SNP dataset, see Methods). We found that the cold-dry tolerant groups showed greater genomic offsets than the warm-humid and warm-dry tolerant groups (two-tailed Wilcoxon rank-sum test, FDR-adjusted Pâ<â0.001, Fig. 4g, h and Supplementary Fig. 4a). In addition, we implemented a separate gradientForest analysis using the outlier SNPs identified by the three genotype-climate association analyses (i.e., outlier SNP dataset), and the analyses yielded similar results (two-tailed Wilcoxon rank-sum test, FDR-adjusted Pâ<â0.05, Fig. 3i, j and Supplementary Fig. 4b). To further validate this, we also carried out a parallel genomic offset analysis using generalised dissimilarity modelling (GDM46), a distance-based method (i.e., FST) that is less sensitive to unequal sampling sizes, to estimate and compare genomic offsets among the three groups. Our GDM analyses using two datasets (R2 positive SNP and outlier SNP datasets, see Methods) show similar patterns to those of the gradientForest analyses (two-tailed Wilcoxon rank-sum test, FDR-adjusted Pâ<â0.001, Supplementary Fig. 5), suggesting that our genomic offset results are robust to different datasets and methods.
Integrating intraspecific variation into ecological niche modelling
We next investigated whether the predicted changes in niche suitability caused by future climate conditions differed among the three climate-tolerant groups. We carried out ecological niche modelling using an ensemble modelling approach that combines the outputs of different modelling algorithms, including maximum entropy, generalised boosted model, generalised additive model and multivariate adaptive regression splines implemented in Biomod247. We projected niche suitability separately for the cold-dry, warm-humid and warm-dry tolerant groups, respectively. All ecological niche models have great discrimination ability (true skill statistics, TSS, 0.71â0.88; area under the receiver operating characteristic curve, AUC, 0.9â0.96, Supplementary Table 5).
Niche modelling inferred under the different future emission scenarios and decades consistently showed that the three climate-tolerant groups had different degrees of niche suitability change (Fig. 5a, b and Supplementary Fig. 6). Using the current niche suitability value as the benchmark, we calculated the change in the niche suitability index between the current and future climate (NSCâ=âniche suitability index in the future climateâââniche suitability index in the current climate). A negative value indicates a decrease, while a positive value shows an increase in niche suitability under future climate conditions. We found that the NSC decreased to a much smaller degree in those areas that harboured cold-dry tolerant groups than in those areas that maintained the warm-dry and warm-humid groups (two-tailed Wilcoxon rank-sum test, FDR-adjusted Pâ<â0.001, Fig. 5a, b and Supplementary Fig. 6). These results suggest that the areas suitable for warm-tolerant individuals are more likely to decrease suitability under future climate conditions than those areas where cold-dry tolerant individuals live.
Identify potential climate-tolerant populations for evolutionary rescue
Strikingly, while the warm-humid and warm-dry tolerant groups are likely to risk a decline in niche suitability, they are predicted to need less genetic change to cope with future climate conditions. Conversely, while the cold-dry tolerant group is likely to maintain most of its current habitats, it needs to undergo a greater genetic change. To better reflect the intraspecific variation in responses to climate change in the two species, we combined the genomic offset and NSC to estimate a genome-niche index as described in ref. 48 (see also Methods). We considered only niches with increasing suitability (positive NSC value) because these areas will remain hospitable for these species, while other areas with decreasing niche suitability (negative NSC value) will be challenged by habitat suitability decline. Under future climate conditions, the populations in the central areas of the Southwest Mountains, i.e., the junction between different ecological zones, are predicted to experience minor interruptions of genomic offset-based habitat suitability and will thus maintain their current status in the distribution ranges (Fig. 5c, d and Supplementary Fig. 7). Therefore, these populations can be considered potential donors for evolutionary rescue under climate change.
Landscape connectivity establishes evolutionary rescue routes
To explore how landscape features would affect the patterns of gene flow between populations, i.e., the potential for setting up evolutionary rescue routes, we implemented linear mixed-effect modelling using maximum likelihood population effects (MLPE49) to fit five landscape features on the matrix of pairwise genetic distance (FST, see Methods). We found that landscape connectivity was strongly related to a combination of habitat suitability and elevation for P. monticolus and habitat suitability for T. elliotii, supported by Akaike weights (\({w}_{i}\)) and a linear mixed modelling approach (Supplementary Table 6). We then extrapolated these relationships to estimate migration possibilities from the donor populations to the populations that are potentially vulnerable to climate change. In both species, individuals from the donor populations could disperse westward to the south Tibetan zone and eastern Himalayan zone, southward to the southwest mountainous zone, northwards to the Loess Plateau zone and east meadow zone, and eastward to the western mountainous plateau zone (Fig. 5e, f). Considering that the ecological niche modelling predicted a substantial decrease in the niche suitability in the southern, northern and eastern parts of the distribution ranges, a plausible evolutionary rescue would be a westward dispersal under the assumption of a westward expansion of suitable niches even beyond the current range (Supplementary Fig. 8). Consequently, the evolutionary rescue of mountainous species depends not only on genomic offset and niche suitability estimates but also on landscape connectivity of the populations.
Demographic model tests the migration possibility among groups
To further explore the possibility of migration between the groups, we compared three models of gene flow (one with continuous gene flow, one with secondary gene flow and one with no gene flow, Supplementary Fig. 9) between groups using Fastsimcoal v2.650. Our model tests and Akaike information criterion (AIC) comparison supported a secondary gene flow (Supplementary Table 7). The demographic parameters generated under the optimal model suggest that migration between groups is likely to have occurred between 11.5â17.5 thousand years ago in the two species (Supplementary Table 8), consistent with a possibility that Pleistocene glaciations have driven range shift and secondary contact between previously isolated groups24,25,26. The gene flow, however, likely occurs along the areas where the three groups meet, as the mixed genotypes are mostly found in individuals from the central areas of the Southwest Mountains (Supplementary Fig. 10). Taken together, the migration estimates demonstrate that the central areas of the Southwest Mountains provide potential corridors and routes of connectivity for the three climate-tolerant groups.
Discussion
Understanding speciesâ responses to climate change plays a vital role in developing effective biodiversity conservation plans51. Intraspecific variation in the genotype-climate association has been documented for many species (e.g. refs. 12, 13, 24,25,26) but has only recently been considered in the vulnerability prediction to climate change11,52. By considering intraspecific variation in climate-associated genotypes, we show that the cold-tolerant populations have a greater genomic offset but risk less niche suitability decline under future climate conditions than warm-tolerant populations. By combining genomic offset and niche suitability change, we consider the populations with minor genome and niche interruptions to be the potential donor populations for evolutionary rescue. We then used landscape genetic analysis to identify potential rescue routes to mitigate the challenge of climate change. Altogether, our study demonstrates a framework of incorporating ecologic genomics, niche modelling and landscape genetics in predicting climate change-driven vulnerability and aiding conservation efforts.
The highly heterogeneous environments in mountainous areas often drive intraspecific genetic divergence and local adaptation of mountainous species to specific climatic conditions24,25,26. T. elliotii and P. monticolus are among the most characteristic passerines of the Sino-Himalayan Mountains. We found that the top climatic variables contributing to genomic variation in the two species were related to seasonal temperature and precipitation (Supplementary Table 2). The amplitudes of climatic variables have long been considered to play an important role in shaping the richness and biodiversity pattern of mountainous birds18,19,20. Our results demonstrate that climatic fluctuation is also the definitive factor in shaping the genomic variation and divergence of these mountainous birds, further supporting the role of topological complexity and climatic heterogeneity in shaping mountainous biodiversity19.
By mapping genotype-climate associations across the distribution ranges of the two birds, we identified that the populations in the cold-dry areas showed distinctly different associations from those in the warm-humid and warm-dry areas. Additionally, our cluster analyses of SNPs that are significantly associated with the top climatic variables separate the cold-dry, warm-humid and warm-dry tolerant individuals. The three genetic groups are geographically congruent with the currently recognised ecological zones. For example, the cold-dry tolerant group of T. elliotii inhabits the southern Tibetan zone, which has an average elevation of 4500âm a.s.l. and is dominated by a plateau with a cold-dry climate and typical alpine meadow and shrubland. The warm-dry groups of T. elliotii and P. monticolus are mostly located in the western mountainous plateau zone, which mostly comprises mid-elevational mountains with a typical temperate climate. However, in the spatially heterogeneous southwest mountainous zone, we found that the populations of T. elliotii in the temperate coniferous forests in the northern part and the tropical broadleaf forests in the southern part are clustered into cold-dry and warm-humid tolerant groups, respectively. The cold-dry tolerant populations showed greater genomic offset values than the warm-humid tolerant populations (Supplementary Fig. 11). Altogether, the intraspecific variation in genotype-climate relationships is likely due to local adaptations to heterogeneous climatic conditions, which drive population-specific responses to future climate change.
Based on the intraspecific variation observed in the two species, we predicted population-specific genomic offsets for future climate change. We found that the populations in the cold areas show higher genomic offsets to future climate change than the populations in the warm areas. When incorporating intraspecific variation in predictions of niche suitability under future climate conditions, we found that the niche suitability for the warm-tolerant populations will decrease to a much greater degree than for the cold-tolerant populations. In addition, we observed an expansion of the potential niches for the cold-tolerant populations but not for the warm-tolerant populations (Supplementary Fig. 8). These results suggest that warm-tolerant populations, despite having relatively minor genomic offsets to climate change, may still need to cope with habitat degradation under future climate conditions.
The cold-tolerant populations are mostly found in the mountains in the southwest mountainous zone and eastern Himalayan zone, which are characterised by drastic elevational variation from mountain tops to valleys. For example, a series of mountains, including Gongga Mountain, reach above 7000âm a.s.l., while the Red River valley lies at elevations below 300âm a.s.l.22. Considering that such a wide elevation gradient allows up- and down-slope dispersal, it seems reasonable to expect that cold-tolerant populations can still maintain their niche suitability even through drastic climate fluctuations. However, the warm-tolerant populations in the western mountainous plateau zone and in the southern part of the southwest mountainous zone are mostly found in the patchily distributed mid-elevational mountains. The summits of these mountains are below 4000âm a.s.l., e.g. 3771âm a.s.l. for the Qinlin Mountains (Mt. Taibai, see ref. 53) and 3105âm a.s.l. for the Shenlongjia Mountains (Mt. Shennongding, see ref. 54). As the species studied herein have almost reached their upper elevational limits, there are few available niches allowing uphill migration55. Consequently, these populations with decreasing niche suitability are more likely to be at risk for habitat decline when the climate becomes warmer.
When species are vulnerable to future climate change, their survival may depend on evolutionary rescue. The populations in the central areas of the Sino-Himalayan Mountains not only have the least genomic offset to climate changes but are also predicted to have only minor niche suitability interruption throughout climate fluctuations. Consequently, they can serve as potential donor populations for evolutionary rescue. In fact, this part of the mountains has long been considered to be 'glacial refugia', where populations were maintained throughout the Pleistocene glaciations despite the drastic climatic fluctuations, while populations elsewhere experienced severe bottlenecks24,25,26.
Evolutionary rescue is also dependent on landscape connectivity between populations. The central areas of the Southwest Mountains are characterised by a series of parallel mountains, rivers and valleys running in a south-north direction56,57, which provide corridors for the exchange of individuals living in different ecological zones22,58. Indeed, our landscape connectivity analyses show that dispersal possibilities are possible for southward, westward, eastward and northward migration via the central areas of the Southwest Mountains. However, as the southern, northern and eastern parts of the distribution ranges are predicted to face a substantial decrease in niche suitability, an evolutionary rescue is possible via westward migration to the southern Tibetan zone (T. elliotii) and the eastern Himalayan zone (P. monticolus), where suitable niches will expand westward. Consequently, conservation measures to preserve biodiversity from future climate change should consider not only the speciesâ genomic offsets and niche suitability but also the landscape connectivity and their ability to disperse.
Species in the mountains often live in spatially heterogeneous environments and are thus locally adapted to different climate conditions. This intraspecific genetic variation is an important feature that should be considered when estimating vulnerability to climate change. By considering population-level genetic variation in mountainous species in response to future climate change, we demonstrate that a combination of genomic offset and niche suitability provides unique insight into climate change-driven vulnerability. We highlight the importance of integrating multiple factors, i.e., genomic offset, niche suitability and landscape connectivity, in estimating climate change-driven vulnerability.
In reality, the actual evolutionary responses of species to climate change are more complex than those explained by the aforementioned three factors. It requires understanding interactions between local adaptation, phenotypic plasticity, evolutionary potential, effective population size, dispersal ability, and interspecific interactions (i.e., refs. 15, 51, 59). Consequently, a combination of multiple predictors will improve our understanding of climate change-driven species vulnerability. Further implementing experiments on the physiological tolerance of species, i.e., common garden and transplant experiments60, as well as a functional test of climate-sensitive genes48, will shed more light on how species respond to future climate change.
Methods
Study areas
The distribution ranges of T. elliotii and P. monticolus cover the Sino-Himalayan Mountains (including the Southwest Mountains and East Himalayas) and a chain of mid-elevational mountains in Central China (Fig. 2a). The Southwest Mountains are located at the junction of the Palaearctic and Oriental biogeographic realms21. This region comprises a cluster of high mountains with a high range in elevation and heterogeneous climates. The northern part of the Southwest Mountains, the southern Tibetan zone, is located in the southeastern part of the Qinghai-Tibetan Plateau, which has an average elevation of 4500âm a.s.l. (STZ in Fig. 2a). This zone is dominated by a cold-dry climate and characterised by typical alpine meadow and shrubland. The southern part, the southwest mountainous zone (SMZ in Fig. 2a), has a vertical climatic zonation ranging from subtropical broadleaf forests and temperate coniferous forests to alpine meadows22. On the southwest side of the Southwest Mountains lies the East Himalayan biodiversity hotspot, which is zoo-geographically attributed to the eastern Himalayan zone (EHZ in Fig. 2a). This region is on the southern margin of the Qinghai-Tibetan Plateau and encompasses a broad range of ecological habitats varying from grassy meadows to dense humid evergreen forest. On the east side of the Southwest Mountains, a chain of mid-elevational mountains spans eastward. These mountains and adjacent lowlands belong to the western mountainous plateau zone (WMPZ in Fig. 2a).
Across their distribution ranges, T. elliotii and P. monticulus are commonly found between 1500 and 45,000âm a.s.l. Because distribution maps obtained from the International Union for Conservation of Nature (IUCN, https://www.iucn.org/) include unsuitable habitats, i.e., areas below 1000âm a.s.l., we restricted our study to areas by integrating distribution records and expert distribution maps based on ecological niche modelling61. This restriction removed parts of unsuitable habitats, i.e., a lowland Sichuan Basin with an average elevation of 500âm a.s.l. from our subsequent analyses.
De novo genome assembly and annotation of T. elliotii
We generated a 10X Genomics linked-read library for muscle from a T. elliotii individual collected from Shenlongjia Mountain, Hubei Province, China (voucher number SNJ08157). Sequencing was carried out on BGI-seq 500 platform with PE150. We cleaned the raw reads using SOAPfilter v2.2 with the following steps; (1) removing reads with >10% of N; (2) removing reads with >60% low-quality bases (Phred score â¤10); (3) removing reads with undersize insert size; (4) removing PCR duplicates. All cleaned reads were used to assemble the genome using Supernova v2.0.162 under the 'pseudohap' mode and the intra-scaffold gaps were filled using GapCloser v1.1263. We measured genome completeness using BUSCO v3.0.227. We applied the homologue-based approach to annotate the protein-coding genes by using the protein sequences of Gallus gallus and Taeniopygia guttata. The protein sequences of these reference genes were aligned to T. elliotii genome using TABASTN v2.2.2664 with an e-value cut off 1eâ5, and multiple adjacent hits of the same query were connected by genBlastA v1.0.465. Homologous blocks with lengths larger than 30% of the query protein length were retained. The connected hit region was later extended to include its 2âKb flanking regions, on which gene structure was predicted by Genewise v2.4.166. We then used Muscle v3.8.3167 to align the annotated proteins with the reference proteins. Predicted proteins with a length â¥30 amino acids and identity value â¥40% were retained.
Sample information for resequencing
We sampled 55 and 58 individuals from the distribution ranges of T. elliotii and P. monticolus, respectively (Fig. 2b and Supplementary Data 1). The tissues of these birds were stored in a â80â°C freezer before transportation to the sequencing centre. The Zoological Museum of Institute of Zoology has the authority for specimen collecting and exemption of export/import of samples for scientific purposes (No. 1999/84, provided by Article VII from CITES). Tissue collecting procedures conform to the regulations of the Animal Experimental and Medical Ethics Committee of the Institute of Zoology, Chinese Academy of Sciences.
Genomic data generation and processing
Genomic DNA was extracted from tissues and subsequently sequenced on the Illumina sequencing platform (NovaSeq 6000) at Berry Genomics Corporation (Beijing, China). DNA libraries were constructed according to the manufacturerâs instructions and subsequently sequenced in PE150. We conducted quality control to filter reads, (1) â¥10% unidentified nucleotides (N); (2) >10 nt aligned to the adaptor sequence, allowing â¤10% mismatches; (3) >50% bases having Phred quality <5; and (4) PCR duplicates.
Qualified reads of P. monticolus and T. elliotii were aligned against the great tit genome (GCA_001522545.328) and the de novo genome of T. elliotii generated in this study using BWA v0.7.1768 We classified variants using the HaplotypeCaller function from GATK v3.7.069 and Samtools v1.3.170 with default settings and then intersected two VCFs for each species to obtain the final dataset. We further used VCFtools v0.1.12b71 and BCFtools70 to remove indels and keep only biallelic SNPs with the following filtering criteria: (1) minQ >30, (2) min-DP >7 and max-DP <1000, (3) max-missing countsâ=â5, (4) SNPs at least 5âbp away from INDEL regions.
Genomic offset modelling with gradientForest
We used gradientForest29 to model compositional genetic turnover (i.e., turnover in allele frequencies) using nonlinear functions of climatic gradients. The turnover function transforms multidimensional climatic variables to multidimensional genetic space while selecting and weighting these variables so that they can best summarise genomic variation11,29. We used only SNPs with a minor allele frequency >0.05 because rare alleles tend to yield false positives13. We collected 19 climatic variables at 2.5-min resolution from WorldClim72 (http://www.worldclim.org). We randomly selected 50,000 SNPs because of constraints in computational time. For each SNP, we used 500 regression trees to build a function for each of 19 climatic variables. Only SNPs with R2â>â0 (a measure of the response of individual SNPs to the gradient of a given climatic variable) were considered to be âpredictiveâ loci and were further used in the aggregate turnover function, accounting for the importance of climatic variables and the goodness of fit for each SNP.
The gradientForest analysis provided a weighted R2 value for each of 19 climatic variables to assess its importance (Supplementary Table 2). We selected the top variables by ranking their importance and discarded those that were highly correlated with a variable with a higher weighted R2 value (Pearsonâs râ>â0.7). We then used turnover functions to examine changes in allele frequencies among the top variables and to transform them into genetic importance values. To visualise the resulting multidimensional genetic patterns in geographic and biological space, we used principal component analysis (PCA) to reduce the transformed climatic variables into three PCs11,29.
We used the gradientForest outputs calculated for the current climate as the baselines to predict genomic variation under future climate conditions. The Euclidean distance between the current and predicted values is referred to as the genomic offset11. To explore a range of potential future climate conditions, we used four climate models (CCSM4, MICRO-5, MPI-ESM-LR and CNRM-CM5-2) and four emission scenarios (RCP4.5 and RPC 8.5 for 2050 and 2070 projections, Supplementary Table 3) for a total of sixteen future climate conditions to predict genomic offset. To estimate the spatial regions where the genotype-climate relationship will be most disrupted under future climate, we transformed the climatic variables from each of the 16 climate conditions into genetic importance using the turnover functions as described above. For each grid, we calculated the Euclidean distance between the current and future genomic importance values11,29, which serves as a metric of genomic offset. We averaged the genomic offset values of the four climate models. For each of the future climate emission scenarios and decades, we compared genomic offset values among different climate-tolerant groups using the two-tailed Wilcoxon rank-sum test. A P value of 0.05 after FDR correction for multiple comparisons was considered to be significant. We used two datasets, i.e., R2 positive SNPs identified in the gradientForest (5446 and 7294 SNPs for T. elliotii and P. monticolus) and outlier SNPs identified in a combination of three genotype-climate association approaches (72 and 798 SNPs for T. elliotii and P. monticolus, see below), to check the consistency of the results.
Genomic offset modelling with GDM
Ideally, the gradientForest analysis should use at least four individuals for each locality (e.g. ref. 13). However, our study areas are in highly heterogeneous mountains, which presents a great challenge for obtaining samples. Although our study, to our knowledge, represents the densest bird sampling in this region, we could only include two or three individuals of T. elliotii for three localities (i.e., Rangtang, Mangkang and Wulong). As these localities are on the junction areas of the three identified genetic groups, we wanted to retain them because they would help to define fine-scale population structure and estimate the extent of genetic admixture. Given this small sampling size, we implemented a deep-sequencing strategy because low precision in statistical inference due to the low sample size can be offset by using a large number of SNPs in the analyses73,74.
We then conducted a parallel genomic offset analysis using R package GDM46, a distance-based method (i.e., FST), to estimate the relationship between genetic variation and climatic variables. It has been shown that FST analysis provides reliable estimates even for populations with small sample sizes (i.e., two individuals) given a large number of SNPs (i.e., >1500)74,75. After removing autocorrelated variables (Pearsonâs râ>â0.7), we kept BIO1, BIO2, BIO4, BIO7 and BIO14 for T. elliotii, and BIO1, BIO2, BIO4, BIO5, BIO12 and BIO14 for P. monticolus in the GDM analyses. We calculated the pairwise FST matrix among sampling localities based on 100,000 SNPs that were randomly extracted from genomes using diveRsity package76, and rescaled the FST values to range between 0 and 1. To estimate the relative important of these selected variables, we rescaled the maximum value of the fitted I-Splines between 0 and 1, which is proportional to variable importance. We used 'gdm.transform' to predict and map the pattern of genotype-climate association across the distribution ranges of the two species11.
We predicted the genomic offsets to future climate conditions following the same procedure as gradientFroest analyses. Specifically, we calculated genomic offset as the extent of mismatch between current and expected future genetic variation based on genotype-climate association modelled by GDM analysis, under the sixteen climatic conditions that were obtained from four climatic models and four emission scenarios and decades (Supplementary Table 3). We predicted the genomic offset by 'predict.gdm' function, and then obtained a metric of genomic offset for each of the gridded climate points. The resultant genomic offsets were mapped with ArcGIS 10.1 to show the geographic distribution of population-level variation. We compared the genomic offset values among groups using the two-tailed Wilcoxon rank-sum test and a P value of 0.05 after FDR correction for multiple comparisons was considered to be significant. In addition, we run the separate GDM analyses using the datasets of 72 and 798 outlier SNPs identified for T. elliotii and P. monticolus, respectively, in the three analyses of genotype-climate associations (see below).
Identify SNPs associated with climatic variables
To control for false positives that are often observed in individual genotype-environment programmes77, we used three approaches to identify SNPs that are highly associated with the top climatic variables identified by gradientForest, i.e., LFMM34, RAD and dbRAD35,36. For the LFMM, we ran five separate MCMC runs with a latent factor of Kâ=â3 (see Results). P values from all five runs were combined and adjusted for multiple tests using a FDR correction of Pâ<â0.05. RAD and dbRAD analyses were carried out using the R package vegan v2.5-7 (https://CRAN.R-project.org/package=vegan). Only the constrained axes with a significance of Pâ<â0.05 were retained and subsequently used to build the loading value ('species scores' in vegan). SNPs were identified as outliers if their loading values were greater than three standard deviations of the average loading values. We regarded SNPs identified by all three methods as climate-associated SNPs. We annotated these SNPs using SnpEff78 and the resultant genes were enriched using Kobas v3.079 with Benjamini and Hochberg adjusted P values.
Subsequently, we used these climate-associated SNPs to define climate-tolerant groups using a model-based clustering algorithm in Admixture v1.345. We run the programme with tenfold cross-validation and 10,000 bootstrapping replicates for the coancestry cluster (K) between 1 and 6. The optimal cluster was selected as the one with the smallest cross-validation error45. The genetic ancestry of the individual is assigned to the group identified in the optimal genetic cluster (Kâ=â3, see Results) if it has a proportion of inferred ancestry >60%.
Ecological niche modelling predicted niche suitability under future climate
We used an ensemble modelling approach in the R package Biomod247 to model niche suitability under current and future climates. After initially including six commonly used models, i.e., maximum entropy, generalised boosted model, generalised additive model, multivariate adaptive regression splines, classification tree analysis and random forest, we combined the first four models due to the poor performance of the last two models. This ensemble framework was applied to capture the variation in different ecological niche model algorithms and enhance the robustness of prediction80,81.
We collected distribution records from museum collections and the Global Biodiversity Information Facility (GBIF, https://www.gbif.org). We first tested how reducing sampling bias affected autocorrelation using a range of thresholds to thin distribution records, i.e., from 5 to 35âkm with 5-km intervals. We used Moranâs I to estimate the degree of spatial autocorrelation in each thinning threshold. Compared to that of the non-thinned records, we found that Moranâs I decreased and reached a stable level at the approximately 10-km and 20-km thresholds (Supplementary Fig. 12). We decided to use a 10-km threshold to minimise sampling bias because this resolution is considered to essentially capture the environmental variation in the mountains and has been used in previous ecological niche modelling in mountainous areas (e.g., 82,83,84). We removed redundant records and kept only one record every ten kilometres using the function âthin_bâ in the R package ENMwizard85. We run separate models for individuals that were identified as tolerant to cold-dry, warm-dry and warm-humid conditions to evaluate whether they would be affected differentially by climate change. Because some individuals of T. elliotii located in the adjoined areas of the three ecological zones, i.e., Mangkang, Hanyuan and Wulong, showed admixed genotypes, we removed these records from modelling. After filtering, a total of 292 records were retained for T. elliotii and 619 for P. monticolus. These records included 175 and 284 cold-dry tolerant individuals, 49 and 155 warm-humid tolerant individuals and 68 and 180 warm-dry tolerant individuals for T. elliotii and P. monticolus, respectively.
We used 19 climatic variables at a 2.5-min resolution from WorldClim72 to model the current niche suitability. After removing autocorrelated variables (Pearsonâs râ>â0.7) using usdm86, we kept the remaining uncorrelated variables for subsequent analyses (Supplementary Table 9). We used ENMeval87 wrapped in ENMwizard85 with different feature classes (linear, quadratic and hinge) and regularisation multipliers (RM, from 0.5 to 5 with an increment of 0.5) and selected the best combinations for the maximum entropy models85. We used clamping to avoid extreme predictions for climatic values falling outside the range88. To enhance the foresting accuracy, we removed distribution records with the 10% lowest probability detected by ENMeval87. We estimated the best fitting models based on the lowest AIC89.
For the four ecological niche models, we randomly generated 10,000 pseudo-absence points (or background points90) across each groupâs range and gave equal weights to presence data and background points (i.e., 50% balancing the weights of presence and background points to a prevalence of 0.5)91,92. We employed cross-validation with five repeats by randomly splitting distribution records into two subsets; 70% of the data were used to calibrate the models, and the remaining 30% were used for testing. To increase prediction accuracy, we excluded the models with AUC <0.8 or TSS <0.6 from the final ensemble prediction92. We assigned weights to these models based on their TSS values and constructed ensemble models by calculating the weighted mean of niche suitability across the predictions92. When ecological niche modelling is projected outside the range of the climatic variables on which models were calibrated, there are usually nonanalogous climates (i.e., areas where the value of at least one predictor variable is outside the training region)93. To minimise such uncertainties, we made conservative predictions and restricted our projections to those analogous climates that can be sampled in the current distribution ranges.
We projected niche suitability under future climate conditions using four climate models (CCSM4, MIROC5, MPI-ESM-LR and CNRM-CM5) under four emission scenarios (RCP4.5 2050, RCP8.5 2050, RCP4.5 2070 and RCP8.5 2070). We run separate models for each dataset and then averaged the projections across the four climate models. For each of the four emission scenarios and decades, we calculated the change in niche suitability between the current and future climate (NSCâ=âniche suitability index in the future climateâââniche suitability index in the current climate)48. A negative value indicates decreasing niche suitability, while a positive value shows increasing niche suitability in future climate conditions. We compared NSCs between the cold-dry, warm-dry and warm-humid tolerant groups in each species using a two-tailed Wilcoxon rank-sum test and an FDR-adjusted Pâ<â0.05 was considered to be significant.
Combining genomic offset and NSC
As genomic offset and NSC revealed different estimates of the climate change-driven vulnerability among the three groups, we used a genome-niche index, which allows incorporating genomic offset and NSC to better predict the responses of two species subject to future climate change48. Only the niches of increasing future suitability were considered because populations in these areas would not be challenged by niche suitability decline, as we aim to find the populations that are least interrupted by future climate change. The genome-niche index (gni) of each grid is calculated as follows: \({{gni}}_{i}={{nsc}}_{i}^{\alpha }\;{{go}}_{i}^{1-\alpha }\), where \({ns}{c}_{i}\) and \({{go}}_{i}\) are the NSC and genomic offset at location i, respectively, and \(\alpha \in \left[{{{{\mathrm{0,1}}}}}\right]\) is the weight of normalised \({{nsc}}_{i}\), iâ=â1, 2,â¦, n. To minimise the total deviation between \({{ngi}}_{i}\) and \({{nsc}}_{i}\), \({{go}}_{i}\), we used the equation below to find the smallest value of \({{\min }}Y:\)48
We used the artificial bee colony (ABC) algorithm to find the smallest value of min Y and estimate the optimal α value using the R package ABCoptim (https://CRAN.R-project.org/package=ABCoptim). We set the population size to 20, the maximum number of iterations to 1000 and the limit to 50. We normalised \({{nsc}}_{i}\) and \({{go}}_{i}\) values between 0.1 and 0.9. After the ABC algorithm converged, we obtained optimal estimates of \({{nsc}}_{i}\), i.e., under the projection RCP8.5 2050, \(\alpha\)â=â0.52725 for T. elliotii and \(\alpha\)â=â0.617253 for P. monticolus, which were further calculated to obtain \({{gni}}_{i}={{nsc}}_{i}^{0.52725}\;{{go}}_{i}^{0.47275}\) for T. elliotii and \({{gni}}_{i}={{nsc}}_{i}^{0.617253}\;{{go}}_{i}^{0.382747}\) for P. monticolus. We used ArcGIS 10.1 to visualise the resultant genome-niche index across each speciesâ distribution range.
Landscape genetic analysis
To investigate the influence of different landscape features on genetic structure, we employed a modelling approach to test a set of alternative hypotheses of landscape connectivity. We generated a resistance surface for each landscape feature, including elevation, slope, a standard deviation of elevation, habitat suitability and land cover using Circuitscape94. For each elevation, slope and standard deviation of elevation raster layers, we extracted values for all distribution records (used for ecological niche modelling) for each species and calculated the average value for each landscape feature. We used this value as the benchmark (i.e., the lowest resistance cost) and regarded the absolute difference between this value and the focal value of each grid as the resistance cost for this given grid. In addition, we normalised the niche suitability feature by subtracting the habitat suitability index obtained from the ecological niche modelling by 1. We rescaled these continuous values between 1 (the lowest resistance) and 10 (the highest resistance) to normalise these landscape features.
According to our field observations and written sources23, T. elliotii and P. monticolus are mainly found in coniferous forests and broadleaf forests (and shrublands in T. elliotii). We thus assigned the lowest resistance (i.e., 1) to these habitat types. We then assigned a resistance cost of 3 to savannas, 5 to croplands, 7 to barren and sparsely vegetated areas and 10 to water bodies (Supplementary Table 10). Considering that this assignment is arbitrary, we carried out a sensitivity test by assigning different resistance costs to these habitat types. Our model test showed that assigning a low resistance cost to forest habitats always gave better model performance, thus validating our resistance cost assignment in land cover (Supplementary Table 10).
We randomly extracted 100,000 SNPs from the genome to calculate the genetic distance (FST) between populations in the R package diveRsity76. We used linear mixed-effect models to estimate the effects of landscape features on the patterns of gene flow in lme495 to account for multiple memberships by MLPE. Here, the term 'population' refers to localities. In all models, pairwise genetic distance was used as the dependent variable, and resistance costs from the landscape features were used as the independent variables. All independent variables were z-transformed to meet the normality assumption. As a high degree of collinearity between landscape features may bias the parameter estimation, we kept only features that had a Pearsonâs râ<â0.8. We evaluated modelling performance by ranking AIC with second-order bias correction (AICc). The best model was selected by computing Akaike weights (\({w}_{i}\))96. As model selection with AICc can be biased by the nonindependence of pairwise datasets97, we estimated model fit with marginal R2 (i.e., landscape features) and conditional R2 (i.e., landscape features and population effect) using MuMIn98.
Demographic model tests the migration possibility between groups
We used a coalescence-based approach to test the presence of migration between the groups using Fastsimcoal v2.650. We first generated a two-dimensional, folded site frequency spectrum (SFS) dataset based on whole-genome-wide SNPs for each pair of groups using easySFS.py (https://github.com/isaacovercast/easySFS). We compared three demographic models for each pair of groups, M1) no occurrence of gene flow, M2) occurrence of a secondary gene flow assuming that the groups have been in secondary contact at some time, and M3) occurrence of a continuous gene flow assuming that the groups frequently exchange migrants (Supplementary Fig. 9). For each model, we run 100 replicates with 100,000 coalescent simulations with a minimum of 20 and a maximum of 40 cycles in a conditional maximisation algorithm. We specified a mutation rate of 3.3eâ9 per site per generation following the estimate for passerine birds99. We used the AIC to evaluate which model had the higher likelihood (Supplementary Table 7). We simulated 100 replicates of the SFS from the *_maxL.par file (i.e., the parameter estimates that produced the maximum likelihood) for the best-fit run (minimised difference between maximum estimated likelihood and maximum observed likelihood) of the best-fit model. We run the 100 bootstrap replicates and ranked their maximum observed likelihoods. The demographic parameters from the top 30 runs were used to estimate confidence intervals100.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Data availability
The genome assembly and sequencing data of the de novo sequenced individual of T. elliotii have been deposited to the National Genomics Data centre (https://db.cngb.org/) under BioProject accession number CNP0003256 and NCBI with BioProject accession number PRJNA860040. The resequencing data from P. monticolus and T. elliotii have been deposited in the National Genomics Data centre (https://db.cngb.org/) under the accession number CNP0002314 and CNP0002315, respectively. VCF datasets, location records used in ecological niche modelling, climatic data used in genotype-climate-association, outputs of ecological niche modelling and genomic offset analyses are available in Dryad (https://doi.org/10.5061/dryad.brv15dvb5).
Code availability
Analysis scripts can be found in Dryad (https://doi.org/10.5061/dryad.brv15dvb5).
References
Wiens, J. J. Climate-related local extinctions are already widespread among plant and animal species. PLoS Biol. 14, e2001104 (2016).
Scheffers, B. R. et al. The broad footprint of climate change from genes to biomes to people. Science 354, aaf7671 (2016).
Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37â42 (2003).
Urban, M. C. Accelerating extinction risk from climate change. Science 348, 571â573 (2015).
Wingfield, J. C. et al. Organism-environment interactions in a changing world: a mechanistic approach. J. Ornithol. 152, 279â288 (2011).
Mendoza-Gonzalez, G., Martinez, M. L., Rojas-Soto, O. R., Vazquez, G. & Gallego-Fernandez, J. B. Ecological niche modeling of coastal dune plants and future potential distribution in response to climate change and sea level rise. Glob. Change Biol. 19, 2524â2535 (2013).
Saunders, S. P. et al. Community science validates climate suitability projections from ecological niche modeling. Ecol. Appl. 30, 17 (2020).
Peterson, A. T., Cobos, M. E. & Jimenez-Garcia, D. Major challenges for correlational ecological niche model projections to future climate conditions. Ann. N. Y. Acad. Sci. 1429, 66â77 (2018).
Mays, H. L. et al. Genomic analysis of demographic history and Ecological niche modeling in the endangered Sumatran Rhinoceros Dicerorhinus sumatrensis. Curr. Biol. 28, 70â76 (2018).
Malcolm, R. J., Liu, C., Neilson, P. R., Hansen, L. & Hannah, L. A. Global warming and extinctions of endemic species from biodiversity hotspots. Conserv. Biol. 20, 538â548 (2005).
Fitzpatrick, M. C. & Keller, S. R. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol. Lett. 18, 1â16 (2015).
Gotelli, J. N. & Stanton-Geddes, J. Climate change, genetic markers and species distribution modelling. J. Biogeogr. 42, 1577â1585 (2015).
Ruegg, K. et al. Ecological genomics predicts climate vulnerability in an endangered southwestern songbird. Ecol. Lett. 21, 1085â1096 (2018).
Razgour, O. et al. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc. Natl Acad. Sci. USA 116, 10418â10423 (2019).
Valladares, F. et al. The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol. Lett. 17, 1351â1364 (2014).
Bay, R. A. et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science 359, 83â86 (2018).
Rhone, B. et al. Pearl millet genomic vulnerability to climate change in West Africa highlights the need for regional collaboration. Nat. Commun. 11, 5274 (2020).
Rahbek, C. et al. Building mountain biodiversity: geological and evolutionary processes. Science 365, 1114â1119 (2019).
FjeldsÃ¥, J., Bowie, R. C. K. & Rahbek, C. The role of mountain ranges in the diversification of birds. Annu. Rev. Ecol. Evol. Syst. 43, 249â265 (2012).
Freeman, B. G., Scholer, M. N., Ruiz-Gutierrez, V. & Fitzpatrick, J. W. Climate change causes upslope shifts and mountaintop extirpations in a tropical bird community. Proc. Natl Acad. Sci. USA 115, 11982â11987 (2018).
He, J. K., Lin, S. L., Li, J. T., Yu, J. H. & Jiang, H. S. Evolutionary history of zoogeographical regions surrounding the Tibetan Plateau. Commun. Biol. 3, 9 (2020).
Wu, Y. J. et al. Explaining the species richness of birds along a subtropical elevational gradient in the Hengduan Mountains. J. Biogeogr. 40, 2310â2323 (2013).
del Hoyo, J., Elliott, A., Sargatal, J. & Christie, D. A. Handbook of the Birds of the World (Lynx Edicions, 2013).
Qu, Y. et al. Lineage diversification and historical demography of a montane bird Garrulax elliotii - implications for the Pleistocene evolutionary history of the eastern Himalayas. BMC Evolut. Biol. 11, 174 (2011).
Qu, Y. et al. Long-term isolation and stability explain high genetic diversity in the Eastern Himalaya. Mol. Ecol. 23, 705â720 (2014).
Wang, W. J. et al. Glacial expansion and diversification of an East Asian montane bird, the green-backed tit (Parus monticolus). J. Biogeogr. 40, 1156â1169 (2013).
Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210â3212 (2015).
Laine, V. N. et al. Evolutionary signals of selection on cognition from the great tit genome and methylome. Nat. Commun. 7, 9 (2016).
Ellis, N., Smith, S. J. & Pitcher, C. R. Gradient forests: calculating importance gradients on physical predictors. Ecology 93, 156â168 (2012).
Giorgetta, M. A. Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project phase 5. J. Adv. Model. Earth Syst. 5, 572â597 (2013).
Gent, P. R. et al. The community climate system model version 4. J. Clim. 24, 4973â4991 (2011).
Watanabe, M. et al. Improved climate simulation by MIROC5: mean states, variability, and climate sensitivity. J. Clim. 23, 6312â6335 (2010).
Voldoire, A. et al. The CNRM-CM5.1 global climate model: description and basic evaluation. Clim. Dyn. 40, 2091â2121 (2013).
Frichot, E., Schoville, S. D., Bouchard, G. & Francois, O. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol. Biol. Evol. 30, 1687â1699 (2013).
Forester, B. R., Jones, M. R., Joost, S., Landguth, E. L. & Lasky, J. R. Detecting spatial genetic signatures of local adaptation in heterogeneous landscapes. Mol. Ecololgy 25, 104â120 (2016).
Forester, B. R., Lasky, J. R., Wagner, H. H. & Urban, D. L. Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. Mol. Ecol. 27, 2215â2233 (2018).
Li, C. et al. Two Antarctic penguin genomes reveal insights into their evolutionary history and molecular changes related to the Antarctic environment. Gigascience 3, 27 (2014).
Pirri, F. et al. Selection-driven adaptation to the extreme Antarctic environment in Emperor penguin. Preprint at bioRxiv https://doi.org/10.1101/2021.12.14.471946 (2021).
Wang, L. C. et al. Involvement of the Arabidopsis HIT1/AtVPS53 tethering protein homologuein the acclimation of the plasma membrane to heat stess.J. Exp. Bot. 62, 3609â3620 (2011).
Piñol, R. A. et al. Preoptic BRS3 neurons increase body temperature and heart rate via multiple pathways. Cell Metab. 33, 1389â1403 (2021).
Guilherme, A. et al. Neuronal modulation of brown adipose activity through perturbation of white adipocyte lipogenesis. Mol. Metab. 16, 116â125 (2018).
Zhang, Y., Guo, W., zhang, Y., Zhang, H. & Wu, C. Insights into hypoxic adaptation in Tibetan chicken embryos from comparative proteomics. Comp. Biochem. Physiol. Part D. 31, 100602 (2019).
Pizzagalli, M. D., Bensimon, A. & Superti-Furga, G. A guide to plasma membrane solute carrier proteins. FEBS J. 288, 2784â2835 (2020).
Qu, Y. et al. Rapid phenotypic evolution with shallow genomic differentiation during early stages of high elevation adaptation in Eurasian Tree Sparrows. Natl Sci. Rev. 7, 113â127 (2020).
Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655â1664 (2009).
Ferrier, S., Manion, G., Elith, J. & Richardson, K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Diversity Distrib. 13, 252â264 (2007).
Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD - a platform for ensemble forecasting of species distributions. Ecography 32, 369â373 (2009).
Chen, Y. et al. Large-scale genome-wide reveals climate adaptive variability in a cosmopolitan pest. Nat. Commun. 12, 7206 (2021).
Clarke, R. T., Rothery, P. & Raybould, A. F. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. J. Agric. Biol. Environ. Stat. 7, 361â372 (2002).
Excoffier, L., Dupanloup, I., Huerta-Sanchez, E., Sousa, V. C. & Foll, M. Robust demographic inference from genomic and SNP data. PLoS Genet. 9, e1003905 (2013).
Foden, W. B. et al. Climate change vulnerability assessment of species. WIREs Clim. Change 10, e551 (2019).
Smith, T. B. et al. Genomic vulnerablity and soci-economic threats under climate change in an African rainforest bird. Evolut. Appl. 14, 1239â1247 (2021).
Liu, B., Liang, E. Y., Liu, K. & Camarero, J. J. Species- and elevation-dependent growth responses to climate warming of mountain forests in the Qinling Mountains, central China. Forests 9, 11 (2018).
Dang, H. S., Zhang, Y. J., Zhang, K. R., Jiang, M. X. & Zhang, Q. F. Climate-growth relationships of subalpine fir (Abies fargesii) across the altitudinal range in the Shennongjia Mountains, central China. Clim. Change 117, 903â917 (2013).
Lingua, E., Cherubini, P., Motta, R. & Nola, P. Spatial structure along an altitudinal gradient in the Italian central Alps suggests competition and facilitation among coniferous species. J. Veg. Sci. 19, 425â436 (2008).
Zhang, D. C., Zhang, Y. H., Boufford, D. E. & Sun, H. Elevational patterns of species richness and endemism for some important taxa in the Hengduan Mountains, southwestern China. Biodivers. Conserv. 18, 699â716 (2009).
Zhang, R. Z., Zheng, D., Yang, Q. Y. & Liu, Y. H. Physical Geography of Hengduan Mountains (Science Press, 1997).
Liu, Y. et al. Sino-Himalayan mountains act as cradles of diversity and immigration centres in the diversification of parrotbills (Paradoxornithidae). J. Biogeogr. 43, 1488â1501 (2016).
Bush, A. et al. Incorporating evolutionary adaptation in species distribution modeling reduces projected vulnerability to climate change. Ecol. Lett. 17, 1468â148 (2016).
Sparks, M. M., Westley, A. A. H., Falke, J. A. & Quinn, T. P. Thermal adaptation and phenotypic plasticity in a warming world: insights from common garden experiments on Alaskan sockeye salmon. Glob. Change Biol. 23, 5203â5217 (2017).
Merow, C., Wilson, A. M. & Jetz, W. Integrating occurrence data and expert maps for improved species range predictions. Glob. Ecol. Biogeogr. 26, 243â258 (2017).
Weisenfeld, N. I., Kumar, V., Shah, P., Church, D. M. & Jaffe, D. B. Direct determination of diploid genome sequences. Genome Res. 27, 757â767 (2017).
Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience 1, 18 (2012).
Altschul, S. F. et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389â3402 (1997).
She, R., Chu, J. S. C., Wang, K., Pei, J. & Chen, N. GenBlastA: enabling BLAST to identify homologous gene sequences. Genome Res. 19, 143â149 (2009).
Birney, E., Clamp, M. & Durbin, R. GeneWise and genomewise. Genome Res. 14, 988â995 (2004).
Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792â1797 (2004).
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754â1760 (2009).
McKenna, A. et al. The genome analysis toolkit: a mapreduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297â1303 (2010).
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078â2079 (2009).
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156â2158 (2011).
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965â1978 (2005).
Robinson, J. D., Bunnefeld, L., Hearn, J., Stone, G. N. & Hickerson, M. J. ABC inference of multi-population divergence with admixture from unphased population genomic data. Mol. Ecol. 23, 4458â4471 (2014).
Nazareno, A. G., Bemmels, J. B., Dick, C. W. & Lohmann, L. G. Minimum sample sizes for population genomics: an empirical study from an Amazonian plant species. Mol. Ecol. Resour. 17, 1136â1147 (2017).
Willing, E. M., Dreyer, C. & van Oosterhout, C. Estimates of genetic differentiation measured by FST do not necessary require large sample size when using many SNP markers. PLoS One 7, e2649 (2012).
Keenan, K., Mcginnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. diveRsity: an Rpackage for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol. Evol. 4, 782â788 (2013).
Rellstab, C., Gugerli, F., Eckert, I. A., Hancock, M. A. & Holderegger, R. A practical guide to environmental assocaition analysis in landscape genomics. Mol. Ecol. 24, 4348â4370 (2015).
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80â92 (2012).
Xie, C. et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316âW322 (2011).
Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42â47 (2007).
Araújo, M. B. et al. Standards for distribution models in biodiversity assessments. Sci. Adv. 5, eaat4858 (2019).
Boria, R. A., Olson, L. E., Goodman, S. M. & Anderson, R. P. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecol. Model. 275, 73â77 (2014).
Anderson, R. P. & Raza, A. The effect of the extent of the study region on GISmodels of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. J. Biogeogr. 37, 1378â1393 (2010).
Pearson, R. G., Raxworthy, C., Nakamura, M. & Peterson, A. T. Predicting species distributions from small numbers of occurrence records: a test case using crypticgeckos in Madagascar. J. Biogeogr. 34, 102â117 (2007).
Heming, N. M., Dambros, C. & Gutiérrez, E. E. ENMwizard: advanced techniques for Ecological Niche Modeling made easy. https://github.com/HemingNM/ENMwizard (2018).
Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling. Ecography 37, 191â203 (2014).
Muscarella, R. et al. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models. Methods Ecol. Evol. 5, 1198â1205 (2014).
Owens, H. L. et al. Constraints on interpretation of ecological niche models by limited environmental ranges on calibration areas. Ecol. Model. 263, 10â18 (2013).
Akaike, H. New look at statistical-model identification. IEEE Trans. Autom. Control AC19, 716â723 (1974).
Phillips, S. J. et al. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol. Appl. 19, 181â197 (2009).
Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol. Evol. 3, 327â338 (2012).
Bellard, C. et al. Will climate change promote future invasions? Glob. Change Biol. 19, 3740â3748 (2013).
Elith, J., Kearney, M. & Phillips, S. The art of modelling range-shifting species. Methods Ecol. Evol. 1, 330â342 (2010).
Anantharaman, R., Hall, K., Shah, V. B. & Edelman, A. Circuitscape in Julia: high performance connectivity modelling to support conservation decisions. Proc. JuliaCon Conf. 1, 58 (2020).
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1â48 (2015).
Anderson, D. R. & Burnham, K. P. Avoiding pitfalls when using information-theoretic methods. J. Wildl. Manag. 66, 912â918 (2002).
Van Strien, M. J., Keller, D. & Holderegger, R. A new analytical approach to landscape genetic modelling: least-cost transect analysis and linear mixed models. Mol. Ecol. 21, 4010â4023 (2012).
BartoÅ, K. MuMIn: multi-model inference, R package version 1.9.13 (2013).
Zhang, G. et al. Comparative genomics reveal insights into avian genome evolution and adaptation. Science 346, 1311â1320 (2014).
Roesti, M., Kueng, B., Moser, D. & Berner, D. The genomics of ecological vicariance in threespine stickleback fish. Nat. Commun. 6, 8767 (2015).
Acknowledgements
We sincerely thank Huijie Qiao and Xuan Liu for their facilitating ecological niche modelling and landscape genetic analyses. This research was funded by the Strategic Priority Research Programme of the CAS (XDA20050204 to Y.Q.), the Second Tibetan Plateau Scientific Expedition and Research (STEP) programme (2019QZKK0501 to Y.Q.), and the National Natural Science Foundation of China (NSFC32020103005 to Y.Q. and 3213000355 to F.L.).
Author information
Authors and Affiliations
Contributions
Y.Q. and P.G.P.E. designed research; G.S., X.L. and F.L. collected samples; Y.C., Z.J., P.F., P.G.P.E. and Y.Q. performed research; Y.Q., Y.C. and P.G.P.E. wrote the paper. Z.J., X.L., G.S. and F.L. commented on the paper.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Nature Communications thanks Orly Razgour and the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer review reports are available.
Additional information
Publisherâs note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the articleâs Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the articleâs Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Chen, Y., Jiang, Z., Fan, P. et al. The combination of genomic offset and niche modelling provides insights into climate change-driven vulnerability. Nat Commun 13, 4821 (2022). https://doi.org/10.1038/s41467-022-32546-z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467-022-32546-z
This article is cited by
-
Haplotype-based pangenomes reveal genetic variations and climate adaptations in moso bamboo populations
Nature Communications (2024)