Skip to main content

Advertisement

To split or not to split? Multilocus phylogeny and molecular species delimitation of southeast Asian toads (family: Bufonidae)

Abstract

Background

Recent studies have demonstrated that Bayesian species delimitation based on the multispecies coalescent model can produce inaccurate results by misinterpreting population splits as species divergences. An approach based on the genealogical divergence index (gdi) was shown to be a viable alternative, especially for delimiting allopatric populations where gene flow is low. We implemented these analyses to assess species boundaries in Southeast Asian toads, a group that is understudied and characterized by numerous unresolved species complexes.

Results

Multilocus phylogenetic analyses showed that deep evolutionary relationships including the genera Sigalegalephrynus, Ghatophryne, Parapelophryne, Leptophryne, Pseudobufo, Rentapia, and Phrynoides remain unresolved. Comparison of genetic divergences revealed that intraspecific divergences among allopatric populations of Pelophyrne signata (Borneo vs. Peninsular Malaysia), Ingerophrynus parvus (Peninsular Malaysia vs. Myanmar), and Leptophryne borbonica (Peninsular Malaysia, Java, Borneo, and Sumatra) are consistent with interspecific divergences of other Southeast Asian bufonid taxa. Conversely, interspecific divergences between Pelophryne guentheri/P. api, Ansonia latiffi/A. leptopus, and I. gollum/I. divergens were low (< 3%) and consistent with intraspecific divergences of other closely related taxa. The BPP analysis produced variable results depending on prior settings and priors estimated from empirical data produced the best results that were also congruent with the gdi analysis.

Conclusions

This study showed that the evolutionary history of Southeast Asian toads is difficult to resolve and numerous relationships remain ambiguous. Although some results from the species delimitation analyses were inconclusive, they were nevertheless efficacious at identifying potential new species and taxonomic incompatibilities for future in-depth investigation. We also demonstrated the sensitivity of BPP to different priors and that careful selection priors based on empirical data can greatly improve the analysis. Finally, the gdi can be a robust tool to complement other species delimitation methods.

Background

The instability of Southeast Asian amphibian systematics cannot be overstated and reflects the rate at which new insights are gained, largely through the use of molecular data [1]. The increased effort to sequence rare, historical, and widespread taxa at a finer genetic and geographical scale has brought upon a renewed understanding of the region’s complex biodiversity patterns and evolutionary history. Consequently, taxonomic nomenclature is constantly in flux [2,3,4,5,6,7,8,9,10,11], as it attempts to keep pace with the increasing frequency at which new species are being described (especially cryptic species) and revised [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Although these insights have improved our understanding of Southeast Asian amphibian biodiversity considerably, many taxa have yet to be sequenced or subjected to molecular analyses, and more fine-scale geographic sampling is still needed to adequately characterize species boundaries and distribution ranges.

True toads (family Bufonidae) have a native distribution throughout Asia, Africa, Europe, and America, but have also been widely introduced to Australia and the Oceanic islands [27]. Southeast Asian taxa represent a relatively small portion of the family’s diversity and consequently, have received comparatively less research attention. However, studies have begun to show that the diversity of Southeast Asian toads is severely underestimated and much of their evolutionary history remains unresolved [13, 20, 21, 26, 28,29,30]. Moreover, most research focused on small subsets relating to specific taxon groups, thereby precluding the detection of broad biodiversity and evolutionary patterns that can only be revealed by more comprehensive studies conducted at a broader scale [5].

Although many Southeast Asian toad species have been sequenced and are publicly available on GenBank, these sequences have not been collated and analysed within a comprehensive phylogenetic framework. This includes widespread taxa that are co-distributed across separate landmasses and distinct biogeographic regions such as Borneo, Sumatra, Java, Peninsular Malaysia, and Indochina. Such taxa are of particular interest because they are phenotypically similar, yet have not been analysed within a broader phylogeographic context. Studies that have included wider geographic sampling have revealed that widespread species were composed of multiple undescribed lineages with more fragmented distributions [7, 16, 21, 31,32,33,34,35,36,37,38,39]. However, no such study has been conducted on Southeast Asian toads, thereby highlighting the urgent need for more in-depth studies to fill the obvious gap in our understanding of the group’s evolution, biodiversity, and distribution. This study aims to provide a comprehensive working phylogeny and preliminary analysis of genetic variation, phylogeographic structure, and species boundaries, with emphasis on cryptic complexes that are co-distributed across multiple biogeographic regions.

The use of molecular data and development of statistical species delimitation methods have become indispensable tools to aid in the detection and delimitation of new species [40, 41]. Among these, methods based on the multispecies coalescent (MSC) model are largely regarded as robust, especially for closely related species. However, studies based on simulations and empirical data have revealed that certain MSC-based methods may over-split by delimiting populations as species [22, 42, 43]. Widely used species delimitation programs such as Bayesian Phylogenetics and Phylogeograpy (BPP; Yang and Rannala, 2010) were found to be unable to differentiate population structure from species cladogenesis when sequence data were simulated under the protracted speciation model [42]. Subsequently, the empirical genealogical divergence index (gdi) was proposed and found to be more reliable than BPP in differentiating population structure from species divergence, especially for allopatric populations [43, 44]. As such, the main goals of this study are to: 1) estimate a comprehensive multilocus phylogeny of Southeast Asian bufonids to elucidate ambiguous evolutionary relationships; and 2) implement a multi-step approach to species discovery and delimitation to screen for potential new species and taxonomic incompatibilities. We first applied mitochondrial divergence thresholds specific to amphibians to identify candidate species for species delimitation analyses [45, 46]. We then used BPP to test putative species boundaries [47] and the gdi analysis to validate the results [43, 44].

Results

Phylogenetic estimation

