Immune Tolerance in Mytilus galloprovincialis Hemocytes After Repeated Contact With Vibrio splendidus.

Mediterranean mussels (Mytilus galloprovincialis) are sessile filter feeders that live in close contact with numerous marine microorganisms. As is the case in all invertebrates, mussels lack an adaptive immune system, but they respond to pathogens, injuries or environmental stress in a very efficient manner. However, it is not known if they are able to modify their immune response when they reencounter the same pathogen. In this work, we studied the transcriptomic response of mussel hemocytes before and after two consecutive sublethal challenges with Vibrio splendidus. The first exposure significantly regulated genes related to inflammation, migration and response to bacteria. However, after the second exposure, the differentially expressed genes were related to the control and inhibition of ROS production and the resolution of the inflammatory response. Our results also show that the second injection with V. splendidus led to changes at the transcriptional (control of the expression of pro-inflammatory transcripts), cellular (shift in the hemocyte population distribution), and functional levels (inhibition of ROS production). These results suggest that a modified immune response after the second challenge allowed the mussels to tolerate rather than fight the infection, which minimized tissue damage.

Introduction

Mytilus galloprovincialis is widely distributed throughout the world and has a high ecological and economic impact (1). Due to their status as sessile and filter-feeding animals, bivalves are exposed to a continuous stream of microorganisms, some of them pathogens (2), and environmental pollutants (3). Consequently, mussels have been used as sensors in ecotoxicological studies to monitor the quality of the marine environment (4, 5). Bivalves, such as Crassostrea gigas, are susceptible to diseases that may cause massive mortality (6, 7). Surprisingly, despite living in the same ecosystem and being exposed to the same pathogens, no significant mortality in M. galloprovincialis has been reported (8, 9).

Although lacking acquired immune response, hemocytes, which are the immune cells in bivalves, respond to pathogens with chemotaxis, encapsulation, phagocytic activity and the release of oxygen and nitrogen radicals (10). Moreover, hemocytes can recognize pathogen-associated molecular patterns (PAMPs) via pattern recognition receptors (PRRs) to activate intracellular signaling pathways to finally trigger the synthesis of antimicrobial effectors (11, 12). In this sense, transcriptomic information regarding the modulation of the hemocyte immune response in bivalves remains scarce (13–17). In mussels, several genes related to key immune functions have been characterized over the past few years. These include different molecules involved in specific pathogen recognition, such as C-type lectins (18), C1q domain-containing proteins (19), and proteins with a fibrinogen-related domain (FReD) (20). Compared to other bivalves, mussels are also particularly rich in antimicrobial peptides (AMPs), and myticin C is an example of an important immune effector with chemotactic, antibacterial and antiviral activities (21, 22). Lysozyme, which is able to hydrolyze the central components of the bacterial wall, is another key player in the effector arm of the mussel immune response (23). Gerdol and Venier (24) have reviewed the presence and the interplay between the different molecular components of the mussel immune defense system by using information found in public sequence databases.

In recent years, some studies have suggested that invertebrates may respond to an infection due to some degree of innate immune memory or “priming.” Because of this, exposure to a non-lethal dose of a pathogen could provide protection against later infection with the same pathogen (10, 25). Several studies have recently reported the protection of oysters (Crassostrea gigas) from subsequent infection with ostreid herpesvirus (OsHV-1) by the use of poly I:C as a priming effector (26–28).

The main objectives of the present study were to characterize the transcriptomic and functional response of mussel hemocytes after injection with Vibrio splendidus, which has been reported to produce mortality in mussels (9), and to analyze whether a different type of response could be elicited after a second interaction with the same pathogen. The outcome of this experimental approach might help to reveal the trainability of the mussel immune response and to identify genes associated with this process.

Materials and Methods

Animals

Adult M. galloprovincialis with shells 8–10 cm in length were obtained from a commercial shellfish farm (Vigo, Galicia, Spain) and maintained in open circuit filtered seawater tanks at 15°C with aeration. The animals were fed daily with Phaeodactylum tricornutum and Isochrysis galbana. Prior to the experiments, the animals were acclimatized to aquarium conditions for 1 week.

Experimental Approach

Twenty mussels were marked and notched in the shell, and hemolymph (500 μl) was withdrawn from the adductor muscle of each mussel with a 0.5 mm diameter (25 G) disposable needle. The hemolymph sampled at time zero (t0) was centrifuged at 4°C at 3,000 g for 10 min, and the pellet was resuspended in 500 μl of TRIzol (Invitrogen), immediately homogenized and stored at −80°C until RNA isolation.

After 1 week, 10 mussels were injected in the adductor muscle with 100 μl of filtered seawater (FSW). The other 10 mussels were injected in the same way with 100 μl of a solution of V. splendidus (reference strain, LGP32) at a non-lethal concentration (1 × 107 UFC/ml). One day post injection (24 hpi) and 7 days post injection, hemolymph (500 μl) was sampled again from individual mussels and centrifuged in the same conditions, and the pellet was resuspended in 500 μl of TRIzol (Invitrogen). Samples were immediately homogenized and kept at −80°C until RNA isolation.

