We collected 12 Kurixalus samples in April and July 2018 in north-western Zhejiang Province, People’s Republic of China (Fig. 1; 31.06°N, 119.85°E). Specimens in April were photographed and had buccal swabs taken for preliminary analysis and subsequently released (Fig. 2), with a follow-up collection expedition taking place in July, pending positive preliminary results for a potential novel species.

Figure 1. 

Map of sampling site and Kurixalus species. The sample for Kurixalus sp. nov. were collected in April and July 2018 in north-western Zhejiang Province, People’s Republic of China. Map generated in ArcMap 10.4.

Figure 2. 

A Kurixalus sp. nov. specimen (HM 323117) in-situ from 26 April 2018 and B dorsal and C ventral view of NJFU20180704005.

Specimens collected in July were humanely euthanised through cooling in line with Shine et al. (2015) and a subsequent application of 20% benzocaine applied to the venter (Torreilles et al. 2009). Specimens were deposited at the Biological Museum at Nanjing Forestry University (institutional code NJFU; Table 1). Genomic materials were collected from buccal swabs for the initial three individuals found in April (photographic vouchers deposited in the repository institution HerpMapper.org (institutional code HM) (HerpMapper 2020): HM 244044, HM 323117 and HM 323118) and thigh muscle tissues for the 11 specimens collected subsequently (Table 1). Genomic DNA was extracted from both swabs and tissues using a Qiagen DNA extraction kit (Blood and Tissue Kit; Qiagen, Germany) according to the manufacturer’s protocol.

For all 11 individuals from which we extracted tissues, we amplified one mitochondrial and one nuclear gene fragment. For the mtDNA, we sequenced 827 bp from a section of the genes 12S rRNA, the complete tRNA-Valine (Val) and 16S rRNA, using the primer pair F0001 (5’-AGA TAC CCC ACT ATG CCT ACC C-3’), R1169 (5’-GTG GCT GCT TTT AGG CCC ACT-3’) (Wilkinson et al. 2002). For the nuclear gene, we sequenced 476 bp of the Tyrosine exon-1 (TYR), using the primer pair L2976 (5’-TGC TGG GCR TCT CTC CAR TCC CA-3’), H2977 (5’-AGG TCC TCY TRA GGA AGG AAT G-3’) (Bossuyt and Milinkovitch 2000).

The Polymerase Chain Reactions (PCR) were carried out in 20 µl reaction with 50 to 100 ng of template DNA, with 1.0 µl of each primer (10 mM). The final concentration of each PCR reaction resulted to 1.5 µl of MgCl₂ (25 mM), 1.6 µl of dNTP (2.5 mM), 2.0 µl of 10× Buffer and 0.1 µl of TaKaRa Taq DNA polymerase (5 unit/µl). PCR amplifications were performed under the following thermal profiles: initial denaturation at 95 °C for 5 min, followed by 35 cycles with denaturation at 94 °C for 1 min, annealing at 55 °C for the mtDNA genes fragment and 54 °C for TYR for 1 min and extension at 72 °C for 1 min. The cycles were followed by a 10 min final extension at 72 °C. The amplified PCR products were sent for purification and sequencing to Cosmo Genetech Co. (Cosmo Genetech, Republic of Korea) on an ABI platform.

To reconstruct the independent and concantenated genes tree, we relied on two different datasets: (i) 827 bp-long fragments of mtDNA 12S rRNA, tRNA-Val and 16S rRNA (n taxa = 98), (ii) 486 bp-long fragments of sequences of protein-coding nuDNA Tyrosinase gene (TYR; n taxa = 110); and, (iii) 80 concatenated sequences of partial 12S rRNA (292 bp) and TYR (479 bp). We trimmed the sequences in each dataset manually and aligned the three sequences datasets indepedently using Clustal Omega (Sievers et al. 2011) in Geneious Prime (Kearse et al. 2012).

