Fifty-one leaf samples were collected from six Taxus species in China and stored at 4°C (S1–15: T. mairei, S16–20: T. chinensis, S21–25: T. yunnanensis, S26–30: T. wallichiana, S31–36: T. cuspidata, S37–42: T. media, S43–51: for the external validation of FDA models). The images and detailed information of the samples are shown in Supplementary Figure S1 and Table 1. The sources of all samples were identified by Professor ZH and MY from the Plant Identification Teaching and Research Office of the Fujian University of Traditional Chinese Medicine, Fujian Province, China. Voucher specimens were kept in the Comprehensive Medical Research Institute of Fujian University of Traditional Chinese Medicine.

Methanol, acetonitrile, and formic acid for UPLC analysis were purchased from Merck (Darmstadt, Germany). Deionized water was prepared daily by using a Millipore Milli-Q purification system (Millipore, Bedford, MA, United States). 10-DAB, baccatin III (BAC), 7-xylosyl-10-deacetyltaxol (7-xyl-10-DAT), 10-DAT, cephalomannine (CE), 7-epi-10-deacetyltaxol (7-epi-10-DAT), paclitaxel (TAXOL), 7-epi-paclitaxel (7-epi-TAXOL), gallocatechin (GC), catechin (C), isoquercitrin (IQC), nicotiflorin (NFR), and triptolide (IS₁) were purchased from Chengdu Mansite Bio-Technology Co., Ltd. (Chengdu, China). Epigallocatechin (EGC), taxifolin (TAX), aromadendrin (ARO), amentoflavone (AF), 7-demethylginkgetin (DGK), GK, isoginkgetin (IGG), and SDN were purchased from Baoji Herbest Bio-Technology Co., Ltd. (Baoji, China). Epicatechin (EC), rutin (RT), QC, quercetin (QR), casticin (IS₂), and liquiritin (IS₃) with purities exceeding 98% (determined by HPLC) were purchased from the China National Institute for Food and Drug Control (Beijing, China), and their purity was more than 98% (determined by HPLC). Figure 1 shows the chemical structures of the 24 analytes and the three internal standards.

Dried Taxus leaf samples were ground into 65 mesh powder, weighed accurately to 0.10 g, suspended in 25 mL of 80% methanol aqueous solution, accurately weighed in a 50 mL capped conical flask, and extracted in an ultrasonic bath with an output power of 300 W (40 kHz) for 30 min. The extract was added with 80% methanol to compensate for weight loss. The extract was centrifuged (12,000 r/min, 10 min) then passed through a 0.22 μm microporous membrane. The sample was diluted 5-fold by adding 800 μL of methanol solution to 200 μL of the sample solution. Then, the diluted sample was added to the internal standard mixed solution at a 1:1 v/v ratio. All samples were kept at 4°C.

The 24 standards were weighed separately and dissolved in methanol (UPLC–MS grade) to obtain single stock solutions with accurate concentrations. Each stock solution was rediluted and mixed with methanol to obtain a series of working standard solutions, and calibration curves were established by using the mixed working standard solution. All standard solutions were stored in brown glass bottles at 4°C.

Internal standards, namely, triptolide (IS₁), casticin (IS₂), and liquiritin (IS₃), were dissolved in methanol to a concentration of approximately 1 mg/mL individually. Each stock solution was rediluted with methanol to prepare a mixed internal standard solution with 50 ng/mL IS₁, 50 ng/mL IS₃, and 100 ng/mL IS₂. A total of 500 μL of the mixed internal standard was added to 500 μL of the mixed working standard solution or sample solution and filtered through a 0.22 μm micropore membrane before use.

Chromatographic analysis was performed on a Waters UPLC system (Waters, United States). Twenty-four analytes were chromatographically separated on a Waters Cortecs C₁₈ column (2.1 × 100 mm, 1.6 μm) through chromatographic separation. The column temperature was 45°C. The mobile phases consisted of 0.1% formic acid in water (A) and acetonitrile (B). The gradient elution program was as follows: 95% A at 0–0.5 min, 95%–55% A at 0.5–2.5 min, 55% A at 2.5–4.5 min, 55%–25% A at 4.5–8.5 min, 25%–95% A at 8.5–8.6 min, and 95% A at 8.6–10.5 min. The flow rate was 0.2 mL/min, and the volume of the injected sample was 2 μL.

