The mitochondrial genome of a sea anemone Bolocera sp. exhibits novel genetic structures potentially involved in adaptation to the deep‐sea environment

Abstract The deep sea is one of the most extensive ecosystems on earth. Organisms living there survive in an extremely harsh environment, and their mitochondrial energy metabolism might be a result of evolution. As one of the most important organelles, mitochondria generate energy through energy metabolism and play an important role in almost all biological activities. In this study, the mitogenome of a deep‐sea sea anemone (Bolocera sp.) was sequenced and characterized. Like other metazoans, it contained 13 energy pathway protein‐coding genes and two ribosomal RNAs. However, it also exhibited some unique features: just two transfer RNA genes, two group I introns, two transposon‐like noncanonical open reading frames (ORFs), and a control region‐like (CR‐like) element. All of the mitochondrial genes were coded by the same strand (the H‐strand). The genetic order and orientation were identical to those of most sequenced actiniarians. Phylogenetic analyses showed that this species was closely related to Bolocera tuediae. Positive selection analysis showed that three residues (31 L and 42 N in ATP6, 570 S in ND5) of Bolocera sp. were positively selected sites. By comparing these features with those of shallow sea anemone species, we deduced that these novel gene features may influence the activity of mitochondrial genes. This study may provide some clues regarding the adaptation of Bolocera sp. to the deep‐sea environment.


| INTRODUCTION
The deep sea is the part of the ocean below the continental shelves, and it is the most extensive ecosystem on earth (Rex, 1981). The organisms living there survive in an extremely harsh environment, tolerating hundreds of bars of pressure, small amounts of oxygen, very little food, constant darkness, and low temperature (Sanders & Hessler, 1969). Because of the sparseness of animal life and the technical difficulties in sampling the deep-sea benthos, our knowledge of deep-sea organisms is based almost entirely on morphological distinctions at the species level (Etter, Rex, Chase, & Quattro, 1999). There is little information about the adaptive molecular mechanisms of the organisms living in the deep-sea environment (Sanders & Hessler, 1969).
As the powerhouse of the cell, mitochondria generate energy by oxidative phosphorylation (OXPHOS) (Luo et al., 2008) and play important roles in energy metabolism and various biosynthetic pathways (Green & Reed, 1998;Newmeyer & Ferguson-Miller, 2003).
Although several complete or partial mitogenomes of sea anemones have been sequenced in recent years, the members of these genera are highly diverse (Emblem et al., 2014), and the information regarding these mitogenomes remains incomplete. For deep-water species from extreme environments in particular, little mitochondrial data have been reported. To evaluate the variation in deep-sea sea anemone gene structures and their adaptations to the deep-sea environment, this study determined the complete mitogenome of a deep-sea sea anemone (Bolocera sp.), identified its mitochondrial gene structures and arrangements, and elucidated their evolutionary characteristics.

| Specimens and DNA extraction
The specimen was collected on December 15, 2014, from a seamount on the Pacific Ocean (137.44˚E/8.52˚N) at a depth of 1106 m, fixed in 99% ethanol, and stored at 4°C. DNA was extracted using a Genomic DNA Kit (Tiangen Co. Beijing, China) according to the manufacturer's instructions. The specimen was not an endangered or protected species, and no specific permits were required for our collection process.

| PCR and sequencing
To identify the subspecies of the specimen, two conserved genes (COI and cytb) were sequenced. The complete mitogenomes of closely related species were downloaded from the NCBI Entrez Database and amplified to search for conserved regions where primers for the complete mitogenome clone were designed. The complete mitogenome was amplified by overlapping PCR. All PCRs were performed in a 50 μl volume, which included 1 μl template DNA (approximately 100 ng), each primer at a concentration of 0.3 μmol/L, 5 μl of 10 × LA Taq buffer (Mg 2+ plus), 5 μl of dNTP Mix (2.5 mmol/L), and 1 U of LA Taq (TaKaRa, Japan). The PCR amplifications used the following procedure: one cycle of denaturation for 5 min at 94°C; 35 cycles of 40 s at 94°C, 40 s at the primer-specific annealing temperature, and 5 min at 72°C; and finally a 10-min extension at 72°C. After purification, the PCR products were directly sequenced in both directions three times with the PCR primers. Sequencing was performed by ThermoFisher Scientific (Guangzhou, China).