We calculated sequences similarity and estimated the genetic distance (or net evolutionary divergence) on the datasets of mtDNA 12S rRNA-trNA-Val-16S rRNA (n sequences = 98) and nuDNA TYR (n sequence = 110) using MEGA X (Kumar et al. 2018). We estimated the net average of evolutionary divergence between groups of sequences in each dataset; hence, we assigned 19 groups of species for 12S rRNA dataset and 16 groups of species for TYR dataset. In MEGA X, we conducted the analyses using the Maximum Compo-site Likelihood algorithm (Tamura et al. 2004) and modelled the rate variation amongst sites with a gamma distribution (shape parameter = 1). We considered differences in the composition bias amongst sequences in our evolutionary comparisons (Tamura and Kumar 2002), thus, all ambiguous positions were removed for each sequence pair using pairwise deletion option. These final datasets resulted in the totality of 116 positions, 301 sites for 12S rRNA and 110 positions, 486 sites for TYR.

For subsequent phylogenetic analyses, we downloaded supplemental sequences data of 98 homologous sequences of Kurixalus and Zhangixalus and other Rhacophoridae genera from Genbank (Wilkinson et al. 2002; Frost et al. 2006; Li et al. 2009; Nguyen et al. 2014; Yu et al. 2017a; Yu et al. 2018). GenBank accession numbers for both the new and previously deposited data are given in Table 1. We then created an initial alignment, based on nucleotide sequences with ClustalW2 (Larkin et al. 2007) and refined it manually. The final trimmed sequences resulted in 771 bp of concatenated 12S rRNA and TYR (n taxa = 79).

We used Partition Finder v. 2.1.1 (Lanfear et al. 2012) to determine the best-fit partitioning of the defined subsets. For the concatenated genes dataset, we defined four subsets by considering a fixed model for non-coding fragment and one subset for every single codon position with respect to the protein coding TYR gene fragments. Based on the Bayesian Information Criterion (BIC) values, we selected the following models for the following gene fragments: non-coding 12S rRNA (fixed subset): 1–292 (SYM+G) and protein coding TYR (subset 1): 293–771/1 (GTR+G); TYR (subset 2): 294–771/2 (K80+I+G) and TYR (subset 3): 295–771/3 (GTR+I+G). We used the models selected as a priori in further phylogeny analyses.

We built phylogenetic trees for all three datasets: mtDNA 12S rRNA-tRNA-Val-16S rRNA, protein coding nuDNA TYR and concatenated 12S rRNA-TYR using Bayesian Inference methodologies with MrBayes v.3.2.6 (Ronquist and Huelsenbeck 2003). For each tree dataset, we performed four separate analyses with 50 million generations of Markov Chain Monte Carlo and discarded the first 20 percent generations as burn-in until a convergence was reached (here we obtained > 0.005 split frequencies).

To test the presence of population differentiation using the TYR marker, we ran an analysis of molecular variance (AMOVA; Excoffier et al. 1992) using Arlequin v.3.5.2.2 (Excoffier and Lischer 2010). Here we used the AMOVA to test the three clades recovered from our phylogeny, based on TYR (Outgroup, Clade A and Clade B; see phylogenetic tree in Suppl. material 1: Fig. S2). In addition, we phased the diploid sequences of TYR gene (486 sites; n = 216) and analysed the haplotype using DnaSP v.5.0 (Librado and Rozas 2009). Before analysing the haplotypes, we assigned each haplotype group to its species, resulting 10 Kurixalus species and four closely-related species: K. inexpectatus sp. nov., K. banaensis, K. baliogaster, K. bisacculus, Kurixalus sp., K. hainanus, K. odontotarsus, K. verrucosus, K. eiffingeri, K. idiootocus and Zhangixalus appendiculatus, R. ocellatus, R. romeri and Orixalus carinensis. Out of the 486 sites, we disregarded missing haplotypes and removed all invariable sites. Then, we converted the haplotype analysis in DnaSP v.5.0 to an RDF input file format and we built the reticulated haplotype network from the phased TYR in NETWORK v.10.2.0 (Fluxus Technology Ltd; UK) using the Median-joining method (Bandelt et al. 1995).