Mass spectrometry analysis was conducted with a Waters TQS triple quadrupole mass spectrometer in the switching mode of electrospray positive- and negative-ion modes of multiple reaction monitoring (MRM). The optimized MS conditions were fixed as follows: capillary voltage, 2.50 kV; desolvent gas flow, 800 L/h (N₂); desolvent gas temperature, 500°C; ion source temperature, 150°C; secondary cone hole extraction voltage, 3.00 V; cone gas flow, 50 L/h (N₂); collision gas, argon. The most suitable collision energy (COE) and cone voltage (CV) for each analyte was optimized. See Table 2 for the specific parameters.

The quantitative data of the analytes and FDA were statistically analyzed with SPSS statistical software 22.0 (SPSS Inc., Chicago, United States), and the difference was significant. Data were analyzed by Masslynx version 4.2. HCA was performed by using Heatmap Illustrator 1.0. PCA and OPLS-DA were conducted with SIMCA 14.1 software (Umetrics, Umea, Sweden). RF was established with the RF package (version 4.6–14) in the R environment (version 3.5.2).

In this experiment, we selected methanol–water, methanol–0.1% formic acid water, acetonitrile–water, and acetonitrile–0.1% formic acid water as the candidate mobile phases and the chromatographic peak separation effect and peak shape as the evaluation indicators. Acetonitrile–0.1% formic acid water was selected on the basis of the optimized result. Our established method separated five isomeric compounds (GC and EGC, C and EC, GK and IGG, 10-DAT and 7-epi-10-DAT, and TAXOL and 7-epi-TAXOL) simultaneously for the first time. To develop a sensitive and accurate quantitative method, we studied the quantitative ion selection of analytes in the positive- and negative-ion ESI modes on the basis of the optimized CV and COE values by using manually optimized MRM parameters, as shown in Table 2. The chromatograms of the representative T. mairei sample (S3, methodological validation sample) and mixed standard solution are provided in Supplementary Figure S2 and Figure 2, respectively. Among the 24 analytes measured, we found that four different classes of flavonoids had the best response in the negative-ion mode. Eight taxoids had the best response in the positive-ion mode. IS₁ was used as the internal standard of taxoids; IS₂ was utilized as the internal standard of flavanols, flavonols, and biflavones; and IS₃ was applied as the internal standard of dihydroflavones. ISs could be separated from the chromatographic peaks of the 24 analytes. This characteristic plays an important role in the accurate determination of complex components to ensure the accuracy of results. ISs and the analytes had similar chemical structures, which are indicated by the same color. IS₁ and taxoids are marked with blue; IS₂, flavonols, and biflavones are marked with black; and IS₃, flavanols, and dihydroflavones are marked with red. Refer to Figure 1 for specific analyte results.