| Complete mitogenome analysis
The sequence alignment was conducted using Clustal X. The proteincoding genes and rRNA genes were determined by BLAST and the NCBI Entrez Database and by comparison with the mitogenome sequences of homologous species. The tRNA genes and their secondary structures were predicted by the Web-based tRNAscan-SE 1.21 (Lowe & Eddy, 1997). The skew in nucleotide composition was calculated by AT skew and GC skew and measured according to the following formulae: AT skew = (A − T)/(A + T) and GC skew = (G − C)/ (G + C) (Perna & Kocher, 1995), where A, T, C, and G are the occurrences of the corresponding bases. Codon usage was calculated by the Codon Usage Database (http://www.kazusa.or.jp/codon/). The gene map of the complete mitogenome was depicted by OGDRAW (http://ogdraw.mpimp-golm.mpg.de/).

| Phylogenetic analysis
To illustrate the phylogenetic relationships among sea anemone, the . Geodia neptuni (AY320032) (Demospongiae) was selected as the out-group. The concatenated nucleotide sequences of 13 energy pathway protein-coding genes were aligned using Clustal X with the default settings. The maximum likelihood (ML) method was employed to analyze the phylogenetic tree. The GTR + I + G model was selected as the best nucleotide substitution model by ModelTest 3.7 (Posada & Crandall, 1998). The ML analysis was performed by MEGA 5.1 with 1,000 bootstrap replicates.

| Positive selection analysis
The selective pressure imposed on the mitogenomes of sea anemones was evaluated using CODEML from the PAML package. Two different tree-building methods were used because the CODEML likelihood analysis is sensitive to tree topology. The two-ratio and free-ratio model (M1 model) was used in the mitogenome analysis.
The branch-site model was used to determine whether these genes have undergone positive selection in the foreground lineage. Bayes Empirical Bayes (BEB) analysis was used to calculate the Bayesian posterior probability of the positively selected sites.

| Genome organization
Similar to most metazoan mitogenomes, the mitogenome of Bolocera sp. is a circular molecule. The complete mitochondrial DNA of Bolocera sp. contained 19,463 bp. It shared the highest overall similarity (96.43%) with the B. tuediae mitochondrial DNA sequence (Emblem et al., 2014). In addition to the common individual base composition differences, the alignment of the two complete mitogenomes also showed several insertions or deletions of genetic fragments. Therefore, we speculated that Bolocera sp. was a new subspecies of Bolocera. Similar to those of other anthozoan species, the mitogenome of Bolocera sp. showed a different evolution pattern than most metazoans (Shearer, Van Oppen, Romano, & Wörheide, 2002). It encoded 16 protein-coding genes (13 energy pathway protein-coding genes, a heg gene, and two unknown ORFs), two tRNA genes, and two rRNA genes ( Figure 1, Table 1). All of the genes were coded by the same strand (the H-strand) and transcribed in the same direction. The gene order and orientation were identical to other sequenced actiniarians. In all the mitogenomes studied, no gene overlaps were observed, and the intergenic spacers varied. Several typical anthozoan mitogenome features (Beagley & Wolstenholme, 2013;Beagley et al., 1998;Emblem et al., 2014) were also observed in Bolocera sp., including the presence of two group I introns, two tRNA genes, and large intergenic spacers. In addition, two noncanonical protein-coding genes (ORFC and ORFD) were observed. A genetic fragment similar to the control region (CR) of metazoans was also observed, which was first reported in sea anemones. The complete mitochondrial DNA sequence was deposited in the GenBank database under the accession number KU507297.
F I G U R E 1 Graphical map of complete mitogenome of Bolocera sp. Different genes are represented by different boxes in different colors. tRNAs are displayed according to the one-letter code. Genes encoded by the heavy strand are shown outside the circle, and those encoded by the light strand are shown inside. The direction of the arrows shows the direction of transcription. The inner ring shows the GC content of the mitogenome

| Group I introns and unknown ORFs
In addition to the standard set of 13 energy pathway protein-coding genes, an additional heg gene and two group I introns as well as two unknown ORFs were identified in the Bolocera sp. mitogenome (Figures 1 and 2, Table 1). Group I introns are genetic insertion elements that are extremely rare in metazoans and have only been identified within the mitogenomes of hexacorals and some sponges (Boore, 1999). In Lophelia pertusa (a scleractinian coral), a complex group I intron is inserted in the ND5 gene to host seven essential mitochondrial protein genes and a rRNA gene (Emblem et al., 2011). In this study, two group I introns were detected (ND5 intron and CO I intron). Both of the group I introns contained the same gene components as those in other hexacorallian species except for scleractinians ( Figure 2). A gene structure analysis of hexacorallians showed that the ND5 intron absorbed more genes but that the CO I intron lost the heg gene. The reason for this difference might be that the scleractinian ancestor had a different evolutionary pathway than other hexacorallian suborder species before significant mitogenome rearrangement, which is supported by the opinion of Mónica (Medina, Collins, Takaoka, Kuehl, & Boore, 2006). Consistent with previous studies, no obvious similarity was detected. Flot & Tillier (2007) strongly suggested that the ORFs are expressed proteins, and Emblem even identified the ORFs to be functional elements (Emblem et al., 2014). The ORFs appear to be evolving under some selection relative to other protein-coding genes (Emblem et al., 2014). The patterns (gained, lost, or truncated) exhibited by these noncanonical ORFs are consistent with transposon-like elements (Winckler, Szafranski, & Glöckner, 2005). The ORFs detected in this study showed the same characteristics as the ORFs above, so they may play similar roles. Transposable elements can influence neighboring genes by altering splicing and polyadenylation patterns, or by functioning as enhancers or promoters (Girard & Freeling, 1999;Slotkin & Martienssen, 2007;Waterland & Jirtle, 2003). They have been demonstrated to play essential roles in the host response to stress and in facilitating the adaptation of populations (Blot, 1994;Casacuberta & González, 2013;Chénais, Caruso, Hiard, & Casse, 2012). Considering the common features of the species that carry noncanonical ORFs and the particular characteristics of the environment in which they live, the ORFs identified in this study may play important roles in adjusting mitochondrial energy metabolism.
T A B L E 1 Gene structure of the mitogenome of Bolocera sp.

| Genome composition and skewness
In metazoan mitogenomes, the frequency of each nucleotide utilized varies among different taxa. The mitogenome of metazoans usually has a strand-specific bias in nucleotide composition (Alexandre, Nelly, & Jean, 2005). The A + T content of the mitogenome is extremely high in insects and nematodes and lower in vertebrates (Saccone, Giorgi, Gissi, Pesole, & Reyes, 1999). In the anthozoan mitogenome, the A + T content ranges from 54.9% to 68.1% (Brugler & France, 2007;Foox et al., 2016). The nucleotide composition of the H-strand in Bolocera sp. was biased toward A and T, and the overall A + T content was 60.51%. Similar results have also been observed in other actiniarian mitogenomes (Foox et al., 2016;Zhang & Zhu, 2016). As is commonly found in metazoan mitogenomes, the A + T content varied in different regions ( Table 2). The CR-like element presented the highest A + T content (66.16%), and the lowest content was found in the 2 tRNA genes (49.65%).
The AT skew was negative (−0.127), whereas the GC skew was positive (0.108) in Bolocera sp. In all sequenced actiniarian mitogenomes, the trends were the same, and the AT-skew and GC-skew values were similar (Table 2). This result indicated that actiniarian mitogenomes favor Ts and Gs. Similar nucleotide skew patterns have also been observed in the mitogenomes of other hexacorallian subclasses (Brugler & France, 2007;van Oppen et al., 2000).

| Protein-coding genes
In this study, 16 protein-coding genes (13 energy pathway proteincoding genes, a heg gene, and two unknown ORFs) were identified ( Table 1). All of these genes were coded by the same strand (the F I G U R E 2 Linearized schemes of mitochondrial gene arrangements in anthozoans. (a) Linearized mitochondrial gene arrangements in actiniarians, Stichopathes lutkeni (Antipatharia) as outgroups. (b) Linearized mitochondrial gene arrangements in different suborders of anthozoa. Lengths of the genes correspond to relative lengths of the genomes in a. tRNAs are displayed according to the one-letter code. Species names and NCBI accession numbers are given under each of the linearized schemes H-strand) and transcribed in the same direction. In metazoans, most of the mitochondrial protein-coding genes start with an ATN codon (Liao et al., 2010;Ma et al., 2014;Wang, Chao, Fang, & Yu, 2016).
In the present study, except for the inserted heg gene, which is considered to be a selfish genetic element (Edgell, 2009), all of the genes were initiated by typical ATG codons. This pattern of initiation codon usage has also been observed in other basal metazoans Boore & Brown, 1995;Pont-Kingdon et al., 1998).
However, in this study, all of the protein-coding genes were terminated by complete TAA (11) and TAG (5) termination codons. This situation has also been observed in some sequenced cnidaria (Flot & Tillier, 2007;Kayal & Lavrov, 2008).

