The S. miltiorrhiza seedlings used in the study were purchased from the plantation of S. miltiorrhiza in Pingyi County, Shandong Province (35°30' N, 117°35' E). The area has a temperate seasonal climate with an average elevation of 87.9 m. The annual mean precipitation, air temperature, and average relative humidity are 836.0 mm, 14.3 °C and 67.1%, respectively (Lu, 2020; Tian et al., 2021). Healthy seedlings with approximately 4 basal leaves were cultivated from the upper and middle parts of the annual root. They were identified as S. miltiorrhiza Bge. By Associate professor Xiaoyun Wang from the Research Center for Traditional Chinese Medicine Resources and Ethnic Minority Medicine of Jiangxi University of Chinese Medicine (JXUCM).

Healthy and disease-free seedlings were selected and planted in Shennong garden of JXUCM. The physical and chemical properties of the topsoil (0–20 cm) in the planting area of S. miltiorrhiza in Shennong garden were as follows: total nitrogen 0.30 g kg⁻¹, total phosphorus 0.28 g kg⁻¹, total potassium 27.15 g kg⁻¹, available nitrogen 0.01 g kg⁻¹, available phosphorus 0.01 g kg⁻¹, available potassium 0.08 g kg⁻¹, organic matter 1.97 g kg⁻¹, pH 4.6, and total Cd 0.92 mg kg⁻¹. The total Cd content was lower than the critical value of Cd (1.0 mg kg⁻¹) in soils, which could ensure the normal growth of plants.

The indicated amounts of cadmium chloride hemi (pentahydrate) (CdCl₂·2.5H₂O) were mixed well with sieved topsoils from Shennong garden. Based on the study of Zhang et al. (2013) and Wei et al. (2020), four Cd treatment levels were set in the present study as follows: 0 (CK), 25 mg kg⁻¹ (T1), 50 mg kg⁻¹ (T2), and 100 mg kg⁻¹ (T3). Subsequently, 2 kg of the soils were placed a flower pot (16 × 17 cm) and incubated for 30 d. Then, the S. miltiorrhiza seedlings, which had grown in Shennong garden for 30 d, were transplanted into the pots, and each pot included one seedling. There were three repetitions at each level and three seedlings in each repetition. After 15 d of treatment, root samples were collected, washed with ultrapure water, quickly frozen with liquid nitrogen, and transferred to −80°C until further metabolomic analysis. Voucher specimens (No. DS-001) were deposited in a public herbarium in the Research Center for Traditional Chinese Medicine Resources and Ethnic Minority Medicine of JXUCM.

A 25 mg fresh sample in a 500 μL mixture of methanol and water (3:1, v/v) (including an isotope-labeled internal standard mixture) was ground at 35 Hz for 4 min and lysed in an ultrasonic water bath for 5 min. After sitting at −40°C for 1 h, the samples were centrifuged at 12000 rpm for 15 min at 4°C. The supernatant was stored at −80°C until liquid chromatography-tandem mass spectrometry (LC–MS/MS) analysis. Quality control (QC) samples were prepared with a mix of the supernatants from all samples.

Ultrahigh-performance liquid chromatography and quadrupole orbital well hybrid TM mass spectrometry coupled with an Acquity UPLC HSS T3 liquid chromatographic column were conducted to separate the target compounds.

The chromatographic conditions were as follows: an Acquity UPLC HSS T3 liquid chromatographic column (2.1 mm × 100 mm, 1.8 μm) was used; the column temperature was 35°C; mobile phase A and phase B were the aqueous phase (containing 5 mmol L⁻¹ ammonium acetate and 5 mmol L⁻¹ acetic acid) and acetonitrile, respectively; gradient elution was carried out (0–0.7 min, 1% B; 0.7–9.5 min, 1–99% B; 9.5–11.8 min, 99% B; 11.8–12.0 min, 99–1% B; 12–14.8 min, 1% B); the flow rate was 0.5 mL min⁻¹ and the injection volume was 3 μL.

A Thermo Q Exactive HFX mass spectrometer was used to collect the primary and secondary mass spectrometry data of the samples under the control of Xcalibur 4.0.27 (Thermo Scientific). The detailed parameters are presented in Supplementary Table S1.