Relying solely on a distance-based method is insufficient. The coalescent-based species delimitation was determined as the most efficient method for comparative study of species delimitation in genus of Kurixalus (Yu et al. 2017a). To test the assumption that the individuals samples belonged to a new species rather than an exotic or invasive clade of K. idiootocus, we employed a topology testing and species delimitation approach using both the coalescent-based methods. First, we designed two competing topology species tree models: model 1 and model 2 with two independent datasets consisting respectively of 79 unlinked sequences of mtDNA 12S rRNA (292 bp) and nuDNA TYR (451 bp). Model 1 designated the new Kurixalus clade, K. inexpectatus sp. nov. as clumped within the clade of most closely related species, Kurixalus idiootocus, whereas Model 2 assigned K. inexpectatus sp. nov. as a new species, split from K. idiootocus. For a comparison between topologies, we ran a nested sampling analyses on species tree Model 1 and Model 2 with NS package implemented in BEAST v.2.6.6 (Bouckaert et al. 2019). We selected MCMC sub-chain length of 10,000 with particle count of 10 and an epsilon of 1.0×10^-9 as parameters for each nested sampling analysis. Then, we evaluated topology of species trees of Model 1 and Model 2 by comparing the tree marginal L estimate (MLE) value and the Bayes factor obtained from the nested samplings. We calculated the Bayes factor with the following formula: Bayes factor = (MLE value of Model 1) – (MLE value of Model 2). We selected the best species tree model through the Bayes factor value, in which a positive Bayes factor is in favour of that particular model. We visualised the most likely species tree with Densitree (Bouckaert 2010).

Additionally, the recent study on the phylogeography of Taiwanese Kurixalus showed that the genus colonised the Island attributes through a land-bridge during the last glacial maxima (Yu et al. 2021). We further inferred the lineage origins and divergence between our focal taxa and K. idiootocus distributed in Taiwan Island. To do so, we estimated the time divergence of Kurixalus lineage by calibrating the species tree using an uncorrelated lognormal relaxed molecular clock with StarBeast RLC v.2.6.6 (Bouckaert et al. 2019). For both Model 1 and 2 datasets, we enforced three similar calibration points. Due to the absence of fossil records of Rhacophoridae in Asian mainland (Yu et al. 2021), we relied on paleogeological events for our primary calibration source. For secondary calibration, we adapted the range of molecular dating estimations of related literature (Pan et al. 2017; Yu et al. 2021). The three calibration points described as: (i) Emergence of Zhangixalus nigropunctatus in Southeast Asian and Chinese mainland ca. 11.39 Mya (High posterior density (HPD) 95%: 8.89–14.16; Pan et al. 2017), (ii) Emergence of stem group of South-eastern Asian clades of Kurixalus involves K. verrucosus group and its representative members ca. 7.4 Mya (HPD 95%: Yu et al. 2021) and (iii) Emergence of stem group of K. eiffingeri and K. idiootocus in Taiwan ca. 5.50 Mya (HPD 95%: 8.75 -3.25; Yu et al. 2021), simultaneously with the earliest island formation after physically separating from the south -astern mainland of China through the formation of the Taiwan strait (ca. 5.0-2.0 Mya modern shape; Teng 1990). We tested both birth-death and Yule priors on our species tree datasets and finally selected Yule as the best tree prior due to a better pattern of bifurcation for each crown node generated in trees. We ran four independent analyses, with MCMC chains of 20 million generations and 1,000 pre-burn-in steps for each dataset of model. We verified the convergence of the generated trees by evaluating the MCMC outputs with Tracer v.1.7.1 (Rambaut et al. 2018). Here, we ensured the values of effective sample size (ESS) for all parameters to be at least more than 1,000. We summarised a maximum clade credibility (MCC) tree for the calibrated tree time using TreeAnnotator, an application attached to BEAST v.2.6.6.

Finally, we projected the possible dispersal pathways, based on the molecular dating estimates focusing on the clade containing our focal species and Taiwanese Kurixalus on paleomaps using QGIS v.2.18.15. The oscillayers used to reconstruct the Plio-Pliocene maps was adapted from datasets provided in Gamisch (2019).

