Untargeted metabolomics approach for unraveling robust biomarkers of nutritional status in fasted gilthead sea bream (Sparus aurata)

Reagents and chemicals

HPLC-grade water was obtained from a Mili-Q water purification system (Millipore Ltd., Bedford, MA, USA). HPLC-grade methanol (MeOH), HPLC-supergradient acetonitrile (ACN), sodium hydroxide (>99%) and reagent grade ammonium acetate (NH₄Ac) were obtained from Scharlab (Barcelona, Spain). Leucine-enkephalin (mass-axis calibration), formic acid (mobile phase modifier) and analytical-grade standards methionine sulfoxide and trimethylamine N-oxide were purchased from Sigma-Aldrich (Saint Louis, MO, USA).

Animal care and sampling

Two-year-old gilthead sea bream of Atlantic origin (average initial weight: 380 g) were reared from early life stages in the indoor experimental facilities of the Institute of Aquaculture Torre de la Sal (IATS), following natural light and temperature conditions at our latitude (40°5′N, 0°10′E). The oxygen content of water was always higher than 85% saturation, unionized ammonia remained below toxic levels (<0.02 mg/l), and rearing density was maintained lower than 15 kg/m³.

At mid-summer (July 2014), 30 fish were randomly allocated in two tanks (500 L). One group continued to be fed with a standard commercial diet (Biomar, EFICO Forte 824) to visual satiety one time per day, whereas the other group remained unfed for a 10-day period. At the end of this period, 10 fish from fasted and fed groups (following overnight fasting) were randomly sampled and anaesthetized with 100 mg/L of aminobenzoic acid ethyl ester (MS-222, Sigma-Aldrich) for blood and tissue sampling. Blood was taken from caudal vessels with vacutainer tubes with a clot activator. Liver and visceral adipose tissue were extracted and weighed. Blood samples were allowed to clot for 30 min at room temperature, and then centrifuged at 1,300 g for 10 min. The obtained samples were stored at −20 °C until analysis.

All procedures were approved by the IATS Ethics and Animal Welfare Committee according to national (Royal Decree RD53/2013) and EU legislation (2010/63/EU) on the handling of animals for experiments.

Sample processing

Serum samples were centrifuged at 12,500 g for 10 min. Supernatant (400 µL) was diluted with 1.2 mL of ACN followed by centrifugation (12,500 g for 10 min). Then, 750 µL of supernatant were stored for hydrophilic interaction liquid chromatography (HILIC), and another 750 µL aliquot was evaporated to dryness by MiVac Duo Concentrator (40 °C, 60 min) and dissolved with MeOH (75 µL) and Mili-Q Water (675 µL) for reversed phase (RP) analysis (details in Fig. S1). Quality control (QC) samples were prepared by pooling 50 µL of each sample extract. All samples were stored at −20 °C until injection.

UHPLC-HRMS

A Waters Acquity UPLC system (Waters, Milford, MA, USA) was coupled to a hybrid quadrupole-TOF mass spectrometer (Xevo G2 QTOF, Waters, Manchester, UK), using a Z-spray-ESI interface operating in positive and negative ionization mode. The UHPLC separation was performed using Acquity UPLC^® BEH C18 1.7 µm particle size analytical column 100 × 2.1 mm (Waters) at 300 µL/min flow rate for RP analysis. An Acquity UPLC^® HILIC 1.7 µm particle size analytical column 100 × 2.1 mm (Waters) at 300 µL/min flow rate was used for hydrophilic interaction phase separations.

Each serum sample was injected four times, depending on the procedure (RP and HILIC) and the ionization mode selected (ESI+ and ESI−). The RP separation was performed using H₂O with 0.01% formic acid (HCOOH) as weak mobile phase (A) and MeOH with 0.01% HCOOH as strong mobile phase (B). The percentage of B was changed from 10% at 0 min, to 90% at 14 min, 90% at 16 min and 10% at 16.01 min, with a total run time of 18 min for both ESI+ and ESI−. For HILIC separation, the weak mobile phase was a mix of ACN:H₂O (95:5, v/v) with 0.01% HCOOH and 10 mM NH₄Ac (A), and the strong mobile phase was H₂O with 0.01% HCOOH and 10 mM NH₄Ac (B). The B percentage was changed as follows: 0 min, 2%; 1.5 min, 2%; 2.5 min, 15%; 6 min, 50%; 7.5 min, 75%; and finally at 7.51 min, 2%, with a total run time of 10 min, for both ESI+ and ESI−. Sample injection volume was 10 µL in all cases. Nitrogen was used as both the desolvation gas and the nebulizing gas. A capillary voltage of 0.7 kV and 1.5 kV for positive and negative ion modes, respectively, and cone voltage of 25 V were used. MS data were acquired over a m∕z range of 50–1,200. TOF-MS resolution was approximately 20,000 at full width half maximum at m∕z 556.2771. Collision gas was argon 99.995% (Praxair, Valencia, Spain). The desolvation gas flow was set at 1,000 L/h, and the cone gas was set at 80 L/h. Desolvation gas temperature was set to 600 °C, source temperature to 130°C and column temperature to 40°C.

