Introduction {#Sec1}
============

In the past years, there has been a considerable effort to understand the molecular basis of many biological processes of non-model organisms, such as bivalves. Although they are cultured worldwide and have an important ecological value, we are still far from understanding how these animals respond to pathogens. Mussels (*Mytilus galloprovincialis*) are present on all the continents, and they are even considered as a powerful invasive species^[@CR1],[@CR2]^. Their ubiquitous presence and their sessile life explain why mussels have been selected and used as ecological markers for pollution for more than 40 years (mussel watch project)^[@CR3]^. One of the most interesting characteristics of these invertebrates is the filtering activity that they have to use to feed themselves. On average, mussels can filter approximately 7.5 liters of water in one hour^[@CR4]^, meaning they are in intimate contact with millions of microorganisms that are potentially pathogenic to them^[@CR5]^. In fact, Stabili *et al*.^[@CR6]^ reported that the abundance of *Vibrio spp*. is higher in mussels than in the surrounding water. Bivalves are susceptible to numerous diseases that compromise their culture, producing economic losses all over the world. As an example, the ostreid herpesvirus 1 (OsHV-1) has caused massive mortalities in oysters (*Crassostrea gigas*) in different parts of the world^[@CR7],[@CR8]^. Mussels that cohabitate in the same areas as oysters and clams have experienced massive mortalities due to viral and bacterial infections. However, high mortalities have never been reported in the field for *M*. *galloprovincialis*^[@CR9],[@CR10]^.

Mussels, as all invertebrates, have only an innate immune system to fight against pathogenic microorganisms. It is well known that invertebrates recognize conserved pathogen-associated molecular patterns (PAMPs)^[@CR11],[@CR12]^, danger signals, and danger/damage-associated molecular patterns (DAMPs)^[@CR13],[@CR14]^, yet it is not clear how bivalves, and mussels in particular, react to external stimuli taking into account their constant contact with microorganisms. To understand the immune reaction of these animals, a considerable effort has been made in recent years to increase the genomic and transcriptomic resources in bivalves and specifically in the Mediterranean mussel, *M*. *galloprovincialis*^[@CR15]--[@CR20]^.

To date, all the genomic studies carried out in mussels have been conducted using biological pooled samples and have not analyzed the response of individuals. However, we know that mussels have immune effectors that are tremendously variable within the population^[@CR21]--[@CR23]^. This is the case of myticins, which are antimicrobial peptides highly expressed in mussels^[@CR24]^. There are three different myticin genes, named A, B and C, with similar DNA sequences and physicochemical properties. Of them, myticin C is the most variable AMP in mussels with the broader biological properties^[@CR25]^.

In this work, we analyzed the transcriptomic response of six individual naïve mussels and how they respond to a simple injury (control animals injected with filtered sea water) or against a bacterial challenge (infected animals with *Vibrio splendidus*). Our results highlight the variability of the response in individual mussels and the importance of the appropriate experimental controls since the injection on the adductor muscle (the usual method to experimentally infect mussels) might be interpreted by the organism as a danger signal that may influence the unexpected regulation of several expressed genes.

Methods {#Sec2}
=======

Animals {#Sec3}
-------

Adult *M*. *galloprovincialis*, 8--10 cm in shell 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 at least one week.

Experimental design {#Sec4}
-------------------

A schematic representation of the experimental design is shown in Fig. [1](#Fig1){ref-type="fig"}. Twenty naïve mussels were marked and notched on 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. This was the time zero (t0) sampling. The hemolymph 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 with syringe and a 0.5-mm-diameter (25 G) disposable needle and stored at −80 °C until RNA isolation.Figure 1Diagram of the experimental design followed by the stimulation of individual mussels. Ten individual mussels were used per condition (control or infection), but only the three selected mussels for sequencing, in each condition, were represented in the experimental design.

After one week, 10 of the 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 *Vibrio splendidus* (reference strain, LGP32) at a non-lethal concentration (1 × 10^7^ CFU/mL). One day after the challenge, 24 hours post injection (24 hpi), hemolymph (500 µl) was sampled again from individual mussels and centrifuged in the same manner described above, and the pellet was resuspended in 500 µl of TRIzol. Samples were immediately homogenized with syringe and 25 G needle and kept at −80 °C until RNA isolation.

RNA isolation, cDNA production and Illumina sequencing {#Sec5}
------------------------------------------------------

RNA isolation was carried out in the 40 samples (n = 20 naïve at t0, n = 10 FSW injected at 24hpi, and n = 10 bacteria injected at 24hpi) using TRIzol and following the manufacturer's protocol. Purification of RNA after DNase I treatment was performed with RNeasy mini (Qiagen). Next, the concentration and purity of the RNA was measured using a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Inc.), and RNA integrity was tested on an Agilent 2100 Bioanalyzer (Agilent Technologies) before producing cDNA libraries for Illumina sequencing. Only the individuals with the best RNA samples (in terms of RNA quantity and quality) from both sampling points were chosen for Illumina sequencing: control n° 2 (C2), control n° 3 (C3) control n° 4 (C4), infected n° 1 (I1), infected n°2 (I2) and infected n° 10 (I10). In total, 12 RNA samples (2 per individual, the first at t0 and the second 24hpi of FSW or bacteria) were sequenced (details in Table [1](#Tab1){ref-type="table"}).Table 1Summary of the transcriptome bioinformatics pipeline.Reads originrawtrimmedC2 t078,426,94899.59%C3 t044,346,85498.03%C4 t0100,814,19899.59%I1t017,696,89498.08%I2 t093,114,09899.60%I10 t096,780,60299.64%C2 24 h95,296,48499.49%C3 24 h51,708,98897.62%C4 24 h92,661,28299.65%I124 h52,102,30298.77%I2 24 h90,965,87599.52%I10 24h99,262,97099.55%**Assembly**Contigs270,324Range contig length200--15,624Average contig length512N50574**Blast**Contigs identified by Uniprot/SwissProt24.97%Contigs identified by molluscs database54.93%**GO analysis**Annotated contigs24.87%**KEGG analysis**Pathway assigned contigs8.03%

The mRNA-Seq sample preparation kit from Illumina was used according to the manufacturer's instructions. mRNA was extracted from total RNA using oligo (dT) magnetic beads and cleaved into short fragments using fragmentation buffer. A cDNA library compatible with the Illumina NGS technology was then prepared from the fragmented mRNA via reverse transcription, second-strand synthesis and 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 Qbit system. Next-generation sequencing was performed using Illumina HiSeq™ 4000 technology at Macrogen Inc. Korea (Seoul, Republic of Korea).

Bioinformatics and RNA-Seq {#Sec6}
--------------------------

CLC Genomics Workbench, v.10.0.1 (CLC Bio; Qiagen) was used to filter, assemble and perform the RNA-Seq and the statistical analyses of individual mussels.

Raw reads were trimmed to remove low quality sequences (quality score limit 0.05 = PHRED 13), adaptor sequences, and sequences shorter than 70 bp. Then, a reference global transcriptome of the six mussels was assembled with an overlap criterion of 70% and a similarity of 90% to exclude paralogous sequence variants. The settings used were a mismatch cost = 2, deletion cost = 3, insert cost = 3, and minimum contig length = 200 base pairs.

Before the expression analysis, a subsampling step of the trimmed reads was carried out. A random subset of sequences was generated to equal the number of reads present in each sample. Then, RNA-Seq analysis of the subsamples (mismatches = 2, length fraction = 0.8, similarity fraction = 0.8, and maximum hits per read = 10) were performed. The expression values were set as transcripts per million (TPM). Finally, a differential expression analysis test (a Robinson and Smyth's Exact Test, which assumes a Negative Binomial distribution of the data and takes into account the overdispersion caused by biological variability) was used to compare expression levels in each sample and to find the differentially expressed genes (DEGs). Transcripts with absolute fold change (FC) values \>2 and Bonferroni corrected p-value \< 0.05 were retained for further analyses.

BLAST annotation, GO assignment, enrichment and KEGG analysis {#Sec7}
-------------------------------------------------------------

Uniprot/Swissprot BLASTx results were used to obtain the Gene Ontology (GO) term assignments of the contig list with the Blast2GO software^[@CR26]^. To improve the percentage of identification from the contig list, it was also annotated with an inhouse built database made with all of the mollusc sequences present in the NBCI nucleotide database. In both blast approaches, the e-value threshold was set at 1e-3. Then, the enrichment analyses of the up and down regulated DEGs (test set) were conducted, including the global mussel transcriptome as the reference set. A Fisher's exact test was run with default values and a p-value cut-off of 0.05. The option to reduce the enriched list to the most specific GO terms was used. Only over-represented biological process (BP) terms were further analyzed. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in which the DEGs were involved were also analyzed with Blast2GO and summarized following the existent categories in the KEGG database (<http://www.genome.jp/kegg/pathway.html>).

Flow cytometry analysis {#Sec8}
-----------------------

The distribution of myticin C in different cell populations was assayed in mussel hemocytes by immunocytochemistry and flow cytometry (FACSCalibur; BD). The experimental design is depicted in Fig. [2](#Fig2){ref-type="fig"}. Briefly, the studied conditions were five: (1) naïve mussels, (2) mussels injected with FSW, (3) mussels injected with *V*. *splendidus*, (4) mussels bled, let rest for 1 week and injected with FSW and (5) mussels bled, let rest for 1 week and injected with *V*. *splendidus*. A single end-point sampling was performed for each condition. Four mussels per condition were analyzed.Figure 2Experimental design for immunocytochemistry FACS detection of myticin C protein in mussel hemocytes. Four individual mussels were used for each experimental condition (naïve, non-bled FSW, non-bled Vibrio, bled FSW and bled Vibrio).

For each condition, 1.5 mL of hemolymph from individual mussels was withdrawn and immediately fixed in a final concentration of 2% paraformaldehyde (PFA). Samples were divided into three aliquots; the first sample was to analyze myticin C production through a custom antibody^[@CR25]^, a control in which prebleed serum was used instead of the primary antibody and an absolute control without primary antibody. The immunocytochemistry protocol was performed as previously described^[@CR25]^. The cell suspension was fixed for 15 min at 4 °C and washed twice for 10 min in phosphate buffered saline buffer (PBS) prior to permeabilization and staining with the anti-myticin C antibody (1:100). After overnight staining, samples were washed once for 10 min (PBS, 0.1% saponin, 0.2% bovine serum albumin, BSA) and incubated in the dark for 50 min at room temperature with secondary antibody (1:500). Finally, samples were washed for 10 min (PBS), and 200 µl of each condition and replicate was dispensed in a 96-well plate to be analyzed by fluorescence-activated cell sorting (FACS). The density plots and histograms were generated using Cell Quest Pro software (BD), and a one-way ANOVA with post-hoc analysis was used to analyze the significance of the results among all the conditions tested.

Validation of RNA-Seq expression in response to tissue injury {#Sec9}
-------------------------------------------------------------

A new experiment, including a new non-injected control, was carried out following the same layout as the RNA-Seq experimental design with new biological samples (Supplementary Fig. [1A](#MOESM1){ref-type="media"}). Briefly, individual mussels were sampled at t0, let them rest for one week and then divided in three tanks: (1) mussels injected with *V*. *splendidus* (1 × 10^7^ CFU/mL), (2) mussels injected with FSW and (3) mussels remained non- injected. One day after the challenge hemolymph was sampled again from individual mussels and RNA was extracted as previously described. cDNA was synthesized from each individual mussel with 200 ng of total RNA using NZY First-Strand cDNA Synthesis Kit (nzytech) following the manufacturer's protocol.

Gene expression of selected genes (Supplementary Fig. [1B](#MOESM1){ref-type="media"}) was analyzed in a 7300 Real Time PCR System (Applied Biosystems). One microliter of fivefold-diluted cDNA template was mixed with 0.5 ml of each primer (10 mM) and 12.5 ml of SYBR Green PCR master mix (Applied Biosystems) in a final volume of 25 ml. 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 18 S as a reference gene following the Pfaffl method and standardized to the normalized expression of the t0 samplig to calculate fold changes. An independent t-test was used to analyze differences among conditions and differences were considered significant with p-value \< 0.05. For the validation of the RNA-Seq vs the qPCR results linear regression and correlation were performed to analyze the studied genes and conditions.

Results {#Sec10}
=======

Assembly and annotation of mussel transcriptome {#Sec11}
-----------------------------------------------

A summary of the sequence origin, assembly, identification, and annotation results is shown in Table [1](#Tab1){ref-type="table"}. An average of 76 million raw reads was obtained from each individual sample of *M*. *galloprovincialis* hemocytes. The CLC Genomics Workbench was used to filter the raw reads, and over 97% of raw reads successfully passed the quality control in all of the samples. The assembly step was performed with all the samples available to obtain a global mussel transcriptome; 270,324 contigs were assembled with an average length of 512 bp. The putative identities of these sequences were obtained by Blast by two different means; Blast2GO software was used to identify the 24.97% of the contigs through a BLASTx approach against Uniprot, and CLC was used to identify 99.94% of the contigs using an inhouse designed database with all the sequences available in NCBI for molluscs. GO terms were assigned to 24.87% of the contigs and enzyme codes to find KEGG pathways to 8.03% of the sequences.

Mussel transcriptome after bacterial or DAMP stimulation {#Sec12}
--------------------------------------------------------

The experimental design allowed us to sample hemolymph from each individual mussel before and after injection with bacteria or FSW; therefore, the real behavior of the modulated genes could be followed in each animal. Figure [3](#Fig3){ref-type="fig"} shows the distribution of the differentially expressed genes (DEGs) in control and infected animals 24 hpi with regard to their own t0 sampling point. An expected response with more DEGs in infected animals (I1: 3,900 DEGs, I2: 2,286 DEGs, I10: 2,514 DEGs) compared to control animals (C2: 1,562 DEGs, C3: 486 DEGs, C4: 751 DEGs) was found. To validate the results of the RNA-Seq, the expression of IRG1, IFI44, MyD88, SOCS2, Myticin C, C1q and apextrin have been studied by qPCR. The statistical approach showed a significant correlation between the fold changes obtained with the two methodologies, validating our results (Supplementary Fig. [1C](#MOESM1){ref-type="media"}).Figure 3General overview of the expression values of mussels after bacterial challenge (infected) or FSW injection (controls). Stacked column chart reflecting the fold change distribution of DEGs in control and infected animals 24 hpi with regard to their own t0 sampling point.

We analyzed the effect of the injection in each animal and, surprisingly, a very small percentage of the modulated genes were shared between the 3 controls (only 2.9%) and the 3 infected animals (only 9.4%) as shown in the Venn diagrams in Fig. [4a,b](#Fig4){ref-type="fig"}. This result could mean that each mussel responds to the same stimulus in a different way. A simple injection affected the gene expression of control mussels in which a mean of 933 DEGs were significantly modulated after this danger signal. It has to be taken into account that a non-injected control after one-week rest is lacking in this experimental design for Illumina sequencing and RNA-Seq analysis. However, it was indeed included in a complementary experimental design and the gene expression of some immune-genes was analyzed by qPCR (Supplementary Fig. [1D](#MOESM1){ref-type="media"}). Results showed no significant differences between the initial and end-point sampling in non-injected controls except for C1q.Figure 4Venn diagrams of DEGs for each condition. (**a**) in blue, control mussels. (**b**) in red, infected mussels. (**c**) in gray, common genes to control and infected mussels.

Of the total genes, 89.5% was modulated by the infection, and 8.8% was regulated after the tissue injury without bacterial challenge. When the common transcripts shared by the 3 controls and the 3 infected animals were compared (Fig. [4c](#Fig4){ref-type="fig"}), there were only 11 common contigs. Of them, five had informative annotations: C1q domain containing protein MgC1q61 (related to recognition in molluscs), chromobox protein homolog 7 (with epigenetic functions), putative gastrointestinal growth factor xP4 (important in mucosal protection and reconstitution), transcription factor AP-1 (a key transcription factor for immunity) and heat shock protein beta-1 (with numerous biological roles including regulation of stress resistance, inflammation and apoptosis). All these genes, regulated both in control and infected mussels, are directly or indirectly related to the defense response.

The enrichment analysis of the GO terms in each compartment of the Venn diagram showed in Fig. [4c](#Fig4){ref-type="fig"} revealed 100 biological processes (BP) over-represented in the exclusive DEGs for controls, 439 for infected mussels and 37 in the common DEGs list for control and infected mussels. Figure [5](#Fig5){ref-type="fig"} represents the top 30 most significant BP in these groups, divided into up or down regulated according to the fold change of the DEGs that they represent. The enrichment analysis revealed that the BP differentially represented in the common DEGs between control and infected mussels were related to immune terms such as "positive regulation of endothelial cell chemotaxis" representing the heat shock protein 27, "positive regulation of monocyte differentiation" or "negative regulation by host of viral transcription", both terms related to the transcription factor AP-1. GO terms related to histone or DNA methylation were also present: "histone H4-K20 trimethylation" and "DNA hypermethylation", representing the chromobox protein homolog 7. This result could indicate that both control and infected animals were subjected to a stimulus with similar capability to trigger epigenetic changes to respond to damage or infection. In the group of exclusive DEGs of control mussels, unexpected results were found, with BPs related to the immune response, antimicrobial peptides, and more specifically, to myticins, such as "cell killing". Also in the control group, there were more processes related to the immune response both in the up and down regulated genes such as, in the up-regulated group, "regulation of lymphocyte mediated immunity", "protection from natural killer cell mediated cytotoxicity" both terms related to serpin A3C, and, in the down-regulated group, "response to external stimulus", "regulation of programmed cell death" terms representing genes such as the cell death specification protein 2, the GIMAP 4 or the superoxide dismutase. In the group of DEGs exclusive to infected mussels, the most interesting processes found were the ones linked to response to hypoxia, glucocorticoid, lipopolysaccharide and those related to cell proliferation, cell cycle or leucocyte activation, an indication of active hemocyte division and proliferation, representing 70 genes such as the suppressor of cytokine signaling 2, the ghrelin receptor, the fos-related antigen 2, the ficolin 2, the fibrinogen C domain-containing protein 1 and several growth and transcription factors.Figure 5Enrichment analysis of DEGs. The proportions of test/reference sequences for up and down regulated genes are described following this legend: (**a)** blue, exclusive genes for control mussels; (**B**) red, exclusive genes for infected mussels; (**C**) gray, common genes to control and infected mussels.

The KEGG reference pathway analysis showed, in the absolute number of DEGs, that infected mussels always had the maximum representation in every pathway. As we wanted to represent these results proportionally, we calculated the percentage of the DEGs ascribed to every pathway with regard to the modulated genes in each Venn category. This analysis is summarized in Fig. [6](#Fig6){ref-type="fig"}, which shows that metabolism was affected in different ways in control and infected animals. The pathways with more representation in control mussels were related to nucleotide and carbohydrate metabolism and to the metabolism of cofactors and vitamins. In infected mussels, the great majority of DEGs were included in metabolism pathways known to be greatly affected by immune challenges and in pathways related to translation, signal transduction (phosphatidylinositol signaling system, the mTOR signaling pathway and the PI3K-Akt signaling pathway), viral infectious diseases and immune system (defense cells differentiation and signaling) (Supplementary File [1](#MOESM1){ref-type="media"}). The PI3K, Akt and mTOR pathways are linked to the JAK-STAT signaling pathway, and their regulation is intimately related to the immune response and the regulation of processes such as cell proliferation, autophagy and apoptosis^[@CR27]^, in line with the enrichment analysis results for infected mussels.Figure 6Summary of the KEGG reference pathway results for the significantly regulated contigs before and after injection: common for control mussels (blue) and common for infected mussels (red).

With the aim of knowing more about the nature of these DEGs, we focused on the most expressed genes in each group (complete data in Supplementary Table [1](#MOESM1){ref-type="media"}). Tables [2](#Tab2){ref-type="table"}, [3](#Tab3){ref-type="table"} and [4](#Tab4){ref-type="table"} show, respectively, the top 25 genes found exclusively in controls, in infected animals and all those regulated genes shared between both of them, with the fold changes for each individual mussel. It is worth mentioning that antimicrobial peptides (AMPs) such as myticins or myticusin were highly up regulated in FSW injected animals (Table [2](#Tab2){ref-type="table"}) and not in the infected animals, most likely reflecting a reaction against an injury or a danger signal without pathogens. Interestingly, in the DEGs found only in controls, genes related to cell proliferation, differentiation or cell activation were also up regulated and include elastin microfibril interface-located protein 2 (EMILIN-2), stathmin, low affinity epsilon Fc receptor, cell division cycle-associated protein 2 (CDCA2), and signal transducer and transcription activator (STAT). Some of these are associated with the cytoskeleton and cell motility (for example stathmin or the myticins that increase chemotaxis). The most downregulated genes in controls were related to the regulation of cell death: cell death specification protein 2, immediate early response gene 5 protein, protocadherin Fat 1, and GTPase IMAP family member 4.Table 2Top 25 up-regulated DEGs associated to tissue injury (controls). FC: fold change.Contig IDDescriptionFC C2FC C3FC C4MeanMean p-valMg_109517Myticin C30.242.215.2512.57\<0.0001Mg_7618Myticin B25.002.092.859.98\<0.0001Mg_127788Myticin C21.312.672.408.79\<0.0001Mg_1042Elastin microfibril interface-located protein 220.092.373.268.57\<0.0001Mg_33725Myticin C15.242.992.086.77\<0.0001Mg_5071Stathmin9.225.353.996.190.0003Mg_7746Myticusin-alpha precursor6.547.314.206.02\<0.0001Mg_150217Myticusin-alpha precursor6.406.714.225.78\<0.0001Mg_37432Heavy metal-binding protein HIP8.393.624.205.40\<0.0001Mg_26476Low affinity immunoglobulin epsilon Fc receptor7.493.562.714.59\<0.0001Mg_27764Cell division cycle-associated protein 26.834.222.564.540.0001Mg_3246Serine protease inhibitor A3C7.603.792.154.51\<0.0001Mg_5822Signal transducer and transcription activator3.892.452.983.110.0001Mg_4540Signal transducer and activator of transcription 24.122.432.332.960.0056Mg_31595Structural maintenance of chromosomes protein 43.062.892.822.930.0018Mg_1195Neuronal PAS domain-containing protein 4−3.25−3.09−3.92−3.42\<0.0001Mg_1629Sterol regulatory element-binding protein 1−3.81−5.01−2.54−3.79\<0.0001Mg_23409Perlucin−5.22−3.84−5.89−4.98\<0.0001Mg_20788GTPase IMAP family member 4−6.73−2.92−6.03−5.23\<0.0001Mg_157047Protocadherin Fat 1−5.63−5.08−5.87−5.52\<0.0001Mg_29303Notch homolog 2 N-terminal-like protein−5.48−5.01−6.44−5.64\<0.0001Mg_6692Immediate early response gene 5 protein−2.45−12.17−10.66−8.43\<0.0001Mg_10957Cell death specification protein 2−3.92−13.93−8.16−8.67\<0.0001Mg_23430Cornifelin homolog B−15.91−3.61−34.40−17.970.0002Mg_247Putative ariadne-like RING finger protein R811−89.09−3.30−27.59−39.99\<0.0001Table 3Top 25 up-regulated DEGs associated exclusively to a *Vibrio* infection. FC: fold change.Contig IDDescriptionFC I1FC I2FC I10MeanMean p-valMg_80300DNA-directed RNA polymerase II subunit RPB1-like441.8589.211085.67538.91\<0.0001Mg_106Cis-aconitate decarboxylase655.70221.92370.69416.10\<0.0001Mg_457Solute carrier family 12 member 825.4166.52563.75218.56\<0.0001Mg_599Neuropeptide Y receptor type 2-like276.9389.07288.15218.05\<0.0001Mg_261923Tropomyosin-126.4469.72153.2083.12\<0.0001Mg_3149Neuropeptide Y receptor type 2-like72.7711.07144.7076.18\<0.0001Mg_1646PDZ and LIM domain protein Zasp50.97117.8432.4467.08\<0.0001Mg_784Furin-like protease kpc-162.9432.0992.9162.65\<0.0001Mg_24087Suppressor of cytokine signaling 2164.0917.374.6962.05\<0.0001Mg_32754Zinc finger protein 2679.1048.1529.3652.20\<0.0001Mg_19011Phenylalanine--tRNA ligase alpha subunit42.2347.9358.5449.57\<0.0001Mg_2877Growth hormone secretagogue receptor type 192.4538.466.3145.740.0003Mg_9258Suppressor of cytokine signaling 2103.6312.172.8539.550.0078Mg_7966Phenylalanine--tRNA ligase alpha subunit B43.1128.6643.7938.52\<0.0001Mg_6277Interferon-induced protein 44-like22.7627.7857.7336.09\<0.0001Mg_3834Myeloid differentiation primary response protein MyD8866.3523.4817.2435.69\<0.0001Mg_10058Phenylalanine--tRNA ligase alpha subunit37.2323.1646.0135.46\<0.0001Mg_6296Cyclic GMP-AMP synthase64.9417.8312.4531.74\<0.0001Mg_17428Complement C1q-like protein 4−40.12−25.04−40.69−35.28\<0.0001Mg_4119Collagen alpha-1(XIV) chain−53.14−30.65−28.95−37.58\<0.0001Mg_21146Heat shock protein−53.09−47.87−25.07−42.01\<0.0001Mg_44317Nacre apextrin-like protein 1−84.47−8.82−39.92−44.400.0140Mg_3338Headcase protein−156.82−16.21−19.39−64.14\<0.0001Mg_11571Protein NDRG1−59.97−87.36−75.87−74.40\<0.0001Mg_24372Complement C1q tumor necrosis factor-related protein 4−214.51−18.16−19.81−84.16\<0.0001Table 4Regulated DEGs commons to controls and infected animals. FC: fold change.Contig IDDescriptionFC C2FC C3FC C4Mean CFC I1FC I2FC I10Mean II vs C p-valueSifnificanceMg_81902C1q domain containing protein MgC1q617.968.546.857.78−2.37−3.77−2.47−2.87\<0.0001\*\*\*\*Mg_257Putative gastrointestinal growth factor xP43.672.282.272.74−7.42−2.16−5.26−4.950.0085\*\*Mg_384Chromobox protein homolog 7−2.34−2.14−3.09−2.52−2.84−2.40−3.48−2.910.4196nsMg_234Transcription factor AP-1−5.12−4.36−3.55−4.34−15.56−8.14−7.55−10.420.0811nsMg_496Heat shock protein beta-1−290.15−25.12−52.75−122.67−82.31−94.11−80.22−85.550.6821nsP-value and significance show the results of an independent t-test between the controls and the infected mussels.

In infected animals (Table [3](#Tab3){ref-type="table"}), some of the exclusive DEGs were directly involved in the immune response, in concordance with the results of the KEGG pathways analysis (Fig. [6](#Fig6){ref-type="fig"}). Some of these genes that were always up regulated were IRG1, SOCS2, IFI44 and Myd88. The cis-aconitate decarboxylase (also known as Immune-responsive gene 1 protein or IRG1) is an enzyme of de tricarboxylic acid cycle which produces itaconate after an immune challenge. Itaconate has been recently studied for its anti-inflammatory and antimicrobial properties^[@CR28],[@CR29]^, making a direct link between basic metabolism and the immune response. The suppressor of cytokine signaling 2 (SOCS2) is a well-known negative regulator of the JAK-STAT signaling cascade, whose function is to control the inflammatory response^[@CR30]^. The myeloid differentiation primary response protein Myd88, an essential mediator of the Toll signaling pathway, has been characterized in bivalves, and it is up regulated after bacterial infections^[@CR31]^, consistent with our results. The interferon-induced protein 44 (IFI44), on the other hand, has been less well studied, but its upregulation has been associated with viral infections also in bivalves^[@CR32]--[@CR34]^, and it is known to have antiviral properties^[@CR35]^. However, its exact function is still unknown, and it could be possible that in bivalves, it also has antibacterial roles, such as its strong upregulation 24 hours after *V*. *splendidus* injection suggested. The presence of IFI44 in the infected DEGs list is a possible explanation for the "viral infectious diseases" category in the KEGG analysis.

Table [4](#Tab4){ref-type="table"} shows the regulated transcripts shared between controls and infected animals. Two of them (C1q domain containing protein MgC1q61 and Putative gastrointestinal growth factor xP4) were up regulated in all the controls but down regulated in infected mussels, and their gene expression was significantly different between control and infected mussels. Individual responses can override experimental design in many genes, but in this case, even despite individual variability, it was possible to detect experimentally induced changes in hemocyte gene expression. In the specific case of the C1q domain containing protein, it could be due to its role as a pathogen recognition protein and its fast upregulation after an infection, with a return to physiological levels within 24 h^[@CR23]^. The putative gastrointestinal growth factor xP4 is a member of the trefoil factor family (TFF), a group of molecules with a pivotal role in maintaining the surface integrity of mucous epithelia *in vivo*, explaining its up-regulation after the injury in control mussels; although the exact role of the TFF peptides is unknown they are known to have antiapoptotic effects and probably modulate inflammatory processes^[@CR36],[@CR37]^ a possible reason of its diverse behavior in control and infected samples. The other 3 annotated genes, always down regulated, were heat shock protein beta-1, associated with the acute phase response and found modulated after a tissue injury in *M*. *galloprovincialis* at the protein level^[@CR38]^; chromobox protein homolog 7, related to epigenetics via chromatin remodeling and modification of histones^[@CR39]^; and the transcription factor AP-1, known to control the expression of genes related to differentiation, proliferation and apoptosis, and intimately linked to NF-kB^[@CR40]^.

Expression of myticins {#Sec13}
----------------------

The presence of myticins in the group of the most significantly up regulated genes after FSW injection led us to further investigate this fact. Taking into account that we had previously found that these molecules had antibacterial and antiviral properties^[@CR25],[@CR41],[@CR42]^, we expected to find them being up regulated after bacterial challenge and not after a tissue injury, as they were in the controls. Figure [7](#Fig7){ref-type="fig"} shows the expression of all the differentially expressed myticins in each individual. Control animals, after the injection of FSW, presented higher expression values than those challenged with *V*. *splendidus*, where a down regulation of myticins could be observed for only one (I2) or none (I10) of the myticin transcripts that were differentially expressed.Figure 7Expression profiles of myticins in control and infected individual mussels. Dots indicate the fold change (FC) of myticin transcripts 24 hpi with regard to their own t0. The line indicates the mean of the expression for each individual. Note the increase in the expression level in controls with regard to infected mussels and the great individual variability in the number of regulated transcripts, especially among infected mussels, with I10 showing no differentially expressed myticins.

To further confirm these unexpected results, we designed a new experimental protocol to determine if myticins had the same behavior at the protein level using flow cytometry. The experiment had the same layout as the one used to study the transcriptome but with more experimental controls: we included naïve and non-bled mussels to confirm the involvement of tissue injury in the expression of myticins (Fig. [2](#Fig2){ref-type="fig"}). Figure [8a](#Fig8){ref-type="fig"} shows the general profile of the hemocyte population, with R1 being the granulocytes and R2 the hyalinocytes. After *V*. *splendidus* infection in bled (injured) and non-bled mussels, there was a change in the hemocyte population structure with a significant reduction in the percentage of granulocytes with regard to their controls (FSW) and a significant increase in the percentage of hyalinocytes. These changes were present in bled and non-bled animals after infection, but the decrease in the granulocytes was less important in bled or injured mussels (Fig. [8b](#Fig8){ref-type="fig"} and Supplementary Table [3](#MOESM1){ref-type="media"}). It is noteworthy that the hemocyte population changed drastically after bacterial injection, but it recovered significantly, with regard to the non-bled mussels, if a previous damage stimulus had been made (Fig. [8b](#Fig8){ref-type="fig"} and Supplementary Fig. [2](#MOESM1){ref-type="media"}).Figure 8Description of hemocyte populations by flow cytometry. (**a**) Dot plot of the representative cell populations in mussel: R1, granulocytes; R2, hyalinocytes. Represented dots are gated for R6 region. (**b**) Description of the variation in hemocyte populations among samples: Orange, granulocytes; brown, hyalinocytes. Standard deviation is indicated for the four replicates. Different letters indicate significant differences (p-value \< 0.05) among groups: low-case letters for granulocytes and capital letters for hyalinocytes.

In Fig. [9a](#Fig9){ref-type="fig"}, the region considered positive for myticin C expression is shown. These positive cells corresponded to the granulocyte (R1) population (Fig. [9b](#Fig9){ref-type="fig"}). When the production of myticin C was evaluated after the different treatments, a significant decrease in the population positive for myticin C was found after *Vibrio* infection (Fig. [9c](#Fig9){ref-type="fig"} and Supplementary Table [3](#MOESM1){ref-type="media"}). However, this reduction was not observed when mussels were bled before, generating a tissue injury. The percentage of myticin C positive granulocytes was also significantly higher in mussels previously bled and then infected with bacteria when compared with mussels only injected with bacteria. Figure [9d-f](#Fig9){ref-type="fig"} illustrate representative R1 FL1 profiles for every sample group. Naïve or control mussels injected with FSW showed a granulocyte population positive for myticin C (Fig. [9d](#Fig9){ref-type="fig"}) that was lost in *Vibrio* challenged mussels without a previous stimulation (Fig. [9e](#Fig9){ref-type="fig"}). However, the bled/injured mussels after a *Vibrio* challenge did not lose the myticin C labeling (Fig. [9f](#Fig9){ref-type="fig"}) probably because of the effect of the previous stimulus that prevented the hemocyte population from decaying after a subsequent bacterial infection.Figure 9FACS results of the immunocytochemistry of myticin C. (**a**) Criteria to choose myticin C positive cells (M3). (**b**) Myticin C positive cells in a representative sample. (**c**) Statistical analysis of the percentage of granulocytes positive for myticin C immunocytochemistry. Standard deviation is indicated for the four replicates. Significant differences (p-value \< 0.05) among groups are indicated by different letters. (**d**--**f**) Representative overlay of histograms for FL1 intensity of granulocytes (R1) population.

Discussion {#Sec14}
==========

Next generation sequencing technologies are particularly valuable in the study of non-model organisms because they supply the scarce availability of reference genomes, and lead the way to understand many biological processes that could not be studied due to the lack of cellular lines or antibodies. Most of the transcriptomic analyses in the field of bivalves have been conducted using pools of animals or tissues^[@CR43],[@CR44]^. When biological materials, such as hemocytes, are scarce, this is a good solution^[@CR45]--[@CR48]^. Biological pools help us to eliminate individual differences and concentrate on the clear patterns of a particular experiment; determining, for example, how bivalves react to a bacterial infection^[@CR45]--[@CR48]^. Although this approach is totally correct, we cannot forget that we are dealing with wild animals. They are not laboratory homogeneous strains, and they differ in their genetic and physiological backgrounds. In this work, we were interested in the specific transcriptomic response of individual mussels before and after a bacterial infection.

The first unexpected result occurred when we realized that the percentage of DEGs shared among 3 different mussels after the same stimulus was very small. Each individual had its own repertoire of genes, and this result occurred both in controls (reaction of 3 individuals against FSW) and in infected mussels (reaction of 3 individuals against *V*. *splendidus*). The high diversity found in mussels for many genes such as antimicrobial peptides (AMPs), fibrinogen-related proteins (FREPs) or C1q domain containing proteins^[@CR21]--[@CR23]^ increased the differences among individual mussels, which express a specific repertoire of transcripts. The interindividual variability in the basal expression of AMPs was previously reported in oysters and mussels^[@CR15],[@CR21],[@CR49]^. The variability of bivalve responses would be related to the antigenic environment of each individual and would also show the genetic diversity of these animals. We wonder if we are losing important information when we analyze the immune response in pools of animals instead of using individuals and whether we are underestimating the possibility of the individual immunity as part of a possible social immunity of mussels^[@CR50]^, as their impressive resilience and survival capacities point out.

The second striking result was that the 3 control mussels showed 1,562, 486 and 751 DEGs (using this statistical threshold: FC \> \|2\| and Bonferroni corrected p-value \< 0.05), which were probably modulated by a simple tissue injury (FSW injection). The hemocytes are the bivalve immunocompetent cells, but they are also involved in other physiological processes such as basic homeostasis and wound and shell repair^[@CR51],[@CR52]^. As mussels host diverse and abundant microbial communities, they probably do not respond to all the microorganisms that they are in contact with. Only if there is tissue damage can they trigger an immune reaction. Cell proliferation and migration related genes were overexpressed in control mussels: stathmin, key protein in the cell cycle regulation^[@CR53]^; low affinity epsilon Fc receptor, with essential roles in defense cell growth and differentiation^[@CR54]^; cell division cycle-associated protein 2, which controls the cell cycle^[@CR55]^; signal transducer and transcription activator, member of the STAT family and regulator of cellular immunity, proliferation and differentiation^[@CR56]^; and structural maintenance of chromosomes protein 4, a protein required to enter into mitotic phase^[@CR57]^. In contrast, genes related to the regulation of cell death were down modulated: cell death specification protein 2, with a direct link to the activation of apoptosis^[@CR58]^; protocadherin Fat 1, a tumor suppressor and inhibitor of cell migration^[@CR59]^; and GTPase IMAP family member 4, member of the GIMAP family which may play a role in defense cells differentiation and apoptosis^[@CR60],[@CR61]^. These results were in agreement with the flow cytometry results. When an initial injury or bleeding occurred, changes in hemocyte structure were milder compared with those observed after a bacterial infection. The response of bivalves to tissue injury or different danger signals has not been explored and deserves to be further investigated. These results also highlight the importance of using appropriate controls in our studies. In fact, the acclimation and rest period is important to rule out possible responses due to stress or handling^[@CR62]^. The usual procedure of injection in the adductor muscle could induce an injury that mussels feel as a danger signal. Therefore, if our aim is to analyze the bivalve immune response against a pathogen, we must clearly identify this specific response without the reaction against the tissue damage produced by the injection. Although our transcriptome study did not take into account a non-injected control, in an additional experiment we could confirm that the initial hemolymph extraction did not induce significant changes in gene expression after one week except for C1q, a pleiotropic gene involved in innate immunity and also in shell repair^[@CR63]^, a process triggered by the notch in the shell for hemolymph sampling.

The third finding emerged when we detected AMPs included in the top 25 most expressed genes from controls (exclusive to controls and not found in infected mussels) (Table [2](#Tab2){ref-type="table"}). To date, the main function recognized for AMPs has been the direct killing of microorganisms^[@CR64]^. We previously reported the biological characteristics of mussel myticins: activity against bacteria^[@CR41]^, activity against molluscan, fish and human viruses^[@CR25],[@CR38]^, chemotactic activity^[@CR25]^, etc. Myticins are very diverse, but they are not usually regulated at the transcript level after infection, they are stored in hemocyte granules ready to act when needed^[@CR42]^. However, even with these premises, it was unexpected that myticins were up regulated in controls. Their expression was increased after a tissue injury, a danger signal. A non-injected control after the rest period in the RNA-Seq could further confirm this results, and although it was not included in the initial experimental design, the qPCR analysis supports our observations that the response to a wound and to a pathogen is quite diverse. The same response is observed with myticusin in that although it has only been described as an antimicrobial peptide^[@CR65]^, it probably can also be involved in cell proliferation and chemotaxis such as myticins are. Therefore, why a tissue injury and not a bacterial infection triggers their transcription? Our hypothesis is that a bacterial signal (PAMP) triggers the release of granules to let the AMPs fight the bacteria, and the tissue injury (DAMP) triggers the mechanisms to produce and store these AMPs. The explanation could be a new and unrecognized function of these AMPs due to exposure to many microorganisms, where wounds matter more than a potential pathogen.

In wild conditions a tissue damage is a possible cause of infection and it makes sense that these animals could be preparing themselves for a possible infection after being wounded. Interestingly, flow cytometry analysis suggests that the myticins could also help in the control of the infection process since the reaction triggered by the bleeding (tissue injury) was enough to avoid the granulocyte and the myticin C+ cell population decrease caused by bacteria, in concordance with the transcriptomic results, which show genes related to cell proliferation. Although further studies will be needed to confirm cell proliferation and hematopoiesis after tissue injury in *M*. *galloprovincialis*, these results open the door for new research topics regarding innate immunity in bivalves.

Conclusions {#Sec15}
===========

In summary, we have raised the question of using an individual approach when facing mussel transcriptomics, which may be true for other organisms as well. We have demonstrated that every single individual expression profile can be very different in terms of the number of genes expressed and in the magnitude of their expression. A careful experimental design should be carried out especially with non-model species, as their individual transcriptomes can be quite diverse. Additionally, some common experimental procedures in bivalves such as shell notching or injection of treatments could trigger unexpected responses such as an unspecific inflammatory response or defense cell proliferation. We found that myticin C may play new roles in preparing the mussels for future pathogenic processes after a danger signal. A tissue injury is a breach of the first defense barrier, which could be easily followed by a subsequent infection.

More studies should be conducted in the future to understand more about these processes. However, currently, we undoubtedly can say that mussels are anything but simple and that more revelations will appear in the study of "non-model mussel immunity".

Supplementary information
=========================

 {#Sec16}

Supplementary Information

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Magalí Rey-Campos and Rebeca Moreira contributed equally.

Supplementary information
=========================

**Supplementary information** accompanies this paper at 10.1038/s41598-019-39870-3.

This work was conducted with the support of the projects AGL2015-65705-R (Ministerio de Economía y Competitividad, Spain), IN607B 2016/12 (Consellería de Economía, Emprego e Industria - GAIN, Xunta de Galicia) and VIVALDI (678589) (EU H2020). MRC thanks additional funding from the Spanish Ministerio de Economía y Competitividad, for her predoctoral contract (BES-2016-076302).

B.N. and A.F. conceived and designed the project. R.M. and M.R.C. did the mussel infections, sampling and RNA extraction. A.F., R.M., V.V.M. and C.G.E. performed the reads trimming, assemply, RNA-Seq and statistical analyses. B.N., R.M. and M.R.C. analyzed the generated data. R.M. did the immunocytochemistry. B.N. wrote the manuscript. All listed authors revised, edited, read and approved the manuscript.

The raw reads of each sampling point and treatment were deposited in the NCBI database with the accession numbers SAMN09096971 to SAMN09096982. Accessible through <https://www.ncbi.nlm.nih.gov/biosample/9096971> to.../909682.

Competing Interests {#FPar1}
===================

The authors declare no competing interests.
