Transcriptome and Metabolome Analyses of the Flowers and Leaves of Chrysanthemum dichrum

Chrysanthemum dichrum is an important wild species in the family Asteraceae. However, because of a lack of genetic information, there has been relatively little research conducted on the molecular mechanisms in C. dichrum. There is no report describing the transcriptome and metabolome of C. dichrum flowers and leaves at different developmental stages. In this study, high-throughput sequencing and RNA-seq analyses were used to investigate the transcriptome of C. dichrum leaves, flower buds, and blooming flowers. Additionally, these three tissues also underwent a metabolomics analysis. A total of 447,313,764 clean reads were assembled into 77,683 unigenes, with an average length of 839 bp. Of the 44,204 annotated unigenes, 42,189, 28,531, 23,420, and 17,599 were annotated using the Nr, Swiss-Prot, KOG, and KEGG databases, respectively. Furthermore, 31,848 differentially expressed genes (DEGs) were detected between the leaves and flower buds, whereas 23,197 DEGs were detected between the leaves and blooming flowers, and 11,240 DEGs were detected between the flower buds and blooming flowers. Finally, a quantitative real-time Polymerase Chain Reaction (qRT-PCR) assay was conducted to validate the identified DEGs. The metabolome data revealed several abundant metabolites in C. dichrum leaves, flower buds, and blooming flowers, including raffinose, 1-kestose, asparagine, glutamine, and other medicinal compounds. The expression patterns of significant DEGs revealed by the transcriptome analysis as well as the data for the differentially abundant metabolites in three C. dichrum tissues provide important genetic and metabolic information relevant for future investigations of the molecular mechanisms in C. dichrum. Moreover, the results of this study may be useful for the molecular breeding, development, and application of C. dichrum resources.


INTRODUCTION
Chrysanthemum dichrum, which is a perennial herb belonging to the genus Chrysanthemum in the family Asteraceae, is a very important germplasm resource. It is a wild relative of the famous ornamental plant Chrysanthemum morifolium. A C. dichrum plant can be more than 30 cm tall. Its main stem, which lies prone or is slanted, is bare, brown, and has multiple branches on the upper part, which contains densely appressed pubescent leaves. It is a highly edible plant with medicinal value and is used for landscaping as a garden greening material. Chrysanthemum dichrum is distributed in southwestern Hebei province in China (Hsianshiu, 1993). However, it can be introduced to other regions relatively easily, making it a wild flower species potentially useful as a functional flower crop.
Current research on C. dichrum mainly focuses on cross breeding and stress resistance, including tolerance to cold and drought. Previous studies proved that C. dichrum, Chrysanthemum indicum "Nankingense, " and other interspecific hybrid progeny are highly tolerant to drought conditions (Li et al., 2020). Chi et al. in order to solve the economic benefit and ornamental quality decrease caused by low temperature in northern China. The genetic control mechanism of interspecific cold tolerance between C. dichrum and C. nankingense was clarified by analyzing the cold-resistant traits such as Low Semi-Lethal Temperature (LT50), basal shoots height and number in the basal shoots stage of the hybrids (Chi et al., 2018). A research team found a 1,474 bp stress-induced CdDREBa promoter in C. dichrum. Identified a stress-inducible CdDREBa promoter from C. dichrum, Clothes several candidate stressrelated cis-acting elements (MYC-box, MYB site, GT-1, And W-box) within it. CdDREBa promoter features a strong low temperature-and right-inducible promoter. This provides an attractive target for engineering inducible promoters in transgenic crops . Additionally, the CdICE1 gene in C. dichrum was cloned and an analysis of its cold tolerance-related function in Arabidopsis thaliana revealed that the miR398-CSD pathway is involved in the induction of freezing resistance (Chen Y. et al., 2013). Earlier research also confirmed that the structural expression of ICE1 gene in C. dichrum in large Chrysanthemum improves the tolerance to low temperature, salt and drought. CdICE1 represents a promising candidate for a biotechnological approach to improve the level of crop abiotic stress tolerance . Other studies applied RNA-seq technology to investigate chrysanthemum and related species (Chen et al., 2009;Wang et al., 2013;Xu et al., 2013;Ren et al., 2014;Hong et al., 2015;Liu et al., 2015), but the C. dichrum transcriptome and metabolome have been rarely examined. Accordingly, research regarding the C. dichrum transcriptome may increase our understanding of the gene regulatory mechanism and biological pathways in this wild species, while also identifying the downstream target genes of key transcription factors to clarify specific processes . Therefore, in this study, we analyzed the transcriptome and metabolome of C. dichrum leaves, flower buds, and blooming flowers. Moreover, we identified the genes differentially expressed in the examined plant tissues. We also determined the metabolite composition and content in the leaves, flower buds, and blooming flowers and identified some highly abundant economically valuable metabolites with substantial developmental potential. The findings of this study lay the foundation for future investigations regarding chrysanthemum resource development as well as molecular research and breeding related to C. dichrum.