After 2 weeks, the 10 mussels injected with FSW were injected again with FSW. The mussels previously exposed to V. splendidus were injected again with a solution of V. splendidus (reference strain, LGP32) at a non-lethal concentration (1 × 108 UFC/ml). One day after the second injection (24 hpi2), hemolymph (500 μl) was sampled again from individual mussels and centrifuged in the same previously described conditions, and the pellet was resuspended in 500 μl of TRIzol (Invitrogen). The samples were immediately homogenized and kept at −80°C until RNA isolation.

Seven days later (7 d), hemolymph (500 μl) was sampled again from the mussels and centrifuged in the same previously described conditions, and the pellet was resuspended in 500 μl of TRIzol (Invitrogen). The samples were immediately homogenized and kept at −80°C until RNA isolation.

Vibrio splendidus Clearance Assessment

The clearance of V. splendidus was assessed to make sure that the second injection was made after a complete overcome of a possible infection. cDNA was synthesized from samples taken at t0, 24 hpi, and 7 days after the first infection with 100 ng of total RNA using an NZY First-Strand cDNA Synthesis Kit (nzytech). Gene expression of V. splendidus 16S and mussel 18S (used as a reference gene) was analyzed in a Stratagene Mx3005P thermal cycler (Agilent Technologies). In addition, ten-fold dilutions of V. splendidus DNA was included in the plate to extrapolate the bacterial load of each sample.

For 16S detection, 5 μl of five-fold-diluted cDNA template was mixed with 0.6 μl of each primer (10 μM), 0.4 μl of 16S probe (10 μM) and 10 μl of Brilliant III Master Mix 2x Ultrafast (Agilent Technologies) in a final volume of 20 μl. For 18S detection 1 μl of five-fold-diluted cDNA template was mixed with 0.5 μl of each primer (10 μM) and 12.5 μl of Brilliant II SYBR Green (Agilent Technologies) in a final volume of 25 μl. The standard cycling conditions were 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s. All reactions were performed as technical triplicates. The relative expression levels of the genes were normalized using 18S as a reference gene following the Pfaffl method. One-way ANOVA was used to analyze differences in normalized gene expression and bacterial load among the studied samples.

RNA Isolation, cDNA Production, and Illumina Sequencing

RNA isolation was carried out using TRIzol (Invitrogen) according to the manufacturer's protocol. RNA purification was performed after DNase I treatment with the RNeasy Mini Kit (Qiagen). Next, the concentration and purity of the RNA were measured using a NanoDrop ND1000 spectrophotometer, and the RNA integrity was verified with an Agilent 2100 Bioanalyzer (Agilent Technologies). Only the best RNA samples (in terms of RNA quantity and quality) from four sampling points (t0, 24 hpi, 24 hpi2, 7 d) were chosen for the preparation of cDNA libraries compatible with Illumina sequencing. The chosen samples were from the mussels control 2, control 3, control 4, infected 1, infected 2 and infected 10. A total of 24 samples were selected for sequencing that consisted of 2 conditions, 4 sampling points and 3 biological replicates. Some of the samples corresponding to 7 days after the first injection were not of sufficient quantity to be sequenced, and therefore this sampling point was excluded from the analysis.

A TruSEq library preparation kit from Illumina was used according to the manufacturer's instructions. Briefly, eukaryotic mRNA was extracted from total RNA using oligo (dT) magnetic beads and cleaved into short fragments using fragmentation buffer. A cDNA library was then prepared from the fragmented mRNA via reverse transcription, second-strand synthesis and the ligation of specific adapters (paired-ends) after cDNA purification using the QIAquick PCR Purification Kit (Qiagen). The amount of cDNA in each library was quantified through spectrofluorometric analysis using the Qubit system. Paired-end sequencing (2 × 100) was performed using an Illumina HiSeq™ 4000 platform by Macrogen Korea.

A schematic representation of the experimental design for the sequenced samples is shown in Figure 1. The raw sequencing data have been deposited in the NCBI Short Read Archive database under the accession ID SRP145077.

Figure 1

Experimental design used for the stimulation of mussels.

Bioinformatics and RNA-Seq

The CLC Genomics Workbench, v.11.0.1 (CLC Bio; Qiagen), was used to process the raw sequencing output for the de novo assembly of the reference transcriptome and to perform the statistical analyses of gene expression by comparing the three biological replicates for the control and infected mussels at different time points. The raw reads were trimmed to remove adaptor sequences, low quality bases (quality score p-value limit of 0.05), and residual sequences shorter than 70 bp. All reads obtained from the 24 libraries were assembled to obtain a complete reference transcriptome with default word size and bubble size parameters. The assembly was cleaned to remove sequences originating from mussel ribosomal RNA and mitochondrial mRNAs, as well as contaminant transcripts from V. splendidus, ciliates and microalgae. These filtering steps were performed with BLASTn analyses (e-value threshold 1e-10) that were carried out in parallel with the reported assembled mussel genome (29) and the targets mentioned above (the Vibrio splendidus genome from strain NCCB 53037, the ciliate Pseudocohnilembus persalinus genome, the Phaeodactylum tricornutum genome and the Isochrysis galbana transcriptome from BioProject PRJNA248394 were used as references). Contigs that produced a more significant hit when compared to the sequences of the putative contaminants than to the mussel genome were discarded. The quality and completeness of the transcriptome were assessed with BUSCO v.3 (30), which was based on the detection of metazoan Benchmarking Universal Single Copy Orthologs (BUSCOs) according to release 9 of OrthoDB.

