Integration of Transcriptomics and Proteomics Improves the Characterization of the Role of Mussel Gills in a Bacterial Waterborne Infection

In recent years, the immune response of mussels (Mytilus galloprovincialis) has been studied at the transcriptomic level against several bacterial infections. As a result, different immune mechanisms have been revealed, including both conserved essential innate pathways and particularities of the mussel immune response according to its nature and environment. However, there is often a lack of functional verification because mussels are a non-model species and because transcriptomic and proteomic information is not always well correlated. In the current study, a high-throughput quantitative proteomics study coupled to LC-MS/MS analysis using isobaric tandem mass tags (TMTs) for protein labeling was employed to study the mussel gill immune response to a Vibrio splendidus bath (waterborne) infection at a functional protein level. A total of 4,242 proteins were identified and quantified, of which 226 were differentially expressed (DEPs) after infection, giving to the study a depth that was lacking in previous proteomic studies of the bivalve immune response. Modulated proteins evidenced an important cytoskeletal disruption caused by bacterial infection. A conserved network of associated proteins was modulated, regulating oxidative stress and NF-kB inflammatory responses and leading to innate immunity effectors. Proteomic results were submitted to an integrated analysis with those obtained in a previous transcriptomic approach with the same infection. Half of all the quantified proteins had a concordant transcriptomic expression trend, but this concordance increased when focusing on the DEPs. The correlation was higher within the immune-related DEPs, and the activation of the conserved NF-kB pro-inflammatory pathway was the main response in both approaches. The results of both techniques could be integrated to obtain a more complete vision of the response.