Illumina Sequencing and Assembly
This study involved analyses of C. dichrum leaves, flower buds, and blooming flowers (Figure 1). The transcriptomes of these plant tissues were sequenced using the Illumina HiSeq TM 4000 platform. A total of 445,396,894 clean reads were obtained and 77,683 unigene sequences were assembled using de novo assembly technology. The N50 was 1,413 bp and the maximum and minimum lengths were 14,459 and 201 bp, respectively, with an average length of 839 bp ( Table 1).

Gene Annotation and Functional Classification
Of the 77,683 unigenes, 44,204 were annotated. More specifically, 42,189 unigenes were annotated using the Nr database, 28,531 unigenes were annotated using the Swiss-Prot database, 23,420 unigenes were annotated using the KOG database, and 17,599 unigenes were annotated using the KEGG database. A total of 33,479 unigenes were not been annotated (Figure 2).
Regarding the GO analysis, 55,748 unigenes were classified into 50 groups in the three main categories. The biological process category had the most unigenes (27,682). "Metabolic process" and "Cellular process" were the main GO terms in this category, with 6,476 and 5,866 unigenes, respectively. A total of 16,466 unigenes were classified in the cell component category, with most annotated with the "Cells" (3,756 unigenes) and "Cell parts" (3,754 unigenes) GO terms. The molecular function category included 11,599 annotated unigenes, with "Catalytic activity" (6,303 unigenes) and "Binding" (4,296 unigenes) representing the main terms. Additionally, of the GO terms in the biological process, cell component, and molecular function categories, "Locomotion, " "Nucleoid, " and "Translation regulator activity" were assigned to the fewest genes, respectively (Figure 4).

Analysis of Differentially Expressed Genes in Chrysanthemum dichrum
We used a false discovery rate (FDR) < 0.05 and a |log 2 (foldchange)| > 1 as the thresholds for identifying differentially expressed genes (DEGs). A total of 31,848 genes were differentially expressed between the C. dichrum leaves and flower buds. Compared with the leaf expression levels, 18,718 and 13,130 DEGs had significantly up-regulated and downregulated expression levels in the flower buds. A total of 23,197 DEGs were identified between the leaves and the blooming flowers. Compared with the corresponding expression during the blooming period, 13,640 and 9,557 DEGs had significantly up-regulated and down-regulated expression levels in the leaves. A comparison between the flower buds and blooming flowers revealed 11,240 DEGs, of which 4,436 and 6,804 DEGs had significantly up-regulated and down-regulated expression levels in the flower buds ( Figure 5). The DEGs underwent GO and KEGG pathway enrichment analyses to determine the differences in biological processes and pathways among the leaves, flower buds, and blooming flowers. The DEGs between the leaves and flower buds were annotated with 47 GO terms. In the biological process category, the main GO terms were "Metabolic process" (GO:0008152), "Cellular process" (GO:0009987), and "Single-organism process" (GO:0044699), which were assigned to 3,146, 2,729, and 940 DEGs, respectively. In the cellular component category, the most represented GO terms were "Cell" (GO:0005623), "Cell part" (GO:0044464), and "Organelle" (GO:0043226), which were assigned to 1,706, 1,706, and 1,056 DEGs, respectively. The most common molecular function GO terms among the DEGs were "Catalytic activity" (GO:0003824), "Binding" (GO:0005488), and "Transporter activity" (GO:0005215), which were assigned to 3,142, 2,008, and 284 DEGs, respectively. A total of 134 KEGG pathways were enriched among 4,253 DEGs. The main enriched KEGG pathways were "Metabolic pathways" (ko01100), "Biosynthesis of secondary metabolites" (ko01110), and "Ribosome" (ko03010), which were associated with 1,853, 1,063, and 364 DEGs, respectively ( Figure 6A).