The reads of each individual mussel and sampling were mapped onto the clean transcriptome with the RNA-Seq tool using the following parameters: mismatch cost = 2, length fraction = 0.8, similarity fraction = 0.8, and maximum hits per read = 10. Differentially expressed genes (DEGs) were identified with a statistical analysis based on generalized linear models and by assuming a negative binomial distribution for the read counts (31). For each set of comparisons, transcripts with an absolute fold change (FC) value > 2 and an FDR-corrected p < 0.05 were considered differentially expressed and retained for further analyses. To find the DEGs at each time point after the Vibrio challenges, the injected samples were compared with their respective controls. To find primed and tolerized genes, the Vibrio challenged samples from the second injection were compared with the challenged samples from the first injection, and the same comparison was made in the controls to confirm that the selected genes were not modulated in control animals.

BLAST Annotation, GO Assignments, and Enrichment Analysis

The transcriptome was functionally annotated with the Blast2GO software (32) by assigning gene ontology (GO) terms based on the significant BLASTx matches found in the UniProt/Swissprot database. To improve the annotation rate, we performed an additional BLASTn analysis against an in-house database, which included all the molluscan sequences present in the NBCI nucleotide database. In both cases, the e-value threshold for annotation was set to 1e-3. Then, functional enrichment analyses of the DEGs (test set) were conducted using the full mussel transcriptome as the reference set. For this purpose, a two-tailed Fisher's exact test was performed with the default parameters and a p-value cut-off of 0.05. The test was performed on the basis of overrepresented biological process (BP) gene ontology terms.

Functional Assays: Hemocyte Distribution, Apoptosis, and ROS Analyses

The previously described experimental design was repeated using eight biological replicates (each replicate included a single mussel) to determine whether functional immune parameters were also affected by a second exposure to the same bacterial pathogen. Hemolymph was collected from the adductor muscle of the eight individual mussels using a disposable syringe, and the cell concentration was adjusted to 106 cells ml−1 with FSW.

The hemocyte populations were evaluated by flow cytometry. Two FSC/SSC gates were created that included both the viable granulocyte and the hyalinocyte populations. Data were acquired using a FACS Calibur flow cytometer (Becton and Dickinson), and the analysis was carried out using CellQuest software (Becton and Dickinson).

To investigate the possibility that the mussels that received a single and two subsequent injections of V. splendidus showed changes in cell death rates, an apoptosis analysis was performed. Hemocytes were centrifuged and resuspended in 1 ml of binding buffer (BB1X). Then, 5 μl of annexin V (Invitrogen) and 5 μl of actinomycin (BD Pharmingen) were added to the cell suspensions. The samples were incubated for 15 min at room temperature in the dark and analyzed by flow cytometry.

The respiratory burst activity of hemocytes was determined by the luminol-enhanced chemiluminescence method (CL) in 96-well-plates. We used 5-amino-2,3-dihydro-1,4-phthalazinedione (luminol, Sigma Aldrich) as a light emitter and phorbol myristate acetate (PMA, Sigma Aldrich) or zymosan A (Sigma Aldrich) to trigger the production of reactive oxygen species (ROS). A stock solution of 0.1 M luminol was prepared in dimethyl sulphoxide (DMSO, Sigma Aldrich) and diluted in FSW to obtain the luminol working solution (final concentration of 10 mM). Zymosan A (20 mg ml−1) was diluted in the luminol working solution to obtain a final concentration of 1 mg ml−1. The PMA stock solution (1 mg ml−1 in ethanol) was also diluted in the luminol working solution to obtain a final concentration of 1 μg ml−1.

One hundred microliters of hemolymph was dispensed into each well of the 96-well-plates. After 30 min of incubation at 15°C, 100 μl of luminol, PMA or zymosan A were added per well. The relative luminescence units (RLU) were measured in a luminometer (Fluoroskan Ascent, Labsystems) six times at intervals of 5 min with an integration time of 1,000 ms for each measurement.

Results

Assembly and Annotation of the Mussel Transcriptome

The sequencing of the individual hemocyte samples yielded an average of 74.11 million raw reads for each of the 24 libraries. The trimming procedure removed, on average, 0.83% of the raw reads, and a total of 1,778 million reads were assembled into a reference mussel transcriptome containing 260,664 contigs with an average length of 512 bp. The reference transcriptome was highly complete, as just 2% of metazoan BUSCOs were absent, and displayed a fragmentation rate equal to 22%, which was in line with what was expected for such a highly heterozygous species (Table 1).

Table 1

Summary of the transcriptome bioinformatics pipeline.