INTRODUCTION
Mussels (Mytilus galloprovincialis) have great importance in scientific research for different reasons. These bivalves are among the most produced species in aquaculture, having a nonnegligible impact on the European economy (Avdelas et al., 2020). Mussels are filter feeders in contact with large amounts and diverse potential pathogens, and despite this, they do not commonly suffer severe mortality episodes. Mussels can use their efficient innate immune system to deal with different pathogenic agents in a specific way (Costa et al., 2009;Balbi et al., 2013). This attribute is supported by the vast genetic variability that characterizes this species, with several expanded immune gene families as immune recognition C1q factors (Gerdol et al., 2011) or pro-inflammatory cytokines as IL-17 (Saco et al., 2021). The recent discovery of mussel genomic organization as a pangenome revealed presence/absence phenomena in the genetic repertoires of these animals, where a core set of genes was accompanied by great amounts of dispensable and variable genes . Gene families with immune importance are enriched in the dispensable and variable genomic fraction, such as the antimicrobial peptide myticins (Rey- .
To reveal how mussels use their immune system to deal with infections, transcriptomic approaches have been quite successful (Figueras et al., 2019). Classical approaches have studied the transcriptome of the main immune cells, hemocytes, against in vitro and in vivo stimulation with different pathogenic bacteria or pathogen-associated molecular patterns (PAMPs) (Rey- Campos et al., 2019a,b;Moreira et al., 2020). In addition to hemocytes, mucosal immunity should have great importance in filtering bivalves such as mussels. Transcriptomic studies indeed revealed interesting immune expression patterns in more tissues beyond hemocytes (Moreira et al., 2015). A transcriptomic study performed in mussel gills against a bacterial bath (waterborne) infection tried to mimic how mussels recognize an infection from the surrounding water as in natural conditions. The modulation of different pathogen recognition and immune triggering genes revealed that gills may play a natural first step in the activation of the immune response in these animals (Saco et al., 2020). These transcriptomic studies described both essential key immune mechanisms that would therefore be conserved during evolution up to higher vertebrates, as well as particularities adapted to the mussel environment, such as the high variability of certain genes (Rey- , to survive in the ever-changing marine environment. In addition to transcriptomic studies, more functional approaches have been used to improve our understanding of the mussel immune system. Traditional proteomic analyses consist of two-dimensional electrophoresis (2-DE)-based methods and have proven to be useful in several analyses on mussels and other non-model species. However, high-throughput and more effective techniques emerged, allowing the identification of many more proteins in each analysis and a more accurate relative quantification, such as iTRAQ (isobaric tag for relative and absolute quantitation) and TMT (tandem mass tag) approaches, both coupled to liquid chromatographytandem mass spectrometry (LC-MS/MS). They allow much larger identification rates, higher precision and accuracy in relative quantification than traditional 2-DE gel assays (Aslam et al., 2017).
Different proteomic studies have been conducted with mussels and other bivalves subjected to different stimuli. Mussel and other bivalve gills have been used in proteomic studies on environmental stress or in studying the bioaccumulation and toxic effects of different water pollutants. Traditional 2-DE analyses identified a limited number of differentially expressed proteins (DEPs) after the analysis of all the protein spots (Manduzio et al., 2005;Ji et al., 2013b;Song et al., 2016;Chen H. et al., 2018). A big leap in the magnitude of the analyses and the number of identified and quantified proteins has been reported with the most recent LC-MS/MS technological advances Sánchez-Marín et al., 2021).
Proteomics has also been employed for the study of the immune response in bivalves, as in the scallop hepatopancreas infected by Vibrio (Huan et al., 2011), for characterization studies of the immune proteins expressed in hemocytes from different bivalves (Campos et al., 2015;Leprêtre et al., 2020), for the study of the immune response in mussel gills to the intramuscular injection of bacteria (Ji et al., 2013a), among other studies (Castellanos-Martínez et al., 2014;Dinguirard et al., 2018;Jiang et al., 2018;Smits et al., 2020). Most likely, the most studied infection model in bivalve proteomics is oyster infection by herpesvirus and virus-like stimuli (Corporeau et al., 2014;Masood et al., 2016;Leprêtre et al., 2021).
Few immune proteomic responses to bacterial infections have been reported in mussels or in other bivalves, and the evolution of proteomic techniques since the 2DE analyses caused some past studies to have severe limitations. The objective of the current work was to carry out a high-throughput proteomic study in gills of the marine mussel Mytilus galloprovincialis subjected to a waterborne infection with Vibrio splendidus, simulating natural infection conditions by a bloom of certain infective bacteria. The infection model followed the same design as a previous transcriptomic experiment, which had revealed key mechanisms in gills in terms of pathogen recognition and activation of the immune response (Saco et al., 2020). Multiplexing of TMTlabeled protein samples coupled to liquid chromatographytandem mass spectrometry (LC-MS/MS) was applied to this proteomic study, and the results were also complemented with those obtained by the previous transcriptomic approach with the objective of gaining a more complete understanding of the marine mussel response to this infection model.

Animals, Experimental Infection, and Sampling
Adult mussels (Mytilus galloprovincialis) were obtained in May from a mussel farm (Vigo, Galicia, NW Spain) and maintained in open-circuit filtered seawater (FSW) tanks at 15 • C with aeration. The mussels were fed daily with Phaeodactylum tricornutum and Isocrhysis galbana. The mussels were used for the experiments after at least 1 week of acclimatization.
Vibrio splendidus (reference strain LGP32) was cultured in tryptic soy agar (TSA) plates at 22 • C for 24 h before use. Bacterial suspensions were prepared in FSW before infection. Serial dilutions were made to calculate the number of colony forming units (CFU)/mL 24 h after inoculation.
Waterborne infections were performed by treating mussels with a concentration of 10 8 CFU/mL V. splendidus in tank water. Control mussels were maintained in another tank with only filtered seawater (FSW). Gill samples from 5 control and 5 infected mussels were taken 24 h after the infection. To perform the sampling, gills were cautiously separated from gonads and mantle and frozen directly at −80 • C.

Protein Quantification
Gill samples were homogenized with several ultrasonication cycles (5-s pulses, with 30-s ice stops and 20% amplitude) in lysis buffer consisting of 7 M urea, 2 M thiourea and 4% CHAPS. Once the tissue was correctly disaggregated and proteins released and solubilized in the buffer, samples were centrifuged (14,000 g, 60 min, 4 • C), and the supernatants were collected in 50-µL aliquots and frozen at −80 • C.
Protein quantification was performed following the Bradford method in microplates. Both BSA standards and samples were analyzed in triplicate, and the absorbance was read at 595 nm in a BIO-RAD iMark TM Microplate Reader. After quantification, samples were diluted to an equal concentration of 0.5 µ g/µ L.

Protein Precipitation
To perform the precipitation step, 600 µL of prechilled acetone (−20 • C) was added to 100 µL of each sample. Samples were incubated for 2 h at −20 • C to ensure precipitation. After the incubation step, the samples were centrifuged (8,000 g, 10 min, 4 • C). The supernatant consisting of acetone was discarded, and the pellet was dried for a maximum of 10 min before adding 100 µL of 100 mM TEAB (triethylammonium bicarbonate). A brief bath sonication was performed to resuspend the pellet. Then, 5 µL of TCEP (Bond-Breaker TCEP Solution) was added, and samples were incubated under soft shaking conditions (300 rpm) for 1 h at 55 • C. To alkylate the samples, 5 µL of 375 mM iodoacetamide (diluted in 100 mM TEAB) was added, and the samples were incubated for 30 min at room temperature and under dark conditions. Finally, 10 µL of trypsin (0.1 µg µL −1 in 100 mM TEAB) was added, and the samples were incubated at 37 • C overnight with soft shaking at 300 rpm. Because we began protein precipitation with 100 µL of homogenized protein sample at 0.5 µg/µL, we expected to precipitate a maximum of 50 µg of proteins per sample.

Peptide Fractionation
The 10 already labeled samples were combined in the same tube at equal concentrations, for which 10.3 µL was added from each sample, and then the final volume was raised to 300 µL with TFA (trifluoroacetic acid) 0, 1%. The mix consisting of 300 µL with an estimated peptide amount of 30 µg was submitted to peptide fractionation using the Pierce High pH Reversed-Phase Peptide Fractionation Kit (Thermo Fisher Scientific) following the kit protocol. After conditioning and washing the column, a 300-µL sample mix was applied to the column. Elution solutions consisted of a series of 8 solutions with increasing acetonitrile concentrations from 12.5 to 50% supplemented with 0.1% triethylamine. Therefore, 8 fractions or elutions were obtained, and each one of them was divided in 2 aliquots of 150 µL with approximately 2 µ g of peptides.

Liquid Chromatography-Tandem Mass Spectrometry Analysis
Each 150 µL aliquot was dried using a SpeedVack (Thermo Fisher Scientific), and 2 µg of peptides from each aliquot was suspended in 13 µL of 0.5% formic acid. A brief bath sonication for 10 min was followed by centrifugation (14,000 rpm, 1 min, 10 • C). From each aliquot, 10 µL corresponding to 1.5 µg of peptides was analyzed with LC-MS/MS using a Proxeon EASY-nLC II liquid chromatography system (Thermo Fisher Scientific) coupled to an LTQ-Orbitrap Elite mass spectrometer (Thermo Fisher Scientific).

Protein Identification
To identify the proteins using the MS/MS spectra, a customized protein database was generated from the six-frame translation of the 162,167 contig sequences from the mussel gill transcriptome performed under the same bacterial bath infection stimulus (Saco et al., 2020). MS/MS spectra were searched against this database using PEAKS Studio v8.0 software (Bioinformatics Solutions Inc., Waterloo, Canada). A contaminant database (CRAPome) was included, as well as a decoy sequence database to calculate the false discovery rate (FDR) in the peptide identification, which was set to 1% FDR for PSMs (Peptide-Spectrum Matches). Trypsin was indicated to be the enzyme for protein cleavage, and only one missed cleavage was allowed. For peptide identification, the precursor and fragment ion mass tolerances were set to 10 ppm and 0.02 Da, respectively. Carbamidomethylation of cysteine and TMT 10 plex were set as fixed modifications, and oxidation on methionine was set as a variable modification. For protein identification, the number of matched peptide sequences was set to ≥ 2, the number of unique peptides to ≥ 1 and the protein identification PEAKS score to ≥ 20. In addition to the mentioned self-built protein database, for comparative purposes, the same procedure (both the identification step and the downstream analyses) was repeated using the NCBI non-redundant protein database restrained to Mollusca (23/05/2019; 749,599 sequences).

Protein Quantification
For protein quantification, only unique peptides were considered. The spectrum filter for the reporter ions was set to FDR 1% and quality ≥ 10. Intra-experiment normalization was performed by calculating an equalizer factor to normalize the reporter ion intensities for each TMT 10plex label using the total ion counts per channel. After intra-normalization, a logarithmic transformation was applied to the normalized protein abundance data to meet the normality and homoscedasticity assumptions necessary to carry out parametric statistics, which is a standard procedure for proteomic datasets. Differential expression analysis was performed with a t-test followed by a multi-test correction. These corrections were calculated using SGoF + software v.3.8 (Carvajal-Rodriguez and de Uña-Alvarez, 2011) following the strategy previously described in Diz et al. (2011). Adjusted p-values for false discovery rate (q-values, Storey, 2003) were obtained and set to ≤ 0.05 as the criterion to define DEPs.

Proteomic Data Analysis
Annotation of the identified and quantified proteins was derived from the annotation of the transcriptome used to construct the identification database. This annotation was performed using a Blastx approach on OmicsBox (BioBam Bioinformatics, Valencia, Spain) (Götz et al., 2008) against UniProt/Swiss-Prot with an e-value of 10 −3 . The annotation of DEPs was manually verified with Blastp against the NCBI non-redundant protein database.
Two enrichment analyses, both of which were performed with generic GO slim terms and with specific GO terms, were performed with Fisher's exact test and FDR ≤ 0.05. The first analysis considered all the identified and quantified proteins as the test set against the reference set (identification database) to analyze which types of proteins were identified in the proteome (a sort of functional clustering). The other enrichment analysis was performed with the modulated proteins against the reference (all the identified proteins) to associate enriched terms with the response to the infection. Enrichment analyses were represented using a ratio between the weight of each term in the test set and the weight in the reference set. Thus, processes with better test/reference ratio will be of interest as they are more enriched in the studied test set than in the reference set.
A protein interaction network was built using the STRING algorithm (Szklarczyk et al., 2015) and the modulated proteins. Default parameters were employed, including the use of human homologs to the mussel proteins as a reference, a common procedure with non-model species that are not included in the server. Only high confidence interactions were shown (0.7 minimum required interaction score) and disconnected nodes were removed. Enriched biological processes in the network were also analyzed with the STRING algorithm. A second STRING network was constructed, following the same criteria (0.7 minimum score and no disconnected nodes), using homolog proteins from a phylogenetically closer species, the marine invertebrate Strongylocentrotus purpuratus. Both analyses were then compared to check which interactions are recovered using the homologs from each species and which ones are shared by both analyses.
Heatmaps and PCA were constructed using ClustVis (Metsalu and Vilo, 2015). For the PCA, singular value decomposition was used, and for the heatmaps, distances were calculated with the unit variance scaling method. Clustering was performed based on average correlation. Schemes were created with BioRender.com.

Proteomic and Transcriptomic Data Analysis
Transcriptome analysis performed on gills from the same marine mussel species (M. galloprovincialis) after waterborne V. splendidus infection (Saco et al., 2020) was employed both as a database for proteomic analysis and to retrieve transcriptomic expression of the detected proteins. The experimental design was the same for both approaches, with the same concentration of bacteria in the infection (10 8 CFU/mL V. splendidus) and the same sampling time (24 h). The transcriptomic analysis was performed with 3 biological replicates for each control and infected group instead of the 5 biological replicates employed in the proteome analysis. Reads accessible from Bioproject PRJNA638821 under the accessions SRR11996464, SRR11996686, and SRR11996723 for control samples and SRR11996734, SRR11996735, and SRR11996743 for infected samples were assembled with the parameters indicated in Saco et al. (2020). Expression of transcriptomic contigs equivalent to the modulated proteins in the proteome was retrieved in transcripts per million (TPMs). Log2 transformation of the fold change values was used for the correlation analyses between proteomic quantification and transcriptomic expression. Linear regression algorithm analyses were made in R.

Characterization of the Mussel Gill Proteome
The use of the specific mussel gill transcriptome (Saco et al., 2020) allowed us to achieve better protein identification rates in the analysis than the NCBI database ( Table 1). After the LC-MS/MS analysis, 245,047 MS/MS spectra were obtained, and using the mentioned mussel gill transcriptome as a database, 35,217 peptide spectrum matches (PSMs) were retrieved. A total of 4,388 proteins were identified from 28,141 peptide sequences based on the parameters used. Considering only unique peptides that contained signals from reporter (TMT) ions, 4,242 proteins could be quantified. Protein relative abundance data could be retrieved for those proteins.
Enrichment analyses were performed with all the identified and quantified proteins to reveal a general characterization of what type of proteins dominate the mussel gill proteome. For this purpose, GO Slim terms were employed. This broad overview of the proteome is shown in Figure 1A for the molecular functions, Figure 1B for the cellular components and Figure 1C for the biological processes. More than half of the proteins were related to binding functions, both inside the cell and in the plasma Information about the transcriptome database is highlighted in gray as this database allowed better results and was the one employed for the whole downstream analysis.
FIGURE 1 | General characterization of gill proteome from the marine mussel (M. galloprovincialis) based on the confidently identified proteins. Identified proteins in the proteomic analysis were submitted to an enrichment analysis with GO-Slim terms to obtain a generic view of the type of proteins that were identified in the gill proteome. The top 15 groups of GO-Slim terms with the highest number of representative proteins are shown for (A) "molecular function"; (B) "cellular component"; and (C) "biological processes." (D) The most present specific biological processes are also represented without GO-Slim generalization.
membrane. These proteins are implicated in functions such as response to stress, signal transduction and immune processes, among others. Biological processes of the general proteome were also studied with specific GO terms, without GO slim generalization, specifying the importance of gill proteins in the recognition of exogenous antigens and in the transduction of C-type lectin, Fc-epsilon receptor, interleukin or Wnt signaling pathways, among others ( Figure 1D).

Modulated Proteomic Response in Mussel Gills Against Bacterial Bath Infection
Based on the quantification data, after applying the statistical methodology explained, a final number of 226 DEPs were obtained (q-value ≤ 0.05). The DEPs and the rest of the quantified proteins are displayed in Supplementary Table 1. Modulated proteins displayed a clear difference between the two conditions (Figure 2A), and no great variability was observed between biological replicates within the control and infected groups of samples ( Figure 2B). The biological processes of the modulated proteins were first revealed with enrichment analyses of the upregulated and downregulated DEPs. Functions such as activation of innate immune cells/response and response to oxidative stress are among the most common upregulated proteins against infection ( Figure 3A). Other apparently less-related immune terms, more focused on regulatory functions, response to ions, and regulation of gene expression, are also directly related to the response to the infection. Cilia movement, as expected in the gills, involved in cell motility was the dominant process in the downregulated proteins ( Figure 3B). Immune processes downregulated with the infection were response to cytokine stimulus, response to external stimulus and oxidation-reduction processes.  Focusing on specific functions of the modulated proteins, many of them were related to the cytoskeleton or to immune pathway regulation, as revealed by the GO Slim biological processes represented in Figure 4A. Table 2 also displays the identity of the DEPs specifically related to the mentioned cytoskeleton-related and immune-related functions Frontiers in Marine Science | www.frontiersin.org that dominated the response. All the cytoskeletal proteins that were modulated (main components such as actin and tubulin, transport proteins such as dynein and many regulators) evidenced general cytoskeletal disruption and reorganization with the infection, which will be discussed in the following sections. DEPs with immune implications were involved mainly in the regulation of the NF-kB pro-inflammatory response and the response to oxidative stress. An interaction net was constructed with human homolog proteins to the mussel DEPs ( Figure 4B). The net was significant (PPI enrichment p-value < 1.0e-16), meaning that proteins modulated in the analysis showed a significantly higher degree of interactions within each other than random proteins. The core of the interaction network was enriched in cytoskeletal and immune DEPs (Figure 4B and Supplementary Table 2). The comodulation of this net of associated proteins in the gill mussel proteome could suggest a conservation of these protein interactions in the mussel innate immune response with respect to the higher vertebrates in which protein interactions are usually described. Using homolog proteins from a phylogenetically closer species, the sea urchin (Strongylocentrotus purpuratus), another significant interaction net was constructed (Supplementary Figure 1A). Interactions detected using sea urchin were mostly shared with the human homologs analysis, since sea urchin interactions information mainly relied on data from homologs in the more studied model species (Supplementary Figure 1B and Supplementary Table 3).
An effort was made to summarize how the modulated proteins regulate the NF-kB immune response according to their function (Figure 5). The NF-kB inflammatory response is initiated after the recognition of pathogen molecular patterns with pro-inflammatory cytokine signaling. Then, an intracellular pathway liberates NF-kB from its inhibitory factor and allows it to express pro-inflammatory and innate immune effectors. The trend was clear: proteins that can positively regulate the pathway were mainly upregulated, while inhibitory proteins were all downregulated. Figure 5 also displays the response related to oxidative stress, comprising ROS-and NOS-related proteins, as well as the calcium-regulatory response that affects them. This redox stress response caused by the infection would also affect the NF-kB regulatory net, as both responses are associated, mainly through the Pirin protein.
Other immune-related proteins were also modulated, affecting different functions than previously seen NF-kB and ROS. These proteins are also listed in Table 2 and display more effector roles against infection: Tubby protein homolog stimulates phagocytosis, alanine and aspartate aminotransferases, which are related to apoptosis under bacterial infection, Kinesin-like proteins related to antibacterial lytic granule secretion, Histone deacetylase 6, and granulin, among others.
Gene expression was retrieved from the equivalent transcriptomic analysis (Saco et al., 2020) for all the proteins detected in the proteome. The expression of the equivalent contigs is represented for immune-related DEPs in Figure 6A. There was a general correlation between these modulated immune proteins and their transcriptomic expression trends. Furthermore, the transcriptomic modulated response agreed with the proteomic study in key processes, as will be discussed later. Transcriptomic and proteomic comparison can shed some light on the putative reasons why some expected terms do not appear in the proteomic study, i.e., the great transcriptomic response related to pathogen recognition and cytokine signaling (Saco et al., 2020) was not modulated in the proteomic experiment. Figure 6B shows how contigs corresponding to the detected proteins have higher average expression than the overall transcriptome. The depth of both approaches is a factor to consider as well since the proteome retrieved 4,388 proteins and the transcriptome retrieved 162,167 contigs.
Correlation analyses between the quantitative proteomic and RNA-Seq transcriptomic experiments revealed that 139 out of 226 DEPs (60%) were concordant gene-protein interactions (gene expression and protein quantification fold change values with the same trend). Fold change values for these concordant DEPs had a significant correlation of R 2 = 0.4192 ( Figure 6C). The correlation was more noticeable within the modulated immune proteins (Figure 6C), which is especially important since DEPs related to other functions, such as the cytoskeleton, showed null correlation (R 2 = 0.0641). When focusing on the entire proteome characterization instead of the modulated response, the agreement decreased, with half the terms concordant (2,137 out of 4,242 quantified proteins) and with less correlation in the magnitudes of their fold changes ( Figure 6D). The complementarity between proteomic and transcriptomic experiments for the immune implicated DEPs supposes a crossvalidation for the implication of certain biological responses to the infection, as NF-kB and the oxidative responses. Other groups of proteins/transcripts did not show complementarity, as the cytoskeletal DEPs in the proteome or the modulated immune recognition response in the transcriptome. Immune recognition proteins of interest were identified and quantified in the proteome but they were not modulated in the study (listed in Table 2). Relative quantification of some of these important but non-modulated proteins is shown in Supplementary Figure 2, along with their corresponding transcriptomic expression: AIF-1, Myd88 and a toll-interacting protein.

DISCUSSION
After the LC-MS/MS analysis, a mussel gill proteome with high protein identification and quantification rates was obtained. Identification rates in proteomic analyses, especially in nonmodel species, are improved by using customized databases from appropriate and specific transcriptomes, as we did, rather than using general databases such as NCBI (Romero et al., 2019). The number of identified proteins is similar to the number of proteins reported in other studies in bivalves following high-throughput proteomic approaches based on ITRAQ or TMT sample labeling Sánchez-Marín et al., 2021)., far above the numbers obtained in traditional 2-DE assays (Manduzio et al., 2005;Ji et al., 2013a;Song et al., 2016;Chen H. et al., 2018). General characterization of the mussel gill proteome revealed functional information on terms that are typically found in gill transcriptomes, in particular the biological processes related to Relevant DEPs related to cystoskeletal functions are represented as this biological process dominated the modulated response (Enrichment Analysis with GO Slim terms performed with the DEPs as test set). Immunerelated DEPs are displayed as well because of their importance in order to understand the proteomic immune response to the studied infection. It was also included a list of non-modulated but significantly identified and quantified proteins of interest for the immune response.
recognition and signaling pathways (Moreira et al., 2015). The proteins identified corroborated and expanded the results of a previous 2-DE proteomic characterization of the M. edulis mussel gill (Rocher et al., 2015). In the current work, the study for the first time of the proteomic modulated response to a bacterial waterborne infection with Vibrio splendidus in mussel gills demonstrated that modulation against the infection occurs not only at the transcriptomic level but also at the protein scale. The study of the same infection model by the combination of proteomic and transcriptomic approaches reported different agents taking part in the same response: the main modulated genes in the transcriptome were focused on pathogen recognition terms and cytokines activating the NF-kB inflammatory response (Saco et al., 2020), while the proteome was dominated by proteins that positively regulate the activation of the NF-kB pathway at an intracellular level.
Agreement between transcriptomic and proteomic analyses is not always obtained. High-throughput proteomic analysis of oyster response to viral-like simulation with dsRNA revealed disparities between transcriptomic and proteomic analyses (Masood et al., 2016). Discrepancies had been observed before in gills of mussels with bacteria injected in the muscle but in a lower number of proteins/genes in the comparison because of the 2D analyses (Ji et al., 2013b). Other studies have reported a generally good correlation between transcriptomic and proteomic results, such as in the genes/proteins expressed in the sperm of different mussel species, even though not for all cases (Romero et al., 2019). In our study, because of using the indicated transcriptome as a database for protein identification, it was possible to recover the transcriptomic expression of the contigs that corresponded exactly to the proteins identified and quantified. There was a greater concordance in terms of expression trends between both techniques when focusing on DEPs (60% concordant terms with R 2 = 0.4) than with the complete proteome (50% concordant terms with R 2 = 0.2). When classifying these DEPs according to their function, the correlation was observed to be much more important for certain functions, such as immune proteins (70% concordant terms with R 2 = 0.5), while it was null for the cytoskeleton-related DEPs (60% concordant terms with R 2 = 0.0641). The fact that the transcriptomic-proteomic complementarity within a certain experiment is better for some specific functions is consistent with what has been observed in other studies (Xiong et al., 2015;Romero et al., 2019) and points to a functional cross-validation of certain biological responses to the stimulus.
The modulation of cytoskeletal proteins was a clear result in the proteomic analysis and was not among the modulated responses in the transcriptomic study of the same infection model. The presence of these proteins in the analysis, according to what the lack of agreement with the transcriptomic results points out, would not be due to a modulation in the gene transcription and subsequent protein production. These proteins would be modulated because of the physical effects of infective bacteria disrupting the cells and defensive cells (hemocytes) migrating to the infection site (Wu et al., 2013;de Souza Santos and Orth, 2015). Cytoskeleton disruption FIGURE 5 | Scheme showing the differentially expressed proteins (DEPs) that regulate the intracellular inflammatory NF-kB pathway and the oxidative and calcium regulation response. The names of the upregulated DEPs are indicated in green, while downregulated DEPs are indicated in red. The activity that each protein exerts is represented with arrows if it is an activation or with an inhibition sign if it is an inhibitory protein. The pathway represented in the scheme follows the same order that is explained in the Discussion. Proteins with functions enhancing the NF-kB pathway were generally upregulated (green), while inhibitory proteins were downregulated in the experiment (red). The scheme shows how the NF-kB response is activated and the oxidative stress is maintained under control.
is a direct consequence of infective bacteria entering host cells (Barbieri et al., 2002;Yoshida and Sasakawa, 2003). The Vibrio splendidus strain employed has been described to use OmpU porins to attach and invade mollusc cells, which would subvert the host cell actin cytoskeleton, causing this general cytoskeleton disruption (Duperthuy et al., 2011). Revising mussel and other bivalve proteomes studying Vibrio and other different infection stimuli, the detection of cytoskeletal proteins among the modulated proteins is recurrent (Huan et al., 2011;Wu et al., 2013;Corporeau et al., 2014;Zhang et al., 2014;Wang et al., 2017). Within the cytoskeletal-related DEPs in our work, associated proteins that could interact with each other were found, as the homolog proteins interaction network revealed. Tubulin is modulated and is known to be the major binding site of bacterial peptidoglycans, which can regulate its polymerization (Dziarski et al., 2000). Tubulin-related proteins such as microtubule organizing center proteins and tubulin polymerization-promoting proteins were also modulated. Modulated dynein proteins have also been described to interact with bacteria, particularly with the Pilin protein of gram-negative bacteria, such as Vibrio splendidus (Kausar et al., 2013).
However, most of the cytoskeletal DEPs were related to actin modulation, i.e., thymosin (a regulator of actin polymerization), nephrocystin-4 (actin network regulation in multi-ciliated epithelial cells, such as mussel gills), Cdc42-interacting protein 4 (which reorganizes the actin cytoskeleton during endocytosis), F actin capping protein and even cytoplasmic actin protein.
Profilin is another DEP regulator of actin polymerization that has been directly related to the actin-based intracellular motility of different bacteria (Geese et al., 2000;Mimuro et al., 2000;Liverman et al., 2007), so it has been typically modulated with bacterial challenges in a skin mucus cod proteome (Rajan et al., 2013) and in bivalve transcriptomes (Somboonwiwat et al., 2006;Araya et al., 2010).
Cytoskeletal disruption caused by infection affects adhesion, phagocytosis and oxidative burst responses in host cells (Tanguy et al., 2013). Modulated proteins related to all these functions were found. In a proteomic study of oysters infected with OsHV-1 herpesvirus, general cytoskeletal disruption (actin, dynein, F actin capping protein) was associated with oxidative stress (modulation of catalase, glutathione-S-transferase, superoxide dismutase). Actin could be related to that oxidative response, FIGURE 6 | Comparative analysis between proteomic quantification and transcriptomic expression for the differentially expressed proteins (DEPs). (A) For the immune-related DEPs, the mRNA expression of their equivalent contigs was retrieved from a previous transcriptomic study in gills against the same bath infection, which was also used as a database for proteomic identification (Saco et al., 2020). The heatmap shows the expression in TPMs (transcripts per million) for these immune DEPs in each individual sample used for the transcriptomic analysis (three biological replicates for each group: control and infected). The legend in the heatmap indicates for every protein if their fold change values (both in the proteome and in the transcriptome) are positive (dark blue) or negative (light blue). (B) Average TPM values are represented for the whole transcriptome and for the differentially expressed contigs, as well as for the equivalent contigs to the proteins quantified in the proteome and the DEPs. How the proteins detected in the proteome correspond to contigs with higher average expression than the overall transcriptome (it is more possible to identify in proteomics a term with high transcriptomic expression than one with low expression) can be seen. (C) Correlation analysis between the log 2 -fold change values of the proteomic quantification and the transcriptomic expression of the DEPs. From the total number of 226 DEPs, 87 were discordant gene-protein pairs, indicated in pink (different trend between both approaches), and 139 (61%) were concordant gene-protein pairs, indicated in blue (same expression trend between both approaches). The R-squared value when considering all the DEPs (gray line) and when considering the concordant gene-protein pairs (blue line) are represented, as well as the p-value of the fold change value correlation for the concordant gene-protein pairs. The same correlation analysis was performed considering only the immune-related DEPs, the same ones that were used in (A). How the correlation increases with certain types of modulated genes (immune functions) can be seen, while it is highly decreased in other functions (cytoskeletal-related DEPs, with R 2 = 0.0641; data not shown). (D) Correlation analysis between the log 2 -fold change values of the proteomic quantification and the transcriptomic expression for all the quantified proteins in the proteome. The regression analysis considering the whole proteome was not significant. From the total number of 4,242 quantified proteins, 2,105 were discordant gene-protein pairs (pink), while 2,137 were concordant gene-protein pairs (blue). All the concordant gene-protein pairs (blue dots) were used for the represented regression (blue line), and the R 2 and p-value are indicated for this analysis with concordant terms. since the carbonylation and glutathionylation of this protein have been related in Mytilus edulis with a recognition system of oxidative stress (McDonagh et al., 2005). It has been reported that "decreased actin dynamics cause depolarization of the mitochondrial membrane and an increase in reactive oxygen species (ROS) production" (Gourlay et al., 2004). In accordance with these previous reports, our proteomic study on infected mussel gills also reported an important response related to oxidative stress along with cytoskeletal modulation. Upmodulated cytochrome c oxidase subunit 12 and methionine-R-sulfoxide reductase B2 proteins could also be implicated in the regulation of ROS. Oxidative proteins such as dual oxidase were downregulated, as were antioxidant proteins such as hydroperoxide glutathione peroxidase. Additionally, related to oxidative stress situations is the alteration of calcium homeostasis (Görlach et al., 2015). The modulation of oxidative/antioxidant and calcium-related proteins is common in previous bivalve proteomic studies, not only under infection conditions but also against stress stimuli of different kinds (Manduzio et al., 2005;Ji et al., 2013b;Zhang et al., 2015;Chen H. et al., 2018;Franco-Martínez et al., 2018). Calmodulin is an upregulated DEP that controls the function of numerous enzymes, ion channels and other proteins through Ca 2+ regulation. The Ca 2+ flux regulated by calmodulin is necessary for NOS production in the immune response (Tagliazucchi and Conte, 2005;Ma et al., 2008). However, Striatin-3, an associated protein that directly binds calmodulin and is responsible for the activation of endothelial NO synthase (Lu et al., 2004), was downregulated. The kinase suppressor of Ras Ksr1 was also downregulated in infected samples. Ksr1 is responsible for the assembly of inducible nitric oxide (NO) synthase (iNOS) upon bacterial infections (Zhang et al., 2011). NO production would therefore be inhibited to avoid increasing oxidative stress conditions. Proteins related to the regulation of STIM1-dependent calcium entry (Surfeit locus protein 4 and EF-hand domain-containing family member B) were downregulated as well.
As mentioned before, NF-kB regulation was the key response detected in the analysis to the infective situation caused by the bacteria. The results showed how many of those modulated proteins could interact with each other in what would be an NF-kB regulatory network, also highlighting associations with other DEPs implicated in different functions. The regulation of the NF-kB inflammatory response is essential for the immune response; therefore, it is conserved from ancestral molluscs as mussels to the innate component of the immune response in superior vertebrates. Several pathogen recognition proteins that could recognize the bacteria and initiate this response were modulated in the transcriptome (Saco et al., 2020) and were identified but not modulated in the proteome: fibrinogen-like proteins or FREPs, different types of lectins, C1-q containing proteins, and allograft inflammatory factor-1 (AIF-1). AIF-1 is a pro-inflammatory cytokine implicated in immune activation in molluscs (De Zoysa et al., 2010;Zhang et al., 2013) and can activate the NF-kB response (Egaña-Gorroño et al., 2019;Zhou et al., 2019). Immune recognition and signaling terms dominated the modulated response in the transcriptomic study of the same infection model, but they were not modulated in the current proteomic approach. The lack of concordance between both approaches concerning this aspect of the response can be due to the different depth of each technique, to the different technique that is used to detect each type of molecule (proteins and transcripts) or to biological reasons.
Despite this lack of concordance, the modulated DEPs completely covered the entire pathway of the NF-kB response. Downmodulated Gelsolin-like protein 1 and Gelsolin-like protein 2 are calcium-regulated and actin-associated proteins. Gelsolin proteins can bind LPS from the gram-negative Vibrio membrane, resulting in neutralization of LPS effects, including cytoskeletal actin remodeling and NF-kB response (Bucki et al., 2005(Bucki et al., , 2008. The clear downmodulation of gelsolin would support the responses observed in the proteome. In an unstimulated context, NF-kB is captured in the cytoplasm by IκB proteins (inhibitors of NF-kB). The NF-kB activation pathway is activated after pathogens are recognized (TLR recognition) or inflammatory signals arrive at the cells (inflammatory cytokines such as AIF-1 and IL-17) (Liu et al., 2017). Then, the NF-kB pathway is activated by different transductor proteins also identified in the proteome, such as MyD88 and TRAF proteins. Scaffold mitochondrial proteins implicated in pathway transduction were modulated, such as NIP-SNAP and MCC1 (Cao et al., 2016;Yamamoto et al., 2017). The phosphorylation of TAK-1 is another signaling step, and it is enhanced by the upregulation of the 14-3-3 protein, which is directly associated with the upregulation of calmodulin (Luk et al., 1999;Somboonwiwat et al., 2010;Zuo et al., 2010).
The transduction pathway results in the phosphorylation and subsequent activation of IκB kinase (IKK). IKK phosphorylates two serine residues of the IκB protein, leading to ubiquitination and degradation by the proteasome and releasing NF-kB. This crucial step is activated in the proteomic response, since proteins that would enhance the process are upregulated (PACRG and the associated proteins STING and TRIM32) (Albor et al., 2006;Dunphy et al., 2018;Meschede et al., 2020) and inhibitory proteins of the phosphorylation and degradation of IκB are all downregulated: SAMHD1, RASSF2, β-Arrestin (Witherow et al., 2004;Song et al., 2012;Chen S. et al., 2018). The upregulated ubiquitin-conjugating enzyme E2 is implicated in the elimination of phosphorylated IκB in the proteasome, and its implication in the NF-kB response has already been reported in bivalves (Hatada et al., 2000;Liu et al., 2019). Released NF-kB is transported to the nucleus, and this process in certain cell types is reported to be carried by a dynein complex (Mikenberg et al., 2007). Translocation of released NF-kB into the nucleus is regulated by importin α, which is also upregulated in the proteome, in a nuclear localization signal (NLS)-dependent process (Fagerlund et al., 2005(Fagerlund et al., , 2008. Importin α is an adapter protein of KPNB1, the major nuclear receptor for NF-kB import (Liang et al., 2013).
Once NF-kB has been translocated into the nucleus, we found modulation of the Pirin protein, a transcriptional coregulator of NF-kB. Pirin facilitates NF-kB binding to target genomic regions in a redox-state-dependent manner: it is more active when oxidized in its ferric form (Liu et al., 2013) under oxidative conditions as the infection would unleash. Modulation of Pirin quantity has also been associated with NF-kB activity inside the nucleus (Brzóska and Kruszewski, 2010;Licciulli et al., 2010). Two proteins implicated in the inhibition of NF-kB transcriptional activity are downregulated in the proteome: the E3 ubiquitin-protein ligase TRIM39 and 5 -3 exoribonuclease 2 (XRN2) (Rother et al., 2016;Suzuki et al., 2016).
Proteins corresponding to typical downstream genes whose expression can be regulated by NF-kB were identified as well: HSPG2, Fibulin-1, IF44L, and GVIN1. These proteins are implicated in effector functions against infection. Myticins were also identified in the proteome. Myticins are antimicrobial peptides showing cytokine-like functions and are the most highly expressed genes in mussel defensive cells and hemocytes (Mitta et al., 1999;Balseiro et al., 2011). Antimicrobial peptide expression in bivalves has been demonstrated to be regulated by NF-kB (Oyanedel et al., 2016). More effector proteins were found to be modulated in the experiment. Adhesion-related proteins, commonly implicated in bivalve immune responses and related to cell migration at the site of infection, were modulated: fibulin-1, hemicentin, cell adhesion molecule 3, lachesin (Gueguen et al., 2003;Xu and Faisal, 2008;Xu et al., 2018) and filamin A, which are also implicated in NF-kB and ROS responses (Muscolini et al., 2011;Uotila et al., 2017). Phagocytosis-related Tubby protein homolog and granulin are also implicated in effector actions against infection (Caberoy et al., 2010;Wu et al., 2018;Hambrook et al., 2019).
The proteomic approach contributed with a complementary vision to the transcriptomic results, adding a functional validation of the processes indicated by the modulation of gene expression. Proteomic results described effects caused by bacterial infection in the host cytoskeleton that were previously undetected in the transcriptomic analysis. Furthermore, a net of proteins acting together was described, linking cytoskeletal disruption with the main modulated responses: the regulation of oxidative stress and the inflammatory NF-kB axis. The induction of the NF-kB response after the pathogen recognition steps is in agreement with previous results at the transcriptomic level in a similar experiment. In fact, the correlation between proteomic quantification and transcriptomic expression was higher within these immune DEPs, validating this response at two different levels of gene expression. Discrepancies were also observed, since the modulated response with each approach was based on different terms that would take part in the same response pathway: the transcriptome reported receptors and cytokines that would activate the NF-kB response, while in the proteome, a network of proteins that would activate the intracellular pathway was obtained. Because of this, an enriched and more complete overall vision of the mussel response to this natural infection model was obtained with the combination of both approaches.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (Deutsch et al., 2020) via the PRIDE (Perez-Riverol et al., 2019) partner repository with the dataset identifier PXD028043 and 10.6019/PXD028043.

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 the 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. APD, AS, SB, and AP performed the proteomic analyses. APD, AS, and AF performed the bioinformatic analyses. AS wrote the manuscript. All authors contributed to the article and approved the submitted version.