Multivariate Statistical Analysis of the Metabolome
A principal component analysis (PCA) can generally reflect the overall metabolic differences between samples in each group and the degree of variability between samples within a group. The first and second principal components (PC1 and PC2, respectively) of the sample group data for the C. dichrum leaves, flower buds, and blooming flowers accounted for 54.8 and 28.4% of the total variability, respectively (Figure 7). These results reflect relatively large overall metabolic differences, but low diversity between samples.
Qualitative and quantitative data were obtained for 568 metabolites in C. dichrum leaves, flower buds, and blooming flowers. Additionally, 216 annotated substances were highly abundant in the three plant tissues, including m-cresol, phytosphingosine, tartaric acid, abietic acid, neohesperidin, oxalic acid, salicylic acid, xylitol, glutamine, asparagine, 1kestose, and other heterochromia-related metabolites. These compounds have diverse uses in medicines, skin care products, and food additives, indicative of the importance and utility of C. dichrum.
We screened 48 metabolites with significant differences between the C. dichrum leaves and flower buds. Of these metabolites, 22 and 26 were, respectively, more and less abundant in flower buds than in leaves. A total of 49 metabolites differed significantly between the leaves and blooming flowers, of which 28 and 21 metabolites were respectively more and less abundant in blooming flowers than in leaves. A comparison between the flower buds and blooming flowers revealed 22 significantly different metabolites, with all but one more abundant in the blooming flowers than in the flower buds (Figure 8).

KEGG Annotation and KO Enrichment Analysis
On the basis of the KEGG pathway enrichment analysis, 296 metabolites were assigned to 17 metabolic pathways. The most represented pathway was the "Global and overview " (51 metabolites), followed by "Carbohydrate metabolism" (32 metabolites) and "Amino acid metabolism" (18 metabolites). The least annotated metabolic pathway was "Metabolism of terpenoids and polyketides" (1 metabolite) (Figure 10).
The four KO pathways with the most differentially abundant metabolites and the smallest P/Q value between the leaves and flower buds were arginine biosynthesis, glycerolipid metabolism, citrate cycle (TCA cycle), and carbon metabolism ( Figure 11A). Additionally, the four KO pathways with the most differentially abundant metabolites and the smallest P/Q value between the leaves and blooming flowers were glycerolipid metabolism, glyoxylate and dicarboxylate metabolism, ABC transporters, and carbon metabolism ( Figure 11B). Furthermore, the four KO pathways with the most differentially abundant metabolites and the smallest P/Q value between the flower buds and blooming flowers were fatty acid elongation; cutin, suberine, and wax biosynthesis; glutathione metabolism; and fatty acid metabolism ( Figure 11C).

Verification of Gene Expression Profiles Using qRT-PCR
To further verify the expression profiles of genes in the Illumina sequencing analyses, 15 ungenes with high expression levels and large differential multiple were selected for qRT-PCR, and leaf (YS1), flower bud (YS2) and blooming flower (YS3) were selected for RNA-seq. According to transcriptome sequencing data, the expression level of unigenes in the former group was significantly higher than that in the latter group in the three control groups. The RT-PCR results showed that the expression patterns of these 15 genes were consistent with the sequencing data (Additional file 10, Figure 12).