Sample	Raw	Trimmed
C2 t0	78,426,948	99.59%
C3 t0	44,346,854	98.03%
C4 t0	100,814,198	99.59%
I1 t0	17,696,894	98.08%
I2 t0	93,114,098	99.60%
I10 t0	96,780,602	99.64%
C2 24h	95,296,484	99.49%
C3 24h	51,708,988	97.62%
C4 24h	92,661,282	99.65%
I1 24h	52,102,302	98.77%
I2 24h	90,965,875	99.52%
I10 24h	99,262,970	99.55%
C2 24h2	93,665,372	99.96%
C3 24h2	47,632,114	99.09%
C4 24h2	83,629,996	99.62%
I1 24h2	25,833,124	98.25%
I2 24h2	90,636,926	99.60%
I10 24h2	101,411,946	99.64%
C2 7d2	78,528,396	99.55%
C3 7d2	46,511,774	99.36%
C4 7d2	95,180,866	99.69%
I1 7d2	13,593,322	97.43%
I2 7d2	90,614,296	99.31%
I10 7d2	98,236,250	99.57%
Assembly statistics
Contigs	260,664
Range contig length	200–15,624
Average contig length	512
N50	576
Complete metazoan BUSCOs	740/978 (75.69%)
Fragmented metazoan BUSCOs	218/978 (22.29%)
Missing metazoan BUSCOs	20/978 (2.04%)
Blast
Contigs with hit in Uniprot/SwissProt	19.93%
Contigs with hit in mollusks database	42.64%
GO analysis
Annotated contigs	23.35%
KEGG analysis
Pathway assigned contigs	7.27%

Two different BLAST approaches were used to annotate the assembled transcriptome. In brief, 42.64% of the contigs were found to have a significant match in the custom database, which included all the mollusk nucleotide sequences available from NCBI, and Blast2GO was used to annotate 19.93% of the contigs through a BLASTx search of UniProt/SwissProt (Supplementary File 1). Based on these results, gene ontology (GO) terms were assigned to 23.35% of the contigs. Table 1 shows the sequencing output of all the samples and the main metrics of the transcriptome assembly and annotation.

Transcriptomic Response After Injection and Reinjection With the Same Pathogen

We carried out a differential gene expression analysis to gain insights into the dynamics of the transcriptional response of mussel hemocytes to an experimental infection with V. splendidus. First, the clearance of V. splendidus by the injected mussels was confirmed 7 days after the first injection: V. splendidus detection increased 24 hpi and was rapidly controlled 7 days after the injection, returning to control levels (Supplementary File 2). To analyze the transcriptomic response at each sampling point, the Vibrio injected animals were compared to control animals (FSW-injected) (Figure 2A). The monitoring of the transcriptional profiles enabled us to assess whether a second interaction with the same pathogen could elicit a different type of response compared to the response elicited by the first injection. A total of 1,216 differentially expressed genes (DEGs) were detected 24 h after the first injection (24 hpi). However, this number dramatically decreased to 236 DEGs 24 h after the second injection (24 hpi2) and dropped further to 80 DEGs when the transcriptional profiles were compared 7 days after reinjection (7 d) (Figure 2B).

Figure 2

(A) Representative scheme of the comparisons made for the differential expression analysis between injected and control individuals at the analyzed sampling points. The thresholds used to detect statistical significance were a fold change (FC) > |2| and an FDR < 0.05. (B) Distribution of the response magnitudes of the DEGs. Statistically significant gene modulation is shown according to intensity (fold-change) and direction (up- and down-regulation).

To detect any significant alteration in biological pathways, a Fisher's exact test was performed. An enrichment analysis of the GO annotations associated with the DEGs at each sampling point (24 hpi, 24 hpi2, and 7 d) was conducted. The 30 most significantly enriched GO terms for each sampling point are shown in Figure 3. After the first injection (24 hpi), GO terms related to the immune system were found, such as those related to the regulation of innate immune response, inflammatory response, cell migration and defense response to bacteria. After the second injection (24 hpi2), genes related to the inflammatory response seemed to also be modulated and were represented in processes involved in the regulation of the NF-kB signaling pathway. Moreover, processes involved in defense response to bacteria and fungi, the negative regulation of ROS, apoptosis and glucose homeostasis appeared to be regulated after reinjection. The last sampling point (7 d) showed a modulation of genes related to GO terms involved in neural processes (long-term memory and learning), tissue regeneration (cell population proliferation and proteoglycan, glycosaminoglycan, mucopolysaccharide, and collagen metabolism) and the resolution of infections (oxidation-reduction processes and defense response to pathogens).

Figure 3

Enrichment analysis of DEGs according to the experimental design. Bars represent the proportions between the percentages of sequences in the test set (our DEGs list) and the reference set (transcriptome). (A) Biological processes overrepresented in infected mussels 24 h after the first infection. (B) Biological processes overrepresented 24 h after the second infection. (C) Biological processes overrepresented 7 d after the second infection.

The most highly expressed genes at each sampling point are shown in Table 2 (complete information in Supplementary File 3). After the first exposure, several genes showed high expression values and were significantly decreased as the experiment progressed (reinjection and 7 d). This was the case for perlucin-like protein, which is directly involved in pathogen recognition, the spore cortex-lytic enzyme, which can destroy the bacterial cell wall, and the henna protein, which is important for melanization; all of these genes could play crucial roles in the killing and sequestration of invading pathogens. These genes reached very high fold-change values after the first injection, and they decreased after reinjection and exhibited their lowest values 7 days after reinjection. Other interesting genes related to recognition (lectin, neurocan, and galaxin), acute phase response (HSP70 and sacsin), antimicrobial response (apextrin) and apoptosis (caspase 3 and the GTPase IMAP family member 4) were up- or down-regulated in a balanced manner (15 up- and 10 down-regulated). However, 1 day after the second challenge, the majority of the most highly regulated genes were down-regulated (23 down- and 2 up-regulated). For example, antimicrobial peptides such as defensin MGD-1 or myticin B were not differently regulated after the first injection, but they were indeed inhibited after reinjection (FC −1,360) and 7 days after the second injection (FC −3,745) (Supplementary File 3).

