2.1 Plant material and pretreatment

A total of 56 batches of CR (Table S1 in Supplementary materials) were from two species, of which 44 batches were purchased from the markets and 12 batches were collected from wild plants. The purchased batches were mainly from Sichuan, Inner Mongolia and the three northeast provinces of China while the collected batches were obtained from Heilongjiang. All samples were identified by Prof Yanxu Chang (Tianjin University of Traditional Chinese Medicine) using morphological authentication. The voucher specimens were deposited in State Key Laboratory of Component-based Chinese Medicine (Tianjin, China). The CR samples were powdered with a pulverizer and filtered through a 50-mesh sieve. Each sample powder (0.100 g) was weighed accurately and extracted by an ultrasonator with 4.0 mL of 50 % methanol (v/v) for 40 min at 50 Hz. The extract was centrifuged at 7300 rpm for 10 min. All the supernatant were filtered through a 0.22 µm filter membrane before UHPLC analysis.

2.2 Chemicals and reagents

HPLC-grade methanol and acetonitrile (Fisher, Pittsburg, PA, USA), HPLC-grade formic acid (FA) (Anaqua™, Wilmington, DE, USA), and ultrapure water was prepared by Milli-Q academic ultra-pure water system (Millipore, Milford, MA, USA). Other reagents were analytical grade. Eleven reference standards, including caffeic acid, ferulic acid, isoferulic acid, cimifugin, cimicifugoside, cimifugin-4′-O-β-D-glucopyranoside, cimigenol xyloside, cimigenol-3-O-α-L-arabinoside, cimicidanol-3-O-α-L-arabinoside, acetylcimigenol-3-O-α-L-arabinopyranside, 26-deoxycimicifugoside were purchased from Chengdu Desite Bio-Technology Co., Ltd (Chengdu, China). 2 × Taq PCR Mix, Plant Genomic DNA kit, ddH₂O were obtained from Tiangen Biochemical Technology Co., Ltd (Beijing, China). DNA marker, 6 × loading buffer and GoldView™ were supplied from Takara Biomedical Technology (Beijing) Co., Ltd. The primer was synthesized by Sangon Co., Ltd (Shanghai, China).

2.3 DNA barcoding analysis

Each sample was sprayed by 75 % ethanol solution and gently wiped with degreased cotton to remove impurities. Each sample (0.080 g) was powdered for 10 min at 70 Hz in a tissue grinder (Servicebio, China). After transferring the sample to a new centrifuge tube, the genomic DNA extraction was isolated by Tiangen Plant Genomic DNA kit (Tiangen Biotech, China) with minor modifications. The extracted samples were quantitatively analyzed by NanoDrop spectrometer (Thermo Fischer Scientific, USA) and stored at −20 °C for later use.

The ITS2 sequences were amplified from genomic DNA by polymerase chain reaction (PCR) using universal primers of ITS2F (5′-GCGATACTTGGTGTGAAT-3′) and ITS3R (5′-GACGCTTCTCCAGACTACAAT-3′) (Ren et al., 2014). The PCR was performed in a total volume of 50 μL, containing approximately 50–200 ng template DNA, forward primer (2.5 mM, 2 μL), reverse primer (2.5 mM, 2 μL), 2 × Taq PCR Mix (25 μL), ddH₂O added to 50 μL. The PCR conditions were as follows: 94 °C for 5 min; 40 cycles at 94 °C for 30 s, 56 °C for 30 s and 72 °C for 45 s; 72 °C for 10 min. PCR products (5 μL of each) were detected by electrophoresis on 1.5 % agarose gel in 1 × TAE buffer for 40 min at 80 V. Purified PCR products were sequenced in both directions.

