Development of a Gill Assay Library for Ecological Proteomics of Threespine Sticklebacks (Gasterosteus aculeatus)

A data-independent acquisition (DIA) assay library for quantitative analyses of proteome dynamics has been developed for gills of threespine sticklebacks (Gasterosteus aculeatus). A raw spectral library was generated by data-dependent acquisition (DDA) and annotation of tryptic peptides to MSMS spectra and protein database identifiers. The assay library was constructed from the raw spectral library by removal of low-quality, ambiguous, and low-signal peptides. Only unique proteins represented by at least two peptides are included in the assay library, which consists of 1506 proteins, 5074 peptides, 5104 precursors, and 25,322 transitions. This assay library was used with DIA data to identify biochemical differences in gill proteomes of four populations representing different eco- and morpho-types of threespine sticklebacks. The assay library revealed unique and reproducible proteome signatures. Warm-adapted, low-plated, brackish-water fish from Laguna de la Bocana del Rosario (Mexico) show elevated HSP47, extracellular matrix, and innate immunity proteins whereas several immunoglobulins, interferon-induced proteins, ubiquitins, proteolytic enzymes, and nucleic acid remodeling proteins are reduced. Fully-plated, brackish-water fish from Westchester Lagoon (Alaska) display elevated ion regulation, GTPase signaling, and contractile cytoskeleton proteins, altered abundances of many ribosomal, calcium signaling and immunity proteins, and depleted transcriptional regulators and metabolic enzymes. Low-plated freshwater fish from Lake Solano (California) have elevated inflammasomes and proteolytic proteins whereas several iron containing and ion regulatory proteins are reduced. Gills of fully-plated, marine fish from Bodega Harbor (California) have elevated oxidative metabolism enzymes and reduced transglutaminase 2, collagens, and clathrin heavy chains. These distinct proteome signatures represent targets for testing ecological and evolutionary influences on molecular mechanisms of gill function in threespine sticklebacks. Furthermore, the gill assay library represents a model for other tissues and paves the way for accurate and reproducible network analyses of environmental context-dependent proteome dynamics in complex organisms.


Introduction
Enormous efforts have been devoted towards development of quantitative proteomics approaches that enable large network-scale analyses of proteome dynamics. As a result, it is now possible to routinely assess systems-scale protein dynamics with superior quantitative accuracy, e.g. for quality control of antibody selectivity and specificity (1). An exciting recent development is data-independent acquisition (DIA) technology, which permits precise and simultaneous quantitation of thousands of proteins (2,3).
This technique is well suited for label-free approaches to be used with virtually any species studied in its natural environment (4,5). However, application of these technical developments for ecological proteomics is not common and its potential is underutilized. One of the main technical challenges and time-consuming aspects for precise DIA quantitative proteomics is the selection and validation of diagnostically reliable transitions that unambiguously and quantitatively represent unique proteins in diverse samples (6)(7)(8). In this study we have addressed this challenge by developing a reliable, standardized, and manually curated DIA assay library assay for gills of threespine sticklebacks (Gasterosteus aculeatus). DIA targeted data extraction has recently been introduced for biomedical applications and shown to be as quantitatively consistent and accurate as selected reaction monitoring (SRM), which is considered the gold standard of quantitative proteomics (3,7,8). DIA-based assays combine advantages of selective/multiple reaction monitoring (SRM/ MRM) approaches on the one hand and large-scale, high-throughput DDA on the other hand. Like SRM/ MRM approaches, DIA methods monitor well-defined transitions representing specific precursors and proteins, although, in DIA the transitions are targeted post-acquisition by utilizing MSMS spectral libraries (9). In addition, DIA approaches are suitable for quantitation of more than a thousand proteins simultaneously, which is why they are sometimes referred to as hyper-reaction monitoring (HRM) (10). The utility of the stickleback gill assay library developed in this study was demonstrated by comparing gill proteomes of four stickleback populations adapted to four different environments. Quantitative reproducibility of this assay library was verified by analyzing four different population samplings for each of the four habitats.