The best-fit substitution model scheme for the ML analysis was GTR + F + R5 for the partitions 12S and 16S; TIM2 + F + I + G4 for the partition CO1; and TN + F + G4 for the partitions CXCR4, NCX1, and RAG-1. Both ML and BEAST phylogenies produced largely congruent results. At the generic level, the ML tree was poorly supported at the basal nodes, including the genera Sigalegalephrynus, Parapelophryne, Leptophryne, Pseudobufo, Ghatophryne, Rentapia, and Phrynoidis. Relationships of the other genera (Sabahphrynus, Ingerophrynus, Pelophryne, and Ansonia) were relatively well resolved. The BEAST tree differed only in the placement of the genus Ghatophryne, which was recovered as the sister lineage to Parapelophryne + Leptophryne, as opposed to the ML tree, where it was the sister lineage to Rentapia + Phrynoidis. Species-level relationships were also largely congruent, with a few exceptions including the placement of Ansonia kraensis, A. inthanon and I. gollum. The phylogenetic placement of Ingerophrynus gollum and its closest relative, I. divergens was not well-resolved in both Bayesian and ML phylogenies. In the Bayesian phylogeny, I. divergens was paraphyletic with respect to I. gollum, however, the clade lacks statistical support (Additional file 2). In the ML phylogeny, I. gollum and I. divergens were recovered as sister species but the entire clade lacked phylogenetic structure and support (Additional file 1). Branch support was relatively high for both phylogenies except for the most basal splits. While the ML phylogeny recovered a polytomy at the base of the tree, the BEAST phylogeny recovered the genus Sigalegalephrynus as the earliest diverging ingroup with high support (PP = 1.0; Fig. 1, Additional files 1, 2).

Fig. 1
figure1

Bayesian phylogeny of derived from 6236 bp comprising three mitochondrial (12S, 16S, CO1) and three nuclear genes (CXCR4, NCX1, RAG-1). Black dots denote branch support with PP ≥0.95, while red dots denote 0.9 ≤ PP < 0.95. All images in the figure were taken by the first author

Identification of candidate species

Below the 3% threshold, three taxon pairs were identified as candidates for lumping (Fig. 2). These include Pelophryne guentheri/P. api (2.2%), Ansonia leptopus/A. latiffi (1.2–2.3%), and Ingerophrynus gollum/I. divergens (1.9–2.7%). Intraspecific divergences among populations of Pelophryne signata from Peninsular Malaysia exceeded this threshold (0.7–3.6%). From 4 to 5%, another taxon pair, Leptophryne cruentata/L. javanica (5%) was identified for lumping, whereas two other taxon pairs were identified for splitting. These include, populations of Leptophryne borbonica from Borneo vs. Sumatra and from Java vs. Peninsular Malaysia (PM). Above the 5% threshold, three taxa were earmarked for splitting including Pelophryne signata from Borneo vs. PM (5.2–6.4%), Ingerophrynus parvus from PM vs. Myanmar (5.1–5.5%), and L. borbonica from Java+PM vs. Sumatra+Borneo (6.1–10.1%; Fig. 2). Individual pairwise distances for all populations are presented in Additional file 5.

Fig. 2
figure2

Threshold-based approach to identify candidate species for lumping and splitting using uncorrected p-distances at the 16S gene. Dotted red lines show divergence thresholds of 3–5%. All images in the figure were taken by the first author

Species delimitation

A comparison between posterior distributions generated from the priors and empirical data showed that there was sufficient information in the data and that results were not driven by the priors (Additional file 3). Overall, the A10 and A11 analyses produced similar results, indicating that the estimated species tree did not differ from the inferred Bayesian phylogeny. The exception was the comparison between Ingerophrynus gollum/divergens where posterior probabilities for the A10 and A11 analyses using empirical priors were markedly different (0.58 vs. 0.86). This was not unexpected considering the uncertain phylogenetic placement of this species (Additional files 1 and 2). Both A10 and A11 analyses were sensitive to prior settings and posterior probabilities based on empirical priors were generally lower than those based on pre-defined priors (Tables 1 and 2). For the candidates for lumping, certain non-empirical prior sets were highly supported but support was weak when empirical priors were employed (Table 1). This was corroborated by the uncertain values in the gdi analysis (Table 3). For the candidates for splitting, BPP analyses based on empirical priors and the gdi highly supported the distinction between Leptophryne borbonica from Java + PM/Sumatra + Borneo and Pelophryne signata/cf. signata as separate species (Tables 2 and 3). The gdi analysis also supported the split between L. borbonica from Sumatra vs. Borneo (mean ± SD = 0.70 ± 0.19; Fig. 3) and between Ingerophrynus parvus from PM vs. Myanmar (mean ± SD = 0.71 ± 0.20; Fig. 4). However, these were weakly and moderately supported in the BPP analyses (pp = 0.68 and 0.87 respectively). Complete results for every A11 analyses performed are provided in Additional file 6.

Table 1 Posterior probabilities from both BPP analyses (A10 | A11) for the candidates for lumping under different prior settings. A10 = species delimitation on a guide tree; A11 = joint species delimitation and species tree estimation
Table 2 Posterior probabilities from both BPP analyses (A10 | A11) for the candidates for splitting under different prior settings. A10 = species delimitation on a guide tree; A11 = joint species delimitation and species tree estimation
Table 3 Summary of the different lines of evidence used for species delimitation. BPP support summarizes results from A10 and A11 analyses using empirical priors. PP ≥ 0.95 denotes high support, 0.90 ≤ pp. < 0.95 denotes moderate support, and pp. < 0.9 denotes weak support
Fig. 3
figure3

Left: subset of the Bayesian phylogeny showing the genera Sigalegalephrynus, Ghatophryne, Parapelophryne, Leptophryne, Pseudobufo, Rentapia, and Phrynoidis. Highlighted tips are candidate taxa included in species delimitation analyses. Top right: boxplots of uncorrected p-distances (16S) for the genus Leptophryne. Bottom right: density plot of gdi values. Populations are considered distinct species when gdi values are > 0.7, while low gdi values < 0.2 indicate that populations belong to the same species. Values of 0.2 > gdi < 0.7 indicate ambiguous species status

Fig. 4
figure4

