Mitochondrial genomes of Japanese Babina frogs ( Ranidae , Anura ) : unique gene arrangements and the phylogenetic position of genus Babina

Ryosuke Kakehashi, Atsushi Kurabayashi, Shohei Oumi, Seiki Katsuren, Masaki Hoso† and Masayuki Sumida* Institute for Amphibian Biology, Graduate School of Science, Hiroshima University, Hiroshima 739-8526, Japan Section of Agriculture and Forest, Amami City Government, Amami, Kagoshima 894-0048, Japan Okinawa Prefectural Institute of Health and Environment, Okinawa 901-1202, Japan Netherlands Centre for Biodiversity Naturalis, A107, Van Steenis, Einsteinweg 2, 2333 CC, Leiden, The Netherlands


INTRODUCTION
Genus Babina is a member of Ranidae, a large family of frogs currently comprising 10 species (see below).Thompson (1912) described the first two species in this genus-B.holsti and B. subaspera.These two original Babina species are famous for their "fifth finger": their first finger is a medially and ventrally fleshy sheath encasing a spine-like elongated 1st metacarpal, giving the appearance of five fingers on a hand.Because this character is unique among amphibians, these species have attracted the attention of many researchers in ecology, developmental biology, and taxonomy fields (e.g., Tokita and Iwai, 2010).B. holsti and B. subaspera are endemic to the small Ryukyu Islands, Okinawa and Amami, respectively.Recent environmental destruction has devastated populations of these species; and as a result, B. holsti and B. subaspera have been listed as a class B1 endangered species in the IUCN Red List of Threatened Species (http://www.iucnredlist.org)and designated as "natural monument species" in Okinawa and Kagoshima Prefectures, respectively.In addition, another Babina species, B. okinavana is listed as endangered.Prompt measures to conserve these species are awaited, and a recent population genetics study has been performed on B. subaspera using microsatellite markers (Iwai and Shoda-Kagaya, 2012), but no usable mitochondrial genomic marker on these species has been obtained.
Since the initial description of genus Babina (Thompson, 1912), several taxonomic changes have been made.Consequently, the original Babina species have been moved to genus Rana, a large ranid taxon, or occasionally to subgenus Babina of Rana (e.g., Maeda and Matsui, 1999).Recently, Frost et al. (2006) re-realized genus Babina and separated it from genus Rana.They also assigned 10 species to this genus: 2 species are the original species (designated here as Babina sensu stricto) while the others (Babina sensu lato) are from genus Nidirana (or subgenus Nidirana of genus Rana).Although the monophyly of these two Babina groups has been suggested (Kurabayashi et al., 2010;Kuraishi and Matsui, unpublished), there is no morphological or ecological foundation (e.g., synapomorphy in Babina) that justifies assigning them to the same genus.Further, which ranid group is the sister taxon of genus Babina is an unresolved matter in ranid phylogeny.Previous molecular phylogenetic studies have resulted in three distinct phylogenetic relationships for Babina.First, Stuart (2008) proposed that Babina is a sister taxon of the Glandirana, Hylarana, and Sanguirana clades based on two nuclear genes.By contrast, Cai et al. (2007) showed the close relationship of Babina with Rana and Lithobates based on mitochondrial ribosomal RNA gene data.From the same mitochondrial gene data, Kurabayashi et al. (2010) suggested that Babina is the sister group of Odorrana.Frost et al. (2006) and Pyron and Wiens (2011) also support the sister relationship of Babina and Odorrana.In every case, the statistical supports of Babina nodes were not high, and the Babina sensu stricto species was not used in these analyses (excluding Kurabayashi et al., 2010).Further, relatively small molecular data sets were applied in phylogenetic analyses in the previous studies regarding Babina.However, no consensus on the phylogenetic position of Babina has been reached.
To perform population genetics studies and resolve the phylogenetic problems related to Babina, mt genomic information seems to be a usable marker.Vertebrate mitochondrial (mt) DNA comprises closed circular molecules.
The gene arrangement of the animal mt genomes tends to be conserved among vertebrates, with all 37 genes and the CR organized relatively in the same gene order in almost all species--from teleost fish to eutherian mammals (Anderson et al., 1981;Roe et al., 1985;Tzeng et al., 1992;Boore, 1999).In anurans, the mt gene arrangement of archaeobatrachians (a paraphyletic assemblage of basal anuran families; e.g., Duellman, 2003) is identical to the typical arrangement in vertebrates (vertebrate-type arrangement, see Fig. 1; e.g., San Mauro et al., 2004).Yet the gene arrangement in the neobatrachian group (phylogenetically nested anurans including family Ranidae) is different.Four trns genes are commonly rearranged in the neobatrachian groups (neobatrachian-type arrangement) (Sumida et al., 2001;Zhang et al., 2005), and further gene arrangements have been reported in some natatanurans.Although a few convergent gene rearrangements have been observed in vertebrate mtDNA, many of these derived genomic rearrangements are clearly useful as phylogenetic markers (thus, synapomorphic characters) for the specification of monophyletic groups in natatanurans (e.g., Kurabayashi et al., 2006Kurabayashi et al., , 2008;;Zhang et al., 2009).Furthermore, and remarkably, Kurabayashi et al. (2010) showed unexpected gene rearrangement in ranid taxa including Babina and suggested the possible phylogenetic utility of the gene arrangement data of these taxa.
In the present study, we determined the complete mtDNA sequences of 3 endangered Babina species in Japan (2 sensu stricto species and 1 sensu lato species, B. okinavana) to identify mt genes with high sequence divergence and to check the unknown mt gene rearrangements.The former information can be applied in interpopulation genetic analyses of these species, and the latter will furnish possible phylogenetic markers to elucidate the phylogenetic problems related to this genus.Miochondrail Genomes of Babina species We also sequenced the mtDNA of the American bullfrog, L. catesbeianus.This species, of which complete mt genomic information has not been obtained, is very common and notorious as invader species in many countries.Furthermore, L. catesbeianus is one of the possible sister taxa of genus Babina (Cai et al., 2007).This paper describes the novel mt genomic structures and estimates the possible rearrangement pathway for the observed gene arrangements in these ranid taxa.Based on this information, we discuss the most related genus of Babina.

MATERIALS AND METHODS
Frog specimens Since the taxonomic revision of amphibians by Frost et al. (2006), the names and taxonomic ranks of anuran taxa, including those for a multitude of ranid groups, have been vertiginously changed.To prevent unnecessary confusion, we follow the newest system for anuran taxonomy from Frost (2011) in this paper.
The experiments in this study were performed using 3 frog specimens of Babina species B. holsti, B. subaspera, and B. okinavana collected in the Ryukyu islands, Japan (Okinawa Island, Amami Island, and Iriomote Island, respectively).Another ranid species, L. catesbeianus, kept by the strain maintenance team of the Institute for Amphibian Biology, Hiroshima University, was also used.PCR, subcloning, and sequencing Total DNA was extracted from a clipped toe of each frog specimen using the DNeasy Tissue Kit (QIAGEN).To determine the complete mt genomic sequences of 4 species, normal or long-and-accurate (LA) PCR using Ex-and LA-Taq (Takara-Bio.),respectively, was used according to the manufacturer's instructions (sequences of amplification primers are available in Supplementary Table S1).The resultant PCR fragments were electrophoresed on 0.7% agarose gel; then target DNAs were purified from excised pieces of gel using GenElute Agarose Spin Column (Sigma).The purified DNA fragments were used for sequencing and subcloning.
The primer-walking method was employed for mtDNA sequencing using an automated DNA sequencer (3100, ABI) with the BigDye Terminator Cycle Sequencing Kit (ABI) (sequences of the walking primers are available in Supplementary Table S1).A PCR fragment (correspond- Babina holsti and B. subaspera

Odorrana tormota
Neobatrachian-type arrangement ing to the CRs) containing long tandem repeats from which walking primers could not be constructed were subcloned into a pCR 2.1 vector using the TOPO TA Cloning Kit (Invitrogen).To determine the precise sequence of the long tandem repeats in the CR of B. okinavana, a series of deleted subclones was constructed from the single resultant clone by the Exonuclease III deletion method (Henikoff, 1987).
Genetic distance Sequence data were analyzed with GENETIX (Software Development).Genes were identified by comparison with corresponding gene sequences from other ranid species, and in the case of trns, by determination of secondary structural features.
The nucleotide sequences of 2 rrns, 13 mt protein-coding genes, and the partial CR were separately aligned using MUSCLE (Edgar, 2004).For the CR, it is known that multiple repetitive sequences exist in neobatrachians, such as Pelophylax nigromaculatus and Buergeria buergeri (Sumida et al., 2001;Sano et al., 2004).Thus, the alignable non-repetitive portions, about 1.6 kbp sequences from the 3′ end of the repeat to around CSB III were used for analyses as the partial CRs.Genetic distances of each gene and the partial CR were calculated among 6 ranid species (B. holsti, B. okinavana, B. subaspera, O. splendida, O. tormota, and L. catesbeianus).Gaps and missing sites were deleted in the calculations.The homogeneity of variance in the distance among different genes and the CR was tested using Bartlett's test.A Kruskal-Wallis test was performed to test for significant differences between the distance of each gene and the CR.

Molecular phylogenetic analyses
The phylogenetic analyses of the higher anuran groups were performed with 34 anurans, 8 other amphibians (6 caudates and 2 caecilians), and 4 other vertebrates (a mammal, bird, reptile, and teleost fish) for which complete mt sequences are available (see Fig. 1 legend).The nucleotide sequences of 2 ribosomal RNA and 13 mt protein-coding genes were separately aligned using MUSCLE (Edgar, 2004).To exclude gaps and ambiguous areas, the alignments were revised using Gblocks 0.91b software (Castresana, 2000) with the default parameters.Then, the corrected alignments from 15 genes were combined (11,345 bp).The phylogeny was analyzed with the combined data by maximum likelihood (ML) and Bayesian inference (BI) methods.The ML and BI analyses were respectively performed using TREEFINDER ver.Mar.2011 (Jobb, 2011) and MrBayes5D (Tanabe, 2008), a modified version of MrBayes3.1.2(Ronquist and Huelsenbeck, 2003).
Different substitution models were applied for each gene partitioned by codon positions in both analyses.The best-fit substitution model for each partition was estimated using Akaike's information criteria implemented in Kakusan4 ver.4.0 software (Tanabe, 2011).In the BI analysis, the following settings were applied: number of Markov chain Monte Carlo generations, 11 million; sampling frequency, 100; burn-in, 1,001.The robustness of the resultant ML tree was evaluated using bootstrap probabilities calculated from nonparametric bootstrap analyses with 1,000 pseudo-replications, and statistical support of the resultant BI tree was determined based on Bayesian posterior probability (BPP).

Mitochondrial genome organization of the Babina species
In this study, we determined the complete nucleotide sequences of the mt genomes from three Babina species (B.holsti, B. okinavana, and B. subaspera) and L. catesbeianus.Sizes were 19,113 bp, 19,959 bp, 18,525 bp, and 17,682 bp, respectively.As in other animals, the mt genomes of the four ranid species contain 37 genes consisting of 2 rrns, 22 trns, and 13 protein-coding genes, all of which are similar in length to their counterparts in other anurans (data not shown).In contrast, however, B. okinavana (sensu lato group) mtDNA has one pseudogene of mt trnH (ps-trnH, whose nucleotide similarity with the corresponding gene is 63.5%; Supplementary Fig. S1).Similarly, the mt genomes of B. holsti and B. subaspera (sensu stricto group) possess ps-trnH (similarities are 75.8% and 76.2% in B. holsti and B. subaspera, respectively; Supplementary Fig. S1).Additionally, the Babina sensu stricto group has ps-trnE (similarities are 72.1% and 70.0% in B. holsti and B. subaspera, respectively; Supplementary Fig. S1).
We also identified a major non-coding region downstream of cob in Babina species and in L. catesbeianus (see Fig. 1).Upon determining that this region contained TAS (termination-associated sequence), O H (Heavystrand replication origin), and CSB I-III (conserved sequence blocks I, II, and III), common features of the CR, we identified it as the CR.Noncoding sequences of 36-37 bp in length were also found.These were deemed to be of L-strand replication origin (O L ) because of their potential to form a hairpin structure characteristic of vertebrate O L (data not shown) and because their positions between trnN and trnC in the WANCY trns cluster correspond to the typical O L position (see Fig. 1).
As shown in Fig. 1, the mt gene arrangement of L. catesbeianus and B. okinavana is similar to that in most vertebrates and consistent with the typical neobatrachiantype arrangement having a LTPF trn-cluster.However, in the B. okinavana mt genome, there is an additional noncoding region of 551 bp between nad5 and nad6; and ps-trnH is additionally found downstream of the CR.In contrast to the former, the mt gene arrangement of the Babina sensu stricto group differs from the neobatrachian-type arrangement and shows the unique position of the ps-trnH-trnS (AGY)-nad5 gene region between the CR and the LTPF trn-cluster.The translations of nad5 from its original position (upstream of the CR) to downstream of the CR were reported in several neobatrachian mt genomes (e.g., dicroglossids, rhacophorids, mantellids, and the ranid Staurois latopalmatus; see Kurabayashi et al., 2010, Alam et al., 2010).However, the CR downstream position of trnS (AGY) and trnH or its pseudogene is found only in Babina and a few ranid species (S. latopalmatus and Odorrana spp.; see Kurabayashi et al., 2010).

Genetic distances of ranids
We compared the currently available nucleotide sequences of all mt genes (2 rrns, 13 protein-coding genes and a concatenated sequence of 22 trns) and the CR among 6 ranid species, including Babina and their related genera Odorrana and Lithobates.The genetic distance for each gene and the CR is shown in Fig. 2. Mean values of nucleotide divergence for each gene and the CR range from 10.6% in 12S rrns to 30.8% in atp8.Bartlett's test showed that the variances of nucleotide divergence are significantly heterogeneous (P < 0.05) among genes (and the CR).The Kruskal-Wallis test also rejected the uniformity of nucleotide divergences among the mt genes and the CR (P < 0.01).In the ranid taxa analyzed, atp8, nad5, and the CR showed the highest nucleotide divergence (30.8%, 26.7%, and 25.8%, respectively), and 12S showed the lowest nucleotide divergence (10.6%).

Phylogenetic analyses
We analyzed the phylogenetic relationship among anurans using the long mitochondrial sequence data (11,345 bp).The resultant Bayesian (BI) tree is shown in Fig. 3.The ML tree based on the same data set has an identical tree topology.In these trees, neobatrachians form a clade (BP = 100) in which the  ranoid families form a clade (BP = 100).Within the ranoid clade, ranid frogs form a monophyletic group (i.e., family Ranidae, BP = 100).These relationships are consistent with recent molecular phylogenetic studies on anurans (e.g., Igawa et al., 2008;Kurabayashi et al., 2010).In the ranid clade, genus Babina forms a monophyletic group (BP = 100).Genus Odorrana becomes the sister taxon of Babina.The BP support in ML is not high for this node, but BPP is high (BP = 67.8;BPP = 98.8).Genus Lithobates becomes the sister taxon to the Babina + Odorrana clade.

DISCUSSION
Nucleotide diversity in Babina mt genes One of the purposes of this study was to find mt genes with high nucleotide divergence, genes which will be useful as genetic markers among populations of Babina, including three endangered species (B. holsti, B. subaspera, and B. okinavana).Among the ranid taxa examined here (Fig. 2), 12S and 16S rrns and trns have small nucleotide divergence (10.6%, 14.2%, and 12.2%, respectively), while those of protein-coding genes and the CR are large.In particular, atp8, nad5, and the CR have large nucleotide divergence (30.8%, 26.7%, and 25.8%, respectively).Given the very short length of atp8 (168 bp), nad5 and the CR seem to be candidates for a suitable polymorphic marker for Babina and related ranid taxa.

2006
).Furthermore, the CR shows low divergence in salamanders and newts (Mueller, 2006;Kurabayashi et al., 2012), yet the CR is the third-most variable region in Odorrana (Kurabayashi et al., 2010).These findings emphasize that performing a preliminary survey of fastevolving genes and the CR is important for determining phylogenetic and/or population genetic markers for closely related amphibian taxa.

