Print
A taxonomic revision of the genus Phrynoglossus in Indochina with the description of a new species and comments on the classification within Occidozyginae (Amphibia, Anura, Dicroglossidae)
expand article infoGunther Köhler, Joseph Vargas, Ni Lar Than§, Tilman Schell|, Axel Janke|, Steffen U. Pauls|#, Panupong Thammachoti¤
‡ Senckenberg Forschungsinstitut und Naturmuseum, Frankfurt am Main, Germany
§ East Yangon University, Yangon, Myanmar
| LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany
¶ Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
# Justus-Liebig-University, Gießen, Germany
¤ Chulalongkorn University, Bangkok, Thailand
Open Access

Abstract

Abstract

We revise the frogs of the genus Phrynoglossus from Indochina based on data of external morphology, bioacoustics and molecular genetics. The results of this integrative study provide evidence for the recognition of three distinct species, one of which we describe as new. Phrynoglossus martensii has a vast geographic distribution from central and southern Thailand across southern China to Vietnam, Laos, and Cambodia. Phrynoglossus myanhessei sp. nov. is distributed in central Myanmar whereas Phrynoglossus magnapustulosus is restricted to the Khorat Plateau, Thailand. These three species occur in allopatry and differ in their mating calls, external morphology, and in genetic distances of the 16S gene of 3.8–5.9%. Finally, we discuss and provide evolutionary evidence for the recognition of Phrynoglossus as a genus distinct from Occidozyga. Members of both genera form reciprocal monophyletic groups in our analyses of mtDNA data (16S) and are well differentiated from each other in morphology and ecology. Furthermore, they differ in the amplexus mode with Phrynoglossus having an inguinal amplexus whereas it is axillary in Occidozyga. We further provide a de novo draft genome of the holotype based on short-read sequencing technology to a coverage of 25-fold. This resource will permanently link the genetic characterization of the species to the name-bearing type specimen.

Key words

amplexus mode, bioacoustics, cryptic species diversity, Dicroglossidae, genome, Khorat Plateau, Myanmar, new species, Occidozyga, Phrynoglossus myanhessei sp. n., Thailand

Introduction

Southeast Asia is recognized as a major global biodiversity hotspot (Myers et al. 2000; Mittermeier 2004; Corlett 2014), and the amphibian order Anura is well represented in this region. However, it seems that the actual diversity of anurans in this part of the world is still grossly underestimated due to poor sampling of areas with difficult or dangerous access and because most of the supposedly wide-spread species actually constitute species complexes of two or more distinct species (e.g., Funk et al. 2011; Hasan et al. 2012; Köhler et al. 2019). During field work in Myanmar and Thailand we came across variation in the advertisement call of the common mud or puddle frogs, in this region all currently referred to as Occidozyga martensii Peters, 1867 (Niyomwan et al. 2019; Frost 2020). These observations prompted us to collect this group of frogs to obtain a reasonable geographic coverage with the goal to evaluate whether these really are a single species occurring in a large geographic area or if these actually represent a complex of several – morphologically similar (cryptic) – species with smaller geographic ranges.

The generic placement of martensii in either Phrynoglossus Peters, 1867 or Occidozyga Kuhl and van Hasselt, 1822 has caused some controversy among recent authors. It seems that the majority of authors have treated Phrynoglossus as a synonym of Occidozyga (e.g., Duellman and Trueb 1994; Nguyen et al. 2009; Frost 2020) whereas others have recognized them as two valid genera (e.g., Manthey and Grossmann 1997; Ohler and Dubois 1999; Ziegler, 2002; Köhler et al. 2018). Here we provide a short summary of the taxonomic history of Phrynoglossus and Occidozyga:

The genus name Occidozyga was introduced by Kuhl and van Hasselt 1822. Its type species (by subsequent designation of Stejneger 1925) was described as Rana lima Gravenhorst 1829 a few years later. In 1859, Günther described Oxyglossus laevis based on material from the Philippines. Oxyglossus Tschudi 1838 is an objective synonym of Occidozyga Kuhl and van Hasselt 1822 according to Ohler and Dubois (1999) whereas Ooeidozyga, Oxydozyga, Occidogyna are incorrect subsequent spellings of Occidozyga (Dubois 1981, Frost 2020). The genus name Oreobatrachus Boulenger 1896 is a synonym of Phrynoglossus Peters 1867 according to Smith (1931).

In 1896, Boulenger introduced his new species Oreobatrachus baluensis from „Mount Kina Balu, North Borneo [= Sabah]“, Malaysia (Borneo). Smith (1931) placed Oreobatrachus Boulenger 1896 in the synonymy of Phrynoglossus Peters 1867. With the introduction of the latter genus, Peters (1867) described the new species P. martensii, the holoype of which he had received from Bangkok, Thailand. In 1877, Peters described Microdiscopus sumatranus from Sumatra, Indonesia. Smith (1916) discussed the taxa laevis, lima and martensii and considered martensii to be a subspecies of laevis. Also, he described and illustrated the tadpoles of lima and martensii and provided natural history notes for the two taxa. He stated (Smith 1916: 175) that “O. lima is strictly aquatic in its habits” whereas martensii, “although never far from water, is seldom to be found in it.” Taylor (1922) introduced the new species Micrixalus diminutiva based on specimens from “near Pasananka, Zamboanga, Mindanao”, Philippines. Inger (1954) transferred the latter taxon to the genus Ooeidozyga. In 1927, Smith described two new species from Sulawesi, Indonesia: Ooeidozyga semipalmata from „Lowah, near Mt. Bonthain“, and Ooeidozyga celebensis from “Djikoro, Mt. Bonthain”. Based on specimens from “Rana Mese, 1200 m. H., West-Flores”, Indonesia, Mertens (1927) described Oxydozyga floresiana. In 1958, Taylor and Elbel introduced their new species Occidozyga magnapustulosa with type locality in “Ban Na Phua (subvillage), Kan Luang (village), Na Kae (district), Nakhon Phanom (province), Thailand, elevation approx. 200 m”. More recently, Iskandar et al. (2011) analyzed the morphological variation in populations of frogs related to Occidozyga semipalmata and recognized the populations from „Bantayan, Mount Tompotika, Balantak Mountains …, Desa (=Village) Bualemo, … Sulawesi Province, Indonesia“ as a distinct species that they named Occidozyga tompotika.

In his treatment of the herpetofauna of Mount Kinabalu, Borneo, Smith (1931) pointed out the unique tongue morphology of Ooeidozyga lima, strikingly different from the other species of that genus, and proposed to retain only O. lima in the genus and to transfer the remaining species to the genus Phrynoglossus. However, this proposal was not widely accepted by subsequent authors (e.g., Inger 1954, 1996) who preferred to maintain all species in a single genus (i.e., Occidozyga respectively its synonym Ooeidozyga). Inger (1996) argued that of the five supposedly diagnostic characters listed for Phrynoglossus, only one, the tongue morphology, would distinguish the species of the latter genus from Occidozyga.

Further evidence for considering Phrynoglossus as a synonym of Occidozyga was provided in the large-scale molecular phylogeny of amphibians of Pyron and Wiens (2011) where the two species of Phrynoglossus were nested within the species of Occidozyga, thereby rendering the latter genus paraphyletic if Phrynoglossus was recognized as a valid genus. Basal to the clade containing the species of Phrynoglossus and Occidozyga is Ingerana in the tree of these authors. However, aside from tongue morphology, there is additional evidence in favor of recognizing these two groups as distinct genera (see our Results section), and we therefore recognize Phrynoglossus as a valid genus distinct from Occidozyga.

Materials and methods

Specimens examined for this study were personally collected by GK, NLT, and PT (see Appendix 1 for specimens examined). Specimens labeled with GK field numbers were deposited in the collections of Senckenberg Forschungsmuseum Frankfurt (SMF) or at East Yangon University (EYU), Thanlyin, Myanmar, and those with PT field numbers were deposited in the collection of the Chulalongkorn University, Museum of Natural History (CUMZ), Bangkok.

Prior to preservation of collected specimens in the field, we took color photographs of each individual in life. We euthanized the frogs with a pericardial injection of T61. We cut tissue samples from one forelimb or from the tongue and preserved these in 98% non-denatured ethanol for DNA extraction. The tissue samples were deposited in the collection of SMF and CUMZ. Specimens were then preserved by injecting a solution of 5 mL absolute (i.e., 36%) formalin in 1 L of 96% ethanol into the body cavity, and stored in 70% ethanol. Coordinates and elevation were recorded using Garmin GPS receivers with built-in altimeters. All coordinates were recorded in decimal degrees, WGS 1984 datum. Capitalized colors and color codes (in parentheses) followed Köhler (2012).

In evaluating species’ boundaries within and among populations, we followed the evolutionary species concept (Simpson 1951; Wiley 1978). As lines of evidence for species delimitation, we applied a phenotypic criterion (external morphology), the genetic distinctness of a mitochondrial genetic marker as well as a criterion for reproductive isolation (bioacoustic data).

Abbreviations used are EYD (eye diameter); FL (foot length); HL (head length); HW (head width); IND (internasal diameter); IOD (interorbital diameter); NED (nostril–eye distance); HNL (hand length); SHL (shank length); SL (snout length); SVL (snout–vent length); TED (tympanum–eye distance); THL (thigh length); TYD (longitudinal tympanum diameter). Webbing formulae follow (Savage and Heyer 1997). Terminology of snout shape follows Heyer et al. (1990).

We recorded anuran vocalizations using a digital audio recorder (Olympus LS-12) with a Sennheiser ME 66 shotgun microphone capsule and a Sennheiser K6 powering module. The microphone was positioned between 0.5 and 1.5 m from the calling frog. Aside from the GPS coordinates and elevation above sea level of the locality we also noted ambient air temperature and determined SVL of the recorded individual. Files were recorded as uncompressed 24-bit WAV files at a sampling frequency of 96 kHz. Audio files were deposited in the Fonoteca Zoológica, Museo Nacional de Ciencias Naturales, Madrid, Spain.

The spectral and temporal parameters were analyzed and the power spectra were calculated in RAVEN PRO 1.4. (Bioacoustics Research Program 2011). Spectrograms were obtained using the Blackman window function at 3db Filter Bandwidth of 141 Hz; grid spacing of 21.5 Hz; overlap 90%. Frequency information was generated through Fast Fourier Transformation (FFT, width 2,048 samples). Temporal measurements of calls such as repetition rates, duration of notes, and number of pulses, were measured on the waveforms. Terminology in call descriptions follows Köhler et al. (2017). The map was created using ArcMap 10.4. Additionally to the specimens examined by us, we also plotted specimens from the Field Museum (FMNH), Chicago, not examined by us.

Marker based analysis