We employed the established analytical method to identify the 24 compounds in six Taxus species. Structures were unambiguously assigned on the basis of the retention times and MS spectra of the reference standards (Supplementary Figure S4). The ESI mass spectra provided the characteristic quasi-molecular ions of GC [M–H]⁻ at m/z 304.90; EGC [M–H]⁻ at m/z 304.90; C [M–H]⁻ at m/z 289.01; EC [M–H]⁻ at m/z 289.01; RT [M–H]⁻ at m/z 608.81; IQC [M–H]⁻ at m/z 463.20; NFR [M–H]⁻ at m/z 593.15; QC [M–H]⁻ at m/z 447.10; QR [M–H]⁻ at m/z 301.08; TAX [M–H]⁻ at m/z 303.03; ARO [M–H]⁻ at m/z 287.05; AF [M–H]⁻ at m/z 536.90; DGK [M–H]⁻ at m/z 551.00; GK [M–H]⁻ at m/z 565.00; IGG [M–H]⁻ at m/z 565.00; SDN [M–H]⁻ at m/z 578.83; 10-DAB [M–H]⁻ at m/z 567.20; BAC [M + Na]⁺ at m/z 609.10; 7-xyl-10-DAT [M + Na]⁺ at m/z 966.20; 10-DAT [M + Na]⁺ at m/z 834.20; CE [M + Na]⁺ at m/z 854.20; 7-epi-10-DAT [M + Na]⁺ at m/z 834.00; TAXOL [M + Na]⁺ at m/z 876.10; and 7-epi-TAXOL [M + Na]⁺ at m/z 876.20. Moreover, the characteristic fragmentation behaviors of the 24 compounds revealed by MS/MS analysis in our present study were identical to those in previous studies (Wang Y.-J. et al., 2019; Gai et al., 2020; Luo et al., 2020; Yi et al., 2020) and our previous pre-experimental HPLC–Q-TOF–MS/MS work (Supplementary Figure S3; Supplementary Table S2): m/z 305→219→179→125 for GC and EGC; 289→245→203→151→125 for C and EC; m/z 609→301→257→151 for RT; m/z 463→300→179→151 for IQC; m/z 593→285→151 for NFR; m/z 447→301→179→151 for QC; m/z 303→285→151→107 for TAX, m/z 287→259→243→201 for ARO; m/z 301→273→179→151→107 for QR; m/z 567→531→445 for 10-DAB; m/z 537→443→417→399→375 for AF; m/z 609→549→427→367 for BAC; m/z 551→519→475→457→389 for DGK; m/z 966→681→308 for 7-xyl-10-DAT; m/z 834→549→308 for 10-DAT and 7-epi-10-DAT; m/z 565→533→389→374 for GK and IGG, m/z 854→591→286 for CE; m/z 876→591→533→308 for TAXAL and 7-epi-TAXOL; m/z 579→547→403→165 for SDN, and GC, RT, TAX, SDN, and TAXOL were used as examples to depict the detailed identification processes of flavanols, flavonols, dihydroflavones, biflavones, and taxoids in Figures 3A–E, respectively.

Eight taxoids, four flavanols, five flavonols, two dihydroflavones, and five biflavones in 42 leaf samples from Taxus species were determined and analyzed by UPLC–MS/MS after method validation (See Figure 4; 1: T. mairei, 2: T. chinensis, 3: T. yunnanensis, 4: T. wallichiana, 5: T. cuspidata, 6: T. media). Under the UPLC–MS/MS condition, all analytes can be detected in the leaves of the different Taxus species. The contents of chemical components in different Taxus species were obviously different. We first discuss toxoids. In T. wallichiana, the amount of 7-Epi-TAXOL was less than the LOQ. Among the eight taxoids, 10-DAB, BAC, 7-xyl-10-DAT, 10-DAT, and TAXOL were higher in T. yunnanensis than in the other five Taxus species; CE, TAXOL, and 7-epi-TAXOL were higher in T. cuspidata than in the other five Taxus species; and 7-epi-10-DAT was higher in T. chinensis than in the other five Taxus species. The eight other taxoids were abundant in T. yunnanensis and had low abundance in T. wallichiana (T. yunnanensis > T. cuspidate > T. chinensis > T. media > T. mairei > T. wallichiana). Therefore, T. yunnanensis was used as a plant material for taxoid extraction. Second, we discuss flavanols. We found that among the four flavanols, GC was higher in T. yunnanensis than in the other Taxus species, and EGC, C, and EC were higher in T. mairei, T. chinensis, and T. media than in the other three Taxus species. Third, we found that among the five flavonols, QR was present at low levels in the six Taxus species. RT, IQC, NFR, and QC were higher in T. maireis, T. chinensis, and T. media than in the other Taxus species. Fourth, we discovered that of the two dihydroflavones, TAX and ARO were higher in T. cuspidata and T. media than in the other Taxus species. Fifth, among the five biflavones, AF was higher in T. mairei, T. chinensis, and T. media than in the other Taxus species. DGK was higher in T. cuspidata, GK was higher in T. chinensis, IGG was higher in T. wallichiana, and SDN was higher in T. mairei than in the other Taxus species. The components in the leaves of different Taxus species were different, and quantitative analysis could aid the effective use of the leaves of different Taxus varieties. Moreover, we comparatively analyzed the quantitative results of 24 analytes from the six Taxus species on the basis of their mean ± SDs in Figure 4 , and the content determination results are provided in Table 7.

