microRNAs Facilitate Comprehensive Responses of Bathymodiolinae Mussel Against Symbiotic and Nonsymbiotic Bacteria Stimulation

Background:As the dominant species inhabiting both cold seeps and hydrothermal vents,Bathymodiolinae mussels are one of the most successful megafauna in the deep sea.They thrive in dark and food-insucient environmentsby harboring sulfur-oxidizing bacteria (SOB)and/or methane-oxidizing bacteria (MOB)ingill bacteriocytesand obtain the majority of their nutrition from them.Many attempts have been made to decode the mechanisms underlying their symbiosis, which yetremained largely undisclosedfor years due to the lack of cultivable symbionts. In the present study,the globalexpression pattern of immune-related genes and miRNAswere surveyed inGigantidasplatifronsduring bacterial challengesusing enriched symbiontsor nonsymbioticVibrio in attempting to reveal the molecular mechanisms underlying chemosynthetic symbiosis. Results: Multiple PRRs such as TLRs, LRRs and C1q were found vigorously modulated during challenges whiledistinctly clusteredbetween symbiotic and nonsymbiotic bacteria stimulation. As downstream of the immune response,dozens of immune effectors including HSP70, P450, CD82 andvacuolar protein sorting-associated proteinwere modulated simultaneously, contributing to the ne tuning of cellular homeostasis, lysosome activity and bacteria engulfment in either symbiotic and nonsymbiotic bacteria challenge.A total of 459 miRNAs were identied in gill tissue of G. platifrons while dozens of themwere differentially expressedduring the challenge.Among these miRNAs, some were also found in differentexpression patternbetween symbiont or nonsymbiontchallenges and targeting apoptosis and phagosome maturation-related genes, including caspase8, inhibitor of apoptosis, cAMP-responsive element-binding protein,IκB, Rab and integrin. Conclusion:It was suggested that G. platifrons PRRs might function cooperativelyto facilitate the specialized immune recognition to MOBs or nonsymbioticbacteria. Meanwhile, a shared expression pattern of immune effectorswas observed between bacterial challenges, indicatingthe conservative response of Bathymodiolinae mussels in promoting the adhesion andengulfment of symbionts and nonsymbiont. Nevertheless, the differentially expressed miRNAs were yet suggested to facilitate specialized modulationinsymbiosis by repressing apoptosis- and phagosome maturation-related genes.With the orchestra of immune-related genes and miRNAs, G. platifronsmussels could therefore maintain arobust immune response against invading pathogens while establishing symbiosis with chemosynthetic bacteria.