Sequences were assembled by Geneious 9.0.2. Then, the complete ITS2 sequences were annotated and cut based on the ITS2 Database (https://its2.bioapps.biozentrum.uni-wuerzburg.de/). After final alignment in MEGA 6.0, 56 ITS2 sequences were imported into NCBI (https://www.ncbi.nlm.nih.gov/) to preliminarily determine the species attribution. Genetic distance was calculated based on Kimura 2-parameter (K2P) model to evaluate intraspecific and interspecific variation. The phylogenetic tree was constructed by Neighbor-joining (NJ) method with 1000 bootstrap replications to summarize the genetic relationships.

2.4 UHPLC-Q-TOF-MS analysis

The metabolic characteristics and composition characterization of herbal samples were collected on a 1290 UHPLC system together with a 6520 Q-TOF mass spectrometer (Agilent, Santa Clara, CA, USA). The samples were separated on a ZORBAX Eclipse Plus C18 column (4.6 × 100 mm, 1.8 μm, Agilent Technologies, MD, USA) at 35 °C. The mobile phase consisted of solvent A (0.1 % FA-water) and solvent B (acetonitrile) with a gradient elution. The elution gradient of metabolomics was as follows: 0–8 min, 5 %–35 % B; 8–20 min, 35 %–57 % B; 20–40 min, 57 %–81 % B; 40–42 min, 81 %–90 % B. The flow rate was 0.3 mL/min and the injection volume was 3 μL. The other ion source parameters were set as follows: source temperature, 550 °C; drying gas temperature, 325 °C; skimmer voltage, 65 V; fragmentor voltage, 120 V; capillary voltage, 3.5 kV; ion spray voltage, −4.5 kV; collision energy (CE), 10, 35 and 40 V; nebulizer gas pressure, 40 psig; drying gas, N₂; gas flow rate, 10 L/min; detection range, m/z 50–1500 in both positive and negative mode.

The DDA mode containing mass range, precursor ion lists (PILs) and static exclusion range lists were used for component characterization. The main different parameters were as follows. According to the organized compounds database, the mass range were set to 100–400 Da, 400–600 Da, 600–700 Da, 700–800 Da, 800–1200 Da in positive mode in order to acquire as much as possible mass spectrometric information, respectively. Furthermore, the molecular weight of the components in CR was mainly concentrated in the range of 400–800 Da, so 400–600 Da, 600–700 Da and 700–800 Da were set to obtain more compound information. The PILs mainly consist of those with low response and the static exclusion range lists were interfering ions (collected in Table S2). The optimal gradient including 0–15 min, 5 %–100 % B; 15–18 min, 100 % B was adopted for component characterization due to its good peak separation effect and more time-saving.

2.5 Quantitative analysis of chemical markers

The quantitative analysis of the herbal samples were acquired by an Ultimate 3000 High Performance Liquid Chromatography with Diode-Array Detection (HPLC-DAD) (Thermo Fisher Scientific, United States). An Agilent Ultimate AQ-C18 (4.6 × 250 mm, 5 μm, Agilent) was used for subsequent analysis. The mobile phase consisted of 0.1 % FA-water (A) and acetonitrile (B), and the optimal gradient conditions were as follows: 0–10 min, 5 %–20 % B; 10–18 min, 20 %–26 % B, 18–28 min, 26 %B; 28–29 min, 26 %–30 % B, 29–35 min, 30 %B; 35–45 min, 30 %–50 %B. The flow rate was set at 1.0 mL/min, the injection volume was 10.0 μL, the column temperature of 30 °C and the detection wavelength was 310 nm.

The quantification method of the potential markers was validated for linearity, precision, repeatability and recovery according to the guiding principle of the validation of analytical methods of Chinese Pharmacopoeia (Yan et al., 2022). The mixed standard solutions of different concentrations were adopted to gain the standard curves. The precision was evaluated by injecting six consecutive needles of the XSM-14 sample. The repeatability of the method was investigated by inspecting six duplicate samples of the XSM-14. The stability of the samples were computed within 24 h. Furthermore, the recovery was measured by adding half amount of the mixed standard solutions into the samples.

2.6 Multivariate statistical analysis

The multivariate statistical analysis mainly included Principal Component Analysis (PCA), Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA), K-Nearest Neighbor (KNN), Adaptive boosting algorithm (AdaBoost) and Fisher discriminant analysis. Firstly, the data processing and analysis of UHPLC-Q-TOF-MS were carried out by the Agilent MassHunter software B.07.00. Agilent Mass Profinder software (version B.10.00) was applied to peak identify, peak match and then metabolomics variable information was obtained. Secondly, total metabolic variants were used for PCA and OPLS-DA analysis by SIMCA (14.1 version) to classify the 56 batches. Further screening the differential components between C. foetida and C. dahurica (Han et al., 2022). Thirdly, AdaBoost and KNN algorithms were served to calculate the grouping accuracy of the selected markers. Finally, these markers were performed by Fisher discriminant analysis using SPSS software. This model would be applied to distinguish and identify C. foetida and C. dahurica with unknown origin.

3.1 Discrimination of C. foetida and C. dahurica using ITS2 barcoding