| Transfer RNA genes
In most cnidarian mitogenomes, only two tRNAs (tRNA Trp and tRNA Met ) were detected (except for octocorallians, in which only one tRNA Met was detected) Beaton, Roge, & Cavalier-Smith, 1998;Kayal & Lavrov, 2008). A study by Beagley and Wolstenholme (2013) showed that nuclear DNA-encoded functional tRNAs were detected in mitochondria, and the missing tRNAs are believed to be encoded in the nucleus and later imported into the mitochondrion (Schneider & Maréchal-Drouard, 2000). In this study, the mitogenome of Bolocera sp. contained two tRNAs (tRNA Trp and tRNA Met ) (Figure 3). Both tRNAs could fold into a clover-leaf secondary structure, and the anticodon usage was identical to most of the observed sea anemone species. One mismatched base pair (C-A) was detected in tRNA Met . Interestingly, the unmatched base pair occurred on the amino acid acceptor arm. Such stem mismatches seem to be a common phenomenon for mitochondrial tRNAs in many species (Jiang et al., 2013;Liao et al., 2010;Miller et al., 2005;Wang et al., 2016) and are probably corrected by a post RNA-editing mechanism (Lavrov, Brown, & Boore, 2000). T A B L E 2 Genomic characteristics of the mitogenome of Bolocera sp.

