Insights into deep‐sea adaptations and host–symbiont interactions: A comparative transcriptome study on Bathymodiolus mussels and their coastal relatives

Mussels (Bivalve: Mytilidae) have adapted to various habitats, from fresh water to the deep sea. To understand their adaptive characteristics in different habitats, particularly in the bathymodiolin mussels in deep‐sea chemosynthetic ecosystems, we conducted a comparative transcriptomic analysis between deep‐sea bathymodiolin mussels and their shallow‐water relatives. A number of gene families related to stress responses were shared across all mussels, without specific or significantly expanded families in deep‐sea species, indicating that all mussels are capable of adapting to diverse harsh environments, but that different members of the same gene family may be preferentially utilized by different species. One of the most extraordinary trait of bathymodiolin mussels is their endosymbiosis. Lineage‐specific and positively selected TLRs and highly expressed C1QDC proteins were identified in the gills of the bathymodiolins, suggesting their possible functions in symbiont recognition. However, pattern recognition receptors of the bathymodiolins were globally reduced, facilitating the invasion and maintenance of the symbionts obtained by either endocytosis or phagocytosis. Additionally, various transporters were positively selected or more highly expressed in the deep‐sea mussels, indicating a means by which necessary materials could be provided for the symbionts. Key genes supporting lysosomal activity were also positively selected or more highly expressed in the deep‐sea mussels, suggesting that nutrition fixed by the symbionts can be absorbed in a “farming” way wherein the symbionts are digested by lysosomes. Regulation of key physiological processes including lysosome activity, apoptosis and immune reactions is needed to maintain a stable host–symbiont relationship, but the mechanisms are still unclear.

Since the hydrothermal vent community was discovered in 1977, which is powered by chemosynthetic bacteria, mussels of subfamily Bathymodiolinae have been found to be one of the dominant species in a number of deep-sea habitats including hydrothermal vents and cold seeps, and to be also common in whale falls and sunken woods (Duperron et al., 2011;Johnson, Won, Harvey, & Vrijenhoek, 2013;Kiel & Hansen, 2015;Lorion et al., 2013;Tyler & Young, 2003). Bathymodiolin mussels represent a diverse macrofaunal group that can survive in many harsh environments. Previous phylogenetic work using both fossil and molecular evidence showed that vent and seep communities were colonized recently by modern mytilid taxa, with whale falls and sunken wood possibly serving as stepstones Lorion et al., 2013;Thubaut, Puillandre, Faure, Cruaud, & Samadi, 2013). After ancestral divergence, the bathymodiolins evolved a range of properties to survive in the new habitats. However, the means by which this colonization took place, including how bathymodiolin mussels have adapted to deep-sea chemosynthetic environments over the long term, and the genetic bases for these adaptations, are still under discussion. Answers to these questions may shed light on the adaptive mechanisms of mussels in diverse environments.
Despite numerous investigations, the nature of the relationship between bathymodiolin mussels and their endosymbionts remains controversial, especially with respect to the host-symbiont recognition and the interactive processes in the onset of endosymbiosis. In cnidarians, the key model organism for studying the host-symbiont relationships, the recognition and the dynamic control of symbionts begins with initial cell contact and persists through the entire symbiotic processes. Such winnowing not only determines the species of the endosymbionts, but also the dynamics and fate of the symbiosis (Dunn & Weis, 2009). Evidence for the activation of tolerogenic host immune pathways has been found in the cnidarian during the dinoflagellate infection, which may function to aid the establishment of the symbiosis (Detournay, Schnitzler, Poole, & Weis, 2012). Meanwhile, in insects, the host immune systems is either simplified or suppressed by the symbionts to facilitate infection by and maintenance of beneficial bacteria (Ratzka, Gross, & Feldhaar, 2012;Richards et al., 2010). Mussels have complete gene sets related to the innate immune responses (Gerdol & Venier, 2015), as also found in Bathymodiolus platifrons (Wong et al., 2015) and B. azoricus (Bettencourt et al., 2010b). However, how the innate immune system of the bathymodiolin interacts with the symbionts is still unclear.
Transcriptomic analysis is a powerful approach to reveal mRNAlevel responses of mussels dwelling in diverse habitats (Wang, Liu, & Wu, 2016;Woo et al., 2011;Xu et al., 2016), and also has been widely applied to study symbiotic relationships (Braendle et al., 2003;Heddi et al., 2005;Nakabachi et al., 2005). However, only a few transcriptomic studies have examined host-symbiont interactions in deep-sea chemosynthetic mussels, and most of them are focused on the physiology of the host (Bettencourt et al., 2010b;Wong et al., 2015). Linking variation in transcriptomic analysis in a comparative framework can offer a nuanced insight into the mechanistic basis of adaption (Pan et al., 2015;Velotta et al., 2016;Yang, Qi, & Fu, 2014). To define the general adaptive genetic traits of the bathymodiolins, we compared B. platifrons mussels from a cold seep, which host methanotrophic symbionts (Barry et al., 2002), and B. manusensis from a hydrothermal vent, which host thiotrophic symbionts (Lorion et al., 2013), to their shallow-water relatives. These were represented by Modiolus kurilensis from cold-water environments off northeast China and the intertidal mussel Perna viridis from warm water. M. kurilensis has been found to be a close relative of the bathymodiolins based on mitochondrial genomic data, and will serve as good candidate for comparison with Bathymodiolus spp. (under review). We aim to (i) provide an outline of the genetic basis for mussel adaptations to diverse habitats, focusing on deep-sea mussels; (ii) elucidate genetic features that permit the linkage between mussels and their endosymbionts, particularly with regard to immune recognition; and (iii) test hypotheses on how deep-sea mussels obtain nutrients from a transcriptomic framework. This work will contribute to our understanding of the mechanism that allows mussels survive in harsh environments.

