Complete mitochondrial genomes of three fairy shrimps from snowmelt pools in Japan
BMC Zoology volume 7, Article number: 11 (2022)
Fairy shrimps belong to order Anostraca, class Branchiopoda, subphylum Crustacea, and phylum Arthropoda. Three fairy shrimp species (Eubranchipus uchidai, E. asanumai, and E. hatanakai) that inhabit snowmelt pools are currently known in Japan. Whole mitochondrial genomes are useful genetic information for conducting phylogenetic analyses. Mitochondrial genome sequences for Branchiopoda members are gradually being collated.
Six whole mitochondrial genomes from the three Eubranchipus species are presented here. Eubranchipus species share the anostracan pattern of gene arrangement in their mitochondrial genomes. The mitochondrial genomes of the Eubranchipus species have a higher GC content than those of other anostracans. Accelerated substitution rates in the lineage of Eubranchipus species were observed.
This study is the first to obtain whole mitochondrial genomes for Far Eastern Eubranchipus species. We show that the nucleotide sequences of cytochrome oxidase subunit I and the 16S ribosomal RNA of E. asanumai presented in a previous study were nuclear mitochondrial DNA segments. Higher GC contents and accelerated substitution rates are specific characteristics of the mitochondrial genomes of Far Eastern Eubranchipus. The results will be useful for further investigations of the evolution of Anostraca as well as Branchiopoda.
The fairy shrimp genus Eubranchipus belongs to family Chirocephalidae, order Anostraca, class Branchiopoda, subphylum Crustacea, and phylum Arthropoda. These shrimps occur in temporary pools formed by snowmelt in the forest groves in northern Japan . They hatch from resting eggs when the snowmelt water appears in early spring and mature in the pool and lay eggs just before the water dries up in late spring. Thus far, studies of these shrimps have been limited [1,2,3,4], because they are found only in inconspicuous locations such as forest bushes and only during a short period . In 2018, Takahashi et al.  described three new species of Eubranchipus for the first time in 62 years, from Far East Asia. They also reviewed the morphological ambiguity of the earlier description of E. uchidai and updated the molecular systematics with the new Asian taxa.
The animal mitochondrial genome is a small, extrachromosomal genome. It has a simple conserved structure that is approximately 16 kb long, with 13 protein-coding genes (COX1, COX2, ATP8, ATP6, COX3, ND3, ND5, ND4, ND4L, ND6, CYTB, ND1, and ND2), two ribosomal RNA genes (16S and 12S), 22 tRNA genes, and a noncoding region known as a control region . Utilization of the complete mitochondrial genome for molecular phylogenetic analysis has two main advantages: one is that more phylogenetic information can be obtained than when using partial sequences, and the other is that NUMTs are avoided. In this study, we determined the complete mitochondrial genomes of three Eubranchipus species in Japan and compared these with those of other anostracans.
Although the animal mitochondrial genome is conserved across the animal kingdom, rearrangements in gene order may occur. In the Branchiopoda, three gene order patterns have been reported [7, 8]. The first is the “ancestral pancrustacean pattern” , which is shared by Daphnia and Triops species. The second, based on the ancestral pancrustacean pattern, is the “anostracan pattern” , which translocates (trnI + trnQ) from between the control region and trnM to between trnW and trnC, together with an inversion of trnI. Artemia, Phallocryptus, and Streptocephalus species share this pattern. The third pattern is observed in Branchinella kugenumaensis, which underwent a further large inversion of the block (trnM + ND2 + trnW + trnI) from the anostracan pattern .
GC content differences have been observed in branchiopod mitochondrial genomes . Anostraca (36.8%) and Onychocaudata (Limnadia lenticularis + Cladocera) (35.0%) have a significantly higher GC content than Notostraca (30.5%), which is explained by the preferential AT to GC substitution bias during the evolution of the Anostraca and Onychocaudata lineages.
In this study, we report on the sequencing and analysis of six mitochondrial genomes from three Eubranchipus species of Chirocephalidae, Anostraca. The newly obtained mitochondrial genomes were analyzed alongside all currently available mitochondrial genomes of Branchiopoda members to gain insight into the evolution of this class of crustaceans.
Complete mitochondrial genomes
We used two individuals of each species (Eubranchipus hatanakai Takahashi & Hamasaki, in Takahashi et al., 2018 [Sample IDs: Eh1, Eh6], Eubranchipus uchidai (Kikuchi, 1957) [Sample IDs: Eu17, Eu36], and Eubranchipus asanumai Takahashi, in Takahashi et al., 2018 [Sample IDs: Ea3, Ea4]), to identify the complete mitochondrial genomes. The estimated respective DINs and concentrations of the extracted genomic DNAs were as follows: Eh1 (7.4, 39.7 ng/μL), Eh6 (7.3, 57.8 ng/μL), Eu17 (7.6, 53.5 ng/μL), Eu36 (7.8, 59.1 ng/μL), Ea3 (7.8, 24.2 ng/μL), and Ea4 (7.5, 34.2 ng/μL). In total, 4.4–7.3 Gb of corrected data were obtained from 5.0–7.8 Gb of raw data (i.e., 90.2–93.7%) for each fastq dataset using the Pollux 1.0.2 program  (Additional file 1: Table S1). The estimated k-mers are shown in Additional file 1: Table S1. De novo genome sequence assembly was performed using the Ray 2.1.0 program , and 329,762–491,652 assembled sequences were obtained. Partial mitochondrial DNA sequences were obtained from these assembled sequences via tBLASTn searches  ([Eh1]: scaffold-325,202 (7647 bp) containing ND2-COX1-COX2-ATP8-ATP6-COX3-ND3, scaffold-21 (4137 bp) containing ND5-ND4-ND4L-ND6-CYTB, scaffold-42 (3350 bp) containing ND1; [Eh6]: scaffold-0 (7649 bp) containing ND2-COX1-COX2-ATP8-ATP6-COX3-ND3, scaffold-328,332 (8201 bp) containing ND5-ND4-ND4L-ND6-CYTB-ND1; [Eu17]: scaffold-96 (3279 bp) containing ND2-COX1, scaffold-31 (3483 bp) containing COX2-ATP8-ATP6-COX3-ND3, scaffold-12 (7853 bp) containing ND5-ND4-ND4L-ND6-CYTB-ND1; [Eu36]: scaffold-34 (6444 bp) containing ND2-COX1-COX2-ATP8-ATP6-COX3-ND3, scaffold-12 (7889 bp) containing ND5-ND4-ND4L-ND6-CYTB-ND1; [Ea3]: scaffold-157 (5954 bp) containing ND2-COX1-COX2-ATP8-ATP6-COX3-ND3, scaffold-12 (8281 bp) containing ND5-ND4-ND4L-ND6-CYTB-ND1; [Ea4]: scaffold-403,293 (6440 bp) containing ND2-COX1-COX2-ATP8-ATP6-COX3-ND3, scaffold-4 (8198 bp) containing ND5-ND4-ND4L-ND6-CYTB-ND1) (Additional file 2: Selected scaffold sequences). We then performed mitochondrial genome sequence assembly in the NOVOPlasty 3.2 program  using the longest mitochondrial DNA sequence obtained as the seed for each dataset (Eh6: scaffold-328,332, Eu17: scaffold-12, Eu36: scaffold-12, Ea3: scaffold-12, Ea4: scaffold-4). Because the complete mitochondrial genome sequence of Eh1 could not be obtained using scaffold-325,202 as the seed in the NOVOPlasty assembly, we used a combined sequence comprising two scaffolds (scaffold-325,202 and scaffold-21) as the seed. Circularized assembly sequences were obtained for all samples. Next, we remapped the fastq reads to the assembled mitochondrial genome sequence to check the sequence read depth (Additional file 3: Fig. S1). We thus obtained complete mitochondrial genomes for E. hatanakai (Eh1 [17,006 bp], Eh6 [17,006 bp]), E. uchidai (Eu17 [15,795 bp], Eu36 [15,795 bp]), and E. asanumai (Ea3 [17,503 bp], Ea4 [17,503 bp]) (Fig. 1). We found three variant sites (172AG, 9010 AC, and 14234CT) in Eh1, one (4804AG) in Eh6, and one (13191AT) in Eu17. Lists of the annotated loci of the six sequences are given in Additional file 4: Table S2.
These Eubranchipus species share the anostracan gene arrangement pattern in their mitochondrial genomes . The nucleotide sequence data obtained in this study were deposited in the DDBJ/EMBL/GenBank International Nucleotide Sequence Database (accession numbers: LC633437–LC633442).
COX1 and 16S sequences
We constructed phylogenetic trees for the COX1 (Fig. 2) and 16S (Fig. 3) sequences of Chirocephalidae species. For E. uchidai, the COX1 sequences of Eu17 and Eu36 were identical to that of the A type of E. uchidai in the DNA databank (LC314408.1, ), but we observed one nucleotide difference in the 16S sequence (LC314409.1, ). For E. hatanakai, the COX1 sequence of Eh6 was identical to that of E. hatanakai in the DNA databank (LC314402.1, ), while there was one nucleotide difference in that of Eh1. The 16S sequences of both Eh1 and Eh6, however, were identical to that of E. hatanakai in the DNA databank (LC314403.1, ).
We observed distinct differences for E. asanumai, however: there were 44 nucleotide differences (0.090 differences per site) between the COX1 sequences of Ea3 and Ea4 and that of E. asanumai in the DNA databank (LC314404.1, ). In addition, there were 13 nucleotide differences (0.027 differences per site) between the 16S sequences of Ea3 and Ea4 and that of E. asanumai in the DNA databank (LC314405.1, ).
We extracted genomic DNA and total RNA samples from one individual of E. asanumai (ST-M2) and conducted PCR and RT-PCR analyses on them using the same primer sets. The COX1 sequence amplified from the genomic DNA sample was identical to that of the E. asanumai sequence in the DNA databank (LC314404.1, ) (Fig. 2), and the 16S sequence had only one nucleotide different from that of the E. asanumai sequence in the DNA databank (LC314405.1, ) (Fig. 3). In contrast, both the COX1 sequence and the 16S sequence amplified from the total RNA sample had only one nucleotide different from Ea3 and Ea4 (Figs. 2 and 3).
Nucleotide differences in mitochondrial DNA within and between Eubranchipus species
We observed four nucleotide differences between the two mitochondrial DNA sequences (Eh1 and Eh6) of E. hatanakai: 172AG nonsynonymous in COX1, 3368CG nonsynonymous in COX3, 7766AG synonymous in ND4, and 8628AG nonsynonymous in ND6. We also observed four nucleotide differences between the two mitochondrial DNA sequences (Eu17 and Eu36) of E. uchidai: 2675AG nonsynonymous in ATP6, 3150CT nonsynonymous in ATP6, 8010AG nonsynonymous in ND4L, and 13,703 AC in the region between rrnS and trnM (probably the D-loop region). However, we found 20 nucleotide differences between the two mitochondrial DNA sequences (Ea3 and Ea4) of E. asanumai, 10 of which were as follows: 1727AG synonymous in COX2, 2959CT synonymous in ATP6, 3243GA synonymous in COX3, 3324AG synonymous in COX3, 5484TC synonymous in ND5, 6989TC nonsynonymous in ND4, 7849AG synonymous in ND4L, 7923AG synonymous in ND4L, 8199CT in trnP, and 16998AG synonymous in ND2, and the remaining 10 were 10542AG, 10974CT, 12940AG, 13547CT, 13793TA, 14906TC, 15843TG, 15899TC, 15957CT, and 16195GA in the region between rrnS and trnM (probably the D-loop region). There were no insertion or deletion differences between the two sequences of each species.
The genetic distances for each locus among the Eubranchipus species are listed in Additional file 5: Table S3. The overall genetic distance between E. hatanakai and E. uchidai was 0.192, ranging from 0.033 for trnH to 0.340 for ATP8. The overall genetic distance between E. hatanakai and E. asanumai was 0.189, ranging from 0.033 for trnL1 to 0.331 for ATP8. In addition, the overall genetic distance between E. uchidai and E. asanumai was 0.041, ranging from 0 for 10 tRNA loci (trnL2, trnK, trnD, trnG, trnR, trnS1, trnH, trnT, trnI, and trnC) to 0.070 for trnV. The overall genetic distances between E. hatanakai and E. grubii, between E. uchidai and E. grubii, and between E. asanumai and E. grubii were 0.269, 0.260, and 0.260, respectively.
GC content of mitochondrial genomes of Branchiopoda species
We estimated the GC contents of the whole mitochondrial genomes, 13 protein-coding genes with first, second, and third codon positions, and two rRNA regions for Branchiopoda (Tables 1 and 2). The GC content of the whole mitochondrial genomes of the Anostraca, Artemiidae, ranged from 35.49% for Artemia sinica to 37.55% for A. urmiana; that of the Chirocephalidae ranged from 33.04% in E. grubii to 38.96% in E. hatanakai; and that of the Thamnocephalidae and Streptocephalidae ranged from 31.83% in B. kugenumaensis (Japan) to 35.40% in Streptocephalus sirindhornae. These were significantly different among taxa (Kruskal–Wallis, p < 0.05). Multiple comparisons indicated that Artemiidae had a significantly higher GC content than the Thamnocephalidae and Streptocephalidae (Steel–Dwass, p < 0.05), but other comparisons did not indicate significant differences. We obtained similar results for the third codon position and two rRNAs. Multiple comparisons of the third codon position (Steel–Dwass, p < 0.05) and the two rRNAs (Steel–Dwass, p < 0.05) also indicated that the Artemiidae had a significantly higher GC content than the Thamnocephalidae and Streptocephalidae, while other comparisons did not indicate any significant differences.
The GC content of the second codon position of the Artemiidae ranged from 36.40% in A. franciscana to 36.96% in A. tibetiana; that of the Chirocephalidae ranged from 37.09% in E. grubii to 38.22% in E. hatanakai; and that of the Thamnocephalidae and Streptocephalidae ranged from 35.59% in B. kugenumaensis (Japan) to 36.83% in Phallocryptus tserensodnomi, and these values differed significantly among taxa (Kruskal–Wallis, p < 0.05). Multiple comparisons indicated that the Chirocephalidae had a significantly higher GC content than the Thamnocephalidae and Streptocephalidae (Steel–Dwass, p < 0.05), while other comparisons did not indicate any significant differences.
Phylogenetic tree of Branchiopoda species
In addition to our three mitochondrial genome sequences from E. hatanakai (Eh1), E. uchidai (Eu36), and E. asanumai (Ea3), we used mitochondrial DNA sequences from ten other Anostraca species, six Notostraca species, and eight Diplostraca species (Table 1). Selected substitution models are listed in Additional file 6: Table S4. For the ML approach , we used the models selected by AIC . The corrected AIC (AICc)  selected the same models as were selected by AIC. For the BS approach , we used the models selected by the BIC . For the NJ method , we selected the Tamura and Nei  model with a γ correction (α = 0.41) as the best available model.
We constructed the NJ tree of Branchiopoda species using the 13 protein-coding mitochondrial genes (Fig. 4). The divergent orders in Anostraca were Artemiidae, Chirocephalidae, Thamnocephalidae and Streptocephalidae. P. tserensodnomi and B. kugenumaensis (Thamnocephalidae) did not form a monophyletic cluster, whereas Cladocera and Spinicaudata (Diplostraca) did form a cluster. Two B. kugenumaensis mitochondrial genome sequences had been deposited in the DNA database, one from Japan and the other from China, and there are large nucleotide differences between them.
We also reconstructed the BS tree (Additional file 7: Fig. S2) and the ML tree (Additional file 8: Fig. S3). The topologies of the NJ (Fig. 4) and BS (Additional file 7: Fig. S2) trees were identical, but there was one difference between the topology of the ML tree and those of the NJ and BS trees: Daphnia galeata and D. laevis formed a cluster with a 64% bootstrap value in the ML tree (Additional file 8: Fig. S3). We compared the two topologies (NJ and BS versus ML) using the likelihood method (Table 3). The likelihood of the topology of the ML tree was higher than that of the NJ and BS trees, but according to three tests, this difference was not statistically significant.
In all three phylogenetic trees (Fig. 4, Additional file 7: Fig. S2, and Additional file 8: Fig. S3), the branches of the Diplostraca, Artemiidae, and Chirocephalidae clusters appeared longer than those of the Notostraca, Thamnocephalidae, and Streptocephalidae clusters. We investigated the constancy of the substitution rate and the variation among lineages (Table 4, Additional file 9: Fig. S4). The “no clock” model means that substitution rates are entirely free to vary from branch to branch, whereas the “global clock” model means that all branches have the same substitution rate. There was a significant difference between the no clock and global clock models, suggesting variations in the substitution rates among lineages. The “local clock1” model assumes that the Artemiidae, Chirocephalidae, and Diplostraca lineages have higher substitution rates than those of the Notostraca, Thamnocephalidae, and Streptocephalidae. There was a significant difference between the global clock and local clock1 models. In the “local clock2” model, the substitution rates of the Artemiidae, Chirocephalidae, and Diplostraca were assumed to be independent, and this was also significantly different from the global clock model. This implies that the acceleration of the substitution rate occurred independently in the Artemiidae, Chirocephalidae, and Diplostraca lineages.
We estimated the overall ratios of nonsynonymous substitutions per nonsynonymous site to synonymous substitutions per synonymous site (ω) to be 0.109, 0.106, 0.121, 0.166, and 0.138 for the Artemiidae, Chirocephalidae, Thamnocephalidae and Streptocephalidae, Diplostraca, and Notostraca clusters, respectively.
Nuclear mitochondrial DNA segments in previous annotations of Eubranchipus asanumai
After the study conducted by Takahashi et al.  was published, different types of sequence of both COX1 and 16S were obtained from different E. asanumai individuals based on PCR and sequencing using genomic DNA. We initially considered the possibility of the presence of another species in the Shiretoko area. However, based on morphological observations, this was not possible . In this study, therefore, we determined the complete mitochondrial genome. This is thus the first study to reveal the whole mitochondrial genomes of Far Eastern Eubranchipus species. At the same time, we performed RT-PCR using extracted total RNA. Because the sequences amplified using this technique would be functional, these sequences could not be NUMTs. It is regrettable that (for unknown reasons) the COX1 and 16S sequences of E. asanumai presented by Takahashi et al.  were incorrect. In addition, separate from the complete mitochondrial genome sequences of E. asanumai, we observed a 244 bp fragment from the Ea3 Ray genome assembly data (scaffold-278,847) and a 387 bp fragment from the Ea4 Ray genome assembly data (scaffold-236,301) (Additional file 2: Selected scaffold sequences). They were almost identical (one nucleotide difference) to the sequence of LC314404.1 (NUMT of COXI), which was published by Takahashi et al. We also identified a 139 bp fragment from the Ea3 Ray genome assembly data (scaffold-458,031), which was almost identical (one nucleotide difference) to the sequence of LC314405.1 (NUMT of 16S), also published by Takahashi et al. (Additional file 2: Selected scaffold sequences). Further, we identified a fragment (around 350 bp) from the Eu17 (scaffold-250,790) and Eu36 (scaffold-122,170) Ray genome assembly data for E. uchidai (Additional file 2: Selected scaffold sequences) that was identical to the nucleotide sequence of LC314404.1 (NUMT of COX1). We were not able to observe the 3′ end of the nucleotide sequence, which contains the reverse primer khCOI-R  region. We were therefore unable to amplify the NUMT of COX1 using the primer set of Takahashi et al. . It is likely that the NUMT of COX1 occurred in the common ancestor of E. asanumai and E. uchidai, and that the 3′ part of the sequence was lost in the E. uchidai genome. We were also not able to observe the NUMT of 16S in the Ray genome assembly data for E. uchidai (Eu17 and Eu36), or that of COX1 or 16S in the data for E. hatanakai (Eh1 and Eh6).
Although Takahashi et al.  checked that there were no gaps in the COX1 or 16S rRNA compared with other species and that all nucleotide changes, which are subject to functional constraints, in COX1 were synonymous, this was not sufficient, since further analysis (such as RT-PCR) would have been needed to avoid NUMTs. The entries in the DDBJ/EMBL/GenBank International Nucleotide Sequence Database (LC314404.2 for COX1 and LC314405.2 for 16S) have been corrected accordingly.
Although both of the corrected nucleotide sequences for E. asanumai were more similar to those of E. uchidai than the previous sequences, they nevertheless formed distinct clusters in the phylogenetic trees (Figs. 2 and 3). Nucleotide differences between E. asanumai and E. uchidai ranged from 5.5 to 6.0% for COX1 (Additional file 10: Table S5) and from 1.8 to 2.0% for 16S (Additional file 11: Table S6).
The high GC content of the mitochondrial genomes of Eubranchipus species
Luchetti et al.  showed that the mitochondrial genomes of the Anostraca and Onychocaudata (L. lenticularis + Cladocera) have a significantly higher GC content than those of Notostraca species. According to these authors, this can be explained by a preferential AT to GC substitution bias during the evolution of the Anostraca and Onychocaudata lineages. In this study, we observed differences in substitution bias among Anostraca species. The Artemiidae and Chirocephalidae tend to have a higher GC content than other Anostraca species. Although Eubranchipus species have a similar or higher GC content than Artemia species, we did not find any significant differences between the Chirocephalidae and Thamnocephalidae + Streptocephalidae (Table 2). For this reason, E. grubii, which has a relatively low GC content, is included in the Chirocephalidae. If E. grubii were eliminated from comparisons, the number of species of Chirocephalidae would decrease, and statistical tests could not be performed. Thus, we suggest that Far Eastern Eubranchipus species also have a relatively high GC content, equivalent to that of Artemia species. We infer that the higher GC contents of Eubranchipus and Artemia species than that of the Thamnocephalidae and Streptocephalidae are also caused by a preferential AT to GC substitution bias. In this study, however, the number of Eubranchipus species analyzed is limited. Further investigation, including other Eubranchipus species, is thus needed to clarify the differences in substitution bias among Anostraca species.
Phylogenetic tree of the mitochondrial genome data for Branchiopoda
The current knowledge on the phylogenetic relationships of Branchiopoda is that Anostraca diverged first, followed by Notostraca and Diplostraca (Cladocera and Spinicaudata) [9, 40,41,42,43,44]. The topologies we obtained (Fig. 4, Additional file 7: Fig. S2, and Additional file 8: Fig. S3) support these relationships.
Luchetti et al.  indicated that the Anostraca and Onychocaudata (Diplostraca) have a significantly higher substitution rate than the Notostraca, similar to what we found (Table 4). Further, we observed differences in substitution rates among the Anostraca: accelerated substitution rates within the Artemiidae and Chirocephalidae lineages. In this study, we assumed an acceleration of the Chirocephalidae lineage, including E. grubii (Additional file 9: Fig. S4); however, there is also the possibility that the acceleration occurred in the common ancestor of the Far Eastern Eubranchipus species. Because we used only one species, E. grubii, as out of Far Eastern Eubranchipus species in the study, further investigation including other Eubranchipus species is needed to clarify whether the acceleration occurred in the common ancestor of the Chirocephalidae or that of the Far Eastern Eubranchipus lineage. It is likely that substitution rate accelerations in the Artemiidae and Chirocephalidae lineages occurred independently, because the Artemiidae and Chirocephalidae do not form a monophyletic cluster (Fig. 4). This phylogenetic relationship is similar to that reported in previous studies [7, 45,46,47].
Accelerations of substitution rates in the Artemiidae, Chirocephalidae, and Diplostraca lineages are not associated with increases in nonsynonymous substitutions, because the overall ω values of the Artemiidae (0.109), Chirocephalidae (0.106), and Diplostraca (0.166) do not differ from those of the Thamnocephalidae and Streptocephalidae (0.121) or the Notostraca (0.138). Lower ω values (< 1) were commonly reported in previous studies [9, 18, 21, 22, 28]. This is due to selective constraints acting on the genes of mitochondrial genomes.
For the Artemiidae, it has been suggested that halophilic habits may be correlated with accelerated substitution rates . This is not the case for the Chirocephalidae because members of this family inhabit freshwater bodies. Other potential reasons for accelerated substitution rates such as body size, temperature, generation time, and population size have been considered [49,50,51,52]. The typical body length of Anostraca is 1–5 cm , and the three Eubranchipus species are within this range [1, 2, 4]. Generation time is difficult to estimate for Branchiopoda. Most branchiopod eggs are drought-resistant and can remain dormant for decades under anoxic conditions, and during pool inundation, only a fraction of each egg clutch hatches . Eubranchipus species also have this “bet-hedging” strategy. It is thus not straightforward to estimate or compare generation times among Anostraca species.
We infer that the effect of population size is the most plausible reason for the accelerated substitution rate of Eubranchipus species. It has been suggested that evolution occurs rapidly in small populations . Habitat for Eubranchipus species is rather limited in Japan. Presumably, a few eggs originally came from another location to the current habitats a long time ago, and acceleration of the substitution rate occurred in these small populations, after which the population sizes increased rapidly over several generations. The genetic diversity of these species would therefore also be expected to be small. In this study, we used mitochondrial genome data for our phylogenetic analyses, but it is also necessary to conduct investigations using nuclear genome data, to better understand the substitution rates and genetic diversity of Eubranchipus species.
In the current study, we present six new mitochondrial genome sequences from three Eubranchipus species. These are the first reports of the entire mitochondrial genome sequences of these Eubranchipus species. We show that these species shared the anostracan pattern of gene arrangements in their mitochondrial genomes. We observed a higher GC substitution bias in Eubranchipus than in other Anostraca species. We noted that the COX1 and 16S sequences presented in Takahashi et al.  were NUMTs, and we have corrected these in the present study. We also conducted a phylogenetic analysis of Branchiopoda species using mitochondrial genome data. We observed accelerations of substitution rates within the lineages of Eubranchipus species. Higher GC content and accelerated substitution rates are the specific characteristics of the mitochondrial genome of Eubranchipus.
Eubranchipus uchidai specimens were collected from a temporary snowmelt pool in Ishikari, Hokkaido, on 29 April 2017. Eubranchipus asanumai specimens were collected from a temporary snowmelt pool in Shiretoko, Hokkaido, on 18 May 2018. Eubranchipus hatanakai specimens were collected from a temporary snowmelt pool in Yuza, Yamagata, on 10 April 2018.
Genomic DNA extraction
Two E. hatanakai individuals (sample IDs: Eh1 [male], Eh6 [male]), two E. uchidai individuals (sample IDs: Eu17 [female], Eu36 [female]), and two E. asanumai individuals (sample IDs: Ea3 [male], Ea4 [male]) were used for genomic DNA extraction and sequencing. Each individual was homogenized using a disposable BioMasher II homogenizer (Nippi, Tokyo, Japan). The genomic DNA was extracted using the conventional sodium dodecyl sulfate lysis and phenol–chloroform method and RNase A (Sigma, St. Louis, MO, USA) was used to digest any contaminated RNA. Then, to remove RNase A proteins, the genomic DNA was purified using a NucleoSpin gDNA Clean-up Kit (TaKaRa Bio, Kusatsu, Japan). The quality of the extracted genomic DNA was checked using an Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA, USA).
A sequencing library was constructed using the Ion Xpress Plus Fragment Library Kit (Thermo Fisher Scientific, Waltham, MA, USA), and sequencing was performed using the Ion Proton System (Thermo Fisher Scientific). Approximately 6 Gb were sequenced per sample. To correct homopolymer errors in the fastq data, we used Pollux 1.0.2  with option “-k 31”. The optimal k-mer length for the corrected fastq data was then estimated using KmerGenie 1.7016 .
The de novo genome sequence assembly was performed using Ray 2.1.0  with option “-k”. The k-mer values used are shown in Additional file 1: Table S1. Next, standalone BLAST (tBLASTn) searches  were performed using 13 amino acid sequences from the mitochondrial genome sequence data of E. grubii (MT410793.1) as queries against the Ray assembly data. The longest mitochondrial DNA sequence obtained for each dataset was used as the seed for the mitochondrial genome sequence assembly in NOVOPlasty 3.2 . To check the sequence read depth, read remapping to the assembled mitochondrial genome sequence was performed using BWA 0.7.12  with the “aln” option and SAMtools 0.1.19 . The remapping was performed using the pipeline of the Management and Analysis System for Enormous Reads (Maser) of the Platform Project for Supporting Drug Discovery and Life Science Research (Platform for Drug Discovery, Informatics, and Structural Life Science) . The Integrative Genomics Viewer 2.8.13  was used to visualize mapped data. Mapped sites were counted using igvtools .
Gene annotations were performed using the MITOS WebServer  and manually corrected by following the annotation of the complete mitochondrial genome of E. grubii (MT410793.1).
Reanalysis of mitochondrial DNA sequences of E. asanumai
We used one E. asanumai individual (ST-M2) from the same locality as the paratypes. It was collected on 28 May 2016 and had been stored in ethanol. Its whole body was cut in half along the median plane; one half was used for genomic DNA extraction, and the other was used for total RNA extraction. Each half was homogenized independently using a disposable BioMasher II homogenizer (Nippi, Tokyo, Japan). Genomic DNA was extracted using the conventional sodium dodecyl sulfate lysis and phenol–chloroform method. Total RNA was extracted using TRIzol Reagent (Thermo Fisher Scientific) and treated with DNase I (Nippon Gene, Tokyo, Japan) to digest any contaminated genomic DNA. Next, to remove DNase I proteins, the RNA sample was extracted using TRIzol LS Reagent (Thermo Fisher Scientific). Extracted DNA and RNA were confirmed using 1% agarose gel electrophoresis.
The PCRs for COX1 and 16S sequences using the genomic DNA were performed using the same primers and PCR conditions as used by Takahashi et al. . The RT-PCR using total RNA was performed in one tube in a 5.0 μL mixture containing approximately 0.5 μg total RNA, 0.2 μM of each primer, 1× KAPA 2G Fast Multiplex PCR Kits HS (KAPA Biosystems, Wilmington, MA, USA), and 10 U of SuperScriptIII Reverse Transcriptase (Thermo Fisher Scientific). The reaction conditions were 50 °C for 30 min and 94 °C for 2 min, followed by 35 cycles of denaturation at 94 °C for 15 s, annealing at 50 °C for 30 s, and extension at 72 °C for 30 s. The same primers for the COX1 and 16S sequences were used for RT-PCR. The PCR products were confirmed via 1% agarose gel electrophoresis and purified using a High Pure PCR Product Purification Kit (Roche Diagnostics, Mannheim, Germany). DNA sequencing was performed on the PCR products using a BigDye Terminator v1.1 Cycle Sequencing Kit and the ABI PRISM 310 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA). To confirm the sequence, both strands of DNA were sequenced.
Phylogenetic analysis of COX1 and 16S was carried out following Takahashi et al. . Multiple alignments were performed using MUSCLE , implemented in MEGA7. The ML method  was used to construct phylogenetic trees from the COX1 and 16S data using RAxML  with the general time reversible (GTR)  + γ model. Bootstrap probabilities  were computed from 500 replicates. Pair-wise p-distances as percentages were also calculated using MEGA7.
Phylogenetic inference among Branchiopoda species
The nucleotide sequences of 13 protein-coding genes were retrieved, concatenated, and used for the analysis. Multiple alignments for the translated amino acid sequences of each gene were performed using MUSCLE , implemented using MEGA7. Sites containing gaps were removed. The translated amino acid sequences were then returned to the nucleotide sequencing data. The GC content was calculated using MEGA7.
The ML method was used to construct a phylogenetic tree using RAxML-NG 1.0.1  with 1000 bootstrap replicates. The BS approach was used to construct a phylogenetic tree using MrBayes version 3.2  with 10,000,000 generations. Nucleotide substitution models selected by ModelTest-NG 0.1.6 were used for both the ML method and the BS approach. The NJ method  was used to construct a phylogenetic tree with 1000 bootstrap replicates in MEGA7. The “Find Best DNA/Protein Models (ML)” option of MEGA7 was used to select the best-fitting nucleotide substitution models for the NJ tree. To compare the topologies obtained from these three methods, we conducted the likelihood analysis using BASEML with GTR (also known as REV) + γ in PAML version 4.4 . Differences in substitution rates among lineages were tested using BASEML with GTR + γ. We used CODEML in PAML version 4.4 to estimate the ω values, and performed the Kishino and Hasegawa test  and the Shimodaira and Hasegawa test .
Availability of data and materials
DNA sequences determined in this study were deposited in the DDBJ/EMBL/GenBank International Nucleotide Sequence Database (accession numbers: LC633437–LC633442).
12S ribosomal RNA
16S ribosomal RNA
Akaike information criterion
ATP synthase F0 subunit 8
ATP synthase F0 subunit 6
Bayesian information criterion
Cytochrome c oxidase subunit I
Cytochrome c oxidase subunit II
Cytochrome c oxidase subunit III
DNA Data Bank of Japan
DNA integrity number
European Molecular Biology Laboratory
General time reversible
NADH dehydrogenase subunit 1
NADH dehydrogenase subunit 2
NADH dehydrogenase subunit 3
NADH dehydrogenase subunit 4
NADH dehydrogenase subunit 4 L
NADH dehydrogenase subunit 5
NADH dehydrogenase subunit 6
Nuclear mitochondrial DNA segments
Polymerase chain reaction
Reverse transcription polymerase chain reaction
Takahashi N, Kitano T, Hatanaka Y, Nagahata Y, Tshistjakov YA, Hamasaki M, et al. Three new species of the fairy shrimp Eubranchipus Verill, 1870 (Branchiopoda: Anostraca) from northern Japan and far eastern Russia. BMC Zool. 2018;3(1):1–16.
Kikuchi H. Occurrence of a new fairy shrimp, Chirocephalopsis uchidai sp. nov., from Hokkaido, Japan (Chirocephalidae Anostraca). J Fac Sci Hokkaido Univ Ser 6 Zool. 1957;13:59–62.
Moriya H. Notes on a fairy shrimp, Eubranchipus uchidai (Kikuchi) (Anostraca), from Japan. Hydrobiologia. 1985;120:97–101.
Ooyagi A. First records of a fairy shrimp, Eubranchipus uchidai Kikuchi (Anostraca), in Honshu, Japan, with some notes on its ecology. J Nat Hist Aomori. 1996;1:25–9 (in Japanese).
Igarashi S, Mikami H. Water environment of a vernal pool in Ishikari coastal area: habitat characteristics of fairy shrimp in northern part of Japan, Eubranchipus uchidai Kikuchi (CRUSTACEA, ANOSTRACA). Report Hokkaido Institute Environ Sci. 1999;26:31–4 (in Japanese).
Boore JL. Animal mitochondrial genomes. Nucleic Acids Res. 1999;27(8):1767–80.
Yang R, Chen Y. The complete mitochondrial genome of the freshwater fairy shrimp Branchinella kugenumaensis Ishikawa 1894 (Crustacea: Anostraca: Thamnocephalidae). Mitochondrial DNA Part B. 2020;5(1):1048–9.
Cook CE, Yue Q, Akam M. Mitochondrial genomes suggest that hexapods and crustaceans are mutually paraphyletic. Proc Biol Sci. 2005;272(1569):1295–304.
Luchetti A, Forni G, Skaist AM, Wheelan SJ, Mantovani B. Mitochondrial genome diversity and evolution in Branchiopoda (Crustacea). Zool Lett. 2019;5(1):1–13.
Marinier E, Brown DG, McConkey BJ. Pollux: platform independent error correction of single and mixed genomes. BMC Bioinformatics. 2015;16(1):10.
Boisvert S, Laviolette F, Corbeil J. Ray: simultaneous assembly of reads from a mix of high-throughput sequencing technologies. J Comput Biol. 2010;17(11):1519–33.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.
Dierckxsens N, Mardulyn P, Smits G. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 2017;45(4):e18.
Fan YP, Lu B, Yang JS. The complete mitogenome of the fairy shrimp Phallocryptus tserensodnomi (Crustacea: Anostraca: Thamnocephalidae). Mitochondrial DNA Part A. 2016;27(5):3113–4.
Tladi M, Dalu T, Rogers DC, Nyamukondiwa C, Parbhu SP, Teske PR, et al. The complete mitogenome of the fairy shrimp Streptocephalus cafer (Lovén, 1847) (Crustacea: Branchiopoda: Anostraca) from an ephemeral pond in Botswana, southern Africa. Mitochondrial DNA Part B. 2020;5(1):623–5.
Liu XC, Li HW, Jermnak U, Yang JS. The complete mitogenome of the freshwater fairy shrimp Streptocephalus sirindhornae (Crustacea: Anostraca: Streptocephalidae). Mitochondrial DNA Part A. 2016;27(5):3189–91.
Asem A, Li W, Wang P, Eimanifar A, De VS, Van SG. The complete mitochondrial genome of Artemia sinica Cai, 1989 (Crustacea: Anostraca) using next-generation sequencing. Mitochondrial DNA Part B. 2019;4(1):746–7.
Zhang H, Luo Q, Sun J, Liu F, Wu G, Yu J, et al. Mitochondrial genome sequences of Artemia tibetiana and Artemia urmiana: assessing molecular changes for high plateau adaptation. Sci China Life Sci. 2013;56(5):440–52.
Perez ML, Valverde JR, Batuecas B, Amat F, Marco R, Garesse R. Speciation in the Artemia genus: mitochondrial DNA analysis of bisexual and parthenogenetic brine shrimps. J Mol Evol. 1994;38(2):156–68.
Ribeiro MM, Facchin S, Pereira AH, Kalapothakis E, Xu S, Han B, et al. Mitogenome of Daphnia laevis (Cladocera, Daphniidae) from Brazil. Mitochondrial DNA Part B. 2019;4(1):194–6.
Fields PD, Obbard DJ, Mctaggart SJ, Galimov Y, Little TJ, Ebert D. Molecular Phylogenetics and evolution Mitogenome phylogeographic analysis of a planktonic crustacean. Mol Phylogenet Evol. 2018;129:138–48.
Tokishita S, Shibuya H, Kobayashi T, Sakamoto M, Ha J. Diversification of mitochondrial genome of Daphnia galeata (Cladocera, Crustacea): comparison with phylogenetic consideration of the complete sequences of clones isolated from five lakes in Japan. Gene. 2017;611:38–46.
Crease TJ. The complete sequence of the mitochondrial genome of Daphnia pulex (Cladocera: Crustacea). Gene. 1999;233:89–99.
Lee J, Kim D, Choi B, Kato Y, Lee J, Kim D, et al. Complete mitochondrial genome of the freshwater water flea Daphnia magna NIES strain (Cladocera, Daphniidae): rearrangement of two ribosomal RNA genes. Mitochondrial DNA Part B. 2020;5(2):1822–3.
Liu P, Xu S, Huang Q, Dumont HJ, Lin Q, Han B. The mitochondrial genome of Diaphanosoma dubium with comparison with Daphnia magna. Mitochondrial DNA Part B. 2017;2(2):926–7.
Bellec L, Debruyne R, Utge J, Rabet N. The first complete mitochondrial genome of Limnadia lenticularis (Branchiopoda, Spinicaudata), with new insights on its phylogeography and on the taxonomy of the genus. Hydrobiologia. 2019;826(1):145–58.
Tladi M, Dalu T, Rogers DC, Nyamukondiwa C, Emami-khoyi A, Oliver JC, et al. The complete mitogenome of an undescribed clam shrimp of the genus Gondwanalimnadia (Branchiopoda: Spinicaudata), from a temporary wetland in Central District. Botswana Mitochondrial DNA Part B. 2020;5(2):1238–40.
Sun X, Cheng J. International journal of biological macromolecules characterization of the complete mitochondrial genome of Chinese Triops granarius and implications for species delimitation. Int J Biol Macromol. 2019;135:734–44.
Umetsu K, Iwabuchi N, Yuasa I, Saitou N, Clark PF, Boxshall G, et al. Complete mitochondrial DNA sequence of a tadpole shrimp (Triops cancriformis) and analysis of museum samples. Electrophoresis. 2002;23:4080–4.
Gan HM, Tan MH, Lee YP, Austin CM. The complete mitogenome of the Australian tadpole shrimp Triops australiensis (Spencer & Hall, 1895) (Crustacea: Branchiopoda: Notostraca). Mitochondrial DNA Part A. 2016;27(3):2028–9.
Felsenstein J. Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Evol. 1981;17:368–76.
Akaike H. A new look at the statistical model identification. IEEE Trans Automatic Control. 1974;19(6):716–23.
Sugiura N. Further analysts of the data by akaike's information criterion and the finite corrections: further analysts of the data by akaike's. Commun Stat - Theory Methods. 1978;7(1):13–26.
Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogeny. Bioinformatics. 2001;17:754–5.
Schwarz G. Estimating the dimension of a model. Ann Stat. 1978;6:461–4.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987;4:406–25.
Tamura K, Nei M. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993;10(3):512–26.
Kishino H, Hasegawa M. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in hominoidea. J Mol Evol. 1989;29:170–9.
Shimodaira H, Hasegawa M. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol Biol Evol. 1999;16:1114–6.
Stenderup JT, Olesen J, Glenner H. Molecular phylogeny of the Branchiopoda (Crustacea)–multiple approaches suggest a ‘diplostracan’ ancestry of the Notostraca. Mol Phylogenet Evol. 2006;41:182–94.
Olesen J. Phylogeny of Branchiopoda (Crustacea)–character evolution and contribution of uniquely preserved fossils. Arthropod Syst Phylogeny. 2009;67(1):3–39.
von Reumont BM, Jenner RA, Wills MA, Ampio ED, Ebersberger I, Meyer B, et al. Pancrustacean phylogeny in the light of new phylogenomic data: support for Remipedia as the possible sister group of Hexapoda. Mol Biol Evol. 2012;29(3):1031–45.
Schwentner M, Richter S, Rogers DC, Giribet G. Tetraconatan phylogeny with special focus on Malacostraca and Branchiopoda: highlighting the strength of taxon-specific matrices in phylogenomics. Proc R Soc B. 2018;285:20181524.
Uozumi T, Ishiwata K, Grygier MJ, Sanoamuang LO, Su ZH. Three nuclear protein-coding genes corroborate a recent phylogenomic model of the Branchiopoda (Crustacea) and provide estimates of the divergence times of the major branchiopodan taxa. Genes Genet Syst. 2021;96(1):13–24.
Remigio EA, Hebert PDN. Affinities among Anostracan (Crustacea: Branchiopoda) families inferred from phylogenetic analyses of multiple gene sequences. Mol Phylogenet Evol. 2000;17(1):117–28.
Weekers PH, Murugan G, Vanfleteren JR, Belk D, Dumont HJ. Phylogenetic analysis of anostracans (Branchiopoda: Anostraca) inferred from nuclear 18S ribosomal DNA (18S rDNA) sequences. Mol Phylogenet Evol. 2002;25:535–44.
deWaard JR, Sacherova V, Cristescu ME, Remigio EA, Crease TJ, Hebert PD. Probing the relationships of the branchiopod crustaceans. Mol Phylogenet Evol. 2006;39:491–502.
Hebert PDN, Remigio EA, Colbourne JK, Taylor DJ, Wilson CC. Accelerated molecular evolution in halophilic crustaceans. Evolution. 2002;56(5):909–26.
Fontanillas E, Welch JJ, Thomas JA, Bromham L. The influence of body size and net diversification rate on molecular evolution during the radiation of animal phyla. BMC Evol Biol. 2007;7:1–12.
Gillooly JF, Allen AP, West GB, Brown JH. The rate of DNA evolution: effects of body size and temperature on the molecular clock. Proc Natl Acad Sci U S A. 2005;102(1):140–5.
Welch JJ, Bininda-Emonds ORP, Bromham L. Correlates of substitution rate variation in mammalian protein-coding sequences. BMC Evol Biol. 2008;8(1):1–12.
Thomas JA, Welch JJ, Lanfear R, Bromham L. A generation time effect on the rate of molecular evolution in invertebrates. Mol Biol Evol. 2010;27(5):1173–80.
Rogers DC. Branchiopoda (Anostraca, Notostraca, Laevicaudata, Spinicaudata, Cyclestherida). In: Likens GF, editor. Encyclopedia of inland waters. 2. Oxford: Academic Press; 2009. p. 242–9.
Ohta T. Slightly deleterious mutant substitutions in evolution. Nature. 1973;246:96–8.
Chikhi R, Medvedev P. Informed and automated k-mer size selection for genome assembly. Bioinformatics. 2014;30(1):31–7.
Li H, Durbin R. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Kinjo S, Monma N, Misu S, Kitamura N, Imoto J, Yoshitake K, et al. Maser: one-stop platform for NGS big data from analysis to visualization. Database. 2018;2018:bay027.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6.
Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69(2):313–9.
Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16(2):111–20.
Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
Tavaré S. Some probabilistic and statistical problems in the analysis of DNA sequences. Lect Math Life Sci. 1986;17:57–86.
Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39:783–91.
Darriba D, Posada D, Kozlov AM, Stamatakis A, Morel B, Flouri T. ModelTest-NG: a new and scalable tool for the selection of DNA and protein evolutionary models. Mol Biol Evol. 2020;37(1):291–4.
Kozlov AM, Darriba D, Flouri T, Morel B, Stamatakis A. RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics. 2019;35(21):4453–5.
Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MRBAYES 3.2: efficient Bayesian phylogenetic inference and model selection across a large model space. Syst Biol. 2012;61:539–42.
Yang Z. PAML 4: a program package for phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24:1586–91.
Computations were partially performed on the NIG supercomputer at ROIS National Institute of Genetics. We thank Akihiko Utsumi for the technical assistance. We are grateful to the journal editor and reviewers for their comments and suggestions, which have helped us greatly in improving the manuscript.
Ethics approval and consent to participate
N. Takahashi collected the fairy shrimp specimens in the special protection zone of Shiretoko National Park, according to permit no. 1,604,271, which was obtained from the Kushiro Nature Conservation Office, Ministry of the Environment. No approval was required for specimens collected from other locations for this study. Ethical approvals are not required at Yamagata University or Ibaraki University for research conducted on invertebrates such as the fairy shrimp species used in this study.
Consent for publication
The authors have no competing interests to declare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
: Table S1. Summary of sequencing statistics.
: Table S2. Annotation of mitochondrial genes in three Eubranchipus species.
: Table S3. Genetic distance among Eubranchipus species.
: Table S4. Selected nucleoitde substitution models.
: Table S5. Pair-wise nucleotide difference per site in % for COX1 data. Compared sites were 487 bp.
: Table S6. Pair-wise nucleotide difference per site in % for 16S rRNA data. Compared sites were 399 bp.
About this article
Cite this article
Kitano, T., Sato, H., Takahashi, N. et al. Complete mitochondrial genomes of three fairy shrimps from snowmelt pools in Japan. BMC Zool 7, 11 (2022). https://doi.org/10.1186/s40850-022-00111-2