Differentially expressed genes(DEGs)analysis
At each sampling time point, genes with the normalized expression values (TPM) ≥ 2-fold change and an adjusted P-value < 0.5, compared with the control group, were considered as differentially expressed genes. Based on these criteria, a total of 748, 228, 239, 434, 887, and 1661 DEGs were identified after LOV treatment for 2, 4, 6, 8, 10, and 12 h, respectively. The sample treated with LOV for 12 hours showed the highest number of up-regulated genes (Fig. 1a). No genes have yet been found to be significantly differentially expressed in all six sampling time points. It has also been observed that when S.miltiorrhiza was first exposed to LOV, there was a period of rapid and significant changes in gene expression within a short timespan of 2 hours. Subsequently, the number of DEGs exhibited a tendency to gradually increase over time.
GO enrichment analysis of DEGs
The DEGs were analyzed by GO enrichment to characterize their functions at each sampling time points. The significantly enriched GO terms were identified based on the criterion of a P-value < 0.05 (Fig. 2a-f). The top 10 significantly enriched GO terms in the biological process (BP) category at each time point are presented in Supplementary Table S1. The results indicated that the BPs varied significantly at different time points. Among these BPs, the terms with the most DEGs at each time point were, in the following order: “carbohydrate metabolic process” (GO:0005975, 2 h, 46 DEGs), “terpene biosynthetic process” (GO:0016114, 4 h, 8 DEGs), “lipid X metabolic process” (GO:2001289, 6 h, 2 DEGs), “RNA processing” (GO:0006396, 8 h, 25 DEGs), “isoprenoid biosynthetic process” (GO:0008299, 10 h, 18 DEGs), “metabolic process” (GO:0008152, 12 h, 504 DEGs).
KEGG enrichment analysis of DEGs
KEGG is a knowledge base for the systematic analysis of gene functions, linking genomic information with higher-order functional information [41]. In this study, the DEGs were further analyzed using KEGG enrichment. The results showed that DEGs in all six sampling time points could be assigned to 26 pathways. Among these pathways, 24 pathways were related to metabolism. It indicated that the inhibitor primarily affected plant metabolism. After the 2 h, 10 h and 12 h of LOV treatment, “photosynthesis (map00195)” was significantly enriched. “Inositol phosphate metabolism (map00562)” was significantly enriched at 2 h and 4 h. “Biosynthesis of various plant secondary metabolites (map00999)” and “tyrosine metabolism (map00350)” were significantly enriched at 2 h and 12 h (Fig. 3).
Notably, DEGs enriched in various terpenoid biosynthesis-related processes at different time points. For instance, “monoterpenoid biosynthesis (map00902)” was significantly enriched in the samples treated with LOV for 2 h; “diterpenoid biosynthesis (map00904)” and “zeatin biosynthesis (map00908)” were significantly enriched at 4 h; “carotenoid biosynthesis (map00906)” was significantly enriched at 6 h; “sesquiterpenoid and triterpenoid biosynthesis (map00909)” and “terpenoid backbone biosynthesis (map00900)” were significantly enriched at 10 h (Fig. 3). This suggests that when plants are disturbed by LOV, the order of influence on terpenoid biosynthesis-related pathways is related to the topological relationship among these pathways. Interestingly, the effect on “terpenoid backbone biosynthesis (map00900)” occurs at the late stage.
Analysis of differentially expressed TFs
TFs play crucial roles in regulating plant growth and development, as well as in responding to environmental stresses. In this experiment, a total of 79 differentially expressed TFs were identified in the samples treated with LOV, involving 18 of the 56 transcription factor families in S.miltiorrhiza, namely bHLH, bZIP, C2H2, HD-ZIP, HSF, MADS, MYB, NF-YA, AT-hook, GRAS, AP2/ERF, NZZ/ SPL, Nin-like, SRS, TCP, WRKY, ZF-HD, and BES1. The one with the highest number of differentially expressed TFs was the AP2/ERF family, comprising a total of 16 differentially expressed TFs (Fig. 4a).
The results also showed that the genes from 18 transcription factor families were expressed differentially at various time points. The sample with the most transcription factor families in differentially expressed genes was taken from the 12th hour (Fig. 4b). AT-hook and Nin-like families were predominantly differentially expressed in the early phase, while bHLH, HSF, AP2/ERF, and WRKY families were predominantly differentially expressed in the late phase. Some transcription factor families are expressed differentially at multiple time points. For example, the WOX family was differentially expressed at 2 h, 4 h, 10 h, and 12 h; the HSF family was differentially expressed at 2 h, 6 h, 8 h, and 12 h; the bZIP family was differentially expressed at 4 h, 8 h, 10 h, and 12 h. SmWOX11-like, a member of the WOX family, was significantly differentially expressed in the samples treated with LOV for 4 h, 10 h, and 12 h. Details of the differentially expressed TFs at different time points are shown in Supplementary Table S2.
Analysis of GCNs
To investigate the co-expression dynamics of genes, a weighted gene co-expression network analysis was conducted. Total 116 co-expressed modules were identified. Among them, the turquoise module was the largest, containing 8,679 genes, while the steel blue and salmon modules showed a strong correlation with HMGR concentration (P-value < 0.05, Fig. 5a-c). The top 10 hub genes in the three modules were identified by analyzing nodes with centrality in the largest clusters (Table S3). Functional enrichment analysis of the three gene modules based on GO and KEGG was conducted to elucidate the underlying biological functions of these genes within their networks.
The genes of these three modules were enriched in various GO terms and metabolic pathways. The turquoise module was found to be significantly enriched in 944 GO terms (P-value < 0.05) (Fig. 6a). The KEGG enrichment analysis revealed that these genes can be categorized into 17 pathways (P-value < 0.05) (Fig. 7a), with 9 of these pathways belonging to the genetic information processing category. Analysis of the expression profiles of the top 10 hub genes in this module revealed that the peak expression of these genes occurred at hour 8 (Fig. S2a). It was indicated that the genes involved in genetic information processing biological processes were actively expressed after inhibitor treatment for 8 hours. A total of 623 and 331 GO terms were significantly enriched in the steel blue and salmon modules (P-value < 0.05) (Fig. 6b-c). The enriched terms indicated that the steel blue and salmon modules were primarily associated with macromolecular biosynthesis and sugar metabolism. The KEGG enrichment analysis revealed that these genes were assigned to 27 and 24 pathways (P-value < 0.05). Among them, 14 and 17 pathways are associated with energy metabolism (Fig. 7b-c). This indicates that many pathways related to energy metabolism responded positively to changes in HMGR activity. Analysis of the expression profiles of the top 10 hub genes in the steel blue and salmon modules revealed that the peak expression of these genes occurred at hour 10 (Fig. S2b-c). The results indicated that the genes in the energy metabolism pathway were actively expressed after 10 hours of inhibitor treatment.
Analysis of GRNs
Unlike the GCN, the GRN can be used to study the regulatory relationships between genes. It is known that gene regulation occurs in the nucleus, and the proteins encoded by these genes play roles in various subcellular locations. GO organizes biological knowledge about cell components hierarchically. This makes GRNs analysis not only interpretable mechanistically but also provides new mechanistic insights [42]. Therefore, in constructing GRNs, DEGs were first classified into five gene sets based on GO subcellular localization annotations. The GO terms of these five gene sets are related to “chloroplast”, “mitochondria”, “endoplasmic reticulum”, “Golgi apparatus”, and “cytoplasmic lysosomal”. Five GRNs were constructed based on these five gene sets and differentially expressed TFs (Fig. 8, Fig. S4), and the top 5 hub genes in each network were identified (Table S5).
The target genes of these specific regulators are involved in various biological pathways (Table S6). For example, in the chloroplast gene regulatory module, SmWRKY48 targeted genes involved in fatty acid biosynthesis, pyruvate metabolism, photosynthesis, arginine biosynthesis, and other pathways. In the mitochondrial gene regulatory module, the target genes of SmAGL16 were involved in various pathways, such as nicotinate and nicotinamide metabolism, pantothenate and CoA biosynthesis, amino acid metabolism and degradation, aminoacyl-tRNA biosynthesis, folate biosynthesis, and other biological processes. In the cytoplasmic lysate gene regulatory module, SmWRKY28 targeted genes involved in glutathione metabolism, pyruvate metabolism, cysteine and methionine metabolism, and glycolysis/gluconeogenesis pathways. Similarly, in A.thaliana, the homologous genes of these TFs also exhibit the behavior of targeting genes involved in various biological processes. For example, AtWRKY48, a homologous gene of SmWRKY48, targeted AtASP, AtIAA2, and AtDLDG1 in A.thaliana acclimation to a combination of high light and heat stress; AtAGL16, a homologous gene of SmAGL16, regulated the expression of AtHKT1;1, AtHsfA6a, and AtMYB102 in A.thaliana salt stress response; AtWRKY28, a homologous gene of SmWRKY28, targeted genes in the salicylic acid and jasmonic acid/ethylene-dependent defense pathways in A.thaliana streptococcus bacteriophage and oxalic acid stress response [43–45].
TFs, as a regulatory center, regulate S.miltiorrhiza responses to MVA pathway interference by integrating various biological pathways. In the cytoplasmic lysate gene regulatory network, SmWRKY48-SmTCP4-SmWRKY28 worked as a regulation hub to integrate various biological processes, including pyruvate metabolism, glycolysis/gluconeogenesis, amino acid metabolism, and ubiquinone and other terpenoid-quinone biosynthesis (Fig. 8). In order to further verify whether the identified TFs have a specific response to LOV, the expression profiles of TFs were analyzed after treatment with different concentrations of LOV (Fig. 9). It was found that the expression of SmWRKY28 and SmWRKY48 fluctuated greatly, while there was no significant change in SmTCP4 expression.
In this study, it was relatively common that a target gene was regulated by several transcription factors. For example, in the chloroplast gene regulatory module, the gene SmPSBR, which is involved in the photosynthesis pathway, was simultaneously targeted by the transcription factors SmCPRF1, SmABR1, and SmWRKY48. In the mitochondrial gene regulatory module, the gene SmPDL2, which is involved in the folate biosynthesis pathway, was targeted by the transcription factors SmWRKY48, SmAGL16, and SmERF010.
Validation of RNA-Seq sequencing data by qPCR analysis
To validate the reliability of transcriptome sequencing data, 10 genes were selected for expression correlation analysis: DXS, HMGS, HMGR, IDI, KSL, GA3, CPRF1, TCP4, WRKY40, and WRKY48. The results indicate that the majority of the selected genes displayed expression patterns consistent with the results of transcriptome sequencing, affirming the credibility of our transcriptome data (Fig. S3).