The raw data were converted into mzXML format with the software ProteoWizard (https://proteowizard.sourceforge.io/), and the recognition, extraction, alignment and integration of the peaks was conducted with the R package (XCMS as the core). Then, the MS2 database was applied for metabolite annotation. The cutoff for annotation was set at 0.3. Data with no definite substance name and no spectral ratio or substances with a missing quantity greater than 50% in the comparison group samples were filtered and removed. For substances with a missing quantity less than 50%, the K-nearest neighbor (KNN) algorithm was used to simulate the missing value. Finally, the total ion current (TIC) or internal standard (IS) of each sample was used to normalize the data, and 305 metabolites were identified for further analysis.

Unsupervised principal component analysis (PCA) was performed to analyze the distribution of root samples in CK, T1, T2 and T3. Supervised orthogonal partial least squares-discriminant analysis (OPLS-DA) was used to distinguish the metabolic profiles of the roots between the CK and each Cd-added group (T1, T2 and T3). To test the OPLS-DA model, a cross-validation residual variance test (CV-ANOVA test) and 200 permutation tests were carried out. The model with p < 0.05 in the CV-ANOVA test or R2 > 0.7 and Q2 > 0.4 in the permutation test was reliable (Yuan et al., 2020). The variable importance in projection (VIP) was vital for explaining the data of the OPLS-DA model. The discriminating metabolites were defined based on a VIP value above 1 and a p value below 0.05 (Yuan et al., 2020). Metabolic pathway analysis of these discriminating metabolites was conducted. Meanwhile, one-way ANOVA was performed to calculate the variability in the relative contents of the metabolites.

Multivariate statistical analysis (PCA and OPLS-DA) and validation of the OPLS-DA model were carried out with SIMCA-P 14.1 (Umetrics, Sweden). Univariate statistical analysis and one-way ANOVA were conducted with IBM SPSS Statistics 20.0 (SPSS, Inc., IBM Corp, New York, USA). The metabolic pathway analysis was carried out with the online software MetaboAnalyst 4.0 (http://www.Metaboanalyst.ca/faces/ModuleView.xhtml). The heatmap was drawn by online software (https://matrix2png.msl.ubc.ca/bin/matrix2png.cgi). The Venn diagrams were drawn by Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

By observing the differences in the peak height of the internal standard between the QC samples and the peak conditions of the internal standard in the blank samples, this study determined whether the instrument was stable and whether there were residues in the detection process. Supplementary Figure S1 shows that both the retention time and the response strength of the internal standard in the QC samples were stable, and the data acquisition stability of the instrument was excellent. Supplementary Figure S2 shows that no obvious peaks other than the internal standards were detected in the blank samples.

The metabolite profiles of the roots were investigated using an untargeted global metabolomic platform with LC–MS/MS. A total of 305 metabolites were annotated and quantified in S. miltiorrhiza roots with different levels of Cd treatment (Supplementary Table S2), which could be categorized into seven major groups based on their molecular structure: amino acids, organic acids, fatty acids, ketones, sugars, amides, and others (Figure 1C; Supplementary Table S3).

There were 136 upregulated metabolites and 169 downregulated metabolites in T1 and 141 upregulated metabolites and 164 downregulated metabolites in both T2 and T3 (Figures 1A,B). Proportions of the metabolites in the roots changed in response to Cd addition (Figure 1C; Supplementary Table S3). Compared to CK, the total content of amino acids increased by 0.2-fold in T1, approximately 0.9-fold in T2 and T3 (p < 0.05), and the total content of fatty acids decreased by 87% in T1, 59% in T2 and 75% in T3 (p < 0.05). There were no significant differences in the total content of organic acids between CK and each Cd-added group (Figure 1; Supplementary Table S3).

According to the results of the PCA, root samples from CK could be separated from the Cd-treated groups (Figure 2). Based on the reliable OPLS-DA model with a p value below 0.05 (Figure 3; Supplementary Figure S3A), S. miltiorrhiza roots in T1 and CK could be distinguished by metabolites with major contributions from amino acids, organic acids, and fatty acids (Figure 4): 174 discriminating metabolites were selected, 56 of which were upregulated and 118 of which were downregulated. Among them, most fatty acids (especially α-linolenic acid, linoleic acid (C18:2), and oleic acid (C18:1)) and organic acids were downregulated, with chorismate and o-tyrosine contributing the most to the upregulation of organic acids and amino acids in S. miltiorrhiza roots, respectively.

Similarly, based on the significant OPLS-DA model (Supplementary Figures S3B,C), the root samples could be significantly distinguished by metabolites with the major contributions from 161 metabolites (81 of them were upregulated and 80 were downregulated) in T2 and CK and from 191 metabolites (120 of them were upregulated and 71 were downregulated) in T3 and CK. Major contributors to the differences were amino acids, organic acids, and fatty acids (Figure 4): compared to CK, in T2, all of these amino acids (especially DL-tryptophan, L-aspartic acid, and L-proline) were upregulated, and most fatty acids (e.g., linoleic acid, oleic acid) were downregulated, with chorismate contributing the most to the upregulation of the organic acids (Figure 4); in T3, most of the amino acids (especially L-proline, L-serine and L-histidine), organic acids (especially eudesmic acid and rosmarinic acid) and all of the saturated fatty acids (e.g., stearic acid (C18:0)) were upregulated, and oleic acid was downregulated (Figure 4).

The discriminating metabolites were annotated into the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and 65 metabolic pathways were obtained (Figure 5A; Supplementary Table S4). Among them, biosynthesis of unsaturated fatty acids, isoquinoline alkaloid biosynthesis, betalain biosynthesis, aminoacyl-tRNA biosynthesis, and tyrosine metabolism were the significantly enriched metabolic pathways (p < 0.05) (Figure 5B; Supplementary Table S4). In addition, a metabolic map of the resistance process was developed based on these results (Figure 6).

Our results clearly revealed that amino acids, organic acids and fatty acids in S. miltiorrhiza roots played essential roles in resisting Cd stress. As the basic components of proteins, amino acids can participate in the regulation of ion transport and nitrogen metabolism and play vital roles in resisting abiotic stress in plants (Sharma and Dietz, 2006; Xu et al., 2012; Tian, 2021). Organic acids are a class of compounds containing carboxyl groups (excluding amino acids), most of which can combine with metal ions or alkaloids (Kuang, 2017). Fatty acids and their derivatives are essential energy stores and the main components of membrane lipids (Chaffai et al., 2009; Liu et al., 2014; Li-Beisson et al., 2016). Due to the important roles of these three classes of metabolites in plant survival and stress resistance (Liu et al., 2014; Panchal et al., 2021; Trovato et al., 2021), in this study, the changes in the proportions of amino acids, organic acids, and fatty acids were larger than those of the other metabolites in S. miltiorrhiza roots under Cd treatment (Figure 1).

In comparison with CK, in T1, T2 and T3, the total relative contents of the amino acids significantly increased (p < 0.05), the total relative contents of the fatty acids significantly decreased (p < 0.05), and the total relative contents of the organic acids presented no significant differences (Figure 1; Supplementary Table S3). In disagreement with these results, in Oryza sativa L. roots, Cd at a low content could promote the production of amino acids and organic acids, and Cd at a high content could inhibit the production of amino acids and organic acids (Tang et al., 2016), while in agreement with our results, there was a significant decline in the production of fatty acids in Sedum plumbizincicola L. roots under Cd stress (Sun et al., 2020). This might result from the fact that the plant metabolic response to abiotic stress is affected by the stress modes, stress intensity, and plant species (Hu and Xu, 2014; Zhang and Chen, 2021).

As reported above, amino acids were the main contributors to the differences in S. miltiorrhiza roots under different levels of Cd stress (Figure 4). Well known as an abiotic and biotic stress indicator, proline can function as a hydroxyl radical scavenger and it plays a vital role in the adjustment to osmotic stresses in plants (Sharma and Dietz, 2006; Zemanov et al., 2017). Tyrosine is the precursor of many metabolites (including tocopherol, plastoquinone and ubiquinone) that are essential to the survival and stress resistance of plants (Kilgore and Kutchan, 2016; Cassels and Sáez-Briones, 2018; Xu et al., 2020). Tryptophan, an aromatic amino acid, plays important roles in the regulation of plant development and it acts as a precursor for the biosynthesis of the hormone auxin (Sanjaya et al., 2008; Liu et al., 2011a). Histidine plays important roles in the regulation of the biosynthesis of other amino acids and in the chelation and transport of metal ions (Stepansky and Leustek, 2006). Aspartic acid can be fed into the synthesis of other amino acids (e.g., lysine, methionine and threonine) (Angelovici et al., 2009). Serine is essential for the synthesis of proteins and other biomolecules, including nucleotides and serine-derived lipids (e.g., phosphatidylserine and sphingolipids), and is involved in the resistance to various stresses (Ho and Saito, 2001; Waditee et al., 2007; Ros et al., 2013). Due to the different roles of these amino acids in the stress resistance of plants, with different levels of Cd treatment, the levels of amino acids playing vital roles in resisting Cd stress were different (Benral and McGrath, 1994; Wang et al., 2015; Cosio and Renault, 2020) (Figures 5, 6).

In the present study, o-tyrosine contributed more to the upregulation of the discriminating amino acids between T1 and CK; DL-tryptophan, L-aspartic acid, and L-proline contributed more to the upregulation of the discriminating amino acids between T2 and CK, and L-proline, L-serine and L-histidine contributed more to the upregulation of the discriminating amino acids between T3 and CK (Figure 4; Supplementary Table S2). This illustrated that there were different Cd stress resistance patterns of amino acids in S. miltiorrhiza roots under different levels of Cd stress. Under moderate Cd (T2) stress, S. miltiorrhiza roots might mainly increase the relative contents of DL-tryptophan, L-aspartic acid, and L-proline to resist Cd, but under high-level Cd (T3) stress, the roots might mainly increase the relative contents of L-proline, L-serine and L-histidine to resist Cd. Inconsistent with this, compared to CK, the proline content declined significantly and the aspartic acid content increased significantly with moderate Cd treatment; and the histidine content increased significantly and the serine content declined significantly with high Cd treatment (p < 0.05) of Crassocephalum crepidioides (Zhu et al., 2018). As the level of Cd increased, the contents of tryptophan and proline increased significantly, and the L-aspartic acid content declined significantly in Solanum nigrum (Xu et al., 2012). These differences might be caused by the differences in the Cd sensitivity of different plants or by the different methods of applying Cd stress (Zhang and Chen, 2021).

Our study revealed that based on the changes in the relative contents of organic acids in S. miltiorrhiza roots with different levels of Cd treatment, chorismate and eudesmic acid in T1, chorismate in T2, and eudesmic acid and rosmarinic acid in T3 were upregulated. In the present study, organic acids were among the main contributors in distinguishing CK samples from Cd treated samples (T1, T2 and T3) (Figure 4; Supplementary Table S2): most organic acids showed a decreasing accumulation, with chorismate and eudesmic acid contributing the most to the upregulation of organic acids in T1; most organic acids showed an increasing accumulation, with chorismate contributing most to the upregulation of organic acids in T2; and most organic acids presented an increasing accumulation (especially eudesmic acid and rosmarinic acid) in T3.

Chorismate could function as a key branching point between primary and secondary metabolisms, as well as being a precursor of aromatic amino acids (e.g., tryptophan, phenylalanine, tyrosine) and hormonal substances (e.g., indoleacetic acid and salicylic acid) that play essential roles in plant metabolism (Kristin and Michael, 2010). Eudesmic acid and rosmarinic acid are polyphenol derivatives and phenol compounds, respectively, and both have strong antioxidation activity and strong hydroxyl radical scavenging activity (Lv et al., 2005; Wang et al., 2014; Korkmaz et al., 2018). Under Cd stress, the elevation of intracellular chorismite would result in some pathways, including the biosynthesis of aromatic amino acids (e.g., phenylalanine and tryptophan) and secondary metabolites (e.g., polyphenols and flavonoids) (Figure 6; Supplementary Table S2) (Malik, 1979), and the activity of the enzymes involved in phenolic compound metabolism might increase (Michalak, 2006). This might lead to an increase in phenolic compounds and polyphenol derivatives, which could bind heavy metals, enhance the activity of antioxidant enzymes, and reduce the harmful effects of heavy metals on plants (Güez et al., 2017; Manquián-Cerdaa et al., 2018) (Figure 6). Therefore, in response to low-level (T1) and high-level (T3) Cd treatment, S. miltiorrhiza roots upregulate their secondary metabolism, mainly by upregulating eudesmic acid and rosmarinic acid to resist Cd stress (Figure 6) (Zoufan et al., 2020).

Our study proved that fatty acids played vital roles in resisting Cd stress. Compared to CK, most discriminating fatty acids (e.g., oleic acid and linoleic acid) decreased significantly both in T1 and T2, and unsaturated fatty acids (e.g., oleic acid and linoleic acid) decreased significantly, with all of the discriminating saturated fatty acids (e.g., stearic acid) being upregulated (Figure 5; Supplementary Table S2) in T3. As discussed above, in plants, fatty acids and their derivatives play vital roles in improving stress tolerance by participating in various defense pathways, including basal, systemic, and effector-triggered immunity (Chaffai et al., 2009; Liu et al., 2014; Li-Beisson et al., 2016). Unsaturated fatty acids (e.g., oleic acid and linoleic acid) are produced by the catalysis of fatty acid desaturases (FADs), which can catalyze the formation of double bonds at specific positions in the chain of saturated fatty acids (e.g., stearic acid) to regulate the response to various stresses (Liu et al., 2011b; Park et al., 2015). It is important to regulate the lipid composition and adjust the unsaturation level of membrane fatty acids to cope with metal stress (Thompson, 1992) (Figures 4, 6; Supplementary Table S2).

Therefore, oleic acid and linoleic acid, which could induce the production of plant reactive oxygen species, showed a decreased accumulation in T1, T2 and T3 to maintain the normal growth of S. miltiorrhiza under Cd stress (Cury-Boaventura and Curi, 2005) (Figure 5), and high-level Cd stress might inhibit the expression of FADs, leading to the enrichment of saturated fatty acids (especially stearic acid) (Figure 5) (Liu et al., 2013). In contrast to our study, other researchers found that Cd stress increased the contents of oleic acid, linoleic acid, and stearic acid in Asterioneilu glacialis (Jones et al., 1987) and decreased the linoleic acid content in Sedum plumbizincicola (Sun et al., 2020) and the oleic acid content in Sedum alfredii (Luo et al., 2014). As mentioned above, these differences might be caused by the differences in the Cd sensitivities of different plants (Zhang and Chen, 2021).