1 INTRODUCTION

More and more investigations suggest that life is abundant in the deep sea (Jamieson, 2015; Downey et al., 2018), which upends the conventional view regarding the deep sea as the 'life desert' (Anderson and Rice, 2006). The deep sea is a unique environment due to its high hydrostatic pressure, low temperature and darkness (Anderson and Rice, 2006). High hydrostatic pressure changes the intra or inter-molecular interactions that tend to induce macromolecular irreversible unfolding and aggregation (Balny et al., 1997; Wilton et al., 2008; Lan et al., 2017). Another limitation is the ice-cold temperature that tends to cause nucleic acids to adopt unfavorable structures and hinders enzyme activities (Anderson and Foote, 1975; Feller and Gerday, 2003).

High throughput sequencing can provide a vast amount of information with just a small piece of tissue. At present, multi-omics investigations have greatly advanced our understanding of deep-sea adaptation, such as adapting vision abilities across three deep-sea fishes (Musilova et al., 2019), enhancing levels of organic osmolytes in hadal (the ocean deeper than 6 000 m) snailfish and hadal amphipods (Simonato et al., 2006; Lan et al., 2017; Wang et al., 2019), adjusting the membrane lipid composition in a hadal snailfish (Wang et al., 2019), chemosymbiotic clams (Lan et al., 2019), and mussels (Sun et al., 2017), positively selecting cytoskeleton related genes in bathyal (the ocean at depth from 1 000 to 4 000 m) fish (Lan et al., 2018), positively selecting the DNA repair genes in hadal amphipod (Lan et al., 2017) and bathyal fish (Lan et al., 2018), and adapting the symbiosis system in chemosynthetic ecosystems, including for animals such as limpets (Liu et al., 2020), snails (Sun et al., 2020), clams (Lan et al., 2019), squat lobsters (Cheng et al., 2019), mussels (Sun et al., 2017; Zheng et al., 2017), shrimps (Zhang et al., 2017a), polynoids (Mehr et al., 2015; Zhang et al., 2017b), and crabs (Hui et al., 2017) from hydrothermal vents or cold seeps.

A hadal trench is the deepest point on the Earth. It is ideal for deep-sea adaptation research because of its unique characters, such as the highest hydrostatic pressure on the Earth and a constant low temperature (Jamieson, 2015). However, sample collection is difficult in such a unique environment. At present, only a hadal snailfish genome (Wang et al., 2019) and a hadal amphipod transcriptome (Lan et al., 2017) have been reported. Scarce omics data in hadal trenches limit our knowledge to further understand the molecular mechanism for deep-sea adaptation.

Sea cucumbers (common name for holothurians) are the dominant invertebrates in hadal trenches (Beliaev and Brueggeman, 1989; Jamieson, 2015). Here, we present a transcriptome of a holothurian Paelopatides sp. (Holothuroidea, Synallactida) from the Mariana Trench, which, to our knowledge, is the deepest in-situ fixed macrofauna specimen from the hadal deep. To investigate the molecular mechanism of adaptation to the hadal environment, a positive selection analysis was performed by comparing the transcriptome of the hadal holothurian with its shallow-water relatives.

2 MATERIAL AND METHOD 2.1 Sample collection, ossicle preparation, and RNA extraction and sequencing

During a cruise by DY37-II Dive116 on 16 June 2016, a sea cucumber Paelopatides sp. was captured by the human-occupied vehicle (HOV) Jiaolong on a muddy bottom of the Mariana Trench (141°56.1719'E, 10°57.1693'N) at a depth of ~6 501 m (Fig. 1a), where the seawater temperature was about 1.68 ℃. The body wall of the hadal specimen was crushed into pieces and fixed with RNAlater in situ at the same time as collection. After landing on the deck, the dissected tissues with RNAlater were stored at -80 ℃ until used for sequencing. Two shallow-water holothurians were also collected in November 2016. Stichopus horrens was sampled on a muddy patch of a seagrass bed in Sanya, China. While Apostichopus japonicus was sampled from a muddy coast in Qingdao, China (Fig. 1a). They were captured at a depth of no deeper than 10 m, where the seawater temperature was about 22 ℃ and 15 ℃ respectively. The body wall of these two shallow-water samples were dissected in RNAlater and then stored at -80 ℃.