5
Threespine sticklebacks are euryhaline fish that represent a long-standing model in behavioral, evolutionary, and ecological physiology research (11). These fish are ancestrally marine but the species is also represented by many anadromous and land-locked populations that have colonized brackish and freshwater habitats following the retreat of Pleistocene glaciers (12)(13)(14). This species is eurythermal but typically restricted to cold-temperate climates throughout coastal areas in the entire Northern hemisphere (15,16). However, some exceptional refuge populations have evolved to colonize warm-water habitats Threespine sticklebacks ideally combine key advantages as a model organism for integrative biology and cross-disciplinary research at the interface of ecophysiology, behavior, evolutionary biology, biochemistry and genetics. They are represented by many populations that have naturally adapted to diverse environments, their complete genome has been sequenced, and their proteome is well annotated (14). A sequenced and well-annotated genome represents a major resource for studying the biology of a particular species (19). Therefore, much attention has been devoted towards developing and applying whole genome sequencing approaches and strategies for genome, transcriptome, and proteome annotation to an increasing number of species (20). The ability to unambiguously map peptide and protein sequences based on a well-annotated species-specific reference database to unique genomic loci represents a critical prerequisite for proper proteome-wide and unambiguous quantitative analyses. In this study we have utilized these unique resources and advantages of threespine sticklebacks to develop a strategy for longterm quantitative ecological proteomics of this species. This strategy is based on the generation of validated assay libraries that are tissue-specific and can be used in a highly consistent manner, i.e.
allowing the quantitation of the same set of proteins based on identical, previously validated precursors 6 and transitions in every sample. To demonstrate its usefulness, the current study has been limited to a single tissue (gill), but the approach can be expanded to include other tissues for integrative studies of environmental effects on whole organism proteomes. The availability of this gill assay library will propel precise, quantitative and network-scale analyses of proteome dynamics in fish experiencing a wide variety of ecological and evolutionary contexts.