Phylogenetic relationship of Babina and related genera
As mentioned above, several alternative hypotheses have proposed a sister taxon for genus Babina.This study is the first attempt to tackle this problem using all mt sequence data, and our data revealed the sister relationship between Babina and Odorrana (Fig. 3).Genus Lithobates, which has been suggested to have a close relationship with Babina (Cai et al., 2007), becomes the sister taxon of the Babina + Odorrana clade in our analysis.These relationships recovered here are consistent with the results of three previous molecular phylogenetic studies (Frost et al., 2006;Kurabayashi et al., 2010;Pyron and Wiens, 2011).Although the BP support of the Babina and Odorrana clade is not sufficient, BPP is high (BP = 67.8;BPP = 98.8).Furthermore, a unique mt gene arrangement seems to support the mono-phyly of Babina and Odorrana (see below).

Evolutionary pathway and systematic implications of mt gene rearrangement
The present study shows that the mt gene arrangements of Babina species differ from those of typical neobatrachians (Fig. 1), and Babina sensu lato and groups possess different gene orders.Furthermore, the unique position of trnH and its pseudogene (downstream of the CR) is observed in all Babina and Odorrana (sister taxon of Babina) species examined so far.By contrast, mtDNA of Lithobates, the sister taxon of the Babina + Odorrana clade, possesses the typical neobatrachian-type gene arrangement.These observations seem to indicate that gene rearrangement events occurred and that the unique trnH (and its pseudogene) emerged in an ancestral lineage which led to Babina and Odorrana after the splitting of Lihobates (see Fig. 4).
Animal mt genomes lack introns and intergenic spacers, and there are no multicopy genes (with a few exceptions; e.g., Kurabayashi et al., 2008).In such genomes, the occurrence of unregulated gene rearrangement will destroy an essential single-copy gene.Thus, gene rearrangements in animal mt genomes are generally explained by the "Duplication and Deletion model" (Moritz and Brown, 1986;Boore, 1999;San Mauro et al., Fig. 4. Mitochondrial gene rearrangement pathway of Babina and related taxa.Possible rearrangement lineages and putative duplication or deletion events in each lineage are illustrated.ML and BI tree topologies (Fig. 3) are used here."X" and "Ψ" under the gene boxes indicate "deleted gene" and "pseudogene" in extant species, respectively.Gene names and the other abbreviations are the same as those in Fig. 1.

2006
).In this model, a duplication including multiple genes is caused by a replication slippage or DNA recombination (see Kurabayashi et al., 2008).Then, one of copied genes is deleted through an accumulation of natural mutations.The duplicated CRs and several remnants of the duplicated genes (i.e., ps-trns) observed in ranoids, including the Babina species examined here, support the occurrence of duplication and deletion events.According to the gene arrangement mechanism, there seems to be only one duplication event that can explain the divergent gene arrangements between Babina sensu lato and sensu strict groups as well as the unique trnH position shared by Babina and Odorrana as follows.First, the gene region from trnH to trnE was duplicated in the common ancestral lineage of Babina and Odorrana (Fig. 4I).From the duplicated genomic condition, the upstream trnH gene and the downstream trnS (AGY) to trnE region would have been deleted in the lineage which previously led to Odorrana; then, the present Odorrana mt gene arrangement occurred (Fig. 4II).In the common ancestor lineage of Babina, the duplicated genomic condition was maintained in the last common ancestor of the sensu stricto and lato groups.In the lineage leading to B. okinavana (sensu lato group), the downstream copied genes were deleted, and only trnH remained as the pseudogene (Fig. 4III).Finally, in the lineage leading to B. holsti and B. subaspera (sensu strict group), the original trnS (AGY), nad5, and trnE were deleted, and the copied trnH remained as the pseudogene (Fig. 4IV).Consequently, these species have the translocated trnS (AGY), nad5, and trnE genes as well as an additional trnH pseudogene.
If we hypothesize a single duplication event, the gene rearrangements observed in Babina and Odorrana can easily be explained by only several deletion events.In this context, the derived position of trnH and ps-trnH (CR downstream) can be regarded as a synapomorphic character of these genera (Fig. 4II-IV).Furthermore, the additional rearranged gene positions (trnS (AGY)-nad5-trnE between the CR and LTPF trns cluster) found in B. holsti and B. subaspera will be synapomorphic with the Babina sensu stricto group (Fig. 4IV).
The mt genome arrangement is thought to be a usable phylogenetic marker, especially for the taxa in which it is difficult for phylogeny to be resolved by normal gene sequence data (e.g., Boore andBrown, 1994, 1998;Boore et al., 1995;Macey et al., 2000;Kurabayashi and Ueshima, 2000;Inoue et al., 2003;Mabuchi et al., 2004).The unique trnH and ps-trnH position shared by Babina and Odorrana seems to be one such example.There remain many unresolved relationships among ranid genera (e.g., the position of genera Glandirana, Amolops, Hylarana; Cai et al., 2007;Che et al., 2007;Pyron and Wiens, 2011), but the mt genomic information may be a powerful tool for resolving these problems.

