Research Article |
Corresponding author: Panupong Thammachoti ( panupong.th@chula.ac.th ) Academic editor: Raffael Ernst
© 2021 Gunther Köhler, Joseph Vargas, Ni Lar Than, Tilman Schell, Axel Janke, Steffen U. Pauls, Panupong Thammachoti.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Köhler G, Vargas J, Than NL, Schell T, Janke A, Pauls SU, Thammachoti P (2021) 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). Vertebrate Zoology 71: 1-26. https://doi.org/10.3897/vz.71.e60312
|
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.
amplexus mode, bioacoustics, cryptic species diversity, Dicroglossidae, genome, Khorat Plateau, Myanmar, new species, Occidozyga, Phrynoglossus myanhessei sp. n., Thailand
Southeast Asia is recognized as a major global biodiversity hotspot (
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.,
The genus name Occidozyga was introduced by Kuhl and van Hasselt 1822. Its type species (by subsequent designation of
In 1896, Boulenger introduced his new species Oreobatrachus baluensis from „Mount Kina Balu, North Borneo [= Sabah]“, Malaysia (Borneo).
In his treatment of the herpetofauna of Mount Kinabalu, Borneo,
Further evidence for considering Phrynoglossus as a synonym of Occidozyga was provided in the large-scale molecular phylogeny of amphibians of
Specimens examined for this study were personally collected by GK, NLT, and PT (see Appendix
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
In evaluating species’ boundaries within and among populations, we followed the evolutionary species concept (
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 (
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. (
We extracted DNA following the protocol of
We aligned the sequences with MUSCLE (
Bayesian Inference analyses (BI) used MrBayes 3.2 (
De novo whole genome sequencing and assembly began with genomic DNA extraction from muscle tissue of the holotype according to standard phenol/chloroform procedures (
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 (
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
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.
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
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) |
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) |
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.
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.
Thus, we formally resurrect Phrynoglossus from the synonymy of Occidozyga and define the two genera as follows (data also from
Rana lima Gravenhorst, 1829.
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.
Occidozyga lima (Gravenhorst, 1829)
Phrynoglossus martensii Peters, 1867
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.
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.”
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).
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
Coloration in life was recorded as follows (Fig.
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).
The paratypes agree well with the holotype in general appearance; morphometrics and coloration (see Table
“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.
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.
As currently known, Phrynoglossus myanhessei is restricted to central and lower Myanmar (Fig.
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
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
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 (
In previous publications, the geographic distribution of Phrynoglossus magnapustulosus was given as “Nakhon Phanon, Ubon, Loei, and Chiang Mai provinces, Thailand“ (e.g.,
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
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
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.
Comparative Specimens Examined
Occidozyga lima—China: Guangdong: no further data:
Phrynoglossus celebensis—Indonesia: Sulawesi Selatan: Djikoro, Mt. Bouthain:
Phrynoglossus floresianus—Indonesia: Flores: Nusa Tenggara Barat, Rana Mese:
Phrynoglossus laevis—Philippines: Laguna: Mt. Makiling, Los Banos:
Phrynoglossus myanhessei—Myanmar: Bago: Bago Yoma, 425 m: GK-7104–05, 7139; Bago, 50 m:
Phrynoglossus magnapustulosus—Thailand: 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 martensii—Thailand: 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 sumatranus—Indonesia: Bali: Batoeriti:
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 |
|
MW217492 |
Occidozyga lima |
|
MW217494 |
Occidozyga lima |
|
MW217493 |
Occidozyga lima | KR827958 | |
Occidozyga lima | KR827959 | |
Occidozyga lima | KR827960 | |
Phrynoglossus myanhessei |
|
MW217501 |
Phrynoglossus myanhessei |
|
MW217502 |
Phrynoglossus myanhessei |
|
MW217503 |
Phrynoglossus myanhessei |
|
MW217499 |
Phrynoglossus myanhessei |
|
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 |
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 (
Low quality bases and adapter sequences were trimmed from the raw reads using Trimmomatic’s 0.39 (
To filter out reads originating from possible contamination, Kraken 2.0.9 (
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 (
Mitochondrial genome assembly and annotation. NOVOplasty 4.2 (
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
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
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.