Table 2

Top 25 DEGs at each sampling point. FC, fold change.

I24h vs. C24h		I24h2 vs. C24h2		I7d vs. C7d
FC	Description	FC	Description	FC	Description
−6068.31	Caspase-3	−20257.80	Vitellogenin	−3745.45	Myticin B
−5928.63	Neuronal acetylcholine receptor subunit alpha-2	−7278.01	TRPM8 channel-associated factor 2	−1031.26	Nacre apextrin-like protein 1
−5633.81	Metalloproteinase inhibitor 3	−5276.43	Latrophilin Cirl	−956.71	Neurocan core protein
−5114.18	GTPase IMAP family member 4	−5076.28	Cytochrome c oxidase subunit 3	−712.63	Putative L-cysteine desulfhydrase 1
−3240.46	Hepatic lectin	−4899.42	Phosphoenolpyruvate carboxykinase	−693.30	Glyoxylate reductase/hydroxypyruvate reductase
−3109.94	Nacre apextrin-like protein 1	−4849.28	Bacterial hemoglobin	−620.87	Complement C1q-like protein 3
−2965.72	Nephrin	−3940.14	Betaine–homocysteine S-methyltransferase 1	−607.11	Metalloproteinase inhibitor 3
−2778.00	Neurocan core protein	−3754.37	Lysosome membrane protein 2	−596.49	Phosphate carrier protein
−2611.24	Sacsin	−3724.18	Multi-CRP-I 3	−562.79	Glyoxylate reductase/hydroxypyruvate reductase
−2395.37	Hemicentin-1	−3590.02	40S ribosomal protein SA	−510.31	Neuronal acetylcholine receptor subunit alpha-2
2410.00	TLD domain-containing protein 1	−3381.63	Venom allergen 5.01	−407.37	Probable cysteine protease RD21B
2436.01	Polyubiquitin	−3007.79	Vitellogenin-2	−246.99	Multi-CRP-I 3
2489.15	Hemicentin-1	−2451.55	Metalloproteinase inhibitor 3	−203.99	Hepatic lectin
2502.58	MAM and LDL-receptor class A domain-containing protein	−2405.10	Nacre apextrin-like protein 1	−198.48	Tubulin alpha chain
2725.96	D-arabinono-1,4-lactone oxidase	−2107.39	Myticin B	165.42	Sushi, nidogen and EGF-like domain-containing protein 1
3084.32	2′-5′-oligoadenylate synthase 1A	−1566.39	Papilin	197.54	C3a anaphylatoxin chemotactic receptor
3102.00	WAP four-disulfide core domain protein 2	−1482.80	Malate dehydrogenase, cytoplasmic	212.68	Mammalian ependymin-related protein 1
3139.69	Shell protein-5	−1466.81	Heme-binding protein 2	317.59	Proprotein convertase subtilisin/kexin type 5
3251.22	Galaxin	−1444.83	Stress-associated endoplasmic reticulum protein 2	336.24	Perlucin-like protein
3255.70	Netrin receptor DCC	−1435.92	Cytochrome c oxidase subunit 2	359.33	Peroxidasin homolog
3294.33	Nephrin	−1360.83	Defensin MGD-1	419.55	Spore cortex-lytic enzyme
4011.09	Spore cortex-lytic enzyme	−1283.57	Cytosolic 10-formyltetrahydrofolate dehydrogenase	427.42	Phenylalanine-4-hydroxylase
4027.72	Protein henna	−999.46	Glyceraldehyde-3-phosphate dehydrogenase	735.42	Protein henna
4927.82	Heat shock 70 kDa protein 12B	1357.30	Allene oxide synthase-lipoxygenase protein	915.43	Collagen alpha-2 (VIII) chain
5793.60	Low-density lipoprotein receptor-related protein 6	1451.53	Spore cortex-lytic enzyme	4103.11	Lysozyme

Changes in the Response After the Second Encounter With the Same Pathogen

Next, we looked for trainable genes with different expression values after reinjection compared to those after the first injection. We compared the transcriptomes of the challenged animals after the first and second exposures and selected primed genes, which were those with increased expression values after the second encounter with the same pathogen, or tolerized genes, which were those that showed decreased expression after reinjection (Figure 4A). A schematic representation of the expression behavior of these genes during the first and second injections with regard to naïve animals is shown in Figure 4B. Thirty-nine genes showed significantly increased expression after previous stimulation with V. splendidus, and 31 showed decreased expression. The expression levels of all these genes did not change in control animals after the first or second stimulation (Figure 4C).

Figure 4