HCA in Heatmap Illustrator 1.0 software was used to analyze the contents of the 24 analytes in the leaves of six Taxus species, and the content data were normalized. The clustering method was the Ward method, and the distance type was square Euclidean distance. The HCA results are shown in Figure 5. The 42 batches of leaf samples from the six Taxus species divided into six categories. Therefore, the established content determination method can distinguish the six different varieties of Taxus species. The differences in the contents of the 24 analytes in the leaves of different Taxus species can be clearly seen in Figure 5. Among the eight taxoids, 10-DAB, BAC, 7-xyl-10-DAT, and 10-DAT were higher in T. yunnanensis; CE, TAXOL, and 7-epi-TAXOL were higher in T. cuspidata; and 7-epi-10-DAT was higher in T. chinensis than in the other species. The contents of four flavanols and five flavonols were higher in T. chinensis and T. media than in the other Taxus species. However, GC was higher in T. yunnanensis than in the other Taxus species. The contents of two dihydroflavones were higher in T. cuspidata and T. media than in the other Taxus species. Among the biflavones, DGK was higher in T. cuspidata than in the other species. GK content was higher in T. chinensis, IGG content was higher in T. wallichiana, and SDN content was higher in T. mairei than in the other Taxus species. HCA revealed the differences in the contents of the 24 analytes in leaf samples from six different Taxus species with high intuitiveness.

PCA is the most commonly data analysis method in the statistical analysis of multiple elements. In this method, the information of original features is retained to the maximum extent without the loss of important information (Jolliffe and Cadima, 2016; Giuliani, 2017). We imported the data into Simca 14.1 with sample batch as the observation value and the 24 analytes as the variables. Unit variance scaling was used for data normalization to generate the PCA model. The scores and accumulated R²X and Q² with different numbers of principal components are displayed in Figures 6A, B. Consistent with the HCA results, the 42 batches of Taxus leaf samples divided into six categories. The samples were well separated, indicating differences in the contents of the 24 analytes in the leaves of different Taxus species (Figure 6A). The main components of the PCA model can predict 92.8% of the changes in the original dataset (R²X [cum] = 0.928), and 73.4% of the cumulative prediction rate was predicted in the model through cross-fold validation (Q² [cum] = 0.734, Figure 6B). Therefore, this model can represent the information of the original data, and the obtained results are scientific and effective. In the next step, we used OPLS-DA to identify the chemical components that distinguish the six Taxus species.

OPLS-DA is a supervised discriminant analysis method that can be used to find variables that cause differences between samples (Li et al., 2021; Wang et al., 2022). By using PCA, we further searched for the chemical components that distinguish the different Taxus species. Once again, we imported the content determination results of the 24 analytes in the 42 batches of Taxus leaf samples into SIMCA 14.1 for OPLS-DA. The OPLS-DA scores and the variable importance predictive values (VIPs) of the 24 analytes are presented in Figures 6C, D. The established OPLS-DA model can predict 97.4% of the information in the X matrix (R²X [cum] = 0.974) and 87.6% of the information in the Y matrix (R²Y [cum] = 0.876). Cross-validation (Q² [cum] = 0.854) revealed that our model is stable and has good prediction ability. The OPLS-DA score plot was similar to the PCA score plot, which divided the six Taxus species into six categories. This result also verified that their compositions were different (Figure 6C). Two hundred permutation tests were used for the internal validation of the model to prevent overfitting. The R ² and Q ² values of the six varieties were greater than 0.9 and were higher than the left (negative) of the R ² and Q ² values, indicating that the model and results are reliable without overfitting (Supplementary Figure S5). Therefore, VIP was used to further analyze the differential components of the six Taxus species (Figure 6D). The VIP indicates the contribution of each analyte to distinguishing samples, and the variables with substantial contribution to grouping and VIPs greater than 1.0 were screened as the index (Giuliani, 2017). The following analytes with VIPs greater than 1.0 can be used to distinguish the six different Taxus species: 10-DAB, 10-DAT, ARO, GK, BAC, GC, TAX, TAXOL, QC, EGC, 7-xyl-10-DAT, CE, C, and RT. The above analytes can be applied as potential index components to distinguish the quality of the six different Taxus species.