Introduction
Symbiosis between microorganisms and animals or plants is considered to be an ingenious innovation of life [1]. Many multicellular organisms can live together with bacteria in symbiotic relationships, either loosely or tightly and in epi-symbiosis or endosymbiosis) [2]. It has been demonstrated that organisms can bene t signi cantly from symbiosis, gaining more metabolic potential, enlarging terrestrial habitats, and even receiving shielding from pathogens or predators [3]. In the wild, there is a diverse range of types of symbiotic relationships between host and bacteria (symbiosis in the present study is de ned exclusively as mutualism rather than commensalism or parasitism). All symbionts evolved from freeliving ancestors before coevolutionary processes occurred that resulted in a mutualistic relationship with a host [4]. The decoding of symbiosis, especially the biological processes underlying the establishment and maintenance of symbiosis, is therefore regarded as crucial for understanding the adaptation and evolution of life and has attracted much attention ever since [1,3].
As the front line in defending the self from non-self, the immune system of a host plays an indispensable role in controlling both the establishment and maintenance of symbiosis [5]. It has been demonstrated that symbionts can be acquired by hosts either horizontally from the environment or vertically from the germ cells of parents (mainly maternal) after immune recognition. Multiple molecular or cellular immune processes, including the excretion of antimicrobial peptides, phagocytosis, and apoptosis, can be vigorously modulated simultaneously, promoting colonization in symbiotic tissue or cells such as the light organ (squid), bacteriome (aphid), and trophosome (tubeworm) [6][7][8]. Once colonized, symbionts are further monitored and controlled by the host immune system, avoiding their overgrowth or drastic decline and maintaining the balance of symbiosis [5,9,10]. It is therefore interesting to know how symbionts were initially discriminated from nonsymbionts and how immune processes were further modulated during establishment or maintenance of symbiosis [11]. With the help of state-of-the-art molecular tools such as genome and transcriptome sequencing, more molecules involved in the immune recognition and signal transduction have now been identi ed and were found diversi ed greatly across species.
As a class of endogenously encoded small non-coding RNAs, microRNAs (miRNAs) are known to play indispensable roles in the post-transcriptional modulation of gene expression [12]. To date, more than 38,000 mature miRNAs have been identi ed across about 271 species (according to miRBase.org) and majority of veri ed miRNAs could repress the translation of target genes after binding with 3'-UTR region [12]. Accordingly, a diversity of biological processes including cell proliferation, growth, differentiation and immune response could be further modulated by miRNAs [13]. Interestingly, some recent studies have also demonstrated the participation of miRNAs in host-symbiont interaction especially in plants. For example, dozens of plant miRNA including miR-171, miR-393 and miR-396 have been found playing crucial role in the symbiosis between root and fungal by targeting nodule signaling pathway, auxin signaling pathway and etc [14]. In contrast with the thorough investigations in plants, few symbiosisrelated miRNAs has been found in animal disregarding the large number of miRNAs identi ed to date.
Recently, some miRNAs in aphid or coral were found highly expressed in symbiont-housing tissue or in response to endosymbiont infection, which shielded new insight in their interaction with symbionts [15,16]. However, how exactly these miRNAs participate in symbiosis has remained largely uninvestigated.
As one of the dominant species in both cold seeps and hydrothermal vents, Bathymodiolinae mussels (Mytilidae: Bathymodiolinae) have been found in symbiosis with bacteria, bearing sulfur-oxidizing bacteria (SOB) and/or methane-oxidizing bacteria (MOB) in specialized epithelium cells of their gill tissue (bacteriocytes) [17]. It has been reported that Bathymodiolinae mussels could acquire symbionts horizontally since their settlement and regain them throughout their life span including the adulthood [18]. Moreover, the symbionts were found rstly distributed in both mantle and gills in juvenile before gradually restricted within gill bacteriocytes [19]. Holobionts of Bathymodiolinae mussel and chemosynthetic bacteria was therefore regarded as an ideal model in investigating both the symbiosis and deep sea adaptation against the extreme environment of seeps and vents (cold, dark, insu cient photosynthesis based organic matter, but rich in methane or H 2 S) [20]. Many studies have been thereafter undertaken to determine the mechanisms beneath their symbiosis [20][21][22][23][24]. For instance, several reports have revealed the participation of PRRs in the innate response of Bathymodiolinae mussels against a Vibrio challenge or long-term acclimatization [25][26][27][28]. However, few studies were conducted with symbiont challenge due to the lack of cultivable symbiotic SOBs or MOBs, leaving the immune recognition and signal transduction underlying the onset and maintenance of symbiosis largely unknown.
Since it was rst discovered in 1987 in Sagami Bay, Gigantidas platifrons (formerly named as Bathymodiolus platifrons) has been found to be dominant in cold seeps and hydrothermal vents of Okinawa Trough and Formosa Ridge of the South China Sea [29][30][31]. It was found that G. platifrons only harbored MOBs in their bacteriocytes, making them an ideal model for investigating the immune response against symbionts. Recent studies found that multiple PRRs, including immunoglobin domain containing proteins, PGRPs, Toll-like receptors (TLRs), and C1qDC proteins were present extensively in the genome of G. platifrons and might play a crucial role in symbiosis [20,32]. Besides, works by our lab have further surveyed the global immune response of G. platifrons after short-time decolonization (symbiont depletion) and found multiple PRRs (such as leucine-rich repeat protein or LRRs) that respond to simultaneous MOBs or nonsymbiotic bacteria challenge [21,23]. However, given that regaining of symbionts in adult Bathymodiolinae was most likely accomplished in symbiotic state instead of aposymbiotic state, it is still necessary to know whether the host immune recognition could differ and whether symbionts could render the host with more robust immune response. Moreover, given the crucial role of miRNAs in host-symbiont interaction across plants and animals, it's also interesting whether Bathymodiolinae mussels could encode miRNAs modulating symbiosis-related process by targeting immune-related genes. In this study, the Bathymodiolinae mussel G. platifrons collected from the Formosa ridge in the South China Sea were challenged with either symbiotic MOBs or nonsymbiotic Vibrio bacteria and subjected to both miRNA and transcriptome sequencing. The aim of the study was to (1) investigate the expression pattern of immune related genes as well as miRNAs in G. platifrons holobionts against both symbionts and nonsymbiotic Vibrio bacteria challenges; (2) decode subsequent immune effects mediated by genes that response to symbiont challenge; and (3) survey potential modulation on symbiosis-related process endowed by Bathymodiolinae miRNAs, hopefully providing more information in the interaction between Bathymodiolinae mussels with their chemosynthetic bacteria.