(A) Representative scheme of the comparisons made to select the primed and tolerized genes. The statistical parameters used as thresholds were an FDR < 0.05 and a fold change (FC) > |2| in the case of Vibrio injected individuals and a FC < |2| in control mussels. (B) Mean of the tendency of the primed and tolerized genes. (C) Fold change of each primed and tolerized gene and their respective controls. The mean and standard error of all represented genes is also shown.

The list of the annotated primed and tolerized DEGs is shown in Table 3. The regulation of these genes suggests an attenuation of inflammation, a decrease in radical oxygen species (ROS) production and the inhibition of apoptosis in the second contact with V. splendidus. Primed genes with increased expression in the second encounter, such as those encoding myomodulin neuropeptides, furin-like protease (KPC-1), and plasminogen/apolipoprotein(a), were directly or indirectly related to the control and inhibition of inflammatory processes. Moreover, there was high expression of genes involved in the inhibition of ROS, such as mitochondrial uncoupling protein 2-like (UCP2), oxidative stress-induced growth inhibitor 1 (OKL38), and NAD kinase. Finally, genes involved in the reduction of cell death (inhibitor of p53-induced apoptosis-beta) and DNA repair (DNA repair protein complementing XP-A cells) were also present in our set of primed genes. However, tolerized genes with a reduced expression level after the second encounter were associated with the activation of apoptosis and inflammatory response and included regulator of nonsense transcripts 1 (UPF1), nephrin, H/ACA ribonucleoprotein complex subunit 4 (DISKERIN), nuclear migration protein (nudC), and the dual serine/threonine and tyrosine protein kinase (RIP5).

Table 3

Identified primed and tolerized DEGs.

Contig name	FC	Description
PRIMED DEGs 24h2 VS. 24h
Mg_contig_5596	423.00	Myomodulin neuropeptides
Mg_contig_14406	230.36	Mitochondrial uncoupling protein 2-like: UCP2
Mg_contig_1335	193.45	Tetraspanin-7
Mg_contig_9808	155.40	Mitochondrial uncoupling protein 2-like: UCP2
Mg_contig_4328	142.49	Furin-like protease: KPC-1
Mg_contig_30535	136.05	Tubulin beta chain
Mg_contig_14134	128.00	Plasminogen
Mg_contig_3665	57.18	Cytochrome P450 3A24
Mg_contig_50120	34.93	DNA repair protein complementing XP-A cells: XPA
Mg_contig_24689	30.88	Solute carrier family 12 member 8
Mg_contig_16026	30.15	Ropporin-1-like protein
Mg_contig_21973	27.24	Inhibitor of p53-induced apoptosis-beta
Mg_contig_13532	26.93	Zinc finger protein Eos
Mg_contig_3349	25.19	Oxidative stress-induced growth inhibitor 1: OKL38
Mg_contig_13357	23.18	Apolipoprotein(a)
Mg_contig_41851	22.52	Tetratricopeptide repeat protein 38
Mg_contig_36957	22.50	Hydrocephalus-inducing protein
Mg_contig_32198	22.28	Solute carrier family 12 member 8
Mg_contig_45429	18.20	Alpha-L-fucosidase
Mg_contig_15455	17.40	Solute carrier family 46 member 3
Mg_contig_13725	14.90	NAD kinase
TOLERIZED DEGs 24h2 VS. 24h
Mg_contig_49351	−157.33	Usherin
Mg_contig_53728	−69.50	NFX1-type zinc finger-containing protein 1
Mg_contig_88014	−40.66	Regulator of nonsense transcripts 1: UPF1
Mg_contig_50460	−22.38	NFX1-type zinc finger-containing protein 1
Mg_contig_38773	−21.16	Nephrin
Mg_contig_1952	−19.36	H/ACA ribonucleoprotein complex subunit 4: DISKERIN
Mg_contig_2166	−18.87	Phenylalanine-tRNA ligase beta subunit
Mg_contig_39783	−17.22	Nuclear migration protein: nudC
Mg_contig_7966	−16.66	Phenylalanine-tRNA ligase alpha subunit B
Mg_contig_13367	−16.20	Dual serine/threonine and tyrosine protein kinase: RIP5
Mg_contig_14736	−13.65	Phenylalanine-tRNA ligase beta subunit

Priming Induces the Modification of Functional Hemocyte Responses

Flow cytometry was used to better understand how two consecutive injections of V. splendidus affected mussel hemocytes. Two cell populations, granulocytes (R2) and hyalinocytes (R3), were well-defined in the control group (mussels injected with FSW at both sampling points) (Figure 5A). When mussels were injected with V. splendidus, the cell population structure was altered, and it was almost impossible to establish two separate populations of granulocytes and hyalinocytes. However, if mussels were stimulated with V. splendidus and received a second injection with the same pathogen, the hemocyte population structure was restored to that found in naïve mussels and showed a distribution similar to that in the controls. Quantitatively, the numbers of granulocytes and hyalinocytes were significantly reduced in mussels injected once with the bacteria compared to control animals. However, when mussels had been previously injected and received a second bacterial challenge, both cell types reached similar values to those found in the control (Figure 5B).

Figure 5