The acoustic recordings of putative new Kurixalus species were recorded between April and July in 2018 at 24 °C with a linear PCM recorder (Tascam DR-40; California, USA) linked to a unidirectional microphone (Unidirectional electret condenser microphone HT-81, HTDZ; Xi’an, China). To determine the relationship with other species, we first compared the number of consecutive calls within a series of calls between the individuals recorded and K. idiootocus, K. eiffingeri, K. berylliniris and K. wangi. We then compared the call properties of the new population with that of K. idiootocus as it was the most closely-related species and the only species with the same number of calls within a series of consecutive calls (see results). The recordings of K. idiootocus were obtained in central Taiwan (23.9240 N, 120.8910 E) in July 2013, using a Tascam DR-70D digital recorder (TEAC Corporation, Tokyo, Japan) and a Sennheiser ME67/K6 directional microphone (Sennheiser Electronic GmbH & Co. KG, Hanover, Germany). All our recordings were recorded at a sampling rate of 44.1 kHz with 16-bit resolution. Temperature was recorded with a Tecpel DIT-517 infrared thermometer (between 22 and 25 °C; TECPEL Corporation, New Taipei, Taiwan). The genus emits a series of continuous notes, pooled in bouts of continuous calls. To compare K. idiootocus and the new Kurixalus population, we selected one entire series of consecutive calls for each individual and analysed 373 advertisement calls in total, including 238 calls for K. idiootocus (9 to 24 calls in a bout from each of 16 males) and 135 calls for the new population (9 to 21 calls in a bout for 9 males).

We used Raven Pro v.1.5 (Cornell Lab of Ornithology 2011) to analyse our recordings. Nine properties, including six temporal and three spectral properties were measured in the two Kurixalus species (K. idiootocus vs. K. inexpectatus sp. nov): number of calls in a bout, bout length (s), call interval (ms), call length (ms), rise time (ms), fall time (ms), max frequency (kHz), 2^nd frequency (kHz), relative amplitude (dB). We measured the number of calls and the length of a series of consecutive calls. Call interval was measured as the duration from the end of a call to the beginning of the next call (Fig. 3A, B).

Figure 3. 

The call property measurements. This figure shows A the waveform of entire series of a consecutive call B the waveform C the spectrogram of two calls and D the spectral power distribution of a single call from the new Kurixalus population from Zhejiang, China. We extracted the number of calls in a bout, call interval, call duration (CD), rise time (RT), fall time (FT), dominant frequency (here also the max frequency), secondary peak frequency and the relative amplitude of two peaks.

Call duration refers to the time between the onset and offset of a call. Rise time refers to the time between the onset of call and the local maximum in the waveform. Fall time is the time between the local maximum in the waveform and the offset of a call (Fig. 3B). Dominant frequency is the strongest frequency in the duration of a call and the secondary frequency is the maximum frequency of the second harmonic. Relative amplitude is the difference in amplitude between the dominant frequency and secondary high (dominant – secondary). In both species, the dominant frequency was the primary harmonic, as indicated through preliminary analyses and later confirmed by determining the fundamental frequency for a subset of calls as the reciprocal of the average period of the quasi-periodic fine-temporal waveform. All frequency measurements were based on 1024-point fast Fourier transformation and Hann windows and were made from the average power spectrum computed over the duration of a call.

We first corrected the calls for temperature variation by adjusting the value of each variable to the average temperature of all recordings using the equation originating from the linear regression of each focal variable in function of temperature. As the contribution of each call property for each individual is not independent of other call properties and consequently correlated, we used a principal component analysis (PCA) to convert those call properties into a set of values of linearly uncorrelated factors. The PCA provided four principal components with Eigenvalues larger than 0.5, explaining 95.7% of the total variance. We used a Discriminant Function Analysis (DFA) to classify the call properties and test for the correctness of group assignment. We then plotted the two significant PCs against each other to illustrate the divergence between the two species. Finally, to determine the differing call variables between the two clades, we used a Multivariate Analysis of Variance (MANOVA) to compare each call property between these two species.

