Protein changes as robust signatures of fish chronic stress: a proteomics approach to fish welfare research

Background Aquaculture is a fast-growing industry and therefore welfare and environmental impact have become of utmost importance. Preventing stress associated to common aquaculture practices and optimizing the fish stress response by quantification of the stress level, are important steps towards the improvement of welfare standards. Stress is characterized by a cascade of physiological responses that, in-turn, induce further changes at the whole-animal level. These can either increase fitness or impair welfare. Nevertheless, monitorization of this dynamic process has, up until now, relied on indicators that are only a snapshot of the stress level experienced. Promising technological tools, such as proteomics, allow an unbiased approach for the discovery of potential biomarkers for stress monitoring. Within this scope, using Gilthead seabream (Sparus aurata) as a model, three chronic stress conditions, namely overcrowding, handling and hypoxia, were employed to evaluate the potential of the fish protein-based adaptations as reliable signatures of chronic stress, in contrast with the commonly used hormonal and metabolic indicators. Results A broad spectrum of biological variation regarding cortisol and glucose levels was observed, the values of which rose higher in net-handled fish. In this sense, a potential pattern of stressor-specificity was clear, as the level of response varied markedly between a persistent (crowding) and a repetitive stressor (handling). Gel-based proteomics analysis of the plasma proteome also revealed that net-handled fish had the highest number of differential proteins, compared to the other trials. Mass spectrometric analysis, followed by gene ontology enrichment and protein-protein interaction analyses, characterized those as humoral components of the innate immune system and key elements of the response to stimulus. Conclusions Overall, this study represents the first screening of more reliable signatures of physiological adaptation to chronic stress in fish, allowing the future development of novel biomarker models to monitor fish welfare.