PCR amplification was performed after extracting DNA fragments. The PCR results showed that the ITS2 regions of 56 samples were successfully amplified by the universal primers ITS2F/ITS3R (Fig. S1) and high-quality bidirectional sequencing trace files were obtained. After removing the 5.8 s and 28 s rRNA gene sequences at both ends, a total of 56 ITS2 sequences were acquired. Among them, 16 of 56 batches of CR were identified as C. foetida and the rest were C. dahurica. All the sequences were 219 bp in length. According to the analysis of variable sites, C. foetida can be classified into four main haplotypes (F1 ∼ F4), while C. dahurica has three main haplotypes (D1 ∼ D3) (Table 1). The GC-content of C. foetida and C. dahurica were 51.8 %∼53.0 % and 50.2 %∼50.7 %, respectively (Table S3).

In order to construct a consensus phylogenic tree with bootstrap percentages, the Neighbor-Joining (NJ) algorithm was applied to the ITS2 regions using MEGA 6.0 software (Fig. S2). C. foetida and C. dahurica could be clearly distinguished into two groups. Meanwhile, K2P Genetic distance within and between species were calculated by this software. The intra specific distance of C. foetida and C. dahurica from different regions were 0 ∼ 0.0046 and 0 ∼ 0.0093, with an average intra specific distance of 0.0013 and 0.0021, respectively. It indicated that there was a little variation in the genetic process between different regions and also proved that ITS2 sequence, as a DNA barcoding of CR, exhibited good stability. The interspecific distance between the two cultivars of CR was 0.0485, which was far greater than the intraspecific difference (Table S3).

This study requires correct species identification of C. foetida and C. dahurica to ensure accurate elucidation of the chemical differences between the two plants and to discover the impact of species on the quality of medicinal materials. It is usually difficult for inexperienced researchers to employ morphological authentication methods to confirm the original plant species of CR. Recently, DNA barcoding was used to identify accurately the plant species unaffected by external conditions (Gao et al., 2019). Here, the ITS2 regions of 56 samples were successfully amplified and the sequences were obtained. The results showed that the ITS2 could authenticate the original plant species of CR with 100 % success rate. The ITS2 region had the potential to be a good DNA barcoding for identification of medicinal species of CR. It could not provide differences in composition and its content between the two species. Subsequently, the differential components were screened and verified to lay the foundation for further improvement on the quality evaluation of CR based on metabolomics and chemometrics.

3.2 Identification of chemical composition in CR

3.2.2

3.2.2 Identification of triterpenoid saponins

Triterpenoid saponins were the primary bioactive components of CR. Up to date, approximately 400 triterpenoid saponins (mostly 9,19-cycloartane type) had been discovered and characterized from the Cimicifuga genus. In our study, total 101 constituents had been identified by means of UHPLC-Q-TOF-MS based on the above approach. Most triterpenoid saponins were liable to form [M + H]⁺ ion and [M + Na]⁺ ion in positive mode. The main cleavage pathways of triterpenoid saponins were prone to lose water, acetyl groups, dimethylethylene oxide and the glycosyl groups, resulting in neutral losses of 18.01 Da, 60.02 Da, 72.01 Da, 132.05 Da, 162.05 Da. The following were examples of the conventional fragmentation pathways for triterpenoid saponins.

The main cleavage pathway of cimigenol-type of triterpenoid saponins was to lose water, the glycosyl groups and was prone to twist to form dimethylethylene oxide. Compound c106 cimigenol-3-O-α-L-arabinoside had a [M + Na]⁺ peak at m/z 643.38 (1.54 ppm). Then m/z 585.37 ([M + H − 2H₂O]⁺) and m/z 453.32 ([M + H − 2H₂O − Ara]⁺) were formed after successively removing two molecules of water (36.02 Da) and arabinose (132.05 Da). And m/z 435.31 ([M + H − 3H₂O − Ara]⁺) was also detected after removing a molecule of water (Fig. 2). In the positive mode, the precursor ion of compound c126 cimigenol-3-O-β-D-xyloside was m/z 643.38 [M + Na]⁺, and the molecular formula was presumed to be C₃₅H₅₆O₉. m/z 567.36, m/z 495.35, m/z 363.26 were generated successively with continuous water loss, dimethylethylene oxide and xylopyranose. Compound c69 12β-hydroxy-7,8-didehydrocimigenol 3-O-α-L-arabinopyranoside had a [M + H]⁺ at m/z 635.37 (2.95 ppm). In the secondary mass spectrometry, fragments of m/z 617.37 [M + H − H₂O]⁺, m/z 599.36 [M + H − 2H₂O]⁺, m/z 545.31 [M + H − H₂O − C₄H₈O]⁺, m/z 467.31 m/z [M + H − 2H₂O − Ara]⁺ and m/z 377.26 [M + H − 3H₂O − Ara − C₄H₈O]⁺ were produced (Pang et al., 2021).

