Genomic insights into local adaptation in the Asiatic toad Bufo gargarizans, and its genomic offset to climate warming

Abstract Genomic signatures of local adaptation have been identified in many species but remain sparsely studied in amphibians. Here, we explored genome‐wide divergence within the Asiatic toad, Bufo gargarizans, to study local adaptation and genomic offset (i.e., the mismatch between current and future genotype‐environment relationships) under climate warming scenarios. We obtained high‐quality SNP data for 94 Asiatic toads from 21 populations in China to study spatial patterns of genomic variation, local adaptation, and genomic offset to warming in this wide‐ranging species. Population structure and genetic diversity analysis based on high‐quality SNPs revealed three clusters of B. gargarizans in the western, central‐eastern, and northeastern portions of the species' range in China. Populations generally dispersed along two migration routes, one from the west to the central‐east and one from the central‐east to the northeast. Both genetic diversity and pairwise F ST were climatically correlated, and pairwise F ST was also correlated with geographic distance. Spatial genomic patterns in B. gargarizans were determined by the local environment and geographic distance. Global warming will increase the extirpation risk of B. gargarizans.

lizards , but the trend of genetic diversity under the background of global change still remains a sparsely studied area.
Most studies of local adaptation at the molecular level have focused on identifying candidate genes and characterizing their functions (Tiffin & Ross-Ibarra, 2014). Nevertheless, a few studies have been done to compare current and future genotype-environment mismatches and thereby explore the extent of genomic offset to climate change (Cummins et al., 2019). Genomic offset is by definition the genetic distance between the present and the predicted genomic composition, with a smaller genomic offset translated into a lower risk of population decline (Rellstab et al., 2021). Many singlenucleotide polymorphism (SNP) loci and candidate genes closely related to environment variables have been identified in a diverse array of taxa including several species of anuran amphibians such as Pseudophryne guentheri (Cummins et al., 2019), Rana temporaria (Bonin et al., 2005), and Nanorana parkeri . These SNPs and genes have potential impact on the survival of organisms in nature. In recent years, increasingly more studies have used SNP datasets to visualize current and future genomic landscapes (Fitzpatrick & Keller, 2015;Ingvarsson & Bernhardsson, 2020).
These studies not only provide key insights into how organisms adapt to their environments at the genomic level but also contribute to better coping with the impact of future environmental change from the perspective of biological conservation.
Amphibians have unique habitat requirements and life-history characteristics, making them more vulnerable to climate change compared with other terrestrial vertebrates (Ficetola et al., 2015).
The effects of climate change on individual species ultimately depend on their physiological tolerance and dispersal ability (Cummins et al., 2019;Lawler et al., 2010;Qu & Wiens, 2020).
China exhibits complex geomorphological and climatic characteristics. The altitudinal gradient (low in the east and high in the west) has a substantial effect on the distribution and genetic patterns of organisms (Huang et al., 2013;Qu et al., 2014;Wu et al., 2016). Numerous mountain ranges, deep valleys, and rivers make the topography of China particularly complex and provide opportunities for populations to become geographically isolated.
For example, geological and climatic events have contributed to the current geographic distributions of frogs of the genera Babina, Nanorana, Pelophylax, Quasipaa, Rana, and Rhacophorus (Che et al., 2010;Li et al., 2012Liu et al., 2015;Yang et al., 2017;Zhang et al., 2008). Geographic barrier has also been shown to be a key factor determining spatial patterns of genetic variation in several species of anurans such as Feirana quadranus (Wang et al., 2012), Quasipaa boulengeri (Yan et al., 2013), Odorrana schmackeri , Nanorana pleskei , and Leptobrachium ailaonicum (Zhang et al., 2009). Studies of the genomic offset to climate warming at the population level can provide key insights into how patterns of biodiversity vary in a region in response to climate change; however, such studies have been rarely conducted.
The Asiatic toad (Bufo gargarizans) is common throughout its range in China, part of Russia, Japan, and the Korean Peninsula (Fu et al., 2005;Zhan & Fu, 2011). Its wide distribution spanning nearly 20 degrees of latitude with high environmental heterogeneity makes it an excellent model for studying the effects of environmental and geographical factors on the genetic structure of populations (Fu et al., 2005;Macey et al., 1998;Pan et al., 2018;Tong et al., 2017;Wu & Hu, 2010;Zhan & Fu, 2011). Previous phylogenetic studies of B. gargarizans using mitochondrial and nuclear markers have suggested that vicariance, repeated range expansions and heterogeneous landscapes have played key roles in shaping the genetic structure of populations (Fu et al., 2005;Pan et al., 2018;Tong et al., 2017;Zhan & Fu, 2011). However, quantitative analyses of the contribution of geographical and environmental factors to the formation of genetic patterns in B. gargarizans have not yet been conducted. Such analyses are of importance in predicting the impact of future environmental change on B. gargarizans populations and developing effective protection strategies for this species that has undergone rapid population decline and even local extinction in China.
Over the past decade, genome-wide SNP datasets have been increasingly used as molecular markers to study the genetic pattern and diversity of organisms for the following reasons. First, SNP datasets are cost-effective, allow sampling schemes to cover a broader geographical distribution range, and more reliably evaluate genetic patterns and variation at different spatial scales (Anderson et al., 2010;Maigret et al., 2020). Second, using SNP datasets can reduce the number of individuals required for quantifying differences among sampling locations (Nazareno et al., 2017). Third, SNP datasets are useful in detecting weak spatial genetic patterns, thus allowing better implementation of genetic conservation practices (Maigret et al., 2020). Here, we conducted a genome-wide analysis of genetic divergence within B. gargarizans using landscape genetic approaches. Specifically, we aimed to (1) identify the climatic and geographic variables contributing to genomic divergence and local adaptation in B. gargarizans, (2) determine how climate change will affect future patterns of genetic diversity among populations, and (3) evaluate the genomic offset of B. gargarizans to climate change across its range.

