Overview of G. platifrons gill tissue and transcriptome/miRNA sequencing
Some of the deep sea mussel G. platifrons were dissected immediately after retrieve by ROV to clarify the fitness of samples. As observed, tissues of fresh collected mussels remained intact while the gill were found composed by numerous homorhabdic filaments (Fig. S1 A). With 4’, 6-diamidino-2-phenylindole (DAPI) staining, we could see that gill filaments were made of monolayer epithelium cells overlying a central lumen containing haemocytes (Fig. S1 B C). The symbiotic MOBs were distributing exclusively in the apical region of bacteriocytes. Transmission electron micrograph further demonstrated that majority of the gill cells were bacteriocytes while other cells such as ciliated cells and mucous cells were also observed. Noticeable, though most symbionts were within membrane-delimited vacuoles, some were engulfed by lysosomes (Fig. S1 D E).
After on board acclimation, Gigantidas mussels were then challenged with sterilized seawater, enriched symbiont MOBs (Fig. S2 A) or V. alginolyticus (Fig. S2 B) for 12 and 24 h correspondingly and designated as CT12, EN12, VA12, CT24, EN24, and VA24 group. Eighteen transcriptome libraries along with miRNA libraries were further constructed with gill tissue of G. platifrons from above six groups correspondingly. A total of 849.40 M clean reads were obtained for transcriptome sequencing with reads lengths over 127.41 Gbp (Supplementary Table 1). After filtration with reads containing adapters or with over 10% unknown nucleotides or more than 50% low quality bases, 118.40 Gbp qualified data were retained and mapped against the Gigantidas genome by HISAT2. Consequently, over 69.17% of the sequencing reads were successfully aligned with the genome, while the mapping rate of each group ranged from 56.27–82.25% (Supplementary Table 1). Comparatively, a total of 331.79 M clean reads from eighteen libraries were obtained for miRNA sequencing and 232.54 M reads were retained after quality control and engaged for genome mapping. As a result, about 189.73 M reads were mapped with genome and suitable for subsequent analysis such as miRNA identification and expressional evaluation (Supplementary Table 1).
Differentially Expressed Genes And Mirnas During Bacterial Challenge
All the aligned reads were then processed for transcripts assembly and expressional evaluation. As a result, 24,595 genes out of 33,962 were found to be expressing among all groups (Fig. 1A, Supplementary Table 2). In comparison with CT12 group where mussels were injected with sterilized seawater for 12 h, a sum of 95 genes were found significantly up-regulated in the EN12 group, where Gigantidas were challenged with enriched MOB symbionts for 12 h (Supplementary Table 3). Meanwhile, a total of 59 genes were also down-regulated in the EN12 group. When Bathymodiolinae mussels were challenged with nonsymbiotic V. alginolyticus, about 182 genes were found significantly increased at 12 h (VA12 group, in comparison with CT12 group) while 61 genes decreased. When mussels were challenged with symbionts for 24 h, only 73 genes were found robustly up-regulated (EN24) while the transcripts of 65 genes were down-regulated (in comparison with the CT24 group). Comparatively, transcripts of 206 genes were promoted at 24 h post V. alginolyticus challenge while that of 82 genes were down-regulated (VA24 group, Supplementary Table 3). A Veen diagram of these DEGs was subsequently constructed (Fig. 1B). It transpired that about 39 genes were responsive in both the EN12 and VA12 groups, while 33 genes were regulated remarkably in both the EN24 and VA24 groups. Only 17 of 269 genes were found to have been vigorously modulated in both the EN12 and EN24 groups, while 51 of 471 genes were found to be responsive in both the VA12 and VA24 groups.
For miRNA sequencing, a total of 459 miRNAs were identified in gill tissue of G. platifrons (Supplementary Table 4). Among these miRNAs, 386 miRNAs were found conserved across species by sharing same seed region and therefore designated as known miRNAs. A total of 73 miRNAs were first reported given the seed region and suggested as novel ones. Moreover, about 105 miRNAs were found with two more precursors in genome (up to six for gpl-miR-544a, See Fig. S3). The overall expression level of all miRNAs in each groups were then compared using box plot given the log10(TPM + 1) values (Fig. 1C). The bottom and top of the box represented the first and third quartiles of in corresponding group while the line insides the box stood for the median value. It transpired that the median quartile in CT12, EN12 and VA24 groups were similar while that in VA12, CT24 and EN24 groups were similar. Noticeably, the third quartile of EN24 group was significantly lower than the rest groups while that in VA24 group were markedly higher.
The differentially expressed miRNAs (DE miRNAs) were then determined (Supplementary Table 5, Fig. S4). Consequently, the expression levels of 30 miRNAs were promoted in EN12 group while that of 31 miRNAs were repressed when compared with CT12 group. Comparatively, a total of 13 miRNAs were up-regulated in VA12 group and 24 were down-regulated. When challenged for 24 h, only 21 miRNAs were differentially expressed in EN24 group, including 11 increased ones and 10 decreased ones. Similarly, 20 miRNAs were vigorously modulated in VA24 group, among which 14 miRNAs were promoted and six miRNAs were repressed. Among these DE miRNAs, 19 miRNAs were responsive to both EN12 and VA12 group, among which three miRNAs were found in opposite pattern (gpl-novel-47, gpl-miR-7538 and gpl-miR-4981, increased in EN12 group yet decreased in VA12 group). For the rest 16 DE miRNAs that share similar pattern, only four of them were up-regulated (gpl-novel-49, gpl-miR-479a, gpl-novel-72 and gpl-miR-3610). At the meantime, about seven miRNAs were found differentially expressed in both EN24 and VA24 group and none was in opposite pattern. In detail, four miRNA including gpl-miR-9570, gpl-miR-9272, gpl-miR-190 and gpl-miR-4981 were up-regulated while gpl-miR-D16, gpl-miR-5324 and gpl-miR-100 were down-regulated.
Functional Annotation Of Degs And Targets Of De Mirnas
GO annotation of all DEGs was subsequently conducted by Blast2GO and visualized by WEGO. As a result, immune-related functions and processes, such as signal transduction, the cellular response to stimuli, immune responses, and cell death, were found and were suggested to have been modulated during both symbiotic and nonsymbiotic bacterial challenges (Fig. S5). Moreover, genes involved in neurotransmitter binding, transcription regulation, cellular communication, and biological adhesion were also vigorously regulated by Bathymodiolinae mussels in an MOB challenge at 24 h (Fig. S5 B). It was also found that more immune-related processes, such as scavenger receptor activity, hormone metabolic processes, and cell killing, were vigorously modulated during a V. alginolyticus challenge (Fig. S5 C, D).
The target genes of all DE miRNAs were then predicted (Supplementary Table 6). Consequently, a total of 744 unique genes were predicted as putative targets of the DE miRNAs. GO distribution analysis further demonstrated that multiple immune-related processes, such as immune system process and response to stimulus could be modulated by above DE miRNAs (Fig. S6).
Distinct expression pattern of PRRs in symbiotic and nonsymbiotic bacterial challenges
As important molecules in immune recognition, 29 PRRs, including 17 C1q proteins, two IL17, three low-density lipoprotein receptor-related protein (LRPs), PGRP_scaffold2290, LRR_Scaffold_175.36, LRR74_Scaffold_93.19, TLR2_scaffold1476, CD209_Scaffold_209.75, and low affinity immunoglobulin epsilon Fc receptor (FCER_Scaffold_21.25) were also identified as immune responsive genes in either MOB or V. alginolyticus challenges (Supplementary Table 7). The expression patterns of the above PRRs were therefore surveyed. It was found that immune responsive PRRs could cluster distinctly between the EN and VA groups given their expression level (Fig. 2). The PRRs in the EN12 and EN24 groups were found to be initially branched together before clustering with the VA12 and VA24 groups.
Many PRRs were found to be responsive only in the EN groups (Fig. S7). For example, multiple C1q proteins, including C1q protein_Scaffold_559.36, C1q protein_scaffold3283, C1q protein_Scaffold_6.43, C1q protein_Scaffold_230.38, and C1q protein_Scaffold_6.44 were found to be up-regulated only in the EN12 group, while TLR2_scaffold1476, C1qTNF3_scaffold2842, C1q protein_Scaffold_507.62, interferon alpha-inducible protein (IFI) 27_scaffold1314, C1q protein_scaffold3227, and C1q protein_scaffold2249 were found to be up-regulated only in the EN24 group. Similarly, LRP4_Scaffold_1957.6, LRP2_Scaffold_836.16, C1q protein_scaffold2097, C1q protein_Scaffold_2534.1, C1q protein_Scaffold_230.37, C1q protein_Scaffold_146.9, and C1q protein_Scaffold_602.35 were found to be responsive only in the VA12 group, while PGRP_scaffold2290, low affinity immunoglobulin epsilon Fc receptor (FCER)_Scaffold_21.25, LRR_Scaffold_175.36, interleukin (IL)17_Scaffold_442.2, IL17_Scaffold_1168.9, LRR74_Scaffold_93.19, C1q protein_scaffold1716, C1q protein_Scaffold_838.19, and C1q protein_Scaffold_308.16 were found to be up-regulated only in the VA24 group. Two PRRs (CD209_Scaffold_297.5 and LRP2_Scaffold_1332.6) were found to be responsive in both EN and VA challenges.
Expression Pattern Of Immune Effectors Responsive To Bacterial Challenges
As described previously, multiple immune-related processes could be modulated in MOB or V. alginolyticus challenges. The expression patterns of immune effectors during challenges were surveyed for further confirmation. As a result, multiple genes involved in immune signaling transduction, cytokine expression, cell migration, and adhesion and oxidation-redox homeostasis were found to be vigorously regulated (Fig. 3, Supplementary Table 8). In detail, three mammalian ependymin-related proteins (EPDRs), two GTPase IMAP family member 4 (GIMA4) genes, two calmodulin (CaM) genes, and two cytochrome P450 genes, along with the caspase8 (Casp8), heat shock protein 70 (HSP70), and the cathepsin L (catL) genes were significantly modulated in the EN12 group. However, only three HSP70 genes and two E3 ubiquitin-protein ligase TRIM genes, along with the protein mab-21, neuronal acetylcholine receptor (nAChR), and the baculoviral IAP repeat-containing protein (BIRC) genes were found to be vigorously modulated at 24 h post MOB challenge. In addition to the genes mentioned above, immune genes such as the inhibitor of apoptosis (IAP), endoplasmic reticulum resident protein (ERP), G-protein coupled receptor (GPR), macrophage migration inhibitory factor (MIF), and lipopolysaccharide-induced TNF-alpha factor (LITAF) were also found responsive at 12 h post V. alginolyticus challenge. Only two P450 genes, three HSP70 genes, and four TRIM genes, along with the nAChR, BIRC, and the superoxide dismutase (SOD) genes were vigorously modulated when Gigantidas was stressed by V. alginolyticus for 24 h.
Diversity of immune-related signal transducers were targeted by DE miRNAs responsive to bacterial challenges
Besides these differentially expressed PRRs and immune effectors, there were also multiple immune-related genes that were targeted by DE miRNAs from either MOB challenge or nonsymbiont challenge. In detail, four PRRs including TLR4 (TLR4_scaffold2249) and LRRs (LRR_Scaffold_832.4, LRR_ Scaffold_405.8 and LRR74_Scaffold_342.14) in addition with two immune effectors including lysosomal protective protein (CSTA) and matrix metalloproteinase-2 (MMP) were found targeted by miRNAs that differentially expressed in EN12 group (Fig. 4A). Meanwhile, phagocytosis-related receptors or signal transducers such as CD82, vacuolar protein sorting-associated protein 33 (VSP33), Ras-related protein Rab-5C (Rab5C) and integrin beta (INTB), along with apoptosis modulators including IAP, Ras-responsive element-binding protein (RREB1) and cAMP-responsive element-binding protein 2 (CREB2) were also suggested as targets of DE miRNAs in EN12 group. When the Bathymodiolinae mussels were challenged by MOB for 24 h, only one PRR (LRR74_Scaffold_342.14) were found being continuously targeted (Fig. 4B). Notwithstanding, diversity of immune-related transducers such as CaM, NF-kappa-B inhibitor alpha (IκB) and TNF receptor-associated factor 6 (TRAF6), along with caspase8, VSP33 were now putatively being modulated.
When the Bathymodiolinae mussels were stimulated by nonsymbiont V. alginolyticus, more immune-related target genes were found. For example, PRRs including TLR2_scaffold1476, TLR4_scaffold2249, LRR74_Scaffold_342.14 and variable lymphocyte receptor (VLR_ Scaffold_1558.11) were suggested as putative targets of miRNAs decreased in VA12 group (Fig. 5A). Meanwhile, phagocytosis-related genes (CD82, VSP33, INTB), apoptosis-related genes (RREB1, CREB2, Baculoviral IAP repeat-containing protein or BIRCs), along with some immune-related transducers or effectors such as CaM, TRAF6, Calcium/calmodulin-dependent protein kinase type 1 (CDPKs), serine/threonine-protein kinase TBK1 and MMP were also putatively modulated by DE miRNAs in VA12 group. When the stimulus continued, however, only TLR4_scaffold2192, TBK1 caspase8, RREB1, CREB-regulated transcription coactivator 1 (CRCT1), MMP and bactericidal permeability increasing protein (BPI) were found being targeted by DE miRNAs in VA24 group (Fig. 5B).