To study the differences in metabolites between wolfberry (Lycium, of the family Solanaceae) and other species, we selected three monocot crops and three dicot crops as reference objects. The monocot crops included three major staple crops maize (Zea mays L.), rice (Oryza sativa L.) and wheat (Triticum aestivum L.), with rice being the model plant for monocots. The dicot crops included legume crop soybean (Glycine max (L.) Merr.) and the fruit crops tomato (Solanum lycopersicum L.) and grape (Vitis vinifera L.). Wolfberry (L.barbarum ‘Ningqi No.1’, L. barbarum var. auranticarpum, L. ruthenicum) was obtained from Yinchuan, Ningxia (E113°42′, N34°48′). Maize (Waxy maize, Red Waxy maize, Fruit maize, Sweet maize), Grape (Red grape, Green grape, Seedless red grape, Jufeng grape), Tomato (Pink tomato, Tomato, Cherry tomato, Millennium cherry tomato) were obtained from Haikou, Hainan (E110°20′, N20°02′). Soybeans (GDC058, GDC062, GDC063) were obtained from Zhengzhou, Henan (E113°42′, N34°48′). Wheat (Lumai 21, Heng 7228, Linmai No. 2, Xinong 529) was obtained from Zhaoxian, Hebei (E114°28′, N38°02′). Rice (Huanghuazhan, Nipponbare, Mimghui 63, Zhengshan 97, Zhonghua 11) was obtained from Hainan University (E110°20′, N20°02′).