We extracted DNA following the protocol of Ivanova et al. (2006). To eliminate potential PCR-inhibiting contaminants, the tissue samples were incubated for 14 hours at 4°C in 200 µL low PBS buffer (20 µL PBS in 180 µL of water) before overnight digestion with the vertebrate lysis buffer at 56°C. After extraction, DNA was eluted in 50 µL TE buffer. Fragments of the mitochondrial 16S rRNA (16S) were amplified in an Eppendorf Mastercycler® Pro using the following protocol: initial denaturation for 2 min at 94°C; followed by 40 cycles with denaturation for 35 s at 94°C, hybridization for 35 s at 48.5°C, and elongation for 60 s at 72°C; final elongation for 10 min at 72°C. The reaction mix for each sample contained 1 µL DNA template, 14 µL water, 2.5 µL PCR-buffer, 1 µL 25 mM MgCl2, 4 µL 2.5 mM dNTPs (Invitrogen), 0.5 µL (containing 2.5 units) Taq Polymerase (PeqLab), and 1 µL of the primer (16S: forward: L2510, 5’–CGCCTGTTTATCAAAAACAT–3’; reverse: H3056, 5’–CCGGTCTGAACTCAGATCACGT–3’; from Eurofins MWG Operon). DNA extraction, PCR, and sequencing were done at SMF for the samples from Myanmar and at Chulalongkorn University for the samples from Thailand. We generated 35 new sequences for this study (see Appendix 2). Additionally, we downloaded relevant 16S sequences from GenBank (Appendix 2). Because the frogs of the genus Ingerana Dubois, 1987 are supposedly the closest relatives of Phrynoglossus and Occidozyga (see Pyron and Wiens 2011) we included sequences of Ingerana tenasserimensis (Sclater 1892) in our analyses. Our dataset contains the type species of Ingerana (i.e., Rana tenasserimensis Sclater, 1892), Occidozyga (i.e., Rana lima Gravenhorst, 1829), and Phrynoglossus (i.e., Phrynoglossus martensii Peters, 1867).

We aligned the sequences with MUSCLE (Edgar 2004) using the default settings in Geneious 6.1.2. (Kearse et al. 2012). For software applications, sequence data formatting was converted using the online server Alter (Glez-Peña et al. 2010). The best substitution model for each gene was identified using PartitionFinder2 (Lanfear et al. 2017), with linked branch lengths via PhyML 3.0 (Guindon et al. 2010). Model selection used the corrected (for finite sample size) Akaike Information Criterion (AICc) (Burnham and Anderson 2002). Limnonectes limborgi (GK-7110) was used as outgroup (Pyron and Wiens 2011).

Bayesian Inference analyses (BI) used MrBayes 3.2 (Ronquist et al. 2012) with five runs with eight chains. The first 25% of trees were discarded as burn-in. MCMC runs used an initial set of 1,000,000 generations with sampling every 500 generations, and adding 500,000 generations until chains reached convergence. Convergence was considered achieved when the standard deviation of split frequencies was 0.015 or less. Additionally, convergence was diagnosed by PRSF (Potential Scale Reduction Factor), which should approach 1.0 as runs converge (Gelman and Rubin 1992). We used the IQTree webserver (Trifinopoulos et al. 2016) to run a Maximum Likelihood (ML) analysis using 10,000 ultrafast Bootstrap approximation (UFBoot) replicates with 10,000 maximum iterations and minimum correlation coefficient of 0.99 (Minh et al. 2013), plus 10,000 replicates of Shimodaira-Hasegawa approximate likelihood ratio (SH-aLRT), which proved to be accurate with a high statistical power (Guindon et al. 2010). We used FigTree 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree) for viewing the trees. We estimated evolutionary genetic divergence, computing uncorrected pairwise distances with MEGA 7.0.26 (Kumar et al. 2016) to assess the degree of intra- and interspecific differences, using a Bootstrap estimation method of 10,000 replications. We built a species tree using BEAST 2.4.7 (Ogilvie et al. 2017; Bouckaert et al. 2018) with 1,000,000 generations for the MCMC model. The resulting tree was visualized using DensiTree 2.2.6 (Bouckaert and Heled 2014). Automatic barcode gap discovery (ABGD) (Puillandre et al. 2012) was run through its webserver (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html), with using Simple Distance and default values for Prior Intraspecific divergence, except for relative gap width (1.5), which did not work for our genes (see also Kekkonen et al. 2015). Because high values in relative gap-width tend to overly split species (Yang et al. 2016), we used an intermediate value of 1.0.

Whole genome sequencing and assembly

De novo whole genome sequencing and assembly began with genomic DNA extraction from muscle tissue of the holotype according to standard phenol/chloroform procedures (Sambrook and Russell 2001). The resulting DNA was resuspended in TE buffer (10mM Tris-Cl, 0.1 mM EDTA) and stored at –20°C. Quality checks for high molecular weight DNA were performed by agarose gel electrophoresis (Sambrook and Russell 2001). The DNA sample was shipped on dry ice to Novogene (UK) for short read Illumina genome sequencing. One genomic library (insert size: 350 bp) was prepared and 150 bp paired-end reads were sequenced on an Illumina NovaSeq 6000 platform (San Diego, CA). Raw reads are deposited under the accession number SRR13288470.