We collected eighteen morphological measurements three times each and averaged the values for further analyses. Morphometric data were taken using digital calipers to the nearest 0.1 mm and included the following characters: snout-vent length (SVL), head width (HDW), distance between left and right articulations of jaw, head length (HDL), from the tip of the snout to the articulation of the jaw, snout length (SNT), from tip of snout to the anterior corner of the eye, horizontal eye diameter (EYE) from the anterior to the posterior corner of the eye, width of the upper eyelid (UEW), the horizontal length of the upper eyelid, internares distance (IND), the distance from nostril to eye (DNE), from the posterior border of nostril to anterior border of the eye, narrowest interorbital distance (IOD), greatest horizontal tympanum diameter (TMP), tympanum-eye distance (TEY) from anterior edge of tympanum to posterior corner of eye, hand length (HND) from distal end of radio-ulna to tip of finger III, radio-ulna length (RAD), forelimb length (FLL), distance from the proximate end of radio-ulna to distal end of finger III, thigh length (THL), distance from vent to distal end of femur, tibia length (TIB), foot length (FL) from proximal end of inner metatarsal tubercle to tip of toe IV and the length of the foot and tarsus (TFL), distance from tibio-tarsal joint to tip of toe IV. All specimens were measured by a single author (YY) to minimise sampling error. The dataset is available Suppl. materials.

To be able to compare with the morphometrics of other clades, we extracted data from the literature for all species available (Suppl. material 1: Table S4): K. berylliniris, K. wangi, K. idiootocus, K. eiffingeri, K. bisacculus, K. lenquanensis, K. odontotarsus. K. yangi, K. naso, K. viridescens and K. ananjevae (Kuramoto and Wang 1987; Matsui and Orlov 2004; Nguyen et al. 2014; Tao et al. 2014; Wu et al. 2016; Yu et al. 2017b; Yu et al. 2018; Zeng et al. 2021). The data were, however, incomplete for the variables TEY, TFL, THL, HND and RAD and these variables were removed from the analyses. We only kept data points that had a full dataset for the remaining 13 variables and were males to enable further morphological comparison without impact of sexual dimorphism. In total, we harvested data for 68 Kurixalus sp. individuals, including our samples. We then removed variation due to size difference between individuals by dividing each of the variables by the SVL of the matching individual. Furthermore and because of the low sample size for most species, we created two categories for the statistical analyses: one including our focal clade and the other one with all other non-focal species. The reasoning behind this segregation being that a clade would have to be extremely divergent to be morphologically different from all other species of the genera. We also created a subsection of the dataset including our focal clade (n = 12) and K. idiootocus as it is the most phylogenetically closely-related clade (n = 6).

As the variables were strongly correlated (Pearson’s correlation; Table 2), we decided to use a factor reduction statistical analysis to identify the independent dimensions of the morphological characters. The principal component analysis was set such that principal components were to be extracted if their Eigen value was > 1, under a varimax rotation. Variables were selected as loading into a PC if loading > 0.58 (Table 3). Once the PCs were extracted, we tested for significant differences between our focal clade and all other species through one-way ANOVA and then between our focal clade and K. idiootocus through a second one-way ANOVA. All analyses were run in SPSS (SPSS, Inc., Chicago, USA). Additionally, after standardising morphological measurements by SVL, we also ran a two-sample t-test comparing the putative new species (n = 12) and K. idiootocus (n = 8) on 12 of the morphological characters.

Our analyses resulted in minor differences in the evolutionary divergence between the 12S rRNA gene fragments of K. inexpectatus sp. nov. and K. idiootocus (mean = 0.0004 SD ± 0.0004). Similarly, the protein coding nuclear TYR between K. inexpectatus sp. nov. and K. idiootocus showed a comparatively smaller mean of substitution rate (mean = 0.0035 ± 0.0004; value marked with double asterisks (**) in Table 4) than that of other species groups (Table 4).

Overall, the Bayesian Inference (BI) trees inferred from both mtDNA 12S rRNA-tRNA-Val-16S rRNA and nuDNA TYR fragments showed strong patterns of genetic structures for the East Asian and Southeast Asian Kurixalus phylogeny relationship, recovering four strongly supported clades (Clades A, B, C and D; see the distributions of the clades and the phylogenetic tree in Suppl. material 1: Fig. S1), including two major clades (Clades A and B; Suppl. material 1: Fig. S2) within the Kurixalus lineage. Although showing a discordant topology for the clades distributed in Southeast Asia (Suppl. material 1: Figs S1 and S2), the mtDNA and nuDNA trees converged towards a similar phylogenetic position for K. inexpectatus sp. nov., highlighting a sister relationship with K. idiootocus (Suppl. material 1: Figs S1, S2).

