3.1 Metabolite content
We detected 1153 features (740 in positive and 413 in negative ionisation mode, respectively). The abundance of each metabolite in each sample was calculated, and an accumulation column chart was used for visualisation so that compositional and structural differences between metabolites in different groups could be intuitively compared. Figure 1 shows the top 20 most abundant metabolites.
Some of the identified metabolites play specific roles in organisms, such as hormones and vitamins. We annotated all metabolites using KEGG database BR08001 to infer the biological roles played by the metabolites. The number of terms related to biological roles was then counted and a column chart was drawn (Fig. 2).
Metabolites detected by positive ion mode were divided into nine categories; control vs. treatment peptides were 46.27% vs. 58.55%, nucleic acids were 43.39% vs. 27.62%, vitamins and cofactors were 4.56% vs. 4.23%, and carbohydrates were 3.89% vs. 7.85%. Metabolites detected by negative ion mode were divided into eight categories; control vs. treatment organic acids were 50.45% vs. 49.76%, carbohydrates were 21.51% vs. 34.13%, nucleic acids were 9.14% vs. 5.40%, lipids were 8.14% vs. 6.32%, and peptides were 6.85% vs. 3.69%.
3.2 Principal component analysis
PCA was performed on the identified metabolites, a scatter diagram was drawn from the results, and structural differences between metabolites in different samples were inferred. Where the distributions of groups are clearly in different areas, this indicates that metabolite structures are quite different. In the PCA diagram (Fig. 3), each point corresponds to a sample, and the distance between two points is proportional to the differences of in metabolite structures between two groups (Euclidean distance). Different groups are marked in different colours, and the area marked by the ellipse is the 95% confidence region. In positive ion mode, the first principal component accounts for 67.30% of the variation, and the second principal component accounts for 8.10%. In negative ion mode, the first principal component accounts for 67.80%, and the second principal component accounts for 10.40%.
3.3 Screening and analysis of differential metabolites
In order to accurately identify differential metabolites, an OPLS-DA model was used to analyse the metabolomic data. The R2Y and Q2 values were close to 1 (in positive ion mode, R2X = 0.744, R2Y = 1, Q2 = 0.995; in negative ion mode, R2X = 0.776, R2Y = 0.999, Q2 = 0.994). This indicates that the OPLS-DA model was stable and highly reliable for this dataset, and could be confidently used to explore differences between treatment and control groups (Fig. 4).
In positive ion mode, 740 differentially abundant metabolites were identified, of which 401 were upregulated and 339 were downregulated; in negative ion mode, 413 differentially abundant metabolites were identified, of which 190 were upregulated and 223 were downregulated (Fig. 5).
3.4 Random forest analysis
In random forest analysis, mean decrease in accuracy and mean decrease in gini are used to measure the importance of discriminant groupings of metabolites. The value of a metabolite is changed to a random number, and the reduction in prediction accuracy of the random forest is the mean decrease accuracy; The mean decrease in gini is the effect of a metabolite on the heterogeneity of observed values at all nodes of the classification tree; the larger these two values are, the greater the importance of metabolites in the random forest. Figure 6 shows the 15 most important metabolites in the random forest, and these metabolites should be significantly different between groups.
In positive ion mode, pc (18:3e /22:6) is the most important metabolite. Compared with control groups, pc (18:3e /22:6), n-tetradecanamide, 9-oxo-ode, n1-(2-amino-2-oxoethyl)-2-phenoxyacetamide, 2-phenylglycine, pc (18:0/19:2) and isoproterenol were upregulated in treatment groups. Meanwhile, pantothenic acid, feruloyl putrescine, pc(16:2e/2:0), TNK, ILK, 3-(3,4-dihydroxyphenyl) propanoic acid, 1-[5-(2-phenyleth-1-ynyl)-2-thienyl]ethan-1-one oxime, and geranyl pp were less abundant. In negative ion mode, 5-(tert-butyl)-2-methyl-n-(5-methyl-3-isoxazolyl)-3-furamide was the most important metabolite. Compared with control groups, melatonin, 2-(formylamino)benzoic acid, undecanedioic acid, palmitic acid, n-acetyl-l-ornithine, d-glucuronic acid, sucrose, and 15(R), 19(R)-hydroxy prostaglandin F1α were upregulated in treatment groups. Meanwhile, 5-(tert-butyl)-2-methyl-n-(5-methyl-3-isoxazolyl)-3-furamide, mag (18:4), 12-hydroxydodecanoic acid, phloretin, sorbitan monostearate, 5'-deoxy-5'-(methylthio) adenosine and 5-s-cysteinyldopaquinone were downregulated.
3.5 Pathway enrichment analysis
The purpose of enrichment analysis is to search for biological pathways that play key roles in a biological process, and thereby reveal the basic molecular mechanisms of the biological processes. Prior to enrichment analysis, we selected certain metabolites, mainly those that differed significantly between groups (t-test, P< 0.05), and explored which metabolic pathways (KEGG species-specific metabolic pathways, slightly different metabolic pathways in different species) these metabolites are related to by calculating the ORA and p-values of the metabolic pathways, In this way, we determined whether the metabolites of concern (significantly differential metabolites) were significantly enriched in these metabolic pathways. Figure 7 shows the metabolic pathways with significant enrichment of differential metabolites. These metabolic pathways may be important in the biological processes studied.
According to the fold enrichment value ordered from higher to lower, in positive ion mode, differential metabolites were significantly enriched in 33 metabolic pathways, the top six of which were phenylalanine, tyrosine and tryptophan biosynthesis, biotin metabolism, phenylalanine metabolism, nitrogen metabolism, pyrimidine metabolism, and pantothenate and CoA biosynthesis. In negative ion mode, differential metabolites were significantly enriched in 52 metabolic pathways, the top six of which were alanine, aspartate and glutamate metabolism, lipoic acid metabolism, phenylalanine metabolism, cyanoamino acid metabolism cyanoamino, D-Arginine and D-ornithine metabolism, and pantothenate and CoA biosynthesis.
3.6 Topology analysis
The effects of metabolites on metabolic pathways can be estimated by topology analysis, and the effects can be measured by impact. For example, if there are no other metabolites or genes downstream of a metabolite in a metabolic pathway, then we can deduce that the effect of this metabolite on the metabolic pathway is negligable. Conversely, if a metabolite is very close upstream and there are many other metabolites and genes downstream, then the metabolite can be considered to have a large effect on the metabolic pathway.
Topological analysis is often combined with enrichment analysis to determine whether a metabolic pathway plays a key role in the biological process under study. The metabolic pathways in blue in Fig. 8 are prominent in the ORA results, and the vertical axis shows the impact of these metabolic pathways in the topological analysis. In positive ion mode, purine metabolism (map00230), phenylalanine metabolism (map00360), phenylalanine, tyrosine and tryptophan biosynthesis (map00400), pyrimidine metabolism (map00240) and nitrogen metabolism (map00910) play key roles. Alanine, aspartate and glutamate metabolism (map00250) and phenylalanine metabolism (map00360) play key roles in negative ion mode.
3.7 Metabolic pathway mapping
Metabolic pathway mapping can intuitively reflect upstream and downstream relationships, modes of action, and the metabolic pathway topological structure of metabolites, as well as identify genes associated with metabolites. Figure 9 shows metabolic pathways containing only metabolites, and red metabolites are the metabolites of concern (i.e. metabolites with significant differences between groups).