Left: subset of the Bayesian phylogeny showing the genera Sabahphrynus and Ingerophrynus. Highlighted tips are candidate taxa included in species delimitation analyses. Top right: boxplots of uncorrected p-distances (16S) for the genus Ingerophrynus. Bottom right: density plot of gdi values

Discussion

Evolutionary relationships

The phylogenetic relationships inferred from this study bear similarities as well as disagreements with previous studies. The most recent and comprehensive phylogenetic study on the family Bufonidae focused on Arabian taxa but included numerous Southeast Asian species [48]. However, the study lacked a number of key genera such as Pseudobufo, Sigalegalephrynus, and Parapelophryne that were unavailable from GenBank at that point in time. Similar to our results, the backbone of their tree was mostly unresolved and poorly supported. Their phylogeny differed in a number of major relationships, including the placement of the Rentapia + Phrynoidis clade, which they recovered as sister to the Pelophryne + Ansonia clade. Additionally, the genus Sabahphrynus (endemic to Borneo) was recovered as the sister lineage of Strauchbufo raddei from north and east Asia [27], albeit with low support. Our results inferred Sabahphrynus as the sister lineage to the Southeast Asian genus Ingerophrynus, which we favour as a more biogeographically plausible scenario that is also strongly supported in both our phylogenies. The morphological similarities between Sabahphrynus and other Southeast Asian taxa and its obvious disparity with Strauchbufo and other Eurasian taxa also supports this hypothesis. Another relevant multilocus study that focused on the description of Sigalegalephrynus from Indonesia [13] included most Southeast Asian genera, but the relationships of major clades were poorly supported and differed drastically from other studies [48, 49], including ours. Overall, our phylogeny bears the most congruence to the phylogeny of Pyron and Wiens [49] with regards to the placements of Phrynoidis, Rentapia, Sabahphrynus, Ingerophrynus, Ansonia, and Pelophryne. Their phylogeny however, recovered the Southeast Asian genus Leptophryne as the sister lineage to the European genus Epidalea, which we consider to be less parsimonious from a biogeographic perspective.

What all Southeast Asian toad phylogenies share in common, is the poor support along the backbone of the tree, possibly due to the lack of informative characters/genes or evolutionary processes such as hybridization, rapid radiation and incomplete lineage sorting that can obfuscate phylogenetic signal [50,51,52,53,54]. Even though the most robust phylogenies were estimated using 12 loci [48, 49], basal nodes remained poorly supported, especially with regard to a number of enigmatic genera such as Pseudobufo, Ghatophryne, Parapelophryne, and Sigalegalephrynus. These taxa are crucial to the understanding of Southeast Asian biodiversity as they occur in the fringes of the region and represent highly specialized, endemic, and possibly relictual lineages. Moreover, these taxa are represented by relatively few specimens and even fewer molecular samples. The inadequacies of current data to resolve the evolutionary history of Southeast Asian toads highlight the urgent need for more intensive collection efforts and high throughput sequencing at the genomic level to bolster exiting data.

Reconciling results: limitations and recommendations

Results from our study demonstrate that cryptic species delimitation can be a non-trivial task, especially for taxa in the grey zone of the speciation continuum (0.5–2% molecular divergence) [55]. In some cases, morphological differences can evolve faster than the formation of a species barrier (depressed hybrid fitness), but in others, morphological differentiation may not accompany genetic divergence [56]. This can lead to taxonomic discrepancies and incongruent species delimitation conclusions, especially if different operational criteria are used for species diagnosis.

The decision to split or lump can be straightforward in cases where multiple independent lines of evidence are in agreement. However, different methods can produce conflicting results that are often difficult to interpret and reconcile. The BPP analysis was sensitive to certain prior settings, especially the age of the root (τ0). Our results showed that increasing τ to larger values such as 0.04 or 0.4 produced in our opinion, spuriously high posterior probabilities (pp = 1.0 for all tested candidate taxa), possibly because large values of τ pushes coalescent events closer to the tips. The effects of changing θ were less obvious. Larger θ increased posterior probabilities in some cases but not others, indicating that population size and mutation rate need to be estimated on a case-by-case basis depending on the taxa. This demonstrates the sensitivity and potentially misleading results that can be produced by BPP when unrealistic priors are used. Results were more reasonable and congruent with the gdi analysis when empirical priors were employed. We therefore caution against using a range of priors that may not accurately represent the diversification and demographic history of the taxa in question, but instead, employ empirical priors that are estimated from the data.

Candidates for lumping

Leptophryne cruentata and L. javanica were relatively well-diverged (5%) but were not highly supported as separate species. However, both these taxa were represented by singletons and this was reflected in the gdi analysis that produced a flat and diffused posterior distribution. A similar result was seen in Pelophryne api/P. guentheri, which were also represented by singletons. These results should therefore be interpreted with caution and demonstrates the potential impact of sampling density on species delimitation analyses. For Ingerophrynus gollum and I. divergens, result of the BPP A10 analysis based on empirical priors was weak (0.6) whereas the A11 analysis was moderate (0.9). The gdi analysis was largely uncertain while the mitochondrial divergence was relatively low (2–3%). This, coupled with the relatively short branch lengths, suggests that these lineages were recently separated on their individual landmasses (PM and Borneo), likely during the Pleistocene glaciations and may still be in the early stages of speciation. Surprisingly, although the genetic divergence between Ansonia latiffi and A. leptopus was low (1–2%), the BPP and gdi analysis provided moderately strong support in favour of their recognition as distinct species. This could be due to the significant genetic substructuring and uncertain phylogenetic relationships among the different geographic populations, indicating that the A. leptopus group could potentially represent a species complex that requires further investigation. In general, our results did not unequivocally support any of the candidates for lumping, indicating that these taxa could potentially be conspecific or that these lineages are in the nascent stages or grey zone of speciation due to their relatively recent divergence.

Candidates for splitting