Close

Fig. 2

Illustration for the structural elucidation of triterpenoid saponins CR.

Export to PPT

The main feature of 16,23-diketo-type is that the C-16 and C-23 positions both are oxidized to carbonyls, and partially dehydrated to form a ternary oxygen ring structure at C-24 and C-25 positions. Therefore, this type of compounds is extremely easy to remove dimethylethylene oxide and produce highly responsive m/z 73 [C₄H₈O + H]⁺. The [M + H]⁺ peak of compound c54 cimicidanol-3-O-α-L-arabinoside was at m/z 617.36 (−1.44 ppm). Aglycones were produced with dimethylethylene oxide, arabinose and continuous water losses to generate fragments of m/z 545.31 [M + H − C₄H₈O]⁺, m/z 467.31 [M + H − H₂O − Ara]⁺, m/z 395.25 [M + H − H₂O − C₄H₈O − Ara]⁺ and m/z 251.17 [M + H − 3H₂O − C₄H₈O − Ara − C₇H₈O]⁺ (Fig. 2). From above, compound c151 was speculated as cimicidanol according to the fragments of m/z 485.32 [M + H]⁺, m/z 413.27 [M + H − C₄H₈O]⁺, m/z 395.25 [M + H − H₂O − C₄H₈O]⁺. Compound c60 was identified as cimicifugoside H-2 (t_R = 6.561, C₃₅H₅₄O₁₀), as lost a xylopyranose (132.05 Da) and a molecule of dimethylethylene oxide (72.01 Da) and continuous water (n∙18.01 Da) at positive conditions (Cao et al., 2005, Li et al., 2007).

Shengmanol-type of Cimicifugae Rhizoma contains acetyl groups and dimethylethylene oxide at the end of carbon chain, which are easily lost in secondary mass spectrometry. The MS² spectrum of 23-O-acetylcimigenol-3-O-α-L-arabinoside (consistent with peak 148, t_R = 11.55 min) was showed the precursor ion at m/z 663.41 [M + H]⁺. The fragments at m/z 585.36 and 435.31 indicated the elimination of C₂H₄O₂, Ara and three molecules of water. A serious of fragments at m/z 529.3, 469.3, 451.3, 397.2 were both found in peak 78 and 152, indicating that those two compounds had the same skeleton and highly similar in structures. Therefore, they were identified as 23-O-aectyl-7,8-didehydroshengmanol 3-O-α-L-arabinopyranoside and 24-O-acetylshengmanol-7(8)-en-isodahurinol, respectively.

3.2.3

3.2.3 Identified of phenylpropanoids

A total of 49 phenylpropanoids were tentatively identified in the positive and negative mode. The compound of c13, c32 and c36 were orderly identified as caffeic acid, ferulic acid and isoferulic acid by comparing the authentic standards, which cleavage pathways were specified in Fig. 3. The fragment of phenolic acids was characterized by neutral losses of 15.02 Da (–CH₃), 18.01 Da (–H₂O) and 44.01 Da (–CO₂). Based on the same fragment ions m/z 179.03,135.04 as caffeic acid, peak 6 (t_R = 2.518 min, C₁₅H₁₈O₉) was characterized as caffeic acid 4-O-β-D-glucopyranoside.

Close

Fig. 3

Illustration for the structural elucidation of phenylpropanoids and chromones from CR.

Export to PPT

Cimicifugic acids mainly generate signal responses in positive ion mode and were prone to neutral losses such as CO₂, CO, H₂O, and CH₃. Due to the carbon chains of these components contained hydroxyl and carboxyl groups, they were easy to fracture on the carboxyl group oxygen in collision energy spectra. A fragment response at m/z 177 was likely to occur when methoxy group was present on the benzene ring at the same time. Peak 19 (t_R = 3.96, C₂₆H₃₀O₁₂) was regarded as cimicifugaside F attributing to the fragment ions m/z 535.18, m/z 517.21, m/z 491.29 and m/z 163.07. According to the ions m/z 433.11, m/z 389.16, m/z 355.14, m/z 177.05, peak 48 (t_R = 5.83, C₂₁H₂₀O₁₀) was identified as 2-feruloyl piscidic acid. Compared with the literature, peaks 1, 2, 39 and 40 were successively identified as fukiic acid, piscidic acid, cimicifugic acid A and cimicifugic acid B (Werner and Petersen, 2019).

