Salvia miltiorrhiza was cultivated at 25°C with a 16-h-light/8-h-dark cycle in MS basal medium containing 3% sugar and 8% agar. CuSO₄·5H₂O (Sinopharm Chemical Reagent Company., Ltd, China) was dissolved in distilled water to make concentrations of 5, 25, and 100 mM. All the prepared CuSO₄·5H₂O were sterilized through 0.22-μm filters and then added to MS basal medium to a final concentration of 5, 25, or 100 μM. S. miltiorrhiza was treated with different concentrations of CuSO₄·5H₂O for 5 days, and then plant matter was collected for the following experiments.

After 5 days of treatment with different concentrations of CuSO4·5H₂O, a photo of S. miltiorrhiza planted in a culture flask was taken. To measure fresh weight, the whole plant of S. miltiorrhiza was taken out from the culture flask, washing away culture medicum, and measured with electronic balance. To accurately reflect fresh weight changes, at least four plants were used each time, and average weight was calculated.

To extract phenolic acids, fresh roots of copper-treated or non-treated S. miltiorrhiza were collected from a culture flask, dried at 55°C until the weight remained constant, and then grounded into powder. Each sample of 0.1 g was added to 1-ml methanol: water (70: 30, v v-1), followed by 1-h ultrasound and centrifugation at 12,000 rpm for 20 min. The supernatant was then removed and filtered through a 0.22-μm filter. HPLC was performed on an Agilent 1290 with a DAD detector (Agilent Technologies, USA), equipped with an Agilent ODS-SP 5-μm 4.6 mm × 250 mm column. The mobile phase was H₂O (0.1% formic acid, A) and acetonitrile (B). The HPLC program was set to the following line gradient: 10–25% of B for 0–10 min, 25–27% of B for 10–15 min, 27–55% of B for 15–20 min, 55–95% of B for 20–21 min, held 5 min, 95–10% of B for 26–27 min, and 10% of B for 27–30 min. The flow rate was 1 ml/min, the column temperature was set to 35°C, sample injection volume was 10 μl, and the wavelength was 280 nm. Sal B and RA standards (Shanghai Yuanye Bio-Technology Co., Ltd, China) were configured as a mixed standard solution with a concentration of 0.25 mg/ml. The standard curve method was used to calculate the content of Sal B and RA. The experiment was performed with four repetitions, and data were expressed as mean ± SD.

For metabolomics analysis, the same samples were analyzed by LC-MS using an Agilent 6545A Q/TOF mass spectrometer equipped with an electrospray ion (ESI) source. In the negative ion mode, the following conditions were used: a drying N2 gas flow rate, 8 L/min; temperature, 320°C; nebulizer, 35 psig; capillary, 3,000 V; skimmer, 65 V; OCT RFV, 750 V. Mass spectra were recorded in a full-scan mode from 100 to 1,200 m/z range. Agilent 1290 Infinity Ultra-performance liquid chromatography (UPLC, Agilent technologies, USA), equipped with an Agilent ODS-SP 5 μm 4.6 mm × 250 mm column, was used to affect the separations. The mobile phase consisted of acetonitrile and 0.03% acetic acid (mobile phase A) and H2O with 0.03% acetic acid (mobile phase B), which were eluted as follows: 10–13% of B at 0–8 min, 13–23% of B at 8–13 min, 23% of B at 25 min, 23–27% of B at 25–45 min, 27–36% of B at 45–54 min, 36–70% of B at 54–65 min, 70–73% of B at 65–75 min, 73–93% of B at 75–85 min, 93–100% of B at 85–87 min, and then held for 5 min. The flow rate was 0.5 ml/min; the column temperature was set at 35°C, sample injection volume was 5 μl, and post-run time was 12 min.

