Chromosome-level genome assembly of hadal snailfish reveals mechanisms of deep-sea adaptation in vertebrates

  1. Wenjie Xu
  2. Chenglong Zhu
  3. Xueli Gao
  4. Baosheng Wu
  5. Han Xu
  6. Mingliang Hu
  7. Honghui Zeng
  8. Xiaoni Gan
  9. Chenguang Feng
  10. Jiangmin Zheng
  11. Jing Bo
  12. Li-Sheng He
  13. Qiang Qiu
  14. Wen Wang
  15. Shunping He  Is a corresponding author
  16. Kun Wang  Is a corresponding author
  1. School of Ecology and Environment, Northwestern Polytechnical University, China
  2. Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences, China
  3. Key Laboratory of Aquatic Biodiversity and Conservation, Institute of Hydrobiology, Chinese Academy of Sciences, China

Abstract

As the deepest vertebrate in the ocean, the hadal snailfish (Pseudoliparis swirei), which lives at a depth of 6,000–8,000 m, is a representative case for studying adaptation to extreme environments. Despite some preliminary studies on this species in recent years, including their loss of pigmentation, visual and skeletal calcification genes, and the role of trimethylamine N-oxide in adaptation to high-hydrostatic pressure, it is still unknown how they evolved and why they are among the few vertebrate species that have successfully adapted to the deep-sea environment. Using genomic data from different trenches, we found that the hadal snailfish may have entered and fully adapted to such extreme environments only in the last few million years. Meanwhile, phylogenetic relationships show that they spread into different trenches in the Pacific Ocean within a million years. Comparative genomic analysis has also revealed that the genes associated with perception, circadian rhythms, and metabolism have been extensively modified in the hadal snailfish to adapt to its unique environment. More importantly, the tandem duplication of a gene encoding ferritin significantly increased their tolerance to reactive oxygen species, which may be one of the important factors in their adaptation to high-hydrostatic pressure.

eLife assessment

This important study advances our understanding of the potential mechanisms of deep-sea adaptation and sheds light on the evolutionary history of hadal snailfish. Through comparative genomic analysis, the authors provide convincing evidence and propose hypotheses on the timing of trench colonization, population structure, and adaptations to the hadal snailfish genome in response to their environment.

https://doi.org/10.7554/eLife.87198.3.sa0

Introduction

Since its capture in 2014, the hadal snailfish (Pseudoliparis swirei) has captured the attention of biologists and public as the deepest known vertebrate on Earth (Fujii et al., 2010; Gerringer et al., 2017a; Gerringer et al., 2021b). The hadal zone, 6 km below the sea surface where this species lives, and after which the species is named, is characterized by high-hydrostatic pressure (HHP), complete darkness, and barrenness (Jamieson, 2015; Jamieson, 2011). This special organism that survives and thrives in the hadal zone provides us with an unusual case of adaptation to an extreme environment. After several years of anatomical and genomic research (Gerringer, 2019; Gerringer et al., 2017b; Wang et al., 2019), it is now known that the hadal snailfish has degraded vision and a lack of melanin in its skin (Wang et al., 2019). Previous studies have also revealed that having more unsaturated fatty acids in the cell membrane (Cossins and MacDonald, 1984; Fang et al., 2000) and high levels of trimethylamine N-oxide (TMAO) in the body (Yancey et al., 2014) may have played important roles in resistance to HHP. However, much remains to be discovered about this species; our lack of knowledge can be categorized into three principal areas.

The first concerns the origin of hadal snailfish. They have been observed in several trenches in the northwest Pacific Ocean, including the Mariana (Wang et al., 2019), Yap (Gerringer et al., 2021a), Kuril–Kamchatka (Gerringer et al., 2017b), and Japan (Gerringer et al., 2021a) trenches. The question arises: Do they have the ability to migrate across trenches, or do they enter the hadal zone independently? Furthermore, it has been shown that the divergence time between hadal snailfish and Tanaka’s snailfish (a closely related species distributed in shallow areas) is about 20 million years ago (Mya) Wang et al., 2019; but when did the hadal snailfish enter the hadal zone and how long did they take to complete their adaptation to this ecological niche? The second aspect concerns its morphological and physiological characteristics. For example, in such a dark environment, what do hadal snailfish rely on to sense the world: is it smell, taste, or something else? Do they still have circadian rhythms in the absence of sunlight? Does darkness and HHP have any effect on their behavior? The last area concerns the mechanisms by which they tolerate HHP. If unsaturated fatty acids and TMAO are common in marine fish, why have only a few species such as hadal snailfish been observed to reach such depths? Some studies suggest that certain genetic alterations may confer tolerance to HHP (Wang et al., 2019), but if so, how do these alterations help hadal snailfish to adapt to this environment, and which alterations are most critical?

Unfortunately, these questions have not been well resolved because it is difficult for us to make long-term observations of deep-sea organisms in situ. During 2018–2019, we collected multiple samples of hadal snailfish from the Mariana Trench and Tanaka’s snailfish from the Yellow Sea. Based on more data and more refined genome, we have been able to trace the genomic signals left by adaptive evolution in an attempt to more fully understand the evolutionary processes and key changes in this special organism.

Results

Improved genome assembly for Mariana hadal snailfish

A total of four hadal snailfish (P. swirei) and four Tanaka’s snailfish (Liparis tanakae) individuals were collected for this study (Supplementary file 1). Using a combination of Oxford Nanopore Technologies (ONT) long reads, Beijing Genomics Institute (BGI) short reads, and Hi-C sequencing technologies, we generated a chromosome-level genome assembly for hadal snailfish (Supplementary file 2). The genome assembly comprised 1,173 contigs (total length = 626.44 Mb, contig N50 = 4.22 Mb), organized into 24 chromosomes with an anchoring rate of 98.24%. The new assembly filled 1.26 Mb of gaps that were present in our previous assembly and have a much higher level of genome continuity and completeness (with complete BUSCOs of 96.0% [Actinopterygii_odb10 database]) than the two previous assemblies (Figure 1—figure supplements 1 and 2; Supplementary file 3 and 4; Mu et al., 2021; Wang et al., 2019). Moreover, the genome redundancy caused by mis-assembly is also largely reduced in the new assembly (Figure 1—figure supplement 3), which ensures the reliability of the subsequent analysis. Meanwhile, we generated a high-quality chromosomal-level genome assembly for Tanaka’s snailfish for a comparative evolutionary study. We noticed that there is no major chromosomal rearrangement between hadal snailfish and Tanaka’s snailfish, and chromosome numbers are consistent with the previously reported MTZ-ancestor (the last common ancestor of medaka, Tetraodon, and zebrafish) (Kasahara et al., 2007), while the stickleback had undergone several independent chromosomal fusion events (Figure 1—figure supplement 4).

Based on the new genome assemblies, we re-examined the genetic changes that occurred in the common ancestor of hadal snailfish in combination with the new resequencing and transcriptome data. After a thorough scan and careful inspection, we identified 51 absent genes, 20 unitary pseudogenes, 21 lineage-specific expanded genes, 33 genes with insertions and deletions (with a length ≥3 amino acids) in coding regions, and 33 de novo-originated new genes (Supplementary file 8-12). Most of them have not been previously reported and we discuss them in the following sections.

Cross-trench distribution and high level of genetic diversity

Combining the eight new sequenced individuals (four hadal snailfish and four Tanaka’s snailfish) with five previously reported individuals (four hadal snailfish and one Tanaka’s snailfish), we have been able to form an initial perspective of the hadal snailfish at the population level. The principal component analysis (PCA), neighbor-joining tree, and genetic clustering analysis show that the eight hadal snailfish individuals can be divided into two populations, the first with seven individuals and the second with one individual (Figure 1A–C; Supplementary file 1). Interestingly, the first population includes samples from both the Mariana and Yap trenches. Using the mitochondrial data, we found that the divergence time of these individuals from different trenches appears to be only about 44,000 years (Figure 1—figure supplement 5). Combined with additional publicly available mitochondrial data, we noticed that the sample from the Kermadec Trench (Gerringer et al., 2017b), about 6400 km away from the Mariana Trench, is also clustered with individuals from the first population, and the divergence time was estimated to be 1.0 Mya (Figure 1D, Figure 1—figure supplement 6). These results suggest that hadal snailfish have successfully spread to multiple trenches in the Pacific Ocean over the course of a million years. And this dispersal may have been caused by population expansion or deep circulation.

Figure 1 with 8 supplements see all
Sampling information of hadal snailfish, and phylogenetic relationships and population structure of resequenced individuals.

(A) Principal component analysis (PCA) of eight hadal snailfish and five Tanaka’s snailfish. PC, principal component; MHS, Mariana hadal snailfish; YHS, Yap hadal snailfish; TS, Tanaka’s snailfish. (B) Neighbor-joining tree analysis of eight hadal snailfish and five Tanaka’s snailfish using SNPs detected in whole-genome resequencing data. (C) Ancestry results from Admixture under the k = 5 model. (D) Maximum likelihood trees constructed with 13 genes encoding mitochondria in these species, where KHS01-KHS03 were constructed using two mitochondrial genes (co1 and cytb) and manually merged with other species. Divergence times are shown in each node, and the color of each branch represents the survival depth of the species. KHS, Notoliparis kermadecensis. (E) Sampling information of hadal snailfish and Tanaka’s snailfish. Blue represents the ocean, the darker the color, the deeper the depth, depth data from GEBCO Compilation Group (2020) GEBCO 2020 Grid (doi:10.5285/a29c5465-b138-234d-e053-6c86abc040b9).

Genetic diversity of hadal snailfish is about 3.48 times higher than Tanaka’s snailfish. The FST between the two species is close to 0.91, indicating a large genetic divergence (Figure 1—figure supplement 7). That most polymorphisms are unique to each species, and that they share only 6% of their SNPs, is consistent with this observation. In addition, we also estimated the demographic history of hadal snailfish at the population level and observed a significant expansion in the last 60,000 years (Figure 1—figure supplement 8).

Preserved rh1 gene suggests rapid adaptation to hadal zone

Based on the mitochondrial data, the closest known species related to the hadal snailfish were found to be from the genera Careproctus, Crystallias, Rhodichthys, and Paraliparis, which contain many species living at approximately 1,000 m depth (Figure 1D; Supplementary file 13). The divergence time between hadal snailfish and these species was estimated to be about 9.9 Mya, close to the formation time (8 Mya) of the deepest trough of the Mariana Trench (Oakley et al., 2009). We might therefore speculate that the ancestor of hadal snailfish adapted to the deep-sea environment around 1,000 m at about 9.9 Mya, and subsequently gradually adapted to greater depths in the formed or forming trench. The upper and lower limits of the time for hadal snailfish to enter the hadal zone were estimated to be 9.9 Mya (divergence time between hadal snailfish and its closest relatives) and 1.0 Mya (divergence time between different hadal snailfish individuals), respectively. Interestingly, the vision-related genes also confirm a rapid adaptation to hadal zone.