Fig. 1 .
Fig. 1.Comparison of mt gene arrangements among ranid taxa.The transcriptional direction of H-strand encoding genes and the upstream and downstream notations used in this paper are shown by a closed arrow and closed arrowheads, respectively.The H-and L-strand encoded genes are denoted at the top and bottom of each gene box, respectively.The sizes of the boxes do not reflect actual gene length.Closed arrows show the rearranged genes and infer the evolutionary directions of the rearrangements (see the text).Gray boxes represent the rearranged genes from the neobatrachian-type positions.The derived trnH position shared by Odorrana tormota and Babina species is also shown.Transfer RNA genes (trns) are designated by single-letter amino acid codes.L1, L2, S1, and S2 indicate trns for Leu (UUR), Leu (CUN), Ser (AGY), and Ser (UCN), respectively.Trn boxes with "ψ" indicate the pseudogenes.The control region is abbreviated as CR.O L indicates the region of the light-strand replication origin including a typical stem-loop structure.Other genes are abbreviated as follows: 12S and 16S, 12S and 16S ribosomal RNA; cox1-3, cytochrome c oxidase subunits 1-3; cob, cytochrome b; nad1-6 and 4L, NADH dehydrogenase subunits 1-6 and 4L. cox3

Fig. 2 .
Fig. 2. Box plot of the nucleotide divergence of each mt gene and the CR in ranid frogs.Each box plot shows the currently available nucleotide divergence of each mt gene or CR among 6 ranid frogs (3 Babina, 2 Odorrana, and 1 Lithobates species) (see Fig. 3).Abbreviations of gene names are provided.Box plots have a line representing median divergence, and the surrounding box contains the middle 50% (25-75%) of the data, with whiskers showing the 1.5 inter-quartile range.Pluses indicate average nucleotide divergence of the mt genes (and the CR).The genes and the CR are arranged from the smallest average divergence to the largest average divergence.Open circles indicate differences of Babina (B.holsti vs. B. okinavana vs. B. subaspera).
Common ancestral lineage of Babina and Odorrana