For MS^E experiments, two acquisition functions with different collision energies were created. The low energy (LE) function, with a fixed collision energy of 4 eV, and the high energy (HE) function, with a collision energy ramp ranging from 15 to 40 eV in order to obtain the (de)protonated ion from LE function and a wide range of fragment ions from the HE function. Both LE and HE functions used a scan time of 0.3 s with an inter-scan delay of 0.05 s. MS/MS experiments were carried out in the same conditions with different collision energies depending on the fragmentation observed for each compound. Calibrations were conducted from m∕z 50 to 1,200 with a 1:1 mixture of 0.05 M NaOH:5% HCOOH diluted (1:25) with H₂O:ACN (20:80), at a flow rate of 10  µL/min. For automated accurate mass measurement, a leucine-enkephalin solution (0.5 µg/mL) in ACN:H₂O (50:50) at 0.1% HCOOH was pumped at 30  µL/min through the lock-spray needle and measured every 30 s, with a scan time of 0.3 s. The (de)protonated molecule of leucine-enkephalin, at m∕z 556.2771 in positive mode and m∕z 554.2615 in negative mode was used for recalibrating the mass axis during the injection and to ensure a robust accurate mass along time.

Data processing

The workflow of data processing is shown in Fig. 1. LC-MS spectral data were converted from proprietary (.raw, Waters Corp.) to generic (.cdf, NetCDF) format using Databridge application (within MassLynx v 4.1; Waters Corporation) and processed using XCMS R package (https://xcmsonline.scripps.edu/) (Smith et al., 2006). Centwave feature detection algorithm was employed for peak picking (peak width from 5 to 20 s, S/N ratio higher than 10 and mass tolerance of 15 ppm) followed by retention time alignment for the detected features. Peak area normalization (mean centering) was applied to each data set in order to minimize instrumental drifts with a final log2 transformation to the area to standardize the range of independent feature variance followed by pareto scaling. ANOVA analysis followed by Benjamini–Hochberg multiple testing correction was applied to the normalized peak areas of all metabolites to assess differences between fed and control groups.

Multivariate analysis of processed metabolomics data was performed by means of the EZ-Info software (Umetrics, Sweden). First, Principal Component Analysis (PCA) was employed to ensure the absence of outliers and the correct classification of QCs after normalization. Partial Least Squares—Discriminant analysis (PLS-DA) was then applied to maximize the separation of fed and fasted individuals (Fonville et al., 2010). Orthogonal PLS-DA (OPLS-DA) was also carried out (Wiklund et al., 2008) with a high threshold (P[corr] > 0.95) for highlighting the most robust biomarkers.

For elucidation, the MS/MS spectra of the most significant metabolites were compared with reference spectra databases (METLIN, http://metlin.scripps.edu; Human Metabolome DataBase, http://www.hmbd.ca; MassBank, http://www.massbank.eu). For unassigned metabolites, in silico fragmentation software (MetFrag, http://msbi.ipb-halle.de/MetFrag) was employed, with subsequent searches through Chemspider (http://www.chemspider.com) and PubChem (https://pubchem.ncbi.nlm.nih.gov) chemical databases. Injection of standards of methionine sulfoxide and trimethylamine N-oxide served to validate the elucidation workflow.

A retrospective analysis of data previously acquired in MS^E mode served for the refined search of additional relevant metabolites. It consisted in the search of the m∕z ratio (parent ions) of the metabolites of interest in the LE function as well as product ions obtained from MS/MS spectrum online databases (METLIN and Human Metabolome DataBase) in the HE function. Integrated areas of each candidate (parent ion) were compared in samples from fed and fasted groups.

Biometric data

At the end of the experimental period, body weight of fed fish was 15% higher than in fasted fish. This fasting protocol reduced the body fat depots, decreasing significantly (P < 0.001) the hepatosomatic index (100 × liver weight/body weight) from 1.3 to 0.9. The similar trend was found for mesenteric fat, although the decrease of mesenteric fat index (100 × mesenteric fat weight/body weight) from 1.9 to 1.6 was not statistically significant (Table 1). The magnitude of these changes was on the range of expected values for one- and two-year-old fish under similar experimental conditions (Benedito-Palos, Ballester-Lozano & Pérez-Sánchez, 2014; Bermejo-Nogales, Calduch-Giner & Pérez-Sánchez, 2015).

Untargeted metabolomics fingerprinting

Despite of the great potential of GC for chromatographic separation, the nature of serum samples, with medium-high polar compounds in a water-based fluid, pointed out to LC as a more convenient separation technique. UHPLC with sub-2 µm particle size was applied due to its high reproducibility and high separation performance in short-run time analyses. The use of different chromatographic techniques is a key issue to achieve a maximum of detected features when dealing with complex matrices like blood. In our case, serum samples were analyzed with two ionization modes and two different chromatographic columns: RP for a better separation of non-polar compounds, and HILIC to best separate the most polar compounds. In the RP analysis, 6,961 and 3,047 features were detected in both positive and negative ionization modes, while 4,820 and 1,015 features were labeled by XCMS using HILIC separation. This high total number of detected features (m∕z values) highlights the huge detection power and sensitivity of HRMS and makes feasible a wide-view of sample composition to discriminate the most robust markers of nutritional conditions. Many features were only observed under a single ionization mode and chromatography type, reinforcing the importance of employing different chromatographic columns. As an example, a single peak was detected by HILIC for the significant feature elucidated as LysoPC(20:5) while RP chromatography was able to separate ω-3 and ω-6 isomers (Fig. S2).

PLS-DA (of RP and HILIC in both positive and negative ionization modes) clearly discriminated the fasted individuals from those of the fed group (Fig. 2). Both groups were separated along the first component of the analysis, which explained 85–97% of the total variance. Individuals of the same group were distributed along the second PLS-DA component (2–10% of total variance). In the case of OPLS-DA, around 850 features from all four datasets were highlighted as discriminatory between fed and fasted fish with a P[corr] > 0.95 and a corrected P-value < 0.05 (see Fig. S3). Among them, up to 45 different compounds were elucidated as amino acids (4), oligopeptides (8), urea cycle-related metabolites (2), acylcarnitines (5), glutathione-related compounds (5), fatty acids (5), 3-hydroxyisovaleric acid, 3-methoxy-4-hydroxy-phenylglycol (MOPEG) sulphate and phospholipids (14), including phosphatidylcholines (PC) and lysoPC (Table 2).

Phospholipids were characterized by the presence of both the protonated molecule and sodium adduct in the positive LE spectra and their acetate adducts in negative LE spectra. As an example, Fig. S4 shows the m∕z 542.3248 (+0.1 mDa mass error) which corresponds to the protonated molecule, and m∕z 564.3074 (+1.2 mDa mass error) corresponds to the sodium adduct in RP+. MS/MS experiments were also carried out obtaining m∕z184.0735 (−0.4 mDa) (main product ions at 30 eV) in positive ionization mode. These m∕z ions were annotated as phosphocholine in agreement with fragmentation pathways of lysoPC (Xu et al., 2009; González-Domínguez, García-Barrera & Gómez-Ariza, 2014).

Carnitine-related compounds were also elucidated after MS/MS experiments by the presence of m∕z 144.1051 (C₇H₁₃NO${}_2^{\mathrm{+}}$, +2.6 mDa), 85.0303 (C₄H₅O${}_2^{\mathrm{+}}$, +1.3 mDa) and 60.0814 (C₃H₁₀N⁺, −0.3 mDa) as characteristic product ions of these compounds (Möder et al., 1997; Luci, Hirche & Eder, 2008) also observed in METLIN spectra for some of these compounds (Fig. S5).

Tandem mass spectrometry also provides relevant information for isomers differentiation. For example, γ-Glu-Ile and γ-Glu-Leu presented the same molecular formula and close retention times (Fig. S6). MS/MS experiments revealed very similar spectra at 10 and 20 eV with the exception of the m∕z 142.0499 which only appeared for the peak at 3.37 min. After comparing both spectra with METLIN, only γ-Glu-Leu spectrum showed this m∕z ion at 20 eV. The formation of C₄H${}_8^{+}$ have been observed and described in the literature much higher in isoleucine than in leucine (Squire, Beranová & Wesdemiotis, 1995), so in leucine spectra the neutral loss of C₄H₈ can be observed while it does not appears in isoleucine spectra. This strategy was followed for the rest of elucidated compounds.

Functional analysis of elucidated compounds

Biological significance of the concurrent up- or down-regulation of most of the elucidated metabolites during fasting clearly stated that food deprivation increased mobilization of body energy stores and improved the oxidative capacity of metabolic fuels, which paralleled the onset of specific changes in the cell redox-balance. In this regard, the increased mobilization of body fat stores in fasted individuals, exemplified by the loss of liver and adipose tissue mass, was linked to the consistent increase of circulating levels of five sub-products of L-carnitine (compounds 1–5 in Table 2), a carrier of fatty acids across the inner mitochondrial membrane for their subsequent beta-oxidation (Luci, Hirche & Eder, 2008; Ball, Urschel & Pencharz, 2007). At the molecular level, this was early substantiated in similar experimental conditions by a marked up-regulated expression of the two carnitine palmitoyltransferase variants (CPT1A, CPT1B) of the skeletal muscle of gilthead sea bream (Benedito-Palos, Ballester-Lozano & Pérez-Sánchez, 2014), which was encompassed by the increased expression of a high representation (25 enzyme subunits) of regulatory and assembly factors of the five enzyme complex units (Complex I–V) of the mitochondrial respiratory chain (Bermejo-Nogales, Calduch-Giner & Pérez-Sánchez, 2015). Microarray gene expression profiling of either glycolytic or aerobic muscle tissues of fish fed to maintenance ration also indicates that nutrient scarcity is by itself a major factor driving switches in muscle protein turnover and mitochondrial activity (Calduch-Giner et al., 2014). In the present study, this was reinforced by the consistent fasting increase of serum concentrations of urea cycle-related compounds (citrulline, ornithine, argininosuccinate and arginine). Of note, the activity of urea cycle enzymes is typically higher in carnivorous fish than in herbivorous and omnivorous fish species (Chiu, Austic & Rumsey, 1986), and our results highlighted that acyl-carnitine and urea cycle metabolites are specially sensitive to fasting-mediated changes in fatty acid and amino acid catabolism during negative energy balance.

Catecholamines are mobilized into fish circulation during a variety of stressful situations which require modulation of cardiorespiratory function or mobilization of energy reserves. The magnitude of change is dependent on the species and the type and intensity of stress imposed, although a wide range of stressors including hypoxia, hypercapnia, exhaustive and violent exercise, air exposure or anemia are considered strong activators of the hypothalamic-pituitary-interrenal (HPI) axis in fish (Reid, Bernier & Perry, 1998). This also applies to short-term fasting (De Pedro et al., 2003), and the observed increase of MOPEG sulphate, a metabolite of norepinephrine degradation, can be viewed as part of the adaptive response of the HPI axis to cope with fasting hypoglycemia through the activation of lipolysis and gluconeogenesis. This notion was supported by the increased levels of noradrenaline, and other catecholamine metabolites 3,4-dihydroxymandelaldehyde and 3,4-dihydroxyphenylethyleneglycol (DOPEGAL and DOPEG), detected in the refined search step of analysis.

The primary enzymatic antioxidant defense system of living organisms is the glutathione (GSH) redox system that reduces hydrogen peroxide and lipid hydroperoxides at the expense of oxidizing GSH to its disulfide form (GSSG). Once oxidized, GSH can be reduced back by glutathione reductase, using NADPH as an electron donor, and previous studies in gilthead seam bream indicate that either absolute GSH levels or the GSH/GSSG ratio are regulated by dietary oils, increasing the total plasma antioxidant capacity with the increased unsaturation index of dietary oils of marine origin (Saera-Vila et al., 2005). Likewise, total plasma antioxidant capacity is increased in hypoxic fish with a switch from oxidative phosphorylation (OXPHOS) to anaerobic glycolysis (Bermejo-Nogales, Calduch-Giner & Pérez-Sánchez, 2014), which results in reduced mitochondria oxygen consumption and enhanced NADH production from glycolysis (Frezza et al., 2011). Importantly, extension of life span is related in mammals and birds to low antioxidant levels and low rates of generation of reactive oxygen species (ROS) (Lykkesfeldt & Svendsen, 2007; Pamplona et al., 2008). Experimental evidence in rats also indicates that intermittent fasting affects redox balance in a tissue selective manner (Chausse et al., 2015), and our fish metabolomic study highlighted that the depletion of serum GSH during short-term fasting was closely related to changes in the Meister’s γ-glutamyl cycle with a key role in the recovery and delivery of cysteine in the body (Griffith, Bridges & Meister, 1978). This was supported by high circulating concentrations of ϒ-Glu-(Leu/Val/Ile) and pyroglutamic acid in the serum of fasted gilthead sea bream, whereas both GSH and ϒ-Glu-Cys were depleted (Fig. 3). This represents a complex trade-off with a reduced risk of oxidative stress, also highlighted by the decreased concentration of methionine sulfoxide, an oxidized form of methionine that is highly correlated with the risk of oxidative stress (Weissbach et al., 2002). In parallel, other short oligopeptides were either increased or decreased in the serum of fasted fish. It remains to be established if they have a physiological significance or are subproducts of protein hydrolysis.

The fatty acid composition of triacylglycerols (TAG) usually clears a close resemblance to dietary lipids (Benedito-Palos et al., 2010), whereas that of phospholipids is highly regulated and influenced in fish by environmental factors, including temperature and osmolarity (Los & Murata, 2004; Ibarz et al., 2005). The specific effects of ration size have been addressed in gilthead sea bream, and maintenance ratio significantly increased the retention of arachidonic acid (ARA) and docosahexanoic acid (DHA) in muscle phospholipids, whereas the fatty acid composition of storage lipids remained almost unchanged (Benedito-Palos et al., 2013). This lean muscle phenotype was linked in the present study to low plasma levels of linoleic acid and eicosapentaenoic acid, the precursors of ARA and DHA, respectively. At the same time, fasting induced an overall decrease of circulating lysoPC and glycerophosphocholine, whereas the effect on phosphatydylcholines was more selective depending of the composition of fatty acids. In any case, phospholipid metabolism is becoming highly regulated by feed intake at either blood or tissue level. Thus, phospholipids of skeletal muscle would act as a reservoir of long chain poly-unsaturated fatty acids with an enhanced expression of lipoprotein lipase-like, a TAG lipase isoform exclusive of fish linage which is highly expressed in muscle tissues and specifically up-regulated by feed restriction (Benedito-Palos et al., 2013; Rimoldi et al., 2015). This in turn would mediate, at least in part, the changes in the blood composition of phospholipid-related metabolites. This is extensive to trimethylamine N-oxide (TMAO), and high TMAO and choline concentrations are associated in humans with diabetes and advanced cardio-metabolic risk profile (Obeid et al., 2016). In agreement with this, the opposite pattern was found herein in fish under a negative energy balance, which reinforces the close metabolic association between interrelated pathways of phospholipid and oxidative metabolism. Recent metabolomics studies have highlighted a concurrent depletion of reduced GSH and glycerophosphocholine in the gills of another fish species, the golden grey mullet (Liza aurata), as a response to mercury toxicity (Cappello et al., 2016a; Cappello et al., 2016b). This finding also reinforces the view that the response to different challenges such as malnutrition or pollutants toxicity and susceptibility (Kokushi et al., 2016; Wang et al., 2016) can be assessed by the analysis of common responsive metabolites, opening the possibility for future screening of the fish general welfare status through selected biomarkers.

Lastly, major changes in vitamin status are related in our experimental model to biotin metabolism. In humans, the impairment of renal reclamation of biotin results in an elevated urine concentration of 3-hydroxyisovaleric acid (Mock et al., 1997). Accordingly, we found that the fasting increase of this metabolite was concurrent with a low biotin availability (highlighted by refining analysis), which reinforces the role of 3-hydroxyisovaleric acid as a biomarker of B7 vitamin deficiency in a wide range of vertebrate species, including fish.