Fish that inhabit different depths of the sea rely on different vision-related genes (Musilova et al., 2019). Since light with longer wavelengths is absorbed more quickly than those with shorter wavelengths (except for the shortest UV wavelengths), high-energy light with shorter wavelengths, such as blue, is able to penetrate to greater depths (Figure 2A). The genes responsible for absorbing these shorter wavelengths (sws2 and rh2) are therefore much more important to deep-sea species in the photic zone (above 200 m) than those that absorb longer wavelengths (lws). Similarly, the genes providing monochromatic vision in very dim light (rh1 and gnat1) have been proven to be important for deep-sea species of the disphotic zone (from 200 m to 1,000 m) (Musilova et al., 2019). We noticed that the lws gene (long wavelength) has been completely lost in both hadal snailfish and Tanaka’s snailfish; rh2 (central wavelength) has been specifically lost in hadal snailfish (Figure 2B and C); sws2 (short wavelength) has undergone pseudogenization in hadal snailfish (Figure 2—figure supplement 1); while rh1 and gnat1 (perception of very dim light) is both still present and expressed in the eyes of hadal snailfish (Figure 2D). A previous study has also proven the existence of rhodopsin protein in the eyes of hadal snailfish using proteome data (Yan et al., 2021). The preservation and expression of genes for the perception of very dim light suggest that they are still subject to natural selection, at least in the recent past.

Figure 2 with 1 supplement see all
Alterations in vision-related genes in hadal snailfish.

(A) Different colors of light penetrate the depth of the open ocean. Longer wavelengths (such as red) are absorbed at shallower depths, while shorter wavelengths (such as blue) can penetrate to deeper depths. (B) Genetic alterations in the genes encoding the four major proteins involved in activating the photoresponse of vertebrate photoreceptors in the cone cell and rod cell of hadal snailfish. Opsion: rhodopsin, or its cone equivalent. G-protein: heterotrimeric G-protein, transducin. (C) Gene loss of pde6c and opn1mw4 in hadal snailfish. (D) Log10-transformation normalized counts for DESeq2 (COUNTDESEQ2) of vision-related genes in the eyes of hadal snailfish and Tanaka’s snailfish. * represents genes significantly downregulated in hadal snailfish (corrected p<0.05).

Highly expressed auditory genes

Do hadal snailfish compensate for the lack of vision when perceiving the external environment? The genes associated with the olfactory and auditory systems were investigated using both comparative genomic and transcriptomic methods. While the number of olfactory receptors was largely reduced (Figure 3—figure supplement 1), we found that the majority of the auditory genes were well preserved in hadal snailfish. Many of the auditory genes also tended to be significantly more upregulated in the brain of hadal snailfish than in Tanaka’s snailfish (Figure 3A; Supplementary file 14). The upregulated genes involve many aspects of the auditory system, including the development and tethering of otoliths (Kang et al., 2008; Stooke-Vaughan et al., 2015), the development (Iyer and Groves, 2021; Kozlowski et al., 2005; Riley, 2021; Wang et al., 2008), maturation and maintenance of inner ear hair cells, the development and mechanosensitivity of stereocilia (Cirilo et al., 2021; Kitajiri et al., 2010), and other factors (Giffen et al., 2019; Verdoodt et al., 2021; Figure 3A). Of these, the most significant upregulated gene is tmc1, which encodes transmembrane channel-like protein 1, involved in the mechanotransduction process in sensory hair cells of the inner ear that facilitates the conversion of mechanical stimuli into electrical signals used for hearing and homeostasis (Maeda et al., 2014), and some mutations in this gene have been found to be associated with hearing loss (Kitajiri et al., 2007; Riahi et al., 2014). Interestingly, tmc1 is also found to be the gene with the longest deletion specific to hadal snailfish (11 amino acids) in the regions that are generally highly conserved across vertebrate’s genomes (Figure 3—figure supplement 2); the functional implications of this alteration need further verification.

Figure 3 with 2 supplements see all
High expression and gene expansion of hearing-related genes in hadal snailfish.

(A) Upregulation of auditory-related genes in hadal snailfish brain. Red represents upregulated genes in hadal snailfish. (B) Increased copy number of cldnj in hadal snailfish. The relative positions of genes on chromosomes are indicated by arrows, with arrows to the right representing the forward strand and arrows to the left representing the reverse strand. (C) Nanopore sequencing read depth for all cldnj in hadal snailfish.

Moreover, the gene involved in lifelong otolith mineralization, cldnj, has three copies in hadal snailfish, but only one copy in other teleost species, encoding a claudin protein that has a role in tight junctions through calcium-independent cell-adhesion activity (Figure 3B and C; Hardison et al., 2005). This may be important for hadal snailfish because calcium carbonate, the inorganic component in otoliths, is thought not to accumulate efficiently below the carbonate compensation depth (CCD; >4,000–5,000 m) (Jamieson, 2015). It should be noted that the hadal snailfish survive at depths far beyond the limits of CCD, but their otoliths still maintain densities similar to those of sea-surface species (Gerringer et al., 2021a). In our investigation, we found that the expression of cldnj was not significantly upregulated in the brain of the hadal snailfish than in Tanaka’s snailfish, which may be related to the fact that cldnj is mainly expressed in the otocyst, while the expression in the brain is lower. However, due to the immense challenge in obtaining samples of hadal snailfish, the expression of cldnj in the otocyst deserves more in-depth study in the future. Expansion of cldnj was observed in all resequenced individuals of the hadal snailfish (Supplementary file 10), which provides an explanation for the hadal snailfish breaks the depth limitation on calcium carbonate deposition and becomes one of the few species of teleost in hadal zone.

Circadian rhythm decoupled from sunlight and dark adaptation

There is growing evidence that persistent darkness challenges the physiology and behavior of animals, leading to disrupted circadian rhythms, neurological damage, and depressive-behavioral phenotypes (Fisk et al., 2018). Consistent with previous research in cavefish (Policarpo et al., 2021), we noticed that many of the circadian rhythm genes (per2a, cry1a, cry3, cry5, and gpr19) are lost or have undergone pseudogenization in the hadal snailfish (Figure 4A and B, Figure 4—figure supplement 1). Despite that, we noticed that the essential clock control genes are present and expressed in the hadal snailfish, indicating that the rhythm cycle is retained, although it is likely to have been largely uncoupled from sunlight. Moreover, gpr19 deficiency has been reported to prolong the cycle of circadian locomotor activity rhythms (Yamaguchi et al., 2021), so hadal snailfish may have an extended rhythm cycle like cavefish (Cavallari et al., 2011; Yamaguchi et al., 2021).

Figure 4 with 5 supplements see all
Genetic variation of dark adaptation in hadal snailfish.

(A) Genetic changes involved in light-mediated regulation of the molecular clock in hadal snailfish suprachiasmatic nucleus (SCN) neurons. (B) Pseudogenization of gpr19 (gray) in hadal snailfish. Gene structure (top), alignment of nucleotide and amino acid sequences (middle), and sequencing read depth (bottom; the numbers along the x-axis represent the position of the base at the scaffold) for the gpr19 gene. The premature termination (colored red) of gpr19 is due to 82 nucleotide variants and 22 nucleotide deletions (blue). (C) The deletion of one copy of grpr and another copy of downregulated expression in hadal snailfish. The relative positions of genes on chromosomes are indicated by arrows, with arrows to the right representing the forward strand and arrows to the left representing the reverse strand. The heatmap presented is the average of the normalized counts for DESeq2 (COUNTDESEQ2) in all replicate samples from each tissue. * represents tissue in which the grpr-1 was significantly downregulated in hadal snailfish (corrected p<0.05). (D) Gene loss of tmem263 in hadal snailfish.

In addition, in the teleosts closely related to hadal snailfish, there are usually two copies of grpr encoding the gastrin-releasing peptide receptor; we noticed that in hadal snailfish one of them is absent and the other is barely expressed in brain (Figure 4C), whereas a previous study found that the grpr gene in the mouse suprachiasmatic nucleus (SCN) did not fluctuate significantly during a 24 hr light/dark cycle and had a relatively stable expression (Pembroke et al., 2015; Figure 4—figure supplement 1). It has been reported that grpr-deficient mice, while exhibiting normal circadian rhythms, show significantly increased locomotor activity in dark conditions (Wada et al., 1997; Zhao et al., 2023). We might therefore speculate that the absence of that gene might in some way benefit the activity of hadal snailfish under complete darkness.

It should be noted that the abovementioned missing genes are not sufficient to exhibit the full range of changes that occur in the nervous system of hadal snailfish. Previous studies suggest that HHP suppressed the compound action potential in nerve trunks of fishes from shallow areas but not from deep areas, and can perturb the function of G protein-coupled receptors (Siebenaller and Murray, 1995). From our transcriptome data, we also observed that the brain is one of the most divergent organs regarding expression levels between hadal snailfish and Tanaka’s snailfish (Figure 4—figure supplement 2). Specifically, there are 3,587 upregulated genes and 3,433 downregulated genes in the brain of hadal snailfish compared to Tanaka snailfish, and Gene Ontology (GO) functional enrichment analyses revealed that upregulated genes in the hadal snailfish are associated with cilium, DNA repair, and microtubule-based movement, while downregulated genes are enriched in membranes, GTP-binding, proton transmembrane transport, and synaptic vesicles (Supplementary file 15). In line with this observation, one of our previous studies showed that zebrafish brains have the highest number of differentially expressed genes than the other investigated organs when exposed to HHP (Hu et al., 2022). We also identified 15 de novo new genes in hadal snailfish that are highly expressed in the brain (Figure 4—figure supplement 2). The adaptation of the nervous system to HHP deserves more in-depth study in the future.

Possible survival strategy of storing energy

In a previous study, it was noticed that the individual hadal snailfish we investigated retained a large amount of intact food in its stomach and had larger eggs than might otherwise be expected (Gerringer et al., 2017b; Wang et al., 2019). It appears that the hadal snailfish have a survival strategy of storing energy, which is often found in species that need to cope with occasional starvation. Here we find another clue that hints at the existence of this possibility: the pseudogenization of the gene gpr27 in hadal snailfish (Figure 4—figure supplement 3). Gpr27 is a G protein-coupled receptor, belonging to the family of cell surface receptors, involved in various physiological processes and expressed in multiple tissues including the brain, heart, kidney, and immune system. It has been reported that the knockout of gpr27 increases the expression of key enzymes in the carnitine shuttle complex (Nath et al., 2020), especially cpt1, which is essential for the β-oxidation of lipid metabolism. The transcriptome data further confirm that the gene cpt1 is significantly upregulated in the liver of hadal snailfish. As lipid mobilization is thought to be a common metabolic response to short-term starvation in fish (Liao et al., 2017), the inactivation of gpr27 could help hadal snailfish to better survive periods of food deficiency. Although previous surveys have shown that various types of organisms live in the hadal zone, and that the hadal snailfish can survive by eating amphipods and occasionally polychaetes and decapod shrimp (Yan et al., 2021), short-term starvation is still possible because although energy sources are limited by complete darkness, but organisms usually tend to ‘over-reproduce.

Reduced bone mineralization

Vitamin D synthesis is dependent on UV light, with phytoplankton being the origin of vitamin D in food (Björn and Wang, 2000). Whether and how vitamin D reaches the hadal zone through various pathways, for instance, as particulate organic matter, is still unknown. By investigating the genes associated with vitamin D metabolic pathways, we found that these genes are well conserved in the genome of hadal snailfish and are similarly expressed in both hadal snailfish and Tanaka’s snailfish (Figure 4—figure supplement 4), suggesting that vitamin D may not be a limiting factor for hadal zone vertebrates.