FACS analysis of the cell population distribution of mussel hemocytes after subsequent FSW/bacteria injection. (A) Dot plot of a total count of 100,000 events for hemolymph. R2, granulocytes; R3, hyalinocytes. (B) Statistical analysis of the hemocyte population in the different experimental conditions. Asterisks show the statistical significance: *p < 0.05 and ***p < 0.0001. Subsequent stimulations are indicated as follows: F-F, FSW and FSW; F-B, FSW and bacteria; B-B, bacteria and bacteria.

To confirm the results of the transcriptomic analysis that suggested the inhibition of respiratory burst activity after the second injection, we looked closely at the expression values (TPM values) of some representative genes (OKL38 and UCP-2) in individual mussels. All challenged animals exhibited a significant increase in the expression of these two antioxidant genes after reinjection (24 hpi2) (Figure 6A). We also analyzed ROS production in hemocytes from treated mussels. Respiratory burst activity when there was no triggering molecule or that was triggered by PMA was notably decreased after the second injection (Figure 6B), supporting the results of the transcriptomic analysis. Respiratory burst activity mediated by zymosan A did not show significant differences among the three groups of mussels (FSW-FSW, FSW-bacteria, and bacteria-bacteria), which was probably due to the strong stimulating effect of zymosan A that masked the natural response (33).

Figure 6

(A) Expression values (TPM) for two representative inhibitors of oxidative stress (OKL38 and UCP-2) in the six studied individuals (C2, C3, C4, I1, I2, and I10) at two experimental time points (24 hpi and 24 hpi2). (B) Respiratory burst assay. Bars represent the mean percentage and standard error of the relative luminescence units (rlu) for the total hemocytes in eight individual mussels. Asterisks show the statistical significance: *p < 0.05; **p < 0.001; and ***p < 0.0001. Subsequent stimulations are indicated as follows: F-F, FSW and FSW; F-B, FSW and bacteria; B-B, bacteria and bacteria.

The transcriptomic analysis showed that apoptosis was another central process in the response that appeared to be strongly inhibited after the second Vibrio injection (24 hpi2); therefore, we analyzed the expression of cell death inducers such as UPF1, RIP5, and nephrin in individual mussels (Figure 7A) and confirmed the results observed in the global analysis. We next performed an experiment to analyze whether two subsequent challenges resulted in changes in cell death. The number of apoptotic granulocytes and hyalinocytes was significantly increased when mussels were injected with V. splendidus. However, if mussels were injected two consecutive times with a resting time in between, the number of apoptotic cells was similar to that detected after the first injection without any further increase in the apoptotic rate (Figure 7B).

Figure 7

(A) Expression values (TPM) for three representative apoptosis inducers (RIP5, UPF1, and nephrin) in the 6 studied individuals (C2, C3, C4, I1, I2, and I10) at two experimental time points (24 hpi and 24 hpi2). (B) Apoptosis analysis. Bars represent the mean percentage and standard error of apoptotic cells in total hemocytes, hyalinocytes and granulocytes in eight individual mussels. Asterisks show the statistical significance: *p < 0.05 and **p < 0.001. Subsequent stimulations are indicated as follows: F-F, FSW and FSW; F-B, FSW and bacteria; B-B, bacteria and bacteria.

Our transcriptomic and functional results suggest that there is a modulation of the immune response after a second encounter with the same pathogen. Primed genes are involved in the resolution of the inflammatory process and the inhibition of ROS; however, repressed transcripts are related to inflammatory reactions and oxidative stress. There is a shift toward an anti-inflammatory response that attempts to minimize the damage caused by the second encounter with V. splendidus (Figure 8).

Figure 8

Summary of the main biological processes regulated during the subsequent Vibrio challenges. Note that the trained genes shift in terms of response from pro- to anti-inflammatory.

Discussion

One of the characteristics of the innate immune system is its lack of immunological memory. However, in recent years, there is increasing evidence that innate immune cells can become reprogrammed to develop immunological memory after previous encounters with non-self-molecules (34–36). Bivalves, like all invertebrates, do not have an adaptive immune system and, due to their status as filter-feeding animals, are constantly in contact with microorganisms and environmental pollutants. Our results suggest that mussels may control the magnitude of their immune response, which allows them to deal with the continuous exposure to potential pathogens. A continuous reaction in these animals against all pathogenic and non-pathogenic microbes would potentially result in a constant state of inflammation that may be detrimental for the organism.

After the first V. splendidus exposure, there was an initial response represented by a high number of modulated genes, which decreased after the second exposure and almost returned to a basal level at the end of the experiment. The decrease in the number of differentially expressed genes after the second injection is in concordance with previous findings in Crassostrea gigas (37), in which specific protection against a viral infection was achieved after poly I:C priming and a later exposure to the pathogen did not trigger an antiviral response. This suggests that in our experimental design, the primary challenge with a pathogen triggered immune processes that could be reprogrammed afterwards.

When comparing the processes significantly enriched during the two subsequent encounters, a shift from an inflammatory to an anti-inflammatory state can be observed. After the first injection, the expected response would involve the positive regulation of the inflammatory response and the migration of hemocytes, as seen in other bivalves (10, 38). However, a strong decrease in the number of DEGs and the strict control of inflammation, which is a potentially harmful process, could be observed after the second challenge, which was possibly a consequence of the inhibition of the NF-kB signaling pathway (39). Taking into account that Vibrio splendidus was used at a sublethal dose and also that the infection was cleared before the second injection, it is unlikely that bacteria virulence factors have been responsible for the regulation of the inflammatory response. Also, we cannot know if the response after the second injection would be the same if a different pathogen was injected. These aspects should be further investigated.