Common Metabolic Pathway Analysis
Due to the few types of differential metabolites among sample groups, and the absence of target metabolites in some differential metabolites, the common pathway analysis of all genes and all metabolites was carried out.
Combined analysis of transcriptome and metabolome sequencing showed that common pathway analysis was conducted for all genes and all metabolites. Due to the small number of differential metabolites among sample groups and the lack of target metabolites in some differential metabolites, we conducted common pathway analysis for all genes and metabolites.
There were 64 metabolic pathways shared by all genes and metabolites, 9016 candidate genes with pathway annotations, and 64 metabolites with pathway annotations. A total of 3,485 candidate genes (38.65%) and 42 metabolites (65.63%) were annotated for Metabolic pathways. A total of 1,917 candidate genes were annotated by Biosynthesis of secondary metabolites pathway, accounting for 21.26% of the total, and 18 metabolites were annotated, accounting for 28.13% of the total. A total of 570 candidate genes were annotated by Carbon Metabolism Pathway, accounting for 6.32% of the total, and 10 metabolites were annotated, accounting for 15.63% of the total. These three paths are the most annotated (Additional File 14).

Correlation Coefficient Model Analysis of Transcriptome and Metabolome
The correlation between genes and metabolites was evaluated according to Pearson's correlation coefficient. The top 250 differential genes and metabolites with the absolute correlation coefficient greater than 0.5 were screened out, and the correlation network diagram of gene expression and metabolite abundance was drawn (Figure 13). The results showed that there were more positive correlations between different genes and different metabolites, and only 12 groups showed negative correlations.

