Genetic diversity of the Pegu Rice Frog, Microhyla berdmorei (Anura: Microhylidae) based on mitochondrial DNA

Abstract

The Pegu Rice Frog, Microhyla berdmorei is distributed across ten Asian countries. However, the DNA barcoding information (COI gene) is restricted to only Southeast Asian countries. Here, we sampled a specimen of M. berdmorei in Mizoram state, northeast India to allow the genetic diversity of the species across its range. We generated both COI and 16S ribosomal RNA sequences of the studied species to check the population genetic diversity. The Bayesian analyses clearly discriminate M. berdmorei from its sister species Microhyla pulchra. The present datasets of M. berdmorei also revealed 11 and 19 haplotypes with high uncorrected pairwise genetic distances in COI (3.8–11.8%) and 16S rRNA (0–4.6%) gene, respectively. Owing to the high intra-species genetic distances and different haplotypes with sufficient mutational steps in both mitochondrial genes, this study affirms the existence of M. berdmorei species complex or cryptic diversity within its range distribution in South and Southeast Asia.

Keywords:

• Biodiversity hotspot
• amphibians
• phylogeny
• haplotype
• zoogeography

Introduction

The genus Microhyla (family: Microhylidae) contains 50 species around the world; among them, 14 species are reported from India (Frost 2020). This microhylid genus is small-sized terrestrial frogs distributed in the Oriental biogeographic region (Poyarkov et al. 2014, 2018a, 2018b, 2019). Northeast India shares two biodiversity hotspots and is regarded as one of the prolific regions (with 153 known species) for global amphibian diversity (Biju et al. 2019). A total of five Microhyla species, M. berdmorei, M. butleri, M. eos, M. mukhlesuri, and M. mymensinghensis are reported from this region. Among them, M. berdmorei was described early and is distributed in Bangladesh, Cambodia, China, India, Indonesia, Laos, Malaysia, Myanmar, Thailand, and Vietnam (Hasan et al. 2014; Biju et al. 2019). This species generally inhabits moist evergreen forest adjacent to hilly regions and is often found near streams (van Dijk et al. 2004). The species is assessed to be ‘Least Concern’ by the International Union for Conservation of Nature Red List of Threatened Species (IUCN 2021); however, immense deforestation and other anthropogenic threats may decline the local population of M. berdmorei throughout its range distribution. Although, the recent discovery of M. eos from northeast India has increased the checklist of Microhyla, this genus contributed to present persistent taxonomic challenges (Gorin et al. 2020). Nevertheless, high homoplasy and cryptic diversity often obscured accurate diversity estimate of microhylids (Rakotoarison et al. 2017).

The integrated approach combining morphology, molecular, and bioacoustics is evidenced to be the best exercise to illuminate the anuran systematics (Frost et al. 2006; Van Bocxlaer et al. 2006; Köhler et al. 2017). The species diversity, patterns of distribution, phylogenetic relationships, biogeographic origin, and possible routes of colonization of microhylids, including the genus Microhyla have been repeatedly discussed in several studies (Matsui et al. 2011; De Sá et al. 2012; Peloso et al. 2016; Seshadri et al. 2016; Feng et al. 2017; Tu et al. 2018; Vineeth et al. 2018; Garg et al. 2019; Garg and Biju 2019; Poyarkov et al. 2019; Gorin et al. 2020). However, most of the studies deal with mitochondrial ribosomal RNA genes (12S rRNA and 16S rRNA) to define their existing biological questions (Vences et al. 2005; Vieites et al. 2009). As a result, the well-recognized and widely used DNA barcoding (mitochondrial Cytochrome Oxidase subunit I gene) information (Hebert et al. 2003a, 2003b; Biju et al. 2014; Chambers and Hebert 2016) of Microhyla is scanty in the global database, especially from northeast India. As of now, a few mtCOI sequences of M. berdmorei have been generated from Southeast Asian countries to improve the biodiversity assessment from different life stages and fill the Barcode Index Numbers (BINs) of the Barcode of Life Data System (BOLD) (Grosjean et al. 2015; Mulcahy et al. 2018). Nevertheless, most microhylid studies were based on stable species systematics, allowing assessment of their population genetic diversity. Considering the inadequacy of DNA sequence information of M. berdmorei, the present study intended to generate both mitochondrial COI and 16S rRNA sequences from northeast India, to aid our understanding of genetic diversity in this species. We used both phylogeny and haplotype analyses to estimate the population structure of this microhylid frog from its known range distribution. This primary genetic information will assist rapid and reliable estimates of population genetic diversity.

