Genetic Basis of Hydrothermal Vent Adaptation in Bythograeidae Crabs: Insights from Adaptive Evolution of Mitochondrial Protein Coding Genes

Deep-sea hydrothermal vents are very unusual environments that have attracted the attention of scientists since their discovery on the Galapagos Ridge in 1976 (Corliss, 1977; Jannasch and Mottl, 1985; Bettencourt et al., 2010; Sun et al., 2017). Hydrothermal vents are considered short-lived habitats because their temperature is much warmer than that of the surrounding deep-sea (Childress and Fisher, 1992). Hydrothermal vents are characterized by darkness, lack of photosynthetically derived nutrients, low oxygen concentrations, high hydrostatic pressure, and high concentrations of heavy metals and other toxic substances (Hourdez and Lallier, 2007; Martin et al., 2008; Sun et al., 2017). Hydrothermal vent communities illustrate an apparent paradox of displaying a very high biomass despite very challenging environmental conditions when compared to the surrounding deep-sea fauna (Hourdez and Lallier, 2007). Hydrothermal vent ecosystems support dense populations of macrobenthos which are fueled by simple reduced molecules, such as methane and hydrogen sulfide, with the help of chemoautotrophic endosymbionts (Hourdez and Lallier, 2007; Sun et al., 2017). Thus, the organisms living in hydrothermal vents have developed adaptive mechanisms to survive in this hostile environment (Min et al., 2017).

Brachyuran crabs (Bythograeidae) are biomass dominant and many are known to inhabit deep-sea hydrothermal vents (Min et al., 2017). These crabs are considered the most common and abundant species of predator at the top of food chain in the vent ecosystem (Kim et al., 2013). Survival in such extreme conditions requires unique anatomical and physiological adaptations (Bettencourt et al., 2010). These crabs are characterized with specific physiological and morphological traits including reduced eyestalks with vestigial cornea and whole white body (Min et al., 2017). However, studies on the genetic basis for adaptation by vent crabs to the harsh hydrothermal vent environment remain insufficient.

Mitochondria consume the greatest amount (85–90%) of oxygen in cells for oxidative phosphorylation (OXPHOS), which is the primary metabolic pathway for ATP production (Solaini et al., 2010). Five protein complexes in mitochondria are prerequisite for OXPHOS, including nicotinamide adenine dinucleotide (NADH) dehydrogenase (complex I), succinate dehydrogenase (complex II), cytochrome bc1 complex (complex III), cytochrome c oxidase (complex IV), and ATP synthase (complex V) (Eubel et al., 2004; Carroll et al., 2009). The enzymes involved in these complexes are encoded by both nuclear and mitochondrial genes (Wang et al., 2017). Among them, 13 enzymes are encoded by mitochondrial DNA (mtDNA), including seven NADPH reductase subunits (ND1, ND2, ND3, ND4, ND4L, ND5, and ND6), three cytochrome oxidase subunits (COX1, COX2, and COX3), two ATP synthetic subunits (ATP6 and ATP8), and cytochrome b (Cytb) (Kulawiec et al., 2008). Hypoxia is a major threat to the OXPHOS pathway, and if oxygen levels are very low, insufficient ATP availability will result in cell death (Santore et al., 2002). Therefore, mitochondrial protein coding genes are typically chosen to explore the molecular basis of energy metabolism and tolerance to hypoxia (Ning et al., 2010; Shen et al., 2010; Yang et al., 2014; Shi et al., 2018; Wang et al., 2017; Guo et al., 2018).

Thus, the objectives of this study were to increase our understanding of the mechanisms of hypoxia adaptability in vent crabs and to confirm if there was a correlation between mitochondrial protein coding genes (PCGs) and adaptation to the extreme hydrothermal vent environment by comparing the sequences of these genes in Bythograeidae with those in other crabs. Our results demonstrate that a series of positive selection events occurred specifically on the mitochondrial PCGs and on the ancestral and interior lineages of Bythograeidae, suggesting that adaptive evolution of these genes was necessary for Bythograeidae crabs to adapt to the hydrothermal vent environment.