The phylogenetic relationship of concatenated gene fragments of partial 12S rRNA and TYR gene fragments recovered the three major clades within the Kurixalus genus with a Bayesian Posterior (BP) support of 90% for clade A, 52% for clade B, 71% for clade C (Fig. 4). Monophyletic clade A contained a Vietnamese Kurixalus (presumably as K. carinensis, but has been considered different from type species of K. carinensis in Myanmar) and K. romeri and K. ocellatus of Chinese mainland (BP = 100%; Fig. 4). Clade B recovered a monophyletic Kurixalus originated from Taiwan and south-eastern China. Here, the clades endemic to Taiwan, K. eiffingeri (Clade B, BP = 52%; Fig. 4) and K. idiootocus (BP = 97%; Fig. 4) were nested to the monophyletic clade of our Kurixalus sp. sampled in Zhejiang, south-eastern China (clade B, BP = 63%; Fig. 4). Clade B therefore supported a divergence between focal Kurixalus clade and K. idiootocus. Clade C comprised of a large nested monophyletic East Asian and Southeast Asian mainland Kurixalus, included clades of K. banaensis that were ranging in Vietnam and K. verrucosus (BP = 100%), K. baliogaster (BP = 100%), Kurixalus sp. (BP = 74%) K. bisacculus (BP = 67%) and K. odontotarsus (BP = 96%) originating from Western China, Tibet, Yunnan and the Eastern Indo-Chinese Peninsula: Vietnam, Thailand and Cambodia (Table 1; Fig. 4).

Figure 4. 

Bayesian Inference (BI) tree inferred from 79 sequences of concatenated 12S rRNA-TYR gene fragments. The three clades (Clades A, B and C) recovered in the phylogenetic tree are labelled accordingly. Clade B included the species K. inexpectatus sp. nov. described in the present study, indicated by the red box. The value of the node represents the Bayesian posterior probability (BPP) for each clade. The clades are marked with a solid bar and labelled in accordance with their specific name.

The AMOVA provided support to the genetic differentiation recorded while using the TYR marker as it identified 21.80% of variance within clades and 78.20% of variance between clades for the three main clades Kurixalus (n = 108; Fig. 4). The results of the AMOVA also provided a significant F_ST value (0.782; p < 0.05), showing that the three clades were significantly variable. The haplotype generated from the phased TYR fragment (n = 216), based in 477 trimmed sites, resulted in 62 haplotypes with a haplotype diversity (Hd) of 0.958 (Fig. 5). The distribution of haplotypes showed that six haplotype groups were representative of Kurixalus distributed in Taiwan Island and south-eastern China (Clade A; Fig. 5; see phylogenetic tree in Suppl. material 1: Fig. S2). Out of these six haplotype groups in Clade A, three of them represented K. inexpectatus sp. nov. (H4 – H6; Fig. 5). The origin of the haplotypes of K. inexpectatus sp. nov. corresponded to the haplotype groups of Taiwanese Kurixalus: K. idiootocus (H59- H60; Fig. 5) and K. eiffingeri (H 58; Fig. 5), whereas, Clade B contained a large portion of Kurixalus haplotype groups originating from south-eastern Asia, including Z. appendiculatus from Borneo and six Kurixalus species distributed across mainland Southeast Asia: Thailand, Laos, Cambodia and Vietnam (Fig. 5 and Suppl. material 1: Fig. S2). Clade B also included the representative haplotypes of K. hainanus distributed in Hainan Island and K. verrucosus haplotypes originating from Yunnan, south-western China.

Figure 5. 