Materials and methods

An adult specimen of a microhylid frog was collected from Reiek Community Reserve Forest (23.75 N 92.67 E, 1190 m ASL), Mamit district, Mizoram state in northeast India (Figure 1(A)). The exploration was accomplished after acquiring the prior permission (No.A.33011/2/99-CWLW/225) from the Chief Wildlife Warden, Govt. of Mizoram, India. The collected specimen was identified as M. berdmorei based on the below morphological characters (Figure 1(A)). A small frog with smooth skin; dorsally varying from light to dark olive, sometimes with dark brownish spots or marbling scattered on dorsal surface of body; characteristic dark sheds present dorsally in between the eyes, running up to trunk region; few minute light brown spots present on lateral parts of hind limbs as well as on flanks; both throat and chest mottled with dark brown; limbs with faint cross bands; snout-vent length (SVL) of 28.6 mm to 35.7 mm. It can be diagnosed by having characteristic features like, head much broader than length, slightly depressed; snout pointed, a little longer than eyes, projecting a little beyond lower jaw; canthus rostralis prominent; nostrils a little closer to tip of snout than to eyes; tympanum hidden; vomerine teeth absent; forelimbs moderately long; fingers slender, free with rounded tips; length of fingers 1 < 2< 4 < 3, the third finger much longer than snout; subarticular tubercles large and prominent; hindlimbs very long; tibiotarsal articulation reaching beyond tip of snout; heels strongly overlapping when hind limbs folded at right angle to body; tibia 3½ to 4 times as long as broad, more than ⅔ the length of snout to vent; tips of toes swollen into rounded tips which are slightly larger than fingertips; digits with a median groove on the upper surface of the tip separating the upper surface from the lower ones; toes fully webbed; inner metatarsal tubercle prominent and oval; a small rounded outer metatarsal tubercle present (Lalremsanga 2011). The specimen was vouchered (MZMU 1613) in the museum of the Department of Zoology, Mizoram University, India. The muscle tissue was aseptically collected from the hind limb and preserved in 70% molecular grade ethanol.

Figure 1. (A) Map showing the collection localities of Microhyla berdmorei and live individual of M. berdmorei collected from northeast India. (B, C) Bayesian phylogenies based on COI and 16S rRNA genes inferred the clustering pattern of M. berdmorei sampled from different geographical locations in South and Southeast Asian countries. (D, E) The TCS haplotype networks based on COI and 16S rRNA genes showed haplotype network of M. berdmorei collected from South and Southeast Asian countries. The estimated haplotypes are shown in different colored frog icons as represented in the collection locality map.

Total genomic DNA was extracted by using the QIAamp DNA Mini Kit (Qiagen, Valencia, CA) with standard manufacturer’s protocols. The genomic DNA was stored at −30 °C in the Center for DNA Taxonomy laboratory, Zoological Survey of India, Kolkata. The published primer pairs were used to amplify the partial fragment of mitochondrial COI and 16S rRNA genes on a Veriti^® Thermal Cycler (Applied Biosystems, Foster City, CA) with standard thermal profiles (Palumbi 1996; Rassmann 1997; Ward et al. 2005). The 25 µl PCR mixture contains 10 pmol of each primer, 100 ng of DNA template, 1 × PCR buffer, 1.0–1.5 mM of MgCl2, 0.25 mM of each dNTPs, and 1 U of Taq polymerase (Takara BIO Inc., Japan). The amplified PCR products were visualized in 1% agarose gel containing Ethidium bromide (10 mg/ml). The PCR products were further purified by using the QIAquickR Gel extraction Kit (Qiagen, Valencia, CA). The cycle sequencing was performed with the purified PCR products (15 ng), BigDye^®Terminator ver. 3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA), and 3.2 picomoles of the same PCR primer pairs. The thermal profile of the cycle sequencing was set as 96 °C for 1 min, followed by 25 cycles of 96 °C for 10 s, 50 °C for 5 s, and a final extension at 60 °C for 1 min 15 s on the same thermal cycler. The cycle sequencing products were cleaned by BigDye X-Terminator Kit (Applied Bio systems, Foster City, CA), and Sanger sequencing was performed on 48 capillary ABI 3730 Genetic analyzer.