The genetic divergences among populations of the Leptophryne borbonica complex range from 3.9–10.1%, which are clearly within the range of interspecific divergences of other bufonid and amphibian taxa [21, 45, 57, 58]. Based on the type locality, we consider the population from Java to be the true L. borbonica, which is the sister lineage to the population from PM. Although the genetic divergence between these two populations was relatively high (4.6%), their recognition as distinct species was not strongly supported by the BPP and gdi analyses. However, both populations were represented by singletons and thus, more sampling will be required before stronger conclusions can be made. The genetic divergence between populations from either Sumatra or Borneo with the true L. borbonica from Java is high (9–10%; Fig. 3), clearly indicating that they could represent distinct species. When Sumatra + Borneo and Java + PM were compared, all analyses strongly supported those lineages as distinct species. These results show that populations from Sumatra and Borneo are distinct from the true L. borbonica from Java and warrants specific recognition. However, it remains unclear whether populations from Sumatra and Borneo, or Java and PM should be lumped as the results showed that they could potentially be distinct from each other as well. Similarly, the genetic divergence between the true Pelophryne signata from Borneo and its sister lineage from PM and Sumatra (cf signata) is 5.2–6.4%, which is consistent with interspecific divergences of other closely related taxa (Fig. 5). Additionally, the BPP and gdi analyses highly supported the recognition of these populations as distinct species.

Fig. 5
figure5

Left: subset of the Bayesian phylogeny showing the genera Pelophryne and Ansonia. Highlighted tips are candidate taxa included in species delimitation analyses. Top right: boxplots of uncorrected p-distances (16S) for the candidate taxa. Bottom right: density plot of gdi values

The type locality of Ingerophrynus parvus is in PM, but has a wide distribution ranging from southern Myanmar, parts of Indochina, along the Malay Peninsula, Sumatra, and western Java [27]. Unfortunately, only samples from Myanmar and PM were available for analysis. Our results showed that these populations are relatively well-diverged (~ 5%), received moderate support in BPP analyses, and high support in the gdi analysis for being distinct species. Denser population-level sampling throughout its wide distribution range will likely reveal this group to be a species complex with considerable genetic and phylogeographic structure.

Conclusions

This study showed that the evolutionary history of Southeast Asian toads is difficult to resolve using traditional Sanger-based loci and that numerous taxon groups require further revision. Although some results from the species delimitation analyses were inconclusive, they were nevertheless efficacious at identifying numerous potential new species and taxonomic incompatibilities for future in-depth investigation. We hold formal taxonomic changes in abeyance, pending the examination of type material and more robust analyses that incorporates more genetic loci and the consideration of additional, independent lines of evidence.

Methods

Taxon sampling

We obtained DNA sequence data from 107 ingroup samples including all nine Southeast Asian genera (Ansonia, Ingerophrynus, Leptophryne, Pelophryne, Phrynoidis, Pseudobufo, Rentapia, Sabahphrynus, Sigalegalephrynus), representing approximately 50 of the 68 species known to occur in Southeast Asia [27]. Sequences for three mitochondrial (12S, 16S, CO1) and three nuclear genes (CXCR4, NCX1, RAG-1) were acquired from GenBank, including six additional newly sequenced samples from Malaysia (Additional file 4). New sequences were obtained from tissue samples loaned from the La Sierra University Herpetological Collection, Riverside, California, USA and no animals/specimens were used in this study. The six genetic markers were selected because they had the highest coverage across the ingroup taxa. DNA for the new samples was extracted with the Promega Maxwell© RSC Instrument using the Maxwell© RSC Tissue DNA Kit. We used the primers 16Sc-L (5′-GTRGGCCTAAAAGCAGCCAC-3′), and 16Sd-H (5′-CTCCGGTCTGAACTCAGATGACGTAG- 3′) to amplify and sequence the 16S gene [59]. Amplification was done using the following PCR thermal profile: 95 °C for 4 min, followed by 35 cycles of 95 °C for 30 s, 52 °C for 30 s, 72 °C for 70 s, and a final extension phase at 72 °C for 7 min [31]. Amplified DNA products were subsequently visualized on 1.0% agarose gels and sequenced at Genewiz, Frederick, MD. Sequences were assembled, aligned (MUSCLE algorithm), and concatenated in Geneious Pro 5.3 [60] prior to phylogenetic estimation. Collection of specimens used in this study complied with institutional guidelines.

Phylogenetic analysis

The concatenated sequence matrix was partitioned by gene and the program IQ-TREE v1.6.8 was used to estimate a maximum likelihood (ML) phylogeny [61, 62]. The best-fit substitution model for each partition under the Bayesian information criterion was inferred using ModelFinder [63]. Partitions were allowed to merge by executing the “–m MFP + MERGE” function and the best partition scheme was obtained by examining the top 10% partition-merging schemes using the “–rcluster 10” command [64]. Branch support was assessed with 5000 bootstrap replicates using Ultrafast Bootstrap Approximation (UFB; Hoang et al., [65]). UFB values of 95 and above were considered well supported. We also estimated a Bayesian phylogeny using BEAST v2.5 [66], which was implemented through the CIPRES portal [67]. The best-fit substitution model for each partition was estimated via model averaging using the bModelTest plugin in BEAST [68]. A relaxed log-normal and Yule model was selected for the molecular clock and tree priors with all other priors set to default values. We executed two separate MCMC chains at 50,000,000 generations each and checked for convergence using the program Tracer v1.7 [69]. Sampled trees from both MCMC chains were combined using logcombiner and a maximum clade credibility tree was constructed using treeannotator with a burn-in of 10%.

Species delimitation

Candidate species

We first grouped samples into populations to form the fundamental units for species delimitation analyses. Based on inferred relationships from the concatenated BEAST phylogeny, populations were defined as monophyletic lineages that also occurred on separate landmasses or biogeographic regions (e.g. PM, Borneo, Java, Sumatra, Sulawesi, Indochina). We then implemented a threshold-based approach to screen for populations with taxonomic incompatibilities (highly diverged conspecific lineages or non-conspecific lineages with low divergences) [5, 45, 57, 70, 71]. To obtain an overview of intra- vs. interspecific genetic variation, pairwise uncorrected p-distances of the 16S mitochondrial gene were calculated using the program PAUP* [72]. We applied thresholds of 3–5%, which have been shown to effectively capture sequence divergences of nascent and fully diverged species in amphibians [45, 46]. Non-conspecific taxon pairs below the thresholds were flagged as candidates for lumping, while conspecific pairs above the thresholds were identified as candidates for splitting.