3.2.5

3.2.5 Identification of compounds in C. foetida and C. dahurica

The chemical compositions identification of C. foetida and C. dahurica was a basis for screening potential differential components. The TIC of the two species both in positive and negative modes were shown in Fig. 4. In this work, a total of 88 chemical components were ultimately authenticated in the methodology of metabolomics referring to the above qualitative analysis, including 65 from C. foetida and 75 from C. dahurica (Table S4). According to the constituents identified from the two species, there were significant differences in composition and content between them.

Close

Fig. 4

The total ion chromatograms (TIC) of the two species in positive and negative modes. (A) TIC of C. dahurica in positive MS mode; (B) TIC of C. dahurica in negative MS mode; (C) TIC of C. foetida in positive MS mode; and (D) TIC of C. foetida in negative MS mode.

Export to PPT

A comprehensive insight into the secondary metabolites of various species in TCMs is vital for further metabolomics analysis. The DDA mode could maximize the MSⁿ information collection, greatly reducing the difficulty of acquiring low abundant components (Zuo et al., 2019). And the CIF platform was established to rapidly analyze and match a compound from a large amount of MSⁿ information during the data processing stage. Based on this, a total of 157 compounds were tentatively characterized in CR, including 101 triterpenoid saponins, 44 phenylpropanoids and 7 chromones. Compared with previous composition analysis methods, this strategy is more rapid, accurate and efficient, which could become the principal method to identify compounds soon. Furthermore, strengthening this method development and application to provide a practical strategy for characterizing non-targeted metabolites in TCMs.

3.3 Metabolomic analysis for discrimination of C. foetida and C. dahurica

The retention time and peak area of ten randomly selected ion pairs in QC samples were acquired to verify the UHPLC-Q-TOF-MS method. The relative standard deviations (RSDs) of them were less than 5.0 %, indicating the accuracy and reliability of this analytical method. All QC samples were closely clustered into a group in PCA, demonstrating the reproducibility of the analytical system (Fig. 5A and B).

Close

Fig. 5

Multivariate statistical analyses of the two species: PCA score plot in positive mode (A); PCA score plot in negative mode (B); OPLS-DA score plot (C: 2955 ionic characteristic variables in positive mode, D: 2976 ionic characteristic variables in negative mode, E: 48 compounds characteristic variables, F: 4 combinatorial discriminatory quality markers); heat-map in two species(G).

Export to PPT

Metabolomics analysis exhibited excellent properties on screening the differential components in natural products. Total 2955 and 2976 ionic characteristic variables were obtained in positive and negative mode, respectively. The unsupervised model PCA and supervised model OPLS-DA were performed to differentiate the two species in metabolite levels with high fitting and prediction degree. The results showed that the samples from the two species were successfully separated into two distinct clusters in PCA with high fitting and predictive abilities both in positive and negative mode (Fig. 5A and B). The 2955 metabolic variables in positive mode were evidently classified the samples into two groups with a goodness-of-fit R²Y = 0.983 and goodness-of-prediction Q² = 0.95 (Fig. 5C). The 2976 metabolic variables in negative mode had the similar result with R²Y = 0.972 and Q² = 0.935 (Fig. 5D). To better classify and account for the two species, variable importance in projection (VIP) combined with t-test were applied in OPLS-DA mode for significance testing. In this work, about 48 distinctive components were screened and identified between the two species based on analyzing the criteria of VIP > 1 and p < 0.05, including triterpenoid saponins, phenylpropanoids and chromones (Table 3). The 48 compounds characteristic variables could distinguish the samples into two groups in OPLS-DA model (Fig. 5E) and the results were visualized in the Heatmap (Fig. 5G). Among these compounds, the content of 12 components, such as cimifugin-4′-O-β-D-glucose (4), cimicifugoside (5), caffeic acid (7), cimifugin (10), norcimifugin (17) and 26-dedoxycimifugoside (31) were generally higher in C. foetida, while others were higher in C. dahurica. All of the above highlighted significant differences in composition and content between C. foetida and C. dahurica. In addition, relevant studies have shown that C. foetida has significant effects in antidiarrheal, anticomplementary and menopausal syndrome effects (Qiu et al., 2006, Zheng et al., 2013, Zhang et al., 2016), while C. dahurica has better antioxidant, neuroprotective and antibacterial effects (Qin et al., 2016, Lee et al., 2020, Li et al., 2023). Therefore, it is necessary to conduct differential analysis between the two species in order to apply it more clearly in clinical practice.

