Dynamic transcriptomics unveils parallel transcriptional regulation in artemisinin and phenylpropanoid biosynthesis pathways under cold stress in Artemisia annua

Dynamic changes of artemisinin and related metabolites under cold stress

The biosynthesis of terpenoids commences with the synthesis of two 5-carbon precursors: isopentenyl diphosphate (IPP) and its isomer, dimethylallyl diphosphate (DMAPP). These precursors are derived through the cytosolic mevalonate (MVA) pathway and the plastidial methylerythritol phosphate (MEP) pathway25. In the subsequent stages of the terpenoid biosynthetic cascade, key compounds such as artemisinic acid, its reduced form dihydroartemisinic acid (DHAA), and arteannuin B are generated. In our investigation, we conducted a comprehensive metabolite profiling of A. annua, focusing on artemisinin and its associated metabolites. The plant specimens were subjected to low-temperature stress at 4 °C and subsequently analyzed at intervals of 6 h, 2 days, and 7 days using high-performance liquid chromatography (HPLC). Standard curves of four metabolites were shown in Supplementary Fig. 1. In the specimens subjected to this treatment for 6 h, artemisinin levels were observed to range from 0.41 to 0.49 mg/g fresh weight (FW), representing a substantial increase of 4.21 to 5.07 times compared with the control (Fig. 1A). However, this enhancement in artemisinin content was not sustained, as levels were observed to taper off at CD2 and CD7 (Fig. 1A). DHAA, an immediate precursor in the biosynthetic pathway to artemisinin, was likewise identified in our analysis. HPLC analysis demonstrated significant increases in DHAA levels at CH6 and CD7, which were 1.81-fold and 2.16-fold higher compared with the control, respectively (Fig. 1B). The contents of arteannuin B and artemisinic acid, metabolites in the competitively branch pathway of artemisinin biosynthesis, were also tested by HPLC. We observed a significant reduction in the content of arteannuin B in all groups subjected to cold treatment (Fig. 1C). In contrast, the levels of artemisinic acid, a direct precursor to arteannuin B, exhibited a significant increase at CD7, reaching 13.67-fold that of the control group. (Fig. 1D). The artemisinic acid content ranged from 0.25 to 0.27 mg/g fresh weight (FW), representing a substantial increase of 13.16 to 14.09 times over that of the control group (Fig. 1D). The findings indicate that A. annua exhibits an initial surge in artemisinin biosynthesis in response to early low-temperature stress. However, under prolonged cold conditions, the synthesis of artemisinin appears to be attenuated, suggesting a limitation in the biosynthetic process.

Fig. 1
figure 1

The contents of artemisinin and related metabolites in A. annua. (AD) The content of artemisinin (A), dihydroartemisinic acid (DHHA) (B), arteannuin B (C), and artemisinic acid (D) * and ** represent significant differences at P < 0.05 and 0.001, respectively. The population size was n = 6.

Transcriptome sequencing and KEGG enrichment analysis

Transcriptome sequencing was conducted on control and cold-stressed A. annua samples. Furthermore, the feasibility of the transcriptomic data was corroborated by quantitative polymerase chain reaction (qPCR) assays in our prior study, thereby validating the expression profiles obtained26. Correlation heatmaps of the transcriptomic data, coupled with principal component analysis (PCA), substantiated the feasibility of both inter- and intra-group sample comparisons (Supplementary Fig. 2A, 2B). Following the identification of differentially expressed genes (DEGs; |Fold change|>2, FDR < 0.05; Supplementary file 1), we conducted KEGG pathway enrichment analysis. In the KEGG enrichment analysis, the ‘Phenylpropanoid biosynthesis’ and ‘Terpenoid backbone biosynthesis’ pathway were significantly enriched at all cold treatment time points (Fig. 2; Supplementary Fig. 3A, 3B). Furthermore, the ‘Flavonoid biosynthesis’ pathway showed significant enrichment at CH6 and CD7 (Fig. 2, Supplementary Fig. 3B). Additionally, other secondary metabolic pathways were also found to be enriched under cold stress, including ‘Stilbenoid, diarylheptanoid and gingerol biosynthesis’ (Fig. 2) and ‘Ubiquinone and other terpenoid-quinone biosynthesis’ (Supplementary Fig. 3A).

Fig. 2
figure 2

KEGG enrichment analysis of the DEGs at CH6 in A. annua leaves.

The gene expression patterns across the artemisinin, jasmonic acid and phenylpropanoid biosynthetic pathways