Haplotype network inferred from 216 phased nuDNA TYR sequences data (486 sites). The haplotype group for the focal Clade A comprised six representative haplotypes of Kurixalus. Clade A included three K. inexpectatus sp. nov. haplotypes. The size of each haplotype marker matches the haplotype scales. The colour coding matches with the name of the taxa in the legend. The colours used for the boundaries of Clade A and Clade B are coded similarly to the colours of their clades in the phylogenetic tree (Suppl. material 1: Fig. S2).

The topology of the coalescent unlinked 12S rRNA and TYR tree supported the divergence of the focal Kurixalus clade from the most closely-related species, K. idiootocus. Nested sampling analyses on both species trees was favoured on the topology proposed by Model 2 (MLE = - 3688.252; Bayes factor: 651.011; Table 5). This topology provided support on the splitting between the lineages of K. inexpectatus sp. nov. distributed on south-eastern mainland and K. idiootocus distributed on Taiwan Island (Fig. 6), more so than a clumping between K. inexpectatus sp. nov. and K. idiootocus (see Model 2; Table 5).

Figure 6. 

Calibrated species tree of Rhacophoridae represented by Kurixalus, Orixalus, Romerus and Zhangixalus distributed over East Asia and Southeast Asia. The species tree reconstructed from unlinked 12S rRNA and TYR. The asterisk (*) symbol indicates Kurixalus inexpectatus sp. nov. The highlighted lineages divergence noted with (a–d) and the time estimates are synchronised with datation in Table 6. Biogeography models C and D hypothesised early colonisation pathway of Kurixalus to Taiwan Island and potential glacial-driven refugia to the south-eastern mainland, projected on the Plio-Pleistocene oscillations models of.

Our calibrated species tree of unlinked 12S rRNA and TYR gene fragments provided support on the earliest split between Asian lineages of Kurixalus and Zhangixalus to be dated in Mid-Miocene, ca. 11.17 Mya (Table 6; Fig. 5). This split-off was subsequently followed by the emergence of basal clade of Kurixalus in eastern Asia ca. 10.48 Mya [95% Highest Posterior Density (HPD): 8.16 – 12.98; node a; Table 6: Fig. 6]. Stem clade of Kurixalus distributed across the eastern Asian mainland, adjacent islands and south-eastern Asian mainland species group, consisting of members, such as K. banaensis, K. verrucosus. K. baliogaster, K. odontotarsus, Kurixalus sp. K. hainanus and K. bisacculus may have emerged ca. 9.14 Mya [6.87–11.12; node b; Table 6; Fig. 6]. Later, a stem clade of Taiwanese Kurixalus may have emerged, initiated by the isolation of K. eiffingeri in Taiwan and Ryukyu Islands ca. 5.66 Mya [3.32–8.07; node c; Table 6; Fig. 6]. Molecular dating estimates the lineage splitting between our proposed species, Kurixalus inexpectatus that distributed in south-eastern China and its sister clade, Taiwanese K. idiootocus to be in Late Pliocene to Pleistocene, ca. 3.06 Mya (5.82-0.01; node d; Table 6; Fig. 6).

As the contribution of each call property for each individual is not independent of other call properties and consequently correlated, we used a principal component analysis (PCA) to convert those call properties into a set of values of linearly uncorrelated factors. We selected the principal components from the results of the PCA to cover as much as possible of the total variance, resulting in four PCs with Eigenvalues larger than 0.5 and explaining 95.7% of the total variance. We used a Discriminant Function Analysis (DFA) to classify the call properties and test for the correctness of group assignment. We then plotted the two significant PCs against each other to illustrate the divergence between the two species. Finally, to determine the differing call variables between the two clades, we used a Multivariate Analysis of Variance (MANOVA) to compare each call property between these two species.

Based on the descriptions of the advertisement call and the number of calls in a series of consecutive calls, we could first segregate the species into two groups matching with the phylogenetic clustering. We grouped the putative new Kurixalus species and K. idiootocus together, while K. eiffingeri, K. berylliniris and K. wangi were grouped together (Suppl. material 1: Tables S1, S2, S3). From here, when then compared the putative new Kurixalus species and K. idiootocus.