BPP analysis

We used the program BPP v4 [47] to perform species delimitation analyses on the candidate groups. This analysis was based on a single 16S gene, which had the best coverage across all taxa. The A10 analysis (species delimitation using a fixed guide tree) was implemented using relationships derived from phylogenetic analyses as a guide tree. The guide tree was represented by species-level relationships and conspecific individuals were combined under the same nominal identifier. The mapping of individuals to their identifiers was done in the imap file supplied to BPP. To explore the sensitivity of the BPP analysis to different prior settings, we tested a series of 10 different prior sets, using α = 3 in the inverse-gamma as a diffuse prior, and α = 21 as an informative prior. The prior mean was then set to vary by two orders of magnitude [73] while the MCMC was set to 100,000 samples with burnin = 10,000 and sample frequency = 5. Convergence was assessed by comparing the consistency of posterior distributions [44, 73, 74]. To accommodate uncertainty in the guide tree, we also performed the A11 analysis (joint species delimitation and species-tree estimation) using the same prior sets. Because these prior settings may not accurately represent the diversification and demographic history of the taxa, we performed separate A10 and A11 analyses using empirical priors calculated from the data. We used α =3 for both θ and τ priors and the corresponding β parameter was adjusted according to the mean (m) estimate of nucleotide diversity (for θ) and node height (for τ) using the equation m = β/(α-1), for α > 2 [73]. The average node height for each candidate group was obtained from the BEAST analysis. Posterior probabilities, pp. ≥ 0.95 were considered highly supported, 0.90 ≤ pp. < 0.95 were considered moderately supported, whereas pp. < 0.9 were considered weakly supported. Finally, we compared posterior distributions generated from priors only, versus those that were generated from empirical data to determine whether the data contained sufficient information for species delimitation and to ensure that the results were not driven by the priors.

Heuristic gdi analysis

We implemented the A00 analysis in BPP to generate posterior distributions for the parameters τ and θ using empirical priors. Four separate runs were performed to ensure convergence and converged runs were combined to generate posterior distributions for the MSC parameters that were subsequently used to calculate the gdi following the equation: gdi = 1 – e-2τ/ θ [43, 44]. Population A is distinguished from population B by using 2τABA, while 2τABB is used to differentiate population B from population A. Populations are considered distinct species when gdi values are > 0.7, while low gdi values < 0.2 indicate that populations belong to the same species. Values of 0.2 > gdi < 0.7 indicate ambiguous species status [43, 75].

Abbreviations

BPP:

Bayesian Phylogenetics and Phylogeography

gdi:

Genealogical divergence index

MCMC:

Markov Chain Monte Carlo

ML:

Maximum likelihood

MSC:

Multispecies Coalescent

PM:

Peninsular Malaysia

PP:

Posterior probability

rjMCMC:

Reversible-jump Markov Chain Monte Carlo

