Four Pueraria varieties were used as test materials. Three cultivated P. thomsonii varieties [Heshui (HS), Fogang (FG), and Guilin (GL)] were collected from the experimental field at Foshan Institute of Agricultural Sciences, Gaoming District, Foshan City, Guangdong Province, China. The P. lobata variety [Zhangjiajie (YS)] was collected from Yanghuping Town, Yongding District, Zhangjiajie, Hunan Province, China. All plant materials were harvested in December 2021.

The samples were prepared and the metabolites were extracted according to the following steps. Roots were freeze-dried and ground into powder. Approximately 0.1 g of powder was dissolved in 1.2 mL of 70% (v/v) methanol. To ensure that the powder was completely dissolved, samples were vortexed for 30 s every 30 min at least 6 times, and then stored at 4°C overnight. After centrifugation at 12,000 rpm for 10 min, the supernatant was collected and filtered (SCAA-104, 0.22 µm pore size; ANPEL, Shanghai, China) before being used for UPLC-MS/MS analysis.

Flavonoids were profiled on a UPLC-ESI-MS/MS system (UPLC, SHIMADZU Nexera X2) with a C18 column. The two solvents for UPLC consisted of 0.04% (v/v) acetic acid in water (A) and 0.04% (v/v) acetic acid in acetonitrile (B). The gradient program was as follows: (A:B) 95:5 (v/v) at 0 min, 5:95 (v/v) at 11.0 min, 5:95 (v/v) at 12.0 min, 95:5 (v/v) at 12.1 min, and 95:5 (v/v) at 15.0 min; flow rate 0.40 mL/min. A 4-µL sample was injected and the column temperature was maintained at 40°C. For MS, the metabolites were detected on a triple quadrupole–linear ion trap mass spectrometer (API 6500 Q TRAP AB Sciex, CA, USA). The ion source gases were injected at a pressure ranging from 25 to 60 psi. Instrument troubleshooting and mass calibration were performed using 10 µM and 100 µM polypropylene glycol. Multiple reaction monitoring (MRM) mode was used to scan the metabolites. Optimized de-clustering potential (DP) and collision energy (CE) values were used for each MRM transition.

The crude protein contents were quantified by the Kjeldahl method with slight modifications (Rizvi et al., 2022). Briefly, 200 mg of each sample was ground into fine powder and placed into a 300-mL digestion tube. Five mL sulfuric acid and 2 g accelerator were added to each tube after wetting the sample. Each mixture was incubated at 250 ° for 30 min on a digester. The temperature was raised to 400 ° when H₂SO₄ decomposed and emitted while smoke. Each tube was removed from the digester when the solution had turned brownish black. Protein contents were analyzed on an automated protein analyzer.

The enzyme gravimetric method (Nishibata et al., 2009) was used to extract dietary fiber from the samples with slight modifications. Briefly, 100 g of each sample was dried and crushed into fine powder. Each powdered sample was transferred into a high-foot bottle designed to detect dietary fiber. Afterwards, 40 mL of 0.05 ml/L MES-Tris buffer was added and stirred until the sample was completely mixed. The samples were treated with different enzymes, and the treated samples were used to detect the soluble and insoluble dietary fiber.

Total starch content was determined according to the acid hydrolysis method (Rose et al., 1991) with minor modifications. Briefly, 200 mg powdered sample was placed into a 15-mL centrifuge tube to remove the fat and soluble sugars from the sample. Initially, the material was washed several times with a total volume of 150 mL absolute ethanol and 10 mL distilled water. The washed material was transferred into a 50-mL conical flask after removal of the ethanol solution. Subsequently, 3 mL of hydrochloric acid was added to the conical flask and the flask was placed in a boiling water bath for 2 h. Each sample was allowed to cool to room temperature before the pH was adjusted to 7.0. The sample was further incubated for 10 min after adding 2 mL of lead acetate solution, transferred to a 50-mL volumetric flask and diluted with water. Afterwards, the material was filtered and the filtrate was used to quantify total starch contents.

The metabolome analysis was performed by Beijing Biomarker Biotechnology Co., Ltd. Compounds were identified by comparing their ionization spectra with a self-built database. Quantitative analysis was performed using the MRM mode. The data were processed using Analyst 1.6.3 software. Unsupervised principal component analysis (PCA) was performed in R-3.1.1. The results of hierarchical cluster analysis (HCA) of the samples and metabolites are presented as heatmaps. HCA was performed in R using normalized signal intensities of metabolites and displayed as histograms. Similarly, supervised multiple regression orthogonal partial least squares discriminant analysis (OPLS-DA) was performed in R using ropls (R-3.1.1). OPLS-DA modeling was validated with a permutation test (200 permutations). The prediction parameters of the OPLS-DA model included R2X, R2Y, and Q2Y; R2X and R2Y represent the interpretation rate of the built model to the X and Y matrices, respectively, and Q2Y represents the prediction ability of the model. When R2Y and Q2Y parameters are closer to 1, the model is more stable and reliable, and can be used to screen differentially accumulated metabolites (DAMs). Based on the results of OPLS-DA, we analyzed the variable importance in projection (VIP) of the OPLS-DA model from the obtained multivariate analysis. DAMs were determined based on VIP value ≥ 1 and p-value < 0.05. The metabolites selected from different comparison groups were presented as volcano plots, plotted in R (R-3.1.1). Each point in the volcano plot represents a metabolite, and the size of the scatter represents the VIP value of the OPLS-DA model. The larger the scatter, the greater the VIP value, and the more reliable the screened differential metabolite is.