The DFA on the four resulting PCs highlighted that only two PC1and PC3 were significantly different between the two species (PC1: Wilks’ Lambda = 0.93, F_(1,19) = 141.10, p < 0.001 ; PC3: Wilks’ Lambda = 0.17, F_(1,19) = 9.75, p = 0.005) and PC2 (Wilks’ Lambda = 0.11, F_(1,19) = 0.61, p = 0.442) and PC4 (Wilks’ Lambda = 0.12, F_(1,19) = 0.93, p = 0.345) were not. When plotting PC1 and PC3 against each other, a clear segregation of data was visible (Fig. 7). When comparing variables one by one for K. idiootocus and the new clade, the model was significant (MANOVA test; Wilks’ value = 0.087, F_9,14 = 16.27, p < 0.001) and numerous variables were different from each other. In detail, the new clade had longer call intervals, longer call duration and dominant frequency (Table 7).

Figure 7. 

The PCA plot of the Kurixalus inexpectatus sp. nov. (KX) and Kurixalus idiootocus (KI).

The unique PCA, used to identify the independent dimensions of the morphological characters between the individuals collected in this study and other Kurixalus sp. individuals, resulted in two PCs, with eigenvalues of 1.92 and 8.25, explaining a cumulated variation of 78.23% (Table 3). Based on the variables loading on to each of the PCs, we assigned PC1 to the general morphology and PC2 to the horizontal head structure (Table 3).

The results of the one-way ANOVA showed that there was no significant difference between the focal and non-focal groups for either of the PCs (Table 3). However, our focal clade and K. idiootocus were significantly different for both of the PCs (Table 3), such as PC1 (general morphology) p = 0.002 and PC2 (horizontal head structure) p = 0.017. When these variables were plotted against each other (Fig. 8), two non-clustering groups were visible, corresponding to variations between the measurements.

Figure 8. 

Plot of PC1 and PC2 resulting from the PCA and showing the non-clustering of morphological features between Kurixalus inexpectatus sp. nov. and Kurixalus idiootocus.

K. inexpectatus is morphologically most similar to K. idiootocus, its closest relative and, after standardising measurements by SVL, K. inexpectatus differs by having a relatively longer head length (34% vs. 33%), significantly shorter snout (13% vs. 15%; p < 0.001), significantly greater internasal distance (11% vs. 10%; p < 0.001), significantly smaller eye diameter (13% vs. 16%; p < 0.001), nearly significant wider tympanum diameter (7% vs. 6%; p = 0.06), significantly greater distance between the eyes and nares (8% vs. 7%, p = 0.03), significantly longer forelimb length (50% vs. 48%; p = 0.03), shorter tibia length (44% vs. 45%) and longer foot length (42% vs. 40%). Additionally, K. inexpectatus is further distinguished from K. idiootocus in having a tibio-tarsal articulation that extends beyond the anterior corner of the eye (versus centre of eye). K. inexpectatus can be differentiated from K. bisacculus, K. hainanus, K. naso, K. odontotarsus, K. raoi, K. silvaenaias, K. verrucosus and K. yangi by having an average adult SVL of less than 30 mm (27.5 – 31.8, × = 29.2) (vs. larger) (Yu et al. 2018; Hou et al. 2021; Zeng et al. 2021). K. inexpectatus can be further differentiated from K. absconditus, K. baliogaster, K. banaensis, K. berylliniris, K. chaseni, K. eiffingeri, K. gracilloides, K. lenquanensis, K. motokawai, K. viridescens and K. wangi by the presence of a pair of large, symmetrical dark blotches on the chest (vs. absent; Yu et al. 2018; Hou et al. 2021; Zeng et al. 2021). Kurixalus inexpectatus is distinguished from K. ananjevae by having limbs with serrated dermal fringes (vs. smooth; Yu et al. 2018; Hou et al. 2021; Zeng et al. 2021).

The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature and, hence, the new names contained herein are available under that Code from the electronic edition of this article. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix “http://zoobank.org/”. The LSID for this publication is: urn:lsid:zoobank.org:pub: 3CCB356B-F075-4EE5-8366-FE96B855F884. The new species name Kurixalus inexpectatus sp. nov. has been registered under LSID: urn:lsid:zoobank.org:act: 02D394DE-BB1C-4C17-BB70-656D68814C8F.