A k-mer profile was generated from the raw reads using Jellyfish 2.3.0 (Marçais and Kingsford 2011) and analyzed in the GenomeScope webserver (Vurture et al. 2017). Raw reads were trimmed for low quality regions and adapter sequences and filtered for possible contamination using Trimmomatic 0.39 (Bolger et al. 2014) and Kraken 2.0.9 (Wood et al. 2019) respectively. Best k-mer length was estimated by KmerGenie 1.7051 (Chikhi and Medvedev 2014). Nuclear genome assembly was conducted in Velvet 1.2.10 (Zerbino and Birney 2008) without reporting contigs smaller than 500 bp. The mitochondrial genome was assembled from the genomic data separately using NOVOplasty 4.2 (Dierckxsens et al. 2017). Annotation of the mitochondrial genome was manually merged and curated in Geneious Prime 2020.2.3 (https://www.geneious.com) from automatic annotations of GeSeq (Tillich et al. 2017) and MITOS2 (Donath et al. 2019) webservers. Scaffolds matching the mitochondrial genome assembly and scaffolds smaller than 500 bp were removed from the nuclear genome assembly. Finally, the mitochondrial genome assembly was added and quality was checked by backmapping the reads to the assembly with bwa mem 0.7.17-r1188 (Li 2013), searching for possible contamination with Blobtools 1.1.1 (Laetsch and Blaxter 2017) and screening for single-copy orthologs with BUSCO 4.1.4 (Simão et al. 2015). For detailed description see Appendix 3.

Results and Conclusions

The final 16S alignment for 69 samples (genera Ingerana, Occidozyga, and Phrynoglossus) was 599 nucleotide positions long. GTR+G was determined to be the best fitting model of sequence evolution. The trees obtained by BI, ML, *BEAST, and ABGD analyses show a high degree of congruence at well-supported nodes, with some differences in branch arrangement at poorly supported nodes (Figs 1, 2 and 3). Our final phylogenetic analyses recover the deep divergences between the genera Ingerana, Occidozyga, and Phrynoglossus (Figs 1, 2 and 3). These three genera form monophyletic groups in all analyses.

Figure 1. 

Phylogenetic tree of frogs of the subfamily Occidozyginae from a Bayesian Inference analysis of the mitochondrial marker 16S. A scale bar of genetic distance is indicated. Numbers at nodes are bootstrap values (left) and Bayesian posterior probabilities (right), but only for nodes with bootstrap values higher than 75. The tree is rooted using Limnonectes limborgi. — Abbreviations: P. = Phrynoglossus; myanh. = myanhessei sp. n.; magn. = magnapustulosus; I. ten. = Ingerana tenasserimensis.

Ingerana is the sister clade to a clade containing Occidozyga plus Phrynoglossus. There is much genetic structure in the Occidozyga clade indicating the possible presence of cryptic species among the populations currently referred to as O. lima (Gravenhorst, 1829). Among the specimens here referred to the genus Phrynoglossus, four major clades are recognizable and also supported by the ABGD analysis (Fig. 3). Clade 1 (= Group 7 in our ABGD analysis) contains specimens from southern Thailand including Bangkok (= type locality of P. martensii) and are therefore considered to represent the “true” P. martensii. Clade 2 (= Group 8 in our ABGD analysis) contains specimens from northern Thailand, southern China, Vietnam, Laos, and Cambodia. Clade 3 (= Group 11 in our ABGD analysis) is geographically restricted to the Khorat Plateau, Thailand, and includes specimens from Ban Kan Luang (= type locality of P. magnapustulosus Taylor and Elbel, 1958) and are therefore considered to represent the latter species. Finally, Clade 4 (= Group 10 in our ABGD analysis) contains specimens from central and lower Myanmar. We consider these four clades as candidate species. The genetic distances among these clades vary from 3.5% to 5.9% (Table 1).

Figure 2. 

Species tree inferred with ∗BEAST showing density of trees proportional to frequency of occurrence drawn in DensiTree. Abbreviations: P. = Phrynoglossus.

Figure 3. 

Species delimitation analysis of frogs of the subfamily Occidozyginae using the Automatic Barcode Gap Discovery (ABGD) approach; initial partition with prior maximal distance at 0.00167. See text for details. Abbreviations as in Fig. 1.

Table 1.

Uncorrected pairwise distances for the Occidozyginae plus Limnonectes (outgroup) included in this study. For details see text.

Ingerana Occidozyga Limnonectes martensii cf._martensii Myanmar
Occidozyga 0.189
Limnonectes 0.192 0.193
martensii 0.198 0.132 0.203
cf._martensii 0.202 0.125 0.213 0.035
„Myanmar“ 0.203 0.139 0.209 0.054 0.051
magnapustulosus 0.208 0.127 0.207 0.045 0.038 0.059

In bioacoustics, the advertisement calls of males from Clades 1, 3, and 4 can be readily distinguished whereas the calls of males from Clades 1 and 2 are exceedingly similar (Table 2; Fig. 4). Males of P. magnapustulosus (our Clade 3) emit an advertisement call that sounded “mourning” to our ears (“maaaaaaa”). It has a duration of 243–437 ms (mostly in the 300 ms range) with its dominant frequency mostly in the range of 3500–3700 Hz (3187–3790 Hz) and a gap duration between calls of 2.3–5.3 s, mostly 3–4 s. On the contrary, the call of P. martensii (Clade 1) is very short (sounding like “bic”) with a duration of 40–48 ms and a dominant frequency of 3467–3596 Hz. The gap duration between calls varies from 2.8–3.9 s. The calls of males from Myanmar (Clade 4) are somewhat similar to the calls of P. martensii but are of longer duration (82–114 ms) and lower pitched (dominant frequency 2454–2885 Hz). The gap duration between calls varies from 1.8–4.8 s. The advertisement call of our “cf. martensii” (Clade 2) is very similar to the call of P. martensii (Clade 1) and differs only in having mostly a shorter gap duration between calls (0.5–3.0 s, mostly < 2.5 s in “cf. martensii” versus 2.8–3.9 in P. martensii).

Figure 4. 

Advertisement calls of male Phrynoglossus. (A) P. magnapustulosus, GK-7395; (B) P. myanhessei n. sp., SMF 103799; (C) P. martensii, PT-2634 (Bangkok); (D) P. martensii, PT-2076 (Chiang Mai).

Table 2.

Selected bioacoustic parameters of the species related to Phrynoglossus martensii. Range is followed by mean value and standard deviation in parentheses. Dom. Freq. = dominant frequency; Freg. = frequency.

Taxon ID Call duration (sec) Dom. Freq (Hz) Gap duration (sec) Freq 5% (Hz) Freq 95% (Hz)
GK-7475_0544 (n=23 calls) Roi Et Prov., Thailand magnapustulosus 0.243–0.354 (0.291±0.024) 3187–3682 (3447±175) 2.25–4.28 (3.15±0.58) 1680–1938 (1790±74) 4199–5383 (4823±388)
GK-7465_0543 (n=15 calls) Roi Et Prov., Thailand magnapustulosus 0.322–0.437 (0.387±0.031) 3725–3790 (3761±20) 2.67–5.31 (3.92±0.67) 1744–2131 (1894±94) 4091–4177 (4150±29)
GK-7476_0546 (n=9 calls) Roi Et Prov., Thailand magnapustulosus 0.262–0.287 (0.277±0.008) 3750–3750 (3750±0) 2.53–3.20 (2.77±0.19) 1640–1734 (1708±32) 4406–4453 (4427±23)
GK-6727_0183 (n=7 calls) EYU, Myanmar n. sp. “Myanmar” 0.085–0.114 (0.104±0.009) 2454–2476 (2467±11) 2.34–2.90 (2.66±0.23) 1357–1809 (1529±156) 2820–2863 (2839±18)
GK-6926_0008 (n=21 calls) EYU, Myanmar n. sp. “Myanmar” 0.079–0.109 (0.098±0.008) 2756–2885 (2839±41) 2.47–4.78 (3.31±0.67) 1766–2347 (2119±168) 3187–3316 (3265±37)
Not collected (n=26 calls) Thanlyin, Myanmar n. sp. “Myanmar” 0.082–0.102 (0.090±0.005) 2541–2584 (2566±18) 1.84–2.80 (2.18±0.32) 1357–2110 (1831±238) 2929–3015 (2975±26)
Goutte_9357 (n=24 calls) Taiwan cf. martensii 0.041–0.051 (0.046±0.003) 3876–3941 (3914±17) 2.12–2.99 (2.39±0.32) 2304–3488 (3027±414) 4264–4737 (4493±128)
Goutte_9358 (n=28 calls) Taiwan cf. martensii 0.038–0.052 (0.043±0.003) 3919–3962 (3934±15) 1.94–2.78 (2.37±0.65) 2864–3553 (3356±414) 4264–4371 (4300±24)
Ziegler TZ-473 (n=24 calls) Vietnam cf. martensii 0.032–0.046 (0.037±0.004) 3338–3553 (3444±70) 0.53–0.90 (0.72±0.09) 1873–3144 (2874±275) 3725–4027 (3830±76)
Chiang Mai (n=10 calls) Thailand cf. martensii 0.033–0.043 (0.037±0.003) 3660–3704 (3689±14) 1.91–2.56 (2.19±0.20) 2670–3165 (2922±130) 4134–4435 (4254±111)
Bangkok (n=14 calls) Thailand martensii 0.040–0.048 (0.043±0.002) 3467–3596 (3544±37) 2.80–3.86 (3.42±0.34) 2261–3165 (2689±257) 3919–4155 (4028±69)

In external morphology, specimens of these four clades are somewhat conservative in external morphology but show differences in some morphometric values (Fig. 5; Table 3) as well as in the amount of toe webbing and in some details of coloration (see respective diagnosis sections below). Also, there is sexual dimorphism evident in several morphometric characters such as SVL, SHL, FL, and HL. Phrynoglossus magnapustulosus have the smallest SVL (males < 20 mm versus > 20 mm in the other clades) and also have a much larger ratio TYD/SVL compared to those from the other clades.

Figure 5 

part 1. Scatter plots illustrating morphological variation in the species of Phrynoglossus. magnap. = P. magnapustulosus; m = male; f = female. For abbreviations of morphological characters see text. Abbreviations of taxon names as in Fig. 1.

Figure 5 

part 2. Scatter plots illustrating morphological variation in the species of Phrynoglossus. magnap. = P. magnapustulosus; m = male; f = female. For abbreviations of morphological characters see text. Abbreviations of taxon names as in Fig. 1.

Table 3.

Selected measurements and proportions of the species of Phrynoglossus. Range is followed by mean value and standard deviation in parentheses. For abbreviations see text.

P. martensii ♂ 7 | ♀ 8 P. magnapustulosus ♂ 8 | ♀ 12 P. myanhessei ♂ 9 | ♀ 6
SVL males 19.77–23.31 (21.42±1.28) 17.21–19.68 (17.91±0.804) 22.00–27.35 (24.06±1.65)
females 25.55–30.35 (27.60±2.061) 16.90–24.30 (21.66±2.26) 28.38–31.27 (29.67±1.15)
SHL / SVL males 0.441–0.509 (0.479±0.021) 0.442–0.557 (0.509±0.036) 0.457–0.571 (0.506±0.032)
females 0.436–0.479 (0.462±0.016) 0.402–0.548 (0.478±0.041) 0.466–0.592 (0.500±0.047)
FL / SVL males 0.471–0.535 (0.491±0.024) 0.473–0.564 (0.527±0.031) 0.433–0.535 (0.476±0.037)
females 0.432–0.488 (0.468±0.021) 0.439–0.567 (0.482±0.045) 0.407–0.491 (0.449±0.033)
HL / SVL males 0.341–0.383 (0.359±0.015) 0.313–0.377 (0.341±0.024) 0.268–0.378 (0.330±0.031)
females 0.311–0.383 (0.344±0.022) 0.272–0.347 (0.304±0.028) 0.272–0.304 (0.289±0.011)
HW / SVL males 0.336–0.401 (0.357±0.023) 0.335–0.394 (0.364±0.021) 0.286–0.369 (0.338±0.029)
females 0.338–0.388 (0.351±0.017) 0.295–0.360 (0.336±0.019) 0.283–0.354 (0.312±0.025)
HL / HW males 0.909–1.077 (1.011±0.063) 0.855–1.030 (0.937±0.058) 0.897–1.191 (0.980±0.092)
females 0.920–1.020 (0.980±0.031) 0.779–1.057 (0.907±0.098) 0.808–1.043 (0.929±0.093)
IOD / SVL males 0.071–0.091 (0.086±0.007) 0.093–0.135 (0.112±0.015) 0.090–0.132 (0.107±0.013)
females 0.073–0.084 (0.077±0.003) 0.084–0.137 (0.106±0.018) 0.077–0.163 (0.104±0.032)
TYD / SVL males 0.042–0.068 (0.052±0.010) 0.055–0.099 (0.082±0.014) 0.029–0.066 (0.050±0.012)
females 0.035–0.052 (0.043±0.006) 0.062–0.102 (0.080±0.012) 0.034–0.064 (0.051±0.010)
EYD / SVL males 0.109–0.138 (0.122±0.010) 0.091–0.143 (0.116±0.017) 0.094–0.140 (0.112±0.015)
females 0.095–0.119 (0.104±0.008) 0.082–0.118 (0.100±0.012) 0.077–0.132 (0.096±0.019)
NED / SVL males 0.056–0.071 (0.064±0.006) 0.064–0.123 (0.086±0.018) 0.061–0.101 (0.083±0.016)
females 0.055–0.066 (0.060±0.004) 0.058–0.086 (0.074±0.009) 0.060–0.145 (0.090±0.033)
TYD / EYD males 0.357–0.512 (0.425±0.063) 0.482–0.920 (0.722±0.147) 0.291–0.653 (0.460±0.146)
females 0.328–0.515 (0.411±0.060) 0.591–0.962 (0.807±0.120) 0.385–0.727 (0.550±0.143)

Thus, based on the results of the analyses of mtDNA data, morphology, and bioacoustics, we recognize three of the four clades as defined above as species level units (i.e., our OTUs 1, 3 and 4) whereas we tentatively consider our OTU 2 to be conspecific with our OTU 1 since we are unable to find any diagnostic differences in morphology or bioacoustics that would separate these two OTUs. Based on the respective type localities, our Clade 1 is assigned to Phrynoglossus martensii whereas our Clade 3 is referred to as P. magnapustulosus. Clade 4 represents an undescribed species for which no name is available. Thus, we describe the Phrynoglossus from central and lower Myanmar as a new species below.

Regarding the recognition of Phrynoglossus as a valid genus distinct from Occidozyga, we find support by the results of our genetic and morphologic analyses as well as our field observations on the amplexus mode of these frogs. Phrynoglossus has an inguinal amplexus whereas Occidozyga has an axillary amplexus (Fig. 6). The two genera differ readily in tongue morphology with Occidozyga having a very slender, worm-like tongue whereas the tongue is fleshy and thick in Phrynoglossus (Fig. 7a,b) Also, species of the two groups differ strikingly in skin texture (Fig. 7c,d) and mucosome, at least judged from touching these frogs with bare hands: species of Phrynoglossus have a smooth skin and are extremely slimy and thus difficult to constrain manually. On the other hand, Occidozyga has a rough skin and is not slimy at all. Finally, the two groups differ in ecology: species of Phrynoglossus are terrestrial frogs that are usually found on mud at the shore of small water bodies whereas Occidozyga is a fully aquatic species that always remains in the water body. For the issue of possible paraphyly of this group of frogs see also our results of the analysis of the mtDNA data used in this study.

Figure 6. 

Frogs of Phrynoglossus and Occidozyga in amplexus. (A) P. myanhessei n. sp., SMF 103797–98; (B) P. magnapustulosus, not collected; (C) O. lima, not collected (Magwe State, Myanmar); (D) O. lima, not collected (Roi Et province, Thailand). Photos by G.K.

Figure 7. 

Morphological differences between Phrynoglossus and Occidozyga. Tongue morphology in (A) P. myanhessei n. sp., SMF 103353; (B) O. lima, GK-7076; dorsal skin texture in (C) P. myanhessei n. sp., SMF 103353; (D) O. lima, GK-7076. Photos by G.K.

Thus, we formally resurrect Phrynoglossus from the synonymy of Occidozyga and define the two genera as follows (data also from Ohler and Dubois 1999):

Occidozyga Kuhl & van Hasselt, 1822

Type species

Rana lima Gravenhorst, 1829.

Diagnosis

A genus of Asian frogs of the subfamily Occidozyginae Fei, Ye, and Huang, 1990 of the family Dicroglossidae Anderson, 1871, that differs from all other genera of its subfamily by having the following combination of characters: (1); tongue slender, worm-like; (2) vomerine teeth absent; (3) tips of fingers and toes pointed; (4) tympanum hidden; (5) skin not covered by extensive mucous, feels dry to touch in life frogs; (6) throat lining whitish with longitudinal brown stripe; (7) life style fully aquatic; (8) amplexus axillary.

Content

Occidozyga lima (Gravenhorst, 1829)

Phrynoglossus Peters, 1867

Type species

Phrynoglossus martensii Peters, 1867

Diagnosis

A genus of Asian frogs of the subfamily Occidozyginae Fei, Ye, and Huang, 1990 of the family Dicroglossidae Anderson, 1871, that differs from all other genera of its subfamily by having the following combination of characters: (1); tongue fleshy, swollen; (2) vomerine teeth absent; (3) tips of fingers and toes slightly swollen; (4) tympanum moderately distinct; (5) skin covered by extensive mucous, feels slimy to touch in life frogs; (6) throat lining uniformly grey; (7) life style semiaquatic; (8) amplexus inguinal.

Content

Phrynoglossus baluensis (Boulenger, 1896), Phrynoglossus celebensis (Smith, 1927), Phrynoglossus diminutivus (Taylor, 1922), Phrynoglossus floresianus (Mertens, 1927), Phrynoglossus laevis (Günther, 1859), Phrynoglossus magnapustulosus (Taylor and Elbel, 1958), Phrynoglossus martensii Peters, 1867, Phrynoglossus semipalmatus (Smith, 1927), Phrynoglossus sumatranus (Peters, 1877), Phrynoglossus tompotika (Iskandar, Arifin, and Rachmanasah, 2011).

This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the International Commission on Zoological Nomenclature (ICZN). The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information can be viewed through any standard web browser by appending the LSID to the prefix http://zoobank.org. The LSID for this publication is as follows: urn:lsid:zoobank.org:pub:A978D105-6D7F-4858-BA27-8639B76A6039. The LSID registration and any associated information can be viewed in a web browser by adding the LSID to the prefix “http://zoobank.org.”

Phrynoglossus myanhessei Köhler, Vargas, Than & Thammachoti, sp. nov.

Figs 8, 9, 10

Holotype

SMF 103841, an adult male collected at East Yangon University (16.77737, 96.24065; 17 m a.s.l.), Thanlyin, Yangon, Myanmar, collected 6 July 2017 by Gunther Köhler and Ni Lar Than. Field tag number GK-6728.

Paratypes

SMF 103840, same collecting data as holotype; SMF 103797–99, same locality data and collectors as holotype but collected 15 June 2018. All paratypes are adult males except for SMF 103798 which is an adult female.

Diagnosis

A species of the genus Phrynoglossus as defined above that differs from all mainland Southeast Asian congeners by having (1) a large body size (males 22.0–27.4 mm; females 28.4–31.3 mm); (2) relatively small tympanum (ratio TYD/SVL 0.48–0.92); and (3) call duration of male advertisement call 85–114 ms. Phrynoglossus myanhessei differs from its congeners in Indochina (i.e., P. martensii and P. magnapustulosus) in the male advertisement call, most obvious in call duration (85–114 ms in P. myanhessei, 32–52 ms in P. martensii, 243–437 ms in P. magnapustulosus) and dominant frequency (2454–2885 Hz in P. myanhessei, 3338–3962 Hz in P. martensii, 3187–3790 Hz in P. magnapustulosus).

Comparisons

Phrynoglossus myanhessei differs from its congeners as follows (condition for P. myanhessei in parentheses). Phrynoglossus floresianus, P. semipalmatus, and P. tompotika all have distinctly enlarged flattened toe and finger disks (tips of toes rounded, only slightly expanded into discs, not distinctly flattened); P. floresianus and P. laevis are much larger frogs, 35–37 mm in males, 40–51 mm in adult females in P. floresianus, 26–42 mm in males, 35–62 mm in females in P. laevis (males 22.0–27.4 mm; females 28.4–31.3 mm); furthermore, the eyes are directed dorsolaterally in P. laevis and P. celebensis (laterally); P. baluensis usually has a large inverted U-shaped ridge on the dorsum (no such ridge present) and reduced toe webbing with at least two phalanges free of webbing along fourth toe (feet almost fully webbed, less than one phalange free of webbing along fourth toe); P. sumatranus has a dark brown band on each side of the cloaca (such band absent), diamond shaped pupil (ovoid), and its eyes are oriented dorsolaterally (laterally); P. diminutivus has reduced toe webbing with two phalanges free of webbing along fourth toe (feet almost fully webbed, less than one phalange free of webbing along fourth toe); P. martensii and P. magnapustulosus are smaller frogs, 20–23 mm in males, 26–30 mm in adult females in P. martensii, 17–20 mm in males, 17–24 mm in females in P. magnapustulosus (males 22.0–27.4 mm; females 28.4–31.3 mm); furthermore, P. magnapustulosus has a relatively larger tympanum, ratio tympanum/SVL 0.055–0.099, mean 0.082 (0.029–0.066, mean 0.050). Also, P. martensii has a relatively broader head, ratio HW/SVL 0.336–0.401, mean 0.357 in males, 0.338–0.388, mean 0.351 in females (0.286–0.369, mean 0.338 in males, 0.283–0.354, mean 0.312 in females).

Description of the holotype

(Figs 8 and 9). Adult male, as indicated by dark colored throat region and presence of vocal slits; SVL 23.56 mm; habitus robust; head broad, about as wide as long, ratio HL/HW 1.06; snout nearly rounded in dorsal view, projecting beyond lower jaw, rounded in profile; nostril dorsolateral, closer to tip of snout than eye; canthus rounded; ratio EYD/SVL 0.12; IOD (2.12) greater than width of upper eyelid (1.94); tympanum concealed, slightly depressed relative to skin of temporal region, tympanic rim weakly elevated relative to tympanum; ratio TYD/EYD 0.42; vomerine teeth absent; tongue fleshy, rounded, without notch; tips of all four fingers rounded, not expanded into discs; relative finger lengths III>I>IV>II; no webbing; distinct subarticular tubercles, palmar tubercle distinct, bifid; thenar tubercle large, about twice the size of palmar tubercle; tips of toes rounded, slightly expanded into discs; relative toe lengths IV>III>V>II>I; feet almost fully webbed, webbing formula I 0.8–0.8 II 0.8–0.8 III 0.8–0.9 IV 0.9–0.8 V; a well-developed flap of skin on postaxial side of Toe V from level of outer metatarsal tubercle to distal subarticular tubercle; strong fold on distal one-half of tarsus; large, flap-like inner metatarsal tubercle; outer metatarsal tubercle not differentiated, but rather two tiny tubercles present in that area; skin on top of head and on dorsum and flank smooth; skin on throat and venter shagreen; skin on upper surface of forelimbs and thigh with low, rounded tubercles, that of dorsal surface of shank with scattered keratinized pointed tubercles; indistinct, glandular supratympanic fold from posterior edge of upper eyelid along upper margin of tympanum and then obliquely down to shoulder; no dorsolateral fold. Measurements (mm) of holotype: SVL 23.56; HL 8.91; HW 8.44; SL 3.60; EYD 2.81; IOD 2.12; TYD 1.17; TED 0.67; SHL 12.51; THL 12.52; HNL 6.04; FL 12.22; NED 2.37; IND 1.92.

Figure 8. 

Holotype of Phrynoglossus myanhessei n. sp. (SMF 103841). Scale bars equal 10.0 mm (A, B) and 5 mm (C–F). Photos by G.K.

Figure 9. 

Holotype of Phrynoglossus myanhessei sp. n. (SMF 103841). Scale bars equal 5.0 mm. Photos by G.K.

Coloration in life was recorded as follows (Fig. 10): Dorsal and lateral ground color of head and body Drab-Gray (256) with indistinct Raw Umber (280) blotches and mottling; dorsal surfaces of forelimbs True Cinnamon (260) with Raw Umber (280) speckles; dorsal surfaces of hind limbs Drab-Gray (256) with Raw Umber (280) transverse broken bars; a Drab-Gray (256) oblique bar from anterior corner of eye to snout; throat region heavily suffused with Vandyke Brown (282); venter Pale Neutral Gray (296) with Smoky White (261) stipples; ventral surfaces of forelimbs Medium Fawn Color (257); ventral surfaces of hind limbs Cream White (52); iris Olive Brown (278) with a suffusion of Orange-Rufous (56) above and whitish below.

Figure 10. 

Phrynoglossus myanhessei sp. n. in life. (A) Male holotype, SMF 103841; (B) male paratype, SMF 103840; (C) male paratype, SMF 103799; (D) female paratype, SMF 103798. Photos by G.K.

Coloration after about three years preservation in 70% ethanol was recorded as follows: Dorsal and lateral ground color of head and body Glaucous (272) with indistinct Sepia (279) blotches and mottling; dorsal surfaces of hind limbs Glaucous (272) with Sepia (279) transverse broken bars; a Drab-Gray (256) oblique bar from anterior corner of eye to snout; supraocular region Dusky Brown (285); throat region Glaucous (289) with Pale Buff (1) stipples; venter Pale Buff (1) with Smoky White (261) stipples; ventral surfaces of forelimbs Glaucous (289); ventral surfaces of hind limbs Light Buff (2).

Variation

The paratypes agree well with the holotype in general appearance; morphometrics and coloration (see Table 3).

Etymology

“Myan” is Myanmar’s abbreviated name and was chosen because this species is endemic to Myanmar as far as we know. “hessei” was chosen in recognition of the long-term support and funding of Senckenberg by the German State of Hesse. In combination, the species name myanhessei reflects the long-term productive collaboration of researchers from Hesse and Myanmar in the field of herpetology.

Natural history notes

At the type locality, the specimens were collected at night in a patch of muddy grass area, partly open, partly covered by bushes and low trees. The frogs were sitting at the edge of small shallow temporary water bodies.

Geographic Distribution and Conservation

As currently known, Phrynoglossus myanhessei is restricted to central and lower Myanmar (Fig. 11). This species was abundant wherever we found it. Thus, we classify it as Least Concern according to the IUCN categories (IUCN, 2012).

Figure 11. 

Map indicating collecting localities of the Phrynoglossus species occurring in Indochina. Each symbol can represent one or more adjacent localities. Black squares: P. myanhessei n. sp.; black circles: P. martensii; green pentagons: P. cf. martensii; black triangles: P. magnapustulosus; white circles: additional specimens of Phrynoglossus from the FMNH, not examined by authors.

Genomic characterization

Whole genome sequencing and assembly: Illumina sequencing yielded 616,151,068 short-reads with a data amount of 92.4 Gb. K-mer analysis estimated the genome size to 2.6 Gb and heterozygosity 0.9%. The mitochondrial genome was assembled into one circular sequence with a length of 18,348 bp (accession number MW405414). All expected 13 protein coding genes, 2 rRNAs, 23 tRNAs, one D-loop region and the origin of replication could be annotated on the mitochondrial genome sequence. A pairwise alignment to the complete mitochondrial genome of Phrynoglossus martensii (GU177877) shows 77.1% pairwise identity. The final nuclear genome assembly contains 1,446,664 contigs with a total length of 1,829,122,027, an N50 of 1,468 and a GC of 41.07%. The backmapping rate is as high as 98.9% and after filtering, the genome sequence coverage of the holotype specimen is uniformly distributed at 25×. The BUSCO search resulted in 20.8% present BUSCOs (C:8.5%[S:8.4%,D:0.1%],F:12.3%,M:79.2%,n:5310) and no contamination could be identified by interpreting the blobplot. Raw reads and the draft-genome assemblies can be found within the BioProject PRJNA687006. For further details see Appendix 3.

Discussion

Our study provides support for the general assumption that most geographically wide-spread species contain cryptic diversity with one or more unrecognized species “hiding” behind another taxon name due to its similarity in external morphology (Hasan et al. 2012; Köhler et al. 2019). Typically, these anuran species differ from one another in bioacoustics, i.e., in the male advertisement call. Often these species are poorly differentiated in external morphology which was the reason for lumping them before in the first place. The reason for the poor degree of morphological differentiation is possibly the lack of ecological differentiation of the members of such a cluster of frogs. In the case of the Phrynoglossus studied here, it is most likely that they went through the process of speciation in allopatry as indicated by their present distribution pattern. A single prezygotal isolating mechanism is enough to finalize the speciation process no matter how similar the diverging populations are in other characters such as external morphology or ecological traits, e.g., habitat preference, diet, or activity patterns (Wilson, 2001). In anurans the male advertisement call serves as a very effective isolating mechanism to avoid hybridization among similar species under natural conditions (Köhler et al. 2017). Thus, the discovery that there are three distinct male advertisement calls among the species referred to as P. martensii in mainland Southeast Asia is strong evidence for the presence of three distinct species that are not compatible reproductively. In this sense it could be expected that the variation in mtDNA data would be congruent with the results of the bioacoustical analyses. The mtDNA data found a fourth clade of frogs (our Clade 2 with specimens from northern Thailand, southern China, Vietnam, Laos, and Cambodia) that differed from the other three clades by a mean genetic distance of 16S of 3.5–5.1% and therefore qualifies as candidate species (sensu Fouquet et al. 2007). However, since we did not find additional support for the recognition of this clade as a distinct evolutionary species, we tentatively assigned it to the species it is genetically most closely related too (i.e., P. martensii). Additional research including nuclear genetic data is needed to evaluate these northeastern populations.

In previous publications, the geographic distribution of Phrynoglossus magnapustulosus was given as “Nakhon Phanon, Ubon, Loei, and Chiang Mai provinces, Thailand“ (e.g., Taylor 1962; Frost 2020) or “northern and northeastern Thailand” (Khonsue and Thirakhupt 2001). According to our molecular genetic data and also based upon our field observations in Chiang Mai, this statement is erroneous and probably based on misleading characters of external morphology that were considered to be diagnostic for this species such as the presence of “pearly tubercles” on the dorsum (Taylor 1962). In our series of topotypic specimens of P. magnapustulosus we found that dorsal skin texture varies considerably. Some specimens virtually have only a few low tubercles on the dorsum whereas other have distinct ones, pearl-tipped or not (see Fig. 12). A similar degree of variation in dorsal skin texture was observed in individuals of P. martensii (Fig. 13). Thus, we dismiss this character as useful for distinguishing among species of Phrynoglossus. As far as known, P. magnapustulosus is restricted to the Khorat Plateau. Other herpetofaunal species endemic to the Khorat Plateau include Enhydris subtaeniata and Malayemys khoratensis (Murphy and Voris 2014; Ihlow et al. 2016).

Figure 12. 

Phrynoglossus magnapustulosus in life. (A) GK-7855; (B) GK-7853; (C) GK-7877; (D) GK-7880. Photos by G.K.

Figure 13. 

Phrynoglossus martensii in life. (A) GK-7350; (B) GK-7349; (C) GK-7235; (D) GK-7237. Photos by G.K.

The species of Phrynoglossus (or Occidozyga) that are supposedly endemic to the islands of the Sunda Archipelago and the Philippines were beyond the scope of this study. Future studies need to address the phylogenetic relationships and verify their generic assignment as proposed here. The integrative taxonomic approach including also bioacoustic and genetic data along with a thorough analysis of the geographic variation of external morphology will be useful in order to clarify the taxonomic status of these populations.

The evidence for recognizing Phrynoglossus and Occidozyga as distinct evolutionary units at genus level is supported by the monophyly of the two taxa based on the analyses of mtDNA in this study. Furthermore, there are several diagnostic morphological characters that define each genus, some of which may turn out to be apomorphic for one or the other genus. The conspicuous skin morphology of Occidozyga as well as its unique tongue shape are such candidates. Regarding the documented mode of amplexus, finding two different modes in these supposedly closely related genera is interesting. The inguinal amplexus mode that we documented in P. magnapustulosus and P. myanhessei is supported as an autapomorphy of the genus Phrynoglossus by the observation of an inguinal amplexus mode also in P. sumatranus by Eto and Matsui (2012). According to a large-scale analysis of the amplexus mode across all groups of anuran amphibians, reveal the inguinal amplexus as the basal state to all Anura (Carvajal-Castro et al. 2020). However, the latter authors gave the axillary amplexus as the only mode found within the family Dicroglossidae obviously not being aware of the inguinal amplexus mode in Phrynoglossus. Given the phylogenetic position of Phrynoglossus in the amphibian tree of life (Pyron and Wiens 2011), we interpret the amplexus mode in this genus as a reversal from the axillary mode, and thus as an autapomorphy of this genus.

With this description of a new species we provide the complete, annotated mitochondrial genome, a draft genome and 25× coverage of short-read genome resources of the holotype. This fundamental genomic characterization of the species based on the name-bearing specimen will be an invaluable genomic resource for future taxonomic and evolutionary studies. The quality of the nuclear genome assembly provided for the holotype of Phrynoglossus myanhessei is sufficient to genetically characterize the new species on the basis of the name-bearing specimen. The genome data describe the entire genetic variation, including heterozygosity of the holotype individual and make it possible to extract any genetic locus for future taxonomic or phylogenetic studies specifically of occidozygine frogs and frogs in general.

Following the first example from Köhler et al. (2021) we reiterate our proposition that wherever feasible, descriptions of new species should be accompanied by genomic data and a draft genome assembly from the respective holotype as an important resource to future biodiversity research. The investment in a draft genome of newly described holotypes ensures that the genetic fingerprint of newly described species is captured and preserved long-term based on the name-bearing specimen. Furthermore, Köhler et al. (2021) suggest to balance between quality and investment regarding draft genome quality. Therefore the quality of the nuclear genome assembly from 25× coverage is not very high by purpose and still provides a basis for further studies. In addition the complete mitochondrial genome is a valuable resource e.g., for further phylogenetic studies. In the near future, adding genomic data to traditional phenotypic studies will likely become a standard, as costs and efforts to produce them drop dramatically. And the time to start is now.

Acknowledgments

We would like to express our gratitude to Dr. Nyi Nyi Kyaw, Director General, Forest Department, Ministry of Natural Resources and Environmental Conservation, Myanmar, for issuing our collecting permits, U Hla Maung Thein, Director General, Environmental Conservation Department, Ministry of Natural Resources and Environmental Conservation, Myanmar, for issuing our permission to access the genetic resources as in compliance with the Nagoya Protocol, and to Dr. Naing Zaw Htun, Director, Nature and Wildlife Conservation Division, Ministry of Natural Resources and Environmental Conservation, Myanmar, for issuing the exportation permits. Our special thanks are due to Dr. Kyaw Kyaw Khaung, Rector, East Yangon University, Thanlyin, for supporting us to obtain all necessary permissions from the Ministry of Natural Resources and Environmental Conservation, Myanmar. We thank Linda Mogk, Senckenberg Research Institute, Frankfurt, for doing molecular lab work that resulted in some of the new DNA sequences included in this contribution. We are grateful to the following colleagues who provided audio files with recordings of Phrynoglossus: Sandra Goutte, Taiwan; Thomas Ziegler, Köln. The audio files were made available to us via the Fonoteca Zoológica, Museo Nacional de Ciencias Naturales, Madrid, Spain. Melanie Müller and Jakob Adam helped with several tasks in the preparation of data and images for this study. This research is funded partly by Chulalongkorn University: U_GR_63_66_23_10 and also partly financially supported by the Sci-Super VI fund from Faculty of Science, Chulalongkorn University. We thank Noppadon Kitana, Wichase Khonsue and Chatchawan Chaisuekul, Department of Biology, Faculty of Science, Chulalongkorn University, for assisting PT and GK on field work in Thailand. And we are thankful to Andreas Schmitz and Thomas Ziegler for helpful comments on an earlier version of the manuscript.

References

  • Anderson J (1871) A list of the reptilian accession to the Indian Museum, Calcutta from 1865 to 1870, with a description of some new species. Journal of the Asiatic Society of Bengal 40: 12–39.
  • Bioacoustics Research Program (2011) Raven Pro: Interactive Sound Analysis Software (Version 1.4). The Cornell Lab of Ornithology, Ithaca, NY, available from http://www.birds.cornell.edu/raven
  • Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114–2120.
  • Bouckaert R, Vaughan TG, Barido-Sottani J, Duchene S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, De Maio N, Matschiner M, Mendes FK, Müller NF, Ogilvie HA, du Plessis L, Popinga A, Rambaut A, Rasmussen D, Siveroni I, Suchard MA, Wu C-H, Xie D, Zhang C, Stadtler T, Drummond AJ (2018) BEAST 2.5: An Advanced Software Platform for Bayesian Evolutionary Analysis. https://doi.org/10.1101/474296
  • Boulenger GA (1896) Descriptions of new batrachians in the British Museum. Annals and Magazine of Natural History, Series 6, 17: 401–406.
  • Burnham KP, Anderson DR (2002) Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Springer-Verlag New York Inc, New York.
  • Carvajal-Castro JD, López-Aguirre Y, Ospina-L AM, Santos JC, Rojas B, Vargas-Salinas F (2020) Much more than a clasp: evolutionary patterns of amplexus diversity in anurans. Biological Journal of the Linnean Society 129(3): 652–663. https://doi.org/10.1093/biolinnean/blaa009
  • Chikhi R, Medvedev P (2014) Informed and automated k-mer size selection for genome assembly. Bioinformatics 30(1): 31–37.
  • Corlett RT (2014) The ecology of tropical East Asia (2. ed.). Oxford University Press, Oxford
  • Dierckxsens N, Mardulyn P, Smits G (2017) NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Research 45(4): e18-e18
  • Donath A, Jühling FJ, Al-Arab M, Bernhart SH, Reinhardt F, Stadler PF, Middendof M, Bernt M (2019) Improved annotation of protein-coding genes boundaries in metazoan mitochondrial genomes. Nucleic Acids Research 47: 10543–10552.
  • Dubois A (1981) Liste des genres et sous-genres nominaux de Ranoidea (Amphibiens Anoures) du monde, avec identification de leurs espèces types; consequences nomenclaturales. Monitore Zoologico Italiano. Nuova Serie, Supplemento. Firenze 15: 225–284.
  • Dubois A (1987) Miscellanea taxinomica batrachologica (I). Alytes 5: 7–95.
  • Duellman WE, Trueb L (1994) Biology of Amphibians. 2nd edition. New McGraw-Hill, Baltimore and London.
  • Eto K, Matsui M (2012) Field observation of egglaying behavior of a Puddle Frog Occidozyga sumatrana from Bali, Indonesia (Anura: Dicroglossidae). Current Herpetology 31(2): 121–124.
  • Fei L, Ye C-y, Huang Y-z (1990) [Key to Chinese Amphibians]. Publishing House for Scientific and Technological Literature, Chongqing.
  • Fouquet A, Gilles A, Vences M, Marty C, Blanc M, Gemmell NJ (2007) Underestimation of species richness in neotropical frogs revealed by mtDNA analyses. PloS One 2(10) e1109. https://doi.org/10.1371/journal.pone.0001109
  • Funk WC, Caminer M, Ron S (2011) High levels of cryptic species diversity uncovered in Amazonian frogs. Proceedings. Biological Sciences / the Royal Society 279: 1806–1814. https://doi.org/10.1098/rspb.2011.1653
  • Glez-Peña D, Gómez-Blanco D, Reboiro-Jato M, Fdez-Riverola F, Posada D (2010) Alter: Program-oriented conversion of DNA and protein alignments. Nucleic Acids Research 38 (Web Server issue): 14–18. https://doi.org/10.1093/nar/gkq321
  • Gravenhorst JLC (1829) Deliciae Musei Zoologici Vratislaviensis. Fasciculus primus. Chelonios et Batrachia. Leopold Voss, Leipzig.
  • Günther ACLG (1858) Neue Batrachier in der Sammlung des britischen Museums. Archiv für Naturgeschichte 24: 319–328.
  • Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O (2010) New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Systematic Biology 59(3): 307–321. https://doi.org/10.1093/sysbio/syq010
  • Hasan M, Islam MM, Khan MR, Alam M S, Kurabayashi A, Igawa T, Kuramoto M, Sumida M (2012) Cryptic anuran biodiversity in Bangladesh revealed by mitochondrial 16S rRNA gene sequences. Zoological Science 29(3): 162–172. https://doi.org/10.2108/zsj.29.162
  • Heyer WR, Rand AS, Goncalves da Cruz CA, Peixoto OL, Nelson CE (1990) Frogs of Boraceia. Arquivos de Zoologia 31(4): 235–410.
  • Ihlow F, Vamberger M, Flecks M, Hartmann T, Cota M, Makchai S, Meewattana P, Dawson JE, Kheng L, Rödder D, Fritz U (2016) Integrative taxonomy of Southeast Asian snail-eating turtles (Geoemydidae: Malayemys) reveals a new species and mitochondrial introgression. PloS One 11(4) e0153108. https://doi.org/10.1371/journal.pone.0153108
  • Inger RF (1954) Systematics and zoogeography of Philippine Amphibia. Fieldiana. Zoology 33: 183–531.
  • Inger RF (1996) Commentary on a proposed classification of the family Ranidae. Herpetologica 52: 241–246.
  • Iskandar DT, Arifin U, Rachmanasah A. (2011) A new frog (Anura, Dicroglossidae), related to Occidozyga semipalmata Smith, 1927, from the eastern peninsula of Sulawesi, Indonesia. Raffles Bulletin of Zoology. Singapore 59: 219–228.
  • Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A (2012) Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28: 1647–1649. https://doi.org/10.1093/bioinformatics/bts199
  • Kekkonen M, Mutanen M, Kaila L, Nieminen M, Hebert PDN (2015) Delineating species with DNA Barcodes: A case of taxon dependent method performance in moths. PloS One Article ID e0122481. https://doi.org/10.1371/journal.pone.0122481
  • Khonsue W, Thirakhupt K (2001) A checklist of the amphibians in Thailand. The Natural History Journal of Chulalongkorn University 1: 69–82.
  • Köhler G (2012) Color Catalogue for Field Biologists. Herpeton, Offenbach.
  • Köhler G, Khaing KPP, Than NL, Baranski D, Schell T, Greve C, Janke A, Pauls SU (2021) A new genus and species of mud snake from Myanmar (Reptilia, Squamata, Homalopsidae). Zootaxa 4915(3): 301–325.
  • Köhler G, Mogk L, Khaing KPP, Than NL (2019) The genera Fejervarya and Minervarya in Myanmar: Description of a new species, new country records, and taxonomic notes (Amphibia, Anura, Dicroglossidae). Vertebrate Zoology 69: 183–226. https://doi.org/10.26049/VZ69-2-2019-05
  • Köhler J, Jansen M, Rodríguez A, Kok PJR, Toledo LF, Emmrich M, Glaw F, Haddad CFB, Rödel M-O, Vences M (2017) The use of bioacoustics in anuran taxonomy: Theory, terminology, methods and recommendations for best practice. Zootaxa 4251(1): 1–124. https://doi.org/10.11646/zootaxa.4251.1.1
  • Kuhl H van, Hasselt JC (1822) Uittreksels uit breieven van de Heeren Kuhl en van Hasselt, aan de Heeren C. J. Temminck, Th. van Swinderen en W. de Haan. Algemeene Konst-En Letter-Bode 7: 99–104.
  • Kumar S, Stecher G, Tamura K (2016) Mega7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Molecular Biology and Evolution 33(7): 1870–1874. https://doi.org/10.1093/molbev/msw054
  • Laetsch DR, Blaxter ML (2017) BlobTools: Interrogation of genome assemblies. F1000Research 6: 1287.
  • Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B (2017) Partitionfinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Molecular Biology and Evolution 34(3): 772–773. https://doi.org/10.1093/molbev/msw260
  • Li H. (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv: 1303.3997v2
  • Manthey U, Grossmann W (1997) Amphibien und Reptilien Südostasiens. Natur und Tier-Verlag, Münster.
  • Mertens R (1927) Neue Amphibien und Reptilien aus dem Indo-Australischen Archipel, gesammelt während der Sunda-Expedition Rensch. Senckenbergiana biologica 9: 234–242.
  • Minh BQ, Nguyen MAT, Haeseler A von (2013) Ultrafast approximation for phylogenetic bootstrap. Molecular Biology and Evolution 30(5): 1188–1195. https://doi.org/10.1093/molbev/mst024
  • Mittermeier RA (2004) Hotspots revisited: Earth’s biologically richest and most endangered terrestrial ecoregions (1. engl. ed.). Mexico City: CEMEX.
  • Murphy JC, Voris HK (2014) A checklist and key to the homalopsid snakes (Reptilia, Squamata, Serpentes), with the description of new genera. Fieldiana Zoology 1567 1–43.
  • Myers N, Mittermeier RA, Mittermeier CG, da Fonseca Gustavo AB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature 403: 853–858. https://doi.org/10.1038/35002501
  • Nguyen VS, Ho TC, Nguyen QT (2009) Herpetofauna of Vietnam. Edition Chimaira, Frankfurt.
  • Ogilvie HA, Bouckaert RR, Drummond AJ (2017) Starbeast2 brings faster species tree inference and accurate estimates of substitution rates. Molecular Biology and Evolution 34(8): 2101–2114. https://doi.org/10.1093/molbev/msx126
  • Ohler A, Dubois A (1999) The identity of Elachyglossa gyldenstolpei Andersson 1916 (Amphibia, Ranidae), with comments on some aspects of statistical support to taxonomy. Zoologica Scripta 28(3–4): 269–279. https://doi.org/10.1046/j.1463-6409.1999.00002.x
  • Peters WCH (1867) Herpetologische Notizen. Monatsberichte der Königlichen Preussischen Akademie der Wissenschaften zu Berlin 1867: 13–37.
  • Peters WCH (1877) Herpetologische Notizen. II. Bemerkungen über neue oder weniger bekannte Amphibien. Monatsberichte der Königlichen Preussischen Akademie der Wissenschaften zu Berlin 1877: 415–423
  • Pyron R, Wiens JJ (2011) A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Molecular Phylogenetics and Evolution 61: 543–583. https://doi.org/10.1016/j.ympev.2011.06.012
  • Ronquist F, Teslenko M, Mark P, Ayres DL, Darling A, Höhna SH, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space. Systematic Biology 61: 539–542. https://doi.org/10.1093/sysbio/sys029
  • Sambrook J, Russell WD (2001) Protocol 1: DNA isolation from mammalian tissue. Pp. 623–627 in: Sambrook J, Russell WD (eds.) Molecular Cloning: A Laboratory Manual. 3rd Ed. Cold Spring Harbor Laboratory Press, New York
  • Savage JM, Heyer WR (1997) Digital webbing formulae for anurans: a refinement. Herpetological Review 28(3): 131.
  • Sclater WL (1892) On some specimens of frogs in the Indian Museum, Calcutta with description of several new species. Proceedings of the Zoological Society of London 1892: 341–348.
  • Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM (2015) BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31: 3210–3212.
  • Simpson GG (1951) The species concept. Evolution 5: 285–298.
  • Smith MA (1916) On the frogs of the genus Oxyglossis. Journal of the Natural History Society of Siam 2: 172–175.
  • Smith MA (1927) Contributions to the herpetology of the Indo Australian region. Proceedings of the Zoological Society of London 1927: 199–225.
  • Smith MA (1931) The herpetology of Mt. Kinabalu, North Borneo, 13,455 ft. Bulletin of the Raffles Museum 5: 3–32.
  • Stejneger L (1925) Chinese amphibians and reptiles in the United States National Museum. Proceedings of the United States National Museum 66: 1–115.
  • Taylor EH (1922) Additions to the herpetological fauna of the Philippine Islands. II. Philippine Journal of Science 21: 257–303.
  • Taylor EH (1962) The amphibian fauna of Thailand. University of Kansas Science Bulletin 43: 265–599.
  • Taylor EH, Elbel RE (1958) Contribution to the herpetology of Thailand. University of Kansas Science Bulletin 38: 1033–1189.
  • Tillich M, Lehwark P, Pellizzer T, Ulbricht-Jones ES, Fischer A, Bock R, Greiner S (2017) GeSeq–versatile and accurate annotation of organelle genomes. Nucleic Acids Research 45(W1): W6-W11
  • Trifinopoulos J, Nguyen L-T, Haeseler A von, Minh BQ (2016) W-IQ-TREE: A fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Research, 44(W1): W232-W235. https://doi.org/10.1093/nar/gkw256
  • Tschudi JJ von (1838) Classification der Batrachier mit Berücksichtigung der fossilen Thiere dieser Abtheilung der Reptilien. Neuchâtel: Petitpierre.
  • Vurture GW, Sedlazeck FJ, Nattestad M, Underwood CJ, Fang H, et al. (2017) GenomeScope: fast reference-free genome profiling from short reads. Bioinformatics 33: 2202–2204.
  • Wiley EO (1978) The Evolutionary Species Concept Reconsidered. Systematic Biology 27: 17–26.
  • Wilson EO (2001) The diversity of life (New ed.). Penguin, London, 406 pp.
  • Wood DE, Lu J, Langmead B (2019) Improved metagenomic analysis with Kraken 2. Genome Biology 20: 257.
  • Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Research 18(5): 821–829.
  • Ziegler T (2002) Die Amphibien und Reptilien eines Tieflandfeuchtwald-Schutzgebietes in Vietnam. Natur und Tier-Verlag, Münster.

Appendix 1

Comparative Specimens Examined

Occidozyga limaChina: Guangdong: no further data: SMF 6571–72; Lo-fou-shan Mountains: SMF 6573; Hainan: SMF 6570. Indonesien: Java: no further data: SMF 6574–75; Weltevreden: SMF 23436. Myanmar: Bago: Bago Yoma, 125 m: GK-7075–77; Ko Pya Gyi, 55 m: GK-7072; Naypyidaw: near Yamethin, 206 m: SMF 103817–19; Yangon: East Yangon University, 65 m: GK-7057, 7733; near Taw Hlan village: 1: SMF 103815–16. Thailand: Nakhon Phanom: Ban Kan Luang, 173 m: GK-7881–82; Nakhon Ratchasima: near Lake Resort Khorat, 142 m: GK-7917, 7920; Rayong: near Rayong, 110 m: GK-7721; Roi Et: near Ban Sa At Na Di, 160 m: GK-7409–10, 7415–17, 7451–52, 7477–79, 7538–40; near Selaphum, 130 m: GK-7794; Sakon Nakhon: Ban Phaeng Yai: 155: GK-7526. Thailand: Yasothon: 13 km NE Ban Sa At Na Di: 202: GK-7469.

Phrynoglossus celebensisIndonesia: Sulawesi Selatan: Djikoro, Mt. Bouthain: SMF 16332.

Phrynoglossus floresianusIndonesia: Flores: Nusa Tenggara Barat, Rana Mese: SMF 23438–45.

Phrynoglossus laevisPhilippines: Laguna: Mt. Makiling, Los Banos: SMF 74674; Leyte: no further data: SMF 6591; Mt. Balocaue, Baybay: SMF 74614; Luzon: Manila: SMF 6578–80, 6587; Central Luzon: SMF 6581–86; Masbate: Panal: SMF 74317–30; Negros Oriental: Lake Balinsassyao: SMF 75006; Palawan: Calauit Island, northwestern peak of the Busunga: SMF 74134, 74437; Tarusan: SMF 74421–33; Quirino: Sierra Madre: SMF 74637–39; Samar: Hinabangen: SMF 75202–04.

Phrynoglossus myanhesseiMyanmar: Bago: Bago Yoma, 425 m: GK-7104–05, 7139; Bago, 50 m: SMF 103351–53; Magwe: near Taungdwingyi, 170 m: SMF 103800–03; Rakhaing: Ngapali Mountains, 50 m: SMF 103843; Ngapali, Dam Lake, 10 m: SMF 103842; Yangon: East Yangon University, 17 m: SMF 103840–41, 103797–99.

Phrynoglossus magnapustulosusThailand: Nakhon Phanom: Ban Kan Luang, 173 m: GK-7853, 7855, 7876–80, 7887; Nakhon Ratchasima: near Lake Resort Khorat, 142 m: GK-7915–16; Roi Et: near Ban Sa At Na Di, 145 m: GK-7361, 7465, 7496–503, 7805–06, 7395–96, 7412–14, 7475–76, 7509, 7513–16, 7537; Sakon Nakhon: Ban Phaeng Yai, 155 m: GK-7532–34; Yasothon: 13 km NE Ban Sa At Na Di, 202 m: GK-747071.

Phrynoglossus martensiiThailand: Bangkok: Bangkok, 3 m: PT-2634; Chiang Mai: Huai Hong Khrai Royal Development Study Centre, 430 m: GK-7234–35, 7237–38, 7249, PT-1543–441771, 2076; Chaiyaphume: Ban Na Si Nuan, 205: PT-0167–68; Chonburi: Ya Teng Homestay, 30 m: GK-7349–50; Kanchanaburi: Klon Do, Dan Makham Tia District, 33 m: PT-2644–45; Nakornratchasima: Wang Nham Khiau, 515 m: PT-0033–36, 0038–39; Nan: near Klang Wiang, 257 m: PT-2669–70, 2672–73; Prachuap Khiri Khan: Ko Thalu, 8 m: PT-0535; Ratchaburi: Damnoen Saduak, 5 m: GK-7365; Trat: Ko Kut Resort, 35 m: GK-7713–14, PT-1025; Trat, 5 m: GK-7695–701; Satun: near Masayit Gi Ma Min Din, 102 m: PT-2852–53; Songkhla: Wang Pha, 98 m: PT-2754–59.

Phrynoglossus sumatranusIndonesia: Bali: Batoeriti: SMF 23454–55; Gitgit: SMF 23447–53; Java: Jawa Barat: Bogor: SMF 6576–77, 23437, 31165–68; Jawa Tengah: Wonosobo: SMF 31229; Jawa Timur: Punten: SMF 31230–34; Bujutan at Ardjasa, northwestern coast of P. Kangean, Kangean Islands: SMF 55307–08; Sulawesi: no further data: SMF 6600; Sulawesi Utara: Minahasa: SMF 6598–99; Sumatra: Stabat, Deli: SMF 6588–90.

Appendix 2

Genbank accession numbers for the 16S sequences used in this study

Species Specimen number GenBank number
Ingerana tenasserimensis USNM 586923 MG935841
Ingerana tenasserimensis USNM 587300 MG935839
Ingerana tenasserimensis USNM 587302 MG935840
Ingerana tenasserimensis USNMFS35684 MG935838
Limnonectes limborgi GK_7110 MW217495
Occidozyga lima AB530619
Occidozyga lima GK_7409 MW217509
Occidozyga lima GK_7451 MW217508
Occidozyga lima GK_7452 MW217507
Occidozyga lima GK_7721 MW217498
Occidozyga lima GK_7733 MW217497
Occidozyga lima GK_7794 MW217496
Occidozyga lima GK_7924 MW217506
Occidozyga lima ROM 25003 AF206497
Occidozyga lima SMF 103815 MW217492
Occidozyga lima SMF 103817 MW217494
Occidozyga lima SMF 103818 MW217493
Occidozyga lima KR827958
Occidozyga lima KR827959
Occidozyga lima KR827960
Phrynoglossus myanhessei SMF 103797 MW217501
Phrynoglossus myanhessei SMF 103798 MW217502
Phrynoglossus myanhessei SMF 103800 MW217503
Phrynoglossus myanhessei SMF 103840 MW217499
Phrynoglossus myanhessei SMF 103841 MW217500
Phrynoglossus myanhessei USNM 587105 MG935916
Phrynoglossus myanhessei USNM 587107 MG935920
Phrynoglossus myanhessei USNM 587395 MG935918
Phrynoglossus myanhessei USNM 587402 MG935917
Phrynoglossus magnapustulosus KR827981
Phrynoglossus magnapustulosus GK_7395 MW217488
Phrynoglossus magnapustulosus GK_7396 MW217487
Phrynoglossus magnapustulosus GK_7532 MW217486
Phrynoglossus magnapustulosus GK_7533 MW217485
Phrynoglossus magnapustulosus GK_7855 MW217490
Phrynoglossus magnapustulosus GK_7916 MW217489
Phrynoglossus martensii AF285214
Phrynoglossus martensii GU177877
Phrynoglossus martensii KP318725
Phrynoglossus martensii KR827982
Phrynoglossus martensii KR827983
Phrynoglossus martensii KR827984
Phrynoglossus martensii KR827985
Phrynoglossus martensii AM07357 NC_014685
Phrynoglossus martensii AMNH A161171 DQ283357
Phrynoglossus martensii GK_7349 MW217491
Phrynoglossus martensii GK_7695 MW217504
Phrynoglossus martensii GK_7713 MW217505
Phrynoglossus martensii PT_0167 MW217484
Phrynoglossus martensii PT_0600 MW217481
Phrynoglossus martensii PT_0942 MW217480
Phrynoglossus martensii PT_0943 MW217479
Phrynoglossus martensii PT_1172 MW217478
Phrynoglossus martensii PT_1543 MW217477
Phrynoglossus martensii PT_1544 MW217476
Phrynoglossus martensii PT_2634 MW217475
Phrynoglossus martensii PT_2754 MW217483
Phrynoglossus martensii PT_2755 MW217482
Phrynoglossus martensii ROM 22222 AF206467
Phrynoglossus martensii SCUM0437980 DQ458254
Phrynoglossus martensii SCUM0437983 DQ458255
Phrynoglossus martensii SCUMH020 DQ458256
Phrynoglossus martensii TAD P324 KR827986
Phrynoglossus martensii USNM 586940 MG935942
Phrynoglossus martensii USNM 586941 MG935929
Phrynoglossus martensii USNM 586942 MG935941
Phrynoglossus martensii USNM 586943 MG935940

Appendix 3

Genomics

Material and Methods. Default parameters are applied, if not stated otherwise.

Raw data and preprocessing. A k-mer profile was created from the raw data using Jellyfish 2.3.0 (Marçais & Kingsford, 2011) including the parameters “-F 2 -C -m 21 -s 1000000000 -t 96” for count and “-t 96” for histo. The resulting histogram was uploaded to the GenomeScope webserver (Vurture et al., 2017) to retrieve certain statistics of the genome.

Low quality bases and adapter sequences were trimmed from the raw reads using Trimmomatic’s 0.39 (Bolger et al., 2014) paired end mode along with the options to create a summary and “-threads 96”. For adapter trimming all adapter sequences provided within Trimmomatic were used. The following trimmers were set: “ILLUMINACLIP:<adapter_all.fa>:2:30:10:8:true SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33”.

To filter out reads originating from possible contamination, Kraken 2.0.9 (Wood et al.,2019) was ran with a standard database built on March 18th, 2020 (including “complete genomes in RefSeq for the bacterial, archaeal, and viral domains, along with the human genome and a collection of known vectors (UniVec_Core)”) and the paired and unpaired trimmed reads as input. For both runs the options “--threads 96 --unclassified-out <out-file(s)>” were set and for the paired run the option “--paired”.

Nuclear genome assembly and quality control. To estimate the best length of k for genome assembly, KmerGenie 1.7051 (Chikhi & Medvedev 2014) was applied, using the raw data as input and the options “--diploid -s 11 -k 141 -t 60”.

The unclassified paired and unpaired reads were assembled using Velvet 1.2.10 (Zerbino & Birney 2008) by first running velveth for a k-mer length of 61 and second running velvetg with the options “-cov_cutoff auto -ins_length 350 -min_contig_lgth 500”.

The quality of the resulting scaffolds was tested by a) mapping the reads used for assembly back to the assembly, b) checking for possible contamination and c) searching for expected orthologous genes. Backmapping was performed with backmap.pl 0.3 (https://github.com/schellt/backmap), which utilized bwa mem 0.7.17-r1188 (Li, 2013), samtools 1.10 (Li et al., 2009), Qualimap 2.2.1 (Okonechnikov et al., 2016), bedtools 2.28.0 (Quinlan & Hall, 2010) and R 4.0.2 (R Core Team, 2020). Contamination screening on the assembly level was performed with blobtools 1.1.1 (Laetsch & Blaxter, 2017). The bam file resulting from the above mentioned backmapping was coverted to a blobtools readable cov file by first indexing the bam file with samtools index and second converting with blobtools map2cov. To assign Taxonomy IDs blastn 2.10.0+ (Camacho et al., 2009) was used to align the scaffolds against the complete nt database (-task megablast -outfmt ‘6 qseqid staxids bitscore’ -evalue 1e-25 -num_threads 96). From the cov and hits file a blobDB was created and plotted. Completeness in terms of core orthologs was screened with BUSCO 4.1.4 (Simão et al., 2015) along with the tetrapoda_odb10 set and the options “-c 8 -o GK_6728_velvet -m geno --long --offline

Mitochondrial genome assembly and annotation. NOVOplasty 4.2 (Dierckxsens et al., 2017) was used along with the longest annotated gene (ND5; ACZ02636.1) from Occidozyga martensii (GU177877.1) as seed. The config file was changed compared to the included one for the options “Max memory = 100; Read Length = 150; Insert size = 300” and the options “Reference sequence, Chloroplast sequence” were left blank. As input reads the untrimmed raw data was used as recommended. The mt genome assembly was checked for consistency by alignment against the closest available mt genome of Occidozyga martensii (GU177877.1). This was computed with clustalo 1.2.3 (Sievers et al., 2011) in Geneious Prime 2020.2.3 (https://www.geneious.com) using Java 11.0.6+10 (clustalo-1.2.3-Ubuntu-x86_64 -i input.fasta -o clustal.aln -v --outfmt=clustal --output-order=tree-order --iter=0 --cluster-size=100 -t DNA). The assembled mt genome sequence was manually cut based on the alignment to fit to the cut site of the reference. The cut mt genome sequence was submitted to two mt genome annotation webservers: GeSeq (Tillich et al., 2017) and MITOS2 (Donath et al., 2019). For GeSeq the options circular and mitochondrial were chosen as well as tRNAscan-SE 2.0.6 (Chan & Lowe 2019; Chan et al., 2019) was enabled with sequence source as “vertebrate mitochondrial tRNAs”. Furthermore the above mentioned mt reference Sequence (GU177877.1) was used as BLAT reference sequence. The MITOS2 settings can be found in Table A1. The annotations of GeSeq and MITOS2 were manually merged and curated in Geneious.

Table A1.

Settings used for MITOS2.

Property Value
Reference RefSeq 63 Metazoa
Genetic Code 2
Proteins TRUE
tRNAs TRUE
rRNAs TRUE
OH TRUE
OL TRUE
Circular TRUE
Use Al Arab et al. FALSE
E-value Exponent 2.0
Final Maximum Overlap 50nt
Fragment Quality Factor 100.0
Standard Code FALSE
Cutoff 50.0%
Clipping Factor 10.0
Fragment Overlap 20.0%
Local only TRUE
Sensitive only FALSE
ncRNA overlap: 50 nt

Assembly finalization. The final mt genome assembly from NOVOplasty was blasted against the assembly from Velvet to remove contigs representing the mt genome with blastn 2.10.0 and the options “-num_threads 32 -outfmt ‘6 std slen’”. Blast reported multiple hits for 3 contigs. Either the single hit of one contig or the longest hit of one contig align completely to the NOVOplasty assembly. The only exception is one contig aligning with 43% of its 609bp length. Finally, all contigs producing blast hits were considered as mitochondiral origin and removed.

Results

Raw data and preprocessing. Illumina NovaSeq 6000 sequencing yielded 616,151,068 reads with a data amount of 92.4Gb. The GenomeScope results are accessible via the permalink http://genomescope.org/analysis.php?code=uOlM17rRo9kkBpguBh0f. The genome size was estimated as 2.6Gb and heterozygosity to 0.9% (See Figure A1).

Figure A1. 

K-mer profile anaylsis and estimates of genome metrics of GenomeScope.

Raw read trimming was survived by 93.23% of the read pairs, 3.12% forward only, 2.17% reverse only and 1.49% reads were dropped. Of the trimmed paired reads 3.43% and 2.89% of the trimmed unpaired reads were classified via Kraken2.

Nuclear genome assembly and quality control. The Velvet assembly resulted in 1,446,672 contigs with a total length of 1,829,139,474 and an N50 of 1,468. The GC is at 41.07% and uniformly distributed. Of all assembled reads 98.9% could be mapped back to the assembly with a uniform distribution at 25× (Figure A2). No clear cluster representing contamination could be identified in the blobplot (Figure A3). Since the assembly is quite fragmented, there are no clear distinguishable clusters. Blast hits other than Chordata might arise from false positive hits and do not show enough evidence for contamination – especially because non-Chordata hits are sparse and coverage and GC of non-Chordata hits is similar to the rest. The BUSCO search resulted in C:8.5%[S:8.4%,D:0.1%],F:12.3%,M:79.2%,n:5310.

Figure A2. 

Coverage distribution of the assembly. A genome size estimation based on mapped nucleotides and coverage (Schell et al., 2017) with backmap.pl results in 2.83Gb, which is in line with the k-mer based result of 2.56Gb.

Figure A3. 

Blobplot of the assembly.

Mitochondrial genome assembly and annotation. The NOVOplasty assembly resulted in 2 contigs which could be circularized in one way. After the clustalo alignment was performed and manually adjusted by cutting the NOVOplasty assembly at the same site as the reference mt genome from Occidozyga martensii (GU177877.1). Since this alignment showed no huge structural differences, the assembly was used for annotation. The GeSeq annotation could not find the two rRNAs, the D-loop region and the origin of replication. The 13 protein coding genes were automatically found but with wrong reading frame and partially unlikely start and end points. In total 40 tRNAs were annotated because GeSeq does not merge the annotation from different evidences (here BLAT of reference tRNAs and tRNAscan). MITOS2 automatically annotated the two rRNAs, 13 protein coding genes, 23 tRNAs, the D-loop region in two fragments and the origin of replication. Manual curation mainly comprises taking and adjusting the rRNA, origin of replication and D-loop region annotations from MITOS2. The both fragments of the D-loop region were merged into one. Annotations of protein coding genes were adjusted from GeSeq and tRNA annotations were taken from tRNAscan.

Assembly finalization. Of all contigs 8 were identified as mitochondrial origin, of which 1 has smaller identity than 90% (minimum 87.5%) to the NOVOplasty mt genome assembly. The removed 8 contigs cover 94.9% of all positions of the NOVOplasty mt genome assembly. After removing contigs of mitochondrial origin, the final assembly contains 1,446,664 contigs with a total length of 1,829,122,027bp, an N50 of 1,468 and an GC of 41.07%.

References

Marçais G, Kingsford C (2011). A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics, 27(6):764–770.

Vurture GW, Sedlazeck FJ, Nattestad M, Underwood CJ, Fang H et al. (2017). GenomeScope: fast reference-free genome profiling from short reads. Bioinformatics, 33(14):2202–2204.

Bolger AM, Lohse M, Usadel B (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 30(15):2114–2120.

Wood, D. E., Lu, J., & Langmead, B. (2019). Improved metagenomic analysis with Kraken 2. Genome biology, 20(1), 257.

Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R (2009). The sequence alignment/map format and SAMtools.Bioinformatics,25(16), 2078–2079.

Okonechnikov K, Conesa A, García-Alcalde F (2016). Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data.Bioinformatics,32(2), 292–294.

Quinlan AR, Hall, IM (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics,26(6), 841–842.

R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org

Meng, G., Li, Y., Yang, C., & Liu, S. (2018). MitoZ: A toolkit for mitochondrial genome assembly, annotation and visualization. bioRxiv, 489955.

Dierckxsens, N., Mardulyn, P., & Smits, G. (2017). NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic acids research, 45(4), e18–e18.

Sievers, F., Wilm, A., Dineen, D., Gibson, T. J., Karplus, K., Li, W., ... & Thompson, J. D. (2011). Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Molecular systems biology, 7(1), 539.

Tillich, M., Lehwark, P., Pellizzer, T., Ulbricht-Jones, E. S., Fischer, A., Bock, R., & Greiner, S. (2017). GeSeq–versatile and accurate annotation of organelle genomes. Nucleic acids research, 45(W1), W6–W11.

Donath, A., Jühling, F., Al-Arab, M., Bernhart, S. H., Reinhardt, F., Stadler, P. F., ... & Bernt, M. (2019). Improved annotation of protein-coding genes boundaries in metazoan mitochondrial genomes. Nucleic acids research, 47(20), 10543–10552.

Chan, P. P., & Lowe, T. M. (2019). tRNAscan-SE: searching for tRNA genes in genomic sequences. In Gene Prediction (pp. 1–14). Humana, New York, NY.

Chan, P. P., Lin, B. Y., Mak, A. J., & Lowe, T. M. (2019). tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes. bioRxiv, 614032.

Laetsch DR & Blaxter ML (2017). BlobTools: Interrogation of genome assemblies.F1000Research,6(1287), 1287.

Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM (2015). BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics,31(19), 3210–3212.

Chikhi, R., & Medvedev, P. (2014). Informed and automated k-mer size selection for genome assembly. Bioinformatics, 30(1), 31–37.

Zerbino, D. R., & Birney, E. (2008). Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome research, 18(5), 821–829.