Length (bp)
A + T The heg gene and two ORFs do not counted in the 13 energy pathway protein-coding genes. As it was not constrained in the same way as the protein-coding genes, the mitochondrial CR is usually considered to be the most variable portion of the mitogenome (Marshall & Baker, 1997), showing the highest variation in the whole mitogenome (Aquadro & Greenberg, 1983). The structure of the CR has been intensively investigated in vertebrates, but not in invertebrates. Generally, in vertebrates, the mitochondrial CR is divided into three domains, including the ex- However, the nucleotide sequence, length, and number of each motif all vary considerably among vertebrate classes and even within a class (Rand, 1993;Ruokonen & Kvist, 2002;Shaffer & McKnight, 1996). In invertebrates, especially in anthozoans, similar CR structures are not clearly defined, but some similar features have been observed. In A. tenuis, the candidate mitochondrial CR contains repetitive elements and has the potential to form the typical secondary structures of vertebrate D-loops (van Oppen et al., 2000). In Pocillopora, the candidate mitochondrial CR exhibits three characteristics: large size, variability in nucleotide composition, and tandemly arranged repeated sequences (Flot & Tillier, 2007). With the exception of the repeated sequences, the rest of these characteristics were all observed in the 789-bp noncoding region identified here. In addition, the special CR "G(A)nT" motif that is present in S. gregaria and C. parallelus (Zhang, Szymura, & Hewitt, 1995) was also observed in the 789-bp noncoding region ( Figure 4). In all organisms except primates, A + T > G + C in the CR domains (Sbisà et al., 1997), which was also observed in the 789-bp noncoding region. Replication of mitogenome has been shown to initiate near hairpin structures (Clayton, 2000). In Drosophila, the origin of replication is located near a conserved stem-loop structure (Saito, Tamura, & Aotsuka, 2005). In this study, the secondary structure of the 789-bp noncoding fragment presented several similar stem-loop structures ( Figure 5), which is a characteristic feature of O L in vertebrates (Clayton, 1991). In short, while no repetitive elements were F I G U R E 3 Inferred secondary structure of tRNAs in mitogenome of Bolocera sp.

