The test design was based on our previous study (Gong et al., 2020). In brief, sixty healthy male Landes geese of similar weight (65 days old) were randomly assigned to one of three groups (20 geese in each treatment), including one control group and two experimental groups. Treatment included: (1) D0 group: geese fed a regular diet, without any special treatment; (2) D7 group: geese overfeeding on d 19 and lasting for 7 days; (3) D25 group: geese overfeeding on d 1 and lasting for 25 days. The experimental diet was composed of 98% corn, 1% vegetable oil, and 1% salt. The geese were fasted for one night (water was provided) at the end of the overfeeding trial. The next morning, 10 geese of similar weight (90 days old) were chosen from the per groups for liver sampling. The selected geese were slaughtered according to the Administration of Affairs Concerning Experimental Animals of Zhejiang Academy of Agricultural Science, and complete liver stripped from the carcass was weighed. Liver samples were collected and some were stored in −80°C for subsequent analyses.

The water content and fatty acid composition of liver segments collected from 10 geese per group at 90 days old were measured as previously described (Li et al., 2015; Liu et al., 2020). Liver water content and fatty acid composition in the liver segments were measured separately. Fresh liver samples were used as an analysis of water content by oven drying at 105°C for 24 h to a constant weight. A stored sample of liver was used to determine the fatty acid composition. Briefly, total liver lipids were extracted with petroleum ether/anhydrous diethyl ether (1:1, v/v). Fatty acid methyl esters (FAMEs) were obtained using saponification with a methanolic solution of potassium hydroxide. The organic layer was analyzed using an Agilent 7,890 N gas chromatograph equipped with a flame ionization detector (FID) and a CP-Sil 88 fused silica open-tube capillary column (100 m × 0.25 nm; Agilent Technologies, United States). The settings were: carrier gas: helium; column flow: 40 ml/min; temperature program: 140°C, 140–220°C, 25 min; 220°C, 40 min. The temperatures of the injector and detector were 240°C and 260°C, respectively. The fatty acid peaks were identified by comparing retention times of FAMEs from the sample with the retention time of the standards (Cat#: 18919-AMP; Sigma Chemicals, St, Louis, MO, United States). The amount of fatty acids was normalized to 100%. The extraction and measurement of total lipids refer to the method of Folch et al. (1957).

Sample preparation was performed according to the methods of Dunn et al. (2011). A total of 50 ± 1 mg of liver tissue was collected from each sample, placed into 2-ml EP tubes, spiked with 600 μl of extracting solution (methanol and chloroform, 3:1, v/v), and 10 μl of L-2-chlorophenylalanine (1 mg/ml stock in dH2O) and vortexed thoroughly for 30 s. All samples were ground using a grinding mill (Jingxin Industrial Development Co., LTD, Shanghai, China JXFSTPRP-24) for 4 min at 45 Hz and treated for 5 min (incubated in ice water) in an ultrasonic apparatus (Radbon Electronics Co. LTD, Shenzhen, China). After vortexing for 15 s and performing centrifugation at 12,000 rpm for 15 min at 4°C, the supernatant (500 μl) w a s transferred into 1.5-ml EP tubes. Meanwhile, 100 μl of supernatant from each sample was taken as a quality control (QC) sample. Subsequently, the supernatant was desiccated completely in a vacuum freeze drier (Huamei Biochemical Instrument Factory, Taicang, China). The dried samples reacted with 40 μl of methoxy amination hydrochloride (20 mg/ml in pyridine) and incubated for 30 min at 80°C, and 60 μl of bis (trimethylsilyl) trifluoroacetamide reagent (BSTFA, 1% TMCS, v/v; REGIS Technologies, United States) was added to each sample and incubated for 1.5 h at 70°C. When the mixture was cooled to room temperature, 5 μl of a standard mixture of fatty acid methyl esters (FAMEs, C8-C16: 1 mg/ml; C18-C24: 0.5 mg/ml in chloroform) was added to the sample. Afterward, all samples were analyzed by a gas chromatography system coupled with Pegasus HT time-of-flight mass spectrometry (GC-TOF/MS; Dunn et al., 2011). Supplementary Table S1 lists and provides references for derivation reagents used in this study.

Analysis refers to the methods of Yang et al. (2017) and Wen et al. (2018). GC-TOF/MS analysis of the extracted samples was performed using an Agilent 7,890 gas chromatograph (Agilent Corporation, Santa Clara, CA, United States) system coupled with a Pegasus HT time-of-flight mass spectrometry (GC-TOF/MS) system. According to the requirements of the test standard, the system was equipped with a DB-5MS capillary column (30 m × 250 μm × 0.25 μm, J & W Scientific, Folsom, CA, United States) coated with 5% diphenyl cross-linked with 95% dimethylpolysiloxane. In splitless mode with helium as the carrier gas, an aliquot (1 μl) of the sample was injected. Helium was used as the carrier gas, the front inlet purge flow was 3 ml/min, and the column flow was set at 1 ml/min. The initial oven temperature was set to 50°C, held for 1 min, raised to 310°C at a rate of 10°C/min, and maintained for 8 min at 310°C. Set the temperatures of the front injection (280°C), transfer line (280°C), and ion source (250°C), respectively. The energy was −70 eV in electron mode with a solvent delay time of 6.1 min. The mass spectrometry data were collected under the condition of in full-scan mode (m/z range: 50–500; rate: 12.5 spectra per second). Analysis of the QC sample was performed at regular intervals (every ten samples) to monitor the stability of the instrument.

Raw GC-TOF/MS data were first analyzed using Chroma TOF (version 4.3X, LECO) software for peak extraction, baseline correction, deconvolution, peak alignment, and peak integration. The GC-TOF/MS data were normalized by the area of the internal standards (IS). Relative quantification was performed according to peak area metabolite/peak area internal standard (Dunn et al., 2011). The LECO-Fiehn Rtx5 database was used for qualitative analysis (Kind et al., 2009). Peaks detected in <50% of QC samples or relative standard deviation >30% in QC samples were removed (Dunn et al., 2011). Center scaling (standardization) was performed prior to PCA. The data file was scaled with unit-variance scaling (standardization) before all the variables were subjected to OPLS-DA. Multivariate statistical analysis was performed to analyze the resulting normalization using SIMCA software (version 14.1, MKS Data Analytics Solutions, Umea, Sweden), including principal component analysis (PCA; the data were pre-processed using logarithmically transformed and Ctr scaling) and orthogonal projections to latent structure-discriminate analysis (OPLS-DA; Wiklund et al., 2008). The score plot of PCA shows the systematic clusters among the observations. High-level group separation can be obtained and the variables responsible for classification can be understood by supervised OPLS-DA. We identified the potential metabolites in the process of fatty liver formation combined with S-plot obtained from OPLS-DA analysis. In S-plot diagram, the farther away metabolite ions from origin represent the higher VIP value of the ions, and the higher VIP value represents the greater contribute to the difference between the two sample groups, we selected the ions with VIP >1 and further confirmed to the potential metabolite biomarkers of fatty liver formation. When the variable importance for the projection (VIP) values >1.0 in the OPLS-DA model and p < 0.05 in Student’s t-test, significantly different metabolites were identified. The heat map of data normalized by Z-score was generated using TBtools. Then, the Kyoto Encyclopedia of Genes and Genomes (KEGG)¹ was used to search for the related KEGG pathways of the metabolites. The high-quality KEGG metabolic pathway database MetaboAnalyst² was used for pathway analysis and visualization.

The experimental data measured included the water contents of the liver and fatty acid composition. This paper analyzed the data using one-way ANOVA in SPSS 19.0 software. Graphpad Prism 7 software was used for plotting the bar chart and data of fatty acid displayed using logarithmic coordinate axis. Tukey’s multiple comparison test was used to declare significant differences among groups. Significance was set at p < 0.05 for all statistical measures. This article displays results in mean ± SEM.

To examine the effect of overfeeding on the fatty acid composition in the livers of Landes geese, we determined the water content and fatty acid composition in the liver segments collected from 10 geese per group at 90 days old. As shown in Figure 2A, compared to the D0 group, the contents of moisture in the liver in the D7 and D25 groups were significantly reduced (p = 0.0067, p < 0.0001). Compared to 7 days of overfeeding, the contents of moisture in the liver were also significantly reduced on the 25th day of overfeeding (p < 0.0001). It was obvious that with the increase of the day for overfeeding, the contents of moisture in the liver showed a decreasing trend. To further probe the changes in liver metabolism, we examined the fatty acid composition in the livers of overfed Landes geese on different days. The results are shown in Figure 2B. The effects of overfeeding on total lipids in the liver of Landes geese are summarized in Supplementary Figure S1. The results indicated that the major fatty acids were oleic acid (C18:1n9), palmitic acid (C16:0), and stearic acid (C18:0). C18:1n9 was the most abundant fatty acid, accounting for 44.94 and 45.19% of the total fatty acids in the livers of Landes geese after 7 and 25 days of overfeeding, respectively. When compared with the D0 group, the levels of C20:2 (p = 0.0243), C20:4n6 (p = 0.0074), C20:5n3 (p < 0.0001), C22:1n9 (p = 0.0277), C22:3n6 (p = 0.0083), C22:5n3 (p = 0.0232), and C22:6n3 (p = 0.0359) increased in the livers of the D25 groups. In contrast, the level of C18:3n3 decreased (p = 0.004). The level of C20:4n6 (p = 0.0176) and C22:6n3 (p = 0.0360) were significantly increased after 7 days of overfeeding. When compared with the D7 group, C20:5n3 (p = 0.0035) levels were significantly higher in the D25 group. The difference was not significant among other groups.

The analysis of liver metabolism using GC-TOF/MS yielded 508 peaks after referring to the LECO-Fiehn Rtx5 database. At the same time, NIST library was used for identification. With one peak per compound. Data were p retreated with Chroma TOF 4.3X software to correct the data for missing values, eliminate noise, and normalizes the raw data signals was carried out by internal standard (IS) method, and then 422 peaks were preserved. Among these peaks, 198 compounds were relatively quantified, 224 uncharacterized peak were tagged with “unknown.” The D0, D7, and D25 groups of the total ion chromatograms (TIC) of goose liver samples are shown in Supplementary Figure S2A. Total ion current (TIC) plots of 5 QC samples were shown in Supplementary Figure S2B. The TIC peak retention time and peak area of QC samples overlapped well, indicating that it had good repeatability and instrumental stability for the metabolite detections.

To further explore the slight differences in metabolic profiles among all groups, we adopted multivariate statistical analysis methods to analyze metabolomics data. The PCA score plots showed the overall changes in metabolics under the effects of overfeeding. In the result of pairwise comparison (D7 vs. D0, D25 vs. D0 and D25 vs. D7), the values of R²X were 0.546, 0.515, and 0.504, respectively. The classification parameters of the three groups were R²X = 0.504. They are both >0.4. The results indicated that the D0, D7 and D25 groups showed a clear separation (Supplementary Figure S3). The supervised method, OPLS-DA, was used to screen the variables responsible for differences among these three groups. The OPLS-DA model displayed a clear separation and discrimination among groups (Figure 3A), which indicated that the metabolite profiles in the livers of geese apparently changed. The OPLS-DA model was tested by further permutation tests. The results showed that the, respective, R2Y and Q2 intercept values w e r e 0.944 and 0.57 in the model of the control group and D7 groups; 0.991 and 0.899 in the model of the D0 and D25 groups; and 0.976 and 0.909 in the model of the D7 and D25 groups (Figure 3B). The Q² intercept values were −0.56, −0.93 and −0.98, respectively. And all of the R²Y values were close to 1 in our tests, proving that OPLS-DA was established successfully. We identified the potential metabolites in the process of fatty liver formation combined with S-plot obtained from OPLS-DA analysis. In S-plot diagram the farther away metabolite ions from origin represent the higher VIP value of the ions, and the higher VIP value represent the greater contribute to the difference between the two sample groups, we selected the ions with VIP >1 and further confirmed to the potential metabolite biomarkers of fatty liver formation (Supplementary Figure S4).

To further probe into the variables in goose livers, the significantly different metabolites were filtered according to the significance test (p < 0.05) and VIP value (VIP > 1.0) from the OPLS-DA model. Of the metabolites shown to differ (VIP > 1.0 and p < 0.05) between the D7 and D0 groups, 30 were annotated and identified based on the metabolome databases of KEGG and 52 significantly different unidentified peaks (named “unknown”), 52 and 43 between the D25 and D0 groups, and 77 and 53 between the D25 and D7 groups (Supplementary Table S2). There were 95 significantly different metabolites between D25 and D0 (65 down-regulated, 30 up-regulated; Figure 4A), 82 between D7 and D0 (26 down-regulated, 56 up-regulated; Figure 4B), 130 between D25 and D7 (113 down-regulated, 17 up-regulated; Figure 4C). Table 1 shows this result. Subsequently, to further investigate the distinct characteristics of the significantly different metabolites, a hierarchical clustering heat map was plotted (Figure 5). The results from the comparison of the three groups indicated that overfeeding resulted in a continuous increase in oleic acid over time. This was the only significantly different metabolite that was continuously upregulated in the three groups. Although palmitic acid and pantothenic acid were found to have elevated concentrations, there was no difference between the D7 group and the D25 group. There were four metabolites (including palmitic acid, pantothenic acid, cholesterol, and 5′-methylthioadenosine) in the D25 group that exhibited higher concentrations in the D0 group and the D7 group. Meanwhile, overfeeding resulted in a continuous reduced effect in terms of the number of days for some significantly different metabolites. The metabolites included phenylalanine, methyl phosphate, sulfuric acid, and 3-hydroxybenzaldehyde. The other eight metabolites (including glutamic acid, O-phosphorylethanolamine, trans-4-hydroxy-L-proline 2,5,6-dihydrouracil, cytosin, 3-methylamino-1,2-propanediol, inosine 5′-monophosphate, and lactamide) had decreased concentrations after overfeeding for 7 and 25 days. However, compared to the D7 group, the D25 group did not have different concentrations. These 20 metabolites showed a downward trend with the increase in overfeeding days (D25 group vs. D7 group), mostly belonging to organic acids and derivatives (such as 1-aminocyclopropanecarboxylic acid, aminomalonic acid, allylmalonic acid, etc.).

To explore the significantly different metabolic pathways in the formation of fatty liver, the significantly different metabolites were imported into KEGG. When we compared the D7 to the D0, the D25 to the D0, and the D25 to the D7, 13, 29, and 38 metabolic pathways were obtained, respectively. They were then filtered based on the -ln value of p and pathway impact scores, and 6, 7, and 10 metabolic pathways were enriched, respectively (Figure 6). These key metabolic pathways include various amino acid metabolism (arginine and proline metabolism; glutathione metabolism; valine, leucine, and isoleucine biosynthesis; glycine, serine, and threonine metabolism; cysteine and methionine metabolism; alanine, aspartate, and glutamate metabolism; cyanoamino acid metabolism; taurine and hypotaurine metabolism; beta-alanine metabolism), lipid metabolism (glycerolipid metabolism; sphingolipid metabolism), nucleotide metabolism (purine metabolism; pyrimidine metabolism), pantothenate and CoA biosynthesis, sulfur metabolism, glyoxylate, and dicarboxylate metabolism and aminoacyl-tRNA biosynthesis pathways. Several distinct metabolites (including D-glyceric acid, L-cysteine, valine, glycine, glycine, sulfuric acid, uracil, isoleucine, glycerol, β-alanine, asparagine, aspartic acid, pantothenic acid, O-phosphorylethanolamine, inosine 5′-monophosphate, etc.) were involved in these metabolic pathways, as shown in Table 2.