Nonetheless, micro-CT scans have revealed shorter bones and reduced bone density in hadal snailfish, from which it has been inferred that this species has reduced bone mineralization (Gerringer et al., 2021a); this may be a result of lowering density by reducing bone mineralization, allowing to maintain neutral buoyancy without expending too much energy, or it may be a result of making its skeleton more flexible and malleable, which is able to better withstand the effects of HHP. The gene bglap, which encodes a highly abundant bone protein secreted by osteoblasts that binds calcium and hydroxyapatite and regulates bone remodeling and energy metabolism, had been found to be a pseudogene in hadal fish (Wang et al., 2019), which may contribute to this phenotype. Here, we found two more lost genes specific to hadal snailfish, tmem251 and tmem263, that contribute to reduced bone mineralization (Figure 4D, Figure 4—figure supplement 5). These two genes encode transmembrane proteins, and loss-of-function mutations have now been found that may affect bone mineral deposition and thus bone development and body growth (Ain et al., 2021; Wu et al., 2018). Furthermore, many genes that determine chondrocyte differentiation and bone mineralization were found to be differentially expressed in the bones of hadal snailfish and Tanaka’s snailfish. However, it should be noted that this result derives from a single bone sample of a hadal snailfish and needs further verification.

HHP adaptation at cellular levels

HHP exerts broad effects upon cells, including cell membrane fluidity (Casadei et al., 2002; Chong et al., 1983; Kato et al., 2002), protein structure stability (Abe, 2021; Gross and Jaenicke, 1994), and oxidative stress (Aertsen et al., 2005; Moserova et al., 2017). In regard to the effect of cell membrane fluidity, relevant genetic alterations had been identified in previous studies, that is, the amplification of acaa1 (encoding acetyl-CoA acetyltransferase 1, a key regulator of fatty acid β-oxidation in the peroxisome, which plays a controlling role in fatty acid elongation and degradation) may increase the ability to synthesize unsaturated fatty acids (Fang et al., 2000; Wang et al., 2019). As for the stability of the protein structure, previous studies have suggested that the high level of TMAO content could help the marine fishes in resistance to the inhibitory effects of high pressure on numerous proteins (Ma et al., 2014; Yancey et al., 2002; Yancey et al., 2014). We also observed another gene duplication event associated with protein stability. The gene vbp1 (Figure 5—figure supplement 1; Vainberg et al., 1998), encoding prefoldin subunit 3 that promotes protein folding, has two copies in hadal snailfish but one copy in other teleost fishes. But unfortunately, although it is widely known that high pressure leads to the accumulation of reactive oxygen species (ROS) (Abe, 2021; Aertsen et al., 2004; Le et al., 2020), it is still unknown how deep-sea fish cope with this challenge.

We further examined the known ROS-related genes in hadal snailfish, but found that they were not significantly altered in sequence or expression (Figure 5—figure supplement 1). Next, we identified 34 genes that are significantly more highly expressed in all organs of hadal snailfish in comparison to Tanaka’s snailfish and zebrafish, while only 7 genes were found to be significantly more highly expressed in Tanaka’s snailfish using the same criterion (Figure 5—figure supplement 1). The 34 genes are enriched in only one GO category, GO:0000077: DNA damage checkpoint (adjusted p-value: 0.0177). Moreover, 5 of the 34 genes are associated with DNA repair. Interestingly, however, when we analyzed the genes that were both expanded and highly expressed in most tissues, we identified only one gene, fthl27 (encoding a ferritin heavy chain-like protein), which has 14 copies (most of which are tandem duplicates) in hadal snailfish as opposed to 3 copies in Tanaka snailfish (Figure 5A, Figure 5—figure supplement 2). It has also been suggested that ferritin helps control ROS (Orino et al., 2001; Salatino et al., 2019). The expansion of fthl27 was validated in all the eight resequencing individuals by reads mapping (Figure 5—figure supplement 3), indicating that the tandem duplication event occurred at least before the differentiation of these individuals. To test whether the fthl27 can resist oxidative stress, we cultured 293T cells with or without fthl27-overexpression plasmid in cell culture medium supplemented with H2O2 or ferric ammonium citrate (FAC) for 4 hr, and subsequently measured intracellular ROS levels as well as cell viability. The results showed that the intracellular ROS levels of fthl27-overexpression cells were significantly lower than that of the control group (Figure 5B, Figure 5—figure supplement 4). Meanwhile, the fthl27-overexpression cells were also found to had significantly higher cell viability (Figure 5C). Therefore, we hypothesize that the expansion and high expression of this gene may be an important mechanism of resistance to HHP induced ROS in hadal snailfish.

Figure 5 with 5 supplements see all
High-hydrostatic pressure adaptation of molecules and cells in hadal snailfish.

(A) The position of the gene copies of fthl27 in Tanaka’s snailfish and hadal snailfish, and the expression of each gene copy in each tissue. The gene expression presented is the average of transcripts per million (TPM) in all replicate samples from each tissue. (B) Reactive oxygen species (ROS) levels were confirmed by redox-sensitive fluorescent probe using DCFH-DA molecular probe in 293T cell culture medium with or without fthl27-overexpression plasmid added with H2O2 or ferric ammonium citrate (FAC) for 4 hr. Images are merged from bright-field images with fluorescent images using ImageJ, while the mean fluorescence intensity (MFI) is also calculated using ImageJ. Green, cellular ROS. Scale bars = 100 μm. (C) Fluorescence intensities of AlamarBlue for 293T cells with or without the fthl27-overexpression plasmid in cell culture medium supplemented with H2O2 or FAC for 4 hr (n=5). Vertical coordinates are relative fluorescence units, and higher relative fluorescence units indicate stronger cell viability. Significantly higher AlamarBlue fluorescence for treated cells compared to control cells (p<0.001) is indicated by an asterisk (***).

Discussion

The more sequenced individuals provide us with more details about the evolutionary history about the hadal snailfish. For example, given that the divergence time of the hadal snailfish and the other species of the family Liparidae living at a depth of 1,000 m was about 9.9 Mya, and the divergence time between different sequenced hadal snailfish individuals was about 1.1 Mya, it is known that the hadal snailfish entered the hadal zone between 1.1 and 9.9 Mya. Then consider the fact that the genes that are responsible for detecting light in dark environment are well preserved in the hadal snailfish, it is likely that this species have only entered a completely light-free environment in the last millions of years, after the full completion of the Mariana Trench (Oakley et al., 2009). In addition, the phylogenetic relationships between different individuals clearly indicate that they have successfully spread to different trenches within 1.0 Mya (Figure 1—figure supplement 6).

The comparative genomic analysis revealed that the complete absence of light had a profound effect on the hadal snailfish. In addition to the substantial loss of visual genes and loss of pigmentation, many rhythm-related genes were also absent, although some rhythm genes were still present. The gene loss may not only come from relaxation of natural selection, but also for better adaptation. For example, the grpr gene copies are absent or downregulated in hadal snailfish, which could in turn increase their activity in the dark, allowing them to survive better in the dark environment (Wada et al., 1997). The loss of gpr27 may also increase the ability of lipid metabolism, which is essential for coping with short-term food deficiencies (Nath et al., 2020).

The most interesting question about the hadal snailfish is why this is currently one of the very few observed vertebrate species capable of surviving and reproducing at such depths. TMAO, which is able to maintain protein function under high pressure, is thought to be a limiting factor in determining the depth at which fish can survive (Yancey et al., 2014). Results from our previous analysis suggested that positive selection on fmo3 may be important in promoting TMAO synthesis; but after an updated genome and rigorous FDR correction, we found that these genes had not undergone any particular or significant alteration in the amino acid sequences, although this does not exclude the possibility of a small number of amino acids (three species-specific mutations in all the five copies) having an effect on enzyme activity (Figure 5—figure supplement 5). Since the expression of all five copies of fmo3 is similar among the various tissues of Tanaka’s snailfish and hadal snailfish (Supplementary file 16), it seems more likely that the mechanism associated with TMAO degradation is altered in hadal snailfish, although we do not yet have any additional evidence to support this because the genes associated with TMAO degradation are still unclear.

However, the levels of TMAO are not sufficient for us to understand why only the hadal snailfish can tolerate such HHP since this substance is widely present in marine fishes. In contrast, the tandem duplication events of two genes may play a more critical role in the adaptation of the hadal snailfish. The first event is the tandem duplication of cldnj, a gene essential for otolith formation (Figure 3B; Hardison et al., 2005). The two more copies may help the hadal snailfish to maintain the densities of their otoliths far beyond the limits of CCD. Since the dissolution rate of calcium carbonate increases with higher pressure, the otoliths stability may be one of the reasons limiting fishes to dive to even deeper regions. The second event is the massive expansion and high expression of fthl27. Our cellular experiments proved that this gene could help cells to resist ROS burst and protect it from various damages caused by oxidative stress. Meanwhile, comparative transcriptomic analysis observed multiple genes associated with DNA repair are significantly more highly expressed in all tissues of hadal snailfish than in other fishes, which also coincides with the presence of oxidative stress (Figure 5—figure supplement 1).

In summary, we provide chromosome-level genomes of hadal snailfish and Tanaka’s snailfish, as well as additional transcriptome and resequencing data. We report here further advances in our understanding of the origin, specific characteristics, and adaptive mechanisms of the hadal snailfish.

Materials and methods

Sample collection and identification

Request a detailed protocol

All the experiments in this study were conducted in accordance with the preapproved guidelines of the Ethics Committee of the Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences (Sanya, China). The hadal snailfish samples were collected form one site in the Mariana Trench (142°26′E, 11°07′N) at depth of 7,254 m using the deep-sea lander Tianya with a surfacing time of 3 hr (Supplementary file 1). These specimens were identified as conspecific with P. swirei by morphological observations. Tanaka’s snailfish specimens were collected in the southern Yellow Sea in 2018 and identified as L. tanakae on the basis of morphological observations.

Genome sequencing and assembly

Request a detailed protocol

Genomic DNA was extracted from the muscle of four hadal snailfish collected from the Mariana Trench and four Tanaka’s snailfish collected from the southern Yellow Sea. We generated a total of 47.8 gigabases (Gb) of Nanopore reads, 148.6 Gb of BGI short reads, and 123.3 Gb of Hi-C reads for hadal snailfish; and 39.0 Gb of Nanopore reads, 130.3 Gb of BGI short reads, and 99.5 Gb of Hi-C reads for Tanaka’s snailfish.

The genome sizes of hadal snailfish and Tanaka’s snailfish were estimated by k-mer distribution analysis (K = 27) of SOApec v2 (Luo et al., 2012) to be 633.2 Mb and 539.9 Mb, respectively. We then assembled the hadal snailfish and Tanaka’s snailfish genomes based on the filtered Nanopore sequencing data using wtdbg2 v2.4.1 (Ruan and Li, 2020) based on default parameters, followed by two rounds of error correction using NextPolish v1.0 (Hu et al., 2020) based on the filtered BGI sequencing data, and finally assembled them into chromosomal versions using 3D-DNA v180114 (Dudchenko et al., 2017) based on Hi-C data. Finally, BUSCO v 4.1.2 (Manni et al., 2021) was used with the library ‘actinopterygii_odb10’ to analyze and evaluate the completeness of the gene set in our draft genome.

Transcriptome sequencing

Request a detailed protocol