The samples were put in a lyophilizer for vacuum freeze-drying, and they were crushed in a mixer mill (MM 400; Retsch, Haan, Germany) for 1 min at 30 Hz. Next, 80 mg of the powdered sample was weighed into a 2 mL centrifuge tube, and 70% aqueous methanol (v/v) with lidocaine internal standard was added to extract the water-soluble metabolites (pure methanol was used to extract the fat-soluble metabolites). Then, the tube was vortexed for 10 seconds, allowed to stand for 10 min, repeated three times, and placed in a refrigerator at 4°C for 10-12 hours. Then, the sample was centrifuged (4°C, 10000 rpm, 10 min), the supernatant was pipetted, the water-soluble and fat-soluble metabolites were mixed 1:1 and filtered by a microporous filtration membrane (SCAA-104, 13 mm, 0.22 μm, Shanghai Anpu Experimental Technology Co., Ltd., http://www.anpel.com.cn/), and the sample was filtered into an injection bottle for storage for UPLC-MS analysis.

The instruments used for LC-MS/MS analysis included UPLC-Q Exactive Plus Orbitrap HRMS and UPLC-Q-Trap 6500+ MS. The analytical column used was a C18 column (Shim-pack GLSS C18, 1.9 μm, 2.1*100, Shimadzu).

UPLC chromatographic conditions: mobile phase A is an aqueous solution containing 0.04% glacial acetic acid, and mobile phase B is a methanol solution containing 0.04% glacial acetic acid. Elution gradient: At 0 min, V_{phase A}: V_{phase B} = 95:5; At 10 min, V_{phase A}: V_{phase B} = 5:95; At 11 min, V_{phase A}: V_{phase B} = 5:95; At 11.1 min, V_{phase A}: V_{phase B} = 95:5, this ratio was continued until the end (duration is 14 minutes). The column temperature was set to 40°C, the injection volume of the injector was 2 μL, and the flow rate was 0.35 mL/min.

Orbitrap HRMS mass spectrometry conditions: ESI ionization method, mass spectrometry scanning mode is Full MS/ddMS2, ion collection mass range is 100~1200 m/z, and the lysis voltage is set to 20, 40, 60 eV; The spray voltage in positive ion mode is 3.5 kV, the capillary temperature is 350°C, the heater temperature is 350°C, the sheath gas (nitrogen) flow rate is 40 arb, and the auxiliary gas (nitrogen) flow rate is 10 arb. The spray voltage in negative ion mode is -3.0 kV, the capillary temperature is 350°C, the heater temperature is 350°C, the sheath flow rate is 30 arb, and the auxiliary air flow rate is 10 arb.

For wide-target detection of metabolites by UPLC-Q-Trap 6500+ multiplex reaction monitoring mode, the MRM detection window was set to 60 s, and the target cycle time was set to 0.8 s. The raw data were integrated by Multi Quant 3.0.3 to accurately obtain the relative content of each substance.

All statistical analyses were carried out using R (4.1.1, http://www.r-project.org). (PCA) was performed using the R package “FactoMineR” (Lê et al., 2008) with 587 metabolite data from 3 biological replicates of mixed samples of each species to evaluate the metabolome differences across seven species. Circular plots were constructed using the R package Circlize, with the raw data normalized and scaled in the R program; the y-axes of the seven circles were set on the same scale. The hierarchical clustering tree using metabolome data of the seven species was constructed using the R package “hclust”. Orthogonal partial least squares discriminant analysis (OPLS-DA) was conducted by the R package “ropls” (Thévenot et al., 2015) to identify the major discriminant metabolite features among different species. A metabolite feature was considered a species-specific metabolic trait when matching the following criteria compared with the other six species: (i) P value of paired t-test ≤ 0.05; (ii) Fold Change ≥ 3. Venn diagrams were generated on the online website jvenn (http://jvenn.toulouse.inra.fr/app/index.html).

Heatmaps were generated using the R package “pheatmap” with data normalization to divide the metabolome into eight clusters, and comparisons were performed via an average linkage method based on the Manhattan distance. Metabolic traits of different clusters were used for the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to elucidate the differential metabolic pathways among these seven different species. Enrichment analysis was carried out using the online platform MetaboAnalyst (https://www.metaboanalyst.ca/) selecting Oryza sativa japonica (Japanese rice) (KEGG) as a reference database. The bubble plot was drawn by the R package “ggplot2” (Ginestet, 2011).

To dissect the metabolome differences between wolfberry and diverse species, six popular species were selected for comparison with wolfberry for analysis in this study. These species include four cereal crops maize, rice, wheat, legume crop soybean and two fruit crops tomato and grape. To reduce the specific error associated with a variety, we collected commercially available varieties of these seven species and mixed them to prepare mixed samples that represented the metabolome of each species. Collectively, a total of 1043 distinct metabolic features were detected and quantified in these seven species, the identification class is shown in Supplementary Figure S1A . Based on the comparison of our local metabolite database and commercial standards (Chen et al., 2014; Li et al., 2022b), 587 metabolites were identified, which could be divided into 10 categories, including 127 phenylpropanoids, 125 amino acids and their derivatives, 116 organic acids and sugars, 66 lipids, 40 nucleic acids and their derivatives, 35 vitamins and cofactor derivatives, 20 terpenoids, 15 alkaloids, 10 phytohormones and 33 other metabolites ( Figure 1A , Supplementary Table S1 ).

Principal component analysis (PCA) was performed to assess the overall metabolome differences under the unsupervised model. The first principal component (PC1) explained the greatest variance (27.7%) of the variation among the seven species. Among the top 50 contributors of PC1, the proportion of amino acid and their derivative was the highest, followed by organic acids and sugars ( Supplementary Table S2 ). The second principal component (PC2) was orthogonal to PC1 and was the largest in the remaining variance, accounting for 19.4% of the variance. Among the top 50 contributors in PC2, phenylpropanoids had the highest proportion, followed by amino acid and their derivative ( Supplementary Table S2 ). For example, ornithine hydrochloride (ms025), the first-ranked identified metabolite, was found to accumulate more in the species with negative PC2 (grape, wolfberry and rice) than in those with positive PC2 (maize, wheat, soybean and tomato) ( Supplementary Figure S1B ). The results showed that the seven species were clearly separated, and three biological replicates of each species were compactly grouped together. This discrete clustering of the respective species indicates the distinct attributes of each species and the high repeatability and reliability of the experimental results ( Figure 1B ).

To visualize the variation in the metabolome among different species, we generated a circular plot for the seven species ( Figure 1C ). In the circular plot, strong variation in metabolic accumulation across seven species could be observed, with phenylpropanoids, nucleic acids and derivatives as well as phytohormones showing the most pronounced variation. To quantify the degree of variation, we used the coefficient of variation (CV) values of metabolites across the seven species. Phenylpropanoids were found to be the most variable metabolites, with a CV range from 11.56% to 244.93%, which is consistent with the results shown in the circular plot ( Supplementary Figure S1C , Supplementary Table S3 ). In contrast, nucleic acids and derivatives have the least variation, probably because they are crucial components of genetic material in all species, some of which play vital roles in influencing the structure and function of RNA or in post-transcriptional gene regulation (Dominissini et al., 2016). Therefore, the conservation of these compounds, which represented the minimum variation in our data, implied their importance in living organisms, illustrating the role of nucleic acids and derivatives as building blocks in the life of organisms.

The overall profile of all identified metabolites in the seven species was analyzed by hierarchical cluster analysis (HCA). As shown in the heatmap, the metabolic diversity in different species was further indicated. Furthermore, all the compounds were hierarchically clustered into eight main clades, the first seven of which specifically accumulated in wolfberry, rice, maize, wheat, tomato, grape and soybean, indicating that each species has its preferential metabolites ( Figure 2A , Supplementary Table S4 ).

To explore the preferred metabolic pathway of each species, we performed pathway enrichment analysis of the metabolites in clade 1 to clade 7 by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database ( Figure 2B , Supplementary Figure S2 ). Here, we focused on the preferred metabolite of wolfberry in clade 1. The 43 substances in clade 1 showed specific accumulation in wolfberry and were classified into 9 different categories. The two largest metabolites among them are organic acids and sugars and amino acids and their derivatives, accounting for 25.58% and 23.26%, respectively. Organic acids and sugars in wolfberry fruit are not only important nutrients but also the main influencing factors of flavor quality, while amino acids and derivatives are the main nutritional and medicinal components in LBFs (Zhao et al., 2020a). In the KEGG pathway enrichment analysis, the 43 metabolites were involved in 19 pathways. The major pathways are presented in the bubble plot ( Figure 2B ). In addition to amino acid metabolism, they also include glyoxylate and dicarboxylate metabolism, thiamine metabolism, aminoacyl-tRNA metabolism, purine metabolism, cyanoamino acid metabolism and anthocyanin biosynthesis. Similarly, 48 metabolites in clade 6 showed more accumulation in grapes and were enriched in the glutathione metabolic pathway. This could explain the antioxidant activity conferred by glutathione in grapes ( Supplementary Figure S2F ). Collectively, these results suggest that the specific accumulation of metabolites in different species can represent and determine their specific nutritional value.

To explore the affinities among the seven species, a neighbor-joining tree was constructed using metabolome data of the seven species ( Figure 3A ), and we also created a phylogenetic tree using the whole genome protein sequence of the seven species Supplementary (Figure S3) . The dicotyledons wolfberry, tomato, soybean and grape were clustered in both metabolome-based and protein-based trees, while the monocotyledons wheat, rice and maize were clustered together. These results indicate that there are differences in the metabolome between monocots and dicots, and the metabolomes could reflect the evolutionary relationship among different species.

To estimate the classification performance of monocots and dicots, we conducted supervised PCA and partial least-squares discriminant analysis (PLS-DA) among the seven species. The R²X, R²Y and Q² are 0.436, 0.998 and 0.849, respectively, which are the prediction parameters of the PLS-DA model. The R²X value reflects the degree to which the model explains the variability of the input data, while the R²Y and Q² values represent the model’s explanatory power for the output variable and its predictive accuracy, respectively. Here, R²Y and Q² were both close to 0.9, with the same order of magnitude, indicating the stability and trustworthiness of this model (Golbraikh and Tropsha, 2002; Blasco et al., 2015; Su et al., 2022). The results showed that the main principal components (PC1 and PC2) explained 46.6% of the variability, with monocots and dicots being well discriminated from each other ( Figure 3B ).

The PLS-DA loading values of these 587 compounds are listed in Supplementary Table S5 to quantify the contribution to the classification of the metabolites. On the basis of these data, two heatmaps were generated to organize the accumulation levels of the most effective contributors within monocots and dicots. We found that the top 50 contributors preferentially accumulated in dicots ( Figure 3C , according to Supplementary Table S5 , indicated in red), and organic acids and sugars (38%) as well as amino acids and derivatives (22%) were predominant among them, which might be identified as the specific accumulated components of dicots. In contrast, the bottom 50 contributors tended to have a preferential accumulation pattern in monocots ( Figure 3D , Supplementary Table S5 , indicated in blue), and it can be inferred that lipids and phytohormones accumulated specifically only in monocots. The aforementioned results reveal that the differentiation between monocots and dicots can mainly be ascribed to the considerable variation in compounds such as organic acids and sugars, amino acids and derivatives, phenylpropanoids and lipids. Together, these results indicate that the metabolome reflects evolutionary relationships between different crops.

To identify those metabolites that could be used to split and distinguish the seven species, we calculated the fold change for each compound among the species. We defined the metabolites in a certain species whose content was more than 3 times higher than that of other species as species-specific metabolites. There were 16, 29, 14, 21, 15, 52 and 23 species-specific metabolites in wolfberry, rice, wheat, maize, soybean, tomato, and grape, respectively ( Figure 4A ). Seven phenylpropanoids (cyanidin chloride, echinacoside, esculoside, isorhamnetin 3-O-neohesperidoside, methyl p-coumarate, narcissoside and scopoletin), two amino acids and their derivatives (L-asparagine and N-acetylneuraminic acid), four organic acids and sugars (2-hydroxyisocaproic acid, 5-hydroxyhexanoic acid, phenyllactate and shikimic acid), one lipid (cholesterol), one vitamin and coenzyme derivative (riboflavin), and one alkaloid (4-amino-5-hydroxymethyl-2-methylpyrimidine) were among the 16 wolfberry-specific metabolites. Twenty-nine rice-specific metabolites comprised six terpenoids, two phytohormones, four amino acids and their derivatives, five organic acids and sugars, three lipids, eight phenylpropanoids and one vitamin and coenzyme derivative. Six phenylpropanoids, two amino acids and their derivatives, one organic acid and sugar, three lipids, and two nucleic acids and their derivatives are among the 14 wheat-specific metabolites. Among the 21 maize-specific metabolites were two phenylpropanoids, three amino acids and their derivatives, three organic acids and sugars, nine lipids, two nucleic acids and their derivatives, one alkaloid and one other compound. Fifteen soybean-specific metabolites included one phenylpropanoid, two amino acids and their derivatives, three organic acids and sugars, five lipids, one nucleic acid and its derivative, one terpenoid, one alkaloid and one other substance. Fifty-two tomato-specific metabolites included nine phenylpropanoids, 16 amino acids and their derivatives, five organic acids and sugars, 12 nucleic acids and their derivatives, four vitamins and coenzyme derivatives, three alkaloids, one phytohormone and two other compounds. In addition, 23 grape-specific metabolites were identified, comprising 15 phenylpropanoids, three amino acids and their derivatives, three organic acids and sugars, one phytohormone and one other compound.

We found that phenylpropanoids, organic acids and sugars accounted for the largest proportion of specific metabolites in wolfberry, which may determine the special nutritional value of wolfberry. To further compare the metabolites specifically accumulated in wolfberry compared to each species, we performed a Venn diagram analysis. The Venn diagram shows the number of wolfberry-specific metabolites compared with other species. Wolfberry had the most differential metabolites compared with grapes, followed by rice, corn and wheat, the three monocotyledon crops ( Figure 4B ). Venn diagrams of specific metabolites of other species are shown in Supplementary Figure S4 . We found that all species had the highest number of specific metabolites relative to grape, while grape had the most specific metabolites compared to wheat.

The copy numbers of the riboflavin and phenyllactate degradation genes riboflavin kinase (RFK) and phenyllactate UDP-glycosyltransferase (UGT1) were lower in wolfberry than in other species, while the copy number of the phenyllactate synthesis gene hydroxyphenyl-pyruvate reductase (HPPR) were higher in wolfberry, suggesting that the copy number variation of these genes among species may be the main reason for the specific accumulation of riboflavin and phenyllactate in wolfberry. To illuminate the genetic basis of the high metabolite content in wolfberry, we investigated the metabolic pathway of riboflavin and phenyllactate, which are more than five times more abundant in wolfberry than in other species. RFK, encoding riboflavin kinase, which is responsible for the degradation of riboflavin, had the lowest number of copies in wolfberry ( Figure 5 ). This could lead to the weakening of riboflavin’s ability to turn into flavin mononucleotide (FMN), thus increasing the accumulation of riboflavin in wolfberry. HPPR, encoding hydroxyphenyl-pyruvate reductase, which is responsible for synthesizing phenyllactate, has the highest copy number in wolfberry when compared with the other six species ( Figure 6 ). As a result, the amount of phenyllactate in wolfberry may dramatically increase. Moreover, UGT1, encoding a key enzyme in the initial step of phenyllactate degradation, had the fewest copies in wolfberry ( Figure 6 ). This will slow down the reaction of phenyllactate to phenyllactylglucose in wolfberry. Therefore, the increased synthesis and decreased degradation of phenyllactate are jointly responsible for the specific accumulation of phenyllactate in wolfberry ( Figure 6 ). Taken together, our results showed that specific accumulation of riboflavin and phenyllactate in wolfberry may be caused by variation in the number of genes involved in riboflavin and phenyllactate synthesis or degradation.