Mitochondrial DNA phylogeography of the Guizhou odorous frog: limited population genetic structure and evidence for recent population size expansion

Abstract The Guizhou odorous frog Odorrana kweichowensis is endemic to Guizhou Province, China. In this study, a comparative analysis of the mitochondrial COI and ND2 gene sequences was performed to examine genetic diversity in 109 individuals from ten localities across the geographic range of the species. Haplotype diversity and nucleotide diversity were 0.576 and 0.00055, respectively. Phylogenetic analyses almost nested all haplotypes into one lineage. AMOVA indicated that total variation was mainly derived from variation within individual populations. Neutral tests indicated that a recent expansion occurred in the total population. Fst estimations indicated that genetic divergence was not correlated with geographic distance. Accordingly, the species probably experienced a recent population expansion, and there no obvious population genetic structure is apparent. The findings provide useful information for the conservation of this species.


Introduction
The odorous frog Odorrana kweichowensis Li et al. 2018 (Amphibia, Anura, Ranidae) was recently described to science from mountainous streams in Guizhou Province, China (Li et al. 2018). Although Li et al. (2018) reported three populations of O. kweichowensis in the northern part of Guizhou Province, China, the species was suggested to occur in parapatry with its phylogenetically closest species, O. schmackeri, and to have a wider distributional range in Guizhou Province and perhaps into northern Guiangxi Province, China (Li et al. , 2018Zhu 2016). The karst and mountainous topography (e.g., Dalou Mountain and Wuling Mountain series), including some deep valleys (e.g., Wujiang River and Chishui River), in this region would be expected to impede gene flow between populations of amphibians restricted to mountainous streams (e.g., Che et al. 2010;Yan et al. 2013).
We conducted comprehensive surveys of O. kweichowensis and collected a series of specimens from the northern part to the southern edge of Guizhou Province, China. These surveys showed that a variety of human-caused threats, such as habitat destruction, overharvesting, and pollution, have likely caused a recent decline in the number of O. kweichowensis, as also reported in Li et al. (2018). Based on this information, the species was classified as Vulnerable in the IUCN Red List of Threatened Species (IUCN 2021). Understanding population history and genetic diversity is fundamental for conserving species (reference), but until now there has been limited attention toward O. kweichowensis.
Mitochondrial DNA (mtDNA) markers are frequently used for inferring levels of population genetic divergence and structure (Liu et al. 2015;Wang et al. 2017) due to their rapid rates of evolution and maternal inheritance (Sun et al. 2012). In this study, the mitochondrial Cytochrome c oxidase subunit I (COI) and NADH dehydrogenase subunit 2 (ND2) genes were used to reveal the population genetic structure and diversification history of O. kweichowensis to provide information that will assist in the conservation of the species.