Sample collection and preparation
Threespine sticklebacks (Gasterosteus aculeatus) were collected from four different habitats using seining and traps as previously described (4 and ecotype are summarized in figure S1. Additional phenotypes for each fish are provided in figure S2. To equally represent both sexes in the spectral library and assay library, three males and three females were included for each population in PS1. Gill dissection, protein extraction, and protein digestion with immobilized trypsin were performed using an optimized sample preparation protocol as described previously (25).

Data-dependent acquisition (DDA)
Tryptic peptides (2 μL, 300 ng/ μL) were separated along a 5 to 30% acetonitrile (ACN) gradient over 120 min using a nanoAcquity (Waters, Milford, MA) UPLC and injected into an ImpactHD UHR qTOF mass spectrometer (Bruker Daltonics, Bremen, Germany). The conditions used for DDA were as follows: Peptides were first trapped for one min (10 μL/ min flow rate) on a trap column (Symmetry, Waters 186003514) and then separated on a 1.

Data-independent acquisition (DIA)
DIA was performed on all samples collected in all four PSs (96 total biological replicates) using the same instrumentation and separation conditions as for DDA, nanoAcquity (Waters) interfaced online to an ImpactHD UHR qTOF (Bruker Daltonics). Slight variations in retention time (RT) were compensated for 9 by application of internal retention time (iRT) standard calibration in Skyline (figure one). DIA conditions reported previously were used (26) except that the isolation width was set to 10 m/z (0.5 m/z overlap) and spectra were collected over 75 isolation width windows ranging from 390 to 1065 m/z at 37.5 Hz and in two sec. intervals. Except for the acquisition frequency and operating the mass spectrometer in MRM mode, all acquisition parameters were identical to those used for DDA. DIA data from PS1 were used to filter the target list and for selecting iRT standards with Skyline 3.6 (21

Spectral library and DIA assay library construction
A G. aculeatus gill spectral library was generated using the peptide-to-spectrum matches and protein annotation information generated from PS1 DDA data. This information was exported in pepxml and Mascot dat formats and then imported into Skyline 3.6 (21). Skyline was used for spectral library construction and for generating an initial target list, which contained 400,412 transitions, 60,132 precursors, 50,501 peptides, and 5,998 proteins (figure two a-d). Transitions for the initial target list were automatically chosen from library spectra using the Skyline transition settings dialogue based on the following criteria: ion 3 to last ion -1; exclusion of the precursor isolation window (m/z width = 10); interferences. The rules followed for manual curation were as follows: 1. If an interference was greater than the actual peak at proper RT then the corresponding transition was removed; 2. If a precursor peak is not clearly discernable in at least half the samples of one population within three min of the predicted RT then the precursor was removed; 3. If precursor peak widths differ by more than two-fold for half of the samples of the same population from the other half, then the precursor was removed; 4. If more than one transition of the five nominal transitions shows an interference that exceeds the intensity of the actual transition peak at the predicted RT then the corresponding precursor was removed. The resulting final target list contains 25,322 transitions, 5,104 precursors, 5,074 peptides, and 1,506 proteins (figure one a-11 d). The complete assay library including all relevant metadata is available at Panorama Public (https://panoramaweb.org/labkey/JoLi-01.url). This assay represents a tier two assay (27).

Statistical analyses and visual data representation
Quantitative DIA data were normalized at the MS2 level by equalizing medians in each sample followed by comparison of protein quantities in each population against all other populations with the MSstats (22) module integrated in Skyline using the following parameters: confidence interval > 99% (FDR < 1%), scope = protein, summary method = Tukey's median polish, use zero for missing peaks, and q value < 0.005. The q value represents the FDR associated with the mProphet peak scoring model applied to the DIA data, i.e. mProphet FDR < 0.5%. The manual curation was the final step (step seven) after mProphet peak scoring during library assay generation. This final step was only performed for PS1 samples, which were used as training samples for steps five to seven of assay library construction (see previous paragraph). Once the assay library had been established using PS1 samples no manual curation was performed for subsequent DIA data sets (PS2, PS3, and PS4). Rather, the same assay library, parameters, and thresholds as used for PS1 were applied in conjunction with iRT calibration and mProphet peak scoring to quantify transition peaks for all subsequent DIA data sets (PS2, PS3, and PS4). Reproducibility of proteins identified as significantly different in abundance was assessed using two common multiple testing-corrected p value thresholds (0.01 and 0.05). This analysis demonstrates that a multiple testingcorrected p value threshold of 0.05 represents the best compromise for balancing false positives with false negatives (figure three). DIA data exported from Skyline were used for generating heat maps with Genesis 1.8.1 (28). The Venny 2.1 tool was used for visual representation of data as Venn diagrams (29).

Pathway and gene ontology analyses
The UniprotKB ID mapping tool was used to annotate stickleback gill proteins with their corresponding KEGG (Kyoto Encyclopedia of Genes and Genomes https://www.kegg.jp/) orthology (KO) identifiers.
KO identifiers were assigned to 1307 of the 1506 proteins contained in the final assay library (

Gill proteome identification and spectra annotation
Unbiased DDA proteome profiling of all samples collected in PS1 followed by PEAKS 8.0 13 The Scaffold file also contains semi-quantitative spectral count analyses that suggest candidate proteins whose abundances are population-specific (MassIVE AC MSV000081795, CAMP proteome AC CAMPDDA00023). The PEAKS Top3 and Scaffold spectral counting results obtained from PS1 indicate two broad trends: First, gill proteomes differ markedly between the four populations. Second, the gill proteomes of the LaBoRo and WesLag populations are most unique compared to the other populations.
For the purpose of this study the DDA data were used for assay library generation and not for systematic comparison of top3 and spectral counting approaches with assay library quantitation, which exceeds the scope of this article.

Assay library development
DIA targeted data extraction was first introduced six years ago (3). This method eliminates common problems with protein quantitation by DDA (stochastic, irreproducible selection and undersampling of precursor ions) while greatly exceeding the scope (number of proteins) that can be quantified by targeted approaches (SRM, MRM, PRM). Since its inception, many variants of DIA technology have been employed. Two major variants include SWATH on triple-quad instruments and DIA on UHR-qTOF instruments (9,(31)(32)(33)(34)(35). Recently, DIA targeted data extraction has been advanced and quantitative assay precision improved by providing tools that permit target list filtering and generation of curated target lists based on stringent quality control criteria, some of which are evaluated using a "training sample set". In this study we have used Skyline (21) for targeted DIA data extraction and filtering the target list. The resulting assay library is much shorter than the initial raw target list generated from spectral libraries.
Quantitative precision of the assay library is greatly improved because protein assignment ambiguities, low-quality and low-signal transitions, and proteins with less than two unique precursors were eliminated during target list filtering. However, target list filtering also eliminates proteins that are potentially biologically interesting and detectable by multiple peptides if not at least two of these peptides meet all filter criteria. Automatic target list filtering with Skyline eliminated the vast majority of low-quality transitions. Moreover, utilizing all 24 samples from PS1 as a "training sample set" for several automatic steps and a final manual step of target list filtering further eliminated unreliable transitions with ragged peak shapes, interferences, inconsistent peak outlines or low signal-to-noise ratio.
The assay library permits large-scale (network-enabling) analyses of proteome dynamics in gills of threespine sticklebacks exposed to variable ecological contexts. A current limitation is that less than 10% of the overall proteome can be assayed with this assay library. Although it is unlikely that all proteins encoded by the genome are expressed in gills, even greater sensitivity and dynamic range of future instrumentation will gradually alleviate this limitation. The conceptual approach for assay library construction used in our study consisted of first acquiring high-quality fragment ion spectra in DDA mode, second, assigning peptide sequences to these spectra and annotating them to specific proteins, third, generating annotated spectral libraries, and fourth, compiling the final assay from the spectral library. This approach follows a recently described protocol (7,8). However, here we use PEAKS Studio rather than the transproteomics pipeline (TPP) for consolidating multiple search engine results to annotate MSMS spectra from DDA data and Skyline rather than OpenSWATH (36) for spectral library construction and assay library filtering. Schubert and colleagues filtered out fragments smaller than 350 m/z or greater than 2,000 m/z as well as precursors that fall within the precursor isolation window from their assay library and they used a mass accuracy threshold within ±0.05 m/z of the expected mass (8).
We used similar criteria for the generation of the initial target list via the transition settings dialogue in Skyline. In step 1 of reducing the initial library to the final assay library the five most intense y and b ions were selected (min. fragment ion = ion 3, e.g. y4 even if y2 and y3 are more intense, max. fragment ion = last ion -1 excluding the precursor isolation window, m/z threshold ±0.035 range) (Figure one a-d).
Schubert and colleagues recommend that all assays in the final assay library have the same number of transitions to avoid mixed distributions that complicate statistical analysis. However, our inclusion of some precursors represented by four (rather than the nominal five) transitions in the assay library is accounted for since MSstats employs mixed linear models for statistical analysis (22). Removal of transitions that produce interferences in most or all training samples improves quantitative precision without impairing protein coverage of the assay library. In addition, four processed transitions still provide the highest separation power of the mProphet score (23).

Quantitative reproducibility of the stickleback gill assay library
The reliability and quantitative reproducibility of the assay library was very high as assessed by analysis of three additional PSs, each consisting of 24 samples that were not included in the "sample training set" derived from PS1. No manual curation of the assay library was performed for any samples in PS2, PS3, and PS4 and all targeted data extraction procedures were done using Skyline automated workflows, including internal retention time alignment, mProphet peak scoring, and MSstats quantitative analysis (21-23). Visual inspection verified that transition peaks were properly assigned and outlined for precursors of proteins that were identified by MSstats analysis as significantly elevated or reduced in a particular population.
Since all raw data, quantitative results, and metadata are conveniently accessible in Panorama public (37) the corresponding Skyline files can be downloaded directly to reproduce all quantitative and statistical results. Although the population-specific abundance patterns of many proteins were very consistent in all four PSs, there were also differences for many proteins. Most of these differences were small. They can be attributed to three causes: First, individual differences in the expression of at least some proteins are expected between different fish. Second, the level of elevation or reduction of some proteins is slightly different, i.e. just above the significance threshold in some and slightly below that threshold in other PSs. DIA datasets analyzed with the gill assay library represent a consistent matrix for powerful network modeling approaches, which can be complemented in the future by datasets for additional samples to enable comprehensive modeling of ecological and evolutionary effects on gill proteome dynamics.
Although limited to a subset of the proteome (1506 proteins), the gill assay library allows robust and reproducible quantitation of a rich repertoire of striking population-specific proteome signatures. In contrast to genomes that are stable within a single individual, the corresponding proteome is extremely dynamic depending on developmental, ecological, and cell/tissue contexts (38). Therefore, species-and tissue-specific proteomes can be considered excellent "reflections" of an organism's life history exposures and experiences that provide information about its exposome (39,40). Biological implications of the gill proteome signatures identified in this study are discussed below. An analysis of whether these proteome signatures represent the result of evolutionary/ genetic adaptations or physiological plasticity/ acclimation exceeds the scope of this study. Nevertheless, the identification of pertinent proteome signatures and the establishment of a gill assay library enables such future studies.

Ecological significance of population-specific proteome signatures
Pathway and GO enrichment approaches represent convenient tools for surveying broad biological implications of large datasets (41). We have utilized such approaches previously for threespine sticklebacks to interpret large proteomics datasets and condense the associated information (4, 5).
However, only known pathways can be analyzed pathway/GO enrichment approaches and representation of pathways and GO terms in existing databases is highly uneven. Protein annotation with GO terms and pathways is very heterogeneous ranging from none to hundreds per protein. Moreover, functional annotations are often based on only a handful of canonical model organisms studied in particular contexts.
Thus, enrichment analyses are biased against poorly annotated pathways/ GO terms that may be more novel and unique in the particular tissue, species, and contexts studied. In addition, because the assay library covers less than 10% of the stickleback proteome the statistical power of enrichment analyses is lower compared to transcriptomics studies that utilize the entire genome for this purpose. Because of these limitations and caveats of such enrichment analyses we have complemented them with manual, literature-based assessment of pertinent functions of highly significant proteins.

Proteome signature of the Bodega Harbor population
Marine sticklebacks from Bodega Harbor have the least unique gill proteome of all populations analyzed.
A possible reason for this observation is the evolutionarily ancestral state of marine threespine stickleback populations (16), which may include most of the traits displayed by more specialized, derived populations at levels high enough to reflect a "species average". Despite the lack of significant enrichment/ depletion of any PANTHER pathway/ GO term, a distinct proteome signature is evident even  (46), and tissue damage repair (47). Overall, the BodHar gill proteome signature is consistent with the eco-and morpho-type of the BodHar population and suggests specific molecular targets that might underly enhanced growth and SW adaptation of BodHar fish.

Proteome signature of the LaBoRo population
The LaBoRo population is unique among threespine sticklebacks since it chronically experiences temperatures of at least 30°C, which are lethal or highly detrimental for most other stickleback populations (48). Elevated LaBoRo proteins cluster in KEGG pathways that indicate enhanced roles of the extracellular matrix, cell adhesion, the plasma-membrane anchored cytoskeleton, and the complement system. Although many isoforms of common HSPs (HSP70, HSP60, HSP90) are included in the assay library none of them are elevated in the LaBoRo population. This unexpected finding may indicate that 24 these molecular chaperones are not limiting for proteostasis at chronically elevated temperature.
However, two HSP47 (serpin 1) isoforms are strikingly elevated in gills of LaBoRo sticklebacks. HSP47 is a molecular chaperone for collagens and some other ECM proteins that promotes their maturation during thermal stress (49). Although collagens are not among the elevated ECM proteins in LaBoRo fish many other ECM and cell adhesion proteins are highly elevated suggesting a major role of extracellular matrix stabilization in LaBoRo stickleback gills. A key role of HSP47 for ECM reorganization has also been demonstrated in animal models of human diseases, including cardiac infarction (50), tumorigenesis (51,52) and renal glomerulonephritis (53). Moreover, HSP47 plays critical roles for skeletal growth, ECM patterning, and fin regeneration in zebrafish (54). Elevated Hsp47 mRNA levels are correlated with increased temperature and heat-shock in many organisms (49,55,56), including fish (57), suggesting that the high habitat temperature causes severe ECM strain in gills of LaBoRo sticklebacks.
KEGG pathways that are most prominent among reduced LaBoRo proteins indicate a suppression of proteasome, lysosome, and ribosome functions as well as protein anabolism, adaptive immunity, and necroptosis. In addition, RNA/ DNA binding proteins and translation regulators/ factors are significantly enriched among reduced LaBoRo proteins supporting the notion of suppressed protein turnover in the LaBoRo population. A striking proteome feature of LaBoRo sticklebacks is the strong depletion of proteins involved in adaptive immunity while several complement proteins involved in innate immunity are elevated. Fish are ectotherms and thermal stress effects on the vertebrate complement system have been documented previously (58,59). It should be pointed out that the fish gill is a highly vascularized tissue and that proteins detected in this study may be localized in branchial epithelial, vascular, and/ or extracellular compartments (60). The depletion of adaptive immunity proteins is exemplified by enormously reduced levels of interferon-induced GTP-binding protein Mx. This protein (Mx) is a very potent antiviral defense protein that protects cells from many RNA and DNA viruses (61,62). Mx interacts with calcium signaling (63) and proteostasis pathways (e.g. the ER stress response) (64).
Although causality between molecular LaBoRo proteome signatures and temperature is currently unclear, this population chronically experiences temperatures of at least 30°C, which are lethal or highly detrimental for most other stickleback populations (48). In addition, strong temperature effects on the immune system have been reported for threespine sticklebacks (65) and other fish species (66). Whether caused by high habitat temperature or other environmental factors, the unique signature of immune system proteins in LaBoRo fish suggests that their adaptive immune functions are compromised.
Induction of innate immunity via the complement system may represent a compensatory adaptive strategy. Alternatively, activation of the complement system may indicate increased levels of cellular and tissue-damage under the stressful habitat conditions experienced by LaBoRo sticklebacks. Such increased damage would exacerbate the demand for complement clearance of cellular debris (66,67).

Proteome signature of the WesLag population
The PANTHER integrin signaling pathway is significantly depleted from the elevated LaBoRo and WesLag proteins. This finding illustrates that this pathway is well represented by the 1506 assay library proteins but not elevated in either population. Significant enrichment of hydrolase and phosphatase activity among elevated WesLag proteins is exemplified by several highly elevated ATPases that are central for ion transport. The ribonucleoprotein complex is also enriched among elevated WesLag proteins, which cluster into calcium and GTPase signaling, actin cytoskeleton regulation, cell adhesion, pattern recognition, and protein translation KEGG pathways. The most highly elevated WesLag protein is chloride channel protein 2 (CLCN2). This ion channel and the elevated cation ATPases are regulated by environmental salinity in euryhaline fish, including sticklebacks (13, 60,65,68). We interpret these findings as evidence that WesLag fish experience frequent salinity fluctuations in their habitat, either because of tidal action or because of migrations between freshwater, brackish, and marine habitats. This anadromous population may have migrated just prior to collection from seawater to Westchester Lagoon to breed. Intracellular calcium, which controls many physiological processes of euryhaline fish including osmoregulation, is also regulated by the Ca 2+ -ATPase (69). Besides calcium signaling, a prominent signal transduction feature of the WesLag-specific gill proteome signature are various GTPases, including several Ras-related small GTPases that are known to regulate structural stability of the cytoskeleton and epithelial homeostasis and barrier function (70)(71)(72). These processes are key targets during cellular osmoregulation (73).
PANTHER enrichment of amino acid metabolism and blood coagulation among reduced WesLag proteins is accompanied by clustering of these proteins within endocannabinoid and PPAR signaling, spliceosome, lysosome, antigen processing and presentation, necroptosis, glutathione metabolism, and pentose phosphate KEGG pathways. Suppression of these processes in the WesLag population is Although not present in the assay library, the NLRP sensor of LakSol stickleback gill inflammasomes is likely NLRP3 because this sensor preferentially interacts with the PYCARD adaptor (75,76). The PYCARD adaptor activates the effector caspase 1 leading to proteolytic cleavage of an inactive cytokine precursor into active, proinflammatory interleukin-1 beta (74,(77)(78)(79). Inflammasome induction represents a response to parasitism in many species, including fishes (80)(81)(82). Freshwater sticklebacks are frequently parasitized (83) and we did observe cestode parasites (Schistocephalus solidus) in the body cavity of some LakSol sticklebacks during dissection, but never in any of the other populations. We excluded heavily parasitized fish during dissection but defense mechanisms against parasitism may be even more effective in fish from the same habitat that lack visible signs of parasitism. Potent inflammasome induction has also been reported in response to mercury exposure (84,85). This finding is significant because Lake Solano is characterized by high mercury concentration resulting from historical gold mining activities in California (86).
Intriguingly, cathepsins, which are also elevated in gills of LakSol sticklebacks, regulate the appearance and severity of mercury-induced inflammation (84). In addition, cathepsins are known to participate in inflammasome pathways in many other contexts (87-89).
Reduced LakSol proteins are depleted in PANTHER nucleic acid and RNA binding proteins suggesting that the levels of these proteins are maintained similar to those in other populations. Instead, reduced LakSol proteins cluster in focal adhesion, ECM-receptor interaction, tight junction, PI3K-Akt signaling, actin cytoskeleton regulation, and oxidative metabolism KEGG pathways. A reduction of oxidative metabolism is supported by a decrease in hemoglobin isoforms. Decreases in hemoglobin and erythrocytes have been observed as a result of parasitism in freshwater fish (90). For instance, the Brazilian cultivated sport fish Leporinus microcephalus has a reduced hematocrit, erythrocyte count, and hemoglobin concentration caused by the parasite Goezia leporini (91). Depletion of hemoglobins, iron transporting ferritins, and other plasma proteins may also be a consequence of parasitism in LakSol sticklebacks, but this hypothesis remains to be tested. As expected, proteins responsible for active NaCl secretion (e.g. Na + /K + /2Clcotransporter) are reduced in LakSol fish indicating that their gills operate in ion absorption mode (60,92) in their FW habitat.

Conclusions and perspective
To

Data Availability
Raw mass spectrometry data, metadata, and processed data are available at the following public repositories: