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)

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.


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. 2019Frost 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 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 1954Inger , 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).
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.
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://www.abi.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 ). Raw reads were trimmed for low quality regions and adapter sequences and filtered for possible contamination using Trimmomatic 0.39 ) and Kraken 2.0.9  respectively. Best k-mer length was estimated by KmerGenie 1.7051 . 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 . 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 ) and MITOS2  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 , searching for possible contamination with Blobtools 1.1.1  and screening for single-copy orthologs with BUSCO 4.1.4 . 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.
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).
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).
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.
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.
Thus, we formally resurrect Phrynoglossus from the synonymy of Occidozyga and define the two genera as follows (data also from Ohler and Dubois 1999):

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. Phrynoglossus baluensis
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."   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)  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). (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 ( Variation. The paratypes agree well with the holotype in general appearance; morphometrics and coloration (see Table 3).

Description of the holotype
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). 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×.

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.    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 be tween 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 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  to retrieve certain statistics of the genome. Low quality bases and adapter sequences were trimmed from the raw reads using Trimmomatic's 0.39  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: "ILLUMINA-CLIP:<adapter_all.fa>:2:30:10:8:true SLIDINGWIN-DOW:4:20 MINLEN:50 TOPHRED33".
To filter out reads originating from possible contamination, Kraken 2.0.9  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  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 , 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 . 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  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  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-or-der=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  and MITOS2 . For GeSeq the options circular and mitochondrial were chosen as well as tRNAscan-SE 2.0.6  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 MI-TOS2 were manually merged and curated in Geneious.
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).
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 frag-    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%.