Bathymodioline mussels dominate deep-sea methane seep and hydrothermal vent habitats and obtain nutrients and energy primarily through chemosynthetic endosymbiotic bacteria in the bacteriocytes of their gill. However, the molecular mechanisms that orchestrate mussel host- symbiont interactions remain unclear. Here, we constructed a comprehensive cell atlas of the gill in the mussel Gigantidas platifrons from the South China Sea methane seeps (1100m depth) using single-nucleus RNA sequencing (snRNA-seq) and whole-mount in situ hybridisation. We identified 13 types of cells, including three previously unknown ones, uncovered unknown tissue heterogeneity. Every cell type has a designated function in supporting the gill’s structure and function, creating an optimal environment for chemosynthesis, and effectively acquiring nutrients from the endosymbiotic bacteria. Analysis of snRNA-seq of in situ transplanted mussels clearly showed the shifts in cell state in response to environmental oscillations. Our findings provide insight into principles of host-symbiont interaction and the bivalves’ environmental adaption mechanisms.
This study provides an important cell atlas of the gill of the mussel Gigantidas platifrons using a single nucleus RNA-seq dataset, a resource for the community of scientists studying deep sea physiology and metabolism and intracellular host-symbiont relationships. The work, which offers solid insights into cellular responses to starvation stress and molecular mechanisms behind deep-sea chemosymbiosis, is of relevance to scientists interested in host-symbiont relationships across ecosystems.
Mutualistic interactions between multicellular animals and their microbiota play a fundamental role in animals’ adaptation, ecology, and evolution 1, 2. By associating with symbionts, host animals benefit from the metabolic capabilities of their symbionts and gain fitness advantages that allow them to thrive in habitats they could not live in on their own1, 3. Prime examples of such symbioses are bathymodioline mussels and gammaproteobacteria endosymbionts4. The bathymodioline mussels occur worldwide at deep-sea chemosynthetic ecosystems, such as cold-seeps, hydrothermal vents, and whale falls4. The host mussel acquired sulfur-oxidising (SOX) and/or methane-oxidising (MOX) symbiont through horizontal transmission at their early life stage3. Their symbionts, which are hosted in a specialised gill epithelial cell, namely, the bacteriocytes5, 6, utilise the chemical energy from the reduced chemical compounds such as CH4 and H2S released from cold seeps or hydrothermal vents to fix carbon and turn into carbon source for the host mussel7, 8. The ecological success of the bathymodioline symbioses is apparent: the bathymodiolin mussels are often among the most dominant species in the deep-sea chemosynthetic ecosystems9, 10. Thus, it is critical to know how the bathymodioline mussels interact with the symbionts and maintain the stability and efficiency of the symbiosis.
The gill of bathymodiolin mussels has adapted to deep-sea chemosynthetic life at the levels of molecule, cell and tissue11, 12. Compared to shallow water mussels13, the gill filaments of bathymodiolin mussels have enlarged surfaces allowing them to hold more symbionts per filament (Supplementary Fig. S1). This adaption requires not only novel cellular and molecular mechanisms to maintain and support the enlarged gill filament structure but also a strong ciliary ventilation system to pump vent/seep fluid to fuel the symbiotic bacteria. Previous studies based on whole-genome sequencing and bulk RNA-seq projects have shown that genes in the categories of nutrient transporters, lysosomal proteins, and immune receptors are either expanded in the host mussel’s genome or up-regulated in the gill12, 14–17, suggesting that these genes involved in host-symbiont interaction. Though providing deep molecular insights, these studies mainly used homogenised tissues that average genes’ expression levels amongst different cell types and eliminate cell and gene spatial distribution information18–21. In addition, the broad expression and function of potentially ‘symbiosis-related’ proteins also greatly limited data interpretation. Therefore, a systemic atlas of gill cell types and the descriptions of cell-type-specific gene expressional profiles are warranted to a better understanding of the host-symbiont interaction and environmental adaptation mechanisms of the bathymodioline symbiosis.
In recent years, single cell/nucleus RNA-sequencing (sc/sn RNA-seq) technologies have become one of the preferred methods for investigating the composition of complex tissue at the transcriptional level22. snRNA-seq has several advantages, such as compatibility with frozen samples, elimination of dissociation-induced transcriptional stress responses, and reduced dissociation bias23. These advantages are significant for deep-sea and other ecological studies because cell and molecular biology facilities are not available in the field. To examine the cellular and molecular mechanisms of the environmental adaptations and host-symbiont interactions in bathymodioline mussels, we conducted a single nucleus RNA-sequencing-based transcriptomic study. We analysed the gill symbiosis using deep-sea musselsGigantidas platifrons24. Our work provides a proof-of-principle for environmental adaptation mechanisms study in non-traditional but ecologically important organisms with single-nucleus RNA- sequencing technologies.
Results and discussion
1. G. platifrons deep-sea in situ transplant experiment and single-cell transcriptomic analysis
We conducted a G. platifrons deep-sea in situ transplant experiment at the “F-Site” Cold- seep (∼1117 m depth) and retrieved three groups of samples (Fig. 1A), as follows: the ‘Fanmao’ (meaning prosperous) group, which comprised the mussels collected from methane-rich (∼40 μM) site of the cold-seep, where the animals thrived; the ‘starvation’ group, which comprised the mussels collected from the methane-rich site and then moved to a methane-limited (∼0.054 μM) starvation site for 14 days before retrieval; and the ‘reconstitution’ group, which comprised the methane-rich site mussels moved to the ‘Starvation’ site for 11 days and then moved back to the methane-rich site for another 3 days.
We next performed single nucleus RNA-sequencing (snRNA-seq) of the gill posterior tip of the three groups of samples using the microwell-based BD Rhapsody platform™ Single-Cell Analysis System (Fig. 1A). After quality control, we obtained 9717 (Fanmao), 21614 (Starvation), and 28928 (Reconstitution) high-quality single nuclei transcriptomic data.
To minimise the dropout effect in snRNAseq, we utilised a reciprocal PCA (RPCA) strategy to project cells in the three samples onto the same space based on conserved expressed genes among them25. This strategy maximises the number of cells per cluster regardless of the organism’s state and therefore maximises genes per cluster26. We revealed 14 cell clusters and their corresponding marker genes (Figs. 1B- D) with Seurat27, 28. Because there are few available canonical marker genes per cell type for G. platifrons and molluscs, we carefully characterised each cell cluster by 1) examining the cluster marker genes’ functions, 2) identifying the expression pattern of cluster marker gene using whole-mount in situ hybridisation (WISH) and/or double-fluorescent in situ hybridisation (FISH) analyses, 3) conducting scanning electron microscopy (SEM) and transmission electron microscopy (TEM) analyses, and 4) evaluating the stability per cluster using a bootstrapping algorithm. We then successfully identified 13 of the 14 cell clusters, which could be generally divided into four groups, including 1) the supportive cells, 2) the ciliary cells, 3) the proliferation cells, and 4) the bacteriocytes. One cell cluster remained ambiguous. Notably, a bootstrapping strategy suggested that all clusters were highly supported when using all three samples (Supplementary Fig. S2A). When using each sample individually, most clusters were also well supported except of the three ciliary cells (Supplementary Fig. S2 B-D). This could be due to that these cells were derived from the same precursor cells, and at least partially because they were represented by fewer numbers of nuclei and thus fewer available genes, especially in the food groove ciliary cell (Supplementary Table S1). The spatial and functional annotation of the cell clusters revealed that the cell clusters worked in close cooperation to support the symbiont and maintain a high-efficiency chemosynthetic system.
2. Cell atlas of G. platifrons gill
The inter lamina cells, marked by high expression of fibrillar-forming collagens, were the cells located between the two layers of basal membranes. These cells were previously identified as amoeboid cells13, and their function has not been explored. The inter lamina cells were densely distributed around the food groove and at the rim of the gill filament. While on the middle part of the gill filament, the inter lamina cells distributed in parallel rows. It might help to connect the two sheets of basal membranes and maintain the spatial integrity of the gill filament. In addition, we observed the enrichment of ribosomal proteins involved in protein synthesis (Supplementary Table S2), indicating a high metabolic rate29 of inter lamina cells. This implies that the inter lamina cells may process the nutrients acquired from the symbiont.
We have also identified two types of BMC populations which have not been recognised before. The BMC1 expressed genes encoding extracellular matrix and adhesive proteins (Fig. 2B and Supplementary Table S2), such as SCO-spondin (Bpl_scaf_27638-4.18), Zonadhesin (Bpl_scaf_1371-2.4), Cartilage matrix proteins (Bpl_scaf_16371-0.14 and Bpl_scaf_3387-5.47), and Collagen alpha-1(XXI) chain protein (Bpl_scaf_51350-0.19). BMC2 (Fig. 2B and Supplementary Table S2) was distinctive from BMC1 in the context of expression (Fig 2A) and marker genes (Fig 2B). Genes encoding hemicentins (Bpl_scaf_38731-0.50 and Bpl_scaf_7293- 0.14), which are known to stabilise the cells’ contact with the basal membrane30, were highly expressed in BMC2. Thus, these two cell types are likely to help building the basal lamina, and stabilising the epithelial-derived cells, such as bacteriocytes and intercalary cells, on the surface of the basal membrane.
Mucus cells are specialised secretory cells with intracellular mucus vacuoles13. In shallow- water filter-feeding bivalves, mucus cells secret mucus, which cooperates with the ciliary ventilation system to capture, process, and transport food particles to their mouths, known as filter feeding mechanism31–34. Deep-sea mussels may not necessarily retain this function since they do not normally acquire food resources such as phytoplankton and planktonic bacteria but obtain most of nutrient through their chemosynthetic symbiont. Thus, it was hypothesised that mucus cells are involved in other biological functions such as immune responses to pathogens18. Herein, genes encoding proteins with microbe-binding functionalities were enriched in mucus cells (Fig. 2B and Supplementary Table S2), such as C-type lectins (Bpl_scaf_58694-0.22 and Bpl_scaf_38133-0.8), C1q domain-containing protein (Bpl_scaf_55216-2.11), and Fibrinogen C domain-containing protein (Bpl_scaf_18519-2.15)35, 36. The expression of similar immune genes was upregulated in a shallow water mussel Mytilus galloprovincialis when challenged by a pathogenic bacteria37. Our WISH and double-FISH analyses showed that mucus cells were embedded within the outer rim cilia and scattered on the gill lamella alongside the bacteriocytes (Figs. 2G and 2K). These data collectively suggest that mucus cells may help mussel maintain the immune homeostasis of gill. Interestingly, WISH and double-FISH analyses showed that mucus cells were also distributed alongside the gill lamella’s food groove and the inner edge, where the density of bacteriocytes is low (Fig. 2G). Because the food groove is the main entry of food practical to the labial palps and mouth38, this distribution pattern implies that mucus cells may be also involved in capturing planktonic bacteria and sending them to the mouth.
Ciliary and smooth muscle cells
A remarkable feature of bathymodioline mussel’s gill is its ciliary ventilation system, which constantly agitates the water and provides the symbiont with the necessary gas39. We identified four types of ciliary cells (Figs. 3A and 3B): apical ridge ciliary cells (ARCCs), food groove ciliary cells (FGCCs), lateral ciliary cells (LCCs), and intercalary cells (ICs), as well as a newly identified type of smooth muscle cells. All ciliary cells were marked by canonical cilium genes, such as genes encoding flagellar proteins, ciliary motor kinesin proteins, ciliary dynein proteins, and ciliary microtubules and were clearly distinguishable by specifically expressed genes (Supplementary Table S3 and Supplementary Fig. S3). ARCCs were characterized by expression of Bpl_scaf_20631-1.16, which encodes homeobox Dlx6a-like protein (Fig 3D and 3F) which is a marker of the apical ectodermal ridge40. FGCCs expressed genes encoding primary cilia development regulator Tubby-related protein 341 (Bpl_scaf_24834-2.3) and primary ciliary cell structural protein TOG array regulator of axonemal microtubules protein 142, 43 (Bpl_scaf_55620-1.3), suggesting the FGCCs are the sensory ciliary cells that gather information from the surrounding environment. Both ARCCs and FGCCs were located around the ventral tip of the gill filament. LCCs were distributed as two parallel rows along the gill lamella’s outer rim and ciliary disks’ outer rim (Figs. 3G and 3H) and had highly expressed genes involved in cilium structure, cilium movement and ATP synthases (Fig. 3B and Supplementary Table S3), indicating that LCCs may have a strong ability to beat their cilia.
Interestingly, we identified a group of previously unreported smooth muscle cells (SMCs) co- localised with the LCCs as showed by WISH (Fig. 3I and Supplementary Fig. S4). SMCs strongly expressed several Low-density lipoprotein receptors (LDL receptor) and LDLR-associated proteins44 (Fig. 3B). In addition, these cells also expressed the angiotensin-converting enzyme- like protein45 (Bpl_scaf-29477-5.10) and the ‘molecular spring’ titin-like protein46 (Bpl_scaf_56354-7.8) (Supplementary Table S3). The expression of these genes could be commonly found in human vascular smooth muscle cells47–49. Collectively, we suspect that SMCs are involved in lateral cilium movement and the gill slice contraction.
ICs are the specialised ciliary cells surrounding the bacteriocytes (Figs. 3B and 3J). WISH analysis showed that the expressions of IC markers had apparent spatial variations (Fig. 3J). The expression of ICs marker rootletin50–52 (Bpl_scarf_22362-6.9) was considerably higher in the ICs close to the frontal edge of the gill filament, and the expression gradually decreased along with the direction of inter lamina water flow (Fig. 3K and Supplementary Fig. S5), implying that ICs ventilate the water flow and the mucus through the gill filaments. Furthermore, compared with the other three types of ciliary cells, the ICs expressed several genes encoding transcription factors involved in determining cell fate (Fig. 3B and Supplementary Table S3), such as transcription factors Sox 8 (Bpl_scarf_40595-5.30) and Wnt pathway cell polarity regulator secreted frizzled-related protein 5 (Bpl_scaf_57424-6.1)53, suggesting that the ICs might also play regulatory roles54–56.
It has long been known that the bathymodioline mussel gill has three types of proliferation cells that are conserved throughout all filibranchia bivalves: the budding zone at the posterior end of the gill where new gill filaments are continuously formed, and the dorsal and ventral ends of each gill filaments57, 58. These “cambial-zone”-like cell populations could continuously proliferate throughout the whole life span of the mussel58. Our snRNA-seq data recognised three types of proliferation cells, which is consistent with previous findings (Fig. 4A, Fig. 4B and Supplementary Table S4). The gill posterior end budding zone cells (PEBZCs) are located on the first few freshly proliferated filaments of the posterior tip of gill (Figs. 4C and 4D). The PBEZCs marker gradually disappeared at around the 11th -12th role of the gill filament, suggesting the maturation of the gill filaments, which is similar to the developmental pattern reported in another deep-sea mussel Bathymodiolus azoricus 59. The dorsal end proliferation cells (DEPCs), which expressed the hallmarks genes of muscular tissue, as well as cell proliferation and differentiation regulators, were the proliferation cells in connective tissue at the dorsal end of the gill slice (Fig. 4A, 4E). The ventral end proliferation cells (VEPCs) were two symmetrical triangle-shaped cell clusters of small symbiont-free cells (Fig. 4F). In VEPCs, genes encoding ribosomal proteins, chromatin proteins, RNA and DNA binding proteins, and cell proliferation markers were all up-regulated (Fig. 4B and Supplementary Table S4), indicating that VEPCs are meristem-like cells59, 60.
It has been hypothesised that new proliferation cells will be colonised by symbionts, serving as the vital mechanism of bacteriocyte recruitment59. We then determined which proliferation cell type gives rise to the mature cells, especially the bacteriocytes, which has little available information regarding their precursor and transmission mode in bathymodioline mussels59, 61. Slingshot trajectory analysis suggested that the PEBZCs might be the origin of all gill epithelial cells, including the other two proliferation cells (VEPC and DEPC) and bacteriocytes (Supplementary Fig. S6). The sole exception was BMC2, which may be derived from VEPC rather than PEBZC. This result is consistent with previous studies which suggested new gill filaments of the filibranch mussels are formed in the gill’s posterior budding zones59. The colonisation by the symbiont might play a crucial role in determining the fate of the bacteriocytes. Noticeably, in G. platifrons, only the pillar-shaped first row of gill filament, comprised of small meristem-like cells, was symbiont-free (Figs. 4D and D’), whereas all the other gill filaments were colonised by symbionts. This pattern of symbiosis establishment is different from that of B. azoricus, in which PEBZCs are symbiont free and are gradually colonised by symbionts released from the bacteriocytes on the adjacent mature gill filaments after maturation59. On mature gill filaments, the DEPCs and VEPCs are seemingly the sources of new cells that sustain the growth of the gill filament from both dorsal and ventral directions, respectively57.
3. Bacteriocytes and host-symbiont interaction
We conducted the whole-mount FISH using a bacterial 16S rRNA probe for symbiont to determine the spatial distribution of bacteriocytes, and their positions relative to the other cell types on the gill filament (Figs. 5A, B). Bacteriocytes covered the majority of the surface of the gill filament, except the ventral tip, the ciliary disk, and the frontal edge (lateral ciliary). The bacteriocytes were surrounded by intercalary cells with microvilli and cilia on the surface (Fig. 5C). Noticeably, the bacteriocytes’ size and FISH signal intensity, which reflected the symbiont number in each bacteriocyte, gradually decreased from the heavily ciliated frontal edge to the inner edge, demonstrating that the symbiont abundance in each bacteriocyte was affected by the environmental gradient11.
The gene expression profile of bacteriocytes aligned well with the morphological analysis, which suggested that the bacteriocytes have structural, metabolic, and regulatory adaptions to cope with the symbiont (Figs. 5D, 5F and 5I). It has been hypothesised that the bacteriocytes extract nutrients from the symbiont through two pathways: lysing the symbiont through the lysosome (the ’farming’ pathway) or directly utilising the nutrient produced by the symbiont (the ’milking’ pathway)62, 63. Our TEM observations clearly detected intracellular vacuole and lysosome system of the bacteriocyte that harbour, transport, and digest the symbionts (Fig. 5D), which is consistent with previous studies18, 64. The snRNA data showed bacteriocytes expressed cellular membrane synthesis enzyme (phosphoethanolamine N-methyltransferase, Bpl_scaf_15282-3.10) and a series of lysosomal proteins such as lysosomal proteases (cathepsins; Bpl_scaf_61711-5.12, Bpl_scaf_46838-6.40 and Bpl_scaf_59648-4.5; lysosomal alpha-glucosidase-like), protease regulators 56 (TMPRSS15 proteins; Bpl_scaf_52188-1.15 and Bpl_scaf_15410-0.9), and lysosomal traffic regulator proteins (rabenosyn-5, Bpl_scaf_54816-0.3 and Bpl_scaf_52809-1.6, lysosomal-trafficking regulator-like isoform X1). Among the proteases, cathepsins were thought to be evolutionary conserved molecular tools that host utilized to control the residence of their symbiont microbes65. They were also highly expressed in symbiotic tissue of other deep-sea chemosynthetic animals, such as vesicomyid clams and vestimentiferan tubeworms66–68. Bacteriocytes also expressed genes encoding cellular vesicle transports (kinesins, Bpl_scaf_14819-0.11, Bpl_scaf_54265-1.13, and Bpl_scaf_4784-1.40)53, potential amino acid transporter57 (Bpl_scaf_36159-5.5), and genes involved in intracellular vesicle transport, such as the FYVE and coiled-coil domain-containing protein 158 (Bpl_scaf_33726-5.7). In addition to form and mobilise early endosomes, the protein products of these genes could transport symbiont-secreted nutrient vesicles to the host (Fig. 5F and Supplementary Table S5), supporting the “milking” pathway.
The symbiont of G. platifrons belongs to type I methanotrophy, of which the core metabolic function is linked with the development of intracytoplasmic membranes leading to a high lipid/biomass content69, 70. Recent lipid biomarker analyses showed that the gill of G. platifrons contains a high amount of bacterial lipids, which are directly utilised by the host to synthesise most of its lipid contents66. Downstream of the bacteriocyte’s metabolic cascade, genes encoding proteins that may be involved in fatty acid/lipid metabolism, such as perilipin2 (Bpl_scaf_27158-3.8), which is the critical protein to form intracellular lipid droplets71, and a variety of fatty acid metabolism enzymes (acetyl carboxylase 2, Bpl_scaf_55250-5.1172; fatty acid desaturase 1-like isoform X1, Bpl_scaf_35916-0.673, 74; long-chain fatty acid-ligase ACSBG2- like isoform X1, Bpl_scaf_28862-1.5 75, 76), were up-regulated, suggesting that the fatty acid could be a major form of nutrients passing from the symbiont to the host mussel.
Additionally, bacteriocytes expressed several solute carriers, including sodium/ascorbate cotransporter (solute carrier family 23 members 1-like, Bpl_scaf_32311-1.19), sodium/potassium/calcium exchanger (sodium potassium calcium exchanger 3-like, Bpl_scaf_14503-1.28), sodium/chloride ion cotransporter (solute carrier family 12 member 3- like isoform X1, Bpl_scaf_44604-3.7), ferrous iron transporter (solute carrier family 40 member 1-like, Bpl_scaf_63447-0.12) and zinc transporter (zinc transporter ZIP14-like, Bpl_scaf_44428- 5.3). The solute carriers are a large family of ATP-dependent transporters that shuttle a variety of small molecules across the cellular membrane77. In several model symbiotic systems, solute carriers play a vital role in host-symbiont interaction by either providing the symbiont with substrates78, 79 or transporting symbiont-produced nutrients to the host 80–82. Previous studies demonstrated that solute carrier genes were expanded in deep-sea chemosynthetic animals’ genome83 and highly expressed in symbiotic organs84 including bathymodioline mussel’s gill12. In bathymodioline mussels, previous bulk RNA-seq studies detected up-regulated expression of a large variety of solute carriers in the gill, which is consistent with the present study, suggesting that solute carriers may play crucial roles in shuttling nutrients in and out of bacteriocytes and in maintaining the suitable intracellular micro-environment for the symbiont85, 86.
4. Cell-type-specific response to environmental stresses
To examine cell-type-specific acclimatisation to environmental changes the expression profiles of differentially expressed genes (DEGs) were compared between the three states (Fanmao, starved and reconstituted animals) of samples collected in our in situ transplantation experiment (full lists of DEGs are shown in Supplementary Table S6-S18). Notably, for each state three mussels were processed but pooled for nuclei extraction before snRNA sequencing (see Method for detail). We tested our hypothesis here that at the starvation site, a relatively low concentration of methane shall upset symbiont metabolism and thus substantially affect symbiont-hosting bacteriocytes by assessing the transcriptional changes per cell type. We calculated the distances between the three sets of samples on the UMAP plot (Fig. 6A) estimating the centroid coordinates per sample per cell type, and calculating the Euclid distance between the centroid coordinates of each pair of states per cell type. The impact of starvation was variable across cell types, as reflected by the cross-state distances (Fig. 6B, green bar).
Starvation resulted in the most significant transcriptional changes in bacteriocytes reflected by a large Fanmao-vs-starvation distance, followed by VEPC and inter lamina cells. On the other hand, starvation had little impact on the transcriptions of ciliary cells and most supportive cells such as the food groove ciliary cell, BMC2, and mucus cells. Fig 6B also shows that after reconstitution (moving the mussels back to the methane-rich site Fanmao for 3 days), the expressional profile of bacteriocytes rapidly changed back, reflected by a large starvation-vs- reconstitution distance and much smaller Fanmao-vs-reconstitution distance. This result coincided with and was supported by our pseudo-time analysis for bacteriocytes (Fig. 6C), showing that bacteriocytes in the reconstitution are in intermediate and transitional states between Fanmao and starvation.
For the bacteriocytes population, we conducted cell trajectory pseudo-timing and detailed DEG analysis to interpret the mechanism of host-symbiont interactions. The branched heatmap showed distinct gene expression patterns amongst the three treatments (Fig. 6D and supplementary Fig. S7). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the bacteriocytes’ DEGs provides an overall view of the pathways enriched in each environmental state (Supplementary Fig. S8 A-C). The genes encoding ribosomal proteins are highly expressed in the bacteriocytes under the methane-rich ’Fanmao’ state (Supplementary Fig. S8A), suggesting an active protein synthesis and cellular metabolism87. Moreover, organic ion transporters (A BCB P-glycol, Bpl_scaf_18613-16.10; canalicular multispecific organic anion transporter 2, Bpl_scaf_52110-4.43; sodium- and chloride- dependent glycine transporter 1-like, Bpl_scaf_13642-5.8), which may be involved in transporting symbiont-produced nutrients, are also highly expressed.
In starved G. platifrons, negative regulators of cell proliferation, such as autocrine proliferation repressors and receptor-type tyrosine phosphatase beta, were up-regulated in bacteriocytes, suggesting repression of cell growth. We also observed the enrichment of genes in the apoptosis pathway (Supplementary Fig. S8B), which is attributed to the upregulation of caspases (Bpl_scaf_25165-3.8 and Bpl_scaf_24225-0.12) and cathepsins (Bpl_scaf_64706-0.16, Bpl_scaf_61711-5.12) which can trigger caspases-dependent cell death, suggesting cellular stress condition. Correspondingly, a gene encoding baculoviral IAP repeat-containing protein (Bpl_scaf_6172-0.3), which can bind caspase and inhibit apoptosis, was up-regulated in bacteriocytes (Fig. 6E). In Drosophila, the baculoviral IAP proteins are important in animal’s response to cellular stress and promote cell survival88, 89. Similarly, the gene encoding MAP3K14 (Bpl_scaf_60908-4.14) and potential E3 ubiquitin ligase RNF213 (Bpl_scaf_25983-1.26, Bpl_scaf_1894-1.45, Supplementary Fig. S9), which were apoptosis suppressor or regulator, were up-regulated90. These results suggest that starved bacteriocytes were carrying out cell- type-specific adjustments to cope with stresses. The E3 ubiquitin ligases could also work as intracellular immune sensors of bacterial lipopolysaccharides91. Thus, the encoded protein may be able to active downstream immunological toolkit to digest the symbiont population for nutrients.
As mentioned above, bacteriocytes obtain nutrients from endosymbionts. KEGG pathway analyses suggested that phagosome, lysosome-related pathways were upregulated in the starvation state, indicating that bacteriocytes were actively digesting the endosymbionts in the starved G. platifrons. Interestingly, protein synthesis activity was more activated in inter lamina cells and VEPCs as supported by the up-regulation of genes encoding ribosomal proteins in both cell populations (Supplementary Tables S6 and S17). This was contrary to the situation in bacteriocytes and may be a consequence of the high activity of “farming” in bacteriocytes.
We anticipated that after moving the starved G. platifrons back to Fanmao site, bacteriocytes and their endosymbionts would “reconstitute” leading to the partly restored “farming” and “milking” pathways. This was confirmed by higher expression of fatty acid metabolic genes in the mussels in the reconstitution state than those in the starvation state, such as long-chain-fatty-acid-ligase ACSBG2-like (Bpl_scaf_28862-1.5), elongation of very long- chain fatty acids 7-like (Bpl_scaf_5959-1.9), and fatty acid desaturase 1-like (Bpl_scaf_35916- 0.6). These findings were also consistent with the result of KEGG analyses. The mitochondrial trifunctional enzyme (Bpl_scaf_28862-1.5) was also highly expressed, suggesting a high level of energy-producing activity. The KEGG analyses showed that glutamatergic synthase was enriched in the reconstitution state in comparison with those in both the Fanmao and starvation states. In the insect aphid–Buchnera endosymbiosis model, the host-produced glutamate could be transported and directly utilised by the symbiont92, 93. Similar glutamate- based host–symbiont metabolic interaction mechanisms were also proposed in the deep-sea mussel Bathymodiolus thermophilus and B. azoricus63, 94. As the pseudotime analyses showed above (Fig. 6C), reconstitution was the intermediate state between Fanmao and starvation. Metabolic and gene regulatory functions of bacteriocytes were still different from that in Fanmao. The KEGG analyses suggested regulatory (Rap1 signalling, retrograde endocannabinoid signalling, chemokine signalling, glucagon signalling, thyroid hormone synthesis, mRNA surveillance pathway, etc.) and metabolic pathways were enriched in the reconstitution state (fatty acid metabolism, salivary secretion, aldosterone synthesis and section, insulin secretion and pancreatic section) (Supplementary Fig. S8C). For example, we detected up-regulation of the genes encoding carbonic anhydrases (Bpl_scaf_33596-7.3 and Bpl_scaf_48274-0.3), electroneutral sodium bicarbonate exchanger 1 (Bpl_scaf_61230-5.11), and globin-like proteins (Bpl_scaf_50392-1.14 and Bpl_scaf_24370-0.5), which could provide the symbiont with carbon dioxide and oxygen necessary for the symbiont’s chemosynthetic metabolism83–94.
5. Summary and outlook
Using the deep-sea mussel G. platifrons as a model organism, we demonstrated the power of integrating snRNA-seq and WISH data in unravelling the mechanism behind animal-microbe symbiosis. The robustness of this strategy showed stable and highly distinguishable expression patterns of each cell type regardless of the different environmental states. We successfully profiled the specific roles of different types of cells, including the previously unknown cell types, in maintaining the structure and function of the gill. The supportive cells that located in between (inter lamina cells) and on the basal membrane’s surface (BMC1 and BMC2) helped maintain the anatomical structure of the basal membrane. Ciliary and smooth muscle cells involved in gill slice contraction and cilium beat allowing the gathering of material and information from the surrounding environment. Proliferation cells (PBEZCs, DEPCs and VEPCs) gave rise to new cells, including bacteriocytes which obtained nutrients from endosymbionts using intracellular vacuoles and lysosomes.
The snRNA analyses also revealed that different cell types collaborated to support bacteriocytes’ functionality. The PBEZCs gave rise to new bacteriocytes that allowed symbiont colonisation. Bacteriocytes attached to the basal membrane and were stabilised by supportive cells. Mucus cells co-localise with bacteriocytes and help to maintain immune homeostasis. The beating ciliary cells controlled the water flow, providing bacteriocytes with necessary inorganic substances from the environment. This new information on cell-cell interaction certainly advanced our overall understanding of how endosymbiotic microbes and host cells communicate and collaborate, which cannot be easily achieved through other methods.
Moreover, the analysis of snRNA data from bacteriocytes has provided insight into the molecular mechanisms employed by the host to maitain and regulate the symbiosis. Notably, the bacteriocytes-enriched transcripts involved in harbouring, digesting symbionts, and transporting nutrients produced by the symbionts were identify and characterised. Ourin situ transplant experiments by moving mussels between methane-rich and methane-limited sites also provided clues of cell-type specific responses to environmental change. Under a methane- limited environment, the staved mussels more actively consumed endosymbionts through the "farming" pathway. After being moved back to the methane-rich environment, mussels produced more glutamates which sustained the regrowth of symbionts. These preliminary findings showed that the deep-sea mussels were able to control their endosymbionts using a set of genes in response to environmental change. Due to the limitation of remotely operated vehicle (ROV) dives and sampling capacity, we had to have pooled samples of each state for nuclei extraction and sequencing. Thus, the cells per cluster could be considered as technical than biological replicates. Although this sampling process strategy has been broadly used in snRNA/scRNA sequencing95, we recognise the possible violation of assumptions in p-value calculation for DEGs between the three states.
Overall, the single-cell spatial and functional atlas developed in the present work will decipher some common principles of symbiosis and environmental adaption mechanisms of animals. The workflow developed in the present study could provide insightful references for researchers focusing on the mechanistic study of the biological adaptation of biologically and ecologically important non-model animals.
This study was supported by the Science and Technology Innovation Project of Laoshan Laboratory (Project Number No. LSKJ202203104), the National Natural Science Foundation of China (Grant No. 42030407), Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou)(HJ202101, SMSEGL20SC01), Major Project of Basic and Applied Basic Research of Guangdong Province (2019B030302004), the Research Grants Council of Hong Kong (C6026-19G-A and 16101822) and the Guangdong Natural Science Funds for Distinguished Young Scholar (2022B1515020033). We appreciate all the assistance provided by the crew of R/V ’Kexue’ and the operation team of ROV ’Faxian’.
W.H., K.C., P.Y.Q. and L.C.L designed the study. Z.H., L. Cao, Z.Z.S., L. Chao and W.M.X conducted the in situ transplant experiments and collected the samples; Z.Q.Y. and W.H. constructed the snRNA-seq library. H.K. performed the bioinformatics analysis. W.H., Z.H., L. Cao, J.L., H.C. and L.Z. performed the experiments. This paper was written by W.H., H.K., K.C., P.Y.Q. and L.C.L. All authors revised the manuscript. All authors read, approved, and contributed to the final manuscript.
The authors declare no competing interests.
Materials and methods
Deep-sea in situ transplant experiment and sample collection
In situ transplant experiment was conducted at the ’F-site’ cold seep during the R/V ’Kexue’ 2020 South China Sea cold seep cruise. The overall design of the in situ transplant experiment and environmental states are shown in Figure 1. First, mussels in the methane-rich ’Fanmao’ site (meaning ’prosperous’ site, 22°06′55.425″N,119°17′08.287″E, depth 1117m) were scoped into three nylon bags with approximately 10 mussels in each bag. Then, two bags of mussels were transplanted to the low methane Starvation site (Figure 1, 22°07′00″N, 119°17′07.02″E, depth 1147.42m). After 11 days of transplantation, one bag of mussels in the Starvation site was moved back to the ’Fanmao’ site. On the 14th day of the transplant, three bags of mussels: one bag from the Fanmao site (designated as the ’Fanmao’ sample), one bag from the starvation site (designated as the ’Starvation’ sample) and one bag of mussels which were first transplanted to the Starvation site for 11 days and moved back to the Fanmao site for 3 days (designated as the ’Reconstitution’ sample) were all retrieved by ROV Fanxian. The mussels were kept in a hydraulic pressure-sealed biobox during the ascending of the ROV. The biobox is made of heat- insulated material, which will prohibit heat exchange with warm surface water. The samples were immediately processed once onboard R/V Kexue. For snRNA-seq, the posterior end tip of the mussel’s gill was dissected (Figure 1C), snap-frozen with liquid nitrogen and then stored at −80 °C until use. For WISH, ISH and FISH, the gills of the mussels were first fixed with 4% paraformaldehyde (PFA, prepared with autoclaved 0.22 μM membrane filtered in situ seawater) at 4 °C overnight. Then, the gill was washed with ice-cold 1× PBS three times, dehydrated and stored in 100% methanol at −20 °C.
Single-nucleus RNA-sequencing of G. platifrons
The gill nucleus was extracted using the Nuclei PURE prep nuclei isolation kit (Sigma-Aldrich). For each sample, posterior end tips from 3 individual mussels were randomly selected and pooled together. The posterior tips were homogenised in 10 mL of ice-cold lysis solution (Nuclei PURE lysis buffer with 0.1% Triton X100 and 1 mM DTT) on ice. The cell nuclei were then separated by sucrose gradient centrifugation according to the manufacturer’s protocol. The cell nuclei pellets were washed by re-suspending in ice-cold DPBS–BSA solution (1× DPBS, 0.04% nuclease-free BSA, 0.01% RNase inhibitor; Takara). The nucleus was spun down by centrifugation at 500×g for 5 min at 4 °C. This step was repeated to remove the contaminants from the cell plasma. Finally, the nucleus was re- suspended in the DPBS–BSA solution. The concentration of the cell nucleus was counted by Cell Countess II. Single-nucleus RNA-seq libraries were then constructed with the BD Rhapsody™ single-cell analysis system using the BD Rhapsody WTA Amplification Kit according to the manufacturer’s protocol. The libraries were then sequenced with the Illumina HiSeq 4000 platform. The clean reads of three datasets were submitted to the National Centre for Biotechnology Information Sequence Read Archive database (Bioproject: PRJNA779258).
Three raw datasets (Fanmao, Starvation and Reconstitution) were processed individually by following the BD single-cell genomics analysis setup user guide (Doc ID: 47383 Rev. 8.0) and BD single-cell genomics bioinformatics handbook (Doc ID: 54169 Rev. 7.0). In brief, a G. platifrons genome reference v1.0 (the genome could be downloaded at Dryad Digital Repository http://dx.doi.org/10.5061/dryad.h9942)15 was indexed and built using STAR 2.5.2b. Gene names were prepared according to the BD-single cell guideline. The sequencing data were then mapped to the indexed genome using the BD Rhapsody single-nucleus pipeline v1.9 with default parameters to generate an expression matrix per dataset.
Each data matrix was converted to a SingleCellExperiment data format, and empty barcodes were removed using the emptyDrops function of DropletUtils v3.1496. Then, we converted and processed the data in Seurat v328. We removed cells that had <100 or > 2500 genes and <100 or > 6000 unique molecular identifiers (UMIs). We also removed genes that had <10 UMIs in each data matrix. Then, we log-normalised the data and used DoubletFinder v2.097 to remove potential doublet, assuming a 7.5% doublet formation rate. The numbers of retained nuclei were 9717, 21614 and 28928 for Fanmao, Starvation and Reconstitution data, respectively (Supplementary Table S19). We used the top 3000 highly variable genes for principal component analysis (PCA) and a reciprocal PCA approach to integrating the three datasets98.
The first 40 principal components (PCs) were used for uniform manifold approximation and projection (UMAP) dimensional reduction and following clustering (Supplementary Fig. S10). We used an empirical parameter, a resolution of 0.2, in the FindClusters function. The average expression of each gene per cluster were also provided as supplementary data S1. Based on the WISH results (see below), we manually combined the clusters with similar gene expressions and WISH patterns. Eventually, we assigned all the cells to 14 reliable clusters (Supplementary Fig. S11). We adopted a bootstrap sampling strategy to evaluate the stability for each cluster99. We ran bootstrapping with 100 iterations using cells in the three samples together and in each sample individually. We determined marker genes per cluster using the FindAllMarkers function based on the Wilcoxon rank sum test. The marker genes were used for KEGG enrichment using ClusterProfiler v4.2100.
Cell trajectory analysis
We conducted slingshot trajectory analyses using Slingshot101 v2.2.1 for a subset of clusters (’intercalary cells,’ ’dorsal-end proliferation cells,’ ’ventral-end proliferation cells,’ ’mucus cells,’ ’basal membrane cells 1,’ ’basal membrane cells 2’ and ’bacteriocytes’) to explore the developmental trajectory of cells, assuming that all these cells were developed from the same precursor (’posterior-end budding zone’). We performed PCA for the subsampled data set and conducted dimensional reduction using phateR v1.0.7 for slingshot analyses. We used the Potential of Heat-diffusion for Affinity-based Trajectory Embedding (PHATE) because it could better reveal developmental branches than other tools102.
We also examined the effect of deep-sea transplant experiments on shaping gene expression patterns by comparing the expression levels amongst the three different states of a given cluster. This comparison was done for the bacteriocytes. The biased number of cells per state could affect the results of the dimensional reduction and calculation of marker genes, and the sequenced nucleus per state was unbalanced; therefore, we first downsampled the cells per cluster per state to about the same number. Thereafter, we performed PCA and dimensional reduction using UMAP and PHATE, calculated marker genes and conducted slingshot trajectory and KEGG enrichment analyses as mentioned above for each cell type.
Phylogenetic estimation was conducted for E3 ubiquitin ligase RNF213 genes. We downloaded RNF213 genes from GenBank for representative vertebrate and invertebrate species across the tree of life. The amino acid sequences were aligned with theG. platifrons sequences annotated as E3 ubiquitin ligase (Bpl_scaf_25983-1.26 and Bpl_scaf_1894- 1.45) using MAFFT v7.450. The alignment was used to estimate a maximum-likelihood best gene tree and to calculate bootstrap values on each node using RAxML v8.2.
G. platifrons gill fixation and storage
The gill tissues of G. platifrons were dissected within minutes after the ROV, and samples were retrieved on board R/V Kexue. The gill tissues were briefly washed with ice-cold filtered and autoclaved in situ seawater (FAISW) and then fixed in 4% PFA prepared in FAISW at 4 °C overnight. The gill tissues were washed three times with ice-cold PBST, dehydrated in 100% methanol and stored at -20 °C until use.
Whole-mount in situ hybridisation
For WISH and double FISH analyses, the DNA fragments (∼1000bp) of the targeted genes were first PCR amplified with gene-specific primers (GSPs) pairs using G. platifrons gill cDNA as template (the sequences of targeted genes and gene- specific primers were provided in Supplementary Data S2). The amplified fragments were ligated into the pMD18-simpleT vector (Takara) and transformed intoE. coli. Individual colonies were picked up, and their plasmids were sequenced to confirm the inserts. The templates for in vitro mRNA transcription were amplified using T7 forward GSP (sense probe control) or Sp6 reversed GSPs (antisense probe) combined with either forward or reversed gene-specific primer. Labelled probes and control probes were generated using digoxigenin (DIG)-12-UTP (Roche) or fluorescein-12-UTP (Roche) according to the protocol described by Thisse103 with Sp6 and T7 RNA polymerase, respectively.
For WISH, the connective region at the end of the W-shaped gill filament was cut off, and each gill slice was carefully peeled off with fine-tip tweezers. We dissected gill tissues from 5 individual mussels and pooled all the gill slices together. The gill slices were then rehydrated in 75%, 50% and 25% methanol-PBST (1×PBS with 0.1% Tween 20) for 15 min each, followed by 3×5 min PBST washes. The gill slices were then permeabilised with 2 μg/mL proteinase K in PBST for 30 min at 37 °C. Post-digestion fixation was conducted by fixing the gill slices with 4% PFA in PBST for 30 min at room temperature (RT). After 3×5 min PBST wash to remove the residual fixative, the gill slices were pre-hybridised with a hybridisation mix (HM; containing 50% formamide, 5×saline‒sodium citrate (SSC), 0.1% Tween 20, 10 μg/mL heparin, 500 μg/mL yeast tRNA) for 1 h at 65 °C. For each hybridisation, 5-10 gill slices were added to 400 μL of fresh HM containing ∼0.5 ng/μL DIG-labelled probe. Hybridisation was conducted in a 55 °C shaking water bath overnight. Post-hybridisation washes were performed according to the following steps: the gill slices were first washed 3×15 min with hybridisation washing buffer (50% formamide, 5×SSC, 0.1% Tween 20), followed by 3×15 min 2×SSC with 0.1% Tween 20 and 3×15 min 0.2×SSC with 0.1% Tween 20. The washings were also conducted in a shaking water batch, and all the washing buffers were pre-heated to 55 °C. The samples were washed three times with PBST and then blocked in blocking buffer (2.5% sheep serum, 2% BSA in PBST) for 1 h at RT. Each sample was incubated with 1:10,000 diluted anti-DIG-AP antibody (Roche) at 4 °C overnight.
The samples were incubated with the antibody, then washed for 6×15 min with PBST, followed by washing in 3×15 min alkaline Tris buffer (100 mM Tris-HCl, pH 9.5; 100 mM NaCl; 50 mM MgCl2). The samples were incubated in nitro blue tetrazolium/5-Bromo-4-chloro-3-indolyl phosphate staining solution (Sangon). After the desired expression pattern was revealed, the staining reaction was stopped by 3×15 min PBST–EDTA wash (PBST, 1 mM ETDA). The gill slices were cleared by incubating in 100% glycerol overnight at 4 °C and then mounted on glass slides. The results of control hybridisations (with sense probes) were provided in Supplementary Fig. S12-S15. WISH analyses were repeated with another batch of G. platifrons gill slices samples collected during the R/V ’Kexue’ 2017 South China Sea cold seep cruise to confirm the consistency in expression patterns.
Paraffin embedding and double fluorescent in situ hybridisation
The methanal-dehydrated gill slices were incubated in 100% ethanol, a 1:1 mixture of 100% ethanol and xylene and Xylene twice for 1 h each at RT. The samples were embedded by incubating in Paraplast Plus (Sigma- Aldrich) for 2 h at 65 °C and then cooled down to RT. Sections with 5 μM thickness were cut using a microtome (Leica).
For double FISH, sections were dewaxed by incubating in xylene twice, a 1:1 mixture of 100% ethanol and xylene, 100% ethanol twice, 95% ethanol, 85% ethanol and 75% ethanol for 15 min each at RT. The sections were washed with PBST three times for 10 min each and then permeabilised by 2 μg/mL proteinase K (NEB) in PBST for 15 min at RT. Post-digestion fixation was conducted by incubating the sections in 4% PFA in PBST for 30 min at RT. The sections were washed three times with PBST for 15 min each. Pre-hybridisation was conducted by incubating the sections in HM for 1 h at 55 °C. Then, the in situ hybridisation was performed by incubating the sections in ∼0.5 ng/μL fluorescein (Roche)-labelled probed prepared in fresh HM overnight at 55 °C. The sections were washed three times with 2×SSC for 15 min each at 55 °C, cooled down to room temperature and washed three times with PBST.
A second-round hybridisation was conducted on the DIG-labelled oligonucleotides to label the symbiont. The slices were hybridised for 1 h at 46 °C with 100 ng DIG-labelledG. platifrons symbiont-specific probe in FISH buffer (0.9 M NaCl, 0.02 M Tris-HCl, 0.01% sodium dodecyl sulphate [SDS] and 30% formamide). The slices were then washed with FISH washing buffer (0.1 M NaCl, 0.02 M Tris HCl, 0.01% SDS and 5 mM EDTA) three times at 5 min each at 48 °C. The slices were washed with PBST three times and then blocked with blocking buffer (2.5% sheep serum and 2% BSA in sheep serum) for 1 h at room temperature.
The slices were then incubated with 1:1000 diluted anti-fluorescein–peroxidase (POD) (Roche) overnight at 4 °C, then washed six times with PBST for 15 min each and three times with TNT buffer (100 mM Tris-HCl, pH 7.5; 100 mM NaCl; 0.1% Tween 20) for 15 min each. Afterwards, the fluorescent signal of the G. platifrons gene expression pattern was developed by the TSA fluorescein kit (Akoya Biosciences) according to the manufacturer’s protocol. The slices were washed three times, and the remaining POD activity was quenched by incubation in 1% hydrogen peroxide solution for 1 h at room temperature. Then, the slices were washed three times with PBS, blocked with blocking buffer for 30 min at room temperature, incubated with 1:2500 diluted anti-DIG–POD (Roche) for 2 h at RT and washed with PBST for six times and TNT for three times. The symbiont FISH signal was developed using the TSA Cy3 kit (Akoya Biosciences). Finally, the slices were washed with PBST, stained with DAPI, and mounted with ProLong Diamond Antifade Mountant (Thermo Fisher).
All the WISH samples and whole-mount 16S FISH images were observed and imaged with a Nikon Eclipse Ni microscope with a DS-Ri2 camera. The double FISH slides were imaged with a Zeiss LSM710 confocal microscope.
Electron microscopy analysis
The gill slices of the G. platifrons were dissected and fixed in electron microscopy fixative (2.5% glutaraldehyde and 2% PFA) at 4 °C. For SEM analysis, the samples were dehydrated in a graded ethanol series and then dried at the critical point. The samples were then coated with gold (sputter/carbon Thread, EM ACE200) and observed under a scanning electron microscope (VEGA3, Tescan). For TEM analysis, the samples were rinsed with double distilled water, post-fixed with 1% osmium tetroxide and then washed with double distilled water. The samples were then rinsed, dehydrated and embedded in Ep812 resin. Ultra-thin sections were obtained with an ultramicrotome (70 nm thickness, Reichert-Jung Ultracut E). The sections were then double-stained with lead citrate and uranyl acetate. The cells were observed under a transmission electron microscope (JEM1200, Jeol) operated under 100 kV.
- 1Initial symbiont contact orchestrates host-organ-wide transcriptional changes that prime tissue colonizationCell Host & Microbe 14:183–194
- 2Metaorganisms in extreme environments: do microbes play a role in organismal adaptation?Zoology (Jena
- 3Coming together-symbiont acquisition and early development in deep-sea bathymodioline musselsProc. R. Soc 288
- 4Symbiotic diversity in marine animals: the art of harnessing chemosynthesisNat. Rev. Microbiol 6
- 5Life in the dark: phylogenetic and physiological diversity of chemosynthetic symbiosesAnnu. Rev. Microbiol 75:695–718
- 6A new species of deep-sea mussel (Bivalvia: Mytilidae: Gigantidas) from the South China Sea: Morphology, phylogenetic position, and gill- associated microbesDeep Sea Res. Part I Oceanogr. Res. Pap 146:79–90https://doi.org/10.1016/j.dsr.2019.03.001
- 7Symbioses of methanotrophs and deep-sea mussels (Mytilidae: Bathymodiolinae)Prog. Mol. Subcell. Biol 41:227–249
- 8Phylogenetic characterization of endosymbionts in three hydrothermal vent mussels influence on host distributionsMar. Ecol. Prog. Ser 208:147–155
- 9The vent and seep biota : aspects from microbes to ecosystemsSpringer https://doi.org/10.1007/978-90-481-9572-5
- 10Genetic diversity and connectivity of deep-sea hydrothermal vent metapopulationsMol. Ecol 19:4391–4411
- 113D FISH for the quantification of methane- and sulphur-oxidizing endosymbionts in bacteriocytes of the hydrothermal vent mussel Bathymodiolus azoricusISME J 2:284–292
- 12Insights into deep-sea adaptations and host-symbiont interactions: A comparative transcriptome study on Bathymodiolus mussels and their coastal relativesMol. Ecol 26:5133–5148
- 13Ultrastructure of the gill of the hydrothermal-vent mytilid Bathymodiolus spMar. Biol 92:65–72https://doi.org/10.1007/BF00392747
- 14High-throughput transcriptome sequencing of the cold seep mussel Bathymodiolus platifronsSci. Rep 5
- 15Adaptation to deep-sea chemosynthetic environments as revealed by mussel genomes. NatEcol. Evol 1
- 16An Insightful Model to Study Innate Immunity and Stress Response in Deep-Sea Vent Animals: Profiling the Mussel Bathymodiolus azoricusOrganismal and Molecular Malacology 8
- 17Post-capture immune gene expression studies in the deep-sea hydrothermal vent mussel Bathymodiolus azoricus acclimatized to atmospheric pressureFish Shellfish Immunol 42:159–170
- 18Molecular analyses of the gill symbiosis of the bathymodiolin mussel Gigantidas platifronsiScience 24
- 19RNA imaging. Spatially resolved, highly multiplexed RNA profiling in single cellsScience 348
- 20Single-cell RNA sequencing technologies and bioinformatics pipelinesExp. Mol. Med 50:1–14
- 21Single-cell RNA-seq: advances and future challengesNucleic Acids Res 42:8845–8860
- 22From tissues to cell types and back: single-cell gene expression analysis of tissue architectureAnnu. Rev. Biomed. Data Sci 1:29–51
- 23Advantages of single-nucleus over single-cell RNA sequencing of adult kidney: rare cell types and novel cell states revealed in fibrosisJ. Am. Soc. Nephrol 30:23–32
- 24Using Bathymodiolus tissue stable carbon, nitrogen and sulfur isotopes to infer biogeochemical process at a cold seep in the South China SeaDeep Sea Res. Part I Oceanogr. Res 104:52–59https://doi.org/10.1016/j.dsr.2015.06.011
- 25Bayesian approach to single-cell differential expression analysisNat. Methods 11:740–742
- 26netNMF-sc: leveraging gene- gene interactions for imputation and dimensionality reduction in single-cell expression analysisGenome. Res 30:195–204
- 27Integrated analysis of multimodal single-cell dataCell 184:3573–3587
- 28Comprehensive integration of single-cell dataCell 177:1888–1902
- 29Regulation of ribosomal protein genes: An ordered anarchyWIREs RNA 12https://doi.org/10.1002/wrna.1632
- 30Hemicentin-1 is an essential extracellular matrix component of the dermal–epidermal and myotendinous junctionsSci. Rep 11
- 31A functional interpretation of cilia and mucocyte distributions on the abfrontal surface of bivalve gillsMar. Biol 138:295–309
- 32Gill anatomy and the evolution of symbiosis in the bivalve family ThyasiridaeBiol. Bull 208:200–212
- 33Characterization of mussel gill cells in vivo and in vitroCell Tissue Res 321:131–140
- 34Mucocyte distribution and relationship to particle transport on the pseudolamellibranch gill of Crassostrea virginica (Bivalvia:Ostreidae)Mar. Ecol. Progr. Ser 137:133–138
- 35The C1q domain containing proteins of the Mediterranean mussel Mytilus galloprovincialis: A widespread and diverse family of immune-related moleculesDev. Comp. Immunol 35:635–643https://doi.org/10.1016/j.dci.2011.01.018
- 36Pathogen-derived carbohydrate recognition in molluscs immune defenseInt. J. Mol. Sci 19
- 37Transcriptomic response of mussel gills after a Vibrio splendidus infection demonstrates their role in the immune responseFront. Immunol 11
- 38Regulation of particle transport within the ventral groove of the mussel (Mytilus edulis) gill in response to environmental conditionsJ. Exp. Mar. Biol. Ecol 260:199–215
- 39Feeding behaviour of the mussel, Mytilus edulis: New observations, with a minireview of current knowledgeJ. Mar. Biol 2011
- 40The dlx5a/dlx6a genes play essential roles in the early development of zebrafish median fin and pectoral structuresPLOS ONE 9
- 41TULP3 is required for localization of membrane-associated proteins ARL13B and INPP5E to primary ciliaBiochem. Biophys. Res. Commun 509:227–234https://doi.org/10.1016/j.bbrc.2018.12.109
- 42Proteins that control the geometry of microtubules at the ends of ciliaJ. Cell Biol 217:4298–4313
- 43Crescerin uses a TOG domain array to regulate microtubules in the primary ciliumMol. Biol. Cell 26:4248–4264
- 44LDL receptor-related protein mediates uptake of aggregated LDL in human vascular smooth muscle cellsArterioscler. Thromb. Vasc. Biol 20:1572–1579
- 45Angiotensin-converting enzyme in smooth muscle cells promotes atherosclerosis-brief reportArterioscler. Thromb. Vasc. Biol 36:1085–1089
- 46Pulling single molecules of titin by AFM—recent advances and physiological implicationsPflug. Arch. Eur 456:101–115
- 47Hypercholesterolemia, and vascular smooth muscle cells: A perfect trio for vascular pathologyInt. J. Mol. Sci 21
- 48Heart angiotensin- converting enzyme and angiotensin-converting enzyme 2 gene expression associated with male sex and salt-sensitive hypertension in the Dahl RatFront. physiol 12
- 49Role of titin in nonmuscle and smooth muscle cellsAdv. Exp. Med. Biol 481:265–277
- 50Rootletin organizes the ciliary rootlet to achieve neuron sensory function in DrosophilaJ. Cell. Biol 211:435–453
- 51The Drosophila homologue of rootletin is required for mechanosensory function and ciliary rootlet formation in chordotonal sensory neuronsCilia 4
- 52Striated rootlet and nonfilamentous forms of rootletin maintain ciliary functionCurr. Biol 23:2016–2022
- 53Secreted Frizzled-related proteins: searching for relationships and patternsBioessays 24:811–820
- 54The evolution of protostome GATA factors: molecular phylogenetics, synteny, and intron/exon structure reveal orthologous relationshipsBMC Evol. Biol 8:112–112
- 55FOS-1 functions as a transcriptional activator downstream of the C. elegans JNK homolog KGB-1Cell. Signal 30:1–8https://doi.org/10.1016/j.cellsig.2016.11.010
- 56No backbone but lots of Sox: Invertebrate Sox genesInt. J. Biochem. Cell Biol 42:453–464https://doi.org/10.1016/j.biocel.2009.06.013
- 57Gill Development and its functional and evolutionary implications in the blue musselMytilus edulis (Bivalvia: Mytilidae)Biol. Bull 217:173–188
- 58Cambial zones in gills of BivalviaMar. Biol 31:175–180
- 59Forever competent: deep-sea bivalves are colonized by their chemosynthetic symbionts throughout their lifetimeEnviron. Microbiol 16:3699–3713
- 60Proteomic identification reveals the role of ciliary extracellular-like vesicle in cardiovascular functionAdv. Sci 7https://doi.org/10.1002/advs.201903140
- 61On the growth of bivalve gills initiated from a lobule-producing budding zoneBiol. Bull 205:73–82
- 62Methanotrophic symbiont location and fate of carbon incorporated from methane in a hydrocarbon seep musselMar. Biol 129:465–476
- 63Metabolic and physiological interdependencies in the Bathymodiolus azoricus symbiosisISME J 11:463–477
- 64Methane-based symbiosis in a mussel, Bathymodiolus platifrons, from cold seeps in Sagami Bay, JapanInvertebr. Biol 121:47–54
- 65Infection dynamic of symbiotic bacteria in the pea Aphid acyrthosiphon pisum gut and host immune response at the early steps in the infection processPlos One 10
- 66Lipid Biomarker Patterns Reflect Nutritional Strategies of Seep-Dwelling Bivalves From the South China SeaFront. Mar. Sci 9
- 67Genomic adaptations to chemosymbiosis in the deep-sea seep-dwelling tubeworm Lamellibrachia luymesiBMC Biol 17
- 68Genomic signatures supporting the symbiosis and formation of chitinous tube in the deep-sea tubeworm Paraescarpia echinospicaMol. Biol. Evol 38:4116–4134
- 69Genomic evidence that methanotrophic endosymbionts likely provide deep-sea bathymodiolus mussels with a sterol intermediate in cholesterol biosynthesisGenome Biol. Evol 9:1148–1160
- 70Fatty acid biosynthesis pathways in Methylomicrobium buryatense 5G(B1)Front. Microbiol 7
- 71Adipose differentiation-related protein is an ubiquitously expressed lipid storage droplet-associated proteinJ. Lipid Res 38:2249–2263https://doi.org/10.1016/S0022-2275(20)34939-7
- 72Expression, purification, and characterization of human and rat acetyl coenzyme A carboxylase (ACC) isozymesProtein Expr. Purif 51:11–21
- 73Desaturases and elongases involved in polyunsaturated fatty acid biosynthesis in aquatic invertebrates: a comprehensive reviewFish. Sci 84:911–928
- 74Unique fatty acid desaturase capacities uncovered in Hediste diversicolor illustrate the roles of aquatic invertebrates in trophic upgradingPhilos. Trans. R. Soc 375
- 75Mammalian long-chain acyl-CoA synthetasesExp. Biol. Med. (Maywood 233:507–521
- 76Long-chain acyl-CoA synthetase 2 knockdown leads to decreased fatty acid oxidation in fat body and reduced reproductive capacity in the insectRhodnius prolixusBiochim. Biophys. Acta 1861:650–662
- 77The solute carrier families have a remarkably long evolutionary history with the majority of the human families present before divergence of bilaterian speciesMol. Biol. Evol 28:1531–1541
- 78Dual RNA-sequencing analyses of a coral and its native symbiont during the establishment of symbiosisMol. Ecol 29:3921–3937
- 79Transcriptomic differences between day and night in Acropora millepora provide new insights into metabolite exchange and light- enhanced calcification in coralsMol. Ecol 24:4489–4504
- 80Metabolic co-dependence drives the evolutionarily ancient Hydra- Chlorella symbiosisElife 7
- 81Gene family expansions in aphids maintained by endosymbiotic and nonsymbiotic traitsGenome Biol. Evol 8:753–764
- 82Trading amino acids at the aphid-Buchnera symbiotic interfaceP. Natl. Acad. Sci. USA 116:16003–16011
- 83Host-endosymbiont genome integration in a deep-sea chemosymbiotic clamMol. Biol. Evol 38:502–518
- 84Expression of genes involved in the uptake of inorganic carbon in the gill of a deep-sea vesicomyid clam harboring intracellular thioautotrophic bacteriaGene 585:228–240https://doi.org/10.1016/j.gene.2016.03.033
- 85Ascorbic-acid transporter Slc23a1 is essential for vitamin C transport into the brain and for perinatal survivalNat. Med 8:514–517
- 86Zinc transporter ZIP14 functions in hepatic zinc, iron and glucose homeostasis during the innate immune response (endotoxemia)PLoS One 7
- 87Impaired ribosome biogenesis: mechanisms and relevance to cancer and agingAging (Albany NY 11:2512–2540
- 88Understanding IAP function and regulation: a view from DrosophilaCell. Death. Differ 7:1045–1056
- 89IAPS : More than just inhibitors of apoptosis proteinsCell Cycle 7:1036–1046
- 90Targeting NF-κB-inducing kinase (NIK) in immunity, inflammation, and cancerInt. J. Mol. Sci 21
- 91Ubiquitylation of lipopolysaccharide by RNF213 during bacterial infectionNature 594:111–116
- 92The carbonic anhydrases: widening perspectives on their evolution, expression and functionBioessays 10:186–192
- 93Aphid amino acid transporter regulates glutamine supply to intracellular bacterial symbiontsProc. Natl. Acad. Sci. USA 111:320–325
- 94Comparative proteomics of related symbiotic mussel species reveals high variability of host–symbiont interactionsISME J 14:649–656
- 95Missing data and technical variability in single-cell RNA-sequencing experimentsBiostatistics 19:562–578
- 96EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing dataGenome Biol 20
- 97DoubletFinder: Doublet detection in single- cell RNA sequencing data using artificial nearest neighborsCell Systems 8:329–337https://doi.org/10.1016/j.cels.2019.03.003
- 98Comprehensive Integration of Single-Cell DataCell 177:1888–1902https://doi.org/10.1016/j.cell.2019.05.031
- 99Deciphering Hematopoiesis at single cell level through the lens of reduced dimensionsbioRxiv
- 100clusterProfiler 4.0: A universal enrichment tool for interpreting omics dataThe Innovation 2https://doi.org/10.1016/j.xinn.2021.100141
- 101Slingshot: cell lineage and pseudotime inference for single-cell transcriptomicsBMC Genom 19
- 102Visualizing structure and transitions in high-dimensional biological dataNat. Biotechnol 37:1482–1492
- 103High-resolution in situ hybridization to whole-mount zebrafish embryosNat. Protoc 3:59–69