We profiled the flavonoid metabolites in one P. lobata variety [Zhangjiajie (YS)] and three P. thomsonii varieties [Heshui (HS), Fogang (FG), and Guilin (GL)] using UPLC-MS/MS. We performed a PCA to better characterize the overall metabolite differences among the four Pueraria varieties and the degree of variability between samples within the same variety. As displayed in the PCA score chart ( Figure 1A ), the contribution by PC1 was 93.10% and that of PC2 was 4.18%. The flavonoid metabolites among the four Pueraria varieties were clearly distinguished on the two-dimensional graph, and replicates within the same group clustered, indicating the repeatability and reliability of our data. Moreover, YS was clearly separated from HS, FG, and GL, indicating that P. thomsonii and P. lobata differ substantially in their metabolite profiles. We identified 279 flavonoid metabolites in the four varieties, consisting of 90 flavonoids, 78 isoflavones, 38 flavonols, 32 dihydroflavones, 15 chalcones, 7 other flavonoids, 6 dihydroflavonols, 6 flavanols, 4 anthocyanins, and 3 proanthocyanidins ( Figure 1B , Supplementary Table S1 ). The flavonoid metabolites in the four varieties are presented as a heatmap following homogenization ( Figure 1C ). The composition of flavonoid metabolites in YS differed notably from those in HS, FG, and GL, indicating that the flavonoid metabolite compositions of the two species differ substantially.

We used supervised orthogonal signal correction to evaluate the differences in metabolite composition between samples within the same variety, and applied an orthogonal partial least squares-discriminant analysis (OPLS-DA) model to emphasize the distinctions among varieties. OPLS-DA is effective at identifying DAMs as it can be used to optimize population differences. R2X, R2Y, and Q2Y are important parameters for evaluating the OPLS-DA model. Values of R2Y and Q2Y closer to 1 reflect a more stable and reliable model. In addition, Q2Y values above 0.5 are indicative of an effective model, while Q2Y values above 0.9 are indicative of an excellent model. We used the OPLS-DA model to compare the flavonoid metabolite composition of YS vs. HS (R2X=0.913, R2Y=1, Q2Y=0.998; Figure 2A ), YS vs. FG (R2X=0.947, R2Y=1, Q2Y=0.999; Figure 2B ), YS vs. GL (R2X=0.920, R2Y=1, Q2Y=0.998; Figure 2C ), HS vs. FG (R2X=0.841, R2Y=1, Q2Y=0.990; Figure 2D ), HS vs. GL (R2X=0.868, R2Y=1, Q2Y=0.989; Figure 2E ), and FG vs. GL (R2X=0.805, R2Y=1, Q2Y=0.981; Figure 2F ). The high values for R2X, R2Y, and Q2Y in all pairwise comparisons indicate that these analyses were repeatable and reliable and were suitable to screen for DAMs.

We used a VIP value ≥ 1 and p < 0.05 as criteria to screen for DAMs. The metabolites that differed between pairs of samples (YS vs. HS, YS vs. FG, YS vs. GL, HS vs. FG, HS vs. GL, and FG vs. GL) are presented as volcano plots in Figure 3 and the quantification are given in Supplementary Table S2 . The blue dots in the figure represent down-regulated DAMs, and the red dots represent up-regulated DAMs. We identified 205 significant differentially accumulated flavonoid metabolites between YS and HS, of which 180 were down-regulated in HS compared with YS and 25 were up-regulated in HS compared with YS; 205 between YS and FG (172 down-regulated in FG compared with YS and 33 up-regulated in FG compared with YS); 210 between YS and GL (176 down-regulated in GL compared with YS and 34 up-regulated in GL compared with YS); 129 between HS and FG (57 down-regulated in FG compared with HS and 72 up-regulated in FG compared with HS); 128 between HS and GL (59 down-regulated in GL compared with HS and 69 up-regulated in GL compared with HS); and 120 between FG and GL (65 down-regulated in GL compared with FG and 55 up-regulated in GL compared with FG). These results verify that the complement of flavonoid metabolites differs substantially between YS, HS, FG, and GL.