Inflammation is critical in the response against infection; however, this process cannot last for a long time, and a return to a non-inflammatory state requires specific suppressor molecules (40). The expression of primed and tolerized genes in the mussel transcriptome suggests a combined response that attempts to control and limit three key immune processes: inflammation, ROS production, and apoptosis. The primed genes are involved in the attenuation of inflammation by regulating the transport of ions [myomodulin neuropeptides; (41, 42)], inhibiting NF-kB via ubiquitination [KPC-1; (43)], or inhibiting inflammatory pathways [plasminogen/apolipoprotein(a); (44–46)]. On the other hand, the tolerized genes (inhibited in the second challenge) are involved in the regulation of cell death [RIP5; (47)].

As is the case for the inflammatory process, an excess of reactive oxygen species can be harmful to the organism. At physiological levels, ROS are involved in intracellular signaling and defense, but uncontrolled production yields oxidative stress. Therefore, ROS are strictly controlled in all organisms to prevent self-inflicted damage (48, 49). The production of ROS is a well-characterized defense process in bivalves (50), but after the second exposure to the same pathogen, the mussels seemed to actively control oxidative stress by inhibiting respiratory burst activity with the expression of genes such as UCP2 (51), OKL38 (52), and NAD kinase (53). It seems that the control of oxidative stress is one of the central modulated processes after repeated encounters. Accordingly, at the functional level, we observed that the ROS levels were reduced after two subsequent V. splendidus challenges, suggesting that hemocytes could prevent an uncontrolled respiratory burst. Moreover, hemocytes seem to avoid cellular impairment caused by DNA damage resulting from previous oxidative stress (54), with the overexpression after the second injection of the DNA repair protein XPA (55).

A strong link between oxidative response and apoptosis, which is the other strongly modulated process in our results, has been reported (56, 57). The tolerized genes include several modulators of apoptosis. These genes, which have a predominantly pro-apoptotic function, are inhibited after the second exposure and include UPF1 (58), nephrin (59), DISKERIN (60), RIP5 (47), and nudC (61). At the functional level, the number of apoptotic cells was significantly increased after a single injection and was maintained, rather than increased, after the second injection.

In addition, the loss of the hemocyte population distribution after the first challenge and the restoration to normal conditions after the second challenge are in concordance with the presence of a priming process and reflect hemocyte recovery after the previous V. splendidus encounter. In our opinion, this evident change in hemocyte populations is very revealing and complements previous priming results in mollusks (62–65), and it should be further explored.

Eventually, 7 days after the second exposure, the resolution of the infection occurred. Processes related to tissue regeneration and learning are represented in our results, including the up-regulation of proliferation (C3a anaphylatoxin chemotactic receptor, skin secretory protein xP2, sushi, nidogen, and EGF-like domain-containing proteins), maintenance of the extracellular matrix (collagen, techylectin, perlucin, neurocan, and glycoproteins), and learning (phenylanaline-4-hydroxylase/henna protein). The presence of neural processes related to the generation of memory at all sampling points is remarkable. Although this might be due to the strong bias in the GO database in terms of model organisms, we cannot discard genuine evidence of a process that occurs during infection. The protein related to these GO terms, which is involved in the melanization cascade in invertebrates, could also be related to ancient cognition and behavior mechanisms, which are known to occur in invertebrates (66). In any case, it seems that certain processes, such as the metabolism of phenylalanine, could be involved in the generation of innate immune memory. Previous studies in mammals have shown the impairment of cognitive function due to phenylalanine hydroxylase deficiency (67–69), showing a possible link to learning processes and the evolutionary conservation of this mechanism.

In summary, the immune responses of M. galloprovincialis after the first and second encounter with V. splendidus were different. The analysis of the differentially expressed genes suggests that, after the second contact with the bacteria, the mussel hemocytes attempted to control and resolve the inflammatory response to avoid subsequent DNA damage and cell death. There appears a tightly regulated response shifting from a pro-inflammatory response to an anti-inflammatory and probably regenerative phenotype. In conclusion, these results indicate the existence of a secondary immune response in mussels oriented to tolerate infection by inducing anti-inflammatory processes to minimize tissue damage.

Data Availability

The datasets generated for this study can be found in the NCBI Short Read Archive database under the accession ID SRP145077.

Ethics Statement

The Mediterranean mussel, M. galloprovincialis, is not considered an endangered or protected species in any international species catalog, including the CITES list (www.cites.org), and it is not included in the list of species regulated by EC Directive 2010/63/EU. Therefore, no specific authorization is required to work on mussel samples.

Author Contributions

BN and AF conceived and designed the project. RM and MR-C performed the mussel infections, sampling, and RNA extraction. AF, AP, MG, and RM performed the trimming, assembly, annotation and the gene expression, and statistical analyses. BN, MG, RM, and MR-C analyzed the generated data. MR-C performed the functional assays and wrote the manuscript. All listed authors revised, edited, read, and approved the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.