| Mussel samples
Bathymodiolus platifrons specimens were obtained from a cold seep at a depth of 1,100 m in the South China Sea (referred to as SS).
B. manusensis specimens were obtained from a hydrothermal vent at Manus Basin, at a depth of 1,520 m in the western Pacific (referred to hereafter as MS). Both were collected using the ROV Faxian, operated from the RV Kexue, during 2014 and 2015, respectively. M. kurilensis was sampled using SCUBA from Xiaoping Island at Dalian (DL) (Figure 1). P. viridis specimens were previously recorded by Leung et al. (2014) and were not collected afresh for this investigation. Gills were dissected out as quickly as possible after collection and flash-frozen in liquid nitrogen and then stored at À80°C for further study.

| Preparation and sequencing
Three specimens were sequenced for each species. Total RNA was separately extracted using E.Z.N.A mollusc RNA Kit (OMEGA, China) following the manufacturer's instructions. Quality control was evaluated using Nanodrop TM spectrophotometer and gel electrophoresis.
RNA concentration was measured using a Qubit â RNA Assay Kit in Qubit â 2.0 Fluorometer (Life Technologies, CA, USA). mRNA capture using poly-T oligo-attached magnetic beads, cDNA synthesis, pairend library preparation, and Illumina sequencing were completed by NovoGene (Beijing, China). Quality controls were carried out using in-house Perl scripts to remove contaminated adaptors and low-quality sequences. Six thousand reads from each Illumina sample were randomly selected to evaluate the rate of possible contamination by BLAST search (e-value = 1e À 1) against the NCBI NO-REDUNDANT (NR) database. RNA-seq data of P. viridis were retrieved from the National Center for Biotechnology Information (NCBI, Accession nos SRX643383 and SRX643382).

| Data mining and gene functional annotation
Transcriptome assembly was accomplished using TRINITY (2.1.1) (Grabherr et al., 2011) with default parameters, with the exception of min_kmer_cov, which was set to 1. Proteins were predicted by FRAMEDP (Gouzy, Carrere, & Schiex, 2009) using a refined database containing proteins sourced from mollusc species whose genomes are publicly available, as well as the latest Swiss-Prot proteins. For transcripts that could not be orientated by FRAMEDP, ESTSCAN (Iseli, Jongeneel, & Bucher, 1999) was performed. The transcripts were first F I G U R E 1 Sampling sites for deep-sea Bathymodiolus mussels and their coastal relative Modiolus kurilensis. The star indicates a cold seep; the dot indicates a hydrothermal vent; and the square indicates a shallow-water site. The map was created with GEOMAPAPP (version 3.6.2) [Colour figure can be viewed at wileyonlinelibrary.com] ZHENG ET AL. | 5135 annotated using BLASTP (2.2.31+; e-value <1e À 5) against the NCBI NR database, KOG/COG, KO and SWISS-PROT databases. In this step, possible contaminated contigs from human sources were ruled out based on species distribution and sequence identity (0.95) of the best hit of BLASTP results against the NR database. Domains of predicted proteins were recovered by  with the default parameter. The NR annotation and INTERPROSCAN results were used to obtain the GO term assignments of the transcripts with the BLAST2GO (3.1.3; e-value <1e À 6) program. CEGMA estimation (Cai et al., 2015) was applied to assess the completeness of the assembly.