The biosynthesis of artemisinin consumes intermediates of the glycolytic cycle, such as glycerol-3-phosphate and acetyl-CoA. According to the KEGG database, we identified genes associated with the artemisinin synthesis pathway that exhibited significant expression changes (|Fold change| > 2, FDR < 0.05). The anabolic process of artemisinin can be compartmentalized into three distinct phases: the mevalonate (MVA) pathway, the methylerythritol phosphate (MEP) pathway, and the artemisinin synthesis pathway. In the MVA pathway, pyruvate and glycerol-3-phosphate are converted into isopentenyl diphosphate through a series of enzymatic reactions. In this study, we observed differential expression of genes encoding key enzymes such as DXS, DXR, HDS, HDR, and IDI under cold stress (Fig. 3A). Specifically, the genes encoding DXS, DXR, HDS, and HDR were significantly downregulated, whereas IDI showed peak expression levels at CD2 (Fig. 3A). The MEP pathway, responsible for the synthesis of the artemisinin precursor isopentenyl diphosphate, utilizes acetyl-CoA. In our research, the majority of differentially expressed genes in this pathway were downregulated under cold stress, including those encoding ACAT, HMGS, and PMK (Fig. 3A). At least, in the artemisinin synthesis pathway, ADS, CPR, and DBR2 were upregulated under cold stress. Notably, the genes encoding ADS and DBR2 exhibited a significant increase in transient expression levels at CH6 (Fig. 3A).

Fig. 3
figure 3

The expression profile of artemisinin and phenylpropanoid biosynthetic genes in A. annua under cold stress. (A) The expression profile of artemisinin biosynthetic genes. (B) The expression profile of phenylpropanoid biosynthetic genes. Significantly induced genes by cold stress (at least one time point) are indicated with red. AACT acetoacetyl-CoA thiolase, HMGS 3-hydroxy-3-methylglutaryl-coenzyme A synthase, HMGR HMG-CoA reductase, MVK mevalonate kinase, PMK phosphomevalonate kinase, PMD diphosphomevalonate decarboxylase, FPS farnesyl pyrophosphate synthase, DXS 1-deoxy-d-xylulose-5-phosphate synthase, DXR 1-deoxy-d-xylulose-5-phosphate reductoisomerase, MCT 2-C-methyl-D-erythritol 4-phosphate cytidyltransferase, CMK 4-diphosphocytidyl-2-C-methyl-d-erythritol kinase, MECPS 2-C-methyl-d-erythritol 2,4-cyclodiphosphate synthase, HDS 1-hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate synthase, HDR 4-hydroxy-3-methylbut-2-enyl diphosphate reductase, IDI isopentenyl-diphosphate Delta-isomerase, HDS 1-hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate synthase, MDS 2 C-methyl-d-erythritol 2,4-cyclodiphosphate synthase, ADS amorpha-4, 11-diene synthase, CYP71AV1 cytochrome P450–dependent hydroxylase, DBR2 double bond reductase 2, ALDH1 aldehyde dehydrogenase 1, CPR cytochrome P450 reductase, PAL phenylalanine ammonia-lyase, 4CL 4-coumarate: CoA ligase, C4H cinnamate 4-hydroxylase, CHS chalcone synthase, CHI chalcone isomerase, F3H flavanone 3-hydroxylase, CYP75B1 flavonoid 3’-hydroxylase, FLS flavonol synthase, DFR bifunctional dihydroflavonol 4-reductase, LDOX leucoanthocyanidin dioxygenase, UGTs coniferyl-alcohol glucosyltransferase, HCT shikimate O-hydroxycinnamoyltransferase, CCoAOMT caffeoyl-CoA O-methyltransferase, CCR cinnamoyl-CoA reductase, CAD cinnamyl-alcohol dehydrogenas.

Lignins and flavonoids are biosynthesized also from glycolytic intermediates through a shared upstream phenylpropanoid pathway, followed by divergent downstream pathways for their respective syntheses. In our study, genes in the upstream phenylpropanoid pathway, including PAL, 4CL, and C4H, were significantly upregulated at CH6 and CD2 (Fig. 3B). Within the downstream lignin biosynthetic pathway, the majority of genes encoding enzymes involved, such as CCR, CAD, CCoAOMT and HCT, were also significantly upregulated at CH6 and CD2 (Fig. 3B). In the flavonoid biosynthetic pathway, key genes including CHS, CHI, F3H, CYP75B1, DFR, LDOX and UGTs exhibited pronounced upregulation (Fig. 3B). Subsequently, we quantified the levels of lignin, total flavonoids, and anthocyanins. Notably, lignin content significantly increased at CH6, CD2 and CD7 (Supplementary Fig. 4A), total flavonoid content exhibited a significant rise at CD7 (Supplementary Fig. 4B), and anthocyanin content markedly ascended at CD2 and CD7 (Supplementary Fig. 4C).