For data analysis, raw LC-MS data were converted to mzData format using DA reprocessor software (Agilent), and the XCMS package in R was used to perform peak finding, filtering, and alignment. Metabioananlyst 4.0 (https://www.metaboanalyst.ca) was used to obtain normalized data. Metabolites with P < 0.05 and fold change (copper treatment/blank) >1.5 or < 0.5 were regarded as differential metabolites. Metlin database, Human Metabolome database, and our in-house laboratory database were used to identify metabolites. The experiment was performed with four repetitions, and data were expressed as mean ± SD.

Total RNA was isolated from fresh root of copper-treated or non-treated S. miltiorrhiza using RNeasy Plus Kit with the genomic DNA removal step. The concentration and the quality of extracted RNA were evaluated. cDNA library construction and sequencing were performed by the Biomarker Technologies Corporation (Beijing, China). The raw sequence data were filtered to obtain clean data, which were then compared with the S. miltiorrhiza reference genome using HISAT2 software (https://ccb.jhu.edu/software/hisat2/manual.shtml). Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) was used to calculate gene expression. Differential gene expression analysis used edgeR (Robinson et al., 2010), and genes with p < 0.05 and fold change >1.5 or <0.5 are regarded as differentially expressed genes (DEGs). The DEGs were classified into functional categories by blasting against the clusters of orthologous groups against the eukaryotic complete genomes (KOG) database (http://www.ncbi.nlm.gov/KOG) and the gene ontology (GO) database (http://geneontology.org/). The experiment was performed with three repetitions, and data were expressed as mean ± SD.

Total RNA extraction was achieved using a plant total RNA extraction kit (Tiangene Biotech Co., Ltd., China). Briefly, a 0.5-g fresh S. miltiorrhiza root was grounded into powder with liquid nitrogen, and extraction was carried out according to the instructions. RNA reverse transcription (Hieff™ First Strand cDNA Synthesis Super Mix, Yeasen, China) and gene quantification (Hieff™ qPCR SYBR Green Master Mix kit, Yeasen, China) were performed following the protocols of the manufacturer. The qPCR conditions were as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 10 s, 60°C for 30 s, and 70°C for 90 s. mRNA expression levels were calculated using the ^2−ΔΔCt method and presented as a ratio to β-actin. Primer sequences are shown in Supplementary Table 4.

A Laccase Activity Detection kit (Solarbio Life Sciences, Beijing, China) was used to detect laccase activity. To detect laccase activity, a 0.1-g fresh S. miltiorrhiza root was added to a 1-ml extraction buffer, the mixture was homogenized on ice, and then centrifuged at 12,000 rpm for 20 min at 4°C. The supernatant was then moved to a new centrifuge tube, and the absorbance was measured at 562 nm to determine protein concentration of supernatant using a BCA kit (Good Laboratory Practice Bioscience, Montclair, CA, USA). Subsequently, 150 μL of supernatant was mixed with 850 μL of a working buffer and heated in an oven at 45°C for 3 min. As negative control, 150 μL of an extraction buffer was mixed with 850 μL of the working buffer and heated in an oven at 45°C for 3 min. Then, the produced ABTS radical was determined through measuring mixture absorbance at 420 nm. Enzyme activity was calculated through measuring ABTS radical content change according to the instructions of the manufacturer.

The whole genome sequence and protein data from S. miltiorrhiza (BioProject: PRJNA287594) were downloaded from the Genome Warehouse in the BIG Data Center (http://gigadb.org/dataset/100164). The whole genome data and protein data of Arabidopsis thaliana were downloaded from the Arabidopsis thaliana information resource website (TAIR, https://www.arabidopsis.org/). In order to select all the laccases in the genomes of Salvia miltiorrhiza and Arabidopsis thaliana, we downloaded the unique HMM model (the hidden Markov model) of laccase from the Pfam database. HMMER 3.0 software with hmmsearch was used to search for laccase-conserved domain models, including PF00394, PF07731, and PF07732, in the protein database of S. miltiorrhiza and Arabidopsis thaliana. A phylogenetic tree was constructed using the maximum likelihood (ML) method with bootstraps in MEGA X software.

FPKM values of genes in copper ion-treated and non-treated groups were used to calculate correlation coefficient. Gene co-expression analysis between the candidate transcription factors and the key genes was performed by the mean of the Pearson correlation test using IBM SPSS Statistics 20. Data are shown as the mean ± SD (n = 3) (^*p < 0.05 and ^**p < 0.01). The correlation coefficient >0.8 was considered to be co-expressed, and the coexpression network was visualized using Cytoscape.

The promoter sequences of analyzed genes (from 2,000 bp, relative to a translation start site) were extracted using gff file of S. miltiorrhiza genome. Core motif of copper-response element (GTAC) (Quinn et al., 2000; Kropat et al., 2005) and metal response-element (TGCxCxC) (Murphy et al., 1999) were predicted in extracted 2,000 bp promoter sequences of laccases and transcription factors using the Plantcare database (Lescot et al., 2002) and visualized by TB tools (Chen et al., 2020). In addition, transcription factors-binding motifs were also predicated in extracted 2,000 bp promoter sequences of laccases.

To explore the effect of copper on the growth of S. miltiorrhiza plants, we exposed the plants to predetermined copper concentrations (0, 5, 25, and 100 μM) and determined the phenotypic changes thereafter. After 5 days of cultivation, the S. miltiorrhiza plants wilted under treatments with a greater copper concentration. Compared with the 0 and 5 μM treatments, S. miltiorrhiza plants grew slowly and the leaves withered and yellowed in the 25- and 100-μM treatments (Figure 1A). It is noteworthy that low concentration of copper (5 μM) significantly promoted plant growth compared with the blank (0 μM) group (Figure 1A). Thus, 5 μM of copper was optimum for the growth of the S. miltiorrhiza plants.

As the major SalAs, Sal B, and RA are important active secondary metabolites in S. miltiorrhiza. To ascertain the most suitable concentration of copper that promotes SalAs biosyntheses, we used HPLC to determine the contents of Sal B and RA in the copper-stressed groups and the blank group. The low concentration of copper treatment increased Sal B by an average of 4.27-fold compared with the control, but the medium and high concentrations of copper treatments decreased Sal B by 2.16- and 2.06-fold, respectively, compared with the control (Figure 1B and Supplementary Figure 1). Compared with the blank group, low concentration of copper increased RA content 1.55-fold while medium and high concentrations of copper decreased its amount 4.22- and 2.37-fold, respectively (Figure 1B and Supplementary Figure 1). Based on these findings, it is evident that low concentration of copper promotes the biosyntheses of Sal B and RA.

To further confirm the holistic effect of low copper on the accumulation of SalAs in S. miltiorrhiza, the plant samples exposed to low concentration of copper and the blank group were analyzed using UPLC-QTOF-MS. In the negative iron model, the total ion chromatogram was different between the control and the low copper treatment (Supplementary Figure 2). A total of 2,041 ions were detected (Figure 2A). Using the selection criteria of change ratio >1.5 (significant increase; P < 0.05) or < 0.5 (significant decrease; P < 0.05), 90 increased and 70 decreased compounds were detected. The principal component analysis (PCA) analysis showed that the copper stress treatment and control clearly clustered into two categories (Figure 2B). Among the different metabolites, 13 SalAs, including Sal B, RA, lithospermic acid, salvianolic acid D, salvianolic acid A, and danshensu, were significantly increased in the low concentration of copper treatment compared with the control (Figure 3A), indicating that low copper stress can positively regulate the synthesis of SalAs. The pathway enrichment analysis of differential metabolites showed that copper ions had significant effects on amino acid synthesis and metabolism, the pentose phosphate pathway, and carbon fixation in photosynthesis in S. miltiorrhiza (Figure 3B).

To explore how the low concentration of copper treatment regulated the biosynthesis of SalAs in S. miltiorrhiza, a transcriptome analysis was performed between the control and low copper treatments. Using the selection criteria of the change ratio P < 0.05 and fold change >1.5 and <0.5, a total of 418 differential genes were obtained, of which 207 were upregulated and 211 were downregulated in the low copper-treated group compared with the blank group (Supplementary Table 2). Gene Ontology (GO) analysis showed that these significantly changed genes were mainly involved in copper ion binding, oxidoreductase, and ROS metabolism (Figure 4A). Pathway enrichment analysis showed that the differentially expressed genes are related to tyrosine and phenylalanine metabolism, chlorophyll biosynthesis, plant pathogene interaction, photosynthesis, and plant hormone signal transduction (Figure 4B).

A further gene expression pattern indicated that SalAs biosynthesis-related genes (laccase 20, laccase 13, laccase 7, and rosmarinate synthase), photosynthesis-related genes (blue copper protein, basic blue protein, and chlorophyll A-B-binding protein), and zinc transporter 5 were upregulated in the copper-treated group compared with the blank group (Figure 5A and Table 1). In contrast, copper transport-related genes, such as copper transporter 6, copper transporter 5, copper transport protein CCH, and copper-transporting ATPase RAN1, were downregulated in the copper-treated group (Figure 5A and Table 1). In addition, transcription factors, such as ethylene-responsive transcription factor 5, ERF012, ethylene-responsive transcription factor RAP2-7, ethylene-responsive transcription factor TINY, transcription factor MYB14, MYB33, MYB73, basic helix-loop-helix transcription factor, and transcription factor bHLH69, were responsive to low copper treatment (Figure 5A and Table 1).

To determine whether the gene expression changes corresponded to physiological changes, chlorophyll content and copper content in S. miltiorrhiza were measured. Chlorophyll, which plays an extremely important role in photosynthesis, is the main component of chloroplasts in leaves. After the copper ion treatment, we detected the difference in chlorophyll content between the 5-μM copper treatment and the control (Figure 5B). The chlorophyll a, chlorophyll b, and total chlorophyll of S. miltiorrhiza in the 5-μM copper treatment were 1.08, 0.31, and 1.39 mg/g, respectively. The chlorophyll a, chlorophyll b, and total chlorophyll of S. miltiorrhiza in the control were 0.49, 0.03, and 0.52 mg/g, respectively. Compared with the control, the 5-μM copper treatment increased chlorophyll a, chlorophyll b, and total chlorophyll by 2.19, 9.98, and 2.66 times, respectively. Copper content in the treatment group was also significantly increased compared with the control (Figure 5C). The greater chlorophyll content increased the photosynthetic ability of the copper treatment group, the growth state of the plant was healthier, the ability of the plant to resist adverse environments was increased, and the treatment was conducive to the accumulation of active products in S. miltiorrhiza.

Based on the transcriptomics result, we verified the mRNA expression of laccase 20, laccase 13, and laccase 7 using RT-qPCR. The results obtained showed 2.3 and 1.4-fold increases in laccase 20 and laccase 13 expressions, respectively, in the low-copper group relative to the control group. Consequently, laccase enzyme activity increased 3.19-fold in the low-copper group compared with the control group (Figure 6A). These findings suggest the pivotal role of laccase in the synthesis of salvianolic acids in the presence of low copper. To explore the evolutionary relationships among laccases, the sequences of laccases were extracted from S. miltiorrhiza and A. thaliana, and a phylogenetic tree was constructed using the ML method with bootstraps in MEGA X software (Figure 6B). Laccases in S. miltiorrhiza and Arabidopsis thaliana were grouped into seven clusters (I, II, III, IV, V, VI, and VII) based on similarity of amino acid sequence. The copper-induced laccases (marked with a red pentagram) were positioned in Cluster II and Cluster V (Figure 6B). Functional domain analysis found that the three key laccases had the same three motifs and conserved domains, which implies that they have functional similarity (Figure 6C).

To identify the upstream transcription factor involved in the regulation of copper-induced biosynthetic enzymes of SalAs, a coexpression analysis of biosynthetic genes related to SalAs (laccases and rosmarinate synthase) and 12 key transcription factors was carried out. As shown in Figure 7, transcription factors were positively correlated with key candidate genes in the biosynthetic pathways of SalAs. Among them, zinc finger protein ZAT10, Myb family transcription factor PHL5, ethylene-responsive transcription factor 5, and bZIP transcription factor 18 were significantly coexpressed with laccase 20 (gene 3); Zinc finger CCCH domain-containing protein 33/29, zinc finger A20 and AN1 domain-containing stress-associated protein 5, transcription factor MYB73, and ethylene-responsive transcription factor ERF012 were significantly coexpressed with laccase 7 (gene 1); Transcription factor MYB33 was significantly coexpressed with laccase 13 (gene 3) (Supplementary Table 3). To further explore the underlying mechanism by which copper induces expression of laccases and transcription factors, a copper-response element and metal-response elements were searched against 2,000 bp promoter sequences of laccases and transcription factors. As shown in Supplementary Figure 3, a core motif of a copper-response element (GTAC) (Quinn et al., 2000; Kropat et al., 2005) and a metal-response element (TGCxCxC) (Murphy et al., 1999) were identified in extracted 2,000 bp promoter sequences of laccases (laccase 20 and laccase 13) and transcription factors (zinc finger protein ZAT10, Zinc finger CCCH domain-containing protein 33, transcription factor MYB73, and transcription factor MYB33). In addition, MYB transcription factor-binding motifs were also identified in promoter sequences of laccase 20 and laccase 13 (Supplementary Figure 3).

Our results revealed that the growth status and fresh weight of S. miltiorrhiza were significantly improved in the low copper concentration treatment compared with the control. The metabolites involved in carbon fixation in photosynthesis were significantly increased by the low-copper concentration treatment. Importantly, chlorophyll a, chlorophyll b, and total chlorophyll, which play an important role in photosynthesis (Baker, 2008), showed an increase in accumulation in low copper-treated plants. Photosynthesis contributes greatly to the biosynthesis of secondary metabolites in plants (Landi et al., 2020). A previous study has shown that copper played a role in photosynthesis through strengthening the photosynthetic capacity of plants (Douglas and Jaqueline, 2019). In our study, blue copper protein and basic blue protein, which are involved in promoting photosynthesis through acting as electron acceptors (Skyes, 1991), were upregulated in low copper-treated plants. Sal B and RA, the two most representative secondary metabolites of phenolic acids in S. miltiorrhiza (Zhao et al., 2015; Shi et al., 2016, 2019; Zhou et al., 2017; Huang et al., 2019), were significantly increased in S. miltiorrhiza treated with a low concentration of copper. These results suggest that low concentrations of copper induce the accumulation of SalAs might be through regulating photosynthesis related genes such as blue copper proteins in S. miltiorrhiza.

Both tyrosine- and phenylalanine-derived pathways are the starting point for the biosynthesis of SalAs s (Hao et al., 2020). The data obtained in the transcriptomic analysis showed that rosmarinate synthase and laccase were upregulated in plants treated with a low copper concentration compared with the control. Rosmarinate synthase is the key enzyme in the pathway of SalAs biosynthesis in S. miltiorrhiza (Zhou et al., 2018). A previous study has reported that phenolic acids, such as rosmarinic acid and lithospermic acid B, were decreased in rosmarinate synthase-knockdown hairy root lines of S. miltiorrhiza (Zhou et al., 2018). However, overexpression of rosmarinate synthase resulted in a higher content of phenolic acids (up to over 3-fold) in transgenic lines compared with the control (Fu et al., 2020). Recently, scientists have demonstrated that laccase plays a critical role in the biosynthesis of Sal B. The laccase gene family in S. bowleyana has undergone expansion, and it might catalyze the oxidative reaction from RA to Sal B (Zheng et al., 2021). At present, 29 laccase candidates have been found in S. miltiorrhiza, and they contained three signature Cu-oxidase domains (Li et al., 2019b). It is known that laccases are copper-containing enzymes, and their catalytic activity requires the presence of copper ions (Arregui et al., 2019). Here, a core motif of a copper-response element (GTAC) and a metal response-element (TGCxCxC) were identified in the promoter region of laccase 20 and laccase 13. Overexpression and silencing assays demonstrated SmLAC7 and SmLAC20 played an important role in the biosynthesis of Sal B (Li et al., 2019b). In our study, laccase 20, laccase 13, laccase 7, and rosmarinate synthase were significantly upregulated in S. miltiorrhiza in the low-copper treatment compared with the control, suggesting that treatment with a low concentration of copper may promote the accumulation of SalAs via the regulation of laccase 20/13 and rosmarinate synthase.

In addition to the upregulation of laccase 20/13 and rosmarinate synthase, transcription factors, such as ethylene-responsive transcription factor (ERF), transcription factor MYBs, and basic helix-loop-helix transcription factor (bHLH), were also responsive to the low-copper treatment. The content of SalAs was increased in SmERF115-overexpressing hairy roots and was decreased in silencing lines (Sun et al., 2019). Previous studies have reported that SmMYB1 was a positive activator that improved the accumulation of phenolic acids in S. miltiorrhiza (Deng et al., 2020b; Zhou et al., 2021). However, SmMYB36 was reported to inhibit the accumulation of phenolic acids (Ding et al., 2017). For the bHLH transcription factor, the overexpression of SmbHLH37 substantially decreased yields of Sal B (Du et al., 2018). Here, zinc finger protein ZAT10, Myb family transcription factor PHL5, ethylene-responsive transcription factor 5, and bZIP transcription factor 18 were significantly upregulated and coexpressed with laccase 20 under low concentration of copper treatment in S. miltiorrhiza. In our study, a core motif of a copper-response element (GTAC) and a metal-response element (TGCxCxC) were identified in the promoter region of transcription factors (zinc finger protein ZAT10, Zinc finger CCCH domain-containing protein 33, transcription factor MYB73, and transcription factor MYB33). In addition, MYB transcription factor-binding motifs were identified in promoter sequences of laccase 20 and laccase 13. In Arabidopsis, copper has been shown to induce expression changes of transcription factors, including bHLH, MYB, and ERF (Jakubowicz et al., 2010; Perea-García et al., 2010; Cai et al., 2021). Combined, these results suggest that transcription factors, such as ERF, MYB, and bHLH, play a role in copper-induced SalA biosynthesis in S. miltiorrhiza.

We systematically studied the effects of copper on SalA biosynthesis in S. miltiorrhiza. We identified a suitable copper ion concentration to promote the biosynthesis of SalA and provide a potential regulatory mechanism by identifying a series of copper-responsive genes. Our work provides an effective approach to induce SalAs and demonstrates a promising future for the metabolic regulation of S. miltiorrhiza.