Phylogeographic investigation on the spiny frog Quasipaa shini (Amphibia: Anura: Dicroglossidae) using mitochondrial DNA: cryptic species and species complex

Abstract

The spiny-flanked frog Quasipaa shini is distributed in southcentral China. Here, phylogenetic relationships and diversification of the species were examined using mitochondrial COI gene sequences of 70 individuals from eight populations. Phylogenetic analyses revealed that the species contained two deeply divergent clades (Clades I and II) and five subclades in Clade II. Haplotype network, AMOVA and genetic distance estimations also supported the separations of the six phylogroups. Based on the results, we proposed that at least the Clade I should represent a cryptic species. Neutrality tests indicated that all populations and clades have been probably suffered bottleneck effects or subdivisions.

Keywords:

• Genetic diversity
• mitochondrial DNA
• cryptic species
• species complex
• frog

Introduction

The spiny-flanked frog Quasipaa shini (Ahl 1930) is a member of the family Dicroglossidae (Amphibia: Anura; Frost 2018). It is commonly found in the mountain streams at 510–1500 m elevation in southcentral China, ranging from the northern Guangxi, Guizhou, extreme southern Hunan, to southern Chongong provinces (Fei et al. 2012). In this region, huge mountains (e.g. Nanling Mountains and Dayao Mountains) and deep valleys (e.g. major distributaries of Zhujiang River) were expected to promote diversification, especially in amphibians (Che et al. 2010). Whether the remarkable barriers promoted diversification, even speciation within the mountain frog Q. shini, is worth studying.

Because of global changes especially human activities, such as habitat destructions, excessive captures, and pollution, the population size of the spiny frogs including Q. shini has been declining sharply (Fei et al. 2009, 2012). In view of the serious threats and their high sensitivity to environmental factors, Q. shini was listed as the vulnerable (VU) species in the Red List of China’s Vertebrates (Jiang et al. 2016). However, information, especially about the diversification of the species, still remains scarce, and therefore, it is obligatory to develop strategies to make the status of the species clear for protecting it as soon as possible. Obviously, studies on diversification of the species would be a benefit for the conservation of its endemic germplasm resources.

In view of its rapid rate of evolution and maternal inheritance (Sun et al. 2012), Mitochondrial DNA (mtDNA) markers have been increasingly used for estimating levels of genetic diversity and recognizing cryptic species (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 examine the level of diversification in the frog Q. shini occurring from the mountains in the southcentral China. Our target is to explore: (1) whether there is significant population genetic structure even there are lineages that could be recognized as cryptic species; (2) whether the population genetic structure parallels with geographical structure.

Materials and methods

A total of 70 specimens of Q. shini were collected from eight localities (P1–P8) scattering in the distributional range of the species in southcentral China (Figure 1(A); Table 1). Tissues were collected and stored separately in 95% ethanol after specimen fixation. All specimens were preserved in the Chengdu Institute of Biology, Chinese Academy of Sciences (CIB, CAS). Total genomic DNA was extracted from the tissues using a standard phenol-chloroform extraction procedure (Sambrook et al. 1989). COI primers and PCR protocols were in accordance with Che et al. (2012). PCR products were sent to Shanghai Generay Biotech Co. Ltd. (Shanghai, China) for sequencing in both the directions using the same primers as used in the PCR. The new sequences were deposited in the GenBank with accession numbers MK093217–MK093237.

Figure 1. Sampling localities in this study and phylogenetic relationships of Quasipaa shini. (A) Sampling localities (P1–P8) in this study. See information of populations P1–P8 in Table 1. Size of pie charts denoted beside localities is in proportion to number of samples, and different colors represent different clades/subclades. (B) Maximum likelihood tree based on COI gene sequences of haplotypes. Bootstrap supports of major nodes were denoted near them. Insert: the lower and upper frog photos are the phenotype of Clades I and II, respectively. (C) Haplotype network. The circle size is proportional to number of samples. Each subclade was surrounded by a dotted line. One white dot means one mutation step.

Table 1. Genetic diversity and neutrality tests of Quasipaa shini.

Pop/clade	Locality	Voucher number	Latitude (N)	Longitude (E)	n	nh	H	h	π	Tajima’s D	Fu’s Fs
P1	Xingan, Guangxi	CIBWB101-118	25.7918°	110.4687°	18	2	H1–H7	0.856	0.03721	1.09080	7.890*
P2	Ziyuan, Guagnxi	CIBWB201-214	25.9322°	110.6409°	14	2	H1, H2, H4–H6, H8, H9	0.901	0.05004	2.13340*	6.867*
P3	Longsheng, Guangxi	CIBWB301-305	25.7508°	110.1351°	5	2	H10–H12	0.800	0.02579	1.71338	4.428
P4	Jinxiu, Guangxi	CIBWB401-412	24.1415°	110.1174°	12	3	H1, H6, H8, H13–H15	0.848	0.03683	−0.37443	5.663*
P5	Huanjiang, Guangxi	CIBWB501-504	24.8599°	108.2639°	4	1	H16	0.000	0.00000	–	–
P6	Daoxian, Hunan	CIBWB601-607	25.3775°	111.4223°	7	1	H17, H18	0.476	0.00094	0.55902	0.589
P7	Libo, Guizhou	CIBWB701-706	25.4265°	108.1031°	6	1	H19–H21	0.733	0.00185	0.31063	−0.304
P8	Xingyi, Guizhou	CIBWB801-804	25.0309°	104.9875°	4	1	H5	0.000	0.00000	–	–
I	–	–	–	–	10	3		0.511	0.00322	−0.32944	1.485
II	–	–	–	–	60	18		0.932	0.02330	0.72698	2.138
II-1	–	–	–	–	6	3		0.733	0.00357	0.14908	0.826
II-2	–	–	–	–	3	2		0.667	0.00265	–	–
II-3	–	–	–	–	7	2		0.476	0.00094	0.55902	0.589
II-4	–	–	–	–	10	4		0.778	0.00450	1.13264	0.875
II-5	–	–	–	–	34	7		0.831	0.00435	2.10044*	0.053
Total	–	–	–	–	70	21		0.941	0.04046	0.87888	5.850*

Locality: County, Province; n: number of samples; nh: number of haplotypes; H: haplotype; h: haplotype diversity; π: nucleotide diversity. Significance level: *p < .05.

Sequences were aligned using MAFFT 6 (Katoh et al. 2002) with default settings and manually checked. Haplotypes were recognized using DnaSP 5 (Librado and Rozas 2009). Phylogenetic relationships of the haplotypes were reconstructed using maximum likelihood (ML), as implemented in the program PHYML 3.0 (Guindon et al. 2010). For the phylogenetic analyses, COI sequences of one Q. spinosa, one Q. jiulongensis, one Q. boulengeri and one Nanorana taihangnica were downloaded from GenBank (the accession numbers of them were FJ432700, KF199149, KC686711, and KJ569109, respectively) for use as outgroup according to Che et al. (2010). For ML analyses, the best-fitting nucleotide substitution model was selected under the Corrected Akaike Information Criterion (AICc) implemented in JMODELTEST 2.1.7 (Guindon and Gascuel 2003; Darriba et al. 2012). The optimal nucleotide substitution model (GTR + I) was selected for these analyses. Non-parametric bootstrapping with heuristic searches of 1000 replicates was used to assess confidences of branches in ML trees.

Pairwise uncorrected p-distance between phylogroups were estimated using MEGA (Tamura et al. 2011). To visualize the relationships between haplotypes, a haplotype network was constructed using the maximum parsimony method in TCS 1.21 (Clement et al. 2000). Haplotype diversity (h) and nucleotide diversity (π) were estimated using DnaSP. Genetic signals of departure from neutrality or potential population expansion were estimated for populations and major clades using Tajima’s D (Tajima 1989) and Fu’ Fs (Fu 1997) statistics, estimated in DnaSP. 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).

Results

Alignments resulted in 504 base pair long sequences with 74 variable sites. A total of 21 haplotypes were recognized from 70 individuals of eight Q. shini populations. In ML tree, all Q. shini specimens were strongly clustered into one monophyly (Figure 1(B)). It contained two divergent clades, i.e. Clades I and II (Figure 1(B)). Clade I contained three haplotypes for nine specimens from three populations (i.e. populations P1, P2, and P4). Clade II contained 18 haplotypes for 61 specimens from all eight populations. So the two clades commonly existed in the populations P1, P2, and P4 (Figure 1(A)). Clade II was further clustered into five subclades (II-1∼II-5) but with unresolved relationships (Figure 1(B)). Subclades II-2 and II-5 commonly existed in population P3, and subclades II-1 and II-5 commonly existed in population P4 (Figure 1(A)). Other populations each have one subclade, i.e. P5 and P7 each only having subclade II-4, P6 only having subclade II-3 and P8 only having subclade II-5 (Figure 1(A)). The overall genetic distance on the uncorrected p-distance model was 4%. Genetic distance between Clades I and II was 9.4%, much higher than that within Clade I (0.3%) and that between five subclades in Clade II (from 2.5% to 4.4%), of which, the largest one was that between subclades II-1 and II-5, and the smallest one was that between II-4 and II-5. Genetic distance between populations was from 0.7% to 5.5%.

In the haplotype network, the two clades and five subclades of Clade I were all separated by estimations with 95% confidence intervals (Figure 1(C)). Clade I linked Clade II in 39 steps, and the links between five subclades in Clade II need three to ten steps (Figure 1(C)). The estimated haplotype diversity and nucleotide diversity were shown in Table 1. The total population had high h (0.941). In populations, population P2 had the highest h (0.901), and P1, P3, P4, and P7 also had relatively high h (>0.75), but P5 and P8 had zero h. In phylogroups, Clade II had a high h with 0.932, much higher than h of Clade I with 0.511. Subclades II-1, II-4 and II-5 had high h (>0.73), while II-2 and II-3 had low h (<0.67). The total population had middle π with 0.04046. In populations, populations P1–P4 had middle π > 0.025 but <0.052, and populations P5–P8 had very low π (<0.002). In phylogroups, Clade II had middle π with 0.02330, much higher than π of Clade I with 0.00322. All subclades of Clade I had very low π (<0.0046). When the total population was divided into two groups corresponding to two clades (Clades I and Clade II), AMOVA presented that 75.98% of the molecular variance was attributed to the differentiation among groups, whereas 21.18% of the molecular variance was derived from that within groups. When the total population was divided into six groups corresponding to six subclades (Clades I and five subclades of Clade II), the analyses showed that 91.68% of the molecular variance was attributed to the differentiation among groups, whereas 3.23% and 5.1% of the molecular variance was derived from that within groups and within populations, respectively.

In the total population, Tajima’s D and Fu’s Fs tests all resulted in positive values, though the Fu’s Fs value was significant (p < .05; Table 1), suggesting that there was no recent population expansion. In populations, only P4 had negative but not significant Tajima’s D value (p > .05; Table 1), and hence, all of the populations did not experience recent population expansion (Table 1). In phylogroups, only Clade I had negative but not significant Tajima’s D value (p > .05; Table 1), and hence, all of clades/subclades also did not experience recent population expansion (Table 1).

Discussion

In this study, we explored the level of diversification of the mountain frog Q. shini through sampling populations almost across its distributional range in southcentral China (Figure 1(A); Fei et al. 2012). Based on the mitochondrial COI gene sequences, usually used as the DNA barcoding (e.g. Smith et al. 2005; Borisenko et al. 2008; Tavares et al. 2011; Braukmann et al. 2017), our phylogenetic analyses resulted in two major clades including six subclades (Figure 1(B)), and similarly, the haplotype network divided all haplotypes into six groups in accordance with the six subclades (Figure 1(C)). These results indicated that there might be a cryptic species in the relative widely-distributed species of Q. shini. In general, integrative taxonomy together with morphological and molecular pieces of evidence advocated for recognizing new species with discontinuous or diagnostic morphological characters (e.g. Roxo et al. 2015; Li et al. 2018). Correspondingly, we found at least one exclusive morphological character for identifying Clades I and II: in Clade I, there was only one black spine on each orbicular wart on the chest of male in breeding seasons (Figure 1(B) insert: the lower frog photo); but in Clade II, there were more than two black spines on each orbicular wart on the chest of male (Figure 1(B) insert: the upper frog photo). Because the specimens of Clade II presented almost consistent morphology with the holotype and topotype of Q. shini in the description by Ahl (1930), we recommended Clade I as a different species (Quasipaa sp.).

Although for the time being, we have not found diagnostic morphological characteristics between five subclades of Clade II, it did not conclude that there were no cryptic species in them. If we apply the 2% threshold of interspecific genetic divergence as a limit among species (as proposed by Smith et al. 2005), the five subclades with the genetic distance of more than 2.5% between each pair should be considered as different species. This fitted the viewpoint for using genetic divergence to separate species that stated that sufficient genetic distance indicated reproductive isolation between possible species or speciation (e.g. Paggi et al. 1998; Baric and Sturmbauer 1999; Ferguson 2002). De Queiroz (2007) proposed that different species delimitations were derived from different methods using different properties (molecular or morphological) and/or focusing on different stages of the speciation process. The molecularly distinguished but morphologically undistinguished subclades (possible species) in this study probably indicated that discontinuous morphological differences between genetically divergent species are accumulating but was not enough for us to detect at present. In addition, the ecological differentiation as another factor for speciation should also be investigated in the future study. Finally, in view of the presence of deep genetic divergence and the lack of phenotypic discontinuities among local populations/subclades of the Q. shini Clade II, we suggested that this group should be a species complex similar to other problematic amphibian species complex or groups, such as Odorrana schmackeri species complex (Li et al. 2015) and Amolops mantzorum species group (Lu et al. 2014).

All detected clades/subclades of the Q. shini species complex should represent long-time evolutionary units of biodiversity. And our neutrality tests (Table 1) indicated that all populations of the group have not experienced recent population expansions and might have experienced bottleneck effects or subdivisions (Miracle and Campton 1995). Furthermore, in the past decades, the population of the group has decreased because of habitat destructions and excessive captures (Fei et al. 2009, 2012), and thus, the species was listed as VU (Jiang et al. 2016). Therefore, it is important to protect each independent clade/subclade of the group as soon as possible to avoid losing the potential species diversity.

Acknowledgements

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

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Key Research and Development Program of China under Grant number [2016YFC1200705], The Liaoning Province Natural Science Foundation of China under Grant number [201602677].