Additionally, the biosynthesis of artemisinin in A. annua is closely associated with the levels of jasmonic acid (JA). Consequently, we analyzed the expression patterns of key genes involved in JA biosynthesis under cold stress. In our study, genes encoding AOS, AOC, OPR, and fadA were significantly upregulated at CH6 (Supplementary Fig. 5A). However, the expression levels of the majority of genes encoding LOX, which is upstream in the JA biosynthetic pathway, were significantly downregulated at CD2 and CD7 (Supplementary Fig. 5A). Interestingly, the quantification of jasmonic acid revealed a biphasic response under cold treatment, with an initial increase followed by a decline, reaching a peak at CH6 that was 1.81-fold higher than that of the control group (Supplementary Fig. 5B).

PLS-DA analysis of secondary metabolites associated with artemisinin and phenylpropanoid biosynthesis

The biosynthesis of secondary metabolites such as terpenoids, flavonoids, and lignin in plants requires a substantial carbon source as precursors. Consequently, we assessed the total organic carbon content in A. annua leaves subjected to cold treatment. Our results indicate that the levels of total organic carbon (TOC) significantly increased at CD2 and CD7, peaking at CD2 to 1.38 times that of the control group (Supplementary Fig. 6).

To visually delineate the variations in metabolite levels and their interrelationships under cold stress, we performed heatmap visualization of the concentration of key compounds, including JA, total flavonoids, lignin, anthocyanin, TOC, and metabolites from the artemisinin biosynthetic pathway in A. annua (Fig. 4A). We observed that the contents of lignin, total flavonoids, anthocyanins, artemisinin acid, and TOC in A. annua leaves exhibited similar trends under cold treatment, with significant increases at CD2 and CD7. Notably, the levels of JA and artemisinin followed a comparable pattern, showing a marked rise at CH6.

Fig. 4
figure 4

Heatmap visualization and PLS-DA analysis of secondary metabolites. (A) The heatmap of secondary metabolite concentrations in A. annua leaves under cold stress. (B) Constructed PLS-DA loadings and scores plot from metabolite concentrations in samples. (C) Permutation test with 200 permutations of PLS-DA model.

Through the construction of the partial least squares-discriminant analysis (PLS-DA) model27, we found that control and cold treatment samples were clustered according to their metabolite types and quantities (Supplementary file 2). As shown in Fig. 4B, the four different groups were completely discriminated, indicating that cold treatment markedly affects artemisinin and phenylpropanoid biosynthesis pathways. The model parameters were R2X = 0.956, R2Y = 0.962, and Q2 = 0.948 (Fig. 4B). The results indicated that the model was stable and had a strong predictive ability. And the validity of the PLS-DA model was accessed with a permutation test, with a result of R2 = 0.079 and Q2 = − 0.549 (Fig. 4C). These results indicated that there was no over fitting phenomenon in the model. In our model, metabolites from these two biosynthetic pathways were crucial for discriminating between cold-treated groups of different durations. As shown in Fig. 4B, lignin, anthocyanins, and flavonoids were more closely associated with CD7 samples, while artemisinin and JA were more closely related to CH6 samples. Notably, dihydroartemisinin and artemisinic acid also clustered with CD7 samples. Extended cold stress may suppress the pathways from artemisinic acid to artemisinin B and from DHAA to artemisinin, leading to increased synthesis of lignin and flavonoids in A. annua.

Cluster analysis of DEGs related to artemisinin and JA biosynthesis under cold stress

Transcription factors (TFs) have been reported to regulate not only the plant’s cold stress response but also the biosynthesis of secondary metabolites15,28. Based on the PlantTFDB database, a total of 763 transcription factor-encoding genes exhibited differential expression (|Log2 Fold change| >1; FDR < 0.05) under cold stress, predominantly including genes that encode transcription factors (TFs) from the WRKY, AP2/ERF, MYB, NAC, bHLH, and bZIP families (Supplementary file 3). To further characterize the TFs related to artemisinin biosynthesis under cold stress, all differentially expressed TFs, JA, artemisinin and phenylpropanoid biosynthesis pathways genes were subjected to the trend clustering analysis based on Mfuzz package using R 4.3.2 software and divided into 10 clusters (Supplementary Fig. 7; Supplementary file 4).

