The first complete mitochondrial genome of the Mariana Trench Freyastera benthophila (Asteroidea: Brisingida: Brisingidae) allows insights into the deep‐sea adaptive evolution of Brisingida

2.1Sample collection and DNA extraction

The specimen was collected at Mariana Trench, in June 2016 (10°51.0971′N, 141°57.2705′E, at 5,463 m depth). The collection was accomplished by deep‐sea human occupied vehicle (HOV) “Jiao Long” during an expedition. The specimen was preserved in 95% ethanol. Total genomic DNA was extracted from ethanol‐fixed tissue with tissue DNA kit (Omega Bio‐Tek) and stored at −20°C.

2.2PCR amplification and DNA sequencing

The universal metazoan primers for mtDNA were used in PCR. Three fragments cox3, cob, and 16S were successfully amplified with the primers cox3F + cox3R (Boore, Macey, & Medina, 2005), cobF424 + cobR876 (Boore & Brown, 2000), and 16SarL + 16SbrH (Boore et al., 2005). In addition, partial sequences of cox1, cox2, nad4L, atp6, nad4, and nad5 were amplified with the degenerate primers from conserved regions of other starfish in GenBank. The remaining gaps were amplified with the species‐specific primers designed according to the obtained sequences. Finally, the whole mitogenome was amplified based on eight pairs of primers (Table 1).

PCRs were performed with a gradient machine (Applied Biosystems Inc.). The cycling was set up with an initial denature step at 94°C for 5 min, followed by 35 cycles (94°C for 30 s, 45–55°C for 1 min, 72°C 1–3 min), and a final extension was executed at 72°C for 10 min. LA‐PCR was carried out in a 20 μl reaction volume containing 12.6 μl ddH₂O, 2 μl 10 × LA‐PCR buffer (Mg²⁺ plus, Takara), 3.2 μl dNTP mix (2.5 mM each), 0.5 μl each primer (10 μM), 0.2 μl LA Taq DNA polymerase (5 U/μl, Takara), and 1 μl DNA template (50 ng/μl). PCR products were electrophoresed on a 1.0% agarose gel and purified with gel extraction kit (Omega Bio‐Tek) and sequenced with ABI 3730x1 DNA analyzer (Applied Biosystems Inc.).

2.3Gene annotation and Sequence analysis