Background
Managing welfare of fish in captivity is of increasing importance, both for productivity and sustainability reasons [1]. There is still no clear consensus on how welfare should be defined or objectively measured, given the complexity and controversy of the concept [2, 3]. Challenges like divergent coping mechanisms, the incomplete knowledge regarding the nociceptive system of fish (e.g. emotional-like states; cognitive abilities, pain, suffering) [4][5][6][7] and the lack of reliable physiological indicators of fish welfare, make its investigation even more difficult [8].
An aquaculture rearing facility deals with multiple stressful situations (stressors) that are inherent to daily routines and can compromise the fish well-being. These situations are usually unpredictable and uncontrollable for the animal and can range in duration and severity [8,9]. Fish launch a physiological response when faced with these threatening situations [9,10]. This adaptive mechanism, known as stress response, involves a cascade of reactions and enables the fish to cope with the stressor. However, when a stressful event is repeated or prolonged, it exceeds the organism's natural regulatory capacity and the fish fails to regain homeostasis, consequently impairing welfare [11,12].
The physiological stress response starts with the immediate activation of the sympathetic response, followed by a slightly delayed activation of the hypothalamopituitary-interrenal (HPI) axis. As a result, catecholamines and corticosteroids (cortisol in teleosts), respectively, are released into the bloodstream [13,14]. These hormones lead to a series of downstream responses involving alterations in the energy metabolism and respiratory and immune functions [15]. The rapid mobilization of energy substrates such as glucose (the fuel needed for the coping mechanisms) is caused by the activation of both the glycogenolysis in the liver or muscle, and the hepatic gluconeogenesis, by the catecholamines and cortisol, respectively [16,17]. Stressful stimuli can also lead to strenuous exercise fuelled by anaerobic glycolysis in the muscle, generating lactate, which is then released into plasma [18,19]. Prolonged exposure to the stressor will inevitably lead to alterations that are reflected in the whole-animal's performance, like perturbations at the reproduction, immunological, growth and behaviour levels [20].
The plasma levels of cortisol, alongside glucose and lactate, are the most commonly used physiological indicators to assess stress in fish [21]. Nevertheless, some inconsistencies have been reported in several experimental studies. This demonstrates the unreliability of these indicators, mainly in cases of chronic stress, which is mostly due to: (i) a high variability of response levels; (ii) a decrease of the cortisol levels to basal levels within minutes/hours following an acute stressor; (iii) the fact that fish can adapt, to certain extent, to chronic stress and the cortisol response is therefore attenuated; and (iv) the intrinsic and extrinsic factors (e.g. age, sexual maturity, social status, level of domestication, prior experience, nutritional status) that can affect cortisol secretion [22][23][24][25][26][27]. In this sense, it is vital to complement the existing behavioural, biochemical and physiological measures for a correct interpretation of the welfare status of the fish. This will be crucial to form a robust welfare assessment and to allow, in the future, the development of targeted recommendations and legislation. With the increasing research into the welfare of cultured fish, more advanced technologies are gaining popularity. Proteomics are promising alternatives for the discovery of candidate molecular markers that can indicate physiological alterations due to stress exposure [28]. Despite the limitations to the use of these technologies in the aquaculture field [29], several studies prove already the huge potential of proteomics for the identification of stress signatures [30][31][32][33][34].
There is very little data available concerning the process of long-term coping with a chronic stressor and indicators used in this case. Considering this gap in research, we aim, in the present study, to comparatively assess the stress responses of fish at different levels (i.e., plasma stress markers, changes in plasma proteins' abundance and muscle biochemistry). Using Gilthead seabream (Sparus aurata) as model, three chronic stress conditions were employed, and proteomics was used to benchmark potential signatures of stress adaptation in the plasma proteome since several proteins resulting from physiological events are released into circulation. Gilthead seabream was the chosen species in this study since it is one of the most important species in European aquaculture with high commercial value. This work aims to pioneer a better understanding of the underlying molecular mechanisms behind the fish physiological adaptation to long-term stress. Additionally, it aims to bridge the gap between the scientific community and the industry by paving the way for the development of novel biomarkers to monitor fish welfare.

Fish general condition
Fish were monitored every day during the trials. The experimental periods reached the end with a 100% survival rate. The overall condition and growth performance of the fish were also monitored (Additional file 1), and initial (IBW) and final body weights (FBW) were recorded for each experiment. The average body weight was reduced by the end of the net handling (NET) and hypoxia (HYP) trials, in all groups, including the control. However, there were no significant differences in final body weights between the control group and any of the stressed groups (P > 0.05), suggesting that weight reductions were unrelated to the stressor.

Plasma stress markers analysis
Circulating cortisol, glucose, and lactate levels were measured in Gilthead seabream submitted to different chronic stressors and in control fish (Fig.1). The overall levels of these metabolites showed a high variability of biological responses in all trials, with several data points considered outliers (outside of the interquartile range). Cortisol levels presented the highest intervals of values. In the OC trial, only lactate plasma levels presented statistically significant differences, both between control and stressed groups (Lactate CTRL -13. 46
The onset and resolution of rigor mortis (Fig.2) showed significant differences between treatments in the NET and HYP trials, specifically at 8 HAD (P CTRL-NET4 < 0.001), and at 8 (P HYP30-HYP15 < 0.001) and 24 HAD (P HYP30-HYP15 = 0.020), respectively. In the OC trial, fish reached an average maximum rigor strength at 24 HAD. In the NET trial, averaged maximum rigor strength was reached at 48 HAD in CTRL and NET2, and at 24 HAD in NET4 group. In the HYP trial, all groups reached averaged maximum rigor strength at 48 HAD.

Plasma proteomics analysis
A comparative proteomics analysis of the Gilthead seabream plasma between the control and the stress treatments detected, 681, 752 and 681 protein spots for the OC, NET and HYP trials, respectively, within the pH Fig. 2 Post-mortem changes in muscle pH and rigor mortis of gilthead seabream (Sparus aurata) submitted to different chronic stressors (aovercrowding, bnet handling, chypoxia), in two intensities, and unstressed fish (control), stored in ice for 72 h. Data points are the mean ± S.D. of n = 9 for each sampling time. Means labelled * are different at P < 0.05 range of 4-7 and a molecular mass range of 11-114 kDa. After statistical analysis, 19, 360 and 34 protein spots within the OC, NET and HYP trials, respectively, were found to present significantly differential abundance (significance threshold at P < 0.05) between experimental conditions. From these, 7, 171 and 12 were manually excised from the 2D gels for MALDI-TOF/ TOF MS analysis. No proteins were identified with significance for the OC trial. For the NET and HYP trials, 107 and 2 differential protein spots, respectively, were successfully identified by a combination of PMF and MS/MS search, with significant scores (protein score > 76, total ion score > 60, P < 0.05). Among the spots identified from the NET trial, 13 showed more than one significant protein identification (202, 326, 521, 559, 586, 604, 677, 877, 950, 959, 990, 996 and 1157), indicating that multiple proteins migrated to the same spots on the gel. The identified proteins are listed in an additional file (see additional file 2). A representative 2D-gel of the Gilthead seabream plasma proteome is shown in Fig. 3.
Hierarchical clustering (HCA) and principal component (PCA) analyses were performed for the identified 107 proteins spots with differential relative abundance across NET groups to check how well the samples Fig. 3 Representative pattern of gilthead seabream (Sparus aurata) blood plasma on a 12.5% polyacrylamide 2D gel. Black circles represent the 107 proteins identified by MALDI-TOF/TOF MS with significant differences in abundance in NET groups and black squares the 2 proteins with significant differences in abundance in HYP groups (P < 0.05) grouped based on the expression patterns of the protein spots. The PCA ( Fig.4-b) showed two main clusters belonging to the control and NET4 samples, while 2 biological samples belonging to the NET2 group clustered together with the control samples and 1 with the NET4 samples. The 107 differential protein spots were centralized into two principal components (PC), PC1 and PC2, which represented the maximum variation (65.6%) and the next highest variation (5.5%), respectively. The HCA ( Fig. 4-c) likewise revealed two main groups regarding the biological replicates, as observed by the top dendrogram. The protein spots were also grouped into two main clusters, one displaying a pattern of higher relative abundance and the other of lower relative abundance in stressed fish, when compared to the control. As described above for the PCA, higher variability in NET2 was also shown in the HCA.
For the network and GO enrichment analyses the subset of 20 single protein identifications mentioned above was blasted against Danio rerio in the UniprotKB database. A PPI network ( Fig.5-a) was generated on the STRING web tool revealing 61 edges among 18 nodes/ proteins (2 proteins had no interaction with the main network), with a clustering coefficient of 0.677 and a very significant enrichment value (P < 1.0e − 16 ). The analysis was performed on Cytoscape and specific topological parameters were selected to demonstrate the importance and distribution of the nodes in the network: a darker colour intensity of the nodes indicates higher degree, while the size was estimated using the variation in protein abundance (fold-change). For every single entry, one protein spot was chosen as the most representative of each protein (Table 1), based mainly on the protein score and experimental molecular weight and pI close to the theoretical ones. From these 18 spots, 11 were down-and 7 were up-regulated, however, these differences in abundance were mostly significant (log-fold change > 1.0 or < − 1.0, q-value < 0.05) for the NET4 treatment (only 2 were exclusively significant for the NET2 treatment and 2 were significant for both treatments). Thus, the fold-change of these 18 spots between NET4 and CTRL groups was used to estimate the size of the nodes on the PPI network, which ranged from − 4.04 to + 2.78. SERPINC1 (antithrombin-III), TFA (transferrin) and FGA (fibrinogen alpha-chain) occupied the most central positions in the network having the highest number of interactions, while APOA1 (apolipoprotein A-I) showed the highest number of experimentally overrepresented (hypergeometric test, FDR < 0.05) GO Biological Process (BP) terms, mostly linked to the immune system and response to stimulus. No annotations were retrieved for alpha-1-antitrypsin, leucine-rich alpha-2glycoprotein-like, apolipoprotein A-I, apolipoprotein B-100-like, haptoglobin, pentraxin and hyaluronic acidbinding protein 2. In the horizontal bar plot ( Fig. 5-b), only the 9 most significant terms are represented. GO Molecular function enrichment analysis accounted for 8 terms with 5 main proteins (alpha-1-antitrypsin, antithrombin-III, inter-alpha-trypsin-inhibitor, kininogen and alpha-2-macroglobulin) while GO Cellular component revealed 4 enriched terms with 2 main proteins (fibrinogen alpha-chain and alpha-2-macroglobulin). A complete list of all GO terms is described on the additional file 3.

Discussion
In this study, the stress response of farmed Gilthead seabream adults to chronic stress conditions was primarily assessed by observing both changes in the concentration of routine plasma stress indicators, namely cortisol, glucose and lactate, and post-mortem biochemical parameters, explicitly pH and rigor mortis. To evaluate the existence of possible unbiased and reliable markers of chronic stress, proteomics was used to verify the potential of fish protein-based adaptations.
Cortisol is the most commonly used physiological indicator of the primary response to stress [21]. However, it has been shown that this corticosteroid is not a reliable biomarker of long-term stress exposure [27,[35][36][37]. In this study, Gilthead seabream that endured high stocking densities over 54 days showed a possible reconfiguration of the cortisol response. This is supported by the observed downward trend of this metabolite, as compared to unstressed fish. This is suggestive of either a non-activated or an altered responsiveness of the HPI axis, which sometimes leads to the hyporeactivity of the corticosteroid response [38]. The same outcome was observed in juvenile Gilthead seabream confined for 14 days at 26 kg/m 3 [39] and in meagre, cultured at different stocking densities for 40 days [40]. In the NET trial, however, apart from the wide dispersion of observations, plasma cortisol levels were significantly higher in handled fish. This result suggests that the fish were not able to appropriately adapt to the handling stressor. Its persistence, unpredictability and severity, could have prevented the possibility of habituation. Regarding the HYP trial, no effect of the 48 h of hypoxia was reflected in these fish. This suggests an acclimation to the low oxygen environment by a possible adjustment of the oxygen requirement (e.g. reduction of high energy behaviours).
Overall, the aforementioned observations suggest that the cortisol response and the capacity of adaptation are modulated by the nature, duration and intensity of the stressor. However, other factors like species, age, sex and individual coping mechanisms seem to be ubiquitous and impact their adaptive processes [24,38,41]. This process of stress habituation was already suggested and demonstrated in other studies [37,42], but this mechanism is not yet well-understood. High individual Fig.5 A -Protein-protein interaction network generated with 18 differential proteins identified in the plasma of fish from NET trial. Nodes represent proteins and edges the functional associations between them. STRING annotations are described in Table 1. Red arrows represent upregulated proteins in both treatments; blue arrows represent down-regulated proteins in both treatments. D -GO Enrichment analysis of the 18 proteins showing significantly differential abundance between control and NET treatments (hypergeometric test, FDR < 0.05) variability was also observed in each trial, most likely due to individual differences in the stress response related to intrinsic factors of the animal (e.g. coping styles, cognitive perception) [39,43,44]. Additionally, values registered for control fish, in every trial, are higher than the reference values reported in the literature for this species [45]. These discrepancies may have several causes, which is why cortisol should be used with caution when evaluating the magnitude of the stress response. Moreover, the difficulty of measuring resting cortisol levels is also acknowledged to be one of the causes of these discrepancies. The lack of proper planning when sampling cortisol, or the manipulation needed to net and anaesthetize the fish, can result in high "control" cortisol levels that do not correspond to the "genuine" basal levels i.e., the non-manipulated fish levels. Also, it is well-established that following the perception of an acute stressor, the levels of circulating stress markers increase within the first minutes or hours of stress response, returning to basal levels as time elapses, usually within 24 h [17,46,47].
Secondary physiological responses are characterized by an increase in glucose and lactate levels in blood plasma in order to satisfy the increased energy expenditure. Changes in glucose usually follow similar trends as cortisol after the stressor [14]. This is corroborated in this study by the levels of plasma glucose registered in all trials (Fig. 1). Glucose levels, besides following the same trend as cortisol levels, are, in general, below the basal values for this species [45]. This could be related to the fish's inability to maintain the same levels of glucose in the blood due to the high demand for glucose mobilization to other tissues. The decrease of plasma glucose levels in OC is in agreement with the decrease in the cortisol levels, supporting the hypothesis of habituation or exhaustion of the endocrine system [27]. The significant increases in the plasma glucose levels of stressed fish from NET and HYP trials are consistent with previous studies. These showed that glucose rises during air exposure or low oxygen levels, due to stimulation of muscle glycogenolysis and hepatic gluconeogenesis, where glucose is synthesized to maintain the Table 1 String annotations and fold-changes of the proteins in the PPI network. Bold lettering in the "FC" column indicates significant fold-changes (> 1.0 and < − 1.0). List is given in ascending order of spot number  (Fig.3), attributed by the SameSpots software b Accession number -NCBI accession number c Protein IDprotein identification by MALDI-TOF/TOF MS d FC -Log2(fold-change) -significant changes in protein abundance (treated/control). Bold lettering indicates significant fold-changes (> 1.0 and < − 1.0) energetic substrates' demand [48]. Similar to cortisol, glucose and lactate circulating levels also return to basal levels within hours post-stressor, which further makes of these metabolites unreliable markers in case of prolonged stressors [49,50]. Additionally, studies also demonstrate that glucose variations in the blood are not only hormonal-induced due to stressful practices. Factors like variations in the water temperature and pH, anaesthesia, diet composition or fasting can also affect plasma glucose levels [51,52]. When insufficient oxygen is available to maintain the aerobic ATP production, fish resort to anaerobic metabolism to meet cellular requirements. This shift consequently leads to lactate accumulation in the muscle [19,53]. In this study, changes in the circulating lactate levels do not follow the same trends of cortisol and glucose variations. Statistically significant differences in the lactate levels were only observed in the OC trial. In this case, if the cortisol response is indeed lower due to HPIaxis acclimation, as suggested before, the lactate recycling rate in the hepatic glycogenolysis is reduced, explaining the significant plasma lactate increase in stressed fish. Additionally, previous studies show that during hypoxia or intense swimming activity, fish produce lactate in the muscle at a higher rate than it can be processed by other tissues [53].
Post-mortem muscle pH and rigor mortis have been used as tissue indicators of ante-mortem stress in numerous fish species [54][55][56]. After the fish death, both blood circulation and oxygen supply cease. The major source of ATP to the muscle is thus lost, since glycogen can no longer be oxidized. However, for a limited time after death, ATP in the muscle is maintained at a definite level by creatine kinase. Consequently, the depletion of ATP reserves stimulates the breakdown of glycogen by anaerobic glycolysis in the muscle, in order to maintain the energy expenditure. This process results in the accumulation of lactic acid, generating H + ions and consequently lowering muscle pH [57]. Glycolysis continues until all glycogen is consumed or the glycolytic enzymatic system is made inactive by the low pH. Hence, the magnitude and rate of this pH fall depend on the fish's energy reserves prior to death. These energy reserves can be influenced by the intensity and duration of the stress while fish is alive. To our knowledge, no studies were performed with Gilthead seabream regarding the effects of long-term chronic stressors on the evolution of post-mortem biochemical processes in the muscle. Results from this study (Fig. 2) followed the same pH trends as previous studies on gilthead seabream [58,59], however, comparing this with the existent studies on pre-slaughter stress [55,60,61], muscle pH values immediately after death are lower than the ones found in this study, suggesting that stress at slaughter was low in our fish. Poli et al. 2005 state that in cases of exposure to a chronic stressor for a long time before death, the lactic acid produced can be gradually cleared from the muscle, but simultaneously the energy sources, like glycogen, will likewise gradually become exhausted. Hence, when the fish is killed, muscle pH does not suffer a dramatic fall due to an early end of post-mortem anaerobic glycolysis caused by energy source scarcity. This might explain the significant differences found in the HYP trial, where the highest pH values were observed in the highly stressed fish (HYP15), suggesting that these fish had lower energy reserves. Nevertheless, pH values registered after the 24 HAD, in every treatment, are in agreement with the reported by previous studies in this species at the same sampling times [58,62].
Rigor mortis is inextricably correlated with muscle ATP and the pH decline. The onset of rigor mortis occurs with ATP depletion. When ATP reaches low levels, actin and myosin in the muscle bind together, forming the actomyosin complex and causing stiffness of the fish body [63]. A strong relationship between low muscle pH immediately after death, and a rapid onset of the rigor state was demonstrated in a range of fish species [57,60]. In this study, the evolution of rigor mortis (Fig. 2) was similar between treatments and significant differences were only found in the NET and HYP trials at 8, and at 8 and 24 HAD, respectively. A delayed onset was observed, starting between 2 and 6 HAD in every trial and reaching the maximum rigor index between 24 and 48 HAD. This delay is in agreement with the high muscle pH registered immediately after death, supporting the hypothesis of low energetic reserves in our fish at the time of death. Measuring glycogen and ATP content in the fish muscle and liver would be a complementary assessment to infer about the energetic reserves and corroborate our hypothesis.
Plasma proteins were evaluated in this study since blood plasma is a very informative biological fluid, acting as a mirror of the physiological condition of the organism. Stress and stress-related hormones are recognized as modulators of the fish immune system [64], however, responses depend on the intensity and duration of the stressor. The innate immune system is a fundamental defence mechanism in fish [65]. The acute phase response is part of this system and it is mainly regulated by cytokines and glucocorticoids [66]. This response is characterized by the release of acute-phase proteins (APP), by the hepatocytes, into circulation [67]. APP can be classified as "positive" or "negative" depending on whether their plasma concentration increases or decreases during activation of this response [68]. The response profile of our fish demonstrated the same tendency of protein changes.
In this study, the pattern of protein changes observed in the plasma indicate that the fish's immune system was affected mainly by net handling and hypoxia stressors. Nevertheless, net handling was shown to be the most impacting. The levels of 20 different plasma proteins (distributed by 56 significantly differential spots), all related with immunological processes, were shown to be modulated by repetitive net handling, compared to two proteins modulated by hypoxia. As mentioned previously, the same proteins were often detected from different spots on the 2D gels. Such a phenomenon may be due to existent isoforms or caused by adaptive changes of the proteome in an attempt to maintain cellular homeostasis. This adaptation may involve changes at the level of protein degradation, localization, function and activityall of which can be modulated by posttranslational modifications (PTMs) [69]. PTMs can regulate fundamental biochemical processes and be more energetically efficient than altering protein abundance, constituting potential interesting signatures of stress. Studies on PTMs in fish are still scarce.
The changes detected in protein abundance (listed in additional file 2), along with the PPI network and GO enrichment analyses (Fig. 5) performed, confirmed the involvement of several components of the innate immune system in the physiological adaptation to these stressors. Proteins considered to be "positive" APP were likewise shown to be increased in abundance in the plasma of fish stressed by net handling (fibrinogen alpha-chain, complement component C3, haptoglobin, complement factor B, warm-temperature acclimation 65 kDa protein, alpha-1-antitrypsin), while proteins considered as "negative" were decreased (transferrin, interalpha-trypsin inhibitor, apolipoprotein A-I) [70]. A diverse number of proteins involved in the APR was also found previously to be modulated in chronically stressed gilthead seabream [71].
Apolipoprotein A-I (Apo-AI) was modulated only by net handling stress. Seventeen proteoforms were identified in the plasma proteome map, being mostly decreased in abundance when comparing with the unstressed fish. Apo-AI is the main protein constituent of the high-density lipoprotein (HDL), playing a role in lipid metabolism and participating in the reverse transport of cholesterol from tissues to the liver [72,73]. Apo-AI was also found to be decreased in abundance in crowded Atlantic salmon [74]. In cod (Gadus morhua) it acted as a negative regulator of the complement system [75]. Other two apolipoproteins were also found to be down-regulated in the plasma of fish from NET2 and NET4 groups (Apolipoprotein Eb and apolipoprotein B-100).
The complement system is an essential part of the innate immune system which can be activated through three pathways: the classical, alternative and lectin pathways [76]. Fish display a plethora of complement components, mainly complement component C3 (C3), which may present around five proteoforms in a single species [77]. C3 is one of the most abundant proteins in the plasma and plays a central role in the innate immune system, supporting the activation of all three pathways [76]. In this study, C3, identified in 5 proteoforms, and complement factor B (Bf), identified in 4, were found to be increased in abundance by net handling. Contrarily, C3 was down-regulated in fish exposed to low oxygen levels. Bf also plays a role in complement activation by acting as the catalytic subunit of C3 convertase, an enzyme responsible for the proteolytic cleavage of C3, in the classical and alternative pathways [76].
Several metal-binding proteins, existent in the plasma of vertebrates, can chelate iron, zinc and copper, which are essential elements for the virulence of bacteria [78]. Alpha-2-macroglobulin (A2M) is a multifunctional protein [79] found to be down-regulated in the plasma of fish submitted to handling stress. It is mostly known to act as a broad range serine proteinase inhibitor and to bind metal ions [78]. Contrarily, haptoglobin, which is also responsible for the sequestration of iron by binding to hemoglobin, was found to be increased in the plasma of handled fish. Similarly, warm-temperature acclimation-related 65 kDa protein (Wap65), which is involved in the scavenging of free heme [80], was increased in abundance by net handling and hypoxia stressors. Wap65 in fish is the homologue of mammalian hemopexin [81] and in most teleosts presents two proteoforms [82]. In this study, two spots were also matched to this protein suggesting the presence of these two proteoforms. Transferrin (Tf) decreased in abundance in the plasma of fish stressed by net handling. Tf is a plasma protein also capable of binding iron and an important constituent of the iron homeostasis [33].
In fish, antiproteases are important participants of the non-specific humoral immune defence mechanism [70]. A2M is an important factor in this mechanism. Alpha-1-antitrypsin is a serine protease inhibitor, upregulated in net-handled fish, which is responsible for negatively regulating blood clotting molecules to prevent thrombosis [83]. Inter-alpha-trypsin inhibitor H3 is also a serine protease inhibitor, which was found to be down-regulated in the plasma of fish from NET groups. The same pattern of protein changes was verified for fetuin-B, a cysteine proteinase inhibitor recently described in teleosts [84]. Finally, fibrinogen alpha-chain, a beta-globulin involved in blood clotting, an integral part of innate immunity [83], was found to be up-regulated in the plasma of fish belonging to NET groups.

Conclusions
In summary, the overall results suggest that physiological changes were higher in fish exposed to repeated handling, while mild and permanent stressors may allow the fish to refine their physiological processes and adapt to certain challenges. The variability in the response levels of cortisol, glucose and lactate, in fish from the same groups, alongside the possible adaptation suggested by the results, demonstrate that these indicators may not be the most robust in case of chronic stress monitoring. On the other hand, plasma proteomics allowed the detection of a cohesive network of protein changes associated with essential immunological pathways in stressed fish. These proteins will be useful in understanding the biological processes behind protein-based stress adaptation in fish and may, therefore, represent the first screening for potential biomarker candidates of chronic stress in gilthead seabream. This work is the first step for a more scientific and reliable assessment of fish welfare. A multidisciplinary approach, and the study of the stress response from the molecular to the behavioural level might just be the holistic approach needed to achieve such a goal.

Animals
Gilthead seabream (Sparus aurata) were obtained from a commercial fish farm (Maresa, Mariscos de Estero S.A., Huelva, Spain) and kept under quarantine conditions for a 2-week period at the Ramalhete Research Station (CCMAR, University of Algarve, Faro, Portugal). The fish were then individually weighed and distributed among conical fiberglass tanks (500 L), according to the density requirements of each trial. The tanks were supplied with natural flow-through seawater from Ria Formosa, and kept under natural temperature (13.4 ± 2.2°C) and photoperiod, salinity at 34.7 ± 0.8 ‰, and artificial aeration (dissolved oxygen above 5 mg. L − 1). Fish were fed by hand once a day, with a diet manufactured by AquaSoja Portugal, following the species' nutritional requirements.

Experimental design
The study was performed in three separate trials: [1] Overcrowding (OC), [2] Net Handling (NET) and [3] Hypoxia (HYP), due to logistic issues. Each trial followed a 2-week acclimation period and the initial rearing density was established at 10 kg/m 3 (except in the experimental groups of high stocking densities). In the OC trial, during the 54 days of experiment, fish (initial body weight (IBW) = 372.33 ± 6.55 g) were stressed using different high stocking densities, by increasing the number of fish in the tanks. Three different experimental groups were tested in triplicate: Control -10 kg/m 3 (OC CTRL ), medium density -30 kg/m 3 (OC 30 ), high density -45 kg/m 3 (OC 45 ). The NET trial lasted for 45 days and the fish (IBW = 375.69 ± 11.88 g) were stressed by 1-min air exposure, using nets designed to fit inside the tanks and to be lifted to perform the stressful event. The experimental groups were established, in triplicate, as follows: Controlundisturbed fish (the net was also placed in the tanks but not lifted) -(NET CTRL ), fish air-exposed twice a week (NET 2x ) and fish air-exposed four-times a week (NET 4x ). In the HYP trial, fish (IBW = 397.99 ± 16.56 g) were subjected to low levels of saturated oxygen, by injection of nitrogen in the water, for 48 h, according to the following experimental groups (in triplicate): Control -100% saturated oxygen -(HYP CTRL ), 30% saturated oxygen (HYP 30 ) and 15% saturated oxygen (HYP 15 ). Different trial times are due to differences in the nature and severity of the stressor, to which rearing protocols had to be adjusted accordingly.

Sampling procedure
Prior to the sampling day, fish were starved for 48 h to clean the digestive tract. Nine random fish per tank were lethally anaesthetized with tricaine methanesulfonate (MS-222; Sigma Aldrich, St. Louis, Missouri, USA) for the following sampling procedures: 3 fish for rigor mortis index assessment, 3 fish for muscle pH measurement and 6 fish for blood collection. Blood samples of approximately 2 ml were collected from the caudal vein with a heparinized syringe and immediately centrifuged at 2000 g for 20 min. Plasma samples were immediately frozen at − 80°C until posterior analyses. Fish for the measurement of post-mortem biochemical changes (pH and rigor mortis) were stored in polystyrene boxes with ice during the sampling period (72 h). All fish were weighed and measured.

Plasma stress indicators' measurement
Plasma cortisol levels were quantified using a commercial Cortisol ELISA kit RE52061 (IBL International, Hamburg, Germany), following the manufacturer's instructions. Measurements were registered at 450 and 620 nm along with a prepared standard curve on a microplate reader Biotek Synergy 4 Hybrid Technology™ (Biotek Instruments Inc., Winooski, USA). Plasma glucose and lactate levels were assessed through commercial colorimetric kits (Spinreact, Girona, Spain), following the manufacturer's instructions.

Biochemical and quality characterization of fish muscle
Muscle pH measurements were performed (n = 3 per tank), using a waterproof pH spear for food testing (Oakton® Instruments, Nijkerk, Netherlands), in the dorsal muscle, at 0, 1, 2, 4, 6, 8, 24, 48 and 72 h after death (HAD), approximately 1-2 cm apart. At the same post-mortem periods, rigor mortis was assessed (n = 3 per tank) by the rigor index (RI), as previously described [85], using the formula: L 0 (cm) refers to the vertical distance between the base of the caudal fin and the table surface (where the anterior half of the fish is placed), measured immediately after death, whereas L t (cm) corresponds to the same distance, however at selected time intervals. Fish were carefully handled during the measurements to avoid any interference with the rigor onset.

Plasma proteomics analysis Protein labelling
Plasma samples were diluted 80x in DIGE buffer (7 M urea, 2 M thiourea, 4% CHAPS, 30 mM Tris pH 8.5) and the protein content measured with Bradford assay using the BioRad Quick Start Bradford Dye Reagent 1X (Bio-Rad Laboratories, Hercules, California, USA) and bovine serum albumin (BSA) as standard, BioRad Bovine Serum Albumin Standard Set (Bio-Rad Laboratories, Hercules, California, USA). Samples' pH was checked with a pHindicator paper, Sigma-P4536 (Sigma Aldrich, St. Louis, Missouri, USA) and adjusted to 8.5 using 0.1 M NaOH. DIGE minimal labelling of 50 μg of protein was carried out using the CyDye™ DIGE fluor minimal labelling kit 5 nmol (GE Healthcare, Little Chalfont, UK), with 400 pmol fluorescent amine reactive cyanine dyes freshly dissolved in anhydrous dimethylformamide (DMF), following the manufacturer's instructions. Labelling was achieved on ice for 30 min, in the dark, and the reaction quenched with 1 mM of lysine for 10 min. For each trial, six samples per experimental condition were labelled with Cy3 and six with Cy5 to reduce the impact of label difference, while an internal standard consisting of a pool of all samples, with equal amounts, was labelled with Cy2. Samples were randomly sorted to avoid labelling bias.

Image acquisition and analysis
CyDye-labeled gels were scanned on a Typhoon™ laser scanner 9400 (GE Healthcare, Little Chalfont, UK) at 100 μm resolution, with the appropriate laser filters for the excitation and emission wavelengths of each dye (i.e., Cy2-488/520 nm; Cy3-532/580 nm; and Cy5-633/670 nm), according to the manufacturer's recommendations. The voltages of the Photo Multiplier Tube (PMT) were adjusted to obtain a maximum image quality with minimal signal saturation and clipping. Gel images were checked for saturation during the acquisition process using the ImageQuant TL software (GE Healthcare, Little Chalfont, UK). The final images were analysed with SameSpots software (Totallab, Newcastle, UK), including background subtraction (average normalized volume ≤ 100,000 and a spot area ≤ 500), filtering, spot detection, spot matching, normalization and statistical analysis. Spot volume ratios that showed a statistically significant difference (abundance variation of at least 1.0-fold, P < 0.05 -one-way ANOVA on log2transformed normalized spot volumes) were processed for further analysis. Protein spots with statistically different intensities were manually excised from preparative gels and identified by matrix-assisted laser desorption/ ionization time-of-flight/time-of-flight mass spectrometry (MALDI-TOF/TOF MS).

Protein identification by MALDI-TOF/TOF MS
Spots from SYPRO® Ruby-stained (InvitrogenTM, Carlsbad, CA, USA) gilthead seabream plasma 2D gels were picked and subjected to in-gel tryptic digestion, similar as reported before [86]. In this study, gel plugs were washed twice with 50 mM ammonium bicarbonate solution in 50% v/v methanol (MeOH) for 20 min and dehydrated twice for 20 min in 75% acetonitrile (ACN). Proteins were then digested with 8 μL of a solution containing 5 ng/μL trypsin (trypsin Gold, Promega) in 20 mM ammonium bicarbonate for 6 h at 37°C. A 0.1% trifluoroacetic acid (TFA) solution in 50% ACN and a solution of 7 mg/mL α-cyano-4-hydroxycinnamic acid (CHCA) in 50% ACN/0.1% TFA were used for peptide extraction and spotting respectively. MALDI TOF/TOF analysis was performed with a TOF/TOF™ 5800 (AB SCIEX, Redwood City, CA, USA) mass spectrometer in MS and MS/MS mode. For each spot, the 10 most intense peaks of the MS spectrum were selected for MS/ MS acquisition. Database interrogation was carried out over with ProteinPilot v4.5 (AB Sciex) on an in-house Mascot server version 2.6.1 (Matrix Science Ltd., London, UK).
Mass lists were searched against NCBInr database restricted to the taxonomy "other Actinopterygii" (tax ID 7898 excluding 31,033 and 7955) with the following parameters: maximum 2 missed cleavages by trypsin, peptide mass tolerance ±100 ppm, fragment mass tolerance set to 0.5 Da, carbamidomethylation of cysteine selected as fixed modification and tryptophan dioxidation, histidine, tryptophan and methionine oxidation, and tryptophan to kynurenine as variable modifications. Protein hits not satisfying a significance threshold (P < 0.05 and a total ion score > 60) were further searched against vertebrate EST (expressed sequence tags) database also restricted to the taxonomy "other Actinopterygii".

Protein-protein interaction (PPI) network and gene ontology (GO) enrichment analyses
The theoretical molecular masses and isoelectric points (pI) of the MS identified proteins were calculated using the amino-acid sequences (in one-letter code) on the ProtParam Tool (http://us.expasy.org/tools/protparam. html). A significance cutoff was applied for the identified proteins at log-fold change ±1.0. Following, the identified proteins were blasted against Danio rerio, on the UniprotKB database, using the FASTA protein sequences as queries. The orthologues were mapped using STRING web tool v11.0 (https://string-db.org/) to screen for protein-protein interactions (PPI). Gene ontology (GO) enrichment analysis and network visualization and analysis were performed on Cytoscape v3.7.1 (http:// www.cytoscape.org/) with the BiNGO plug-in. Important hub proteins were screened by counting the degree of connectivity of each node in the network. Overrepresented GO terms were identified, using B. rerio as reference, by selecting the hypergeometric test with a significance threshold of 0.05 after Benjamini & Hochberg FDR correction.

Statistical analyses
All univariate and multivariate statistical analyses were performed using R v3.5.3 for MacOSX (https://www.rproject.org). Statistical analyses of the plasma parameters and the post-mortem muscle biochemical changes were performed using plasma cortisol, glucose and lactate levels, muscle pH and rigor index as dependent variables, and the stress treatment as factor. Statistical differences between treatments were analysed independently for each trial (OC, NET and HYP). For rigor index and muscle pH, data were processed separately for each sampling time. Differences in plasma and muscle parameters between treatments were assessed by a one-way analysis of variance (one-way ANOVA) on log10-transformed data, except for rigor mortis data, which was transformed by arcsine square root. Multiple comparisons were carried out by the post-hoc Tukey HSD test. When transformed data failed the Shapiro-Wilk normality test, the non-parametric Kruskal-Wallis on ranks was used, followed by Dunn's test. When transformed data did not verified homoscedasticity assumption by Levene's test, statistical significance was analysed by Welch's ANOVA, followed by Games-Howell. A significance level of α = 0.05 was used in all tests performed. Experimental data is expressed as mean ± standard deviation (SD). Principal component analysis (PCA) and hierarchical clustering analysis of the identified proteins were performed on the log2-transformed normalized spot volumes obtained from SameSpots software, with autoscaling. Heatmap was generated by comparing Z-scores of normalized spot volumes and hierarchical clustering of samples and protein spots was performed using the Euclidean distance and the maximum cluster agglomeration method as distance metrics.
Additional file 2. Protein spots, with statistically different relative abundances (P < 0.05), identified in gilthead seabream (Sparus aurata) blood plasma proteome from NET and HYP trials, by MALDI-TOF/TOF MS after separation by 2D-DIGE. List is given in ascending order of spot number.
Additional file 3. List of the overrepresented terms in the GO Enrichment analysis of the 18 proteins showing significantly differential abundance between control and NET treatments (hypergeometric test, FDR < 0.05).