Twelve unnamed specimens including four adult males, two adult females, and six tadpoles were collected from Wuyishan National Park, Fujian Province, China (Table 1, Fig. 1, Suppl. material 1). For comparisons, 39 specimens of the subgenus Rana were collected, i.e., six R. zhengi from Gulin County, Sichuan Province, China; five adult males, two adult females and six tadpoles of R. sangzhiensis from its type locality, Sangzhi County, Hunan Province, China; two adult males, one female, and one larval of R. johnsi from northern Vietnam; two larval of R. johnsi from Jinxiu County, Guangxi Province, China; eight adult males and six tadpoles of R. johnsi from Shiwandashan Mountains, Guangxi Province, China; and one adult male of R. weiningensis from its type locality, Weining City, Guizhou Province, China (Table 1, Fig. 1, Suppl. material 1). In the field, the frogs and tadpoles were euthanized using isoflurane, and the specimens were fixed in 75% ethanol. Muscle tissue samples were taken and preserved separately in 95% ethanol prior to fixation. The specimens were deposited in Chengdu Institute of Biology (CIB), Chinese Academy of Sciences, Nanning Normal University (NNU), and Institute of Ecology and Biological Resources (IEBR), Vietnam (for voucher numbers see Table 1 and Suppl. material 1).

Figure 1. 

Locations for specimens used in this study. 1. the type locality of Rana wuyiensis sp. nov., Wuyi Mountain, Fujian Province, China; 2. the type locality of R. sangzhiensis, Sangzhi County, Hunan Province, China; 3. the type locality of R. zhengi, Hongya County, Sichuan Province, China; 4. another locality of R. zhengi, Gulin County, Sichuan Province, China; 5. the locality for R. johnsi in Caobang Province, Vietnam; 6. the locality for R. johnsi in Jinxiu County, Guangxi Province, China; 7. the locality for R. johnsi in Shiwandashan Mountains, Guangxi Province, China; 8. the type locality of R. weiningensis, Weining County, Guizhou Province, China.

A total of 40 samples collected in this study was used in molecular analyses, encompassing twelve unnamed specimens from Wuyi Mountain, six R. sangzhiensis, six R. zhengi, 15 R. johnsi, and one R. weiningensis (Table 1). Total DNA was extracted using a standard phenol-chloroform extraction protocol (Sambrook et al. 1989). Three mitochondrial genes (16S rRNA, ND2, and Cyt b) and three nuclear DNA markers (Tyr, BDNF, and RAG1) were amplified and sequenced for the samples. Primer sequences used for PCR are presented in Table 2. 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, 40 s at appropriate annealing temperature (Table 2); and extending at 72 °C for 70 s. PCR products were sequenced with primers same as used in PCR. Sequencing was conducted using an ABI3730 automated DNA sequencer. New sequences were deposited in GenBank (Table 1).

For phylogenetic analyses, the corresponding sequences for congeners especially for the topotypes of species in the subgenus Rana were downloaded from GenBank (Table 1), mainly derived from previous studies (Yuan et al. 2016; Wang et al. 2017; Wan et al. 2020). For phylogenetic analyses, corresponding sequences of one Odorrana versabilis (Liu & Hu, 1962) and one Pelophylax nigromaculatus (Hallowell, 1861) were also downloaded (Table 1), and used as outgroups according to Yuan et al. (2016).