Considering the above content differences and the screening principle of biomarkers: 1) the markers are convenient to obtain and quantify; 2) the markers can clearly distinguish between the two original CR; 3) the specific components in CR (Cui et al., 2022, Lu et al., 2022). Caffeic acid, cimifugin, ferulic acid and isoferulic acid were ultimately screened as potential combinatorial discriminatory quality markers. The two species were significantly divided into two groups by the four markers in OPLS-DA model with high fitting and predictive abilities (Fig. 5F). Furthermore, the grouping accuracy were verified based on AdaBoost and KNN algorithms by Matlab R2022A software. The sample SM-1, 4–14, XSM-1, 2, 7–25 were screened as the training set while others were the testing set at random. The accuracy of variables in different groups were above 83.3 %, which proved the correctness of OPLS-DA model grouping and verified the existence of differences between the two Cimicifuga species (Table 4).

3.4 The verification of the potential bioactive markers

Four combinatorial discriminatory quality markers were quantitatively analyzed by HPLC-DAD with high efficiency and generality. It was indicated that the contents of the four markers from 56 batches were totally different, which may contribute to the differentiation of the species. The repeatability, stability and precision of the four compounds were less than 5.0 % and recoveries were between 96.1 % and 103 % with all RSD values less than or equal to 5.55 %. The correlation coefficients of the linear equations higher than 0.99 and the lower limit of quantitation (LOQ) were 0.08 μg/mL of caffeic acid and cimifugin, 0.048 μg/mL of ferulic acid and 0.16 μg/mL of isoferulic acid, respectively (Table S5). It demonstrated a good linear relationship between these four compounds within their respective concentration ranges. The content of four combinatorial discriminatory quality markers of 56 batches were obtained (Table S6) and the intuitive chart was shown in Fig. 6. In summary, the established HPLC-DAD quantitative analysis for four markers was unequivocally accurate and reliable.

Close

Fig. 6

The contents chart of four markers in all batches of samples.

Export to PPT

Fisher discriminant model was established to classify unknown samples using SPSS software in terms of the above contents of four combinatorial discriminatory quality markers. The twenty batches (including SM-2, 4, 5, 6, 7, 8, 9, 11, 13, 14, XSM-2, 10, 14, 15, 17, 18, 21, 22, 24, 28) were randomly chosen as training sets and the fifteen batches (SM-1, 10, 15, XSM-1, 3, 4, 5, 6, 7, 8, 9, 16, 19, 23, 25) were test sets. The discriminant equation was derived as follows: Y = 7.029X₁ + 1.637X₂ − 0.193X₃ − 0.125X₄ − 3.613 (Y: discriminant score; X₁: caffeic acid; X₂: cimifugin; X₃: ferulic acid; X₄: isoferulic acid). The unknown samples will be classified as C. foetida when the discriminant score is higher than the determination value 0 (the average of 3.188 and −3.188 at the group centroids), if not it will be considered as C. dahurica. In addition, the other 21 batches were utilized to the test sets in order to verify the discrimination ability. The results displayed that the accuracy of classification was 94.4 % in cross-validation group demonstrating high feasibility of this model. The batches were correctly grouped by the discriminant model except two samples. Apart from a little higher content of isoferulic acid, the other three components are generally lower in the two batches.

It was suggested that the four combinatorial discriminatory quality markers could distinguish C. foetida and C. dahurica with high accuracy. The content of ferulic acid in C. foetida was higher than that in C. dahurica, while isoferulic acid was the opposite. Caffeic acid was three times higher in C. foetida than in C. dahurica. Moreover, cimifugin was significantly different between the two species, and it was nearly 16 times higher in C. foetida than in C. dahurica. These differences in the composition of secondary metabolites might be attributed to genetic nuances. Consequently, the contents of four markers will be measured and substituted into the discrimination model so as to differentiate the unknown samples. All above indicated that metabolomics techniques can be applied to distinguish different species from the perspective of compositional differences.