Raw sequencing reads were first processed using Phred with the quality score 20 and assembled in Phrap with default parameters (Ewing & Green, 1998; Ewing, Hillier, Wendl, & Green, 1998). Then, all assemblies and sequence quality were verified manually in Consed (Gordon, Abajian, & Green, 1998). DOGMA (Wyman, Jansen, & Boore, 2004), ORFfinder (http://www.ncbi.nlm.nih.gov/projects/gorf/orfig.cgi), and BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) were used to identify protein‐encoding genes and rRNA genes. The tRNA genes were identified by tRNAscan‐SE 1.21 (Lowe & Eddy, 1997) and ARWEN 1.2.3.c (Laslett & Canbäck, 2008). Secondary structures for tRNAs were drawn using MITOS Web server (Bernt et al., 2013). Codon usage analysis was estimated with CodonW 1.4.4 (Peden, 1999). The mitochondrial gene map was drawn with GenomeVx (Conant & Wolfe, 2008).

2.4Phylogenetic analysis

Twenty echinoderm mt genomes including the one obtained in this study were used for phylogenetic analysis. All complete mtDNA sequences vailable in GenBank are listed in Table 2. Crinoidea is generally considered as the earliest diverged group of echinoderms (Scouras & Smith, 2006). In this study, Antedon mediterranea (Crinoidea) was rooted as the out‐group. The amino acid sequence from each of 13 protein‐coding genes was aligned separately using Clustal ×2.0 (Larkin et al., 2007), and then, the relatively poor homologous sequence was eliminated. The aligned amino acid sequences were concatenated into a single dataset. The phylogenetic reconstruction approach was performed using neighbor joining (NJ) and maximum likelihood (ML) with MEGA 5.0 (Tamura et al., 2011). The assessment of node reliability was performed using 1,000 bootstrap replicates.

2.5Positive selection analysis

To evaluate the variation in selective pressure between deep‐sea F. benthophila and other eight shallow sea starfish, we used a codon‐based likelihood approach implemented in the CODEML program of the pamlX package (Xu & Yang, 2013; Yang, 2007). All models correct the transition/transversion rate and codon usage biases (F3 × 4). The branch model tests were used to analyze the difference of selective pressure between the deep‐sea and shallow sea starfish. The “one‐ratio” model (model 0), “free‐ratio” model (model 1), and “two‐ratio” model were used in the combined dataset of 13 protein‐coding genes. Considering that positive selection may occur in some amino acids during the evolution of a protein, we used two branch site models (A and A null). Bayes empirical Bayes (Yang, Wong, & Nielsen, 2005) analysis was used to calculate the posterior probabilities of a specific codon site.

3.1General features

The mitogenome of the F. benthophila is a 16,175‐bp circular molecule (Figure 1) with a nucleotide composition of 34.70% A, 21.13% C, 10.65% G, and 33.53% T bases. The genome has an overall A + T content of 68.23%, which appears to be high for Asteroidea. Among the eight species in Asteroidea, the lowest A + T content is 56.34% in Acanthaster planci (Table 3). Freyastera benthophila has the smallest complete mitogenome found in Asteroidea thus far. The size of Asteroidea mitogenomes ranged from 16,524 bp in Luidia quinaria to 16,175 bp in F. benthophila (Table 3). The synteny and identity level between F. benthophila and each of the other seven starfish mitogenomes is shown in Figure 2. The lack of similarity between F. benthophila and L. quinaria is the most obvious feature in the plot.

Figure 1

Mitochondrial gene map of Freyastera benthophila. All of 37 genes are encoded on the both strands. Genes for proteins andrRNAs are shown with standard abbreviation. Genes fortRNAs are designated by a single letter for the corresponding amino acid with two leucinetRNAs and two serinetRNAs differentiated by numerals

Figure 2

The synteny and identity level of Freyastera benthophila mitogenome against each of the other seven starfish mitogenomes. Ideograms and ribbons represent the similarity pairwise blastn searches. In F. benthophila ideogram, the 13 coding genes are colored in blue, control regions are colored in black, andrRNAs are colored in green. The figure was produced using Circoletto (Darzentas, 2010).FB(F. benthophila),AB(Acanthaster brevispinus),AP1 (Acanthaster planci),AJ(Aphelasterias japonica),AA(Asterias amurensis),AP2 (Astropecten polyacanthus),LQ(Luidia quinaria),PP(Patiria pectinifera)

The genome encodes 37 genes including 13 protein‐coding genes (PCGs), two rRNA genes, and 22 tRNA genes (duplication of two tRNAs: trnL and trnS) on both strands. Fifteen of the genes are encoded on the negative strand, while the other 22 are encoded on the positive strand. A total of 22 noncoding regions were found, with the largest continuous region (284 bp, A + T = 67.25%) located between trnT and 16S. Due to its AT richness, we predict that this part is mitochondrial control region. Furthermore, we found four overlaps: trnC/trnV, trnA/trnL₁, atp8/atp6, and cox3/trnS₂. Table 4 presents a summary of the organization of F. benthophila mitogenome. The complete mitochondrial DNA sequence has been deposited in GenBank (Accession Number: MG563681).

3.2Protein‐coding genes

With regard to PCGs, nine (cox1‐cox3, nad3‐nad5, nad4L, cob, atp6, and atp8) are encoded by the positive strand, and the remaining three (nad1, nad2, and nad6) are encoded by the negative strand. These features have been observed in all Asteroidea mitogenomes published so far. Thirteen PCGs initiate with the standard start codon ATG. Most of PCGs terminate with the stop codon TAA (9 of 13), and three genes terminate with the stop codon TAG. Incomplete termination codon T is used by cob. However, mitogenomes often use a variety of nonstandard initiation codons (Wolstenholme, 1992). Nonstandard initiation codon GTG and incomplete termination codon TA are also used in other starfish (Table 3). The lengths of PCGs are 11,506 bp, and the A + T content is 67.15% higher than that of other Asteroidea species (Table 3).

The codon usage of F. benthophila is shown in Figure 3. Among PCGs, leucine (15.85%) and cysteine (0.99%) are the most and the least frequently used amino acids, respectively. Codons, UUA (leucine 6.67%) and ACG (threonine 0.08%), are the most and the least frequently used, respectively. We predict that the richness of A and T occurrence frequency of the mitogenome caused the corresponding amino acid bias to some extent. It is obvious that the A + T content of the third codon position (74.10%) is higher than that of the first (63.43%) and second positions (63.67%).

Figure 3

Codon usage in Freyastera benthophila. All codons for amino acids have been classified. Each amino acid is designated by a single letter for the corresponding codon. x‐axis and y‐axis represent the used times of each codon

3.3Ribosomal RNA and transfer RNA genes

Boundaries of both the small and the large ribosomal genes were determined by BLAST and DOGMA. The 16S and 12S genes of F. benthophila are 1,602 bp (A + T = 72.28%) and 891 bp (A + T = 65.66%) in length, respectively. These lengths are typical for Asteroidea, whereas the AT contents are higher than those of other starfish (Table 3).

We analyzed the entire mitogenome sequence of F. benthophila and successfully identified 22 tRNA genes based on their potential secondary structures using the tRNAscan‐SE, ARWEN, and MITOS Web server (Table 4, Supporting Information Figure S1). The length of these tRNA genes ranged from 68 bp (trnS₁ and trnE) to 75 bp (trnK). Twenty‐one of these genes displayed a common cloverleaf secondary structure, and the remaining one lacked a DHU arm from trnS₁. The D‐stem absence has been found in many other starfish, such as Acanthaster brevispinus, Acanthaster planci, Aphelasterias japonica, Asterias amurensis, L. quinaria, and Patiria pectinifera (Asakawa, Himeno, Miura, & Watanabe, 1995; Matsubara et al., 2005; Tang et al., 2014; Yasuda et al., 2006).

3.4Gene arrangement

Mitochondrial gene arrangement has been demonstrated to be an effective means to solve the deep phylogenetic studies (Boore, 1999; Boore & Brown, 1998). In recent years, some research on mt gene arrangement of echinoderms has been reported (Arndt & Smith, 1998; Perseke et al., 2008, 2010; Scouras et al., 2004).

In this study, mitochondrial gene order of echinoderm was compared among species within classes Asteroidea, Echinoidea, Holothuroidea, Ophiuroidea, and Crinoidea (Figure 4). We expected that the mt gene order of starfish may reveal some phylogenetically information. However, the gene component and gene order of eight species of Asteroidea are completely identical to each other. This phenomenon also happened in the class Echinoidea. We obtained 27 complete mt genomes of Echinoidea from NCBI genebank, and the gene component and gene order of 27 species of Echinoidea are also completely identical to each other. Because Strongylocentrotus purpuratus has been considered as a model for developmental and systems biology, we took S. purpuratus as a representative for Echinoidea in Figure 4 (Sodergren et al., 2006). However, mitochondrial gene order has undergone significant changes in the classes of Holothuroidea, Ophiuroidea, and Crinoidea. Scouras et al. (2004) suggested that it is difficult to resolve the echinoderm phylogeny using the mitochondrial gene rearrangement.

Figure 4

Comparison of mitochondrial gene arrangement in Echinodermata. The bars show identical gene blocks. The noncoding regions are not presented, and gene segments are not drawn to scale

It is interesting that mt gene order of the species in the classes of Asteroidea and Echinoidea is completely identical to each other. If the tRNA is not considered, gene order of PCGs in species within the class Holothuroidea is also the same. This raises the questions: As these species are distributed throughout the world's oceans, why had the mt gene order not been changed and how do they evolve over time. More studies of mt genome species are needed to further investigate whether this pattern is common among starfish, sea urchins, and sea cucumbers.

3.5Phylogenetic analysis

The gene order and transcriptional orientation of the eight Asteroidea species are completely identical to each other, so the mt genome structures would not provide the phylogenetic information. Thus, we performed the phylogenetic analysis using all amino acids of mt protein‐coding genes (Figure 5). Almost all the phylogenetic relationships are supported with high values (NJ/ML bootstraps 99–100). Acanthaster brevispinus is first clustered with A. planci and then united with P. pectinifera; meanwhile, Astropecten polyacanthus and L. quinaria formed a clade. And these five starfish formed the Valvatacea clade. Then, A. japonica is first clustered with Pisaster ochraceus and then united with A. amurensis. Finally, F. benthophila with these three species formed a Forcipulatacea clade. Blake (1987) recognized that Brisingida and Forcipulatida are the two orders within the Forcipulatacea and suggested that they were the most primitive asteroids (Blake, 1987, 1988). Mah and Foltz (2011) described that the largest clade within the Forcipulatacea is formed by the Brisingida and Asteriidae, which forms a clade of deep‐sea and Southern Hemisphere taxa. In the present study, the results supported the deep‐sea Brisingida as a sister taxon to the traditional members of the Asteriidae, and the branch support values are higher than those in previous studies (Glover et al., 2016; Mah & Foltz, 2011). However, the number of Brisingida species with complete mitogenome is still limited, and more mitogenomes and analysis are necessary to determine the phylogenetic relationship among members of Brisingida.

Figure 5

Phylogenetic trees based on the concatenated amino acids of 13 protein‐coding genes. The branch length is determined withNJanalysis. Antedon mediterranea was used as out‐group.NJ(left number) andML(right number) bootstrap values are given for each branch. The red dot highlights the species sequenced in this study

3.6Positive selection analysis

We examined the potential positive selection in Brisingida lineage because of the colonization of deep‐sea environments which may affect the function of mitochondrial genes. The results of selective pressure analyses are shown in Table 5. When the ω ratios for the 13 concatenated mitochondrial protein‐coding genes were tested between the deep‐sea F. benthophila and other eight shallow sea starfish, we failed to find a significant difference in their ω ratios, which may be due to the large bias of sample sizes (p >0.05) (Table 5). In addition, in the analyses of individual genes, we found five residues with high posterior probabilities in the atp8 (8 N, 16 I), nad2 (47 D, 196 V), and nad5 (599 N), respectively (Table 5). Similar results have been observed in deep‐sea animals, and the authors concluded that it may be related to the adaptation to environment (Sun et al., 2018; Zhang et al., 2017). Under the deep‐sea extreme environment, survival may require a modified and adapted energy metabolism (Sun et al., 2018).

Because ATP synthase directly produces ATP, variation in ATPase protein sequence should influence ATP production (Mishmar et al., 2003; Wallace, 2007). Amino acid variations have been widely reported in the ATPase proteins (da Fonseca et al., 2008; Mishmar et al., 2003; Zhang et al., 2017; Zhou et al., 2014). Nad2, nad4, and nad5 are suggested to act as proton‐pumping devices (Brandt, 2006; da Fonseca et al., 2008); thus, mutations in these proteins should influence metabolic efficiency (da Fonseca et al., 2008; Hassanin, Ropiquet, Couloux, & Cruaud, 2009; Zhang et al., 2017). Therefore, we predict that mitochondrial protein‐coding genes, specifically atp8, nad2, and nad5, may play an important role in F. benthophila's adaptation to deep‐sea environment.