References

  1. 1.

    Brown RM, Stuart BL. Patterns of biodiversity discovery through time: an historical analysis of amphibian species discoveries in the southeast Asian mainland and adjacent island archipelagos. In: Gower DJ, Johnson K, Richardson J, Rosen B, Ruber L, Williams S, editors. Biotic evolution and environmental change in Southeast Asia. Cambridge: Cambridge University Press; 2012. p. 348–89.

  2. 2.

    Frost DR, Grant T, Faivovich J, Bain RH, Haas A, Haddad CFB, et al. The amphibian tree of life. Bull Am Museum Nat Hist. 2006;297:1–291.

  3. 3.

    Oliver LA, Prendini E, Kraus F, Raxworthy CJ. Systematics and biogeography of the Hylarana frog (Anura: Ranidae) radiation across tropical Australasia, Southeast Asia, and Africa. Mol Phylogenet Evol. 2015;90:176–92. https://doi.org/10.1016/j.ympev.2015.05.001.

  4. 4.

    Chan KO, Grismer LL, Zachariah A, Brown RM, Abraham RK. Polyphyly of Asian tree toads, genus Pedostibes Günther, 1876 (Anura: Bufonidae), and the description of a new genus from Southeast Asia. PLoS One. 2016;11:e0145903.

  5. 5.

    Chan KO, Grismer LL, Brown RM. Comprehensive multi-locus phylogeny of Old World tree frogs (Anura: Rhacophoridae) reveals taxonomic uncertainties and potential cases of over- and underestimation of species diversity. Mol Phylogenet Evol. 2018;127:1010–9.

  6. 6.

    Chen L, Murphy RW, Lathrop A, Ngo A, Orlov NL, Cuc TH, et al. Taxonomic chaos in Asian ranid frogs: an initial phylogenetic resolution. Herpetol J. 2005;15:231–43.

  7. 7.

    Chen J-M, Poyarkov NA, Suwannapoom C, Lathrop A, Wu Y-H, Zhou W-W, et al. Large-scale phylogenetic analyses provide insights into unrecognized diversity and historical biogeography of Asian leaf-litter frogs, genus Leptolalax (Anura: Megophryidae). Mol Phylogenet Evol. 2018;124 March:162–71. https://doi.org/10.1016/j.ympev.2018.02.020.

  8. 8.

    Chen JM, Zhou WW, Poyarkov NA, Stuart BL, Brown RM, Lathrop A, et al. A novel multilocus phylogenetic estimation reveals unrecognized diversity in Asian horned toads, genus Megophrys sensu lato (Anura: Megophryidae). Mol Phylogenet Evol. 2017;109:466. https://doi.org/10.1016/j.ympev.2016.11.019.

  9. 9.

    Li JT, Che J, Bain RH, Zhao EM, Zhang YP. Molecular phylogeny of Rhacophoridae (Anura): a framework of taxonomic reassignment of species within the genera Aquixalus, Chiromantis, Rhacophorus, and Philautus. Mol Phylogenet Evol. 2008;48:302–12.

  10. 10.

    Stuart BL. The phylogenetic problem of Huia (Amphibia: Ranidae). Mol Phylogenet Evol. 2008;46:49–60.

  11. 11.

    Li JT, Che J, Murphy RW, Zhao H, Zhao EM, Rao DQ, et al. New insights to the molecular phylogenetics and generic assessment in the Rhacophoridae (Amphibia: Anura) based on five nuclear and three mitochondrial genes, with comments on the evolution of reproduction. Mol Phylogenet Evol. 2009;53:509–22. https://doi.org/10.1016/j.ympev.2009.06.023.

  12. 12.

    Arifin U, Smart U, Hertwig ST, Smith EN, Iskandar DT, Haas A. Molecular phylogenetic analysis of a taxonomically unstable ranid from Sumatra, Indonesia, reveals a new genus with gastromyzophorous tadpoles and two new species. Zoosystematics Evol. 2018;94(1):163–93.

  13. 13.

    Smart U, Sarker GC, Arifin U, Harvey MB, Sidik I, Hamidy A, et al. A new genus and two new species of arboreal toads from the highlands of Sumatra with a phylogeny of Sundaland toad genera. Herpetologica. 2017;73:63–75. https://doi.org/10.1655/Herpetologica-D-16-00041.

  14. 14.

    Wostl E, Riyanto A, Hamidy A, Kurniawan N, Smith EN, Harvey MB. A taxonomic revision of the Philautus (Anura: Rhacophoridae) of Sumatra with the description of four new species. Herpetol Monogr. 2017;31:70–113. https://doi.org/10.1655/HERPMONOGRAPHS-D-16-00007.

  15. 15.

    Quah ESH, Anuar S, Grismer LL, Wood PL, Siti Azizah MNS, Muin MA. A new species of frog of the genus Abavorana Oliver, Prendini, Kraus & Raxworthy 2015 (Anura: Ranidae) from Gunung Jerai, Kedah, northwestern peninsular Malaysia. Zootaxa. 2017;4320:272–88.

  16. 16.

    Matsui M, Belabut DM, Ahmad N. Two new species of fanged frogs from peninsular Malaysia (Anura: Dicroglossidae). Zootaxa. 2014;3881:75–93.

  17. 17.

    Dehling JM, Dehling DM. A new wide-headed fanged frog of the Limnonectes kuhlii group (Anura: Dicroglossidae) from western Borneo with a redescription of Rana conspicillata Günther, 1872. Zootaxa. 2017;4317:291–309.

  18. 18.

    Matsui M, Nishikawa K. Description of a new species of Limnonectes from Sarawak, Malaysian Borneo (Dicroglossidae, Anura). Curr Herpetol. 2014;33:135–47. https://doi.org/10.5358/hsj.33.135.

  19. 19.

    Matsui M, Nishikawa K, Eto K. A new burrow-utilizing fanged frog from Sarawak, Malaysian Borneo (Anura, Dicroglossidae, Limnonectes). Raffles Bull Zool. 2014;62:679–87.

  20. 20.

    Davis HR, Grismer LL, Klabacka RL, Muin MA, Quah ESH, Anuar S, et al. The phylogenetic relationships of a new stream toad of the genus Ansonia Stoliczka, 1870 (Anura: Bufonidae) from a montane region in peninsular Malaysia. Zootaxa. 2016;4103:137–53.

  21. 21.

    Chan KO, Wood PLJ, Anuar S, Muin MA, Quah ESH, Sumarli AXY, et al. A new species of upland stream toad of the genus Ansonia Stoliczka, 1870 (Anura: Bufonidae) from northeastern peninsular Malaysia. Zootaxa. 2014;3764:427–40.

  22. 22.

    Chan KO, Alexander AM, Grismer LL, Su Y-C, Grismer JL, Quah ESH, et al. Species delimitation with gene flow: a methodological comparison and population genomics approach to elucidate cryptic species boundaries in Malaysian torrent frogs. Mol Ecol. 2017;26:5435–50.

  23. 23.

    Dever JA. A new cryptic species of the Theloderma asperum Complex (Anura: Rhacophoridae) from Myanmar. J Herpetol. 2017;51:425–36. https://doi.org/10.1670/17-026.

  24. 24.

    Poyarkov NA, Orlov NL, Moiseeva AV, Galoyan EA, Nguyen TT, Gogoleva SS. Sorting out moss frogs: mtDNA data on taxonomic diversity and phylogenetic relationships of the Indochinese species of the genus Theloderma (Anura, Rhacophoridae). Russ J Herpetol. 2015;22 December:241–80.

  25. 25.

    Meegaskumbura M, Meegaskumbura S, Bowatte G, Manamendra-Arachchi K, Pethiyagoda R, Hanken J, et al. Taruga (Anura: Rhacophoridae), a new genus of foam-nesting tree frogs endemic to Sri Lanka. Ceylon J Sci (Biological Sci). 2011;39:75–94.

  26. 26.

    Matsui M, Yambun P, Sudin A. Taxonomic relationships of Ansonia anotis Inger, tan, and Yambun, 2001 and Pedostibes maculatus (Mocquard, 1890), with a description of a new genus (Amphibia, Bufonidae). Zool Sci. 2007;24:1159–66.

  27. 27.

    Frost DR. Amphibian Species of the World: an Online Reference. Version 6.0 (accessed 28 Aug 2018). Electronic Database accessible at http://research.amnh.org/herpetology/amphibia/index.html. American Museum of Natural History, New York, USA. 2018.

  28. 28.

    Grismer LL, Wood PL Jr, Aowphol A, Cota M, Grismer MS, Murdoch ML, et al. Out of Borneo, again and again: biogeography of the Stream Toad genus Ansonia Stoliczka (Anura: Bufonidae) and the discovery of the first limestone cave-dwelling species. Biol J Linn Soc. 2017;120:371–95.

  29. 29.

    Matsui M, Tominaga A, Liu W, Khonsue W, Grismer LL, Diesmos AC, et al. Phylogenetic relationships of Ansonia from Southeast Asia inferred from mitochondrial DNA sequences: systematic and biogeographic implications (Anura: Bufonidae). Mol Phylogenet Evol. 2010;54:561–70. https://doi.org/10.1016/j.ympev.2009.08.003.

  30. 30.

    Matsui M, Nishikawa K, Eto K, Hossman MY. A new species of Pelophryne from Western Sarawak, Malaysian Borneo (Anura, Bufonidae). Zool Sci. 2017;34:345–50.

  31. 31.

    McLeod DS. Of least concern? Systematics of a cryptic species complex: Limnonectes kuhlii (Amphibia: Anura: Dicroglossidae). Mol Phylogenet Evol. 2010;56:991–1000. https://doi.org/10.1016/j.ympev.2010.04.004.

  32. 32.

    Peng Z, Pan T, Wan Y, Qian L, Wu J, Zhang B. Phylogeny of Asian Bufonids inferred from mitochondrial DNA sequences (Anura: Amphibia): implication for the speciation of east Asian Bufonids. Mitochondrial DNA Part A DNA Mapping, Seq Anal. 2017;28:358–60.

  33. 33.

    Wogan GOU, Stuart BL, Iskandar DT, Mcguire JA. Deep genetic structure and ecological divergence in a widespread human commensal toad. Biol Lett. 2016;12:20150807.

  34. 34.

    Inger RF, Stuart BL, Iskandar DT. Systematics of a widespread southeast Asian frog, Rana chalconota (Amphibia: Anura: Ranidae). Zool J Linnean Soc. 2009;155:123–47.

  35. 35.

    Liu Z, Chen G, Zhu T, Zeng Z, Lyu Z, Wang J, et al. Prevalence of cryptic species in morphologically uniform taxa – fast speciation and evolutionary radiation in Asian toads. Mol Phylogenet Evol. 2018;127 June:723–31. https://doi.org/10.1016/j.ympev.2018.06.020.

  36. 36.

    Dever JA, Fuiten AM, Konu Ö, Wilkinson JA. Cryptic torrent frogs of Myanmar: an examination of the Amolops marmoratus species complex with the resurrection of Amolops afghanus and the identification of a new species. Copeia. 2012;2012:57–76.

  37. 37.

    Stuart BL, Inger RF, Voris HK. High level of cryptic species diversity revealed by sympatric lineages of southeast Asian forest frogs. Biol Lett. 2006;2:470–4. https://doi.org/10.1098/rsbl.2006.0505.

  38. 38.

    Brown RM, Linkem CW, Siler CD, Sukumaran J, Esselstyn JA, Diesmos AC, et al. Phylogeography and historical demography of Polypedates leucomystax in the islands of Indonesia and the Philippines: evidence for recent human-mediated range expansion? Mol Phylogenet Evol. 2010;57:598–619. https://doi.org/10.1016/j.ympev.2010.06.015.

  39. 39.

    Lu B, Bi K, Fu J. A phylogeographic evaluation of the Amolops mantzorum species group: cryptic species and plateau uplift. Mol Phylogenet Evol. 2014;73:40–52.

  40. 40.

    Camargo A, Sites JJ. Species delimitation: A decade after the renaissance. species Probl - ongoing issues; 2013. p. 225–47. https://doi.org/10.5772/52664.

  41. 41.

    Luo A, Ling C, Ho SYW, Zhu C-D. Comparison of methods for molecular species delimitation across a range of speciation scenarios. Syst Biol. 2018;67:830–46.

  42. 42.

    Sukumaran J, Knowles LL. Multispecies coalescent delimits structure, not species. Proc Natl Acad Sci. 2017;114:1607–12.

  43. 43.

    Jackson ND, Carstens BC, Morales AE, O’Meara BC. Species delimitation with gene flow. Syst Biol. 2017;66:799–812.

  44. 44.

    Leaché AD, Zhu T, Rannala B, Yang Z. The spectre of too many species. Syst Biol. 2018;0:1–14. https://doi.org/10.1093/sysbio/syy051.

  45. 45.

    Fouquet A, Gilles A, Vences M, Marty C, Blanc M, Gemmell NJ. Underestimation of species richness in neotropical frogs revealed by mtDNA analyses. PLos One. 2007;2:e1109. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0001109.

  46. 46.

    Vences M, Thomas M, Bonett RM, Vieites DR. Deciphering amphibian diversity through DNA barcoding: chances and challenges. Philos Trans R Soc Lond Ser B Biol Sci. 2005;360:1859–68. https://doi.org/10.1098/rstb.2005.1717.

  47. 47.

    Yang Z, Rannala B. Bayesian species delimitation using multilocus sequence data. Proc Natl Acad Sci U S A. 2010;107:9264–9. https://doi.org/10.1073/pnas.0913022107.

  48. 48.

    Portik DM, Papenfuss TJ. Historical biogeography resolves the origins of endemic Arabian toad lineages (Anura: Bufonidae): evidence for ancient vicariance and dispersal events with the horn of Africa and South Asia. BMC Evol Biol. 2015;15:1–19. https://doi.org/10.1186/s12862-015-0417-y.

  49. 49.

    Pyron AR, Wiens JJ. A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Mol Phylogenet Evol. 2011;61:543–83.

  50. 50.

    Eaton DAR, Hipp AL, González-Rodríguez A, Cavender-Bares J. Historical introgression among the American live oaks and the comparative nature of tests for introgression. Evolution (N Y). 2015;69:2587–601.

  51. 51.

    Léveillé-Bourret É, Starr JR, Ford BA, Moriarty Lemmon E, Lemmon AR. Resolving rapid radiations within angiosperm families using anchored phylogenomics. Syst Biol. 2018;67:94–112.

  52. 52.

    Meiklejohn KA, Faircloth BC, Glenn TC, Kimball RT, Braun EL. Analysis of a rapid evolutionary radiation using ultraconserved elements: evidence for a bias in some multispecies coalescent methods. Syst Biol. 2016;65:612–27.

  53. 53.

    Tarver JE, Dos Reis M, Mirarab S, Moran RJ, Parker S, O’Reilly JE, et al. The interrelationships of placental mammals and the limits of phylogenetic inference. Genome Biol Evol. 2016;8:330–44.

  54. 54.

    Whitfield JB, Lockhart PJ. Deciphering ancient rapid radiations. Trends Ecol Evol. 2007;22:258–65.

  55. 55.

    Roux C, Fraïsse C, Romiguier J, Anciaux Y, Galtier N, Bierne N. Shedding light on the Grey zone of speciation along a continuum of genomic divergence. PLoS Biol. 2016;14:1–22.

  56. 56.

    Barley AJ, White J, Diesmos AC, Brown RM. The challenge of species delimitation at the extremes: diversification without morphological change in Philippine sun skinks. Evolution (N Y). 2013;67:3556–72. https://doi.org/10.1111/evo.12219.

  57. 57.

    Vences M, Thomas M, van der Meijden A, Chiari Y, Vieites DR. Comparative performance of the 16S rRNA gene in DNA barcoding of amphibians. Front Zool. 2005;2:5.

  58. 58.

    Vieites DR, Wollenberg KC, Andreone F, Kohler J, Glaw F, Vences M. Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proc Natl Acad Sci. 2009;106:8267–72. https://doi.org/10.1073/pnas.0810821106.

  59. 59.

    Evans BJ, Brown RM, JA MG, Supriatna J, Andayani N, Diesmos A, et al. Phylogenetics of fanged frogs: testing biogeographical hypotheses at the interface of the Asian and Australian faunal zones. Syst Biol. 2003;52:794–819.

  60. 60.

    Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28:1647–9. https://doi.org/10.1093/bioinformatics/bts199.

  61. 61.

    Nguyen LT, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74. https://doi.org/10.1093/molbev/msu300.

  62. 62.

    Chernomor O, Von Haeseler A, Minh BQ. Terrace aware data structure for Phylogenomic inference from Supermatrices. Syst Biol. 2016;65:997–1008.

  63. 63.

    Kalyaanamoorthy S, Minh BQ, Wong TKF, Von Haeseler A, Jermiin LS. ModelFinder: Fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9. https://www.nature.com/articles/nmeth.4285.

  64. 64.

    Lanfear R, Calcott B, Kainer D, Mayer C, Stamatakis A. Selecting optimal partitioning schemes for phylogenomic datasets. BMC Evol Biol. 2014;14:1–14.

  65. 65.

    Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Le SV. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol. 2017;35:518–22. https://doi.org/10.1093/molbev/msx281.

  66. 66.

    Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:1–6.

  67. 67.

    Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. 2010 Gatew Comput Environ Work GCE 2010. 2010;:1–8.

  68. 68.

    Bouckaert R, Drummond A. bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC Evol Biol. 2017;17:1–11. https://doi.org/10.1101/020792.

  69. 69.

    Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67:901–904. https://academic.oup.com/sysbio/article/67/5/901/4989127.

  70. 70.

    Padial JM, De La Riva I. Integrative taxonomists should use and produce DNA barcodes. Zootaxa. 2007.

  71. 71.

    Vieites DR, Min M-S, Wake DB. Rapid diversification and dispersal during periods of global warming by plethodontid salamanders. Proc Natl Acad Sci U S A. 2007;104:19903–7.

  72. 72.

    Swofford DL. PAUP*. Phylogenetic analysis using parsimony (*and other methods). Sunderland: Sinauer Associates; 2002.

  73. 73.

    Flouri T, Jiao X, Rannala B, Yang Z, Yoder A. Species tree inference with BPP using genomic sequences and the multispecies coalescent. Mol Biol Evol. 2018:1–9. https://doi.org/10.1093/molbev/msy147.

  74. 74.

    Yang Z. A tutorial of BPP for species tree estimation and species delimitation. Curr Zool. 2015;61:854–65.

  75. 75.

    Pinho C, Hey J. Divergence with gene flow: models and data. Annu Rev Ecol Evol Syst. 2010;41:215–30. https://doi.org/10.1146/annurev-ecolsys-102209-144644.

