Mitochondrial DNA revealed the validation of Quasipaa robertingeri (Amphibia: Anura: Dicroglossidae) and its population genetic diversity

Abstract The spiny frog Quasipaa robertingeri is endemic to a narrow region of southwest China and its taxonomic validation is still controversial. Based on COI gene sequences of 110 individuals from seven populations of Q. robertingeri and its related species, we investigated the phylogenetic position and population genetic structure of the species. Phylogenetic analyses indicated that Q. robertingeri was deeply genetically separated from its closely related species Q. boulengeri, indicating the validation of the species. All samples of Q. robertingeri were clustered into two divergent lineages. Haplotype network, AMOVA, and genetic distance estimations also supported the separations of the two groups. Neutrality tests indicated that one lineage has been likely independently experienced a recent population expansion, leading to a secondary contact area between the two lineages.


Introduction
The spiny frog Quasipaa robertingeri belongs to the family Dicroglossidae (Amphibia, Anura; Frost 2020), and is commonly found in the mountain streams at elevations ca. 400-2000 m in a narrow region, bordering Sichuan and Guizhou provinces and Chongqing City in southwest China (Fei et al. 2009(Fei et al. , 2012. The validation of the species is still on debates. Some molecular phylogenetic studies suggested that it should be synonym with Q. boulengeri, which is broadly sympatric with Q. robertingeri (Che et al. 2009;Zhang et al. 2018). But Fei et al. (2009Fei et al. ( , 2012 still recognized Q. robertingeri as a valid species based on morphology. Whatever, there has been no work with comprehensive sampling on the population level of the species and its related species to resolve the systematic problems. In recent years, a variety of human-caused threats especially excessive captures and habitat destructions lead to the dramatic decline of Q. robertingeri (Jiang et al. 2016). In view of the serious threats and their high sensitivity to environmental factors, this species was listed as the VU (vulnerable) species in the red list of China's vertebrates (Jiang et al. 2016). However, knowledges especially on the population diversification of this species remains scarce, and it is obligatory for us to put strategies to resolve the status of the species. Investigations on population diversification of the VU species would promote conservation of its endemic germplasm resources and genetic diversity.
In light of its rapid rate of evolution and maternal inheritance (Sun et al. 2012), mitochondrial DNA markers have been often used to evaluate genetic diversity, detect phylogenetic relationships, and recognize phylogeographic clusters (Weiss et al. 2011;Liu et al. 2015;Wang et al. 2017). In this study, the mitochondrial cytochrome oxidase subunit I (COI) gene was used to investigate phylogenetic position of Q. robertingeri and its population genetic diversity.

Materials and methods
A total of 110 specimens of Q. robertingeri were collected from seven localities (P1-P7) scattered in the boundary of Sichuan and Guizhou provinces and Chongqing City, China ( Figure 1(A,B); Table 1). In addition, 10 specimens of Q. boulengeri were also collected in six same places together with Q. robertingeri (Figure 1(A)). Tissues were collected and stored separately in 95% ethanol before specimen fixation. Specimens were preserved in the Chengdu Institute of Biology, Chinese Academy of Sciences (CIB, CAS). Total DNA was extracted from the tissues using a standard phenol-chloroform extraction procedure (Sambrook et al. 1989). COI primers and PCR protocols follow Che et al. (2012). Sequencing was conducted in both directions using the same primers as used in the PCR using ABI3730 sequencer in Shanghai Generay Biotech Co. Ltd. (Shanghai, China). The new sequences were deposited in the GenBank with  Sequences were assembled and aligned using BioEdit 7.0.9.0 (Hall 1999) with default settings. Haplotypes were recognized using DnaSP 5 (Librado and Rozas 2009). Phylogenetic relationships of haplotypes of Q. robertingeri and its related species were reconstructed using maximum likelihood (ML), as implemented in the program PHYML 3.0 (Guindon et al. 2010). Nanorana parkeri was used as outgroup according to Che et al. (2009). For ML analyses, the best-fitting nucleotide substitution model was selected under the corrected Akaike information criterion (AICc) using JMODELTEST 2.1.7 (Guindon and Gascuel 2003;Darriba et al. 2012). The optimal nucleotide substitution model (GTR þ IþG) was selected for the analyses. Non-parametric bootstrapping with heuristic searches of 1000 replicates was used to assess confidences of branches in ML trees. Pairwise uncorrected pdistances between clades and populations were estimated using MEGA 6 (Tamura et al. 2013). In addition, a haplotype network of Q. robertingeri was constructed using maximum parsimony method in TCS 1.21 (Clement et al. 2000).
Haplotype diversity (h) and nucleotide diversity (p) were estimated using DnaSP. 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.

Results
Alignment resulted in 504 base pair long sequences. Sixteen haplotypes were recognized from 110 individuals from seven populations of Q. robertingeri (Table 1). Phylogenetic analyses clustered the 16 haplotypes of Q. robertingeri into a clade (Figure 1(C,D)). Ten specimens of Q. boulengeri collected in this study were nested with five Q. boulengeri sequences and two sequences named as 'Q. robertingeri' downloaded from GenBank in a clade (Figure 1(C)). The Q. robertingeri clade, the Q. boulengeri clade, and one Q. verrucospinosa were clustered into a big clade (bootstrap support ¼ 70%) separating from other congeners (Figure 1(C)) though the relationships of the three species in this clade were not resolved. The Q. robertingeri clade contains two divergent lineages, lineages A and B. Lineage A contained nine haplotypes (Haps 1-9), and lineage B contained seven haplotypes (Haps 10-16) ( Figure  1(C)). Lineage A contained samples from populations P1, P2, P4, and P5, and lineage B contained samples from populations P3, P5, P6, and P7 (Figure 1(B)).
The least mean genetic distance between Q. robertingeri and its congeners was 11.5% (vs. Q. boulengeri). Genetic distance between lineages A and B of Q. robertingeri was 3.4%.
Genetic distances between populations fall in the range of 0.0-4.8%, with an overall average at 2.3%.
Haplotype diversity (h) and nucleotide diversity (p) are presented in Table 1. The total haplotype diversity of Q. robertingeri was 0.755. Lineage A had higher haplotype diversity (0.736) than B (0.354). In populations, the population P4 showed the highest haplotype diversity (0.670), followed by the P1 population (0.667) and P2 (0.518), and other populations showed low haplotype diversity (h < 0.5). The total nucleotide diversity was 0.01771, but two major lineages have low nucleotide diversity (p < 0.045). The population P4 showed the highest nucleotide diversity (0.01259) and other populations exhibited low nucleotide diversity (p < 0.0039). In the total population and lineage A, Tajima's D, and Fu's Fs tests all resulted in no-significant positive values (p>.05), and as well, all populations had no significant negative values (p>.05; Table 1), suggesting that there was no recent population expansion in these groups. As note, in lineage B, Tajima's D was significant negative (p<.05; Table 1), indicating a signal for a recent expansion on population size of this lineage.

Discussion
In this study, we explored the phylogenetic status and diversification of the mountain frog Q. robertingeri through sampling populations across its narrow distributional range in southwestern China (Figure 1(C)). Che et al. (2009) andZhang et al. (2018) all based on sequences of a very little number of samples of 'Q. robertingeri' proposed that 'Q. robertingeri' was invalid. However, in our tree based on mitochondrial COI gene sequences showed that the samples recognized as 'Q. robertingeri' in Che et al. (2012) and Zhang et al. (2018) were nested into the Q. boulengeri clade, but the true Q. robertingeri samples was divergent from its closely related species Q. boulengeri and Q. verrucospinosa ( Figure  1(D)). Moreover, the smallest genetic distance on COI gene was 11.5% between Q. robertingeri and Q. boulengeri, being much higher and near four times of 3% as interspecific genetic distance proposed by Zhang et al. (2018) and Vences et al. (2005). Hence, all COI sequences recognized as 'Q. robertingeri' on GenBank were misidentified, and Q. robertingeri should be a valid species. Obviously, Q. robertingeri and Q. boulengeri are sympatric in the distributional range of Q. robertingeri (Figure 1(C)), easily preventing the classifications of the two species also with the reason on the superficially morphological similarity of them. Our result indicated that to evaluate the validity of a species should first take a comprehensive sampling on population level of it and its closely related species.
All 110 samples of Q. robertingeri were clustered into two major divergent lineages (Figure 1(D,E)). They should represent long-time evolutionary units of biodiversity with mean genetic distance of 3.4%. But within each lineage, the divergence between haplotypes was shallow indicating that populations in each lineage were probably young. Neutrality tests suggested that lineage B has been experienced a recent population expansion. Accordingly, the locality for the population P4 located between the distributional ranges of lineages A and B might be the secondary contact area. Based on the results of this study, we at least need to protect each independent lineage of the endemic species as soon as possible to avoid losing important genetic diversity.