| Orthology determination and positive selection test
To reconstruct a phylogeny of mussels, orthologous groups were predicted from the BLASTP result with ORTHOMCL version 2.0.3 (Li, Stoeckert, & Roos, 2003) using default settings for two data sets. were downloaded from NCBI. These three transcriptome sequences were assembled with TRINITY as previously described. One-to-one orthologs were determined using an original Perl script. To determine the lineage-specific, positively selected codons, we used an optimized branch-site model on the basis of the best maximum-likelihood (ML) tree. A likelihood ratio test was used to compare a model that allowed sites under positive selection in the foreground branch with the null model, in which the omega value was fixed to 1. pvalues were computed based on the chi-square statistic adjusted by the FDR method; genes with adjusted p values <.05 were treated as candidates for positive selection, which are termed as positively selected genes (PSGs) in the following analysis.
2.5 | Quality assessment and differential expression analysis Clean reads were mapped to the assembly using BOWTIE2 (Langmead & Salzberg, 2012). The differentially expressed genes (DEGs) were screened using DESEQ2 (Love, Huber, & Anders, 2014; p adj < .05). The expression levels of gills in each mussel were presented as FPKM, which is generated by DESEQ2.

| Phylogenetic analyses and statistics
Proteins were aligned using MSAPROBS (Liu & Schmidt, 2014). Corresponding CDSs (coding sequence) were aligned based on the peptide alignments using an original Perl script. Phylogenetic inferences were conducted on the protein alignments under maximum likelihood (ML) using RAXML version 8.2 (Stamatakis, 2006). Enrichment analysis (p adj < .05) on the selected data sets (differentially expressed genes and positively selected genes) was performed using ERICHPIPELINE developed by NOVOGENE (Beijing, China).
Heatmaps and other diagrams were made in R. The degree of enrichment was indicated as a Rich factor.

| RESULTS
3.1 | Illumina sequencing, de novo assembly and sequence validation RNA-seq sequencing was conducted on gill samples from three randomly collected individuals (biological repeats) for each species by Novogene Biotechnology, Beijing. More than 10 Gb clean data were generated for each sample after trimming and quality control. After assembly, mean read length of contigs was longer than 572 bp for all the samples. We also assembled the RNA-seq data of P. viridis under NCBI Accession nos SRX643383 and SRX643382. For B. manusensis, B. platifrons, M. kurilensis and P. viridis, 259,949, 137,591, 136,955 and 120,750 unigenes were obtained, respectively, via de novo assembly of the merged clean data for each species. The completeness of the transcriptome was assessed by CEGMA estimation and by the recovery rate of mitochondrial genome coding protein genes (mtPCGs). Complete sequences were obtained for all the mtPCGs. Of 458 highly conserved eukaryotic genes from the CEGMA set, more than 96% were identified in the assembled transcriptomes, with the exception of B. manusensis (which was still higher than 83%). This suggests a high data set quality (Table 1). *"# Prots" indicates the number of predicted proteins in the set of 248 core eukaryotic genes (CEGs). "Complete" refers to those predicted proteins in the set of 248 CEGs that when aligned to the HMM for the KOG for that protein family give an alignment length that is 70% of the protein length. If a protein is not complete, but still exceeds a precomputed minimum alignment score, then the protein is called "partial." The assessment was carried out using the official pipeline under default parameters.

| Gene annotation
respectively. Similar to other de novo annotations, less than half (22.25%-30.24%) of the unigenes were annotated in at least one database. In NR annotation species distribution results, Crassostrea gigas was ranked highest, followed by Lottia gigantea and Aplysia californica ( Fig. S1). There were 39,687, 27,260, 30,590 and 24,594 transcripts assigned to three major Gene Ontology (GO, 1st level): "cellular components," "molecular functions" and "biological processes" (Fig. S2), showing no obvious differences between species.
The category distribution of the functional annotation in these databases of the four mussels is given in Fig. S3. Pattern recognition receptors (PRRs) were annotated with the combining evidence from BLAST (NR, KEGG and SWISS-PROT databases) searches and INTERPROSCAN domain annotations. Only transcripts with the presence of both conserved domains (according to Gerdol & Venier, 2015) and the corresponding annotations were counted, and the numbers of PRRs are shown in Sheet S1. The representative PRRs are identified in all of the four mussels, large numbers of C1qDC proteins, C-type lectins (CTLs), fibrinogen-related proteins (FREPs) and Toll-like receptors (TLRs) have been found.