A total of 11 transcriptomes from 6 tissues (eye, stomach, heart, liver, muscle, skin) were extracted from three hadal snailfish, while a total of 26 transcriptomes from 10 tissues (brain, spinal cord, eye, bone, cholecyst, stomach, heart, liver, muscle, skin) were extracted from three Tanaka’s snailfish. RNA was subsequently extracted using TRIzol (Invitrogen) and purified using the RNeasy Mini Kit (QIAGEN). Transcriptome reads were obtained from the Illumina HiSeq 2000 sequencing platform. The RNA sequences were filtered using Fastp v0.20 (Chen et al., 2018) and assembled without reference using SPAdes (Bushmanova et al., 2019) with default parameters. Subsequently, TransDecoder(RRID:SCR_017647) v5.5.0 was used to identify coding regions of the transcripts.

Genome annotation

Request a detailed protocol

Both de novo and homology-based predictions were used to identify repetitive elements in hadal snailfish and Tanaka’s snailfish. First, we constructed a de novo transposable element library using RepeatModeler v1.0.11 (Saha et al., 2008), and then used RepeatMasker v4.0.7 (Chen, 2004) to detect repeats. For homologous annotations, the genome sequences were compared with data from Repbase using RepeatMasker v4.0.7 and RepeatProteinMask v1.36 to predict transposable elements. For tandem repeat sequences, we used Tandem Repeats Finder v4.07 (Benson, 1999) to make predictions.

