Genomic evidence for adaptive differentiation among Microhyla fissipes populations: Implications for conservation

Amphibians require both terrestrial and aquatic environments to complete their life cycles. Thus, they are subject to complex selection pressures stemming from different environments, and these selection pressures are likely to vary geographically with variation in temperature and precipitation. Studies of genetic differentiation along geographical clines allow identification of footprints of these selection pressures.

threatened than mammals and birds (Stuart et al., 2004) and, thus, have been receiving increased attention in the midst of the current biodiversity crisis (Adams et al., 2013;Houlahan et al., 2000;Stuart et al., 2004;Wake, 2012;Wake & Vredenburg, 2008). Habitat loss resulting from environmental alterations has been implicated as one of the primary factors affecting amphibian diversity and numbers (Blaustein et al., 1994(Blaustein et al., , 2010Cushman, 2006). For example, natural habitats have been fragmented and degraded through the continuous modification of landscapes by humans (Lindenmayer & Fischer, 2006). This has had negative effects on demographics and genetic characteristics of amphibian populations (Beebee, 2005), including their adaptation to local environments (Johansson et al., 2007).
In recent two decades, genetic data have become increasingly useful in informing amphibian conservation (reviewed in: . The importance of preserving population genetic diversity has long been recognized to influence long-term viability of population as it influences populations' ability to adapt to changes in environmental conditions (Allendorf & Leary, 1986;Chapman et al., 2009;Reed & Frankham, 2003).
Moreover, inbreeding can influence fitness of amphibians by reducing genetic diversity within populations (Allentoft & O'Brien, 2010;Rowe & Beebee, 2003;Rowe et al., 1999). In particular, small effective population sizes (N e ) can lead to loss of diversity and inbreeding, which can also decrease the populations' adaptive potential.
Levels of genetic diversity and patterns of population differentiation in many amphibians have now been successfully assessed using genetic markers, typically microsatellites (Andersen et al., 2004;Beebee, 2005;Burns et al., 2004;Jehle et al., 2005;Yang et al., 2016;Zheng et al., 2021). As such, these studies have typically been limited to a few genetic loci, although studies based on a larger number of markers have started to emerge (e.g., Funk et al., 2018;Guo, Lu, et al., 2016;Hardy et al., 2021;Thörn et al., 2021). Yet, the genetic underpinnings of adaptations in natural populations of amphibians are largely unknown, although quantitative genetic methods have revealed adaptive differentiation in important life history traits (e.g., Berven, 1982;Laugen et al., 2003;Laurila et al., 2002;Palo et al., 2003).
Many amphibian species have become endangered in China due to rapid socio-economic development, population growth, urbanization and climate change. Although Chinese amphibians face the serious threats due to habitat loss, pollution, alien species invasions and overharvesting (Liao et al., 2015;Liu et al., 2018;Xie et al., 2007), little is known about their genetic structuring (but see: Guo, Lu, et al., 2016;Wei et al., 2020;Yang et al., 2016). The ornate chorus frog (Microhyla fissipes) is a case in point. It is a small-sized anuran, and its wide distribution area in China covers two biodiversity hotspots: Indo-Burman and mountains of south central China. The species occurs at altitudes ranging from 0 to 1,400 m and is an ideal model species to study local adaptation across environmental gradients. The species is commonly found on land around rice fields, bogs, ponds and ditches (Fei et al., 2009). It emerges from hibernation at early April and breeding commences in mid-May in mainland China.
In recent years, rapid urbanization following socio-economic development and climate change have reduced their habitats and population sizes (Fei et al., 2009). Hence, from the conservation point of view, it would be useful to know whether this widely distributed species comprises genetically differentiated, locally adapted populations, and whether there is any evidence of reduced genetic diversity in its population isolates.
The aim of this study was to investigate genome-wide genetic variation and differentiation and their possible association with variation in temperature and precipitation among M. fissipes populations to see whether there is genetic evidence for local adaptation.
To this end, the restriction site-associated DNA sequencing (RADseq) approach (Bell et al., 2015;Brelsford, Dufresnes, et al., 2016;Guo, Lu, et al., 2016;Streicher et al., 2014) with the pool-seq strategy (Kofler, Orozco-terWengel, et al., 2011;Schlötterer et al., 2014), was used to identify SNPs from 10 localities spanning 1,398 km long latitudinal, temperature and precipitation gradient (Figure 1). We first used a suite of different outlier tests to identify the putative signatures of directional and/or balancing selection on the SNPs.
We then looked for associations with key environmental factors (viz., temperature and precipitation) using landscape genomic methods.
In addition, we explored the patterns of genome-wide genetic differentiation and variability associated with geography to get insights into the levels of genetic differentiation on different geographic scales. Finally, we evaluated the relative importance of environmental factors and geography in explaining the distribution of genetic variation across M. fissipes populations.

| Specimen sampling
The specimens used in this study were collected with permis-  Table S1 in Appendix S1). All sampling was conducted in breeding seasons of 2014 and 2015 by collecting an arbitrary sample of adults in breeding condition by hand. All animals were sacrificed using the single-pithing method, and a piece of muscle tissue was collected and preserved in 95% ethanol (Lüpold et al., 2017;Yang et al., 2018;Yu et al., 2018). Based on the findings of a previous study, accurate allele frequency estimation should include at least 10 or more individuals (Mita et al., 2013). Hence, 200 specimens, 20 individuals (10 males and 10 females) from each breeding site were used for sequencing in this study.

| DNA extraction, RAD library construction and sequencing
A standard phenol-chloroform method (Sambrook et al., 1989) was used to extract whole genomic DNA from muscle tissues. DNA quality was assessed by visualizing it on 1% agarose gels. To avoid possible unequal representation of individuals within sequencing pools caused by unequal DNA quality, DNA samples exhibiting degradation on agarose gels were re-extracted and re-examined until high quality DNA was obtained. All non-degraded samples were further quantified using both a Qubit ® fluorometer and NanoDrop ® ND-1000 spectrophotometer (Thermo Fisher Scientific). Thereafter, individual samples were diluted to 10 ng/μl and pooled within each of the ten sampling locations. Pooled samples were re-quantified twice using both the fluorometer and spectrophotometer and equalized to 10 ng/μl for constructing of the RAD-seq libraries.
RAD library construction and sequencing were conducted by Gene Denovo Biotechnology Co., Ltd, following the protocol of Guo, Lu, et al. (2016). In brief, restriction enzyme EcoRI was used to fragment DNA followed by adding a P1 adapter with sequencing primer, forward amplification primer and barcode to each of digested DNA pools. Barcoded samples were pooled and randomly sheared.
After this, a P2 adapter was added to the sheared DNA fragments, followed by enrichment of DNA with a P1 adapter by PCR amplification. Finally, DNA fragments at a size range of 300-500 bp were gel-purified, and the barcoded RAD samples were sequenced on the Illumina HiSeq X platform with 150-bp paired-end strategy. All DNA pools were sequenced on a single sequencing lane.

| Data processing, assembly and alignment
We first used an in-house Perl script which allowed for no mismatches between the sequenced barcode and its sequence to de-multiplexed sequences. To reduce the sequencing error rate, we end-trimmed raw reads to a length of 90 bp using a FASTX toolkit (see: http://hanno nlab.cshl.edu/fastx_toolk it/). We then discarded reads containing one or more bases with a Phred quality score below 10 or if more than 5% of the positions had score below 20.
We de novo-assembled the samples from all 10 sampling sites using Velvet 1.2.02 (Zerbino & Birney, 2008). After testing odd kmer lengths from 49 to 73, 100 being the minimum length parameter and 400 being the insert length parameter, we used K-mer of 55 as the optimum. After this, the de novo-assembled sequences from each sampling site were further assembled using CAP3 with default parameters (Huang & Madan, 1999). In the following analyses, we used the singleton and contig sequences based on the assembly with CAP3 as the reference sequences.
We also used BWA 0.6.1 ) to align qualityfiltered reads from all sampling sites to the reference sequences (Xenopus tropicalis UCB_Xtro_10.0; NCBI). The maximum edited distance for each read was two, and the maximum number of alignments for reads paired properly was one. The maximum insert size for a read pair was 500. The maximum insert size was 500 for a read pair. We used SAMtools 0.1.18  to convert the mapping results in SAM format into BAM format and filtered them for a minimum mapping quality of 20. We then removed duplicate reads using Picard Tools v1.133 (http://broad insti tute.github.io/picar d/).
Finally, BAM files were converted into "mpileup" format with maximum 1,000 reads at a position per BAM file using SAMtools 0.1.18 ).

| Estimation of genome-wide genetic variation and differentiation
Nucleotide diversity (Tajima's π), Tajima's D (Tajima, 1989) and population mutation rate (Watterson's Theta, θ W ; Watterson, 1975) were estimated to characterize patterns of genome-wide genetic variation with PoPoolation 1.2.2 (Kofler, Orozco-terWengel, et al., 2011). We firstly defined SNPs across the entire genome based on several stringent criteria, namely, as sequence coverage has a strong effect on the accuracy of allele frequency estimates in sequencing of pooled samples, we sought to improve the accuracy of these parameter estimates using large sliding windows and high sequence coverage . To estimate π and θ W , we required the identified SNPs to possess a minimum minor allele count of six and coverage of 24 to 1,000 reads for each sampling site. Because Tajima's D is sensitive to variation in coverage, it was estimated using only SNPs with a coverage of 48 reads in all sampling sites. A non-overlapping 200-bp sliding window was used to estimate Tajima's π, Tajima's D and θ W across the reference sequences using minimum base Phred quality of 20 for the analysed sites. We calculated F ST values for each pairwise population comparison to estimate of the degree of population differentiation. For these calculations, a minimum count of the minor allele six across all the 10 populations and coverage between 24 and 1,000 reads for each sampling site were analysed using PoPoolation2 .

| Detection of selection footprints
We identified SNPs being likely differentiated as result selection using PoPoolation2 following Guo, Lu, et al. (2016). We employed two independent methods to identify selection. Firstly, we used  (Foll & Gaggiotti, 2008) to estimate the posterior probability that balancing selection/purifying selection (negative alpha coefficient) or directional selection (positive alpha coefficient) had an effect on a given locus. Briefly, we identified the top candidates of selected loci using prior odds of 100 and ran the reversible-jump MCMC chains of 55,000 with a thinning interval of 10, following 20 pilot runs of 5,000 iterations each and burn-in length of 50,000. SNPs identified by both of the two approaches outlined before were considered as truly differentiated loci. Finally, we also performed the BayeScan analyses within each of the population clusters separately.

| Detecting genetic differentiation associated with environmental parameters
Using all SNPs identified by PoPoolation2, we examined the relationships between genetic differentiation and environmental parameters using a Bayesian approach as implemented in Bayenv (Coop et al., 2010). Bayenv accounts for demographic effects when evaluating the relationship between genetic differentiation and environmental parameters. We first estimated neutral covariance matrix (estimated from markers declared as neutral by BayeScan) and then evaluated the relationship between genetic differentiation and environmental parameters (viz. latitude, average annual temperature and precipitation; Table S1). Annual temperature and precipitation data were retrieved from https://www. meteo blue.com. Each environmental parameter was standardized by subtracting the mean from daily values and dividing them by the standard deviation of the parameter among populations (Coop et al., 2010). To confirm that the results were not sensitive to stochastic errors, we ran three independent runs with different random seeds.

| SNP annotation and gene ontology analysis
To examine genetic differentiation across genome features, the reference sequences were annotated by searching matching proteincoding sequences of Xenopus tropicalis from the Ensembl database, with an E-value cut-off of 1 × 10 -5 using BLAST (Altschul et al., 1997).
Gene ontology (GO) terms of the reference sequences were retrieved based on their orthologs using BioMart (Kasprzyk, 2011). We also performed the GO term enrichment analysis using g: profiler (Uku et al., 2019).

| Detection of population structure
We characterized the population structure on the basis of matrices org/~tryph on/popul ation s/). In both cases, we estimated tree topology and branch-support values using 1,000 bootstrap replicates.

| Relationships among population differentiation, geographic distance and environmental parameters
Isolation by distance (IBD; Wright, 1943) is expected to be observed in most amphibians because of their poor dispersal capabili- (distance) is considered a powerful test of IBD (Rousset, 1997). Therefore, we used the Mantel test to test for IBD for all populations in the "VEGAN" R package (Oksanen et al., 2015). Geographic distances between populations were estimated using Google Earth (Table S2 in Appendix S1). Moreover, we carried out separate analyses to test whether IBD could be observed at a finer geographic scale (≤70 km) than that covering the total sampling area (~ 1,400 km) for each of the population clusters ( Figure 1).
Besides the poor dispersal capabilities resulting in low rates of gene flow among populations, environmental variables (i.e., average annual temperature and precipitation) can potentially act as drivers of population differentiation in M. fissipes. Therefore, to compare pairwise genetic distance matrices and pairwise environmental variables, we performed a partial Mantel test to analyse isolation by environment (IBE; Wang & Bradburd, 2014) using the R package in "VEGAN" (Oksanen et al., 2015). We used partial Mantel tests to rule out false-positive relationships stemming from effects of IBD by repeating IBE tests holding geographic distances constant. As the partial Mantel tests are known to have high Type I error rates (Guillot & Rousset, 2013), we also conducted a conditional independence test to identify cause-effect associations between genetic distances, geographic distances and environmental variables (i.e., average annual temperature and precipitation) in the R package "pcalg" (Kalisch et al., 2012). We performed IBE tests for all 10 populations together, as well as for each of the three population clusters separately.
To detect signals of IBD and IBE, we used six different SNP data sets both across all 10 populations, as well as within each of the (f) the neutral SNPs on the basis of the fact that BayeScan did not identify these as outliers, and environmental parameters in Bayenv analyses were not associated with allele frequency variation in these SNPs. Evidence for IBE in the first five data sets could arguably be regarded dubious and even circular; a signal of IBE in the sixth data set should give the most honest and also a conservative test of IBE.

| RAD-seq data set and the assembled reference sequences
We obtained 275.3 million reads, with their numbers ranging from 19.1 to 41.9 million for each population (Table 1). A total of 1,931,696 unique sequences >200 bp were produced by the de novo assembly, and these sequences were used as references for alignment. The mean size of reference sequences was 269 bp (median =273 bp), ranging from 200 to 1,509 bp ( Figure S1 in Appendix S1). A total of 258.1 million reads were aligned to the reference sequences, with the number of aligned reads were ranging from 17.8 to 39.6 million for each population (Table 1). The number of mapped reads was strongly and positively correlated with the number of raw reads within each population (r s ≥0.9636, p < .0001).

| Genome-wide genetic variation
The number of SNPs in a given population identified by PoPoolation ranged from 213,934 to 567,341 (Table 1) (Table 1). Genome-wide Tajima's π ranged from 0 to 0.50, and the genome-wide θ W was between 0 and 0.26 (Figure 2a,b, Figure S2 in Appendix S1). The island populations tended to exhibit less genetic variability than mainland populations (

| Genome-wide genetic differentiation
A total of 20,572 SNPs across the ten populations were identified However, as 15 and 28 of these SNPs were correlated with variation in latitude and altitude, respectively, we excluded them from the further analyses reducing the data set to 317 SNPs.

| Genetic differentiation associated with environmental parameters
The global Bayenv analyses based on the subset of 317 SNPs across 10 populations revealed that allele frequencies in the number of SNPs were strongly correlated with variations in annual average temperature and/or precipitation. Of these, variation in an allele frequency of 69 SNPs was associated with variation in annual average temperature (Figure 2e), and variation in an allele frequency of 248 SNPs was associated with variation in annual average precipitation ( Figure 2f). The variation in allele frequency in seven SNPs was associated with variations in both annual average temperature and precipitation ( Figure 3a).
Of the 69 SNPs with allele frequency variation associated with annual average temperature in global analysis, 54 did so also in analysis restricted to the Sichuan cluster (Table S3 in Appendix S1).
Corresponding figures were 61 in the Guangxi cluster, 60 in the eastern Hainan cluster and 65 in the western Hainan cluster ( Table S3 in Appendix S1). Of the 248 SNPs globally related to precipitation, 206 did so also in analysis restricted to the Sichuan cluster (Table S3 in Appendix S1). Corresponding figures were 235 in the Guangxi cluster, 229 in the eastern Hainan cluster and 234 in the western Hainan cluster ( Table S3 in Appendix S1).
The proportion of SNPs related to annual average temperature or precipitation was significantly (χ 2 -test: p < .05) higher among the 360 outlier SNPs under diversifying selection than that among all 20,572 SNPs. Of the 360 outlier SNPs, variation in allele frequency in three (0.83%, out of 360) was significantly correlated with temperature, but none with precipitation. Further, the difference in corresponding frequencies among all SNPs was statistically significant were associated with precipitation (χ 2 -test: p < .001).

| Candidate genes for local adaptation
A BLAST search against the genome of X. tropicalis found 23 candidate genes for local adaptation using the 360 SNPs identified as outliers (Figure 3b; Table S4 in Appendix S1). These genes were classified under 60 GO terms and were markedly enriched in functional categories "biological processes" and "molecular functions" (Benjamini-Hochberg FDR, p < .05) as compared with all genes in X. tropicalis (Figure 3b). TA B L E 1 Summary statistics of RAD data, including the number of SNPs identified and their average heterozygosity in each population

| Population structure
Average pairwise F ST among the 10 populations based on the 20,572 SNPs was 0.219, ranging from 0.136 (between BH and NB) to 0.443 (between GH and XD; Table S5 in Appendix S1). The principal component analysis (PCoA) plot revealed clear geographically ordered population structuring (Figure 2g). The first PCA-axis differentiated

| Isolation by distance (IBD) and environment (IBE)
No IBD was detected in analysis of all populations using all 20,572 SNPs (Mantel test, r s = −0.101, p = .514). The IBD was detected  Table 3). However, we found that the IBD was detectable across all populations in the three SNP data sets (viz., the 360 SNPs under diversifying selection, SNPs associated with latitude and temperature; Table 3). IBD was detectable in two SNP data sets across local populations from the Sichuan cluster (SNPs associated with latitude and precipitation) and the western Hainan cluster (SNPs associated with precipitation and neutral SNPs; Table 3). For the Hainan population cluster, we only detected IBD in one of the SNP data sets (SNPs associated with precipitation; Table 3).
Population differentiation was associated with environmental temperature across all populations if all 20,572 SNPs were used, but this was not the case with environmental precipitation (Table 4).
Population differentiation was significantly associated with temperature within the Hainan population cluster but not in others ( Table 4).
The relationship between population differentiation and environmental temperature tended to approach significance in the other five different SNP data sets too (Mantel tests, r s ≥0.118, p < .062; Table 4). The relationship between population differentiation and environmental precipitation was significant in the SNPs identified to be associated with latitude (Bayenv analysis: partial Mantel test, r s =0.292, p = .001; Table 4). Geographic distance displayed a sig-

| D ISCUSS I ON
We found that allele frequency variations in numerous SNPs were significantly associated with variation in average annual temperature and precipitation among M. fissipes populations. This suggests that these two environmental factors, or factors in strong association with them, are likely to have been important shaping forces behind the observed population genetic differentiation providing evidence for local adaptation. The findings further indicate moderate levels of genome-wide genetic differentiation across M. fissipes populations both at short (<100 km) and long (>500 km) geographic distances. To this end, our results align with the earlier findings from low coverage population genetic studies, suggesting that anuran amphibians typically exhibit substantial population structuring that often exceeds that observed in other classes of vertebrates (Ward et al., 1992).
Further evidence for adaptive nature of this differentiation was provided by the fact that the patterns of differentiation in a number of outlier loci aligned with variation in environmental temperature and precipitation both in Bayenv and IBE analyses. Such associations are not surprising given that amphibians are ectothermic animals whose biology is strongly affected by moisture and temperature. Note: r s =Spearman rank correlation coefficient, p refers to empirical significance level from 1,000 permutations.

| Genetic differentiation among M. fissipes populations
The patterns of genetic differentiation were heterogeneous across the ten study populations, suggesting a varying but generally fairly likely to result from their limited dispersal ability and high degree of philopatry (Smith & Green, 2005;Ward et al., 1992). One manifestation of limited dispersal ability of amphibians is the generally high IBD among amphibian populations (Brelsford, Rodrigues, et al., 2016;Dodd, 2009;Guo, Lu, et al., 2016). The current study failed to find evidence for IBD at the global scale of sampling. However, we found significant IBD within the local population clusters. The observed discontinuity in IBD at different geographical scales could owe, for instance, to historical discontinuities among different population clusters or to dispersal barriers of natural and/or anthropogenic origin. Whatever the reason, dispersal among local (and global) populations of M. fissipes seems to be limited. Note: The outcomes NaN (not a number) and NA (not applicable) in the partial Mantel tests occurred because variance in some tests was effectively zero, and test statistics could not be estimated.
Given are simple (Mantel test) and tests accounting geographic distance among localities (partial Mantel). r s =Spearman rank correlation coefficient. p refers to empirical significance level obtained from 1,000 permutations.
The IBE analyses revealed some positive associations between genetic and environmental distances (temperature and precipitation) independent of geographic distance, providing further evidence for the role of environmental factors driving patterns of genetic differentiation in certain genomic regions. In our tests, IBE was detected mainly in cases where the analyses were restricted to outlier loci indicated to be associated with variation in temperature and precipitation, corroborating the conjecture that these environmental factors have been drivers of genetic differentiation among M. fissipes populations.

| Genetic variability within M. fissipes populations
The levels of genetic variability within local population of M. fissipes, as assessed by different measures of diversity, were generally moderate to high. However, the island populations tended to exhibit less genetic diversity than the mainland populations. This aligns with the results of numerous previous studies which have found reduced genetic diversity in insular as compared with mainland populations (Bromham & Woolfit, 2004;Frankham, 1997;Garcia-Verdugo et al., 2009;Mason et al., 2011;Yamada & Maki, 2012). Reduced genetic diversity in island populations can be understood based on the fact that they are closed populations subject to stronger effects of inbreeding and genetic drift than that experienced by outbred mainland populations. However, in our case, the differences in genetic diversity between island and mainland populations were not significant.
The moderate levels of genome-wide genetic differentiation among M. fissipes populations, and in particular, the evidence for adaptive differentiation suggests that different populations of this species cannot be considered a homogenous entity. The significant adaptive differentiation among local populations should be taken into consideration by the conservation and management actions pertinent to them. The adaptive differentiation among local populations suggests that possible conservation or management efforts should target local populations, rather than the entire species, and the island populations should be of particular concern as they are likely to be first ones to suffer negative effects of loss of genetic variability and suffer reduced ability to adapt to future environmental changes (Chapman et al., 2009).

| Methodological considerations
It is now widely recognized that many commonly used genome scan approaches are prone to false-positive findings of local adaptation (reviewed in: Hoban et al., 2016). Approaches to reduce the rate false-positive findings include the use of null models incorporating information on the demographic history of the focal populations, as well as the use of relatedness among populations for correcting for a neutral population structure (Hoban et al., 2016). Similarly, access to a reference genome may markedly improve the inference of genetic of local adaptation by reducing errors due to incomplete information McKenna et al., 2010). One can also summarize correlated environmental parameters with a dimensionality reduction analysis such as principal component analysis (Lasky et al., 2012), and correlate allele frequencies with the obtained principal components. Moreover, composite measures of selection can improve one's ability to filter selection signals, but they are not a panacea and univariate statistics they summarize can limit their power (Lotterhos et al., 2017). Unfortunately, most of these approaches are not implementable with data consisting of pooled DNA samples.
Theoretical studies have suggested that the pool-seq strategy can be more effective in SNP discovery and can provide more accurate estimates in allele frequency than individual sequencing (Guo, Lu, et al., 2016;Schlötterer et al., 2014). Nevertheless, other studies have pointed out shortcomings in the pool-seq strategy (Anderson et al., 2014;Lynch et al., 2014;Schlötterer et al., 2014).
These shortcomings are related to the facts that (a) pipetting or errors of DNA quantification may lead to differential representation of individuals and bias allele frequency estimation; (b) it is more difficult to identify improper alignments of short reads based on the reference sequences in pool-seq than in individual sequencing. Despite these possible shortcomings, pool-seq has been recognized as a valid approach to study population adaptive differentiation among animal populations (Corander et al., 2013;Guo et al., 2015;Guo, Li, et al., 2016;Guo, Lu, et al., 2016;Rellstab et al., 2013). Furthermore, by following the fairly stringent protocol of Schlötterer et al. (2014), we believe that our results and inferences should be reliable and not subject to biases stemming from the issues listed before.

| CON CLUS ION
By utilizing 20,572 SNP loci, we explored genome-wide genetic variability and differentiation within and among M. fissipes populations covering a large geographic area. The overall degree of genetic differentiation was moderate to high. Numerous SNPs were found to be in strong association with variation in average annual temperature and/or precipitation across the sampled locations, suggesting that ecological factors play an important role in driving genetic differentiation among local M. fissipes populations. In general, our results provide strong genomic support for the view that adaptive differentiation among amphibian populations can be considerable, and we suggest that this heterogeneity should be factored into planning conservation and management of amphibian species and populations. We further note that although the Indo-Burma and the mountains of south central China have been already recognized as important biodiversity hotspots, our results show that they may contain yet unrecognized genetic biodiversity as evidenced by intraspecific genetic heterogeneity uncovered in this study.

ACK N OWLED G EM ENTS
We

CO N FLI C T O F I NTE R E S T
The authors claim no conflict of interests.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13433.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences are available from GenBank (SRR16095725-SRR16095734).