Sequences were assembled and aligned using the ClustalW module in BioEdit v.7.0.9.0 (Hall 1999) with default settings. The protein-coding gene (Cytb, ND2, BDNF, RAG1, and Tyr1) sequences were translated to amino acid sequences in MEGA v. 6.0 (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. For phylogenetic analyses based on mitochondrial DNA, the dataset was concatenated with mitochondrial gene sequences. The best partition scheme and the best evolutionary model for each partition were chosen for the phylogenetic analyses using PARTITIONFINDER v. 1.1.1 (Robert et al. 2012). In this analysis, 16S gene and each codon position of protein-coding mitochondrial gene were defined, and Bayesian Inference Criteria was used. As a result, the analysis suggested that the best partition scheme is 16S gene/each codon position of protein-coding gene, and selected GTR + G + I model as the best model for each partition. Phylogenetic analyses were conducted using maximum likelihood (ML) and Bayesian Inference (BI) methods, implemented in PhyML v. 3.0 (Guindon et al. 2010) and MrBayes v. 3.2 (Ronquist et al. 2012), respectively. For the ML tree, branch supports (bs) were drawn from 10,000 nonparametric bootstrap replicates. In BI, two runs each with four Markov chains were simultaneously run for 50 million generations with sampling every 1,000 generations. The first 25% trees were removed as the “burn-in” stage followed by calculations of Bayesian posterior probabilities (bpp), and the 50% majority-rule consensus of the post burn-in trees were sampled at stationarity. In addition, to access the genetic isolation between the undescribed species and its closely related species on nuclear DNA, one haplotype network for each nuclear gene dataset was constructed, using the maximum parsimony method in TCS v. 1.21 (Clement et al. 2000). Genetic distance of uncorrected-p-distance model on 16S gene sequences between the new species and its closely related species were estimated using MEGA.

All six adult specimens of the undescribed species were measured (Suppl. material 1). For comparisons, five adult male specimens of R. sangzhiensis, eleven adult male specimens of R. johnsi, and 22 adult male specimens of R. zhengi used in Jiang et al. (1997) were also measured (Suppl. material 1). The terminology and methods followed Fei et al. (2009). Measurements were taken with a dial caliper to 0.1 mm. Twenty-two morphometric characters of adult specimens were measured:

ED eye diameter (distance from the anterior corner to the posterior corner of the eye);

FIIIL third finger length (distance from base to tip of finger III);

FIIL second finger length (distance from base to tip of finger II);

FIL first finger length (distance from base to tip of finger I);

FIVL fourth finger length (distance from base to tip of finger IV);

FL foot length (distance from tarsus to the tip of fourth toe);

HAL hand length (distance from tip of third digit to proximal edge of inner palmar tubercle);

HDL head length (distance from the tip of the snout to the articulation of jaw);

HDW maximum head width (greatest width between the left and right articulations of jaw);

IND internasal distance (minimum distance between the inner margins of the external nares);

IOD interorbital distance (minimum distance between the inner edges of the upper eyelids);

LAL length of lower arm and hand (distance from the elbow to the distal end of the Finger IV);

LW lower arm width (maximum width of the lower arm);

SL snout length (distance from the tip of the snout to the anterior corner of the eye);

SNT distance between the nasal the posterior edge of the vent;

SVL snout-vent length (distance from the tip of the snout to the posterior edge of the vent);

TFL length of foot and tarsus (distance from the tibiotarsal articulation to the distal end of the Toe IV);

THL thigh length (distance from vent to knee);

TL tibia length (distance from knee to tarsus);

TW maximal tibia width;

TYD maximal tympanum diameter;

UEW upper eyelid width (greatest width of the upper eyelid margins measured perpendicular to the anterior-posterior axis).

To reduce the impact of allometry in adults, the correct value from the ratio of each character to SVL was calculated, and then was log-transformed for the following morphometric analyses. One-way ANOVA tests were conducted to test the significance of differences on morphometric characters between the undescribed species and its closely related species. The significance level was set at 0.05.

The morphological description follows the definition in Fei et al. (2009). Sex was determined by examining the gonads. The description of toe webbing followed Savage (1975). The description of digital pad followed Ohler (1995). Comparison characters of known congeners were obtained from the literature (Stejneger 1898; Liu 1946; Fei et al. 1990, 2005, 2009, 2012; Liu et al. 1993; Ye et al. 1993, 1995; Lu et al. 2007; Shen et al. 2007; Li et al. 2008; Yan et al. 2011; Wang et al. 2017; Zhao et al. 2017; Wan et al. 2020). We also examined a series of specimens of Rana (Suppl. material 1).

ML and BI trees of the mitochondrial DNA dataset presented almost consistent topology (Fig. 2A, B). In mitochondrial DNA trees, all samples of the undescribed species were strongly nested into one clade (all supports = 100; Fig. 2B). The R. johnsi group was strongly supported as a monophyletic group containing all samples of R. johnsi, R. sangzhiensis, R. zhengi, and the undescribed species (all supports = 100; Fig. 2B). The R. johnsi group was clustered into the clade corresponding to the subgenus Rana (Fig. 2A). The R. johnsi group contained two clades. In the first clade, samples of R. sangzhiensis and R. zhengi were nested into a clade (all supports = 100), which was weakly clustered as the sister of the undescribed species clade (bs = 52; bpp = 0.80; Fig. 2B). In the second clade, three R. johnsi samples from Vietnam were clustered into one clade, which was sister to the clade containing samples of R. johnsi from two localities of Guangxi Province, China (Figs 1, 2B; Table 1). In addition, the topotype specimen of R. weiningensis was phylogenetically far from the R. johnsi group, and clustered as the basal clade of the genus Rana (Fig. 2A). Haplotype networks based on three nuclear genes all indicated that the undescribed species did not share haplotype with its closely related species R. johnsi, R. sangzhiensis, and R. zhengi (Fig. 2C–E), further indicating the genetic divergence between the undescribed species and its closely related species. As note, on each gene, samples of R. sangzhiensis and R. zhengi massively shared common haplotypes (Fig. 2C–E), indicating their very shallow genetic divergence. The genetic distance on 16S between all samples of undescribed species is less than 0.2% (range 0.0%-0.2%). The genetic distance between the species and its closely related species were as following: vs. R. johnsi from Vietnam 1.3% (range 1.1%–1.7%), vs. R. johnsi from Guangxi, China 0.8% (range 0.8%–0.9%), vs. R. zhengi 1.0% (range 0.9%-1.1%), and vs. R. sangzhiensis 0.9% (range 0.8%–1.1%), being similar to that between the latter four groups (range 0.8%–1.4%). As note, the genetic distance between R. zhengi and R. sangzhiensis was 0.2% (range 0.0%–0.4%), and that between R. johnsi from Vietnam and R. johnsi from Guangxi, China was 0.5% (range 0.4%–1.1%).

Figure 2. 

Phylogenetic relationships of Rana wuyiensis sp. nov. and its relatives A maximum likelihood (ML) tree reconstructed based on the 16S, ND2 and Cyt b gene sequences B a part of the tree highlighting the relationships of the R. johnsi group. ML bootstrap supports/Bayesian posterior probability was denoted beside each node. Sample 1–44 refer to Table 1 C–E haplotype networks constructed based on sequences of nuclear genes BDNF, RAG1, and Tyr, respectively. Different species of the R. johnsi group were denoted as different colors.

The R. johnsi group is phylogenetically clustered into the subgenus Rana, but this group could be identified from other species of the subgenus Rana by the tip of toes with lateroventral grooves (vs. absent in other species of subgenus Rana). The undescribed species could be assigned to this species group by a series of morphological characters: tip of toes flat with lateroventral grooves; body size medium, SVL 41.4–45.6 mm (42.9 ± 1.9 mm, n = 4) in adult males and 47.6–50.3 mm (n = 2) mm in adult females; dorsolateral fold distinct and thin, extending straight from posterior margin of upper eyelid to above groin; tympanum distinct, oval; tibio-tarsal articulation reaching forward beyond tip of snout when leg starched forward; skin ridges distinctly arranged on the dorsal surface of thighs and tibias; adult males with a pair of internal subgular vocal sacs; breeding males possess creamy white nuptial pad with tiny hoar velvety spines on the dorsal surface of the first finger, divided into three parts.

Figure 5. 

Photos of the holotype specimen CIB WY20200913003 of Rana wuyiensis sp. nov. A dorsal view B ventral view C ventral view of hand D ventral view of foot E ventral view of the toe highlighting the lateroventral grooves (arrows).

Although the R. johnsi group and R. weiningensis both have lateroventral grooves on the tip of toes, the undescribed species in the R. johnsi group could be easily distinguished from R. weiningensis by the following characters: males with internal subgular vocal sacs (vs. absent in the latter); males with lager body size (41.4–45.6 mm, n = 4 vs. 32.8–37.4 mm, n = 3 in the latter); and more developed webbing between toes (webbing on fourth toes reaching tip of toe vs. reaching distal subarticular tubercle in the latter).

In the R. johnsi group, the undescribed species could be identified from its closely related species on morphology. ANOVA tests indicated that on the number of transverse skin ridges on the dorsal surface of thighs and tibias, the undescribed species significantly differs from its closely related species (all p-values < 0.01; Table 3; Fig. 3). The undescribed species has larger number of transverse skin ridges either on thighs (mean 14.0 ± 1.7, range 12–16, n = 6), on tibias (mean 12.5 ± 2.0, range 9–15, n = 6), and totally on the two body parts (mean 26.5 ± 2.7, range 22–29, n = 6) than R. sangzhiensis (on thighs mean 9.7 ± 1.3, range 7–11, n = 7; on tibias mean 10.1 ± 1.1, range 8–11, n = 7; and totally on the two parts mean 19.9 ± 1.8, range 17–22, n = 7), R. zhengi (on thighs mean 10.0 ± 1.8, range 7–15, n = 22; on tibias mean 8.1 ± 1.3, range 6–12, n = 22; and totally on the two parts mean 18.1 ± 2.7, range 15–22, n = 22), R. johnsi from Vietnam (on thighs mean 9.3 ± 2.3, range 8–12, n = 3; on tibias mean 9.0 ± 1.0, range 8–10, n = 3; and totally on the two parts mean 18.3 ± 3.2, range 16–22, n = 3), and R. johnsi from Guangxi, China (on thighs mean 10.3 ± 0.9, range 9–12, n = 9; on tibias mean 8.8 ± 1.8, range 6–12, n = 9; and totally on the two parts mean 19.1 ± 2.0, range 16–22, n = 9).

Figure 3. 

Box-plot showing the difference on the number of transverse skin ridges on the dorsal surface of thighs and tibias between different species. Specimens of different species: the holotype CIB WY20200913003 of Rana wuyiensis sp. nov., the topotype specimen CIB SZ2012061203 of R. sangzhiensis, the topotype specimen CIB 950300 of R. zhengi, the specimen NNU 1910009 of R. johnsi from Shiwandashan, Guangxi Province, China, and the specimen IEBR.A 4848 of R. johnsi from Vietnam. Abbreviation: TSR, transverse skin ridges.

On morphometric characters, the results of One-way ANOVA showed that in male, the undescribed species was significantly different from R. sangzhiensis on SVL, HDL, SNT, IOD, UEW, ED, TYD, LAL, HAL, LW, THL, TL, TW, TFL, and FL (all p-values < 0.05), significantly different from R. zhengi on HDL, HDW, SL, IND, IOD, UEW, TYD, LAL, LW, THL, TL, TW, and FL (all p-values < 0.05), significantly different from R. johnsi from Vietnam on SVL, SNT, IND, ED, TYD, and LAL (all p-values < 0.05), and significantly different from R. johnsi from Guangxi, China on HDL, HDW, SNT, IND, UEW, ED, LAL, LW, FIL, FIIIL, and TW (all p-values < 0.05; Table 4).

In total, molecular phylogenetic analyses and morphological comparisons indicated that our specimens from Wuyi Mountain, Fujian Province, China should be classified into the R. johnsi group, and are significantly divergent from its closely related species. The specimens should represent a new species which is described as following section.