For scanning electron microscope (SEM) examinations of ossicles, standard protocols of Smirnov et al. (2000) were followed: first, small pieces of the body wall were dissolved in 6% sodium hypochlorite solution, then rinsed five times in ultrapure water, finally the water with ossicles were transferred to aluminum stubs and dried in a draught drying cabinet. The dried ossicles were observed under SEM Phenom ProX (Thermo Fisher Scientific, Waltham, MA, USA).

The total RNA of each specimen (Paelopatides sp., Stichopus horrens, and Apostichopus japonicus) was separately extracted from its inner body wall using a RNeasy Plus Universal Kit (QIAGEN, Hilden, Germany). The quantity and quality of the RNA were examined by agarose gel electrophoresis and a NanoDrop 2000 (Thermo Fisher Scientific) as well as Liang et al. (2020). A pair-end library (150 bp) was first prepared with an NEBNext Ultra RNA Library Prep Kit for Illumina (New England BioLabs Inc, Ipswich, MA, USA) following the manufacturer's recommendations, then purified with an AMPure XP system (Beckman Coulter, Brea, CA, USA), at last assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Inc., Santa Clara, CA, USA). The library was clustered with a TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer's protocol, and sequenced on an Illumina HiSeq 4000 (Illumina). Library preparation and sequencing were completed by NovoGene (Beijing, China).

2.2 Data filtering, assembly, and coding genes prediction and annotation

The raw sequencing reads were evaluated by FASTQC (v0.11.6; Andrews, 2010) and cleaned by fastp (v0.20; Chen et al., 2018) with default parameters. The clean reads were de novo assembled using TRINITY (v2.4.1; Haas et al., 2013) with min_kmer_cov 2 to get the assembled transcriptome. RSEM (Li and Dewey, 2011) was employed to estimate the gene/isoform expression. Only the isoforms with the highest abundance within the same gene were retained (Li and Dewey, 2011). The CD-HIT-EST (v4.6.8; Li and Godzik, 2006) with a threshold of at least 95% similarity was employed for further cleaning (Lan et al., 2017). Finally, only genes longer than 200 bp were regarded as the valuable non-redundant transcripts for subsequent analyses. Benchmarking Universal Single Copy Orthologs (BUSCO v3.0.2b, metazoa_odb9 database; Waterhouse et al., 2018) was employed to evaluate the quality of the transcripts. Coding sequences of the non-redundant transcripts were predicted and translated with TransDecoder (v5.0.2; Haas et al., 2013) with a minimum amino acid length of 50 bp. The non-redundant transcripts were first aligned to the SWISS-PROT database using blastp in BLAST (v2.2.28; Altschul et al., 1997) and the PFAM-A database using the HMMER 3.0 package Hmmscan (Finn et al., 2011) with an E-value < 0.01. Then, the non-redundant transcripts were annotated using DIAMOND (v0.8.22; Buchfink et al., 2015) with an E-value < 1e-5 against the NCBI NR, NT and KOG database. Finally, the domains of the predicted proteins were recovered using InterProScan 5 (v 5.8-51; Jones et al., 2014) with the default parameters.

2.3 Reference species determination