The repeat masked genome was used for the gene annotation. We used a combination of ab initio gene predictions, homologous gene predictions, and direct gene models produced by transcriptome assembly to identify protein-coding genes structure on the genome as follows:

  • Step 1: Augustus v3.2.1 (Stanke et al., 2008) was used to generate ab initio predictions with internal gene models.

  • Step 2: The protein sequences from seven species – medaka, Atlantic cod, flatfish, stickleback, zebrafish, turbot, and fugu – and the transcriptome-predicted protein sequences were used to align genomic sequences with BLAT v. 35 (Supplementary file 6; Kent, 2002).

  • Step 3: The psl files obtained in the previous step were integrated and the protein sequences that were aligned to the overlapping region of the genome were scored and sorted based on the alignment results using a custom script to filter out the best aligned protein sequences in this region. Then, GeneWise v2.4.1 (Birney et al., 2004) was used to predict gene models with the aligned sequences as well as the corresponding query proteins. The custom scripts have been deposited in GitHub (https://github.com/wk8910/bio_tools/tree/master/42.prediction copy archived at Wang, 2021).

  • Step 4: The Evidence Modeler (EVM) v1.1.1 (Haas et al., 2008) was used to integrate the prediction results with different weights for each.

The integrated gene set was translated into amino acid sequences using InterProScan v5 (Jones et al., 2014) to annotate motifs and domains in protein sequences by searching publicly available databases (including Pfam, PRINTS, PANTHER, ProDom, and SMART), and the genes were further annotated using the KEGG databases.

Variant calling using resequencing data

Request a detailed protocol

Short reads of seven Mariana hadal snailfish, one Yap hadal snailfish, and five Tanaka’s snailfish (Supplementary file 1) were mapped to the hadal snailfish genome assembled in this study with BWA v0.7.12-r1039 (Li, 2013); then SAMtools v1.4 (Li et al., 2009) was used to sort and obtain BAM files. To analyze population genetics, we focused on SNPs and small indels (1–10  bp) (Zhang et al., 2021). The SNPs were called using FreeBayes v0.9.10-3-g47a713e (Garrison and Marth, 2012) with parameters ‘--gvcf --min-coverage 5 --limit-coverage 200’ and filtered by following three thresholds: (1) SNPs with missing rate ≤30%; (2) the highest sequencing depth of SNP position <200×; and (3) the lowest sequencing depth for each allele ≥5. Subsequently, we calculated the distribution of heterozygosity in genomewide regions with 500 kb nonoverlapping sliding windows.

Inference of phylogeny history

SNP tree, PCA, and diversity statistics

Request a detailed protocol

PLINK v1.90b6.6 (Chen et al., 2019) was used to perform PCA and other population divergency statistics, including nucleotide diversity and genetic differentiation (FST). A neighbor-joining tree was constructed with PHYLIP v3.697 (Felsenstein, 1993) for paired genetic distance matrices.

Admixture analysis

Request a detailed protocol

Different K values (from 1 to 5) were tested using Admixture v1.3.0 (Alexander et al., 2009) to infer ancestral populations in all hadal snailfish and Tanaka’s snailfish individuals accessions.

Demographic analysis

Request a detailed protocol

The demographic history of hadal snailfish and Tanaka’s snailfish was inferred with pairwise sequential Markovian coalescent (PSMC) (Li and Durbin, 2011) analysis, based on a substitution rate of 1.9174e-09 per generation for hadal snailfish and 5.6790e-09 per generation for Tanaka’s snailfish. The analysis was performed using the following parameters: −N25 −t15 −r5 −p ‘4+25 × 2+4 + 6’. These mutation rates were estimated using r8s v1.81. The generation time is 1 y for Tanaka snailfish and 3 y for hadal snailfish.

Mitochondrial genome phylogenetic reconstruction and divergence time estimation

Request a detailed protocol

The mitochondria of eight hadal snailfish and five Tanaka’s snailfish were assembled using NOVOPlasty v4.3.1 (Dierckxsens et al., 2017) with default parameters and annotated using MITOS (http://mitos2.bioinf.uni-leipzig.de/index.py). Subsequently, mitochondrial data from currently published species of the Liparidae were combined, and nucleic acid sequences of 13 coding genes on mitochondria were aligned with MUSCLE v3.8.425 (Edgar, 2021) using default parameters, and alignments of the coding sequences were generated with pal2nal v14 using default parameters. The maximum likelihood (ML) tree was constructed with RAxML-8.2.12 (Stamatakis, 2014) using the following parameters: -f a -m GTRGAMMA -p 15256 -x 271828 -N 100. Finally, divergence times were estimated using MCMCtree v4.9j (Yang, 2007) with one soft-bound calibration timepoint (snailfish-stickleback: ~32–73 Ma) based on previous studies. For Notoliparis kermadecensis, we combined all the above mitochondrial data and performed the same above analysis based on co1 and cytb gene sequences to obtain the ML tree and divergence times.

Gene loss and duplication

Request a detailed protocol

Here, we applied an improved read mapping-based method to identify gene loss and duplication, which is effective in reducing false positives and false negatives caused by genome assembly and annotation errors as well as multispecies sequence alignments. The custom scripts have been deposited in GitHub (https://github.com/wenjie-xu-nwpu/hadal_snailfish copy archived at Xu, 2023). Although this method may have limitations for identifying gene loss and duplication in species with long divergence times, the divergence times of hadal snailfish and Tanaka’s snailfish are about 20 million years (Wang et al., 2019), and at least 88% of the reads in all hadal snailfish individuals can be well compared to Tanaka’s snailfish genome, indicating that this method is applicable to this study.

For gene loss, the following methods were used for identification. (1) Short reads of eight hadal snailfish and five Tanaka’s snailfish (~30×) were compared to the stickleback and Tanaka’s snailfish genome using BWA v0.7.12-r1039 (Li, 2013) and subsequently sorted using SAMtools v1.4 (Li et al., 2009) to obtain the BAM files. (2) We obtained the reads depth for each sites in the gene coding region based on the annotation information of the reference genome and subsequently classified the depths we had on individual loci into three types (‘HIGH’ for greater than half of the average coverage, ‘LOW’ for less than 3, and ‘MID’ for the rest). We defined sites with ‘HIGH’ for Tanaka’s snailfish and ‘LOW’ for hadal snailfish as hadal snailfish-specific lost sites (SLSs). Then, the genes with SLSs accounting for at least 40% of the coding sequence length were selected as the candidate specific loss genes. (3) The protein sequences of the genes selected in the previous step were used as a reference to search through the genome of hadal snailfish using BLAT v. 35 (Kent, 2002) and predict the gene structure using GeneWise v2.4.1 (Birney et al., 2004) to determine the genes that were completely lost or partially lost in this species. (4) The synteny alignment between the hadal snailfish, Tanaka’s snailfish, and stickleback was plotted for partial or fully lost of the gene.

For gene duplication, the following methods were used for identification. (1) Short reads of eight hadal snailfish and five Tanaka’s snailfish were compared to the stickleback and Tanaka’s snailfish genome using BWA v0.7.12-r1039 (Li, 2013) and subsequently sorted using SAMtools v1.4 (Li et al., 2009) to obtain the BAM files. (2) The homologous sites whose average value of reads depth of all hadal snailfish individuals were greater than 1.5 the average value of the Tanaka’s snailfish individuals were retained and defined as hadal snailfish specific high-copy sites (HCSs). Then, the genes with HCSs accounting for at least 50% of the coding sequence length were selected as the candidate high-copy genes. (3) We searched for the location of this gene on the hadal snailfish genome using BLAT v. 35 (Kent, 2002) and predicted the gene structure using GeneWise v2.4.1 (Birney et al., 2004) to determine its copy number. (4) Finally, the expansion of this gene was determined by constructing a gene tree of the protein sequences of this gene family from nine species: hadal snailfish, Tanaka’s snailfish, medaka, Atlantic cod, flatfish, stickleback, zebrafish, turbot, and fugu.

Identification of unitary pseudogenes

Request a detailed protocol

Unitary pseudogenes are nonfunctional genes that decay at their original location (Tutar, 2012), and we suggest that some missing homologs will exist in hadal snailfish genome as unitary pseudogenes during their adaptation to the special environment of the hadal zone.

We obtained pseudogenes in hadal snailfish by following five steps. (1) Using the stickleback genome sequence as a reference, we performed synteny alignment for three species (hadal snailfish, Tanaka’s snailfish, and stickleback) with Last v956 (Kiełbasa et al., 2011) using the parameters ‘-E 0.05', generating a total of 382  Mb (of which 290 Mb was informative for all species) of one-to-one alignment sequences with Multiz v1 (Blanchette et al., 2004) using the default parameters. (2) Genes with at least 70% of the coding sequences of stickleback or Tanaka’s snailfish present in the MAF and not present in the corresponding regions of hadal snailfish were selected as alternative unitary pseudogene datasets. (3) We used BLAST v2.9.0 (Altschul et al., 1990) to determine if this gene was present in other regions of the hadal snailfish genome. (4) The hadal snailfish corresponding region was extended left and right by 10 kb, and the genes of stickleback and Tanaka’s snailfish were used as references for predict the gene structure using GeneWise v2.4.1 (Birney et al., 2004). (5) Screening for pseudogenes that were consistent in all hadal snailfish individuals.

De novo-originated new genes

Request a detailed protocol

First, the short reads of eight hadal snailfish and five Tanaka’s snailfish were compared to the hadal snailfish genome using BWA v0.7.12-r1039 (Li, 2013) and subsequently sorted using SAMtools v1.4 (Li et al., 2009) to obtain the BAM files. In the second step, we defined a single-sequenced sample with reads depths <10 at a single locus as a deletion locus. Based on the annotation file of hadal snailfish, we screened all Tanaka’s snailfish individuals for genes with deletions >50%. Next, for the genes specifically present in hadal snailfish selected in the previous step, we used BLAST v2.9.0 (Altschul et al., 1990) to align them with the genomes of eight other fishes (Tanaka’s snailfish, medaka, Atlantic cod, flatfish, stickleback, zebrafish, turbot, and fugu) and screened for genes with a matching region <0.4. Genes with transcripts per million (TPM) maxima less than 1 in each tissue of hadal snailfish were filtered out. The fully annotated genes (presence of start and stop codons) in the results were defined as novel genes of hadal snailfish.

Lineage-specific changes in amino acid sequences

Request a detailed protocol

For 17 species – Tanaka’s snailfish, stickleback, pacific bluefin tuna, medaka, platy fish, Atlantic cod, flatfish, zebrafish, turbot, fugu, spotted gar, coelacanth, chicken, mouse, human, brownbanded bamboo shark, and elephant shark (Supplementary file 6) – we identified one-to-one orthologs for each species and hadal snailfish by the Reciprocal Best-Hits (RBH) method, and subsequently selected genes present in 15 species, including hadal snailfish, as conserved gene sets. Next, the protein sequences of the selected genes were aligned using MAFFT v7.471 (Katoh and Standley, 2013), and a custom script was used to select regions that were consistent in other species and had contiguous specificity at sites greater than 3 bp in hadal snailfish, and that had at least 90% sequence identity for each 5 bp region before and after this variant region (Wu et al., 2021). Finally, genes with consistent variants in all hadal snailfish individuals were selected.

We performed protein structure simulation using AlphaFold2 (Cramer, 2021) for the amino acid sequences of target genes in hadal snailfish and Tanaka’s snailfish. Finally, the highest scoring prediction was selected as the best structure and visualized using UCSF Chimera (Pettersen et al., 2004).

Comparative transcriptome analysis

Request a detailed protocol

For the RNA sequences of hadal snailfish and Tanaka’s snailfish, we used Fastp v0.20.0 (Chen et al., 2018) to filter out low-quality and contaminated reads, and then used Hisat2 v 2.1.0 (Kim et al., 2019) to align them to the respective reference genomes. StringTie v1.3.6 (Pertea et al., 2016) was then used to count the number of reads paired for each gene with the help of gene annotation information of the species, and then TPM values were calculated for each gene in both species. Next, we identified 17,281 one-to-one orthologs of hadal snailfish and Tanaka’s snailfish using the RBH method. Subsequently, we identified the genes that were differentially expressed (DEGs) between the same tissues of two species using the R package DESeq2 with |log2 (foldchange)| ≥ 1 and corrected p<0.05. For genes that are upregulated or downregulated in multiple tissues, we first found by stochastic simulation that a gene is differentially expressed between two species in one organ does not affect the probability that this gene is differentially expressed in any other organ. Subsequently, we counted the genes that were upregulated or downregulated in each tissue to obtain a list of genes that were co-altered in multiple tissues.

Cell lines

Request a detailed protocol

We selected human embryonic kidney (HEK) 293T cells as an in vitro model. HEK293T cells were provided by Fourth Military Medical University (Xi'an, China). The cell line was validated by short tandem repeat analysis and validated as negative for mycoplasma. The cells were maintained in DMEM (Gibco, USA) supplemented with 10% FBS (Gibco) and 1% antibiotic antifungal (Gibco) at 37°C, 5% CO2.

Cell culture, transfection, and ROS detection

Request a detailed protocol

HEK293T cells were inoculated in 6-well plates at a density of 4.0 × 105 cells/well. After a day when the cell density reached 50–60%, pcDNA3.1 and pcDNA3.1-fthl27 were transfected into cells with polyethylenimine (PEI), respectively. After 24 hr of transfection, 6-well plates (2 ml total volume of medium per well) were treated with 200 μl of 0.3 mM H2O2 or 0.4 mM FAC, both for 4 hr. Subsequently, the ROS Reactive Oxygen Species Kit (Aladdin R272916-1000T) was used to measure ROS. ROS levels were measured using a DCFH-DA molecular probe, and fluorescence was observed through a fluorescence microscope with an optional FITC filter, with the background removed to observe changes in fluorescence. Bright-field photos and fluorescent photos were merged using ImageJ (RRID:SCR_003070) v1.53t. Subsequently, we added AlamarBlue reagent to the cells in complete medium, incubated them for 2 hr, and then assayed them with a fluorescence zymograph with excitation light wavelength between 530 and 560 nm and emission light wavelength of 590 nm, and recorded the relative fluorescence units.

Data availability

The sequence data files of the hadal snailfish have been deposited in the NCBI BioProject database with accession numbers PRJNA852951 (genome data) and PRJNA855356 (transcriptome data). The genome assembly file is under accession number JANBZZ000000000.

The following data sets were generated
    1. Northwestern Polytechnical University
    (2023) NCBI BioProject
    ID PRJNA855356. Transcriptome sequencing data of Pseudoliparis swirei.
    1. Northwestern Polytechnical University
    (2023) NCBI BioProject
    ID PRJNA852951. Pseudoliparis swirei isolate:HS2019 Genome sequencing and assembly.
The following previously published data sets were used
    1. Northwestern Polytechnical University
    (2018) NCBI BioProject
    ID PRJNA472845. Pseudoliparis amblystomopsis Raw sequence reads.
    1. BGI-SZ
    (2019) NCBI BioProject
    ID PRJNA512070. Pseudoliparis sp. Yap Trench Genome sequencing and assembly.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/gasterosteus_aculeatus/. Stickleback assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/danio_rerio/. Zebrafish assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/lepisosteus_oculatus/. Spotted gar assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/homo_sapiens/. Human assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/oryzias_latipes/. Medaka assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/xiphophorus_maculatus/. Platyfish assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/gadus_morhua/. Atlantic cod assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/scophthalmus_maximus/. Turbot assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/takifugu_rubripes/. Fugu assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/latimeria_chalumnae/. Coelacanth assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/gallus_gallus/. Chicken assembly and gene annotation.
    1. Genome Reference Consortium
    (2019) Ensembl
    ID release-97/fasta/callorhinchus_milii/. Elephant shark assembly and gene annotation.
    1. Tsinghua University
    (2017) NCBI Assembly
    ID GCF_001970005.1. Genome assembly Flounder_ref_guided_V1.0.
    1. Phyloinformatics Unit
    (2018) NCBI Assembly
    ID GCA_003427335.1. Genome assembly Cpunctatum_v1.0.

References

    1. Björn LO
    2. Wang T
    (2000)
    Vitamin D in an ecological context
    International Journal of Circumpolar Health 59:26–32.
  1. Book
    1. Felsenstein J
    (1993)
    PHYLIP (Phylogeny Inference Package)
    Joseph Felsenstein.
    1. Tutar Y
    (2012) Pseudogenes
    Comparative and Functional Genomics 2012:424526.
    https://doi.org/10.1155/2012/424526

Peer review

Reviewer #1 (Public Review):

This manuscript provides an important case study for in-depth research on the adaptability of vertebrates in deep-sea environments. Through analysis of the genomic data of the hadal snailfish, the authors found that this species may have entered and fully adapted to extreme environments only in the last few million years. Additionally, the study revealed the adaptive features of hadal snailfish in terms of perceptions, circadian rhythms and metabolisms, and the role of ferritin in high-hydrostatic pressure adaptation. Besides, the reads mapping method used to identify events such as gene loss and duplication avoids false positives caused by genome assembly and annotation. This ensures the reliability of the results presented in this manuscript. Overall, these findings provide important clues for a better understanding of deep-sea ecosystems and vertebrate evolution.

https://doi.org/10.7554/eLife.87198.3.sa1

Reviewer #2 (Public Review):

This paper presents improved, chromosome level assemblies of the hadal snailfish and Tanaka's snailfish. This is an extension and update of previous work from the group on the hadal snailfish genome. The chromosomal assemblies allow comparisons of genome architecture between a shallow water snailfish and the hadal snailfish to aid inference on timing of colonization of trenches and genomic changes that may have been adaptive for that move.

The comparisons in genomic architecture are compelling: genes present in Tanaka's snailfish that are lost in hadal snailfish that involve whole regions of the genome that no longer map even though adjacent regions do map between the species and across a large evolutionary distance to stickleback. Or genes that are duplicated in hadal snailfish but only appear as single copy in other fishes. The paper focuses on genes in the eye, in hearing, in circadian rythms, and in ROS scavaging. These are all functions that could play a role in adapting to the hadal environment.

The genomic comparisons all seem sound. Stylistically I would prefer if the authors could introduce the gene product and protein function every time they introduce a gene locus. They introduce a gene and general function, but don't usually note what the protein encoded by the gene is and what it's specific function is.

I found the paper generally well written, and the data compelling and creatively displayed.

Upon revision, the authors have commendably addressed all reviewer comments and added a slew of additional analyses. I find the paper stronger, better argued and have no further questions or comments.

https://doi.org/10.7554/eLife.87198.3.sa2

Author response

The following is the authors’ response to the original reviews.

We would like to thank the reviewers for their thoughtful comments and constructive suggestions. Point-by-point responses to comments are given below:

Reviewer #1 (Recommendations For The Authors):

This manuscript provides an important case study for in-depth research on the adaptability of vertebrates in deep-sea environments. Through analysis of the genomic data of the hadal snailfish, the authors found that this species may have entered and fully adapted to extreme environments only in the last few million years. Additionally, the study revealed the adaptive features of hadal snailfish in terms of perceptions, circadian rhythms and metabolisms, and the role of ferritin in high-hydrostatic pressure adaptation. Besides, the reads mapping method used to identify events such as gene loss and duplication avoids false positives caused by genome assembly and annotation. This ensures the reliability of the results presented in this manuscript. Overall, these findings provide important clues for a better understanding of deep-sea ecosystems and vertebrate evolution.

Reply: Thank you very much for your positive comments and encouragement.

However, there are some issues that need to be further addressed.

1. L119: Please indicate the source of any data used.

Reply: Thank you very much for the suggestion. All data sources used are indicated in Supplementary file 1.

1. L138: The demographic history of hadal snailfish suggests a significant expansion in population size over the last 60,000 years, but the results only show some species, do the results for all individuals support this conclusion?

Reply: Thank you for this suggestion. The estimated demographic history of the hadal snailfish reveals a significant population increase over the past 60,000 years for all individuals. The corresponding results have been incorporated into Figure 1-figure supplements 8B.

Author response image 1
Demographic history for 5 hadal snailfish individuals and 2 Tanaka’s snailfish individuals inferred by PSMC.

The generation time of one year for Tanaka snailfish and three years for hadal snailfish.

1. Figure 1-figure supplements 8: Is there a clear source of evidence for the generation time of 1 year chosen for the PSMC analysis?

Reply: We apologize for the inclusion of an incorrect generation time in Figure 1-figure supplements 8. It is important to note that different generation times do not change the shape of the PSMC curve, they only shift the curve along the axis. Due to the absence of definitive evidence regarding the generation time of the hadal snailfish, we have referred to Wang et al., 2019, assuming a generation time of one year for Tanaka snailfish and three years for hadal snailfish. The generation time has been incorporated into the main text (lines 516-517): “The generation time of one year for Tanaka snailfish and three years for hadal snailfish.”.

1. L237: Transcriptomic data suggest that the greatest changes in the brain of hadal snailfish compared to Tanaka's snailfish, what functions these changes are specifically associated with, and how these functions relate to deep-sea adaptation.

Reply: Thank you for this suggestion. Through comparative transcriptome analysis, we identified 3,587 up-regulated genes and 3,433 down-regulated genes in the brains of hadal snailfish compared to Tanaka's snailfish. Subsequently, we conducted Gene Ontology (GO) functional enrichment analysis on the differentially expressed genes, revealing that the up-regulated genes were primarily associated with cilium, DNA repair, protein binding, ATP binding, and microtubule-based movement. Conversely, the down-regulated genes were associated with membranes, GTP-binding, proton transmembrane transport, and synaptic vesicles, as shown in following table (Supplementary file 15). Previous studies have shown that high hydrostatic pressure induces DNA strand breaks and damage, and that DNA repair-related genes upregulated in the brain may help hadal snailfish overcome these challenges.

Author response table 1
GO enrichment of expression up-regulated and down-regulated genes in hadal snailfish brain.
GO termTypeFunctionAdjusted
p-value
Up_-_
regulated
GO:0005929Cellular Componentcilium2.162E-04
GO:0006281Biological ProcessDNA repair0.001
GO:0005515Molecular Functionprotein binding0.002
GO:0005524Molecular FunctionATP binding0.005
GO:0007018Biological Processmicrotubule-based
movement
0.029
Down-
regulated
GO:0016020Cellular Componentmembrane8.799E-13
GO:0016021Cellular Componentintegral component of
membrane
2.504E-06
GO:0005886Cellular Componentplasma membrane3.248E-05
GO:0005525Molecular FunctionGTP binding1.089E-03
GO:0003924Molecular FunctionGTPase activity0.001
GO: 1902600Biological Processproton transmembrane
transport
0.003
GO:0005739Cellular Componentmitochondrion0.009
GO:0005743Cellular Componentmitochondrial inner
membrane
0.010
GO:0008021Cellular Componentsynaptic vesicle0.049
GO:0007399Biological Processnervous system
development
0.049

We have added new results (Supplementary file 15) and descriptions to show the changes in the brains of hadal snailfish (lines 250-255): “Specifically, there are 3,587 up-regulated genes and 3,433 down-regulated genes in the brain of hadal snailfish compared to Tanaka snailfish, and Gene Ontology (GO) functional enrichment analyses revealed that up-regulated genes in the hadal snailfish are associated with cilium, DNA repair, and microtubule-based movement, while down-regulated genes are enriched in membranes, GTP-binding, proton transmembrane transport, and synaptic vesicles (Supplementary file 15).”

1. L276: What is the relationship between low bone mineralization and deep-sea adaptation, and can low mineralization help deep-sea fish better adapt to the deep sea?

Reply: Thank you for this suggestion. The hadal snailfish exhibits lower bone mineralization compared to Tanaka's snailfish, which may have facilitated its adaptation to the deep sea. On one hand, this reduced bone mineralization could have contributed to the hadal snailfish's ability to maintain neutral buoyancy without excessive energy expenditure. On the other hand, the lower bone mineralization may have also rendered their skeleton more flexible and malleable, enhancing their resilience to high hydrostatic pressure. Accordingly, we added the following new descriptions (lines 295-300): “Nonetheless, micro-CT scans have revealed shorter bones and reduced bone density in hadal snailfish, from which it has been inferred that this species has reduced bone mineralization (M. E. Gerringer et al., 2021); this may be a result of lowering density by reducing bone mineralization, allowing to maintain neutral buoyancy without expending too much energy, or it may be a result of making its skeleton more flexible and malleable, which is able to better withstand the effects of HHP.”

1. L293: The abbreviation HHP was mentioned earlier in the article and does not need to be abbreviated here.

Reply: Thank you for the correction. We have corrected the word. Line 315.

1. L345: It should be "In addition, the phylogenetic relationships between different individuals clearly indicate that they have successfully spread to different trenches about 1.0 Mya".

Reply: Thank you for the correction. We have corrected the word. Line 374.

1. It is curious what functions are associated with the up-regulated and down-regulated genes in all tissues of hadal snailfish compared to Tanaka's snailfish, and what functions have hadal snailfish lost in order to adapt to the deep sea?

Reply: Thank you for this suggestion. We added a description of this finding in the results section (lines 337-343): “Next, we identified 34 genes that are significantly more highly expressed in all organs of hadal snailfish in comparison to Tanaka’s snailfish and zebrafish, while only seven genes were found to be significantly more highly expressed in Tanaka’s snailfish using the same criterion (Figure 5-figure supplements 1). The 34 genes are enriched in only one GO category, GO:0000077: DNA damage checkpoint (Adjusted P-value: 0.0177). Moreover, five of the 34 genes are associated with DNA repair.” This suggests that up-regulated genes in all tissues in hadal snailfish are associated with DNA repair in response to DNA damage caused by high hydrostatic pressure, whereas down-regulated genes do not show enrichment for a particular function.

Overall, the functions lost in hadal snailfish adapted to the deep sea are mainly related to the effects of the dark environment, which can be summarized as follows (lines 375-383): “The comparative genomic analysis revealed that the complete absence of light had a profound effect on the hadal snailfish. In addition to the substantial loss of visual genes and loss of pigmentation, many rhythm-related genes were also absent, although some rhythm genes were still present. The gene loss may not only come from relaxation of natural selection, but also for better adaptation. For example, the grpr gene copies are absent or down-regulated in hadal snailfish, which could in turn increased their activity in the dark, allowing them to survive better in the dark environment (Wada et al., 1997). The loss of gpr27 may also increase the ability of lipid metabolism, which is essential for coping with short-term food deficiencies (Nath et al., 2020).”

Reviewer #2 (Recommendations For The Authors):

I have pointed out some of the examples that struck me as worthy of additional thought/writing/comments from the authors. Any changes/comments are relatively minor.

Reply: Thank you very much for your positive comments on this work.

For comparative transcriptome analyses, reads were mapped back to reference genomes and TPM values were obtained for gene-level count analyses. 1:1 orthologs were used for differential expression analyses. This is indeed the only way to normalize counts across species, by comparing the same gene set in each species. Differential expression statistics were run in DEseq2. This is a robust way to compare gene expression across species and where fold-change values are reported (e.g. Fig 3, creatively by coloring the gene name) the values are best-practice.

In other places, TPM values are reported (e.g. Fig 2D, Fig 4C, Fig 5A, Fig 4-Fig supp 4) to illustrate expression differences within a tissue across species. The comparisons look robust, although it is not made clear how the values were obtained in all cases. For example, in Fig 2D the TPM values appear to be from eyes of individual fish, but in Fig 4C and 5A they must be some kind of average? I think that information should be added to the figure legends.

Of note: TPM values are sensitive to the shape of the RNA abundance distribution from a given sample: A small number of very highly expressed genes might bias TPM values downward for other genes. From one individual to another or from one species to another, it is not obvious to me that we should expect the same TPM distribution from the same tissues, making it a challenging metric for comparison across samples, and especially across species. An alternative measure of RNA abundance is normalized counts that can be output from DEseq2. See:

Zhao, Y., Li, M.C., Konaté, M.M., Chen, L., Das, B., Karlovich, C., Williams, P.M., Evrard, Y.A., Doroshow, J.H. and McShane, L.M., 2021. TPM, FPKM, or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI patient-derived models repository. Journal of translational medicine, 19(1), pp.1-15.

If the authors would like to keep the TPM values, I think it would be useful for them to visualize the TPM value distribution that the numbers were derived from. One way to do this would be to make a violin plot for species/tissue and plot the TPM values of interest on that. That would give a visualization of the ranked value of the gene within the context of all other TPM values. A more highly expressed gene would presumably have a higher rank in context of the specific tissue/species and be more towards the upper tail of the distribution. An example violin plot can be found in Fig 6 of:

Burns, J.A., Gruber, D.F., Gaffney, J.P., Sparks, J.S. and Brugler, M.R., 2022. Transcriptomics of a Greenlandic Snailfish Reveals Exceptionally High Expression of Antifreeze Protein Transcripts. Evolutionary Bioinformatics, 18, p.11769343221118347.

Alternatively, a comparison of TPM and normalized count data (heatmaps?) would be of use for at least some of the reported TPM values to show whether the different normalization methods give comparable outputs in terms of differential expression. One reason for these questions is that DEseq2 uses normalized counts for statistical analyses, but values are expressed as TPM in the noted figures (yes, TPM accounts for transcript length, but can still be subject to distribution biases).

Reply: Thank you for your suggestions. Following your suggestions, we modified Fig 2D, Fig 4C, Fig 4-Fig supp 4, and Fig 5-Fig supp 1, respectively. In the differential expression analyses, only one-to-one orthologues of hadal snailfish and Tanaka's snailfish can get the normalized counts output by DEseq2, so we showed the normalized counts by DEseq2 output for Fig 2D, Fig 4C, Fig 4-Fig supp 4, Fig 5-Fig supp 1, and for Fig 5A, since the copy number of fthl27 genes undergoes specific expansion in hadal snailfish, we visualized the ranking of all fthl27 genes across tissues by plotting violins in Fig 5-Fig supp 2.

Author response image 2
Log10-transformation normalized counts for DESeq2 (COUNTDESEQ2) of vision-related genes in the eyes of hadal snailfish and Tanka's snailfish.

* represents genes significantly downregulated in hadal snailfish (corrected P < 0.05).

Author response image 3
The deletion of one copy of grpr and another copy of down-regulated expression in hadal snailfish.

The relative positions of genes on chromosomes are indicated by arrows, with arrows to the right representing the forward strand and arrows to the left representing the reverse strand. The heatmap presented is the average of the normalized counts for DESeq2 (COUNTDESEQ2) in all replicate samples from each tissue. * represents tissue in which the grpr-1 was significantly down-regulated in hadal snailfish (corrected P < 0.05).

Author response image 4
Expression of the vitamin D related genes in various tissues of hadal snailfish and Tanaka's snailfish.

The heatmap presented is the average of the normalized counts for DESeq2 (COUNTDESEQ2) in all replicate samples from each tissue.

Author response image 5
Expression of the ROS-related genes in different tissues of hadal snailfish and Tanaka's snailfish.

The heatmap presented is the average of the normalized counts for DESeq2 (COUNTDESEQ2) in all replicate samples from each tissue.

Author response image 6
Ranking of the expression of individual copies of fthl27 gene in hadal snailfish and Tanaka's snailfish in various tissues showed that all copies of fthl27 in hadal snailfish have high expression.

The gene expression presented is the average of TPM in all replicate samples from each tissue.

Line 96: Which BUSCOs? In the methods it is noted that the actinopterygii_odb10 BUSCO set was used. I think it should also be noted here so that it is clear which BUSCO set was used for completeness analysis. It could even be informally the ray-finned fish BUSCOs or Actinopterygii BUSCOs.

Reply: Thank you for this suggestion. We used Actinopterygii_odb10 database and we added the BUSCO set to the main text as follows (lines 92-95): “The new assembly filled 1.26 Mb of gaps that were present in our previous assembly and have a much higher level of genome continuity and completeness (with complete BUSCOs of 96.0 % [Actinopterygii_odb10 database]) than the two previous assemblies.”

Lines 102-105: The medaka genome paper proposes the notion that the ancestral chromosome number between medaka, tetraodon, and zebrafish is 24. There may be other evidence of that too. Some of that evidence should be cited here to support the notion that sticklebacks had chromosome fusions to get to 21 chromosomes rather than scorpionfish having chromosome fissions to get to 24. Here's the medaka genome paper:

Kasahara, M., Naruse, K., Sasaki, S., Nakatani, Y., Qu, W., Ahsan, B., Yamada, T., Nagayasu, Y., Doi, K., Kasai, Y. and Jindo, T., 2007. The medaka draft genome and insights into vertebrate genome evolution. Nature, 447(7145), pp.714-719.

Reply: Thank you for your great suggestion. Accordingly, we modified the sentence and added the citation as follows (lines 100-105): “We noticed that there is no major chromosomal rearrangement between hadal snailfish and Tanaka’s snailfish, and chromosome numbers are consistent with the previously reported MTZ-ancestor (the last common ancestor of medaka, Tetraodon, and zebrafish) (Kasahara et al., 2007), while the stickleback had undergone several independent chromosomal fusion events (Figure 1-figure supplements 4).”

Line 161-173: "Along with the expression data, we noticed that these genes exhibit a different level of relaxation of natural selection in hadal snailfish (Figure 2B; Figure 2-figure supplements 1)." With the above statment and evidence, the authors are presumably referring to gene losses and differences in expression levels. I think that since gene expression was not measured in a controlled way it may not be a good measure of selection throughout. The reported genes could be highly expressed under some other condition, selection intact. I find Fig2-Fig supp 1 difficult to interpret. I assume I am looking for regions where Tanaka’s snailfish reads map and Hadal snailfish reads do not, but it is not abundantly clear. Also, other measures of selection might be good to investigate: accumulation of mutations in the region could be evidence of relaxed selection, for example, where essential genes will accumulate fewer mutations than conditional genes or (presumably) genes that are not needed at all. The authors could complete a mutational/SNP analysis using their genome data on the discussed genes if they want to strengthen their case for relaxed selection. Here is a reference (from Arabidopsis) showing these kinds of effects:

Monroe, J.G., Srikant, T., Carbonell-Bejerano, P., Becker, C., Lensink, M., Exposito-Alonso, M., Klein, M., Hildebrandt, J., Neumann, M., Kliebenstein, D. and Weng, M.L., 2022. Mutation bias reflects natural selection in Arabidopsis thaliana. Nature, 602(7895), pp.101-105.

Reply: Thank you for pointing out this important issue. Following your suggestion, we have removed the mention of the down-regulation of some visual genes in the eyes of hadal snailfish and the results of the original Fig2-Fig supp 1 that were based on reads mapping to confirm whether the genes were lost or not. To investigate the potential relaxation of natural selection in the opn1sw2 gene in hadal snailfish, we conducted precise gene structure annotation. Our findings revealed that the opn1sw2 gene is pseudogenized in hadal snailfish, indicating a relaxation of natural selection. We have included this result in Figure 2-figure supplements 1.

Author response image 7
Pseudogenization of opn1sw2 in hadal snailfish.

The deletion changed the protein’s sequence, causing its premature termination.

Accordingly, we have toned down the related conclusions in the main text as follows (lines 164-173): “We noticed that the lws gene (long wavelength) has been completely lost in both hadal snailfish and Tanaka’s snailfish; rh2 (central wavelength) has been specifically lost in hadal snailfish (Figure 2B and 2C); sws2 (short wavelength) has undergone pseudogenization in hadal snailfish (Figure 2-figure supplements 1); while rh1 and gnat1 (perception of very dim light) is both still present and expressed in the eyes of hadal snailfish (Figure 2D). A previous study has also proven the existence of rhodopsin protein in the eyes of hadal snailfish using proteome data (Yan, Lian, Lan, Qian, & He, 2021). The preservation and expression of genes for the perception of very dim light suggests that they are still subject to natural selection, at least in the recent past.”

Line 161-170: What tissue were the transcripts derived from for looking at expression level of opsins? Eyes?

Reply: Thank you for your suggestions. The transcripts used to observe the expression levels of optic proteins were obtained from the eye.

Line 191: What does tmc1 do specifically?

Reply: Thank you for this suggestion. The tmc1 gene encodes transmembrane channel-like protein 1, involved in the mechanotransduction process in sensory hair cells of the inner ear that facilitates the conversion of mechanical stimuli into electrical signals used for hearing and homeostasis. We added functional annotations for the tmc1 in the main text (lines 190-196): “Of these, the most significant upregulated gene is tmc1, which encodes transmembrane channel-like protein 1, involved in the mechanotransduction process in sensory hair cells of the inner ear that facilitates the conversion of mechanical stimuli into electrical signals used for hearing and homeostasis (Maeda et al., 2014), and some mutations in this gene have been found to be associated with hearing loss (Kitajiri, Makishima, Friedman, & Griffith, 2007; Riahi et al., 2014).”

Line 208: "it is likely" is a bit proscriptive

Reply: Thank you for this suggestion. We rephrased the sentence as follows (lines 213-215): “Expansion of cldnj was observed in all resequenced individuals of the hadal snailfish (Supplementary file 10), which provides an explanation for the hadal snailfish breaks the depth limitation on calcium carbonate deposition and becomes one of the few species of teleost in hadal zone.”

Line 199: maybe give a little more info on exactly what cldnj does? e.g. "cldnj encodes a claudin protein that has a role in tight junctions through calcium independent cell-adhesion activity" or something like that.

Reply: Thank you for this suggestion. We have added functional annotations for the cldnj to the main text (lines 200-204): “Moreover, the gene involved in lifelong otolith mineralization, cldnj, has three copies in hadal snailfish, but only one copy in other teleost species, encodes a claudin protein that has a role in tight junctions through calcium independent cell-adhesion activity (Figure 3B, Figure 3C) (Hardison, Lichten, Banerjee-Basu, Becker, & Burgess, 2005).”

Lines 199-210: Paragraph on cldnj: there are extra cldnj genes in the hadal snailfish, but no apparent extra expression. Could the authors mention that in their analysis/discussion of the data?

Reply: Thank you for your suggestions. Despite not observing significant changes in cldnj expression in the brain tissue of hadal snailfish compared to Tanaka's snailfish, it is important to consider that the brain may not be the primary site of cldnj expression. Previous studies in zebrafish have consistently shown expression of cldnj in the otocyst during the critical early growth phase of the otolith, with a lower level of expression observed in the zebrafish brain. However, due to the unavailability of otocyst samples from hadal snailfish in our current study, our findings do not provide confirmation of any additional expression changes resulting from cldnj amplification. Consequently, it is crucial to conduct future comprehensive investigations to explore the expression patterns of cldnj specifically in the otocyst of hadal snailfish. Accordingly, we added a discussion of this result in the main text (lines 209-214): “In our investigation, we found that the expression of cldnj was not significantly up-regulated in the brain of the hadal snailfish than in Tanaka’s snailfish, which may be related to the fact that cldnj is mainly expressed in the otocyst, while the expression in the brain is lower. However, due to the immense challenge in obtaining samples of hadal snailfish, the expression of cldnj in the otocyst deserves more in-depth study in the future.”

Lines 225-231: I wonder whether low expression of a circadian gene might be a time of day effect rather than an evolutionary trait. Could the authors comment?

Reply: Thank you for your suggestions. Previous studies have shown that the grpr gene is expressed relatively consistently in mouse suprachiasmatic nucleus (SCN) throughout the day (Figure 4-figure supplements 1) and we hypothesize that the low expression of grpr-1 gene expression in hadal snailfish is an evolutionary trait. We have modified this result in the main text (lines 232-242): “In addition, in the teleosts closely related to hadal snailfish, there are usually two copies of grpr encoding the gastrin-releasing peptide receptor; we noticed that in hadal snailfish one of them is absent and the other is barely expressed in brain (Figure 4C), whereas a previous study found that the grpr gene in the mouse suprachiasmatic nucleus (SCN) did not fluctuate significantly during a 24-hour light/dark cycle and had a relatively stable expression (Pembroke, Babbs, Davies, Ponting, & Oliver, 2015) (Figure 4-figure supplements 1). It has been reported that grpr deficient mice, while exhibiting normal circadian rhythms, show significantly increased locomotor activity in dark conditions (Wada et al., 1997; Zhao et al., 2023). We might therefore speculate that the absence of that gene might in some way benefit the activity of hadal snailfish under complete darkness.”

Author response image 8
Expression of the grpr in a 24-hour light/dark cycle in the mouse suprachiasmatic nucleus (SCN).

Line 253: What is gpr27? G protein coupled receptor?

Reply: We apologize for the ambiguous description. Gpr27 is a G protein-coupled receptor, belonging to the family of cell surface receptors. We introduced gpr27 in the main text as follows (lines 270-273): “Gpr27 is a G protein-coupled receptor, belonging to the family of cell surface receptors, involved in various physiological processes and expressed in multiple tissues including the brain, heart, kidney, and immune system.”

Line 253: Fig4 Fig supp 3 is a good example of pseudogenization!

Reply: Thank you very much for your recognition.

Line 279: What is bglap? It regulates bone mineralization, but what specifically does that gene do?

Reply: We apologize for the ambiguous description. The bglap gene encodes a highly abundant bone protein secreted by osteoblasts that binds calcium and hydroxyapatite and regulates bone remodeling and energy metabolism. We introduced bglap in the main text as follows (lines 300-304): “The gene bglap, which encodes a highly abundant bone protein secreted by osteoblasts that binds calcium and hydroxyapatite and regulates bone remodeling and energy metabolism, had been found to be a pseudogene in hadal fish (K. Wang et al., 2019), which may contribute to this phenotype.”

Line 299: Introduction of another gene without providing an exact function: acaa1.

Reply: We apologize for the ambiguous description. The acaa1 gene encodes acetyl-CoA acetyltransferase 1, a key regulator of fatty acid β-oxidation in the peroxisome, which plays a controlling role in fatty acid elongation and degradation. We introduced acaa1 in the main text as follows (lines 319-324): “In regard to the effect of cell membrane fluidity, relevant genetic alterations had been identified in previous studies, i.e., the amplification of acaa1 (encoding acetyl-CoA acetyltransferase 1, a key regulator of fatty acid β-oxidation in the peroxisome, which plays a controlling role in fatty acid elongation and degradation) may increase the ability to synthesize unsaturated fatty acids (Fang et al., 2000; K. Wang et al., 2019).”

Fig 5 legend: The DCFH-DA experiment is not an immunofluorescence assay. It is better described as a redox-sensitive fluorescent probe. Please take note throughout.

Reply: Thank you for pointing out our mistakes. We corrected the word. Line 1048 and 1151 as follows: “ROS levels were confirmed by redox-sensitive fluorescent probe using DCFH-DA molecular probe in 293T cell culture medium with or without fthl27-overexpression plasmid added with H2O2 or FAC for 4 hours.”

Line 326: Manuscript notes that ROS levels in transfected cells are "significantly lower" than the control group, but there is no quantification or statistical analysis of ROS levels. In the methods, I noticed the mention of flow cytometry, but do not see any data from that experiment. Proportion of cells with DCFH-DA fluorescence above a threshold would be a good statistic for the experiment... Another could be average fluorescence per cell. Figure 5B shows some images with green dots and it looks like more green in the "control" (which could better be labeled as "mock-transfection") than in the fthl27 overexpression, but this could certainly be quantified by flow cytometry. I recommend that data be added.

Reply: Thank you for your suggestions. We apologize for the error in the main text, we used a fluorescence microscope to observe fluorescence in our experiments, not a flow cytometer. We have corrected it in the methods section as follows (lines 651-653): “ROS levels were measured using a DCFH-DA molecular probe, and fluorescence was observed through a fluorescence microscope with an optional FITC filter, with the background removed to observe changes in fluorescence.” Meanwhile, we processed the images with ImageJ to obtain the respective mean fluorescence intensities (MFI) and found that the MFI of the fthl27-overexpression cells were lower than the control group, which indicated that the ROS levels of the fthl27-overexpression cells were significantly lower than the control group. MFI has been added to Figure 5B.

Author response image 9
ROS levels were confirmed by redox-sensitive fluorescent probe using DCFH-DA molecular probe in 293T cell culture medium with or without fthl27-overexpression plasmid added with H2O2 or FAC for 4 hours.

Images are merged from bright field images with fluorescent images using ImageJ, while the mean fluorescence intensity (MFI) is also calculated using ImageJ. Green, cellular ROS. Scale bars equal 100 μm.

Regarding the ROS experiment: Transfection of HEK293T cells should be reasonably straightforward, and the experiment was controlled appropriately with a mock transfection, but some additional parameters are still needed to help interpret the results. Those include: Direct evidence that the transfection worked, like qPCR, western blots (is the fthl27 tagged with an antigen?), coexpression of a fluorescent protein. Then transfection efficiency should be calculated and reported.

Reply: Thank you for your suggestions. To assess the success of the transfection, we randomly selected a subset of fthl27-transfected HEK293T cells for transcriptome sequencing. This approach allowed us to examine the gene expression profiles and confirm the efficacy of the transfection process. As control samples, we obtained transcriptome data from two untreated HEK293T cells (SRR24835259 and SRR24835265) from NCBI. Subsequently, we extracted the fthl27 gene sequence of the hadal snailfish, along with 1,000 bp upstream and downstream regions, as a separate scaffold. This scaffold was then merged with the human genome to assess the expression levels of each gene in the three transcriptome datasets. The results demonstrated that the fthl27 gene exhibited the highest expression in fthl27-transfected HEK293T cells, while in the control group, the expression of the fthl27 gene was negligible (TPM = 0). Additionally, the expression patterns of other highly expressed genes were similar to those observed in the control group, confirming the successful fthl27 transfection. These findings have been incorporated into Figure 5-figure supplements 3.

Author response image 10
(B) Reads depth of fthl27 gene in fthl27-transfected HEK293T cells and 2 untreated HEK293T cells (SRR24835259 and SRR24835265) transcriptome data.

(C) Expression of each gene in the transcriptome data of fthl27-transfected HEK293T cells and 2 untreated HEK293T cells (SRR24835259 and SRR24835265), where the genes shown are the 4 most highly expressed genes in each sample.

Lines 383-386: expression of DNA repair genes is mentioned, but not shown anywhere in the results?

Reply: Thank you for your suggestions. Accordingly, we added a description of this finding in the results section (lines 337-343): “Next, we identified 34 genes that are significantly more highly expressed in all organs of hadal snailfish in comparison to Tanaka’s snailfish and zebrafish, while only seven genes were found to be significantly more highly expressed in Tanaka’s snailfish using the same criterion (Figure 5-figure supplements 1). The 34 genes are enriched in only one GO category, GO:0000077: DNA damage checkpoint (Adjusted P-value: 0.0177). Moreover, five of the 34 genes are associated with DNA repair.”. And we added the information in the Figure 5-figure supplements 1C.

Author response image 11
Genes were significantly more highly expressed in all tissues of the hadal snailfish compared to Tanaka's snailfish, and 5 genes (purple) were associated with DNA repair.
https://doi.org/10.7554/eLife.87198.3.sa3

Article and author information

Author details

  1. Wenjie Xu

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Chenglong Zhu, Xueli Gao, Baosheng Wu, Han Xu, Mingliang Hu and Honghui Zeng
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6240-8472
  2. Chenglong Zhu

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Visualization, Writing – original draft
    Contributed equally with
    Wenjie Xu, Xueli Gao, Baosheng Wu, Han Xu, Mingliang Hu and Honghui Zeng
    Competing interests
    No competing interests declared
  3. Xueli Gao

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Validation, Methodology, X.G. performed experiments with fthl27
    Contributed equally with
    Wenjie Xu, Chenglong Zhu, Baosheng Wu, Han Xu, Mingliang Hu and Honghui Zeng
    Competing interests
    No competing interests declared
  4. Baosheng Wu

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Resources, Investigation
    Contributed equally with
    Wenjie Xu, Chenglong Zhu, Xueli Gao, Han Xu, Mingliang Hu and Honghui Zeng
    Competing interests
    No competing interests declared
  5. Han Xu

    Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences, Sanya, China
    Contribution
    Resources, Investigation
    Contributed equally with
    Wenjie Xu, Chenglong Zhu, Xueli Gao, Baosheng Wu, Mingliang Hu and Honghui Zeng
    Competing interests
    No competing interests declared
  6. Mingliang Hu

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Visualization
    Contributed equally with
    Wenjie Xu, Chenglong Zhu, Xueli Gao, Baosheng Wu, Han Xu and Honghui Zeng
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9018-6715
  7. Honghui Zeng

    Key Laboratory of Aquatic Biodiversity and Conservation, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China
    Contribution
    Resources
    Contributed equally with
    Wenjie Xu, Chenglong Zhu, Xueli Gao, Baosheng Wu, Han Xu and Mingliang Hu
    Competing interests
    No competing interests declared
  8. Xiaoni Gan

    Key Laboratory of Aquatic Biodiversity and Conservation, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China
    Contribution
    Resources
    Competing interests
    No competing interests declared
  9. Chenguang Feng

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Visualization, Writing – original draft
    Competing interests
    No competing interests declared
  10. Jiangmin Zheng

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Formal analysis, Visualization
    Competing interests
    No competing interests declared
  11. Jing Bo

    Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences, Sanya, China
    Contribution
    Resources
    Competing interests
    No competing interests declared
  12. Li-Sheng He

    Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences, Sanya, China
    Contribution
    Resources
    Competing interests
    No competing interests declared
  13. Qiang Qiu

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Funding acquisition, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  14. Wen Wang

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  15. Shunping He

    1. School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    2. Institute of Deep-Sea Science and Engineering, Chinese Academy of Sciences, Sanya, China
    3. Key Laboratory of Aquatic Biodiversity and Conservation, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China
    Contribution
    Supervision, Funding acquisition, Writing – original draft, Writing – review and editing
    For correspondence
    clad@idsse.ac.cn
    Competing interests
    No competing interests declared
  16. Kun Wang

    School of Ecology and Environment, Northwestern Polytechnical University, Xi'an, China
    Contribution
    Supervision, Funding acquisition, Visualization, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    wk8910@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6059-6529

Funding

National Key Research and Development Program of China

  • Kun Wang

National Natural Science Foundation of China

  • Shunping He

the 1000 Talent Project Shaanxi Province

  • Qiang Qiu

Fundamental Research Funds of Northwestern Polytechnic University

  • Kun Wang

Strategic Priority Research Program of Chinese Academy of Sciences

  • Shunping He

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The project was supported by the National Key R&D Program of China (2022YFC3400300); the National Natural Science Foundation of China (32122021 to KW and 41876179 to SH); the 1000 Talent Project of Shaanxi Province to QQ, KW, and SH; Fundamental Research Funds of Northwestern Polytechnic University; Strategic Priority Research Program of Chinese Academy of Sciences (grant no. XDB42000000), Open Foundation from Marine Sciences in the First-Class Subjects of Zhejiang (No.OFMS011). The authors thank Dr. Yang Zhou, Dr. Yuan Yuan, and Dr. Tao Qin for their advice and discussions during this project.

Senior Editor

  1. Detlef Weigel, Max Planck Institute for Biology Tübingen, Germany

Reviewing Editor

  1. David A Paz-Garcia, Centro de Investigaciones Biológicas del Noroeste (CIBNOR), Mexico

Version history

  1. Sent for peer review: March 9, 2023
  2. Preprint posted: April 17, 2023 (view preprint)
  3. Preprint posted: June 20, 2023 (view preprint)
  4. Preprint posted: December 7, 2023 (view preprint)
  5. Version of Record published: December 22, 2023 (version 1)

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.87198. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2023, Xu, Zhu, Gao et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 451
    Page views
  • 105
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Wenjie Xu
  2. Chenglong Zhu
  3. Xueli Gao
  4. Baosheng Wu
  5. Han Xu
  6. Mingliang Hu
  7. Honghui Zeng
  8. Xiaoni Gan
  9. Chenguang Feng
  10. Jiangmin Zheng
  11. Jing Bo
  12. Li-Sheng He
  13. Qiang Qiu
  14. Wen Wang
  15. Shunping He
  16. Kun Wang
(2023)
Chromosome-level genome assembly of hadal snailfish reveals mechanisms of deep-sea adaptation in vertebrates
eLife 12:RP87198.
https://doi.org/10.7554/eLife.87198.3

Share this article

https://doi.org/10.7554/eLife.87198

Further reading

    1. Computational and Systems Biology
    2. Evolutionary Biology
    Roee Ben Nissan, Eliya Milshtein ... Ron Milo
    Research Article

    Synthetic autotrophy is a promising avenue to sustainable bioproduction from CO2. Here, we use iterative laboratory evolution to generate several distinct autotrophic strains. Utilising this genetic diversity, we identify that just three mutations are sufficient for Escherichia coli to grow autotrophically, when introduced alongside non-native energy (formate dehydrogenase) and carbon-fixing (RuBisCO, phosphoribulokinase, carbonic anhydrase) modules. The mutated genes are involved in glycolysis (pgi), central-carbon regulation (crp), and RNA transcription (rpoB). The pgi mutation reduces the enzyme’s activity, thereby stabilising the carbon-fixing cycle by capping a major branching flux. For the other two mutations, we observe down-regulation of several metabolic pathways and increased expression of native genes associated with the carbon-fixing module (rpiB) and the energy module (fdoGH), as well as an increased ratio of NADH/NAD+ - the cycle’s electron-donor. This study demonstrates the malleability of metabolism and its capacity to switch trophic modes using only a small number of genetic changes and could facilitate transforming other heterotrophic organisms into autotrophs.

    1. Evolutionary Biology
    2. Genetics and Genomics
    Thomas A Sasani, Aaron R Quinlan, Kelley Harris
    Research Article

    Maintaining germline genome integrity is essential and enormously complex. Although many proteins are involved in DNA replication, proofreading, and repair, mutator alleles have largely eluded detection in mammals. DNA replication and repair proteins often recognize sequence motifs or excise lesions at specific nucleotides. Thus, we might expect that the spectrum of de novo mutations – the frequencies of C>T, A>G, etc. – will differ between genomes that harbor either a mutator or wild-type allele. Previously, we used quantitative trait locus mapping to discover candidate mutator alleles in the DNA repair gene Mutyh that increased the C>A germline mutation rate in a family of inbred mice known as the BXDs (Sasani et al., 2022, Ashbrook et al., 2021). In this study we developed a new method to detect alleles associated with mutation spectrum variation and applied it to mutation data from the BXDs. We discovered an additional C>A mutator locus on chromosome 6 that overlaps Ogg1, a DNA glycosylase involved in the same base-excision repair network as Mutyh (David et al., 2007). Its effect depends on the presence of a mutator allele near Mutyh, and BXDs with mutator alleles at both loci have greater numbers of C>A mutations than those with mutator alleles at either locus alone. Our new methods for analyzing mutation spectra reveal evidence of epistasis between germline mutator alleles and may be applicable to mutation data from humans and other model organisms.