Download references

Acknowledgements

We thank the two anonymous reviewers who provided valuable comments on the manuscript.

Funding

KOC’s molecular laboratory work was supported by U.S. National Science Foundation (DEB 1702036) and the National Geographic Society (9722–15). These funding agencies played no role in the design of the study and collection, analysis, interpretation of data, and writing of the manuscript.

Availability of data and materials

All data generated or analysed during this study are included in this published article [and its supplementary information files].

Author information

KOC designed the project, performed the research, and wrote the manuscript; LLG provided genetic samples, manuscript editing, and project conceptualization. All authors have read and approved the manuscript.

Correspondence to Kin Onn Chan.

Ethics declarations

Ethics approval and consent to participate

No animals/specimens were used in this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Maximum likelihood phylogeny derived from 6236 bp comprising three mitochondrial (12S, 16S, CO1) and three nuclear genes (CXCR4, NCX1, RAG-1). Node values denote Ultrafast Bootstrap support values. (PDF 8 kb)

Additional file 2:

Bayesian phylogeny derived from 6236 bp comprising three mitochondrial (12S, 16S, CO1) and three nuclear genes (CXCR4, NCX1, RAG-1). Node values denote Bayesian posterior probabilities. (PDF 7 kb)

Additional file 3:

Posterior distributions of tau generated from priors versus empirical data. (PDF 3012 kb)

Additional file 4:

Molecular samples used in this study and their associated GenBank accession number. The asterisk (*) denotes newly sequenced samples for this study. (XLSX 11 kb)

Additional file 5:

Pairwise uncorrected p-distances of the 16S gene. (CSV 6 kb)

Additional file 6:

Results of the BPP A11 analysis. (ZIP 216 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • BPP
  • Genealogical divergence index
  • Gdi
  • Systematics
  • Species boundaries