All available genomic and transcriptomic data of the sea cucumbers in Synallactida were screened, including the three species (Paelopatides sp., Stichopus horrens, and Apostichopus japonicus) sequenced in this study, and five reference species from previously reported studies (Table 1): Synallactes chuni (SRR2895367 in NCBI SRA; Janies et al., 2016), Stichopus chloronotus (SRR2846098; Janies et al., 2016), Parastichopus californicus (SRR1695477; Cannon et al., 2014), Parastichopus parvimensis (SRR2484238; Reich et al., 2015) with transcriptome data, and Australostichopus mollis (PRJEB10682; Long et al., 2016) with genome data. The proteins and the coding sequences of Australostichopus mollis genome were adopted from http://ryanlab.whitney.ufl.edu/genomes/Amol/. The transcriptome of Holothuria leucospilota (Holothuroidea, Holothuriida; DRR023763; Chieu et al., 2018) was used as the outgroup in the further phylogenetic analysis. All reference transcriptomes were re-assembled, translated, and assessed with the same pipelines as described above.

2.4 Identification of orthologs and phylogenetic analysis

Single copy orthologs of nine species (Paelopatides sp., Stichopus horrens, Apostichopus japonicus, Parastichopus parvimensis, Synallactes chuni, Stichopus chloronotus, Parastichopus californicus, Australostichopus mollis, and Holothuria leucospilota) were first inferred by orthomcl (v2.0.9; Li et al., 2003), and then were multiple aligned with MAFFT (v7.221; Katoh and Standley, 2013). The poorly aligned positions were trimmed by Gblocks (v0.91b; Talavera and Castresana, 2007). ProtTest (v3.4; Darriba et al., 2011) was used to determine the best model. RaxML (v8.2.11; Stamatakis, 2014) was employed to construct a maximum likelihood (ML) tree with 1 000 bootstraps. The divergence time among species was estimated via r8s (v1.7; Sanderson, 2003). Four calibration nodes based on the Miller et al. (2017) were used as time priors: the most recent common ancestor of Pneumonophora (243-309 million years ago, [Ma]), the most recent common ancestor of Synallactida (117-143 Ma), and Parastichopus californicus - Parastichopus parvimensis (5-7 Ma).

2.5 Positive selection analysis

Transcriptome data combined with a branch site model can effectively identify genome-wide positive selection (Yang and dos Reis, 2011; Yang et al., 2015), and therefore were popular tools for molecular adaptation research (Lan et al., 2017, 2018, 2019). The quality of all the transcripts (Fig. 1b) were screened in BUSCO v3.0.2b (metazoa_odb9 database; Waterhouse et al., 2018), and the species with BUSCO missing values higher than that of the hadal species were not used for the positive selection analysis. Finally, three shallow-water species, Stichopus horrens, Apostichopus japonicus, and Parastichopus parvimensis, as well as Paelopatides sp. were retained. The phylogenetic relationship across these four species was reconstructed as mentioned above. ParaAT (v2.0; Zhang et al., 2012) with a -g parameter was used to align the coding DNA sequences of each ortholog according to their amino acid sequence alignment. A modified branch-site model in the codeml module of the phylogenetic analysis by maximum likelihood (PAML, v4.9h; Yang, 2007) was employed to detect the positively selected genes (PSGs) in deep-sea animals as previously reported (Hui et al., 2017; Lan et al., 2017, 2018, 2019; Sun et al., 2017; Zheng et al., 2017; Cheng et al., 2019). The hadal Paelopatides sp. was designated as the foreground phylogeny; the remaining three shallow-water holothurians were designated as the background phylogeny. An alternative branch site model (Model=2, NSsites=2, and fix_omega=0) and a neutral branch site model (Model=2, NSsites=2, fix_ omega=1 and omega=1) were configured. P-values were first computed by a chi-squared test (Zhang et al., 2005) and then corrected by a multiple testing correction (Benjamini and Hochberg, 1995). Genes with Bayesian Empirical Bayes (BEB) sites exceeding 90% and corrected P-values lower than 0.05 were considered as PSGs (Lan et al., 2017). In addition to the annotation information of the seven databases mentioned above, the PSGs were also imported into InterProScan 5 (Jones et al., 2014) to obtain more annotation information.