| Sampling and DNA extraction
We obtained muscle tissue samples from 94 adult Asian toads col-  Table S1 for detailed population information) in western China are in close geographic proximity, the altitudinal difference between populations is not smaller than 500 m (Figure 1a). The sampling area is characterized by a high degree of environmental heterogeneity, with annual mean temperature and annual precipitation both decreasing from south to north and annual precipitation also decreasing from east to west. Our sampling localities included but were not limited to those reported in previous studies of the genetic diversity, population genetic structure, and phylogeography of B. gargarizans based on mitochondrial DNA (Fu et al., 2005;Hu et al., 2007), microsatellite (Pan et al., 2018;Wu & Hu, 2010), or nuclear and mitochondrial DNA data (Wen et al., 2015;Zhan & Fu, 2011).
Tissue samples were brought to our laboratory in Nanjing, where we used dBIOZOL Genomic DNA Extraction Reagent (Bioer Technology) to extract total genomic DNA from each sample. The quality of genomic DNA was assessed by 1% agarose gel electrophoresis. DNA templates were stored at −80°C for later use.

| Sequencing, quality control, and SNP-calling
We used the double-digest restriction site-associated DNA sequencing method to develop genome-wide SNP markers (Peterson et al., 2012). Libraries were constructed for each sample per the steps provided by Sangon Biotech Co., Ltd. First, genomic DNA The Illumina sequencing data were evaluated using FastQC 0.11.8 (Kircher, 2012), and trimmed by Trimmonatic 0.39 (Bolger et al., 2014) to obtain accurate and valid data. We filtered the original data by removing sequences with N bases, adapter sequences, low-quality bases from the forward and reverse sequences (Q value < 20), and reads with lengths less than or equal to 35 nt. A sliding window trimming was performed when the average quality was below a threshold of window size (5) and required quality (20). Gb (NCBI with a BioProject Accession number of PRJNA628553).
We used SAMtools 1.11 to convert SAM files to BAM files (Danecek et al., 2021) and Picard 2.8.3 to remove repeated reads using "MarkDuplicates" for all individuals (Koike & Shiga, 2007)