We screened for DAMs that were shared among the different varieties, as shown in the Venn diagram in Figure 4 . The number of differentially accumulated flavonoid metabolites in YS vs. HS, YS vs. FG, and YS vs. GL was 205, 205, and 210, respectively; however, we identified 167 DAMs in common between each of these comparisons, reflecting the difference between the one P. lobata variety (YS) and the three P. thomsonii varieties (HS, FG, and GL) ( Figure 4A ). We detected only 54 common DAMs between HS vs. FG, HS vs.GL, and FG vs. GL ( Figure 4B ). These results show that the flavonoid metabolites that underlie the differences between YS vs. HS, YS vs. FG, and YS vs. GL are essentially identical. Furthermore, fewer DAMs were common in comparisons of the three P. thomsonii varieties (i.e., 54) than in comparisons of P. lobata and each P. thomsonii variety (i.e., 167). We next classified the 167 common DAMs into subgroups ( Supplementary Table S3 ), of which 59 isoflavones and 55 flavonoids were the major components ( Figure 5 ). Compared to P. thomsonii, most of the more abundant isoflavones and flavonoids in P. lobata were puerarin, daidzein, genistein, biochanin A, acacetin, apigenin, and tricin, and their glycosyl and methyl derivatives. However, the isoflavones and flavonoids mentioned above displayed no clear difference among the three P. thomsonii varieties. We further analyzed DAMs among the three P. thomsonii varieties. Among these metabolites, 12 metabolites including garbanzol, licoflavonol*, 8-Prenylkaempferol, neobavaisoflavone, and pratensein showed drastic differences in the pairwise comparisons between P. thomsonii varieties ( Table 1 ), thus defining metabolite signatures that can be used as markers to distinguish the three P. thomsonii varieties.

Apart from isoflavones and flavonoids, we detected eight chalcones, four anthocyanins, three flavanols, and 18 dihydroflavones that were differentially accumulated in comparisons of P. lobata and each P. thomsonii variety ( Figure 6 ). Among these metabolites, the four anthocyanins peonidin-3-O-glucoside, delphinidin-3-O-sambubioside, pelargonidin-3-O-glucoside, and rosinidin-3-O-glucoside were more abundant in P. lobata than in P. thomsonii ( Figure 6B ). These results indicate that P. lobata accumulated higher levels of flavonoids than P. thomsonii. We also detected high amounts of anthocyanin in P. lobata. In addition, of the top 10 least abundant DAMs identified in the YS vs. HS, YS vs. FG, and YS vs. GL pairwise comparisons, seven were shared ( Figure 7 ): 6-C-methylkaempferol-3-glucoside, kaempferol-3-O-(2’’-O-acetyl) glucuronide, 4,6-dimethoxyisoflavone-7-O-glucoside, biochanin A-7-O-glucoside (Sissotrin), acacetin-7-O-glucoside (tilianin)*, prunetin-5-O-glucoside, and acacetin-7-O-galactoside*. These substances can therefore be used as markers to distinguish between P. lobata and P. thomsonii.

In addition to flavonoids, nutrients such as starch, dietary fiber, and crude protein contribute to the nutritional value of P. thomsonii and P. lobata. We quantified the dietary fiber, total starch, and crude protein contents of P. thomsonii and P. lobata roots using enzyme gravimetry, acid hydrolysis, and the Kjeldahl method, respectively. The P. thomsonii varieties contained more starch than the P. lobata variety. We observed the lowest total starch content in YS (449.85 ± 2.61 mg/g sample), and the highest starch content in GL (629.93 ± 3.73 mg/g), followed by FG (618.24 ± 3.54 mg/g) and HS (613.50 ± 1.93 mg/g) ( Table 2 ).

P. thomsonii is known for its low fiber content. Indeed, per 100 g of sample, we measured the highest insoluble dietary fiber content in YS (10.23 ± 0.28 g), followed by HS (4.26 ± 0.14 g), FG (3.75 ± 0.35 g), and GL (3.24 ± 0.25 g). Per 100 g of sample, we detected the highest soluble dietary fiber content in YS (3.81 ± 0.30 g), followed by GL (0.88 ± 0.03 g), HS (0.29 ± 0.05 g), and FG (0.17 ± 0.03 g) ( Table 2 ).

The P. lobata variety contained the highest crude protein content of the four varieties tested. Per kilogram of sample, the crude protein content was highest in YS (106.24 ± 0.70 g), followed by GL (81.40 ± 0.59 g), HS (72.82 ± 0.91g), and FG (70.95 ± 0.41 g) ( Table 2 ). These results demonstrate that the nutritional compositions of P. thomsonii and P. lobata differ substantially, while those of the three P. thomsonii varieties are similar to each other.