The analysis was comprised of 48 Brachyura species from 19 families (i.e., Homolidae, Leucosiidae, Portunidae, Mithracidae, Majidae, Xanthidae, Menippidae, Bythograeidae, Macrophthalmidae, Varunidae, Grapsidae, Ocypodidae, Dotillidae, Xenograpsidae, Sesarmidae, Mictyridae, Parathelphusidae, Potamidae, and Raninidae), and one Diogenidae species (Clibanarius infraspinatus) as the outgroup (Supplementary Table I). We downloaded the full-length coding sequence (CDS) of 13 mitochondrial PCGs from the National Center for Biotechnology Information (NCBI: https://www.ncbi.nlm.nih.gov/). We used two alignment methods (i.e. CLUSTAL and MUSCLE) as implemented in MEGA 7.0 (Kumar et al., 2016) to align the nucleotide sequences of each mitochondrial PCG and verified them by visual inspection.

A phylogenetic tree of 48 Brachyura species and one outgroup C. infraspinatus (Supplementary Table I) was reconstructed from concatenated 13 mitochondrial PCGs using the Bayesian inference (BI) and maximum likelihood (ML) algorithms. The BI and ML trees were reconstructed using MrBayes 3.2.6 (Huelsenbeck and Ronquist, 2001; Ronquist et al., 2012) and RaxML (Stamatakis, 2014), respectively. jModeltest (Darriba et al., 2012) was used to select the best model, and the SYM+I+G model was optimized to analyze nucleotide alignment. MtArt+I+G was the appropriate model for the amino acid sequence dataset according to ProtTest 3.4 (Abascal et al., 2005; Darriba et al., 2011) based on Akaike’s information criterion. Two simultaneous runs of 10,000,000 generations were conducted for the matrix in the BI analysis. We made two simultaneous runs, sampling trees every 1,000 generations, with three heated and one cold chain to encourage swapping among the Markov-chain Monte Carlo chains. Convergence of the sampled parameters and potential autocorrelation (effective sampling size/ESS for all parameters > 200) were investigated in Tracer 1.6 (http://tree.bio.ed.ac.uk/software/tracer/). Additionally, the average standard deviation of split frequencies between both runs was checked (<0.01). The Bayesian posterior probabilities were obtained from the 50% majority rule consensus of the post-burn-in trees sampled at stationarity, after removing the first 25% of trees as a “burn-in” stage. The resulting phylogenetic trees were visualized in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

The non-synonymous (dN) to synonymous (dS) rate ratio ω is a measure of selective pressure, with values of ω = 1, > 1, and < 1 indicating neutral selection, positive selection, and purifying selection, respectively (Ohta, 1992). The PAML package (Yang, 2007) was used to determine whether adaptive evolution might have occurred in the mitochondrial PCGs of vent crabs. Alignment and consensus trees were used for the posterior molecular evolutionary analysis (Fig. 1).

We used a pair of site models for positive selection at individual codons in vent crabs for each gene, i.e., M8 and M8a (Yang and Nielsen, 2000). The M8a model only allows codons to evolve neutrally or under purifying selection (ω < 1), whereas the M8 model includes a class of sites with ω > 1 (Swanson et al., 2003). We used the free-ratio model (Yang, 1998; Yang and Nielsen, 1998) and branch-site model implemented in CODEML to evaluate whether positive selection was restricted to specific Bythograeidae lineages (Yang, 2007). The positive selection analysis was restricted to the branch of interest (branch leading to the most recently reconstructed ancestor of Bythograeidae, marked a in Figure 1; or the internal branches of Bythograeidae, marked b–i in Fig. 1). The improved branch-site model A (test 2) was performed for every gene in each foreground lineage, which facilitated analysis of the datasets, including all Brachyura species (branch a in Fig. 1) and Bythograeidae only (branches b–i in Fig. 1). The likelihood ratio test (LRT) statistic (2ΔL) approximates a chi-square distribution and was used to compare nested likelihood models. The Bayes empirical Bayes (BEB) approach (Yang et al., 2005) was used to identify amino acids under selection for PAML. We considered candidates sites with posterior probability ≥ 80%.

Positive selected sites were further employed in a series of ML methods implemented in the Datamonkey web server (http://www.datamonkey.org), which has the advantage of improving the dN/dS ratio estimate by incorporating variation in the rate of synonymous substitution (Pond and Frost, 2005; Poon et al., 2009). The fixed-effect likelihood (FEL) calculates site-by-site dN/dS without assuming a prior distribution, and random-effect likelihood (REL) assumes a prior distribution across sites (Poon et al., 2009). Sites with p-values < 0.1 for FEL, and Bayes factor > 50 for REL were considered candidates under positive selection. We then detected selective pressure using the TreeSAAP program (Woolley et al., 2003), which further supports PAML and Datamonkey at the protein physicochemical level. TreeSAAP detected selection based on 31 physicochemical amino acid properties, which were all magnitude category 6–8 changes, and p-values ≤ 0.01 were used as an index for the degree of radical amino acid substitution and positive selection. Finally, to provide further insight into the underlying effects of these positively selected sites, we mapped them onto the protein secondary and three-dimensional structures.

In this study, the 13 PCGs from the mitochondrial genomes of 48 Brachyura species and one outgroup (C. infraspinatus) were examined. We constructed phylogenetic trees using two different methods (BI and ML) with the concatenated 13 mitochondrial PCGs. The relationships of the gene trees obtained from BI and ML were similar to those previously estimated with morphological and molecular data (Dixon et al., 2003; Shen et al., 2013). We used the consensus tree (Fig. 1) that included all species employed in the present study as the working topology in subsequent analyses.

We first used the M0 (one-ratio) model that only allows a single ω ratio for all Brachyura branches. The ω values obtained for all 13 PCGs ranged from 0.00895 to 0.11106 and were significantly < 1 (Table I), suggesting that strong purifying selection played a central role in the evolution of Brachyura mitochondrial PCGs to maintain their important functions in OXPHOS. However, the M1 (free-ratio) model of PAML estimated independent ω along all branches of the phylogeny, which was significantly better than the M0 model (p < 0.05, Table I) for all 13 mitochondrial PCGs, suggesting heterogeneous selective pressures on the different lineages.

Several models were compared to further test whether the evidence for positive selection was restricted to the common ancestor and internal branches of Bythograeidae. In the entire Brachyura dataset, the branch-site models were used to test for positive selection in individual codons for the common ancestor of the Bythograeidae (branch a in Fig. 1). Only the ND1 gene was under positive selection in this lineage, and three codons were identified to be under selection by the BEB approach with posterior probabilities ≥ 0.80 (Table II). Next, to evaluate the diversified selection pressures on various lineages within Bythograeidae crabs, branch-site models were used to test for positive selection in individual codons in each lineage across the Bythograeidae phylogeny (branches b–i in Fig. 1). The results show evidence of positive selection in three Bythograeidae-specific lineages; that is, the lineage leading to the last common ancestor of Gandalfus and Austinograea (branch b), the terminal branch of the Segonzacia mesatlantica (branch c) for COX3, and the common ancestor of Gandalfus (branch d) for ND2 and ND5 (Table III; Supplementary Table II).

Table II.- Selective pressure analyses (Branch-site model) of mitochondrial PCGs (13 genes) in all crabs dataset and evidence of positive selection of ND1 in the ancestral of Bythograeidae.

*Position of the site in the amino acid sequence (Integer); followed by the posterior probabilities for site (Decimals).

*Position of the site in the amino acid sequence (Integer); followed by the posterior probabilities for site (Decimals).

The Bythograeidae-only datasets were also used to test whether specific codons in the mitochondrial PCGs had been subjected to positive selection. The site model (M8 vs. M8a) results detected no significant differences in the estimates using the model incorporating selection (M8) and the neutral model (M8a) for the mitochondrial PCGs (Supplementary Table III). The FEL and REL implemented in Datamonkey which have the advantage of improving the estimate of the ω ratio incorporated variation in the rate of synonymous substitution (Pond and Frost, 2005). A total of 12 sites (FEL: 1 in ND1, 1 in ND4, 2 in ND5; REL: 4 in ATP6, 1 in ATP8, 1 in COX1, 1 in CYTB, 1 in ND4) were identified at seven genes in the Bythograeidae crabs. In addition, 43.75% (14/32) of the putative positively selected sites from seven genes (three in ATP6, one in COX1, three in COX3, one in ND1, two in ND2, one in ND4, and three in ND5) were undergoing radical changes using the protein-level approach performed in TreeSAAP (Supplementary Table IV), which provided additional evidence for positive selection on Bythograeidae crabs.

A positively selected gene is usually more reliable if it can be supported by two or more different methods (Wang et al., 2015). In the present study, COX3 and ND2 were detected using PAML and Datamonkey, whereas five genes (ATP6, ATP8, COX1, CYTB, and ND4) underwent positive selection according to Datamonkey and TreeSAAP. ND1 and ND5 were detected to have undergone strong positive selection according to the three methods (i.e., PAML, Datamonkey, and TreeSAAP). Thus, nine mitochondrial PCGs (ATP6, ATP8, COX1, COX3, CYTB, ND1, ND2, ND4, and ND5) were validated by at least two methods (Supplementary Table IV) and were regarded as positively selected genes.

We mapped all radical amino acid sites onto secondary and tertiary structures to gain insight into the functional significance of the putatively selected sites. Twenty-seven sites (three in ATP6, two in ATP8, one in COX1, three in COX3, 3 in CYTB, three in ND1, three in ND2, two in ND4, and seven in ND5) in nine positively selected genes were situated in the protein binding region (Supplementary Table IV). In addition, a large number of radical changes in amino acids of the positively selected genes were located primarily in the helical transmembrane region (Supplementary Table IV). These results show that most positively selected and radical change sites were located within or close to the functional domain (Supplementary Fig. S1).

Deep-sea vent crabs (Crustacea, Decapoda, Bythograeidae) are a good model to study physiological adaptation to extreme physical and chemical conditions (Kadar and Powell, 2006). Vent crabs are characterized by specific morphological traits, such as a white body and reduced eyestalks with vestigial corneas to survive in the extreme environment (Min et al., 2017). According to previous studies, organisms living in hydrothermal vents have developed enhanced oxygen extraction capacity and specific ways to process sulfide (Williams, 1980; Sanders et al., 1988; Segonzac et al., 1993;Hourdez and Lallier, 2007). However, the genetic basis of these adaptations to the hydrothermal vent environment remains poorly explored. Therefore, the present study is the first systematic survey of all 13 mitochondrial PCGs in vent crabs to assess signs of positive selection and to determine whether environmental tolerance has influenced the evolution of these genes.

The ω values that we obtained from the one-ratio model were significantly less than 1, providing support for the expected presence of purifying selection acting on all mitochondrial genes to maintain important functions during energy metabolism (Mamirova et al., 2007; Popadin et al., 2007; Stewart et al., 2008). The signal for positive selection is usually swamped by continuous negative selection that occurs on most sites in a gene sequence because positive selection mainly acts on only a few sites and for a short period of evolutionary time (Zhang et al., 2005; Shen et al., 2010). These reasons may partly explain why we did not detect positive selection using the site models.

Table IV.- Positive selection in mitochondrial PCGs based on the analysis of the three methods (PAML, DATAMONKEY, or TreeSAAP).

+, positive selection was detected. —, no positive selection was detected. *results are listed in Supplementary Table IV.

Branch-site models are powerful for distinguishing positive selection from the relaxation of purifying selection (Zhang et al., 2005). Thus, four genes (COX3, ND1, ND2, and ND5) showed evidence for significant positive selection according to the LRTs in the branch-site model, and a series of codons were identified as candidate sites that had undergone positive selection (posterior probability ≥ 80%; Tables II, III). According to physiological and morphological researches, Segonzacia, Gandalfus and Austinograea species haemolymph possess higher buffer capacity than littoral species (Tsuchida and Fujikura, 2000; Chausson et al., 2004; Hamasaki et al., 2010). Positive selection of these genes in three Bythograeidae-specific lineages (Table III; Supplementary Table II) might therefore play an important role in promoting their adaptation to hydrothermal vent environment. In addition, adaptive evolution was also supported by evidence that the positively selected sites were detected in Datamonkey and TreeSAAP. Nine genes were determined to have undergone positive selection (Table IV). Particularly, the high proportion of putatively selected codons was localized on or near to the important functional regions (protein binding region and helical transmembrane region) in the mitochondrial protein structure. These nucleotide- and protein-level results indicate extensive adaptive evolution of the mitochondrial PCGs, and suggest that positive selection may be the major driving force for the evolution of mitochondrial PCGs in vent crabs.

Mitochondrial OXPHOS sustains organelle function and plays a central role in cellular energy metabolism (Koopman et al., 2013). The OXPHOS system consists of five multisubunit complexes (Complex I–V), All OXPHOS complexes except complex II, which is exclusively derived from the nuclear DNA, contain subunits that are encoded by both nuclear and mitochondrial genes (Distelmaier et al., 2009; Koopman et al., 2010, 2013). Complex I is the largest OXPHOS enzyme and is central to energy transformation in many prokaryotes and most eukaryotes (Brandt, 2006; Treberg and Brand, 2011). Complex I requires a set of 14 evolutionary conserved “core subunits”, consisting of seven mtDNA-encoded ND subunits and seven nDNA-encoded subunits to oxidize NADH to NAD+ and donates the released electrons to the electron carrier coenzyme Q10 (Koopman et al., 2010; Hirst, 2011). Therefore, the observed positive selection in the ND1, ND2, ND3, and ND5 genes in vent crabs suggests enhanced capability for catalysis using complex I. Complex III contains 11 subunits, one of which is encoded by the mtDNA (CYTB) (Koopman et al., 2013). Complex IV consists of 14 subunits, three of which are mtDNA-encoded (COX1, COX2, and COX3) and catalyze electrons donated to molecular oxygen (O2) to form water (Koopman et al., 2013). About 95% of the O2 organisms breathe is consumed by this complex (Fergusonmiller et al., 2012). The positive selection findings at these genes (CYTB, COX1, and COX3) suggest that vent crabs may have adaptively enhanced oxygen utilization during the OXPHOS process. Complex V (ATP synthase) couples proton flow from the intermembrane space back to the matrix by converting ADP and inorganic phosphate to ATP (Koopman et al., 2013). Complex V is comprised of 19 subunits, two of which are encoded by mtDNA (ATP6 and ATP8). A series of positively selected sites observed in the ATP6 and ATP8 genes was localized in or near the functional regions on the protein structure (Supplementary Table IV), indicating that ATP synthesis was advanced to a certain degree. In summary, these results explain why Bythograeidae crabs can live in hydrothermal vents, and suggest that these crabs might have acquired an enhanced capacity for energy metabolism in an extreme hypoxic environment. Of course, further functional assays are needed to confirm these associations in the future.

Sulfide combines with cytochrome c to suffocate animals to death at a low concentration. Thus, the high concentration of sulfide in the hydrothermal vents is another challenge for organisms. According to previous studies, release of cytochrome c from the mitochondria is a key initial step in the apoptotic process, and cytochrome c oxidase (COX) of the respiratory chain plays a key role in reducing toxicity (Ott et al., 2002; Min and Xu, 2007). Considering the function of COX1 and COX3 in reducing toxicity, positive selection of these genes suggests that vent crabs have evolved an enhanced capacity for inhibiting damage from toxic substances. Similar results were also found in a study in which the COX1 and COX2 mutations helped alvinocaridid shrimp resist sulfide deposits and toxic substances (Wang et al., 2017). Therefore, parallel/convergent studies of these toxic-related genes are necessary in these two species groups in the future.

The present study comprehensively investigated 13 mitochondrial PCGs in Bythograeidae crabs. Significant positive selection was examined at nine PCGs, specific for ND1 and ND5, and strong positive selection was detected by the three methods (i.e., PAML, Datamonkey, and TreeSAAP). These results are well matched with the complex adaptation of Bythograeidae crabs to the hydrothermal vent environment. We anticipate that future studies will elucidate the physiology and biochemistry of these PCGs and seek whether a similar evolutionary mechanism underlies animals inhabiting extreme environments.

This study was funded by the National Natural Science Foundation of China (Grant Number 31702014), and Doctoral Scientific Research Foundation of Yancheng Teachers University to ZFW, and Open Foundation of Jiangsu Key Laboratory for Bioresources of Saline Soils (Grant Number JKLBS2016007).

There is supplementary material associated with this article. Access the material online at: http://dx.doi.org/10.17582/journal.pjz/2019.51.5.1721.1731

The authors declare that there is no conflict of interests regarding the publication of this article.