The generated bi-directional sequences were checked through SeqScanner V1.0 (Applied Biosystems Inc., CA) and nucleotide BLAST search platform (https://blast.ncbi.nlm.nih.gov/). The consensus sequences were submitted to the GenBank to acquire the accession numbers. The publicly available 12 COI sequences and 31 16S rRNA sequences of M. berdmorei were acquired from GenBank database. The sequence of closest species, Microhyla pulchra (COI: KR087825 and 16S rRNA: MN534610) were also acquired from GenBank and used as an outgroup in the Bayesian analysis (BA). The final datasets (16S rRNA: 564 bp and COI: 645 bp) were aligned separately through ClustalX (Thompson et al. 1997) and uncorrected pairwise (p) genetic distances were calculated through the MEGAX program with complete deletion of gap/missing data treatment setting (Kumar et al. 2018). The BA phylogenies were constructed in Mr. Bayes 3.1.2 (Ronquist and Huelsenbeck 2003) after estimating the best fit models (16S rRNA: T92 + G with lowest BIC score 2510.578 and COI: HKY + G with lowest BIC score 3530.804) through Mr.Modeltest v2 (Nylander 2004). The Markov chain Monte Carlo (one cold and three hot chains) was run for 1,000,000 generations with 25% burn-in and trees saving at each 100 generation, until the standard deviation of split frequencies reached 0.01 and the potential scale reduction factor for all parameters bordered on 1.0. Both phylogenies were illustrated in the iTOL tool for better representation (https://itol.embl. de/) (Letunic and Bork 2007). The haplotype networks for both genes were built to elucidate the haplotype diversity within the different population of M. berdmorei. The numbers of haplotypes, haplotype diversity (Hd), and the number of polymorphic sites were estimated through DnaSP v6 (Rozas et al. 2017). The haplotype networks was built through PopART (http://popart.otago.ac.nz) (Leigh and Bryant 2015) using the standard Templeton, Crandall and Sing (TCS) method (Clement et al. 2000).

Results and discussion

The plethora of biodiversity in northeast India has gained worldwide attention due to its unique biogeography, representing two global biodiversity hotspots (Eastern Himalaya and Indo-Burma) (Tordoff et al. 2011). Remarkably, the majority of the faunal composition of this region shares their range distribution with neighboring South and Southeast Asian countries due to the transcontinental dispersion (De Bruyn et al. 2014; Das et al. 2019). To date, this region supports 153 anuran species, among them 51 were described since 2001 (Kamei et al. 2012; Biju et al. 2019). Owing to the lack of accessibility exacerbated by difficult terrains and low-lying river basins, communication gap with the variety of ethnic communities and their indigenous languages, the exploration of anurans has been limited. Besides the few opportunistic sampling and subsequent taxonomic research, the large-scale study dealing with the integrated approach is meagerly attempted to explore the anuran diversity in this region.

The collected frog specimen was identified as M. berdmorei based on the morphological characters and substantiated their adaptation in high altitude as described earlier (van Dijk et al. 2004). The integrated approaches were rigorously adopted to delimitate the microhylid frogs and other closely related species throughout South and Southeast Asia and recognizing their phylogeny, origin, diversification, and colonization (Garg and Biju 2019; Gorin et al. 2020). In the recent past, combination with mitochondrial DNA (12S rRNA + 16S rRNA) and nuclear DNA (BDNF gene) has been extensively used to estimate the time-calibrated phylogeny of Microhyla frogs and discuss their taxonomy, biogeography, and evolution (Gorin et al. 2020). However, the mtCOI sequence information of microhylid frogs is still lacking in northeast India.

We generated both COI and 16S rRNA genetic information to elucidate the population genetic diversity of M. berdmorei among the different populations. The study provides the first DNA barcode data of the studied species from northeast India. The generated COI (MT755953) and 16S rRNA (MW578946) sequences depicted 96.43 and 99.64% similarity in nucleotide BLAST search with the published database sequences of M. berdmorei (MG935596 and MN534602) collected from Myanmar and Bangladesh respectively (Mulcahy et al. 2018; Gorin et al. 2020). The COI barcode and 16S rRNA sequences of M. berdmorei generated from distant localities were cohesively clustered in the BA phylogenies (Figure 1(B,C)).

For COI dataset, generated and database sequences revealed a total of 11 haplotypes with 92 mutational steps, and haplotype diversity (Hd) = 0.97 (Figure 1(D)). Two COI sequences generated from Myanmar (MG935594 and MG935597) and two sequences generated from Malaysia (MG386471 and MG386471) were represented by a single haplotype. For 16S rRNA dataset, generated and database sequences depicted a total of 19 haplotypes with 38 mutational steps, and haplotype diversity (Hd) = 0.93 (Figure 1(E)). The 16S rRNA sequence of the studied sample and six sequences (MH807388, MH807387, MH807386, MH807385, MH807384, MN534602) collected from northeast India and Bandarban, Bangladesh formed a single haplotype. Further, four sequences generated from Myanmar (MG935891, MG935890, MG935889, MG935888), two sequences generated from Thailand (MN534604, KR827906), two sequences generated from Thailand and Cambodia (KR827908, KU840572) and three sequences from Malaysia (MN534600, AB530638, AB598338) were represented by a single haplotype. The remaining 16S rRNA sequences were formed distinct haplotypes with sufficient mutational steps.

The studied sample collected from northeast India revealed genetically close with the Burmese collection with 3.8% genetic distance in COI gene. The northeast Indian sample also depicted high genetic distance (7.1–7.8% from Thailand and Laos samples and 11.2–11.8% from Malaysia samples in COI gene. Further, in 16S rRNA gene analysis, the studied sample collected from Mizoram state showed less genetic distances (0–0.9%) compared with the samples collected from different states (Assam and Tripura) in northeast India, Myanmar, and Bangladesh. The studied sample also depicted 1.7–4.6% genetic distance from other samples collected from different localities in Southeast Asian countries (Thailand, Vietnam, Laos, Cambodia, and Malaysia). A recent study reported a number of deep lineages (with >1.5% uncorrected genetic p-distances in 16S rRNA) within the different population of M. berdmorei collected from Southeast Asian countries and indicated potentially constituting cryptic species diversity (Gorin et al. 2020). Hence, considering the high intra-species genetic distances 3.8–11.8% in COI gene and 0–4.6% in 16S rRNA gene, this study corroborates the previous hypothesis and assumed the existence of possible cryptic diversity and species complex of M. berdmorei within its range distribution. However, this study hinted an additional sampling from board geographical regions and their genetic information would be useful to estimate the population genetics of this species precisely.

A recent study argued the confusing status of M. berdmorei types, and depicted the type locality in Bago, Myanmar (Garg et al. 2019). Besides, the species Microhyla fowleri was first described from Chieng Mai, Thailand, and later on synonymized and resurrected from M. berdmorei. Subsequently, the integrated taxonomic approaches neither elucidate the distinctiveness of M. fowleri, nor synonymized with M. berdmorei (Poyarkov et al. 2014; Garg et al. 2019). The most recent study illuminates the evolutionary position of M. fowleri under M. berdmorei group and noted as a complex of cryptic species (Gorin et al. 2020). Hence, the present study suggests the genetic information of both M. berdmorei and M. fowleri from their typical population which would be useful to discriminate the genetic distinctiveness and actual diversity in South and Southeast Asian countries.

Author's contribution

Conceptualization: SK, HTL; Sampling: LB, HD, LM; Data curation: SK, KT; Formal analysis: SK, KT, VK; Funding acquisition: KC, VK, HTL; Investigation: SK, HTL; Methodology: SK, LB; Project administration: KC, VK, HTL; Resources: KC; Software: SK, VK; Supervision: VK, KC, HTL; Validation: SK, VK; Visualization: SK, HTL; Writing – original draft: SK, HTL, KT; Writing – review and editing: SK, VK, HTL.

Acknowledgments

We are thankful to the Director of Zoological Survey of India (ZSI), Ministry of Environment, Forest and Climate Change (MoEF&CC), Govt. of India for providing necessary facilities and support for the study. The first author (SK) acknowledges the fellowship grant received from the Council of Scientific and Industrial Research (CSIR) Senior Research Associateship (Scientists’ Pool Scheme) Pool No. 9072-A. Authors would like to express their sincere thanks to the Chief Wildlife Warden, Department of Environment, Forest and Climate Change, Govt. of Mizoram for the collection permission, Mathipi Vabeiryureilai for helping in specimen collection.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are openly available in NCBI GenBank database at (https://www.ncbi.nlm.nih.gov) with the accession numbers (MT755953 and MW578946) which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Additional information

Funding

The research is funded by the Core Funding of Zoological Survey of India (ZSI), Kolkata, Ministry of Environment, Forest and Climate Change (MoEF&CC), New Delhi, National Mission for Himalayan Studies (NMHS) Uttarakhand, and Department of Biotechnology (DBT), New Delhi, India.