| Population genetic structure and diversity
We used two independent multivariate and model-based methods to evaluate the genetic structure of populations. First, we used ADMIXTURE 1.4, a model-based clustering program, to characterize population structure (Alexander et al., 2009). The ADMIXTURE analysis was conducted with K values ranging from 2 to 10, and the optimal number of groups was the K value with the lowest crossvalidation error. We then used the package Plink 1.9 to perform principal components analysis (PCA) based on the high-quality SNP data and thereby analyzed population genetic structure (Chang et al., 2015). The number of clusters among populations was assessed using the K-means clustering method for ADMIXTURE analysis and the number of PC axes that minimized the Bayesian information criterion score for PCA. We used the package ggplot2 in R (R Development Core Team, 2020) to visualize the PCA results.
The GTRGAMMA nucleotide substitution model was used in RAxML 8.2.9 (Stamatakis, 2014) to construct a maximum likelihood (ML) phylogenetic tree based on high quality SNPs. The pairwise mean population differentiation (F ST ) was calculated using the highquality SNPs in VCFtools. Genetic diversity, including nucleotide diversity (π), average expected heterozygosity (He), average observed heterozygosity (Ho), and Wright's inbreeding coefficient (F IS ), was measured using the program Populations in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). Genetic diversity at the population level was assessed using the 21 sampled populations (Table S1).