F I G U R E 4
The CR-like sequences of Bolocera sp. The T + A-rich regions were underlined, and the "G(A)nT" motifs were marked with box observed, the 789-bp noncoding region exhibited typical CR structures observed in other invertebrates. Considering the primitiveness of sea anemones in evolution, there must be some primitive features maintained. Therefore, we concluded that the 789-bp noncoding sequence was a candidate CR (that we defined as "CR-like") that has not been reported in other actiniarians. This would be a unique and/ or primitive mitogenomic feature for Bolocera sp., which is in agreement with the CR differences observed within invertebrates. As no other obvious adjusting elements were detected in the Bolocera sp.
mitogenome, this special CR-like element may play an important role in regulating the transcription and replication of the mitogenome in the extreme environment of the deep sea.

| Phylogenetic analyses
To investigate the relationships among anthozoan, we performed phylogenetic analysis based on the nucleotide datasets of the 13 mitochondrial energy pathway protein-coding genes ( Figure 6). All of the topologies showed a high support value. The phylogenetic tree of Anthozoa in this study showed that Bolocera sp. was clustered in the Actiniaria clade and had the closest relationship with Bolocera tuediae. The Actiniaria represented a monophyletic group and together with the sister groups Zoanthidea, Antipatharia, and Scleractinia clustered with the Hexacorallia group, which supported the classification of anthozoan (Daly, Fautin, & Cappola, 2003). Based on seven mitochondrial protein-coding genes, Brugler and France (2007)  zoanthid was at the base of the hexacoral clade, and the antipatharian clade had high support as a sister-taxon to the scleractinian clade.
However, in this study, based on the 13 protein-coding genes, the antipatharians were a sister-group to zoanthids, with a closer relationship to actiniarians, while scleractinians were located at the base of the hexacoral clade, which was identical to results obtained based on 18 S rRNA (Berntson, France, & Mullineaux, 1999). These findings also corroborated earlier studies of sea anemone phylogenetic relationships based on short mitochondrial and nucleotide sequences (Daly, Chaudhuri, Gusmão, & Rodríguez, 2008;Emblem et al., 2014).

| Positive selection analysis
The selective pressures imposed on the mitogenomes of sea anemones were evaluated using CODEML from the PAML package (Table 3).
Two different tree-building methods were used because the CODEML likelihood analysis is sensitive to tree topology. There were no significantly different ω ratios between branches of genus Bolocera and other species when we set the genus Bolocera as a foreground branch using the two-ratio model (p > .05). In the analyses of individual genes, we found that three positive selection sites (31 L and 42 N in ATP6,

S in ND5)
showed BEB values >0.95 using branch-site models.
The NADH dehydrogenase complex, which likely functions as a proton pump (da Fonseca et al., 2008) and influences metabolic performance (Hassanin et al., 2009), has been considered important in the adaptive evolution of the mammalian mitogenome (da Fonseca et al., 2008;Xu et al., 2007). ND2 and ND6 were found to be under positive selection pressure in a mitogenome analysis of Chinese snubnosed monkeys, which is suggestive of adaptive changes related to high altitude and cold weather stress (Yu et al., 2011). ND6 was found to be under positive selection pressure in Tibetan horses living at high altitude (Ning et al., 2010). ATP synthase is directly associated with the produce of ATP (Mishmar et al., 2003;Weiss, Friedrich, Hofhaus, & Preis, 1991;Zhou, Shen, Irwin, Shen, & Zhang, 2014). It has been suggested that variation in ATPase proteins could result in significant variation in mitochondrial adaptation to different environments (Mishmar et al., 2003;Wallace, 2007). Members of the Caprini tribe that live in high-altitude mountain regions show higher levels of adaptive evolution in the ATP synthase complex (Hassanin et al., 2009). The positive selection of ATP6 and ND5 observed in the present study could help us to better understand the adaptation of organisms to the deep-sea environment.

| CONCLUSION
This study characterized the complete mitogenome of a deep-sea benthos species, Bolocera sp., which was speculated to be a new species of Bolocera. The study provided the following conclusions about deep-sea organisms: (1) These organisms have monophyletic genome characteristics similar to those of shallow sea organisms: The basic gene content, order, and orientation of the species were identical to most of those reported homologous species, and phylogenetic analyses indicated that Bolocera sp. is closely related to Bolocera tuediae and belongs to the Actiniidae family.