Results
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 tness of samples. As observed, tissues of fresh collected mussels remained intact while the gill were found composed by numerous homorhabdic laments (Fig. S1 A). With 4', 6-diamidino-2-phenylindole (DAPI) staining, we could see that gill laments 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 (  Table 1). After ltration with reads containing adapters or with over 10% unknown nucleotides or more than 50% low quality bases, 118.40 Gbp quali ed 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 identi cation 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 signi cantly 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 signi cantly 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 identi ed 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 rst 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 log 10 (TPM + 1) values ( Fig. 1C). The bottom and top of the box represented the rst 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 signi cantly 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 upregulated 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 lowdensity lipoprotein receptor-related protein (LRPs), PGRP_scaffold2290, LRR_Scaffold_175.36, LRR74_Scaffold_93.19, TLR2_scaffold1476, CD209_Scaffold_209.75, and low a nity immunoglobulin epsilon Fc receptor (FCER_Scaffold_21.25) were also identi ed 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.
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 con rmation. 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 signi cantly 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), Gprotein coupled receptor (GPR), macrophage migration inhibitory factor (MIF), and lipopolysaccharideinduced 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 immunerelated 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.

Global immune response of Gigantidas against MOBs and nonsymbiotic bacteria
It has been demonstrated that all G. platifrons are in a tight association with type I methanotrophs in their bacteriocytes and can obtain nutrition directly from them [33]. This close relationship between the host and symbiont makes them an ideal model for understanding how organisms recognize their chemosynthetic symbionts [20]. However, the mechanisms controlling the symbiosis between G. platifrons and symbiotic MOBs still remain largely unknown due to the unavailability of cultivable symbionts and accessible mussels. Several methods have been used to harvest symbionts from Bathymodiolinae mussels to date, including enrichment by differential centrifugation and density gradient centrifugation [33][34][35][36]. It has been found that far fewer MOBs can be yielded from density gradient centrifugation compared to differential centrifugation, although their purity is better. In the present study, a modi ed method based on differential centrifugation was applied to obtain symbiotic MOBs, improving the purity with little loss of yield. Enlighted by other immunological studies, an extra step of heating at 56℃ for 30 min was applied for MOBs before they were used for a challenge [37,38]. This procedure deactivated the host proteins without denaturing their tertiary structure, which could minimize in uences brought by byproducts of the MOBs enrichment, such as cytokines and complements while maximize immune response induced solely by MOBs. MOBs successively harvested as described above, along with heat treated V. alginolyticus, were then quanti ed and subjected to injection (Fig. S1, S2).
Recent studies have investigated the expression pattern of the immune-related genes of Bathymodioline during bacterial challenge by qRT-PCR, demonstrating the robust response of the host immune system [25,27,28,39]. However, these studies failed to show the host response globally without state-of-the-art molecular tools. The successful application of next generation sequencing in deep sea mussels now provides a better solution [26,32]. In the present study, expressional alternations of Gigantidas genes during either symbiotic MOB or nonsymbiotic bacterial challenges were surveyed globally. Interestingly, it was found that overall immune responsive genes of Gigantidas mussels against symbiont challenge was far less than that during nonsymbiotic bacterial challenge or in symbionts-depleted Gigantidas mussels or in immune response of shallow mussels such as Mytilus coruscus [23,40]. Notwithstanding, similar phenomenon were also observed in other holobionts such as coral where only dozens of DEG were annotated after symbiont challenge [41,42]. As suggested by Gross et al., the interaction between host and symbionts could undergo pathogenic colonizing stage at rst and then a bene cial stage [5]. Meanwhile, unlike the pathogens, host immune response against symbionts could be highly adapted to protect symbionts rather than eliminating them, which therefore might result in minimized immune response observed here. The mild response caused by symbionts could also be energy-saving as the main purpose of symbiosis is to improve the nutritional state of the two partners. Interestingly, about 5%-15% Gigantidas miRNAs were responsive to either MOB or V. alginolyticus challenges (Fig. 1D). miRNAs are known as crucial modulators for gene expression at post-transcriptional level. These DE miRNAs could also strengthen the immune response of host. Moreover, though only hundreds of genes or miRNAs were found responsive to bacteria, most of them displayed a spatiotemporal-speci c expression pattern between groups. For example, only 39 out of 358 DEGs and 19 out of 79 DE miRNAs were the same at 12 h post MOBs and V. alginolyticus challenges (the number altered to 33 out of 393 genes and 7 out of 34 miRNA at 24 h post challenge) (Fig. 1B,D). It was therefore concluded that Gigantidas might respond against MOBs and V. alginolyticus in two distinct ways. Nevertheless, GO analysis of DEGs and targets of DE miRNAs demonstrated that multiple immune-related processes, such as signal transduction, cellular response to stimuli, immune responses, and cell death, were modulated in both the MOB and V. alginolyticus challenge groups. Noticeably, some immune-related processes and functions were only found in certain groups. For example, immune processes such as neurotransmitter binding, transcription regulation, cellular communication, and biological adhesion, were only regulated after the MOB challenge, while scavenger receptor activity, hormone metabolic processes, and cell killing were vigorously modulated after the V. alginolyticus challenge. These ndings con rmed the involvement of cell communication and cell adhesion in Bathymodiolinae symbiosis, which could also be observed in holobionts such as coral and squid [2,[43][44][45]. Comparatively, scavenger receptor activity and cell killing were also well known in pathogen induced immune response with indispensable role in the elimination of pathogens [44][45][46][47].

Distinct expressional pattern of Gigantidas PRRs against symbiotic MOBs
While diverse immune processes were modulated after MOB challenge, how were they trigged yet remained largely unknown. It is well known that deep sea mussel Gigantidas can only rely on innate immunity for either symbiosis or pathogen elimination while PRRs could play an irreplaceable role by recognizing symbionts or pathogens and further activating the subsequent immune processes [48]. Without immunoglobulins or acquired immunity, how Gigantidas discriminate MOBs and non-symbionts remained largely unknown. The expression pattern of PRRs during both the MOB and V. alginolyticus challenges was then investigated. Consequently, it was found that Gigantidas PRRs were differentially modulated between challenges either at transcriptional level or post-transcriptional level, while the PRRs in the EN12 and EN24 groups shared a more similar pattern than those in the VA12 and VA24 groups (Figs. 2, S7). It seems that different combination of PRRs might function cooperatively as "immunoglobulins" to speci cally recognize different bacteria. Similar results have been reported in other molluscs such as the oyster Crassostrea gigas, where some PRRs were found to be responsive against multiple PAMPs, while others was merely responsive to certain PAMPs or pathogens [49][50][51]. More interestingly, there were 33 PRRs (four PRRs were potentially up-regulated by miRNAs) found to be dramatically up-regulated during the MOB and V. alginolyticus challenges. Among these PRRs, multiple C1q proteins, TLR2_scaffold1476, along with LRR74_Scaffold_342.14 were found to be vigorously modulated after the MOB challenge, while PGRP_scaffold2290, LRR_Scaffold_175.36, LRR74_Scaffold_93.19, TLR2_scaffold1476, TLR4_scaffold2249, and VLR_ Scaffold_1558.11 as well as the remaining C1q proteins were only responsive against the V. alginolyticus challenge. It was suggested that these PRRs might facilitate the specialized immune recognition to MOBs or nonsymbiotic bacteria correspondingly. Though few of them were veri ed in vivo or in vitro, our previous research has found the participation of both C1q, TLR2 and LRR in decolonization of Gigantidas or bacteria challenge, recon rming their role in symbionts recognition [23]. Comparatively, though previously found involved in symbiont recognition of B. septemdierum, B. azoricus, Hydra spp and E. scolopes, LRRs and PGRPs were more likely involved in the recognition of nonsymbionts in G. platifrons given their expression pattern [22,24,52,53]. What's more, C1q proteins were found ubiquitously modulated during the immune response. As reported, C1q proteins were widely expressed and massively expanded in mollucs including mussels [20,54,55]. Given that C1q proteins could bind with diversity of immune-related proteins and acting in concert triggering subsequent immune processes, it was suggested that C1q proteins could function as scaffold of PRRs and dedicate to the immune recognition of either MOBs or nonsymbiotic bacteria correspondingly [56,57].
It is well known that the interaction between PRRs and PAMPs relies greatly on their spatial structure and could therefore vary considerably [58,59]. All of the up-regulated PRRs were clustered according to their protein sequence. Interestingly, all C1q proteins were found to be divided into two clusters, while the majority of proteins in the upper cluster (six of eight proteins) were MOB-responsive. On the other hand, about ve of eight C1q proteins in another cluster were V. alginolyticus-responsive (Fig. 2). Given aforesaid speculation that C1q proteins might act as scaffold of other PRRs and contribute to the recognition of symbionts or nonsymbionts, it was therefore of interest to know how these PRRs were modulated and how the structural divergence in uence the recognition with symbionts. However, the present study failed to answer these questions due to the limited investigation. Nevertheless, the differentially expressed PRRs would undoubtedly modulate the expression of immune effectors and modulators. miRNAs facilitate the host with comprehensive modulation network in symbiosis-related process As mentioned previously, multiple immune processes could be vigorously regulated by Gigantidas either during the symbiosis or pathogen elimination. The expressional alternations of immune-related genes were surveyed manually to recon rm the above conclusion. As a result, multiple immune effectors and modulators, including Casp8, IAP, BIRC, LITAF, MIF, CaM, nAChR, ERP, and catL, were found vigorously modulated at transcriptional level after the challenge. Interestingly, though distinct expressional pattern were found in PRR genes, multiple immune effectors with similar function were found responsive to both MOB and V. alginolyticus challenges. For example, several P450 and HSP70 genes, along with BIRC genes, were down-regulated in both challenges while EPDR genes and nAChR genes were up-regulated during the stress period. Given their conserved function across species, it was therefore suggested that that the robust increase of P450 and HSP70 could dedicate synergistically to the maintenance of homeostasis of Gigantidas during symbiosis or pathogen elimination [60][61][62]. Similarly, EPDR genes are known to play an indispensable role in promoting matrix-mediated cell adhesion while nAChR has a crucial role in the ACh-mediated neuroendocrine-immune system of either vertebrates or invertebrates [63][64][65][66]. The shared expression pattern of the above genes indicated the conservative response of Bathymodiolinae mussels, which could promote the adhesion to MOB and nonsymbiotic V. alginolyticus, facilitate the subsequent symbiosis or elimination. Noticeably, some miRNAs were also responsive to both MOB and V. alginolyticus challenge while putatively promoting the expression level of HSPs (gpl-miR-8386, gpl-novel-28, gpl-miR-4045, gpl-miR-2320), CD82 (gpl-miR-4045), VSP33 (gpl-miR-8386), CDPKs (gpl-miR-2469, gpl-miR-4045, gpl-novel-49) and CREB2 (gpl-miR-4627). Among these genes CD82 and VSP33 are known as crucial receptors or regulators in phagocytosis, and are therefore suggested to promote the engulfment of either symbionts or nonsymbionts.
Besides these shared genes, expressional discriminations of multiple positive regulators of cytokines, such as the MIF, LITAF, and some apoptosis-related genes such as caspase8 could also be observed in different stimulus groups. More interestingly, it was found that multiple Gigantidas miRNAs that targeting above process or genes demonstrated opposite expression pattern in MOB and V. alginolyticus challenges, resulting in a more distinct immune response. For instance, apoptosis has been suggested as an effective way of eliminating pathogens during massive infection, sacri cing the minority and protecting the majority [67][68][69]. Comparatively, repression of apoptosis after recognizing symbionts were therefore expected as symbionts should be unthreatening, if not bene cial, to host cells. Consequently, it was found that Caspase8, an upstream protease that activates the cascade of caspases responsible for cell death, was repressed at transcriptional level after the MOB challenge, while IAPs, which are crucial negative regulators of apoptosis, were repressed in the V. alginolyticus challenge [67,70,71]. In consistency with transcriptome results, it was found that miRNAs targeting anti-apoptotic genes such as IAP (gpl-miR-5887) and CREB2 were repressed after MOB challenge while miRNAs targeting pro-apoptotic genes including caspase8 (gpl-miR-9272) and IκB (gpl-miR-27) were promoted [72]. Considering that most miRNAs were negative regulators in gene translation, apoptosis of host cell were therefore suggested to be repressed when stimulated by MOB instead of nonsymbionts. In addition, genes involved in phagosome localization and lysosome-mediated degradation were also differentially modulated between challenges. For instance, Rab5C and INTB are important molecules in the maturation and translocation phagosome and were suggested repressed by Gigatidas miRNAs (gpl-novel-12, gpl-novel-49) that up-regulated in MOB challenge [73,74]. In comparison, INTB are suggested promoted in nonsymbiont challenge. On the other hand, though phagosome maturation and translocation were speculated repressed by MOB challenge, one catL gene was found signi cantly up-regulated simultaneously. As an important cysteine protease, catL has been found to have a crucial function in lysosome-mediated degradation [75,76]. It has been suggested that lysosome-mediated degradation could play an indispensable role in the nutrition acquisition of Bathymodiolinae mussels as well as the population control of symbionts [20,77,78]. The drastic increase in catL transcripts could therefore enhance the bioactivity of lysosomes in host gills after contacting MOBs and promote above process directly. Actually, modulations on cell apoptosis and lysosome-mediated degradation could also be observed in our previous study during decolonization [23]. It was found that massive IAPs were upregulated while multiple catL genes were down-regulated. Their expression pattern indicated repressed cell apoptosis and lysosome activity during symbiont-depletion and further con rmed our ndings herein.

Conclusions
The present study has demonstrated how Gigantidas respond to symbionts or nonsymbionts by investigating the expression pattern of either protein-coding genes or miRNAs. It is worthy to note that Gigantidas PRRs were differentially modulated in responding to symbiotic MOBs while multiple immunerelated transducers or effectors could be recruited promoting the homeostasis and lysosome activity of host and engulfment of symbionts. Notwithstanding, diversity of immune-related pathways were shared between symbiont-induced or nonsymbiont-induced responses. However, the distinct expression pattern of symbiont-induced miRNAs could further facilitate a more comprehensive modulation network for symbiosis by repressing apoptosis and phagosome maturation. Though the interaction between miRNAs and Gigantidas genes were insu ciently veri ed, the present results have demonstrated the complexity in the symbiosis between Gigantidas mussels and MOBs.

Materials And Methods
Animal collection, maintenance, and bacterial challenge The G. platifrons specimens used in the study were collected from cold seeps in the Formosa Ridge in the South China Sea (22°06'N, 119°17'E) during an expedition by the R/V Kexue in 2017. After retrieval on deck, following collection by a remote operated vehicle (ROV) Faxian, mussels were transferred immediately into the onboard aquarium and maintained at atmospheric pressure in ltered circulating seawater (4℃), with a CH 4 supplement [21]. After acclimation for 48 h, 54 similarly sized G. platifrons mussels (length ranging from 70 to 100 mm) were randomly collected and designed as CT, EN, and VA group correspondingly for subsequent bacterial challenges. Some mussels were subjected to a tissue dissection soon after retrieval and gill tissue was collected for subsequent para n sections, transmission electron microscopy, and MOB puri cation.
Symbiotic MOBs were puri ed using a method reported previously with some modi cations [33]. In short, gill tissue was homogenized on ice and ltered using sterilized gauze and nylon lters, with successive pore sizes of 10, 5, and 3 µm. Flow-though was initially centrifuged under conditions of 300 g at 4℃ for 5 min to discard host cells and then 4,000 g at 4℃ for 15 min to collect symbiotic MOBs. After washing three times using sterilized seawater, enriched MOBs were suspended in sterilized seawater for use. The purity of enriched MOBs was further analyzed by scanning electron microscopy. The genomic DNA of the MOBs was also extracted using a Mollusc DNA kit (Omega) and subjected to fragment cloning and a quantitative real time polymerase chain reaction (qRT-PCR) of the pmoA gene. A standard curve of the pmoA gene was generated simultaneously, given a copy number and Ct value, and used for the quanti cation of enriched MOBs. The nonsymbiotic bacteria, Vibrio alginolyticus (isolated from the macro fauna of Formosa ridge cold seep and kindly provided by Dr. Li Sun from the Institute of Oceanology, Chinese Academy of Sciences), was cultured overnight using 2216E medium at 18℃ and collected by centrifugation. As demonstrated, majority of proteins could be neutralized in activity by heating at 56℃ for 30 min without collapsing its tertiary structure. Therefore, both MOBs and V. alginolyticus were further subjected to heat treatment to deactivate the host protein or extracellular products produced during enrichment or culture. Then, the bacteria (MOB and V. alginolyticus) were diluted to a concentration of 1 × 10 6 copy/mL using ltered seawater before use. Mussels in the CT, EN or VA groups were then challenged with sterilized seawater, MOB and V. alginolyticus (100 µL per individual) respectively and maintained in independent tanks (10 L) before sampling. No mortality was observed during the experiment while gill tissues from three random mussels in each group were then collected at 12 and 24 h post injection. All samples employed for mRNA or small RNA extraction were stored with Trizol reagent (Invitrogen) or liquid nitrogen. Each trial was conducted with three replicates.

Rna Extraction, Library Construction, And Rna-seq Of All Samples
Total RNA for transcriptome sequencing was extracted with Trizol reagent (Invitrogen). Meanwhile, small RNA extraction from gill tissue was conducted using Purelink miRNA isolation kit (Invitrogen) according to the manual. The integrity of total RNA was rst con rmed by both agarose gels and a Bioanalyzer 2100 (Agilent). The purity and concentration of total RNA were then determined using a NanoPhotometer spectrophotometer (Implen) and Qubit Fluorometer (Invitrogen). Quali ed RNA samples were subsequently used for library construction, following the instructions for the Illumina HiSeq 2500 platform. In brief, total RNA for transcriptome sequencing was initially enriched by Oligo(dT) beads for mRNA and subjected to fragmentation afterward. After in vitro transcription for rst-stand cDNA and synthesis for the second-strand, the products were then ligated with sequencing adapters. After PCR ampli cation, all cDNA libraries were nally sequenced by the Illumina HiSeq 2500 platform with pair-end reads. Comparatively, quali ed small RNA was rst subjected to 3' and 5' adapter ligation and ampli ed by PCR. After size selection, the puri ed products were also sequenced using Illumina HiSeq 2500 platform. The resulting sequencing data were then uploaded and deposited at The National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/, BioProject NO. PRJNA540074, PRJNA613553).
Genome mapping and identi cation of differentially expressed genes (DEGs) or miRNA (DE miRNAs) For transcriptome sequencing, the quality control procedure was rst conducted using the raw data by FASTP (https://github.com/OpenGene/fastp). Reads with adapters or containing more than 10% of unknown nucleotides or more than 50% of bases that Q-value ≤ 20 were suggested as low quality and removed automatically. The remaining reads were then mapped with G. platifrons genome using HISAT2 with default parameters. The genome was originally reported by Sun et al and updated by our lab [20]. All mapped reads were subsequently assembled for transcripts by Cu inks while the expression levels of all Gigantidas genes were estimated after normalization by the fragments per kilobase of transcript per million (FPKM) value. The differentially expressed genes (DEGs) between the stimulation and control groups were nally determined by Cuffdiff, with fold changes ≥ 2 and a false discovery rate-adjusted P < 0.05.
For miRNA sequencing, the raw data were also ltered by FASTP toolkit to remove low quality reads (reads containing more than one base that Q-value ≤ 20 or containing unknown nucleotides) or reads with 5' adapter/polyA or without 3' adapter or shorter than 18 nt. The clean tags were then aligned with small RNA deposited in GeneBank database along with Rfam database to remove rRNA, scRNA, snoRNA, snRNA and tRNA. The remaining reads were further mapped with the reference genome of G. platifrons using bowtie-1.00 software to discard these located at exon or intron region. The resulting reads were nally subjected by miRDeep2 software for miRNA identi cation. Noticeably, mature miRNAs and precursor sequences from other species deposited in miRBase were employed for miRDeep2 as references to identify the known and novel miRNAs in G. platifrons. Meanwhile, miRNA candidates with raw count number less than 18 were regarded as low abundance and excluded from subsequent expression analysis. The expression level of all miRNAs were calculated and normalized by transcripts per million (TPM) value. Differentially expressed miRNA (DE miRNAs) between groups were further determined if fold change ≥ 2 and a false discovery rate-adjusted P value ≤ 0.05.

Target prediction of DE miRNAs and Gene ontology (GO) analysis of miRNA targets or DEGs
The target genes of all DE miRNAs were predicated by miRanda using 3'-UTR information of all G. platifrons genes given the genome annotation. Given our experience, stringent parameters with the score threshold raised to 155 and energy threshold adjusted to -23 were set for miRanda software when conducting the prediction.
A GO distribution analysis of all candidate targets or DEGs was conducted by Blast2GO (https://www.blast2go.com/) and further visualized by WEGO (http://wego.genomics.org.cn/). The annotation le of the Gigantidas genome was deposited in Figshare (https:// gshare.com/) under the le name of bapl_v4_annt_with_gene_ID_txt. Immune-related genes were selected manually and subjected to expressional clustering by the pheatmap package (https://cran.rproject.org/web/packages/pheatmap/index.html). Full-length protein sequences of PRR genes were then retrieved based on genome information and subjected to a phylogenetic analysis by Seaview Competing interests the study, and assisted with the analysis and interpretation of the results. ZSZ helped with sample collection and morphological analysis of both mussels and bacteria. CL and LC conducted the mussel sampling during the cruise. LCL conceived the study, coordinatedthe experiment, and helped draft the manuscript. All authorsgave their nal approval for publication.

Figure 1
Overview of differentially expressed genes (DEGs)or miRNAs across samples. (A) A total of 95 and 59 genes were found to be either up-or down-regulated in the EN12 group where Gigantidasplatifronswere challenged with enriched methane-oxidizing bacteria (MOB) symbionts for 12 h and compared to a control (CT12) group where mussels were injected with sterilized seawater for 12 h. Atotal of 182 and 61 genes were found to be signi cantly increased or decreased in the VA12 group where Bathymodiolinae mussels were challenged with nonsymbioticV. alginolyticus when compared to the CT12 group. When mussels were challenged with symbionts for 24 h, the expression levels of 73 and 65 genes were found to be robustly up-or down-regulated. A total of 206 and 82 genes were signi cantly increased and decreased, respectively, in the VA24 group. (B) A veen plot of the above DEGs was subsequently constructed. A total of 36 genes were found to be responsive in both the EN12 and VA12 groups, while 33 genes were regulated vigorously in both the EN24 and VA24 groups. Sixteen of 269 genes were vigorously modulated in both the EN12 and EN24 groups, while51 of 471 genes were found to be responsive in both the VA12 and VA24 groups.  Phylogenetic analysis of differentially expressed pattern recognition receptors (PRRs). All differentially expressed PRRs were subjected to a phylogenetic analysis using their protein sequences by Seaview (maximum likelihood method and LG model with 1000 bootstrap samples). PRRs with similar protein sequence could cluster rstly. The expression pattern of PRRs was also illustrated with colored markings when they were vigorously modulated in corresponding group in comparison with CT12/24 group. As found, many PRRs wereexclusively responsive in the EN groups or VA groups while two PRRs were responsive in both EN and VA challenges. Besides, some PRRs with similar protein sequences could respond synchronously to either MOB or nonsymbioticbacteria challenge.

Figure 3
Heat map and hierarchical clustering of differentially expressed immune effectors. The differentially expressed immune effectors in each group were clustered according to their expression pattern. The alternations of transcripts in each group were also marked with different colors (green if decreased and red if increased).  Schematic diagram of miRNA-mediated immunomodulation network of G.platifronsagainst Vibrio alginolyticus challenge. (A) When the Bathymodiolinae mussels were stimulated by V. alginolyticus for 12 h, more immune-related target genes were found. In detail, PRRs including TLR2, TLR4, LRR74 and VLR were suggested as putative targets of miRNAs decreased in VA12 group. Meanwhile, phagocytosisrelated genes, apoptosis-related genes, along with some immune-related transducers or effectors were also putatively modulated by DE miRNAs in VA12 group. (B) When the stimulus continued to 24 h, only TLR4, TBK1, caspase8, RREB1, CRCT1, MMP and BPI were found being targeted by DE miRNAs in VA24 group. Consistently, target genes were marked in red arrow if miRNAs were down-regulated in VA12 or VA24 group while marked in green when miRNAs were up-regulated.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.