As to the identified PSGs in the hadal Paelopatides sp., we also compared their expression level with their two shallow-water relatives (Stichopus horrens and Apostichopus japonicus). The transcripts per million (TPM) was used as the measure of expression, as it has been recommended to be more comparable across samples than fragments per kilobase of transcript per million mapped reads (FPKM; Li et al., 2010). The heatmap was generated by pheatmap packages (https://cran.r-project.org/web/packages/ pheatmap) in R (v3.2; The R Core Team, 2016).

3 RESULT

In this study, 129 376 486 raw paired end reads (150 bp) of Paelopatides sp. were generated. After trimming, 125 583 692 clean reads (97.07%) were retained and used for the de novo assembly (Table 1). After removing redundant isoforms, 127 210 non-redundant contigs were produced. The non-redundant contigs ranged from 201 to 15 346 bp with a total size of 115 979 417 bp. The contig N50 of the assembled transcripts was 1 260 bp (Supplementary Fig.S1a & b; Table 1). These assembled transcripts hit 935 (95.60%) of the single copy orthologs in the BUSCO metazoan database (Fig. 1b).

In addition, 62 230 protein-coding sequences were obtained from the non-redundant transcripts of Paelopatides sp. (Table 1). The annotation rate of the hadal Paelopatides sp. was slightly lower when compared to that of the two shallow-water holothurian species (Supplementary Fig.S1c; Supplementary Table S1). Interestingly, in the annotation of GO, KEGG, and KOG databases, the annotation percentage of genes related to some pathways was consistently higher in the hadal Paelopatides sp. when compared with their two shallow-water holothurians, including replication and repair; translation; folding, sorting and degradation; cell cycle control, cell division, and chromosome partitioning; chromatin structure and dynamics; and energy production and conversion. Furthermore, the percentage of genes related to rhythmic process and extracellular transport was consistently lower in the hadal Paelopatides sp. (Supplementary Figs.S2-S4; Supplementary Tables S2-S4).

In this study, 706 single-copy orthologs were identified across Paelopatides sp., Stichopus horrens, Apostichopus japonicus, Parastichopus parvimensis, Synallactes chuni, Stichopus chloronotus, Parastichopus californicus, Australostichopus mollis, and Holothuria leucospilota. An ML tree (Fig. 1c) was constructed based on these single copy orthologs with 39 930 amino acids. Synallactes chuni showed the farthest phylogenetic relationship to Paelopatides sp. (Fig. 1c). Stichopus chloronotus, Parastichopus californicus, and Australostichopus mollis had a higher missing number of the single copy orthologs than Paelopatides sp. (Fig. 1b). We did not use these four species in the following positive selection analysis to increase the numbers of single copy orthologs. Moreover, 24 466 paired-pared and shared orthologs were identified across Paelopatides sp., Stichopus horrens, Apostichopus japonicus, and Parastichopus parvimensis. As shown in Fig. 1e & f, 9 341 gene families were shared across all four species. Among them, 6 934 were the single copy orthologs (Fig. 1e & f). In these single copy orthologs, 113 genes in Paelopatides sp. were detected to be positively selected (Supplementary Table S5). According to KEGG ortholog (KO, Supplementary Tables S6) and enrichment analysis, among these 113 PSGs, at least 10 PSGs (such as microphthalmia-associated transcription factor [MITF], transcription factor 7-like 2 [TCF7L2], tartrate-resistant acid phosphatase type 5 [ACP5], lysosomal acid phosphatase [ACP2], disks large protein 1 [DLG1], putative fibrillin-2-like isoform X1 [FBN2], presenilin 1 [PSEN1], large subunit ribosomal protein L38 [RPL38], collagen type II alpha and type V/XI/XXIV/ XXVII alpha [COL2A & COL5AS]) are related to the distinct phenotypic characteristics of Paelopatides sp. And at least 24 PSGs (ETS-related transcription factor Elf-2-like [ELF2], PERQ amino acid-rich with GYF domain-containing protein [GIGYF], cytochrome c oxidase subunit 7a [COX7A], SWI/ SNF related-matrix-associated actin-dependent regulator of chromatin subfamily C [SMARCC], SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily E [member 1 [SMARCE1], type I thyroxine 5'-deiodinase [DIO1], translation factor GUF1 [mitochondrial [GUF1], cell cycle checkpoint control protein RAD9A [RAD9A], replication factor A3 [RPA3], DNA-directed RNA polymerases I/II/III subunit RPABC1 [POLR2E], putative TAR DNA-binding protein 43 isoform X2 [TARDBP], ribonucleoside-diphosphate reductase subunit M1 [RRM1], putative serine/threonine-protein kinase SMG1 [SMG1], transcriptional regulator ATRX [ATRX], alkylated DNA repair protein alkB homolog 6 [ALKBH6], PLAC8 motif-containing protein [PLAC8], putative 46 kDa FK506-binding nuclear protein-like [FKBP], peptide-methionine (S)-S-oxide reductase [MSRA], prefoldin alpha subunit [PFDN5], BCL2-associated athanogene 2 [BAG2], abhydrolase domain-containing protein 6 [ABHD6], beta-1, 4-galactosyltransferase 6 [B4GALT6], cholesterol 24-hydroxylase [CYP46A1], and AT-rich interactive domain-containing protein 2 [ARID2]) may contribute to the deep-sea adaptation, especially for the cold adaptation and cell homeostasis maintenance that mainly include repairing the damaged DNA induced by high hydrostatic pressure, ice-cold temperature, and some heavy metal such as cadmium (Cd), and sustaining structure of protein, membrane and chromosome (Table 2).

4 DISCUSSION

Investigations into deep-sea adaptation have drawn considerable attention for decades, and comparative omics studies have been proven a powerful way to reveal adaptive mechanism (Foote et al., 2015; Hui et al., 2017; Zheng et al., 2017; Liu et al., 2020). Genome-wide positive selection analysis is an effective tool for adaptation research as it can connect ecological change with molecular change (Yang et al., 2015; Lan et al., 2017, 2018, 2019; Wang et al., 2019). In this study, a positive selection analysis was performed across four closely related holothurians. Compared with its three shallow-water relatives, 113 PSGs (Supplementary Table S5) were identified in the hadal holothurian Paelopatides sp. Combined with environmental factors, some PSGs are supposedly related to the distinct phenotypic characteristics of Paelopatides sp., and some supposedly contributing to cold adaptation and cell homeostasis maintenance (Table 2, Supplementary Table S6) in the hadal deep surroundings.

4.1 Translucent white body

Many animals that dwell in dark surroundings are translucent white, such as some cave animals (McGaugh et al., 2014), hadal snailfish (Wang et al., 2019), and some hadal holothurians (Jamieson et al., 2011; Martinez et al., 2019). Our hadal holothurian Paelopatides sp. is also translucent white (Fig. 1c). Compared with its colored shallow-water relatives, the genes MITF and TCF7L2 are positively selected in Paelopatides sp. It is well known that mutations of MITF may cause the white phenotype in many organisms (Schmutz et al., 2009; Hauswirth et al., 2012). In a translucent white hadal snailfish, MITFA is also thought responsible for its unpigmented skin (Wang et al., 2019). Consistently, previous investigations have also shown less melanocytes and lower expression levels of MITF in the transparent white holothurians (Xing et al., 2018). Regarding the MITF in Paelopatides sp., positively selected amino acid sites mainly distribute in and around the active site in the transactivation region of MITF (Fig. 2a). These mutations might influence the activity of MITF and affect the expression profiles of its target genes (Fig. 2e-g). Furthermore, the positively selected amino acids of MITF tend to adopt longer branched amino acids (Fig. 2b). Higher hydrostatic pressure tends to reduce biomacromolecules into smaller volumes and tighter structures (Mozhaev et al., 1996). The PSGs in deep-sea creature substitute with longer branched amino acids that might be involved in sustaining protein structure and coping with high hydrostatic pressure. Additionally, the gene TCF7L2 has been reported to show a lower expression level in non-melanoma skin cells (Dehcheshmeh et al., 2018). These two genes may not only participate in the regulation of melanocytes to produce melanin that gives skin its coloration (Moore, 1995), but also help to protect cells from DNA damage induced by UV (Nguyen and Fisher, 2019). The hadal zone is far from sunlight and lack of UV radiation, there is no need to keep melanogenesis as active as the shallow-water species. Here, two melanogenesis PSGs (MITF and TCF7L2) were also expressed at a lower level in hadal Paelopatides sp. compared to their shallow-water relatives (Fig. 3), which might be an economic strategy to live in darkness. Therefore, it suggests that these two PSGs might be involved in darkness adaptation besides may contribute to the pigment loss and translucent white skin.

4.2 Ossicle degeneration

The endoskeleton ossicles dispersed in the dermis of the body walls could protect sea cucumbers. In addition, the characteristics of ossicles also serves as a taxonomic index in Holothuroidea (Liao, 1997; Miller et al., 2017). Based on our observations (Fig. 1c & d), the body of the hadal holothurian Paelopatides sp. was smooth in low density, small, and concise ossicles in its dermis, and we found no plate ossicle. In contrast, the shallow-water holothurians were covered with rough papillae where high density, large and complex ossicles distributed. Many large plates were found inside their skin. Ossification degeneration phenomena had also been found in some other deep-sea holothurians (Lee Hufford, 1968; Martinez et al., 2019), and ossicle degeneration has been reported in hadal snailfish (Wang et al., 2019). In Paelopatides sp., eight PSGs (ACP5, ACP2, COL2A, COL5AS, DLG1, FBN2, PSEN1, and RPL38) were enriched in the biological process of skeletal system development (Table 2, Supplementary Table S6), and previous investigations had shown that lacking MITF (Lu et al., 2010, 2014) and ACP5 (Hayman et al., 1996) disrupt the ossification. In addition, nine skeleton related PSGs (MITF, ACP5, ACP2, COL2A, COL5AS, DLG1, FBN2, PSEN1 and RPL38) were also consistently expressed at a lower level in the hadal Paelopatides sp. (Fig. 3), which suggests that these PSGs may be the molecular reminders for ossification degeneration of hadal Paelopatides sp.

4.3 Cold adaptation

Temperature is an important environmental driver for species distribution (Danovaro et al., 2004). The temperature of the habitat where the Paelopatides sp. was collected was about 1.68 ℃ (record from HOV Jiaolong). Sea cucumbers are poikilothermal benthos, which cannot regulate their body temperature by themselves (Yang et al., 2006). That means the biochemical systems in the hadal holothurians have to operate in such frozen temperature. In the ice-cold temperature, some enzymes lose their activity, nucleic acids tend to adopt unfavorable structure, and some proteins are out of their normal function when involved in genetic information expression (Feller and Gerday, 2003; Lan et al., 2017, 2018). Moreover, chemical reaction rate have been proofed to be decreased in the low temperature (Carney, 2005; Jamieson, 2015). Even living in such frozen temperature, Paelopatides sp. still adapt and thrive in hadal trenches (observation from Jiaolong; Jamieson, 2015; Martinez et al., 2019). Therefore, how Paelopatides sp. adapt to their ice-cold surroundings is a fantastic point to be further explored.

It is reported that some cold shock domain containing proteins can be involved in cold adaptation by unwinding the unfavorable secondary structures and by facilitating the genetic processes (Thieringer et al., 1998; Lim et al., 2000; Lan et al., 2018). Cold shock proteins are synthesized to overcome the deleterious effects of cold temperatures (Phadtare et al., 1999; Kang et al., 2015; Lan et al., 2017). In this study, ETS domain (TCF7L2, ELF2) and helicase conserved C-terminal domain (ATRX) are positively selected in hadal Paelopatides sp., which are important cold shock domains (Gualerzi et al., 2003; Fry et al., 2018). The homologs of TCF7L2 and ELF2, TCFIIA and ELF1, have been reported to have experienced significant expansion in the hadal amphipod Hirondellea gigas (Lan et al., 2017) and the Antarctic amphipod Gondogeneia antarctica (Kang et al., 2015), and they are suggested to be involved in cold adaptation (Tang et al., 2015; Nie et al., 2019). The Helicase conserved C-terminal domain containing proteins are commonly positively selected across hadal Hirondellea gigas (EIF4G; Lan et al., 2017), bathyal Aldrovandia affinis (EIF4B; Lan et al., 2018), hadal Pseudoliparis swirei (IFIH1; Wang et al., 2019) and the hadal Paelopatides sp. (ATRX). Meanwhile, the gene GIGYF can regulate mRNA surveillance for cold adaptation by selective interaction with the helicase conserved C-terminal domain containing protein EIF4F (Peter et al., 2017), which was positively selected in Paelopatides sp.

In addition, three thermogenesis related genes (COX7A, SMARCC, and SMARCE1) are positively selected in the hadal Paelopatides sp. (Table 2, Supplementary Table S6). COX7A is an adipocyte marker gene, while adipocytes play important roles in adaptive thermogenesis (Fisher et al., 2012; Maurer et al., 2015). SMARCC and SMARCE1 are molecular determinants in response to temperature (Jin, 2018; Leal et al., 2018). The core thermogenesis PSGs (COX7A) had a consistent higher expression in the hadal Paelopatides sp. than its shallow-water relatives (Fig. 3). Moreover, the gene DIO1 and mitochondrial translation factor GUF1 are also positively selected in Paelopatides sp. DIO1 is a deiodinase for thyroid hormone that plays a critical role in cold adaptation in invertebrate (Heyland and Moroz, 2005; Bianco and Kim, 2006; Gereben et al., 2008; Panicker et al., 2008). Mutation of GUF1 would be more sensitive to temperature change and diminished rates of protein synthesis in low temperatures (Bauerschmitt et al., 2008). Therefore, positively selected GUF1 may promote mitochondrial protein synthesis under hadal cold adversity.

4.4 Cell homeostasis maintenance

DNA is one of the most important biomacromolecules carrying genetic information for most organisms. Deep-sea animals are susceptible to DNA damage from high hydrostatic pressure (Dixon et al., 2004; Lan et al., 2017) and ice-cold temperature (Anderson and Foote, 1975; Gualerzi et al., 2003). To ensure genetic fidelity, organisms in the deep sea have evolved a variety of mechanisms to efficiently repair DNA lesions (Dixon et al., 2004). At least eight out of all the PSGs in the hadal Paelopatides sp. are involved in DNA repair including RAD9A, RPA3, POLR2E, TARDBP, RRM1, SMG1, ATRX, and ALKBH6 (Table 2). RPA3 is involved in repairing genetic materials in the replication process (Zou et al., 2006). POL2RE (Bhagavan and Ha, 2015), RAD9A (Lieberman et al., 2017) and ATRX (De La Fuente et al., 2011; Han et al., 2018) are involved in the transcription-coupled repair. TARDBP assembles DNA repair endonucleases to seal DNA breaks (Mitra et al., 2019). RRM1 is necessary for DNA repair especially for the uracil-DNA repair (Ladner, 2001). SMG1 is one of the stress responsive PIKK family members, which serves as the central regulator of DNA damage responses (DDR; Lloyd et al., 2018). In this study, SMG1 had a consistent higher expression in the hadal Paelopatides sp. than its shallow-water relatives (Fig. 3). In addition, ALKBH6 may mediate alkylation damage repair by similarity (Fedeles et al., 2015). These PSGs related to DNA repair may play an important role in maintaining the fidelity of genetic materials in deep-sea surroundings.

In addition to high hydrostatic pressure and low temperature, high concentrations of heavy metals may also cause DNA damages (Lin et al., 2007; Abbà et al., 2011). The upper Mariana Trench (6 500- 7 626 m) is reported to be exposed to the higher concentrations of cadmium (Welty et al., 2018). Cadmium is genotoxic and can induce DNA damage (Lin et al., 2007), while a previous experiment proved that PLAC8 motif-containing protein (OmFCR) and cell cycle checkpoint control protein (RAD9) are involved in the pathway for the response to cadmium exposure (Abbà et al., 2011). In the hadal Paelopatides sp., PLAC8 and RAD9A are positively selected, which might be involved in the adaptation to the high concentration of heavy-metal surroundings in the upper Mariana Trench.

Moreover, many of the rest of the PSGs (Supplementary Table S6) in Paelopatides sp. are enriched in protein folding, lipid metabolism and chromosome regulated, including FKBP (Kuzuhara and Horikoshi, 2004), MSRA (Stadtman et al., 2005), PFDN5, and BAG2 (Supplementary Table S6), which are dedicated to protein folding. ABHD6, B4GALT6, and CYP46A1 are involved in lipid metabolism (Supplementary Table S6). H3 (Liu et al., 2012) and ARID2 (Zhao et al., 2011) are involved in the chromatin-remodeling complex. High hydrostatic pressure may also result in cellular and macromolecular structural changes including causing protein denaturation (Balny et al., 1997; Somero, 2003; Simonato et al., 2006), reducing the fluidity of cell membranes (Yano et al., 1998; Simonato et al., 2006), and changing the characters of chromosome (Pease, 1946). These structural related PSGs might contribute to the adaptation to the hadal environment with extra-high hydrostatic pressure.

5 CONCLUSION

In this study, the deep-sea adaptation of a hadal holothurian Paelopatides sp. was investigated. A high-quality transcript and a phylogenetic tree based on omics data were obtained. Based on the positive selection analysis, 113 PSGs were identified in Paelopatides sp. Some PSGs (such as MITF) may contribute to its translucent white body and degenerated ossicles. Some are thermogenesis related genes (COX7A, SMARCC, SMARCE1, DIO1, and GUF1) or cold shock genes (TCF7L2, ELF2, ATRX, and GIGYF) that may be dedicated to cold adaptation. Many of the other PSGs are mainly involved in cell homeostasis maintenance: at least nine PSGs (RAD9A, RPA3, POLR2E, TARDBP, RRM1, SMG1, ATRX, ALKBH6, and PLAC8) may repair the DNA damage induced by the high hydrostatic pressure and ice-cold temperature. In addition to the DNA repair, PLAC8 may facilitate the adaptation to the high concentration of cadmium in the upper Mariana Trench. The valuable in-situ fixed sample from the Mariana Trench, the high-quality assembly, and the positive selection analysis on the hadal holothurian in this study shed light on deep-sea adaptation research.

6 DATA AVAILABILITY STATEMENT

All transcriptome data generated from this study were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive database (PRJNA613016, SRA Accession Nos.: SRR11535194-SRR11535196).

7 ACKNOWLEDGMENT

We would like to thank Huaining GAO, Zhiqiang WANG, Yanying YE, Yang YANG, Baosheng LI for their helps in in-situ sampler preparation. The authors also want to thank the captains and crews of the R/V Xiangyanghong 09 and the pilots of HOV Jiaolong for their technical support.

8 ABBREVIATION

PSG(s) abbreviate for positively selected gene(s), and full description of abbreviations for each PSG can be found in Table 2 and Supplementary Table S5.

Electronic supplementary material

Supplementary material (Supplementary Figs.S1-S4 and Tables S1-S6) is available in the online version of this article at https://doi.org/10.1007/s00343-020-0241-0.