The 14 analytes selected by OPLS-DA can be utilized as the chemical components distinguishing the six Taxus species. However, the detection of numerous analytes requires great cost and time. In addition, the data scaling method has a strong effect on the model. RF will not excessively scale the data under the condition of mutual comparison (Zhang et al., 2020). Additional stable components were simplified to simplify the model and distinguish the six Taxus species. The RF model was selected for 100 iterations, ranking the importance of the 14 analytes each time. The cumulative times of each variable in the top Nth positions were obtained in accordance with the importance ranking of the variables in 100 RF model iterations. In Figure 7A, each line represents an analyte, and the line on the left has high characteristic importance, whereas the line on the right has the opposite characteristics. Afterward, the frequency of the variable that is the most important parameter in the 100-iteration modeling process was calculated (Figure 7B). The variable with the highest frequency was EGC, which accumulated 48 times, indicating that this analyte has a remarkable role in distinguishing the six Taxus species. Among other analytes, ARO accumulated 23 times, 7-xyl-10-DAT accumulated nine times, 10-DAB accumulated eight times, GK accumulated seven times, and RT accumulated five times. Therefore, we selected EGC, ARO, 7-xyl-10-DAT, 10-DAB, GK and RT to establish the model and evaluate the feasibility of simplifying the model through FDA.

FDA is a typical pattern discrimination method. We employed FDA to verify the feasibility of classifying six Taxus varieties by using the six analytes, namely, EGC, ARO, 7-xyl-10-DAT, 10-DAB, GK, and RT, selected by RF. The discriminant function was produced in accordance with the content data of the six analytes in the 42 sample batches, as follows: Y(T. mairei) = 0.338 7-xyl-10-DAT + 0.064 GK + 0.027 RT −0.008 ARO + 0.130 10-DAB + 0.006 EGC −68.371; Y(T. chinensis) = 0.317 7-xyl-10-DAT + 0.263 GK + 0.065 RT −0.072 ARO + 0.192 10-DAB + 0.007 EGC −389.484; Y(T. yunnanensis) = 1.065 7-xyl-10-DAT + 0.003 GK + 0.044 RT −0.080 ARO + 0.790 10-DAB −0.016 EGC −512.021; Y(T. wallichiana) = 0.005 7-xyl-10-DAT + 0.002 GK + 0.001 RT + 0.000 ARO + 0.002 10-DAB + 0.001 EGC −1.885; Y(T. cuspidata) = 0.010 7-xyl-10-DAT + 0.043 GK − 0.002 RT + 0.683 ARO + 0.005 10-DAB + 0.024 EGC −187.611; Y(T. media) = 0.042 7-xyl-10-DAT + 0.057 GK + 0.015 RT + 0.480 ARO + 0.044 10-DAB + 0.065 EGC −172.673. Additionally, we selected samples S43–S51 outside the model and determined only the content of the six analytes (the content data results are shown in Table 8) for external validation. The samples of the established model and those outside the model based on the discriminant function of the linear combination of predictive variables were correctly discriminated and classified (Table 9), thus verifying the accuracy and practicality of the simplified model. We compared the differences in the chemical structures of the six analytes from Taxus species. 7-Xyl-10-DAT and 10-DAB were taxoids, and their distinguishing radical pharmacophores were hydroxyl and xylose groups; EGC was a flavanol, and its distinguishing radical pharmacophore was hydroxyl group; ARO was a dihydroflavone, and its distinguishing radical pharmacophore was C-ring α-hydroxyl and B-ring hydroxyl groups; GK was a biflavone, and its distinguishing radical pharmacophores were methoxy and hydroxyl groups; RT was a flavonol, and its distinguishing radical pharmacophores was rutinoside group. In Figure 1, circles mark the distinguishing radical pharmacophores of the compounds.