| Identification of orthologs
A total of 226,303 putative orthologs were identified in four mussels using ORTHOMCL (Figure 2). There were 10,730 orthologs shared by all four species, constituting the mussel core gene set. The core genes participate in a wide range of fundamental activities and were distributed in almost every GO function; 6,064 orthologs (deep-set) were only found in both deep-sea mussels. As shown in pfam enrichment analyses, many members of the deep-set possessed domains homologous to transposase, for example transposase DDE domain, DDE superfamily endonuclease, helix-turn-helix of resolvase and others. Meanwhile, 1,714 orthologs were specifically identified from two shallow-water mussels (shallow-set) with domains enriched in autophagy protein Apg6, UV radiation resistance protein, autophagy-related subunit 14 and others (Figure 3a,b). According to the unigene numbers of each gene family, genes that may accelerate genome reorganization such as RNA-directed DNA polymerase from mobile element jockey-like, piggyBac transposable element-derived protein 4-like, and transposable element Tcb2 transposase, transposase Int-Tn were expanded in the deep-sea mussels. Meanwhile, gene families involved in immune recognition such as C1q domaincontaining proteins were reduced in deep-sea mussels (Table 2).

| Phylogenetic analysis
We reconstructed a phylogeny of mussels, as well as other Lophotrochozoa, under a RNA-seq framework using 656 single-copy orthologs with a total length of 80,043 aa. All mussels formed a monophyletic clade, with Crassostrea gigas and Pinctada fucata as a sister group. The Bathymodiolus mussels grouped together as a monophyletic clade, with M. kurilensis as a sister group, while the other mytiliform mussels were grouped in another cluster. Bathymodiolus was not located at a basal position in the Mytilidae clade ( Figure 4). Our results are consistent with the hypothesis that the deep-sea bathymodiolin mussels are recent colonizer to vents and seeps Lorion et al., 2013). The close affinity of the bathymodiolin and modiolin mussels will lead to a more meaningful comparative result for adaptive evolution.

| Positively selected genes (PSGs)
After their divergence from shallow-water species, Bathymodiolus spp. resided in a habitat considerably different and evolved a specific series of physiological features. To detect genes that might directly contribute to the adaption to a deep-sea chemosynthetic lifestyle, 9,652 1:1 orthologs were used to screen the positively selected genes (PSGs). In total, we identified 1,314 PSGs from the two deepsea mussels (D-PSGs) and 228 PSGs from the two shallow-water mussels (S-PSGs), demonstrating that mussels in the deep sea have more genes under positive selection.
Pfam analysis showed that domains related to protein kinase With regard to immune response, there were differences between mussel groups related to extracellular matrix, immune     In addition, a number of genes involved in transport, calcium homeostasis (mucolipin, polycystin, calbindin) and cytoskeleton reorganization (dynein, myosin and some other actin and tubulin related proteins) were both found in D-PSGs and S-PSGs. Among these groups, many more transporters types and quantities were found in D-PSGs, including diverse iron transporters (ferroportin1, coppertransporting ATPase 1, zinc transporter ZIP1, ferrous iron transport protein B, taurine transporter, ammonium transporter Rh type B-A), and many other nutrient, vitamin and coenzyme transporters in the ATP-binding cassette transporter superfamily (ABCs), the permease family and solute carrier family (SLCs) (Table 3; Figure 5).
Furthermore, many lineage-specific PSGs were only identified in D-PSGs, which might be related to the adaptive evolution of the bathymodiolin mussels to the deeper and chemosynthetic environment over the course of their colonization. These PSGs include mucin-5AC-like, caveolin, early endosome antigen 1, vacuolar protein sorting-associated protein 18 (VPS18), vesicle-associated membrane protein 4 (VAMP4), V-type proton ATPases (ATPeV), lysosome-associated membrane glycoprotein (LAMP) and others ( Figure 5).
3.6 | Differences in expression pattern between deep-sea and shallow-sea mussels Among all the significant differentially expressed genes (DEGs), 304 DEGs were expressed at a higher level in deep-sea mussels (D-DEGs), while the shallow-water mussels had 436 upregulated DEGs (S-DEGs). The functional enrichment results indicated that DEGs were mainly involved in the maintenance of cellular homeostasis, transport and immune reactions. However, differences were obvious between the deep-sea and shallow-water mussels. Such differences between the two habitats may reflect innate responses to different environmental conditions, but may also be caused by stress during the sampling process, during which almost all abiotic conditions change dramatically for the deep-sea mussels. To elucidate this, we compared RNA expression in samples used here and mussels fixed in situ following the procedures of Sanders, Beinart, Stewart, Delong, and Girguis (2013). As shown in Fig. S4, the samples were highly correlated, with only 205 DEGs found (Sheet S2). For most of genes discussed here, there was no significant difference in expression (Fig. S4). Thus, we believe our sampling strategy to be reflective of natural expression status.
With regard to cellular homeostasis, different types of molecular chaperones (HSP70, DnaJ homologs) and cytochrome P450s were identified in D-DEGs and S-DEGs, but thioredoxin and various caspases were only in D-DEGs while glutathione peroxidase (GPx) and peroxidases were only found in S-DEGs (Figure 3e). In addition, diverse transporters were identified in both D-DEGs and S-DEGs with different types, while sulphate transporters were specific to D-DEGs.
In terms of immune-related DEGs, C1q domain-containing proteins (e.g., MgC1q13, MgC1q32, MgC1q51 and MgC1q89) were only found in D-DEGs, while lectins (such as C-type lectin and galactoside-binding lectin) were specific to S-DEGs. Meanwhile, different types of TLRs were identified in D-DEGs (TLR2) and S-DEGs (TLR4 and another two TLRs), respectively.
There were also many lineage-specific DEGs, which might be important for environmental adaptation. DEGs such as lysosomal cathepsins, V-type H + -transporting ATPases, caveolin-1, actin beta/ gamma 1, coronin-7 and cation-dependent mannose-6-phosphate receptor were identified in D-DEGs, which might be involved in host-symbiont interactions ( Figure 5). On the other hand, two DEGs, cryptochrome and period circadian protein, were found only in T A B L E 2 The significantly reduced pattern recognition receptor (PRR) gene families in deep-sea mussels compared with their shallow-sea relatives. Gene numbers of each gene family is given F I G U R E 3 The statistics of Pfam enrichment (a-d) and heatmap of consistently differentially expressed orthologs between deep-sea and shallow-sea mussels (DEGs) discussed in this study (e). D-SSGs: special orthologs shared only in two deep-sea mussels; S-SSGs: special orthologs shared only in two shallow-sea mussels; D-PSGs: orthologs which were consistently positively selected in deep-sea mussels; S-PSGs, orthologs which were consistently positively selected in shallow-sea mussels; CR, circadian rhythm; BAMA, B. manusensis; BAPL, B. platifrons; MOKU, M. kurilensis; PEVI, P. viridis. Legends were given on the right side of each chart S-DEGs. They are likely related to the regulation of circadian rhythm in shallow-water mussels.
In addition, many other lysosomal activity-related ortholog (such as glycosidases and syntaxins) also showed higher expression in deep-sea mussels, while some autophagy-related orthologs (such as ATG5, ATG10 and ATG14) were expressed consistently more highly in shallow-water mussels (see Figure 3e), although not significantly so (p adj > .05).