Materials and methods
A total of 109 specimens were collected from ten localities (P1-P10) in Guizhou Province, China that span the known geographical range of O. kweichowensis (Figure 1(A); Table 1 muscle tissue. Muscle tissue samples were taken and preserved separately in 99% ethanol prior to fixation of the voucher specimen in 10% formalin. Preserved specimens were deposited in the Moutai Institute (voucher numbers in Table 1). Total DNA was extracted using a standard phenolchloroform extraction protocol (Sambrook et al. 1989). Two fragments of the mitochondrial COI and ND2 genes were amplified. For COI, Chmf4 (5 0 -TYTCWACWAAYCAYAAAGAY ATCGG-3 0 ) and Chmr4 (5 0 -ACYTCRGGRTGRCCRAARAATCA-3 0 ) were used following Che et al. (2012) and for ND2, Ile-LND2(ATAGGGAGACTTATAGGGGTTC) and Asn-HDN2 (CTAAGTCATTACGGGATCGAGGCC) were used following . Gene fragments were amplified under the following conditions: an initial denaturing step at 95 C for 4 min; 36 cycles of denaturing at 95 C for 30 s, annealing at 46 C (for COI)/57 C (for ND2) for 40 s and extending at 72 C for 70 s. Sequencing was conducted using an ABI3730 automated DNA sequencer at Shanghai DNA BioTechnologies Co., Ltd. (Shanghai, China).
Sequences were assembled and aligned using the ClustalW module in BioEdit v.7.0.9.0 (Hall 1999) under the default settings. The two gene fragments were translated to amino acid sequences in MEGA 6.06 (Tamura et al. 2013), adjusted for open reading frames, and checked to ensure absence of premature stop codons. No-sequenced fragments were treated as missing data. All haplotype sequences were submitted to GenBank under accession numbers MW067106-MW067120 for COI and MW071122-MW071136 for ND2. The two gene fragments were concatenated into a single fragment for each individual for analyses.
Variable sites, conserved sites, and nucleotide composition were estimated using MEGA 6.06. Haplotype diversity (h) and nucleotide diversity (p) were estimated using DnaSP v.5 (Librado and Rozas 2009). Genetic signals of departure from neutrality or potential population expansion were estimated for populations using Tajima's D (Tajima 1989) and Fu' Fs (Fu 1997) statistics, estimated in DnaSP. Pairwise uncorrected pdistances between populations were estimated using MEGA 6.06. The hierarchical distribution of overall diversity was determined using an Analysis of Molecular Variation (AMOVA), as implemented in Arlequin 3.11 (Excoffier et al. 2005) and F-statistics (Fst) among the populations was estimated in DnaSP. For the phylogenetic analyses based on COI þ ND2 genes, corresponding sequences of one O. schmackeri, one O. tianmuii, one O. huanggangensis, one O. yizhangensis, one O. lungshengensis, and one Rana temporaria were downloaded from GenBank (GenBank accession Nos.: MH193609, KX233866, MH193595, MH193616, MH193608 and MH536744) and aligned to the O. kweichowensis dataset as described above. Phylogenetic relationships were reconstructed using Maximum Likelihood (ML) as implemented in the program PHYML v. 3.0 (Guindon et al. 2010). For ML analyses, we ran JMODELTEST v. 2.1.2 (Darriba et al. 2012) with Akaike information criteria on the alignment, resulting in the best-fitting nucleotide substitution models of GTR þ I for the data used in ML. For the ML analysis, branch supports were drawn from 10,000 nonparametric bootstrap replicates. Finally, haplotype networks were constructed using the maximum parsimony method in TCS v. 1.21 (Clement et al. 2000).
In the total population, Tajima's D values and Fu's Fs values were significantly negative (p < 0.01; Table 1), suggesting a recent population expansion. P2, P4, P6, P9 and P10 showed a negative Tajima's D value with no significance and P8 showed a positive Tajima's D value with no significance. In Fu's Fs tests, P2, P4, P9 and P10 showed negative values with no significance, P6 showed a negative value with significance, and P5 and P8 showed positive values with no significance (Table 1).
Pairwise genetic distances among the populations ranged from 0.1% to 0.2% with an average value of 0.2%. AMOVA showed that only 37.4% of the molecular variance was attributed to differentiation among populations, whereas 62.6% of the molecular variance was derived from within populations. The highest Fst value among the populations was 0.750, and the lowest was 0.000 (Table 2).
In the ML tree (Figure 1(B)), all samples of O. kweichowensis clustered into one clade (bootstrap supports, bs ¼ 100) that was recovered as sister to O. schmackeri. The haplotype network showed that H2 was occupied by 70 samples, H1-H11, H14 and H15 each independently linked to H2, and H12 and H13 linked to H11 (Figure 1(C)).

Discussion
The geographical topology within the distributional range of O. kweichowensis in southwestern China is complex and has been proposed to present mountain or river barriers for gene flow that promote speciation (Che et al. 2010;Yan et al. 2013). Mitochondrial genes have often been used to investigate population genetic structure within frog species in southwestern China and have often revealed considerable divergence between populations (e.g., Wang et al. 2012Wang et al. , 2013Wang et al. , 2017. However, our results based on the mitochondrial COI and ND2 gene sequences found low genetic diversity and did not find obvious population genetic structure in the frog O. kweichowensis. In this study, only 16 haplotypes were found among 109 specimens from 10 populations across the range of the species, and one haplotype (H2) occurred in all populations, indicating that all examined populations of the species might have a recent common origin. The genetic distances between populations was very low (average of 0.2% between samples). Further, AMOVA also suggested that more than 60% molecular variance was attributed to the differentiation within populations rather than between populations. This point was also evidenced by neutrality tests that indicated a significant recent population size expansion in the total population of the species. Obviously, the recent population expansion was mainly derived from the ubiquitous distribution of haplotype H2. However, most populations of the frog have been experienced some minimal divergence, for example, the haplotype network indicated that in P4, haplotypes H12 and H13 were probably derived from H11 rather than the common H2. These results indicate that the frogs of the species have probably experienced shallow population divergences only very recently. Nevertheless, based on current results, we could not deduce the ancestral region or expansion center due to seven populations exhibiting unique haplotypes with low genetic diversity.
In summary, genetic diversity and structure in O. kweichowensis were relatively low, indicating that it might be a young species or has experienced bottleneck effects (Miracle and Campton 1995). Future work with more markers might clarify the cause of this pattern of relatively low genetic diversity and structure. However, the existing genetic diversity within this Guizhou-endemic frog species should be protected in light of threats from human activities. Populations P4, P5, P6 and P9 harbor most of the genetic diversity, and so should be the highest priority for conservation.

Acknowledgments
We would like to express our gratitude to the editors and the anonymous reviewers for their critical comments on the manuscript.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability statement
The data of this study are openly available in figshare at http://doi.org/ 10.6084/m9.figshare.13152869. All haplotype sequences were submitted to GenBank under accession numbers MW067106-MW067120 for COI and MW071122-MW071136 for ND2.