In our study, artemisinin and JA biosynthesis genes were primarily enriched in cluster 2,4,5,7 and 9 (Supplementary Fig. 7). Strikingly, ADS, DRB2 and one AACT from artemisinin biosynthesis-related pathway, and AOS and AOC from JA biosynthesis-related pathway were also divided into cluster 2 (Supplementary Fig. 7; Supplementary file 4). Then, we conducted a screen to identify transcription factors that JA affects the biosynthesis of artemisinin, including ORCA29,30,31, ERF132, WRK133, HY534 and MYC210. However, we did not detect any genes encoding the ORCA transcription factors. Interestingly, a subset of AP2-domain transcription factors was also found in cluster 2 (Supplementary Fig. 7; Supplementary file 4), leading us to hypothesize that these factors may serve as candidate regulators analogous to ORCA (Fig. 5A). In cluster 2 (Supplementary Fig. 7; Supplementary file 4), we identified a gene encoding ERF1, along with two genes encoding bZIP transcription factors HY5 and four genes encoding WRKY1 (Fig. 5B and C). Moreover, we also identified a gene encoding for MYC2 under cold stress in cluster 4 (Fig. 5C; Supplementary Fig. 7).

Fig. 5
figure 5

Transcription factors related to A. annua synthesis. (A) Box plots illustrating the expression of candidate ORCA transcription factor genes. (B) Box plots representing the expression profiles of genes encoding ERF1 and WRKY1 transcription factors. (C) Box plots representing the expression profiles of genes encoding HY5, MYC2 and NAC1 transcription factors. (D) Box plots illustrating the expression profiles of genes encoding MYB24, MYB8 and MYB14 transcription factors.

Parallel transcriptional regulation of artemisinin and phenylpropanoid biosynthesis pathways

Lignin, a complex phenolic polymer, serves a crucial role in maintaining the structural integrity of plant cells under low-temperature stress. Our study observed a significant increase in lignin content at CD2 and CD7 (Fig. 4A), which aligns with the concurrent elevation in flavonoid levels, suggesting a coordinated upregulation of these secondary metabolites in response to cold stress. In the trend clustering analysis, 4CL and C4H from the upstream of phenylpropane biosynthesis pathway, HCT, CCoAOMT, CCR, HCT and CAD in the lignin synthesis pathway, and CHS, CHI, F3H, CYP75B1, DFR, LDOX and UGTs in the flavonoid synthesis pathway were mainly enriched in cluster 1,3,4,6,9 and 10 (Supplementary Fig. 7; Supplementary file 4). In cluster 3 and 10 (Supplementary Fig. 7; Supplementary file 4), we observed a significant upregulation of genes encoding for NAC1 and MYB24, respectively (Fig. 5C, D). In contrast, the expression of a gene encoding the MYB8 transcription factor was found to be downregulated under cold stress in cluster 4 (Fig. 5D; Supplementary file 4).

Flavonoids and terpenoids, constituting the two most abundant classes of specialized plant metabolites, originate from distinct biosynthetic routes. However, research conducted at the molecular level, encompassing the roles of transcription factors (TFs), structural genes, and the interplay of biochemical compounds, has indicated potential cross-talk and interconnections between these seemingly divergent pathways. In our study, we also found that some genes encoding PAL, CCoAOMT, C3H, and UGTs are enriched in cluster 2 (Supplementary Fig. 7; Supplementary file 4). Then, TFs related to regulate artemisinin and phenylpropanoid biosynthesis pathways were screened, including HY535,36 and MYB1437. In cluster 2 and cluster 10 (Supplementary Fig. 7; Supplementary file 4), in addition to HY5, we discovered significant upregulation of two genes encoding MYB14 specifically at CH6 and CD2 (Fig. 5D). Beyond transcription factors, evidence suggests that genes encoding enzymes may also serve as crucial connectors between these two pathways. In our study, differentially expressed genes encoding PAL, F3H, FLS, ADS, and DBR2 exhibited a convergent expression pattern under cold stress, specifically enriched in cluster 2, 3, 6 and 10 (Supplementary Fig. 7; Supplementary file 4).

Identification of RNAseq data using quantitative polymerase chain reaction

To verify the accuracy of the RNAseq data, DEGs involved in such pathways as artemisinin, phenylpropanoid biosynthesis and TFS were selected and examined by quantitative polymerase chain reaction (qPCR). The up- or downregulation of each DEG in qPCR was consistent with the results of the RNAseq analysis (Supplementary Fig. 8).