Identification of key genes involved in secondary metabolite biosynthesis in Digitalis purpurea

Abstract

The medicinal plant Digitalis purpurea produces cardiac glycosides that are useful in the pharmaceutical industry. These bioactive compounds are in high demand due to ethnobotany's application to therapeutic procedures. Recent studies have investigated the role of integrative analysis of multi-omics data in understanding cellular metabolic status through systems metabolic engineering approach, as well as its application to genetically engineering metabolic pathways. In spite of numerous omics experiments, most molecular mechanisms involved in metabolic pathways biosynthesis in D. purpurea remain unclear. Using R Package Weighted Gene Co-expression Network Analysis, co-expression analysis was performed on the transcriptome and metabolome data. As a result of our study, we identified transcription factors, transcriptional regulators, protein kinases, transporters, non-coding RNAs, and hub genes that are involved in the production of secondary metabolites. Since jasmonates are involved in the biosynthesis of cardiac glycosides, the candidate genes for Scarecrow-Like Protein 14 (SCL14), Delta24-sterol reductase (DWF1), HYDRA1 (HYD1), and Jasmonate-ZIM domain3 (JAZ3) were validated under methyl jasmonate treatment (MeJA, 100 μM). Despite early induction of JAZ3, which affected downstream genes, it was dramatically suppressed after 48 hours. SCL14, which targets DWF1, and HYD1, which induces cholesterol and cardiac glycoside biosynthesis, were both promoted. The correlation between key genes and main metabolites and validation of expression patterns provide a unique insight into the biosynthesis mechanisms of cardiac glycosides in D. purpurea.

Fig 1. Flowchart of gene co-expression network analysis in D. purpurea.

Data collection and analysis towards downstream analyses are shown.

Fig 2. Selection of an appropriate soft threshold power of β.

Left; the scale-free fit index of power (β) was estimated to be 16 based on the threshold limit of 0.9. Right; the mean connectivity versus soft-thresholding power.

Fig 3. Clustering of 180 modules.

The horizontal red line shows the height cut of 0.25, which corresponds to the correlation of 0.75, to merging the modules.

Fig 4. Clustering of genes and modules.

The cluster dendrogram at the top of the plot shows co-expressed genes. The branches and color bands at the bottom of the plot represent the assigned module.

Fig 5. Correlation of modules and secondary metabolites.

Module Eigengenes (MEs) and secondary metabolites are respectively represented by each row and column. Each cell contains the corresponding correlation at the top and P-value at the bottom. The positive correlation of the module with secondary metabolites and the negative correlation are respectively shown in red and green. The white spectrum indicates the inexistence of modules and secondary metabolites correlations.

Fig 6. The scatterplot of Gene Significance (GS) for digitoxigenin bis-digitoxoside vs. Module Membership (MM) in chocolate3.

The high significant correlation between GS and MM for Digitoxigenin bis-digitoxoside in the chocolate3 module.

Fig 7. Eigengene dendrogram of the modules and digitoxigenin bis-digitoxoside.

A hierarchical clustering dendrogram of eigengenes. The dissimilarity of E_I and E_J is shown by 1- cor (E_I; E_J).

Fig 8. The heatmap plot of adjacencies in the eigengene network, including digitoxigenin bis-digitoxoside.

Each row and column in the heatmap corresponds to one Module Eigengene (labeled by color) and or Digitoxigenin bis-digitoxoside. Low and high adjacencies (negative and positive correlations) are shown in blue and red. The connection strength (adjacency) between eigengenes I and J are defined as A _{I J} = (1 + cor (E_I; E_J))/2.

Fig 9. The gene ontology analysis of candidate modules.

Gene ontology analysis of candidate modules classified into three functional categories including Biological Process, Molecular Function, and, Cellular Component.

Fig 10. The KEGG analysis of candidate modules.

KEGG analysis of candidate modules is performed by the DAVID database.

Fig 11. A total number of 140 hub genes in seven modules.

The top 20 hub genes in each selected module, including blue2 (A), chocolate3 (B), coral3 (C), coral4 (D), darkorange2 (E), lightpink4 (F), and lightsteelblue (G) modules. Nodes represent genes in the network and red nodes indicate the hub genes. The gray line connecting two nodes indicates their connection.

Fig 12. The protein-protein interactions of hub proteins.

Response to the stimulus (GO:0050896), cellular macromolecule metabolic process (GO:0044260), metabolic process (GO:0008152), and cellular process (GO:0009987) are respectively shown in red, blue, green, and yellow.

Fig 13. The relative expression of candidate genes.

The expression pattern of JAZ3, SCL14, DWF1, and HYD1 was evaluated under 100 μM MeJA treatment in Digitalis purpurea L. X axis represents the fold change of the expression value. Y axis represents time points. Vertical bars indicate ± SE of the mean (n = 3). The different letters on columns represent the significant difference given by Duncan’s multiple range tests (P-value ≤ 0.05). There were no significant differences between equal letters (P-value ≤ 0.05). Although the expression of JAZ3 is early induced, a significant suppression is observed after 48 h. Other key genes particularly SCL14, targeting DWF1, and HYD1, inducing cholesterol biosynthesis and cardiac glycoside content, showed a significant increase in the expression.

Fig 14. The comparative relative expression of the candidate genes under 100 μM MeJA treatment in Digitalis purpurea L.

Here, the expression pattern of JAZ3, SCL14, DWF1, and HYD1 is compared. X axis represents the fold change of the expression of candidate genes. Y axis represents time points. Vertical bars indicate ± SE of the mean (n = 3). The different letters on columns represent the significant difference given by Duncan’s multiple range tests (P-value ≤ 0.05). There were no significant differences between equal letters (P-value ≤ 0.05). The interesting point is that the expression of HYD1, involved in the biosynthesis of cholesterol and subsequently the cardiac glycosides, shows a stable increase after 48 h. In contrast, other genes are suppressing and JAZ3, inducing the downstream genes, shows a significant suppression.