| Cell homeostasis maintenance under stress conditions
Both deep-sea mussels and shallow-water mussels might encounter . Such gene families were all found in our results. However, many of them, including genes annotated as ubiquitination related, HSPs (heat-shock proteins) and HIP (heavy metal-binding protein), were reduced in the deep-sea mussels (Sheets S3 and S4).
In addition, there were fewer genes with C1q domains in the deepsea mussels. C1qDC proteins may participate in abiotic challenges, as has been reported in oysters . In contrast to F I G U R E 4 Phylogenetic position of the mussels in Lophotrochozoa. Numbers at the nodes show the support values by 200 bootstrap replications. RAXML was used for the phylogenetic analysis under a GTRCATI model traditional hypotheses, we did not find evidence from deep-sea mussels supporting the hypothesis that bivalves adapt to their environments by the duplicating stress-related genes or obtaining new genes (Zhang et al., 2012). Rather, deep-sea mussels appear to utilize gene sets homologous to their coastal relatives for stress response (Hardivillier et al., 2004;Leignel, Hardivillier, & Laulier, 2005;Serafim, Lopes, Company, Ferreira, & Bebianno, 2008). Based on the fossil records, a modioliform mussel has existed in vents and seeps long before the speciation of Bathymodiolus (Distel, 2000). The fossil species from the ancient chemosynthetic ecosystem suggest that mussel ancestors may already possess the tools for dealing with various environmental stresses. That is, in accordance with previous understanding, mussels are capable of tolerating a broad range of environmental stressors and adapt to diverse habitats.
Although almost all the antistress genes were identified in both deep-sea and shallow-water mussels, we found that most DEGs were not enriched in certain gene families. While some members were upregulated in the deep-sea mussels, other members in the same family were expressed at the same level or even higher in the shallow-water mussels (Figure 3e).
Hypoxia has drawn much attention as a key environmental challenge, as mussels in both coastal zones and deep-sea vents can survive periodic oxygen deficiency (Guezi, Boutet, Andersen, Lallier, & Tanguy, 2014;Hourdez & Lallier, 2007;Huhn, Zamani, von Juterzenka, & Lenz, 2016;Nedoncelle et al., 2015;Roeselers & Newton, 2012;Sui et al., 2016). Globins are small respiratory proteins that bind an oxygen molecule and are found in all domains of life. In vertebrates, different globin types have been identified including haemoglobin, myoglobin, neuroglobin and cytoglobin (Burmester et al., 2004). Thus, they are likely to be involved in adaption to hypoxia stress. Mussels are believed to lack such oxygen-binding proteins, and obtain oxygen by physical gradient exchange. However, an intracellular globin may be present in the gills of the mussel Bathymodiolus heckerae, which might facilitate the adaption to the hypoxia environment (Wittenberg, 1985). In our study, globin-like transcripts were also identified in all four mussel transcriptomes and highly expressed in gill tissues (Sheet S5 and Fig. S5). This suggests a mechanism for hypoxic tolerance, in accordance with previous investigations.
Based on our results, it appears that both deep-sea mussels and shallow-sea mussels can adapt to some extent to diverse harsh environments, including hypoxia, but might respond differently in form and degree.