DISCUSSION
On the basis of a metabolomic analysis of C. dichrum tissue samples from three different developmental periods, we identified several differentially abundant functional metabolites among the 568 total metabolites. Two metabolites associated with food additives, raffinose and 1-kestose, were detected in leaves.
Raffinose, which is also called melitriose, is a trisaccharide comprising galactose, fructose, and glucose. Melitriose is a functional oligosaccharide that can induce substantial bifidobacterial proliferation. It can simultaneously promote the reproduction and growth of beneficial bacteria, such as bifidobacteria and lactobacilli, while also inhibiting the reproduction of harmful bacteria in the intestines, ultimately helping to establish a healthy intestinal microbial flora (Cox, 1966;Mussatto and Mancilha, 2007). Previous studies confirmed that melitriose can also prevent constipation and inhibit diarrhea; function as a two-way regulator to detoxify and protect the liver; inhibit the production of toxins in the body; minimize the burden on the liver; regulate the human immune system; enhance immunity; decrease allergic reactions following an oral administration (Zhang et al., 2004(Zhang et al., , 2008Zhang R. et al., 2013). It also can improve lipid metabolism, lower blood lipid and cholesterol levels, and prevent dental caries; however, it is not effective against oral cariogenic bacteria. Additionally, even in the presence of sucrose, melitriose can decrease tartar formation and protect regions where oral microorganisms are deposited and produce acid, thereby preventing corrosion. It can also whiten and strengthen teeth. Because melitriose is a low-calorie trisaccharide, it does not affect the blood sugar level in the human body and can be consumed by people with diabetes (Yazawa and Tamura, 1982;Muthukumaran et al., 2018).
As the smallest oligofructose (i.e., fructooligosaccharide) compound, 1-kestose is an active ingredient naturally found in diverse foods, including fruits, vegetables, and honey, and it is an excellent water-soluble dietary fiber. Fructooligosaccharides, which comprise 1-kestose (GF2), nystose (GF3), and 1Ffructofruranosy 1 nystose (GF4) (Jayalakshmi et al., 2021), can activate bifidobacteria, while also regulating the intestinal microbial flora balance, ensuring the bowel is appropriately hydrated, regulating blood lipid levels, lowering cholesterol and blood sugar levels, and enhancing immunity (Li et al., 2004;Probert et al., 2004;Mabel et al., 2008;Hidaka et al., 2010;Delgado et al., 2012). Owing to its superior physiological functions, oligofructose has become a widely popular functional food component in the international food market during the past 10 years. It has been used in more than 500 food and health products and medicines. Moreover, it has been described as a healthy sugar for the 21st century (Yun, 1996).
Metabolites (e.g., glutamine) in edible medicinal plants have been detected in C. dichrum flower buds. Basic and clinical experiments demonstrated that glutamine can decrease catabolic activities, promote protein synthesis, improve immune functions, protect the intestinal mucosal barrier, and accelerate wound healing. It can also be used to improve gastrointestinal health by preventing gastric and duodenal ulcers, gastritis, and hyperacidity (Taxuke et al., 2003). Asparagine is a medically relevant metabolite that can lower blood pressure, expand bronchial tubes (useful for treating asthmatics), and prevent peptic ulcers and gastric dysfunction. It can also be used for cultivating microorganisms and treating acrylonitrile wastewater.
In this study, raffinose and asparagine contents were higher in C. dichrum blooming flowers than in leaves. Additionally, another metabolite, xylose, was also highly abundant in blooming flowers. This sugar has been used as a food additive. As a component of xylan, it is widely distributed in plants. Xylose is not digested and absorbed, making it useful as a sweetener that will not lead to weight gain. It can also activate and promote the growth of bifidobacteria in the human intestine, with beneficial effects on health. Ingesting foods containing xylose can improve the microbial environment of the human body and enhance immunity. It cannot be used by microorganisms in the oral cavity. Additionally, it has some physiological functions of dietary fiber and it can decrease serum cholesterol levels and prevent colon cancer. Hence, adding a small amount of xylose to food can lead to increased positive health effects (Liang and Zhao, 2011). Xylose is widely present in plant hemicellulose as macromolecular xylan and can be obtained by degrading xylan with acids or enzymes (Harner et al., 2015). Xylooligosaccharides are an important component of functional foods. Because they are stable and nontoxic sugars that can be efficiently purified, xylooligosaccharides are increasingly being used in food, feed, and medicine, while also being applied in other fields (Zhang Y.H. et al., 2013;Baker et al., 2021).
In the current study, the 1-kestose and xylose contents were high in blooming flowers. Similarly, many other functional metabolites were more abundant in blooming flowers than in flower buds and leaves, including tartaric acid and salicylic acid. Tartaric acid is an antioxidant that has been added to foods. It is most commonly used as a beverage additive and a raw material for the pharmaceutical industry (Kroschwitz and Seidel, 2004). Earlier research proved that salicylic acid is useful for treating acne and for lightening post-acne pigmentation, minimizing pore size, removing fine wrinkles, and preventing sun-induced skin aging. The renewal of skin cells in response to salicylic acid will result in smoother skin (Kligman and Kligman, 1998). In small amounts, salicylic acid is also applicable as a food preservative (Chrubasik et al., 2000). It is also an important raw material for the pharmaceutical industry.
The data presented herein indicate that C. dichrum flowers and leaves are rich in metabolites with important functions. Consequently, C. dichrum is an important plant source of useful compounds. In terms of the extraction of functional metabolites from C. dichrum, the flowers at full bloom are likely the most appropriate plant tissues because they contain large amounts of diverse metabolites.

Plant Materials and RNA Extraction
The C. dichrum plants used in this study were cultivated at 23 • C in a greenhouse at the Beijing Academy of Agriculture and Forestry Sciences (116.3 • E, 39.9 • N) with an 8-h light/16-h dark cycle. Leaves, flower buds, and blooming flowers were collected from C. dichrum plants, with three biological replicates per sample. This plant samples was used for both transcriptome analysis and Metabolomic    Pathways with the smallest P value (or Q value) were used for mapping. The ordinate was Pathway, and the abscissa was enrichment factor (the number of differences in this pathway divided by all numbers). The size represented the number, and the redder the color, the smaller the P/Q value.
analysis. The collected samples were quickly frozen in liquid nitrogen and stored at −80 • C. The RNeasy Plant Mini Kit (Qiagen, China) was used to extract total RNA from the frozen samples. The concentration of the RNA was determined using the NanoDrop ND2000 spectrophotometer (Thermo Scientific).

Library Construction and Sequencing
To construct a sequencing library, rRNA was eliminated from the total RNA extracted from the leaves, flower buds, and blooming flowers using the Ribo-Zero TM Magnetic Kit (Epicenter), after which the mRNA was enriched using oligo-(dT) beads. Fragmentation buffer was used to produce short mRNA fragments, which were then reverse transcribed into cDNA. The second cDNA strand was synthesized in buffer containing DNA polymerase I, RNase H, and dNTP. The cDNA fragments were purified using the QiaQuick PCR extraction kit. After repairing the ends and adding a poly(A) tag, an Illumina sequencing adapter was ligated. The size of the ligation products was determined by Gene de novo using the Illumina HiSeq TM 4000 system (Illumina, San Diego, CA, United States) for the subsequent agarose gel electrophoresis, PCR amplification, and sequencing. The PacBio Sequel system (PacBio, CA, United States) was used for sequencing. To improve the accuracy of the PacBio data, the reads were filtered by deleting reads with adapters, reads with more than 10% unknown nucleotides (N), and reads with more than 40% low-quality nucleotides (Q value ≤ 20). The CD-HIT (version 4.6.7) software (with a sequence identity threshold of 0.99) was used to remove redundant sequences to obtain the final unigene sequences.

Basic Annotation of Unigenes
Unigenes were annotated using the BLASTX algorithm and the following four databases (E value of 1.00E-5): Nr (NCBI), Swiss-Prot 1 , KEGG 2 , and KOG 3 . The sequence direction of each gene was determined according to the best alignments. Unigenes were annotated with GO terms using the Blast2GO program (Conesa et al., 2005). On the basis of the Blast2GO analysis, we selected high-quality unigenes with a high hit rate. These unigenes were functionally classified using the WEGO software (Ye, 2006

Analysis of Chrysanthemum dichrum Transcriptome Sequencing Results
The RPKM values derived from the RNA-seq data were used to represent unigene expression levels (Fourquin et al., 2005). The DEGs in the chrysanthemum transcriptome were screened and analyzed as previously described (Samarskiǐ and Claverie, 1997). To verify the P-value threshold, the FDR was used in multiple tests and analyses (Endress, 2006). The threshold for determining a significant difference in gene expression was set as follows: absolute value of the log 2 (ratio) ≥ 2 and FDR < 0.05 (Mortazavi et al., 2008). Differential gene expression was analyzed for the DEGs with expression levels that differed by at least 2-times between samples.

Alternative Splicing Detection
To analyze alternative splicing events in a unigene, transcripts were divided into gene families according to k-mer similarity. This process uses the coding genome reconstruction tool (Cogent) and then recombines each family into a coding reference genome. This process uses the De Bruijn diagram (Li et al., 2017). The SUPPA tool was used to analyze the alternative splicing events of a unigene (Alamancos et al., 2015).

Gene Expression Analysis Based on qRT-PCR
To conduct qRT-PCR assays, the total RNA extracted from the leaves, flower buds, and blooming flowers were treated with DNase (Promega, United States) to eliminate any residual genomic DNA. The purified RNA served as the template for synthesizing cDNA using a commercial reverse transcription kit (Tsingke, China). The qRT-PCR analysis was completed using the PikoReal system (Thermo Fisher Scientific, Germany). The 20-µL reaction solution comprised 1 µL reverse transcribed cDNA as the template. The PCR program was as follows: 95 • C for 30 s; 40 cycles of 95 • C for 5 s and 60 • C for 30 s. Details regarding the qRT-PCR primers used to determine the relative expression levels of specific genes are listed in Additional file 13. Each sample was analyzed in triplicate, with three biological replicates. The 2 − Ct method was used to calculate relative gene expression levels. The reference control was the Aspergillus gene encoding protein phosphatase 2A (PP2Ac) (Xue et al., 2014).

Chemicals and Reagents
The chemicals and reagents used in this study, such as methanol, ethanol, and acetonitrile, were purchased from Merck (Germany; 4 ). The Milli-Q system (Millipore, Bedford, MA, United States) was used to produce ultrapure water. Authentic standards were purchased from BioBioPha Co., Ltd. 5 and Sigma-Aldrich 6 . All chemicals and reagents were of analytical grade.

Sample Preparation and Extraction
Freeze-dried samples were ground to a powder using zirconia beads and a mixing mill (MM 400, Retsch) set at 30 Hz for 1.5 min. The powdered material was weighed, after which 100 mg was mixed with 1.0 mL 70% methanol aqueous solution (containing 0.1 mg/L lidocaine as an internal standard) overnight at 4 • C. After centrifuging the mixture at 10,000 g for 10 min, the supernatant was collected and filtered (SCAA-104; pore size 0.22 µm; ANPEL, Shanghai, China; 7 ) for the subsequent LC-MS/MS analysis. To assess the repeatability of the entire experiment, quality control samples were mixed with all samples.

AB Sciex QTRAP6500 (UPLC) Analysis
The extracted compounds were analyzed using an LC-ESI-MS/MS system [UPLC: Shim-pack UFLC SHIMADZU CBM30A 8 ; MS/MS: Applied Biosystems 6500 QTRAP 9 ] 7 www.anpel.com.cn/ 8 http://www.shimadzu.com.cn/ 9 http://www.appliedbiosystems.com.cn/ (Chen W. et al., 2013;Wishart et al., 2013). For each sample, a 2-µL aliquot was injected into the Waters ACQUITY UPLC HSS T3 C18 column (2.1 mm × 100 mm, 1.8 µm) operating at 40 • C with a flow rate of 0.4 mL/min. The mobile phases used were acidified water (0.04 % acetic acid) (Phase A) and acidified acetonitrile (0.04 % acetic acid) (Phase B). Compounds were separated using the following gradient: 95:5 Phase A/Phase B at 0 min; 5:95 Phase A/Phase B at 11.0 min; 5:95 Phase A/Phase B at 12.0 min; 95:5 Phase A/Phase B at 12.1 min; and 95:5 Phase A/Phase B at 15.0 min. The eluent was analyzed by an ESI-triple quadrupole-linear ion trap (QTRAP) mass spectrometer. The LIT and triple quadrupole (QQQ) scans were acquired using a triple quadrupole-linear ion trap mass spectrometer (QTRAP; AB Sciex QTRAP6500 System) equipped with an ESI-Turbo Ion-Spray interface. The system was operated in the positive ion mode and controlled with the Analyst 1.6.1 software (AB Sciex). The operating parameters were as follows: ESI source temperature, 500 • C; ion spray voltage, 5,500 V; curtain gas, 25 psi; and collision-activated dissociation, highest setting. The QQQ scans were acquired as MRM experiments with optimized declustering potential and collision energy for each MRM transition. The m/z range was set between 50 and 1,000.

Data Preprocessing and Metabolite Identification
To generate a matrix containing relatively little offset and redundant data, the peak signal/noise was manually checked (>10) and an internal Perl software was used to remove redundant signals resulting from different isotopes and to analyze intra-source fragmentation and K + , Na + , and NH 4 + adducts and dimers. The accurate m/z of each Q1 was obtained, after which the area of each chromatographic peak was calculated. The peaks for the different samples were aligned according to the spectrum and retention time. Metabolites were identified by searching internal and public databases [MassBank, KNApSAcK, HMDB (Zhu et al., 2013), MoTo DB, and METLIN (Worley and Powers, 2013)]. The m/z values, RT, and lysis patterns were compared with standards. The data filtering as well as the peak detection, comparison, and calculation were performed using the Analyst 1.6.1 software.

Multivariate Statistical Analysis
For all samples, the R package model 10 was used to perform the unsupervised dimensionality reduction PCA to initially visualize the differences between different sample groups (Kanehisa et al., 2008).

Differentially Abundant Metabolite Analysis
The T test and the variable importance in the projection (VIP) score of the (O)PLS model were used to rank the metabolites with significant differences in abundance between two groups. The P-value threshold for the T test was set at < 0.05 and the VIP score threshold was set at ≥ 1.

KEGG Pathway Analysis
The enriched KEGG metabolic pathways among the metabolites were determined. The metabolic and signal transduction pathways significantly enriched among the differentially abundant metabolites were identified using the following formula: where N is the number of all metabolites annotated using the KEGG database, n is the number of different metabolites in N, M is the number of all metabolites assigned to a specific pathway, and m is the number of different metabolites in M. The P-value was adjusted using an FDR ≤ 0.05. The pathways satisfying this condition were identified as significantly enriched among the differentially abundant metabolites.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in NCBI, accession number PRJNA744998.