| Habitat suitability
We used Maxent to calculate resistance surface parameters to map the habitat suitability of B. gargarizans in our study landscape (Phillips et al., 2006). A total of 356 occurrence data were obtained from the 21 populations in this study, and 335 georeferenced records downloaded from the Global Biodiversity Information Facility (Table S2; GBIF.org, 2020). We extracted the 23 environmental variables for GIS data layers obtained from related websites, including 19 bioclimatic variables, solar radiation (SR), mean annual actual evapotranspiration (AET), and altitude from the Worldclim data website (Fick & Hijmans, 2017), normalized difference vegetation index (NDVI, data source: https://earth data.nasa.gov), and mean annual solar radiation (Wm-2, SR, data source from 1961 to 1990: https://earth data. nasa.gov). These environmental variables were selected based both on their perceived biological relevance for the distribution of this toad species and the low degree of collinearity between variables (Giovannini et al., 2014). We retained environmental variables with correlation coefficients <0.8. After threshold-based pre-selection, eight environmental variables were included in the habitat suitability model, including mean diurnal range (Bio2), maximum temperature of the warmest month (Bio5), minimum temperature of warmest month (Bio6), temperature annual range (Bio7) annual precipitation (Bio12), precipitation of wettest month (Bio13), precipitation seasonality (Bio15), and NDVI based on the clustering relation of Pearson coefficient ( Figure S1).
We used the eight environmental variables and 356 occurrence records as input data to predict the habitat suitability across the range of B. gargarizans. The predictive power of the habitat suitability model was evaluated using cross-validation and receiver operating characteristic analysis (ROC) based on 20% of the occurrence localities (Phillips et al., 2006) and showed a good fit (AUC > 0.863).
The habitat suitability values ranged from 0 to 1, with the highest values indicating the greatest habitat suitability. Figure S2 shows the predicted range of B. gargarizans. We divided landscape into a 500 × 500 m grid, and generated the resistance surface layer by assigning a habitat suitability value ranging from 0 to 1 to each grid cell in ArcMap 10.4 (ESRI). We used the resistance surface layer to calculate the least cost-weighed distance (LCD) and circuit distance (CD).

| Geographic distance
We used three pairwise distance matrices (Euclidean distance, LCD, and CD) to estimate the degree of geographic isolation between populations. Euclidean distance is by definition the straight-line distance between localities without accounting for the effect of landscape characteristics. LCD and CD account for the heterogeneity of landscapes. LCD is the path that minimizes the total cumulative cost between two points through heterogeneous landscapes (Wang et al., 2009). CD is calculated by summarizing the costs of all possible paths between two points (McRae & Beier, 2007). CD accounts for factors such as species distributions and migration rates in heterogeneous landscapes by considering the relative resistance to different landscape features (McRae et al., 2008). We used the "spDists" function in the R package "sp" to calculate the Euclidean distance between populations (Pebesma & Bivand, 2005). We used Pathmatrix

| Environmental distance
We used the same environmental variables (Bio2, Bio5, Bio6, Bio7, Bio12, Bio13, Bio15, and NDVI) in the habitat suitability model to estimate the environmental distance between populations. We used ArcMap 10.4 to extract values for every environmental variable of each population. We used R to calculate environmental distance by determining the absolute differences between populations and thereby produce eight pairwise environmental distance matrices .

| Causes of genetic divergence
Genetic divergence was estimated by pairwise F ST based on highquality SNPs across all populations in VCFtools. Genetic, environmental, and geographic distance matrices were obtained from the above methods. The Mantel test was used to determine the relationships between environmental/geographic and genetic distance matrices in the R package vegan. Spearman correlation coefficients were calculated, and two-sided p-values were calculated after performing 999 randomizations of the genetic distance matrix using the "mantel. test" function in the R package vegan.
We used reciprocal causal modeling (RCM) with partial Mantel tests, generalized dissimilarity modeling (GDM), and linear mixed effects models with maximum likelihood population effects (MLPE) to evaluate the effects of geographic distance and environmental dissimilarity on genetic differentiation. RCM was used to evaluate relationships between genetic distance and geographic distance and between environmental distance and geographic distance matrices while accounting for either environmental or geographic dis- GDM is a nonlinear extension of matrix regression that can model spatial variation in pairwise F ST between localities caused by pairwise geographic and environmental differences. We performed GDM using the pairwise F ST matrix for high quality SNP dataset as the response variable in the R package "gdm" 1.5.0 (Fitzpatrick et al., 2021). The maximum value of the fitted I-splines was adjusted between 0 and 1 to assess the relative importance of different geographic and environmental variables used in the generalized dissimilarity model. MLPE linear mixed effects can quantitatively analyze relationships between distance matrices while controlling for nonindependence among pairwise data (Shirk et al., 2017). We used the "MLPE.
lmm" function within the R package "ResistanceGA" (Peterman, 2018) to perform MLPE linear mixed effects models to evaluate relationships between genetic distance and each of the 11 geographic and environmental variables. AIC scores were used to assess relative model support by setting REML = FALSE (Peterman, 2018;Shirk et al., 2017).

| Genomic offset to future climate change
We used seven bioclimatic variables (Bio2, Bio5, Bio6, Bio7, Bio12,  (Fick & Hijmans, 2017). GDM was used to predict the genetic change (genetic offset sensu Fitzpatrick & Keller, 2015) needed to track a changing climate relative to future climate conditions under the SSP245 scenarios (Fick & Hijmans, 2017). We mapped the genomic offset to visualize the spatial distribution of population-level vulnerability to climate change.

| SNP genotyping, population genetic structure and diversity
We obtained 308,551 high-quality SNPs from the 94 toads filtered by VCFtools. The number of SNPs was evenly distributed on the different chromosomes. The number of SNPs was lowest (5751) on chromosome 11 and highest (62,353) on chromosome 1 ( Figure S3).
A ML tree constructed using this subset of SNPs revealed that toads in our sample could be divided into three genetically distinct clusters: the western, northeastern, and central-eastern clusters ( Figure S4).
Pairwise F ST values differed significantly from zero in most population pairs and ranged from 0.020 to 0.598, with an overall mean of 0.311 (Table S3). The results of the ADMIXTURE analysis were consistent with the topology of the ML tree; for K values from 2 to 10, the optimal K value was 3 for our 21 populations ( Figure S5).     GDM explained 69.6% of the deviance in turnover in genetic composition (p < 0.0001). GDM showed that the best supported variables affecting genetic distances ranked in order of importance were Bio13 (p = 0.040), Euclidean distance (p < 0.0001) and

| Environmental and geographic factors contributing to local adaptation
Bio5 (p = 0.040; Figure 3). Results of MLPE linear mixed effects models showed that the best supported variable affecting genetic distance was the Bio12 matrix, followed by Euclidean distance, Bio5 and Bio13 matrixes (

| Genomic offset to future climate warming
The

F I G U R E 4
Climate change-related genetic offset predicted from generalized dissimilarity modeling based on an average scenario of SSP245 in the 2100's climate. The color scale indicates the magnitude of the mismatch between current and future climatedriven turnover of alleles, the more genomic offset to future climate conditions and higher risks of population decline.

F I G U R E 3
Fitted I-splines of the generalized dissimilarity modelling for the relationship between observed compositional dissimilarity and predicted ecological distance (a), and the environmental correlates of partial ecological distance based on three important variables, maximum temperature of the warmest month (b), Euclidean distance (c) and precipitation of wettest month (d). Black dots indicate F ST values between populations. The black line represents the model fit, and the shaded area is the error bands (±one standard deviation). Each curve indicates the relative importance of an environmental variable in explaining changes in allele frequency, holding all other variables constant.
The shape of each curve shows the rate at which allele frequencies change along the gradient.
central populations that have often been reported in earlier studies ( Figure S6; Fu et al., 2005). Our results also indicated that Asiatic toads expanded both from west to east and from south to north. This dispersal pattern is supported by the topology of the ML tree, which shows that the western and central-southeastern clusters are ancestral to the northeastern cluster. Migration from west to northeast was not detected in this study. The mito-nuclear discordance can be caused by various mechanisms, including nuclear introgression, sexbiased asymmetries and irregular gene flow between populations (Funk & Omland, 2003;Stobie et al., 2018;Toews & Brelsford, 2012).
Taken together, we tend to conclude that the complex topography, nuclear introgression and/or sex-biased asymmetries all can be the potential main causes of mito-nuclear discordance in B. gargarizans.

| Population genetic diversity
The SNP-based F ST values among the populations in this study were higher than microsatellite-based F ST values among nine populations estimated using 111 samples (Wu & Hu, 2010), but lower than mtDNA control region-based F ST values among 33 populations estimated using 90 samples (Table S3; Fu et al., 2005). Differences in F ST between mitochondrial and nuclear genes are common in the animal species, mainly due to their differences in evolutionary rates (Allio et al., 2017).  (Luquet et al., 2015), and B. andrewsi (Guo et al., 2016 (Reed & Frankham, 2003). Thus, populations of B. gargarizans that have experienced high degrees of climatic heterogeneity are expected to have increased survival potential compared with populations occurring in areas with more homogeneous climates.
F IS was inversely related to Bio7. High F IS might result in genetic drift and ultimately local differentiation (Rai & Jain, 1982) and is associated with harsh environments (Cummins et al., 2019). The high F IS of B. gargarizans populations in suitable environments suggests that the benefits of local adaptation, or possibly the costs of dispersal to or seeking mates in other environments, exceed the costs of inbreeding in B. gargarizans. Therefore, both climatic and geographic factors contribute to genetic diversity in the toad. Genetic diversity within species is of fundamental importance for organisms' adaptation to future environmental changes and, consequently, for their long-term survival.

| Geographic and environmental isolation
The outputs of RCM and GDM models were consistent: the maximum temperature of the warmest month and Euclidian distance were the main drivers of the current pattern of genomic variation. Geographic and ecological factors have been shown to contribute to spatial genomic patterns in several species. For example, environmental variability has played a major role in shaping patterns of genetic differentiation in the long-lived subalpine conifers (Pinus cembra; Tóth et al., 2019), terrestrial-breeding frogs (Cummins et al., 2019), rainbowfishes (Melanotaenia duboulayi; Smith et al., 2020), and song sparrows (Melospiza melodia; Mikles et al., 2020).
IBD and IBE accounted for a large proportion of the spatial genomic variation in B. gargarizans, which is consistent with the results of an earlier study (Pan et al., 2018). Environmental and climatic variables can explain 69.6% of the deviance in turnover in genetic compo-  Sexton et al., 2014). In our study, both IBE and IBD were the main drivers of spatial genomic differences, which were associated with adaptive differentiation and vicariance.
Previous studies on several species in East Asia have shown that the genomic patterns are driven by distinct mechanisms. For example, long-term gene flow and population history have played more important roles than Pleistocene climate fluctuations in shaping spatial genetic structure in the wild boar Sus scrofa (Hu et al., 2020). However, the geological history and environmental heterogeneity of subtropical China have been the major drivers of genomic diversification and local adaptation in the rapid racerunner Eremias velox (Liu et al., 2019).
Therefore, quantifying the relative contributions of IBD and IBE to the organism genomic evolution is helpful to determine whether organisms are impacted by local adaptation and to predict the effects of climate change on the distribution and survival of organisms in the future.

| Adaptation to local environmental conditions
Both genomic diversity (e.g., H E and π) and F IS related to climatic variables were significantly correlated with climate, suggesting that the direction and strength of climate-driven selection on populations of B. gargarizans vary throughout its range. Environmental variables have been reported to drive local adaptation in amphibians. For example, altitude drives variation in genomic diversity among Tibetan frog (Nanorana parkeri) populations . The unique climatic conditions of East Asia, including oxygen partial pressure (Niu et al., 2021), UV radiation (Norsang et al., 2011), and temperature (Muir et al., 2014), have favored the evolution of genes mediating local adaptation in many species. In this study, local environment significantly impacted the genetic diversity and shaped the genetic structure in B. gargarizans. The effects of local temperature and precipitation on amphibians have been demonstrated in many species, such as P. guentheri (Cummins et al., 2019) and R. temporaria (Muir et al., 2014).
Furthermore, temperature and precipitation are the most significant factors that limit the distribution range of amphibians and drive the formation of genetic patterns of amphibians (Lawler et al., 2010;Muir et al., 2014). For example, changes in precipitation and temperature might affect the phenology, abundance, and performance of amphibians (Ficetola & Maiorano, 2016). Therefore, the influence of local environmental variables contributes to local adaptation in B. gargarizans.

| Genomic offset to future climate change
This study and a previous one (Pan et al., 2018)  populations in the East Asian Tertiary relict Euptelea might be more vulnerable to extinction due to climate warming (Cao et al., 2020).
GDM model analysis based on climate data at the end of the century and the genomic variances of B. gargarizans predicted that populations along the Yangtze River and in the western, northern and eastern coastal parts of the species' range might face higher climate change-induced extirpation risk compared with other populations ( Figure 4). Previous studies suggest that climate change (i.e., precipitation and temperature) increases the risk of local extirpation (Cummins et al., 2019;Mi et al., 2022). The slow rate at which physiological tolerances to heat evolve suggests that high temperature will be the main factor restricting the continued survival of species (Qu & Wiens, 2020). Complex landscapes mitigate the effects of higher temperatures to some extent but also result in increased fragmentation of populations in some areas (e.g., western populations). Thus, mitigating carbon emissions and slowing the pace of climate warming are important for protecting many species, including the Asiatic toad.

| CON CLUS ION
Population structure and genetic diversity analysis based on highquality SNPs revealed three clusters of B. gargarizans in the western, central-eastern, and northeastern parts of the species' range in China. Asiatic toads may have experienced an expansion both from west to east and from south to north. Spatial genomic patterns in B. gargarizans were determined by the local environment and geographic distance, as revealed by the fact that both IBE and IBD were the main drivers of spatial genomic differences in B. gargarizans.
Future global warming will increase the risk of extirpation in populations along the Yangtze River and in the western, northern, and eastern coastal parts of the range of B. gargarizans in China.

ACK N OWLED G M ENTS
We thank Chen Chen, Yi-Jing Chen, Xiao-Li Fan, Zi-Yang Gao, Jian-Li Hehuang Pharmaceutical Co., Ltd. Permissions to collect toads used in this study were approved by local authorities.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw sequences obtained in this study have been deposited in the