| Host and symbiont interactions
Marine molluscs are constantly subject to microbes in their natural habitats, including beneficial bacteria (Gerdol & Venier, 2015). Hostmicrobe interactions have been confirmed as important for the health or even survival of many animal species, including mussels (Desriac et al., 2014), and especially for deep-sea mussels. Endosymbiosis inside gill epithelial cells differentiates deep-sea mussels from their coastal counterparts and provides nutrition to the hosts (Barry et al., 2002;Fialamedioni et al., 1994;Streams et al., 1997). We focus on the processes that may determine the fate of the endosymbiotic relationship, including immune response and nutrient T A B L E 3 List of transporters that are positively selected or consistently more highly expressed in deep-sea mussels | 5141 transport. A conceptual scheme on the interactions of hosts and symbiosis is given ( Figure 5).
In bathymodiolin mussels, PRR-coding genes are highly expressed in the gills compared with other tissues, likely due to their mucin-5AC caveolin-1

CYPs
Galectins F I G U R E 5 A conceptual scheme on the interactions of host-symbiont in the bacteriocyte. The symbionts are located inside vacuoles (green) surrounded by the bacteriocyte cytosol (pink). An overview of major processes discussed in this article is shown. The underlined bold indicates the orthologs, which are positively selected in deep-sea mussels; the red up arrow indicates the orthologs, which are expressed consistently more highly in deep-sea mussels; the blue down arrow indicates the orthologs which are expressed consistently more highly in shallow-sea mussels; the dots inside lysosomes indicate lysosomal acid hydrolase, some of which are positively selected and possess higher expression level in deep-sea mussels (e.g., cathepsins). Key genes for the regulation of lysosome activities (mTOR and TFEB) are neither positively selected nor more highly expressed in the deep-sea mussels and are coloured grey. TLRs, Toll-like receptors; SvRs, scavenger receptors; SLCs, solute carrier family; ABCs, ATP-binding cassette transporters; SOD, manganese superoxide dismutase; GST, glutathione S-transferase; CYPs, cytochrome P450s; ATGs, autophagy-related proteins; M6P, mannose-6-phosphate; M6PR, cation-independent mannose-6-phosphate receptor-like; ATPeV, vacuolar H + -ATPase; VPSs, vacuolar protein sorting proteins; VAMPs, vesicle-associated membrane proteins; SNAPs, synaptosome-associated proteins; HGSNAT, heparan-a glucosaminide N-acetyltransferase; mTOR, mammalian target of rapamycin; TFEB, transcription factor EB involvement in symbiont maintenance (Ponnudurai et al., 2016;Wong et al., 2015). However, similarly higher expression of immunerelated genes in gills was also found in shallow-sea mussels, which may be explained by the fact that gills are in close contact with surrounding water and accompanied microbes (Philipp et al., 2012).
Thus, caution should be taken when screening genes related to endosymbiosis using comparative transcriptomes from different tissues. In this work, we screened candidate genes mediating the symbiosis-host interaction by comparing the expression and sequences between the deep-sea and shallow-water groups.
TLRs are highly conserved and capable of binding a very broad range of PAMPs, thus potentially acting as sensors of bacteria, fungi and viruses . A great diversity of TLRs was identified in transcriptomes of the mussels. However, no novel domain organization pattern was identified in the deep-sea mussels. One TLR family is significantly expanded in the deep-sea mussels ( Fig. S6), and most highly expressed in the gills. In addition, three TLRs are positively selected and another one is upregulated in the deep-sea mussels. These TLR genes are plausible candidates for symbiont recognition. The involvement of TLR and downstream pathways in symbiont regulation has been reported in other marine invertebrates (Wippler et al., 2016) including Bathymodiolus spp. mussels (Wong et al., 2015). Furthermore, though reduced as a whole, four C1qDC proteins are still significantly upregulated in the deep-sea mussels. These C1q domain-containing C1qDC proteins might participate in immune recognition with bacteria and parasites and trigger phagocytosis (Gerdol et al., 2011). The involvement of complementary systems in symbiotic interactions have previously been found in cnidarian species (Poole, Kitchen, & Weis, 2016).
It is noteworthy that immune-related PPRs seem to be reduced in the deep-sea mussels. Numbers of C1qDC proteins were significantly reduced in deep-sea mussels (Table 2). Lectins, including C-type lectins and galactoside-binding lectins, were positively selected or more highly expressed only in the shallow-water mussels.
More TLRs were expanded or more highly expressed in shallowwater mussels. A comparative experiment in B. azoricus and Mytilus sp. found that after being challenged by Vibrio spp., fewer immunerelated genes were upregulated in deep-sea mussels (Martins et al., 2014). Expressions of some immune genes are even inhibited after infection (Martins, Queiroz, Santos, & Bettencourt, 2013), and the repression of complementary systems has been shown to be required for maintaining symbiosis in other marine invertebrates (Poole et al., 2016). Similar results were also observed in B. azoricus: in the absence of symbionts, lysozyme expression was upregulated, and declined after re-acquisition of bacteria (Detree et al., 2016).
PRRs appear to be different in both types and expression levels between deep-sea and shallow-water mussels. This polymorphism might somehow generate different immune functions and pathogen resistance to support individual adaptation. Summarily, the chemosynthetic deep-sea mussels might have evolved with reduced PRRs pattern and some specific PRRs (e.g., TLRs, C1qDC proteins) to allow symbionts entry and to successfully maintain them in the bacteriocytes.

