ORIGINAL RESEARCH article
Identification of Plasmatic Biomarkers of Foie Gras Qualities in Duck by Metabolomics
- 1GenPhySE, Université de Toulouse, INRAE, ENVT, Castanet Tolosan, France
- 2Department of Animal Science, Faculty of Agriculture, Tarbiat Modares University, Tehran, Iran
- 3Toxalim, Université de Toulouse, INRA, ENVT, INP-Purpan, UPS, Toulouse, France
- 4Axiom Platform, MetaToul-MetaboHUB, National Infrastructure for Metabolomics and Fluxomics, Toulouse, France
- 5ASSELDOR, Station d’expérimentation appliquée et de démonstration sur l’oie et le canard, La Tour de Glane, Coulaures, France
The foie gras is an emblematic product of French gastronomy composed of waterfowl fatty liver. The organoleptic qualities of this product depend on the liver characteristics such as liver weight (LW) and technological yield (TY) at cooking. One of the main issues for producers is to classify the foie gras with high or low technological quality before cooking them. Thus the study aims at identifying biomarkers of these characteristics with non-invasive biomarkers in duck. 1H-NMR (nuclear magnetic resonance of the proton) analyses were performed on plasma of male mule ducks at different time points during the overfeeding period to obtain a large range of liver characteristics so as to identify plasmatic biomarkers of foie gras. We used two methods, one based on bucket data from the 1H-NMR spectra and another one based on the fingerprints of several metabolites. PLS analyses and Linear models were performed to identify biomarkers. We identified 18 biomarkers of liver weight and 15 biomarkers of technological yield. As these two quality parameters were strongly correlated (−0.82), 13 biomarkers were common. The lactate was the most important biomarker, the other were mainly amino acids. Contrary to the amino acids, the lactate increased with the liver weight and decreased with the technological yield. We also identified 5 biomarkers specific to LW (3 carbohydrates: glucuronic acid, mannose, sorbitol and 2 amino acids: glutamic acid and methionine) that were negatively correlated to liver weight. It was of main interest to identify 2 biomarkers specific to the technological yield. Contrary to the isovaleric acid, the valine was negatively correlated to the technological yield.
The foie gras of duck is a traditional product in France. It corresponds to a liver weighing more than 300 g (JORF, 1993) that is composed of 50 to 60% of lipids. During the cooking process, the lipids can melt resulting in a low cooking technological yield (TY). In addition to deteriorating the sensory qualities of the products, TY also directly influences the economic valuation of the foie gras. TY is controlled by the French legislation as it must exceed 70% in products labeled Foie gras (JORF, 1993).
The manufacturers classify the livers in function of their potential cooking TY using the liver weight (LW) and/or the liver texture. They can adapt their cooking protocol from pasteurization to sterilization to avoid the melting process and to improve TY. However the prediction of the potential cooking TY is not satisfying. As a result, a part of the cooked foie gras products are downgraded into products with a lower added value. One of the major challenges for manufacturers is therefore to predict TY of foie gras and to choose the best food processing. In addition, TY has a −0.80 genetic correlation with LW (Marie-Etancelin et al., 2011) and a −0.83 phenotypic correlation with LW (Bonnefont et al., 2019). But few authors tried to predict it before animal slaughtering. To predict phenotypes by identifying biomarkers, the plasma in an interesting fluid because plasma sampling is minimally invasive. Hermier et al. (1999) already analyzed some phospholipids of plasmatic very low density lipoprotein (VLDL) and high density lipoprotein (HDL) in control and overfed Landes goose. They showed a stronger quantity of phosphatidylcholine in VLDL and a stronger quantity of phosphatidylethanolamine in HDL in overfed goose (Hermier et al., 1999). Furthermore Tavernier et al. (2017) measured the triglycerides (TG) and non-esterified-fatty-acid (NEFA) concentrations in plasma of mule ducks during the overfeeding period. They showed a non-linear increase of TG and NEFA in parallel with LW. Recently, Pioche et al. (2019) measured the total cholesterol, the TG and the glucose in plasma of mule ducks during the overfeeding period. They highlighted a correlation between the total cholesterol and LW (0.88), between TG and LW (0.57) and between the glucose and LW (0.27) (Pioche et al., 2019). Thus these plasmatic molecules could be potential biomarkers of TY. The current study aims at identifying new biomarkers of foie gras LW and TY in mule duck.
The concentration of these metabolites directly reflects the biochemical activity and state of cells and tissues. All these metabolites can be potential biomarkers of TY that can be identify by untargeted metabolomics. Therefore, metabolomics has been strongly used to identify biomarkers (Quinones and Kaddurah-Daouk, 2009).
In the current study, plasmatic biomarkers of LW and TY were investigated by metabolomics approach using proton nuclear magnetic resonance (1H-NMR). The identification of LW and TY biomarkers could be used for both grading the liver of the mule ducks and increasing the added value of foie gras and improving knowledge in human hepatic steatosis.
Materials and Methods
Animal Experimental Design and Characteristics
The experimental protocol is clearly described in the previous paper (Bonnefont et al., 2019). Briefly, 120 male mule ducks (female Cairina moschata x male Anas platyrhynchos) were reared collectively until 12 weeks. They were overfed twice a day during 12 days with two different programs. But the overfeeding programs did not impact duck performances (Bonnefont et al., 2019) and were not taken into account in this current study. A total of 30 ducks were slaughtered at day 6, day 8, day 10 and day 12, 11 h after their last meal (Figure 1A). After 6 days of overfeeding the liver weight is expected to be superior to 300 g that corresponds to French definition of foie gras. After slaughtering the ducks, the liver was weighed. TY was determined after foie gras cooking as described in Rémignon et al. (2018), as following:
Figure 1. Description of the duck samples. Experimental design (A) Evolution of liver weight (B) and technological yield (C) during the overfeeding period (n = 16 or 17 at each time point).
TY (%) = [crude weight of liver in g – (cooked weight of liver in g - visible melted lipids in g)] × 100/[crude weight of liver in g]. At each time point, 16 or 17 ducks were selected among the 30 ducks for plasma analyses (n = 65 in total; Figure 1A). The mean and the variability of LW and TY in the group of 30 birds and in the subgroup were equivalent. LW and TY during the overfeeding period are presented in the Figures 1B,C.
Plasma Sampling and 1H-NMR Analysis
1 h before slaughtering blood was sampled from the venous occipital sinus in heparin-lithium tubes. The tubes were centrifuged at 3,400 g at 4°C during 10 min. The plasmas were collected and stored at −80°C for further analyses.
All plasma samples were extracted according to the procedure previously described to precipitate proteins and to avoid degraded NMR signal due to high lipids levels in fatty liver samples (Nagana Gowda et al., 2015). Briefly, 300 μL of plasma and 600 μL of methanol (MeOH) were mixed and incubated at −20°C for 20 min. After centrifugation (30 min at 11,000 g at 4°C), 800 μL of supernatant was isolated and evaporated using a vacuum concentrator (Concentrator plus, Eppenforf, Germany). Then, dried plasma extracts were diluted in 700 μL of NMR deuterium oxide (D2O) phosphate buffer (pH 7), containing sodium trimethylsilyl propionate (TSP) as an internal standard. The samples were then centrifuged at 4,600 g during 15 min at 4°C and 600 μl of supernatant were transferred into 5 mm NMR tubes. 1H-NMR analyses were performed on a Bruker Avance III HD spectrometer (Bruker Biospin, Rheinstetten, Germany) operating at a proton frequency of 600.13 MHz with an inverse detection 5-mm 1H-13C-15N cryoprobe.
1H-NMR spectra were acquired at 300 K with the nuclear Overhauser effect spectroscopy (NOESY) pulse sequence. A total of 128 transients were collected into 32,000 data points with a spectral width of 20 ppm and an acquisition time of 2.7 s. Prior to the Fourier transform procedure an exponential line-broadening of 0.3 Hz was applied to the free induction decay (FID).
Spectra Pre-processing and Statistical Analysis
1H-NMR spectra were analyzed by two methods: (i) a bucket method and (ii) a metabolite method.
(i) The bucket method is traditionally used to analyze 1H-NMR data. It consisted in converting the 1H-NMR spectra into a bucket intensity table with the Workflow4Metabolomics 3.3 online platform (1 Tremblay et al., 2014). (1) The spectra pre-processing included the following steps: solvent suppression (exclusion of the 5.1 to 4.5 ppm region corresponding to water signal), zero-filling, apodization, Fourier transform, phasing, baseline correction and calibration with TSP at 0.0 ppm. (2) The spectra alignment was performed. (3) The bucketing was defined with a 0.01 ppm interval from 0.5 to 10 ppm. The buckets corresponding to residual methanol from extraction (3.38 to 3.33 ppm) were removed. The intensity of a bucket corresponded to the intensity of the spectrum curve for that bucket. (4) The normalization was done with the whole spectrum intensity method as following:
Normalized bucket intensity = Raw bucket intensity/Whole spectrum intensity.
A table of bucket intensities was obtained with 65 rows corresponding to the animals and 758 columns corresponding to the buckets identified by their chemical shifts.
(ii) The recent metabolite method consists in converting the 1H-NMR spectra into a metabolite relative concentration table with the ASICS R package (R package version 4.0.22) that contains an automatic approach to identify and quantify metabolites in a complex 1H-NMR spectrum from their unique peak pattern (fingerprint) (Tardivel et al., 2017; Lefort et al., 2019). It is a relative quantification in function of the whole spectrum and the metabolite relative concentration has no unit. The metabolite database used consisted in spectra of 176 pure metabolites described in Tardivel et al. (2017). A total of 67 metabolites were identified and quantified in at least one sample. The methanol was removed as it was used to extract the metabolites. Then only 43 metabolites were kept for further analyses as they were present in at least 50% of the samples at least at one time point. Thus the final table of metabolite relative concentration contained 65 rows corresponding to the animals and 43 columns corresponding to the metabolites.
The bucket and metabolite relative concentration tables were analyzed with SIMCA P + software (version 12, Umetrics, AB, Umea, Sweden) for carrying out multivariate statistical analysis. First, the variables were pre-processed with a Pareto normalization. A principal component analysis (PCA) was performed for finding outliers. Then partial least square analyses (PLS) were performed to explain Y variables (LW and TY) by the X variables (bucket or metabolite data). The PLS scatter plots were drawn. The latent variables that corresponded to the scores t1, t2… were new variables summarizing the X variables. The first latent variable t1 explained the largest variation of the X space. Usually the scatter plot of t1 vs t2 is a window in the X space, displaying how the X observations are situated with respect to each other. Here, as only one latent variable was created, the PLS scatter plot represented t1 on the vertical axis vs sample identification on the horizontal axis. The intensity of the Y-values were indicated by the colors of the samples. The goodness-of-fits of the models were estimated by the proportion of cumulative explained variance (R2) for both the X variables (X = buckets or metabolites) and the Y variable (Y = LW or TY) and by the predictive ability of the model (Q2). Q2 is calculated by a cross-validation. The data were divided into 7 parts and each 1/7th in turn was removed. A model was built on the 6/7th of the data until all the data have been predicted. The predictive data are then compared with the original data and the sum of squared errors calculated for the whole dataset. The Root Mean Square Error of Estimation (RMSEE) was computed and indicated the fit of the observations to the model. The Root Mean Square Error after cross validation (RMSECv) was computed. The RMSECv has to be close to the RMSEE. The plot of the Y observed vs Y predicted values were drawn for each PLS model. With a good model all the points fall close to the 45 degree line. The validation of the PLS model was evaluated by comparing the goodness of fit (R2Y and Q2) of the original model with the goodness of fit of 500 models based on data where the order of the Y-observations has been randomly permuted, while the X-matrix (bucket or metabolite) has been kept intact. The permutation plot shows, for a selected Y-variable (LW or TY), on the vertical axis the values of R2Y and Q2 for the original model and for the Y-permuted models. The horizontal axis shows the correlation between the permuted Y-vectors and the original Y-vector for the selected Y (LW or TY). The original Y has a correlation 1.0 with itself. The criteria for validity of the original model is, all permuted Q2-values are lower than the original points and the regression line of the Q2-points intersects the vertical axis at, or below zero. The latent variables associated with interesting axes were analyzed by the variable importance in projection (VIP) method. The variables (buckets or metabolites) with a VIP superior to 1 were considered as “important”. Then, a one-by-one regression with either the LW effect or the TY effect was performed on the whole datasets. The P-values were corrected for multiple tests with the Benjamini-Hochberg correction using the R software (version 3.6.1) and named “BH P-values”. A variable was considered as “significant” when the BH P-value was inferior to 0.050 and “tended to be significant” when the BH P-Value was between 0.050 and 0.100. For the buckets with VIP superior to 1, the corresponding metabolites were identified manually by importing the chemical shift lists into the Human Metabolome Database3 (Wishart et al., 2009). All carbohydrates identified were D-carbohydrates and all amino acids were L-amino acids. To simplify the names of the metabolites, the “D-” and the “L-” were removed before the names of the carbohydrates and the amino acids, respectively. To confirm the identification of the metabolites, the 1H-NMR peaks of these metabolites were manually checked on the sample spectra of some plasma samples with TopSpin software (version 4.0, Bruker BioSpin, Germany). For each metabolite all 1H-NMR peaks were listed. For each 1H-NMR peak the VIP values and the BH P-values of the corresponding buckets were summarized by the number of buckets with VIP superior to 1 and by the range of BH P-values, respectively. The relative concentrations (RC) of a metabolite with the bucket data were estimated with a method adapted from Kostidis et al., 2017, by the following formula:
Metabolite RC j = mean [(intensity Peak ij)/(H number Peak ij)]
Where j represents a specific metabolite.
Where i represents each proton peak of the 1H-NMR spectrum of the j metabolite.
Where “intensity Peak ij” was computed as the sum of the bucket intensity of the i peak for the j metabolite.
Where “H number Peak ij” corresponds to the number of protons that corresponds to the i peak for the j metabolite.
Then the lists of the biomarkers obtained by the bucket method and the metabolite method were compared with a Venn diagram4. All the biomarkers identified by the bucket method and/or the metabolite method were considered as biomarkers. A network analysis based on the correlation of the biomarker RC and the Y-variable (LW or TY) was performed with the functions pls and network of the MixOmics R package (Lê Cao et al., 2009; Rohart et al., 2017) (R package version 18.104.22.168).
The rearing of male mule ducks that were overfed during 6 to 12 days enabled to obtain animals with large variability of performances. The liver of ducks weighed between 302.3 and 914.9 g and the technological yield was between 54.8 and 99.5% (Figures 1B,C). This experimental protocol mimicked the high variability of the liver characteristics that are present in the foie gras industry and was suitable to identify plasmatic biomarkers of foie gras quality by 1H-NMR analysis.
1-Identification of Plasmatic Biomarkers of Liver Weight of Foie Gras
First plasmatic biomarkers of LW were identified with the bucket method. After a PCA no outlier was detected (not shown). The PLS scatter plot had only one latent variable and the parameters were cumulative R2X = 0.304 and R2Y = 0.376. The projection of the samples highlighted an evolution of LW with the first latent variable on the vertical axis (Figure 2A). The prediction of the model was Q2 = 0.283. The original R2Y and Q2-values were higher than the ones obtained after permutation and the regression line of Q2-points intersected the vertical axis below zero (Figure 3A). The RMSEE and the RMSECv were close (145.3 and 153.3, respectively, Figure 4A). A group of 124 buckets with a VIP > 1 explained the latent variable (Supplementary Data 1). For the buckets with VIP > 1, the involved metabolites were identified. They corresponded to 14 metabolites summarized in Table 1. The relative concentrations of the metabolites were computed and the BH P-values of the metabolites were presented in bold in Table 1. In total, 13 out of the 14 metabolites were statistically significant (BH P-Value < 0.05) and were further considered as biomarkers of LW. For each biomarker, the numbers of important peaks that contained at least one bucket with VIP > 1 in comparison to the numbers of 1H-NMR peaks were indicated in parenthesis. The biomarkers of LW were 3 carbohydrates: lactate (HMDB0000190, with 2/2 important peaks), mannose (HMDB0000169 with 8/9) and sorbitol (HMDB0000247 with 5/5), 8 amino acids: alanine (HMDB0000161 with 2/2), arginine (HMDB0000517 with 4/4), glutamic acid (HMDB0000148 with 8/8), glutamine (HMDB0000641 with 3/3), isoleucine (HMDB0000172 with 6/6), lysine (HMDB0000182 with 6/12), methionine (HMDB0000696 with 3/3) and pyroglutamic acid (HMDB0000267 with 4/4) and also glycerol (HMDB0000131 with 3/3) and methylmalonic acid (HMDB00202 with 2/2; Table 1).
Figure 2. PLS score plots for foie gras (A) for liver weight with the bucket method (R2X = 0.304, R2Y = 0.376, Q2 = 0.283) and (B) for liver weight with the metabolite method (R2X = 0.339, R2Y = 0.323, Q2 = 0.269), (C) for technological yield with the bucket method (R2X = 0.313, R2Y = 0.298, Q2 = 0.175) and (D) for technological yield with the metabolite method (R2X = 0.326, R2Y = 0.254, Q2 = 0.159). The numbers correspond to the identification of the samples and the colors to the liver weight values. The legend is indicated on the right of the figures.
Figure 3. Permutation plots. (A) for liver weight with the bucket method, (B) for liver weight with the metabolite method, (C) for technological yield with the bucket method, (D) for technological yield with the metabolite method. 500 permutations were performed. The permutation plot shows, for a selected Y-variable, on the vertical axis the values of R2Y and Q2 for the original model and for the Y-permuted models. The horizontal axis shows the correlation between the permuted Y-vectors and the original Y-vector for the selected Y. The original Y has a correlation 1.0 with itself.
Figure 4. Plots of predicted vs observed data. (A) for liver weight with the bucket method, (B) for liver weight with the metabolite method, (C) for technological yield with the bucket method, and (D) for technological yield with the metabolite method. The Root Mean Square Error of Estimation (RMSEE) and the Root Mean Square Error after cross validation (RMSECv) are indicated.
In parallel, the metabolite method was applied and the 65 spectra were converted into a table of 43 metabolite intensities with the ASICS R package. No outlier was detected by PCA (not shown). The PLS scatter plot had only one latent variable and the parameters were cumulative R2X = 0.339 and R2Y = 0.323 (Figure 2B). The prediction of the model was Q2 = 0.269. The original R2Y and Q2-values were higher than the ones obtained after permutation and the regression line of Q2-points intersected the vertical axis below zero (Figure 3B). The RMSEE and the RMSECv were close (10.7 and 11.4, respectively, Figure 4B). The evolution of LW was well represented on the vertical axis corresponding to the first latent variable (Figure 2B). Only 9 metabolites had a VIP > 1 and 7 metabolites were significant (BH P-value < 0.05) and 1 metabolite tended to be significant (BH P-value < 0.1) (Table 2). Only these 8 metabolites were considered as biomarkers of LW, including 2 carbohydrates: glucuronic acid (VIP = 1.03, BH P-value = 0.006) and lactate (VIP = 4.21, BH P-value = 0.008) and 6 amino acids: alanine (VIP = 1.30, BH P-value = 0.024), glutamine (VIP = 1.10, BH P-value = 0.075), glycine (VIP = 1.21, BH P-value = 0.007), leucine (VIP = 1.12, BH P-value = 0.002), proline (VIP = 2.50, BH P-value < 0.001) and serine (VIP = 1.08, BH P-value = 0.033; Table 2).
In conclusion, 18 biomarkers were identified for LW whose 3 biomarkers identified by both the bucket method and the metabolite method (lactate, alanine and glutamine), 5 ones identified only by the metabolite method (glucuronic acid, glycine, leucine, proline and serine) and 10 ones identified only by the bucket method (mannose, sorbitol, arginine, glutamic acid, isoleucine, lysine, methionine, pyroglutamic acid, glycerol and methylmalonic acid; Figure 5A). For all the 18 biomarkers, their RC were computed with the bucket data and the plots of their RC in function of LW were presented in Figure 6A. The correlation network between LW and the 18 biomarkers was presented in Figure 7A. LW was positively correlated to lactate (0.62) and negatively correlated to other carbohydrates: mannose (−0.78), sorbitol (−0.76) and glucuronic acid (−0.63) and to all the amino acids: alanine (−0.97), glutamic acid (−0.96), glutamine (−0.93), arginine (−0.92), glycine (−0.87), lysine (−0.84), serine (−0.82), methionine (−0.79), leucine (−0.65), proline (−0.64), pyroglutamic acid (−0.32), isoleucine (−0.29) and to glycerol (−0.96) and methylmalonic acid (−0.29; Figure 6A).
Figure 5. Comparisons of biomarker lists with Venn diagram. (A) Biomarkers of liver weight (LW) identified by the bucket method and by the metabolite method (with VIP > 1 and BH P-value < 0.1), (B). Biomarkers of technological yield (TY) identified by the bucket method and by the metabolite method. (C) Biomarkers of liver weight and technological yield identified by at least one method.
Figure 6. Plots of biomarker relative concentration in function of liver weight (LW) (A) or technological yield (TY) (B). The relative concentrations are computed with the bucket data and have no unit. The regression curves are in red.
Figure 7. Correlation networks of foie gras biomarkers (A) of liver weight (Y represents the liver weight) and (B) of technological yield (Y represents the technological yield). The relative concentration used for the correlations are calculated with the bucket data.
2- Identification of Plasmatic Biomarkers of Technological Yield of Foie Gras
A PCA was first performed but no outlier was detected (not shown). The PLS scatter plot to explain TY had only one latent variable and the parameters were cumulative R2X = 0.313 and R2Y = 0.298. The projection of the samples highlighted an evolution of TY with the first latent variable on the vertical axis (Figure 2C). The prediction of the model was Q2 = 0.175. The original R2Y and Q2-values were higher than the ones obtained after permutation and the regression line of Q2-points intersected the vertical axis below zero (Figure 3C). The RMSEE and the RMSECv were close (151.3 and 154.8, respectively, Figure 4C). A group of 128 buckets with a VIP > 1 explained the latent variable (Supplementary Data 2). They corresponded to 16 metabolites (Table 3). In total, 9 out of the 16 metabolites were statistically significant (BH P-Value < 0.05) and 1 tended to be significant (BH P-value = 0.070). These 10 biomarkers were further considered as biomarkers of TY. For each biomarker, the numbers of important peaks that contained at least one bucket with VIP > 1 and the total numbers of 1H-NMR peaks were indicated in parenthesis. The biomarkers of TY were 1 carbohydrate: lactate (HMDB0000190 with 2/2 important peaks), 6 amino acids: alanine (HMDB0000161 with 2/2), arginine (HMDB0000517 with 4/4), glutamine (HMDB0000641 with 3/3), isoleucine (HMDB0000172 with 6/6), lysine (HMDB0000182 with 8/12) and pyroglutamic acid (HMDB0000267 with 4/4) and 3 other organic compounds: glycerol (HMDB0000131 with 3/3), isovaleric acid (HMDB0000718 with 3/3) and methylmalonic acid (HMDB00202 with 2/2; Table 3).
Table 3. List of the 15 biomarkers of foie gras technological yield identified with the bucket method.
In parallel, the metabolite method was performed. No outlier was detected by PCA (not shown). The PLS scatter plot had only one latent variable and the parameters were cumulative R2X = 0.326 and R2Y = 0.254 (Figure 2D). The prediction of the model was Q2 = 0.159. The original R2Y and Q2-values were higher than the ones obtained after permutation and the regression line of Q2-points intersected the vertical axis below zero (Figure 3D). The RMSEE and the RMSECv were close (11.0 and 11.5, respectively, Figure 4D). The latent variable on the vertical axis explained the evolution of TY (Figure 2D) and 9 metabolites had a VIP values superior to 1 (Table 4). But only 7 metabolites were significant (BH P-value < 0.05) and 1 metabolite tended to be significant (BH P-value = 0.080). Only these 8 metabolites were considered as biomarkers of TY. There were only 1 carbohydrate: lactate (VIP = 4.03, BH P-value = 0.046) and 7 amino acids: alanine (VIP = 1.26, BH P-value = 0.080), glutamine (VIP = 1.65, BH P-value = 0.043), glycine (VIP = 1.39, BH P-value = 0.010), leucine (VIP = 1.38, BH P-value < 0.001), proline (VIP = 2.68, BH P-value < 0.001), serine (VIP = 1.33, BH P-value = 0.043) and L-valine (VIP = 1.13, BH P-value = 0.010; Table 4).
Table 4. List of the 9 biomarkers of foie gras technological yield identified with the metabolite method.
In conclusion, 15 biomarkers were identified for TY: 3 biomarkers were commonly identified by the bucket and the metabolite methods (lactate, alanine and glutamine), 5 ones were identified only by the metabolite method (glycine, leucine, proline, serine and valine) and 7 ones were identified by the bucket method (arginine, isoleucine, lysine, pyroglutamic acid, glycerol, isovaleric acid and methylmalonic acid; Figure 5B). For the 15 biomarkers of TY, their RC were computed with the bucket data and the plots of their RC in function of TY were presented in Figure 6B. The correlation network between TY and the biomarker RC was presented in Figure 7B. TY was negatively correlated to lactate (−0.58) and valine (−0.58) and positively correlated to glutamine (0.91), alanine (0.85), proline (0,83), lysine (0.81), arginine (0.73), glycine (0.65), pyroglutamic acid (0.61), isoleucine (0.58), serine (0.57) and leucine (0.35; Figure 6B).
Consequently, the results of the 1H-NMR analysis identified 18 biomarkers for the liver weight of foie gras and 15 for its technological yield (Table 5). As the phenotypic correlation between LW and TY was strong (−0.82, P-value < 0.001), 13 biomarkers were common to LW and TY, of which 1 carbohydrate: lactate and 10 amino acids: alanine, arginine, glutamine, glycine, isoleucine, leucine, lysine, proline, pyroglutamic acid, serine and also glycerol and methylmalonic acid (Figure 5C and Table 5). All these biomarkers were negatively correlated to LW and positively correlated to TY, except the lactate. Thus the small livers with high technological yield were characterized by a low concentration of lactate and a strong concentration of the other biomarkers and it was the contrary for the heavy livers with low technological yield. Moreover, 5 biomarkers were specific to LW: glucuronic acid, mannose, sorbitol, glutamic acid and methionine. All were negatively correlated to LW. Thus their concentrations were higher in small livers than in heavy livers. In addition, two biomarkers were specific to technological yield: the valine was negatively correlated to TY and the isovaleric acid was positively correlated to TY. As a result, the livers with high TY were characterized by a high concentration of isovaleric acid and a small concentration of valine and vice versa (Table 5).
NMR Methodology Discussion
A constant volume of 300 μL was sampled and prepared for further 1H-NMR analysis. Thus the metabolite quantity was supposed to be equivalent between all samples. To build the bucket intensity table, the intensity of each bucket was normalized by the intensity of the whole spectrum that was assumed to be equivalent in all samples, that may insert an error in the bucket intensity values. Moreover, to compute the relative concentration of a metabolite, the intensity of each 1H-NMR peak was divided by the number of proton for this peak as it was proposed by Kostidis et al. (2017). Then, the mean of the intensity of all 1H-NMR peaks for this metabolite was computed. But contrary to Kostidis et al. (2017), the concentration of a metabolite was not divided by the concentration of a standard, as the normalization of the buckets was done by the intensity of the whole spectrum, thus the bucket intensity was already a relative bucket intensity.
In addition, one of the bias of the 1H-NMR analysis is that several 1H-NMR peaks of different metabolites can have the same chemical shift. Thus there is an error in the estimation of the relative concentration of a metabolite, as the relative concentration of a peak of a metabolite can be increased by the presence of a peak of another metabolite at the same chemical shift. The fact that the relative concentration of a metabolite was computed as the mean of all the 1H-NMR of this metabolite may reduce this bias as all the 1H-NMR peaks could not be affected by the peaks of the same metabolites. Moreover, the comparison of the bucket method and metabolite method obtained with the ASICS package (Lê Cao et al., 2009; Rohart et al., 2017) enabled us to be more confident in the current results.
First O-PLS models were developed to analyze 1H-NMR spectra with both the bucket method and the metabolite method (not shown) but as it is known to lead to a severe overfitting of the data, we decided to use PLS models. The R2Y and Q2-values of our PLS models were small. Thus to validate the models, we decided to perform 500 permutation tests to compare the original R2Y and Q2-values to the ones obtained after permutations and we computed the RMSEE and RMSECv parameters. Moreover, the biomarkers that we identified were consistent with the literature of hepatic steatosis of other animal models, giving insights on our results on metabolic mechanisms of hepatic steatosis.
In a study to search for biomarkers of non-alcoholic fatty liver disease (NAFLD) in human, the plasmatic content of glucose was equivalent between the control and the steatosis groups but it was increased in a non-alcoholic steatosis hepatitis (NASH) group (Kalhan et al., 2011). In our study, the plasmatic glucose content was negatively correlated with LW (−0.66) but it was not significant. This difference was probably due to the fact that glucose has many 1H-NMR peaks and is not easy to estimate its relative concentration in plasma with 1H-NMR method.
The liver of mule ducks responds to overfeeding by increasing its glucose uptake capacity (Pioche et al., 2019). The glucose may be converted into sorbitol by the aldose reductase an enzyme that can reduce carbonyl function into alcohol function. The up-regulation of this enzyme induced from high glucose intake led to a strong increase of sorbitol in hepatocytes, resulting in the elevation of intracellular triglycerides (Hotta et al., 2019). In our study, the plasmatic sorbitol was negatively correlated with LW (−0.76). This result confirms the study of a rabbit model of hepatic steatosis in which a high cholesterol diet increased the presence of sorbitol and the activity of sorbitol dehydrogenase enzyme in the blood (Birkner et al., 2007). The sorbitol was already identified as a potential biomarker of alcoholic steatosis in vivo and in vitro in mice (Guo et al., 2018).
In our study, the plasmatic lactate content obtained the highest VIP values (4.21 and 4.03 for LW and TY respectively) therefore it was the metabolite with the strongest importance to draw the first latent variable that separate the livers in function of their LW or TY. It was enhanced when LW increased and reduced when TY increased (Pearson correlation of 0.62 and −0.58, respectively). Similar results were highlighted in the liver, because the liver lactate content was low in low fat loss livers that corresponded to high TY livers (Bonnefont et al., 2014). In addition, the lactate content was already identified as biomarker of liver steatosis in serum for mice (Li et al., 2011), in plasma of human (Kalhan et al., 2011) and in blood of human (Lin et al., 2020).
The amino acid metabolism was strongly impacted by the evolution of fatty livers as many amino acids were identified as biomarkers (10 of LW and TY, 2 of LW only and 1 of TY only). For all the amino acid biomarkers of LW when LW increased their plasmatic contents decreased as their correlations with LW were negative. Thus, the plasmatic amino acid concentration was higher in ducks with small livers with high TY. This corroborates with previous results as the metabolism of livers with low fat loss was oriented toward anabolism contrary to the livers with high fat loss (Theron et al., 2011; Bonnefont et al., 2014).
Many previous studies highlighted the reduction of the amino acid metabolism in animals or humans with steatosis hepatisis. For instance, in obese mice, the role of leucine on the regulation of energy metabolism was shown (Bruckbauer and Zemel, 2014; Fu et al., 2015). The role of leucine, isoleucine and valine deficiency revealed a decreasing in lipogenic gene expression in the liver in mice (Du et al., 2012). Contrary to our results, Kalhan et al. (2011) showed an increase in plasma content for glutamate and isoleucine and an insignificant evolution of leucine and valine in steatosis patients, but these evolutions were significant in NASH patients when compared to controls (Kalhan et al., 2011). Human NAFLD patients had much higher concentrations of proline in urines (Dong et al., 2017) whereas we found a higher proline concentration in plasma of ducks with low LW.
Plasmatic biomarkers of NAFLD have already been searched in many researches. For instance, Li et al. (2011) identified 3 serum biomarkers (glucose, lactate, glutamate/glutamine) to diagnose NAFLD at various stages in mice (Li et al., 2011). High plasmatic levels of isoleucine, alanine and glutamate were also associated with NAFLD severity and glutamate was the top plasmatic amino acid biomarker of obesity (Sookoian and Pirola, 2018). Glutamine supplementation reduced oxidative stress and NAFLD, and increased glucose metabolism in insulin resistant Ob/Ob mice (Leite et al., 2018). Thus the negative correlation of glutamic acid with LW (-0.96) was coherent because the heavy livers were associated to oxidative stress (Bonnefont et al., 2014). The anti-diabetic effect of pyroglutamic acid that is obtained from the glutamic acid were reported in type 2 diabetic rats and mice (Yoshinari and Igarashi, 2011). Moreover, the plasmatic contents of valine, leucine and isoleucine were significantly reduced in animals and patients with hepatic encephalopathy (Fischer et al., 1975). Our results were consistent for leucine and isoleucine (correlation with LW −0.65 and −0.29, respectively) and we did not identify valine as a biomarker of LW, but as a biomarker of TY. A lower circulating glycine levels was observed in association with hepatic insulin resistance in NAFLD patients (Alves et al., 2019) in accordance with our findings where glycine was negatively correlated with LW (−0.87) and positively correlated with TY (0.65). Moreover, the quantity of some metabolites impact directly the development of liver steatosis. Diets deficient in labile methyl groups (choline, methionine, betaine, folate) produced fatty liver (Mato et al., 2008) which validates the negative correlation of methionine with LW (−0.79), and a 5% L-Lysine diet developed fatty livers in rats (Hevia and Visek, 1980). In addition, the evolution of amino acid metabolism that we observed was consistent with the utilization by several authors of plasmatic aminotransferases as biomarkers of liver metabolic functioning (Sookoian and Pirola, 2015).
To our knowledge, there is no reference of methylmalonic acid implication in hepatic steatosis. The role of the isovaleric acid in the metabolism of hepatic steatosis is not clearly describe in the literature. However, in a study about the impact of specific starch to reverse the weight gain and hepatic steatosis induced by high fat diet in mice, the colonic isobutyric acid and isovaleric acid levels were decreased by half compared to the high-fat group of mice (Zhang et al., 2020). Even if we studied the plasma and not the colon, the high concentration of plasmatic isovaleric acid in high technological livers seemed to be consistent.
Conclusion – Perspectives
This study presents the first analysis of plasmatic biomarkers of duck foie gras qualities with a large approach. We identified 18 biomarkers of liver weight and 15 biomarkers of technological yield of foie gras that were mainly lactate and amino acids. As these two quality parameters were strongly correlated, 13 metabolites were biomarkers of both LW and TY. The lactate was the most important biomarker, it increased with LW and decreased with TY. On the contrary the other biomarkers that were mainly amino acids were negatively correlated to LW and positively correlated to TY. We also identified 5 biomarkers specific to LW (3 carbohydrates and 2 amino acids) that were negatively correlated to LW. It was of main interest to identify 2 biomarkers specific to the technological yield. Contrary to the isovaleric acid, the valine was negatively correlated to TY. To predict the technological yield, these metabolites could be measured in order to optimize the valorization of foie gras. But further studies are required to analyze the robustness of the model based on these biomarkers.
Moreover, it was previously shown that plasmatic NEFA, TG and cholesterol were correlated to liver weight (Tavernier et al., 2017; Pioche et al., 2019). In other animal models of hepatic steatosis, lipophilic biomarkers were identified in rats (Goda et al., 2018) and specifically ceramids in mice (Dong et al., 2017) and adolescents (Maldonado-Hernández et al., 2017). Another perspective of the present study could be to search for biomarkers in the lipid fraction of the plasma.
Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
The animal study was reviewed and approved by Comité d’éthique en expérimentation animale - 073 Comité d’éthique Aquitaine Poissons Oiseaux.
ZM performed the statistical analyses and wrote the manuscript. JA made the project to be funding and supervised the animal rearing, overfeeding, and slaughtering. NM-G and CC realized the NMR analyses and the pre-processing of the data. BL helped a lot to analyze the NMR data. MM and AB helped for the analysis and the interpretation of the results. AM and CB supervised the study. All authors contributed to the article and approved the submitted version.
Interprofession Foie gras (CIFOG), and Aquitaine Regional and Dordogne Departmental Councils for funding the project.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors thank the staff of the Station d’expérimentation appliquée et de démonstration sur l’oie et le canard (Dordogne, France) for their excellent supervision of this study, especially F. Lavigne and C. Mondoux. We also thank the staff of the UMR 1388 Génétique Physiologie et Systèmes d’Elevage, in particular M. Bouillier-Oudot and C. Molette for their help in designing the experimental project and getting funding, H. Manse and S. Seidlinger for liver sampling and E. Cobo for blood sampling. The authors also acknowledge the Interprofession Foie gras (CIFOG), the Aquitaine Regional and Dordogne Departmental Councils for funding the project.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.628264/full#supplementary-material
Supplementary Data 1 | Identification of biomarkers of liver weight of duck foie gras by PLS and linear models applied on buckets of 1H-NMR spectra. εThe PLS model to describe the liver weight with bucket data was plotted. The first latent enabled to separate the fatty livers in function of their liver weight. The VIP of the buckets involved in the first latent are presented. ζFor each bucket the effect of the bucket intensity on the liver weight was tested by a linear model with R software. ηThe P-Values were corrected with the Benjamini-Hochberg procedure and named BH P-Values.
Supplementary Data 2 | Identification of biomarkers of technological yield of duck foie gras by PLS and linear models applied on buckets of 1H-NMR spectra. εThe PLS model to describe the technological yield with bucket data was plotted. The first latent enabled to separate the fatty livers in function of their technological yield. The VIP of the buckets involved in the first latent are presented. ζFor each bucket the effect of the bucket intensity on the technological yield was tested by a linear model with R software. ηThe P-Values were corrected with the Benjamini-Hochberg procedure and named BH P-Values.
- ^ https://workflow4metabolomics.org/
- ^ https://bioconductor.org/packages/ASICS/
- ^ http://hmdb.ca/
- ^ https://bioinfogp.cnb.csic.es/tools/venny/index.html
- ^ http://www.bioconductor.org/packages/release/bioc/html/mixOmics.html
Birkner, E., Zalejska-Fiolka, J., Kasperczyk, A., Kasperczyk, S., Grucka-Mamczar, E., Stawiarska-Piêta, B., et al. (2007). The influence of methionine, selenomethionine, and vitamin E on liver metabolic pathways and steatosis in high-cholesterol fed rabbits. Biol. Trace Elem. Res. 120, 179–194. doi: 10.1007/s12011-007-0070-4
Bonnefont, M. D. C., Guerra, A., Théron, L., Molette, C., Canlet, C., and Fernandez, X. (2014). Metabolomic study of fatty livers in ducks: identification by 1H-NMR of metabolic markers associated with technological quality. Poultry Sci. 93, 1542–1552. doi: 10.3382/ps.2013-03546
Bonnefont, M. D. C., Molette, C., Lavigne, F., Manse, H., Bravo, C., Lo, B., et al. (2019). Evolution of liver fattening and foie gras technological yield during the overfeeding period in mule duck. Poultry Sci. 98, 5724–5733. doi: 10.3382/ps/pez359
Bruckbauer, A., and Zemel, M. B. (2014). Synergistic effects of polyphenols and methylxanthines with leucine on AMPK/Sirtuin-mediated metabolism in muscle cells and adipocytes. PLoS One 9:e89166. doi: 10.1371/journal.pone.0089166
Dong, S., Zhan, Z. Y., Cao, H. Y., Wu, C., Bian, Y. Q., Li, J. Y., et al. (2017). Urinary metabolomics analysis identifies key biomarkers of different stages of nonalcoholic fatty liver disease. World J. Gastroenterol. 23, 2771–2784. doi: 10.3748/wjg.v23.i15.2771
Du, Y., Meng, Q., Zhang, Q., and Guo, F. (2012). Isoleucine or valine deprivation stimulates fat loss via increasing energy expenditure and regulating lipid metabolism in WAT. Amino Acids 43, 725–734. doi: 10.1007/s00726-011-1123-8
Fischer, J. E., Funovics, J. M., Aguirre, A., James, J. H., Keane, J. M., Wesdorp, R. I., et al. (1975). The role of plasma amino acids in hepatic encephalopathy. Surgery 78, 276–290. doi: 10.5555/uri:pii:0039606075902147
Fu, L., Li, F., Bruckbauer, A., Cao, Q., Cui, X., Wu, R., et al. (2015). Interaction between leucine and phosphodiesterase 5 inhibition in modulating insulin sensitivity and lipid metabolism. Diabetes Metab. Syndr. Obes. 8, 227–239. doi: 10.2147/DMSO.S82338
Goda, K., Saito, K., Muta, K., Kobayashi, A., Saito, Y., and Sugai, S. (2018). Ether-phosphatidylcholine characterized by consolidated plasma and liver lipidomics is a predictive biomarker for valproic acid-induced hepatic steatosis. J. Toxicol. Sci. 43, 395–405. doi: 10.2131/jts.43.395
Guo, C., Huang, J., Wang, Y., Shi, C., Gao, J., Hong, Y., et al. (2018). Aldose reductase inhibitor protects mice from alcoholic steatosis by repressing saturated fatty acid biosynthesis. Chem. Biol. Interact. 287, 41–48. doi: 10.1016/j.cbi.2018.04.002
Hermier, D., Salichon, M., Guy, G., and Peresson, R. (1999). Differential channelling of liver lipids in relation to susceptibility to hepatic steatosis in the goose. Poultry Sci. 78, 1398–1406. doi: 10.1093/ps/78.10.1398
Hotta, N., Kawamura, T., and Umemura, T. (2019). Are the polyol pathway and hyperuricemia partners in the development of non−alcoholic fatty liver disease in diabetes? J. Diabetes Investig. 11, 786–788. doi: 10.1111/jdi.13190
Kalhan, S. C., Guo, L., Edmison, J., Dasarathy, S., McCullough, A. J., Hanson, R. W., et al. (2011). Plasma metabolomic profile in nonalcoholic fatty liver disease. Metabolism 60, 404–413. doi: 10.1016/j.metabol.2010.03.006
Kostidis, S., Addie, R. D., Morreau, H., Mayboroda, O. A., and Giera, M. (2017). Quantitative NMR analysis of intra-and extracellular metabolism of mammalian cells: a tutorial. Anal. Chim. Acta 980, 1–24. doi: 10.1016/j.aca.2017.05.011
Lefort, G., Liaubet, L., Canlet, C., Tardivel, P., Père, M. C., Quesnel, H., et al. (2019). ASICS: an R package for a whole analysis workflow of 1D 1H NMR spectra. Bioinformatics 35, 4356–4363. doi: 10.1093/bioinformatics/btz248
Leite, J. S. M., Takahashi, H., Junior, J. D., Cruzat, V. F., and Carpinelli, A. R. (2018). “Free and dipeptide forms of L-glutamine supplementation attenuate parameters of oxidative stress and nonalcoholic fatty liver disease (NAFLD), and improve glucose metabolism in insulin resistant Ob/Ob mice,” in Proceedings of the 20th European Congress of Endocrinology, Vol. 56, (Bristol: BioScientifica), doi: 10.1530/endoabs.56.P579
Li, H., Wang, L., Yan, X., Liu, Q., Yu, C., Wei, H., et al. (2011). A proton nuclear magnetic resonance metabonomics approach for biomarker discovery in nonalcoholic fatty liver disease. J. Proteome Res. 10, 2797–2806. doi: 10.1021/pr200047c
Lin, H., Zhu, L., Baker, S. S., Baker, R. D., and Lee, T. (2020). Secreted phosphoglucose isomerase is a novel biomarker of nonalcoholic fatty liver in mice and humans. Biochem. Biophys. Res. Commun. 529, 1101–1105.
Maldonado-Hernández, J., Saldaña-Dávila, G. E., Piña-Aguero, M. I., Núñez-García, B. A., and López-Alarcón, M. G. (2017). Association between plasmatic ceramides profile and AST/ALT ratio: C14: 0 ceramide as predictor of hepatic steatosis in adolescents independently of obesity. Can. J. Gastroenterol. Hepatol. 2017:3689375. doi: 10.1155/2017/3689375
Marie-Etancelin, C., Basso, B., Davail, S., Gontier, K., Fernandez, X., Vitezica, Z. G., et al. (2011). Genetic parameters of product quality and hepatic metabolism in fattened mule ducks. J. Anim. Sci. 89, 669–679. doi: 10.2527/jas.2010-3091
Pioche, T., Skiba, F., Bernadet, M. D., Seiliez, I., Massimino, W., Houssier, M., et al. (2019). Kinetic study of the expression of genes related to hepatic steatosis, global intermediate metabolism and cellular stress during overfeeding in mule ducks. bioRxiv [preprint] doi: 10.1101/690156
Rémignon, H., Yahia, R. B. H., Marty-Gasset, N., and Wilkesman, J. (2018). Apoptosis during the development of the hepatic steatosis in force-fed ducks and cooking yield implications. Poultry Sci. 97, 2211–2217. doi: 10.3382/ps/pey054
Rohart, F., Gautier, B., Singh, A. L., and Lê Cao, KA. (2017). mixOmics: an R package for ‘omics feature selection and multiple data integration. PLoS Comput. Biol. 13:e1005752. doi: 10.1371/journal.pcbi.1005752
Sookoian, S., and Pirola, C. J. (2015). Liver enzymes, metabolomics and genome-wide association studies: From systems biology to the personalized medicine. World J. Gastroenterol. 21, 711–725. doi: 10.3748/wjg.v21.i3.711
Sookoian, S. C., and Pirola, C. J. (2018). The nonalcoholic steatohepatitis metabotype: imbalance of circulating amino acids and transamination reactions reflect impaired mitochondrial function. Hepatology 67, 1177–1178. doi: 10.1002/hep.29705
Tardivel, P. J., Canlet, C., Lefort, G., Tremblay-Franco, M., Debrauwer, L., Concordet, D., et al. (2017). ASICS: an automatic method for identification and quantification of metabolites in complex 1D 1 H NMR spectra. Metabolomics 13:109. doi: 10.1007/s11306-017-1244-5
Tavernier, A., Ricaud, K., Bernadet, M. D., Davail, S., and Gontier, K. (2017). Kinetics of expression of genes involved in glucose metabolism after the last meal in overfed mule ducks. Mol. Cell. Biochem. 430, 127–137. doi: 10.1007/s11010-017-2960-x
Theron, L., Fernandez, X., Marty-Gasset, N., Pichereaux, C., Rossignol, M., Chambon, C., et al. (2011). Identification by proteomic analysis of early post-mortem markers involved in the variability in fat loss during cooking of mule duck “foie gras”. J. Agric. food Chem. 59, 12617–12628. doi: 10.1021/jf203058x
Tremblay, M., Martin, J. F., and Goulitquer, S. (2014). Workflow4Metabolomics: a collaborative research infrastructure for computational metabolomics. Bioinformatics 31, 1493–1495. doi: 10.1093/bioinformatics/btu813
Zhang, Y., Chen, L., Hu, M., Kim, J. J., Lin, R., Xu, J., et al. (2020). Dietary type 2 resistant starch improves systemic inflammation and intestinal permeability by modulating microbiota and metabolites in aged mice on high-fat diet. Aging (Albany N. Y.) 12:9173. doi: 10.18632/aging.103187
Keywords: liver, foie gras, quality, biomarker, plasma, metabolomics
Citation: Mozduri Z, Marty-Gasset N, Lo B, Masoudi AA, Morisson M, Canlet C, Arroyo J, Bonnet A and Bonnefont CMD (2021) Identification of Plasmatic Biomarkers of Foie Gras Qualities in Duck by Metabolomics. Front. Physiol. 12:628264. doi: 10.3389/fphys.2021.628264
Received: 11 November 2020; Accepted: 04 January 2021;
Published: 12 February 2021.
Edited by:Massimiliano Petracci, University of Bologna, Italy
Reviewed by:Carlo Mengucci, University of Bologna, Italy
Yingping Xiao, ZheJiang Academy of Agricultural Sciences, China
Copyright © 2021 Mozduri, Marty-Gasset, Lo, Masoudi, Morisson, Canlet, Arroyo, Bonnet and Bonnefont. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.