| Acquisition of symbionts
Mucin is composed of large glycoproteins on the epithelia and serves as a first shield against pathogens (Chugh et al., 2015;Johansson, Sjovall, & Hansson, 2013). A secreted-type mucin-5AClike gene was positively selected in deep-sea mussels. Interestingly, prtC (Rieder, Fischer, & Haas, 2005), a mucin digestive enzyme expressed by the symbionts, was also positively selected and highly expressed in the methanotroph symbiont of B. platifrons (L. Li, M. X. Wang, C. L. Li, L. F. Li, unpublished data). Through the co-evolution of mucin and the prtC gene, the symbionts may infect the host through a shortcut by restructuring the mucin in advance. After passing through the mucins, the symbionts may invade the hosts in two ways, either via endocytosis (Won et al., 2003) or phagocytosis (Kadar et al., 2005) according to transmission electron scanning results of gill filaments. Based on our results, these two routes are not mutually exclusive and might coexist. Many genes participating in endocytosis and phagocytosis were identified as positively selected or expressed consistently more highly in deep-sea mussels.

| Mutualism between endosymbionts and hosts
Bathymodiolin mussels are studied as model species of the symbiosis in deep-sea chemosynthetic systems. As reported by numerous researchers, nutrient interactions link the mussel and symbionts. In our results, various transport-related genes are positively selected or more highly expressed in the deep-sea mussels (Table 3). As a result, substrates necessary for endosymbionts' metabolisms (e.g., Cu, Fe, sulphate and bicarbonate) can be transported to the bacteriocytes actively or more efficiently. Carbon fixed by the endosymbionts can be captured by two possible routes: "milking" and "farming" (Barry et al., 2002;Kadar et al., 2008;Streams et al., 1997), as mentioned previously. Although the "farming" mode, involving intracellular digestion of symbionts, is well supported based on ultrastructural and enzymatic evidence (Fialamedioni et al., 1994;Fiala-Medioni et al., 2002;Nelson et al., 1995), counterarguments also exist. For example, Kadar et al. (2008) argued that senescent symbionts might be autolysed by bacterial acid phosphatase cell cycling rather than by lysosomes. Our study provides the first transcriptomic evidence from mussels to support nutrient acquisition from symbionts by lysosomal digestion.
The lysosomal system occupies a central and crucial role in degrading macromolecules delivered by endocytosis, phagocytosis and autophagy (Gray et al., 2016;Krystev & Popov, 1976). In our results, genes involved in almost all key lysosomal activities are identified as being positively selected or possessing higher expression trends in deep-sea mussels than in shallow-water mussels ( Figure 5, Sheet S6). Among them, vacuolar H + -ATPase (ATPeV) mediates the pH in the lysosomal lumen (Rothman, Yamashiro, Raymond, Kane, & Stevens, 1989), providing an acidic environment necessary for an activated catalytic activity of cathepsin (Rossi, Deveraux, Turk, & Sali, 2004). Mannose-6-phosphate receptor (M6PR), coronin-7, syntaxins, vesicle-associated membrane proteins (VAMP), vacuolar protein sorting proteins (VPSs) and synaptosome-associated proteins (SNAP) are needed in the transporting, fusion and assemblage of lysosomal enzymes and structural components (ion channels, transporters and catabolic enzymes) (Luzio, Pryor, & Bright, 2007). Our results suggest a specialized function of lysosomes in deep-sea mussels. Considering the concurrence of the upregulated expression of lysosomal genes and the high endosymbiosis biomass, it is reasonable to suggest a key role of lysosomes in the regulation and digestion of endosymbiotic products.
Meantime, diverse transporters are also identified as under positive selection and consistently higher expression in deep-sea mussels in comparison with shallow-water mussels (Table 3). That is, nutrient transfer in deep-sea mussels is more likely accomplished by lysosomal intracellular digestion, but this does not exclude the host mussels from obtaining metabolites continuously from the symbiont as an auxiliary strategy for nutrition acquisition.
Mussel hosts depend on lysosomal activities to meet their nutritional needs, and the balance of hosts and symbionts must be sustained during the life span of Bathymodiolus individuals. Therefore, the regulation of lysosomal digestion of endosymbionts is necessary.
Transcription factor EB (TFEB), a master regulator of lysosomal biogenesis and function (Magini et al., 2013), and the master growth regulator mammalian target of rapamycin complex 1 (mTORC1) (Martina & Puertollano, 2013;Pena-Llopis & Brugarolas, 2011;Settembre et al., 2012) are the most well-studied regulators of lysosomal activity. Although the transcripts of TFEB and mTOR were identified in all of the transcriptomes, none of these two orthologs have shown any significant differences between habitats, indicating that the regulation of lysosomal systems in deep-sea mussels might depend on other mechanisms. In our data set, there were fixed differences in the expression of caspases and ATG proteins between deep-sea mussels and shallow-water mussels (shown in Figure 3e).
As shown, apoptosis-related genes were upregulated while autophagy-related genes were downregulated in deep-sea mussels. Apoptosis might act as a postphagocytic winnowing mechanism in the symbiotic system (Dunn & Weis, 2009). Numerous apoptotic cells were reported in B. thermophilus with high symbiont content compared to mussels harbouring low symbiont content, which suggest that the mussel host might control their symbiont population through apoptotic processes (Guezi, Boutet, Tanguy, & Lallier, 2013).
However, autophagy is an important cellular process of removal and degradation of microbial invaders (Davy, Allemand, & Weis, 2012). It is reasonable to consider that the suppressed expression of autophagy-related genes may help to maintain endosymbionts. To further elucidate this, a loss/reinfection experiment and a more systematic genome analysis involving host and symbiont are needed.

| CONCLUSION
Over the past few decades, numerous investigations have been conducted into bathymodiolin mussels, but controversy remains in many aspects. Here, our comparative transcriptomic analysis provides an outline of the genetic basis for the adaptations of mussels in diverse habitats, focusing on bathymodiolin mussels' adaptations to their deep-sea chemosynthetic environment. According to our phylogenetic results, bathymodiolins have a coastal origin. Almost all the antistress genes were identified in both deep-sea and shallow-water mussels, indicating that all mussels are capable of adapting to various habitats, but different species might prefer to utilize different members of the same gene family. Compared with shallow-water mussels, one of the most extraordinary trait of the bathymodiolins is their endosymbiosis. Together with the highly expressed C1qDC proteins, lineage-specific or positively selected TLRs were identified in the gills of bathymodiolin mussels, suggesting their possible functions in the symbionts recognition. However, pattern recognition receptors were globally reduced, facilitating the invasion and maintenance of the symbiotic bacteria. Additionally, key genes to support lysosomal activity were all positively selected or more highly expressed in the deep-sea mussels in comparison with shallow-sea mussels. Meanwhile, diverse transporters were also identified, suggesting that nutrient transfer in deep-sea mussels may more likely be accomplished by lysosomal intracellular digestion, but this does not exclude the possibility of obtaining metabolites continuously from the symbiont to host as an auxiliary strategy for nutrition acquisition.
Equilibrium relationship between host and symbionts are crucial for a long-term persistence and stability of the mutualism. For deepsea mussels, the autophagic and apoptotic as well as lysosomal activities might be involved in the coordination of the symbiosis. However, the regulation mechanism is still unclear. To resolve this question, a symbiont loss and reinfection experiment and a more systematic genome analysis are needed.

DATA ACCESSIBI LITY
The data have been deposited with links to BioProject Accession no.

CONFLI CT OF INTEREST
The authors declare no competing interests.

AUTHOR CONTRI BUTION
W.M.X., S.S. and L.C.L. designed the study; Z.P. and W.X.C. collected the samples; Z.P. and S.Y. finished the laboratory work. W.M.X. and S.X.Q. performed the bioinformatics analysis. Z.P. drafted the manuscripts. The manuscript was revised by all of the authors. All authors read, approved and contributed to the final manuscript.