Integrated analysis of lncRNAs, mRNAs, and TFs to identify network modules underlying diterpenoid biosynthesis in Salvia miltiorrhiza

Plant materials and MeJA treatment

S. miltiorrhiza seedlings were cultured in the greenhouse under 23 °C–−25 °C. The plants were sprayed with 200 μM MeJA solution as mentioned in a previous report (Luo et al., 2014). After being treated with MeJA solution for 6, 12, 24, and 48 h, the treated and the non-treated (0 h) roots were collected and rinsed with water. These roots were dried gently and quickly with absorbent paper. The cleaned roots were immediately frozen in liquid nitrogen and stored at −80 °C until RNA isolation. Three biological replicates were carried out for each experiment.

S. miltiorrhiza genomic and transcriptomic data

The S. miltiorrhiza genome data were downloaded from NCBI (National Center for Biotechnology Information) Sequence Read Archive database (SRP051524, SRP051564, SRP028388) (Xu et al., 2016) and NCBI (National Center for Biotechnology Information) BioProject: PRJNA682867 (Ma et al., 2021). The transcriptomic data were obtained from our previous study (Zhou et al., 2021) with an accession number assigned to PRJNA712174. The transcriptomic data were gathered from four varieties of S. miltiorrhiza root tissues during the tanshinone accumulation stage and included Fragments Per Kilobase Million (FPKM) expression value and annotation information.

Pipeline for lncRNA identification

The following four filters were used to shortlist the bona fide lncRNAs from the transcriptomic data: (1) transcripts with annotation information in one of the Nr (Non-RedundantProtein Sequence), Swiss-Prot, and COG/KOG databases were removed; (2) transcripts shorter than 200 nt with an open reading frame (ORF) longer than 100 aa were discarded, found and extracted ORFs on getorf (http://emboss.bioinformatics.nl/cgi-bin/emboss/getorf), selection of ORF length <100 aa by R (version 4.2.1) script; (3) transcripts were searched against the Pfam database (Punta et al., 2012) (http://pfam.xfam.org/) by HMMER to remove transcripts possibly containing known protein domains; and (4) the protein-coding potential of each transcript was calculated using PLEK (Li, Zhang & Zhou, 2014) and Coding Potential Calculator 2 (CPC2, http://cpc2.cbi.pku.edu.cn) (Kang et al., 2017), and these with PLEK and CPC2 scores >0 were discarded. Through the above process, identified lncRNAs were obtained. FPKM value less than 0.05 was used as the standard for low expression levels (Li et al., 2016). The transcripts that remained were regarded as candidate lncRNAs.

Characterization and conservation analysis of S. miltiorrhiza lncRNAs

To gain more understanding of these lncRNAs in S. miltiorrhiza, we compared several different features of the lncRNAs and mRNAs: GC content and transcript length. GC content and transcript length were determined by R (version 4.2.1), and statistical analysis was carried out by Excel.

We aligned the lncRNA sequences identified here with BLAST+ (Camacho et al., 2009) (blast−2.11.0+) against the genome sequences of the Lamiaceae family: Salvia splendens, Salvia hispanica, Mentha longifolia, Scutellaria baicalensis, Pogostemon cablin, and Sesamum indicum, of which Salvia splendens and Salvia hispanica both belong to the Salvia genus; Solanaceae family: Nicotiana tabacum, Brassicaceae family: Brassica napus and Arabidopsis thaliana; and Selaginellaceae family: Selaginella moellendorffii. A cutoff E-value <1e−10 was used. The genomes were downloaded from the NCBI databases: GCF_004379255.1 (SspV2), GCF_023119035.1 (UniMelb_Shisp_WGS_1.0), GCA_001642375.2 (Mlong_CMEN585_v), GCA_005771605.1 (ASM577160v1), GCA_023678885.1 (GZUCM_PCab_1.0), GCF_000512975.1 (S_indicum_v1.0), GCF_000715135.1 (Ntab-TN90), GCF_020379485.1 (Da-Ae), GCF_000001735.4 (TAIR10.1), and GCF_000143415.4 (https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000143415.4/). The lncRNAs that had coverage of >20% of matched regions were defined as conserved lncRNAs.

Precursors of miRNA and miRNA target prediction in S. miltiorrhiza lncRNAs

The candidate lncRNAs may act as the precursors of miRNAs, the S. miltiorrhiza miRNAs in miRBase (Kozomara, Birgaoanu & Griffiths-Jones, 2019) (Release 22.1, http://www.mirbase.org/) and PmiREN2.0 (Guo et al., 2022) (https://pmiren.com/) were aligned to the sequences of the candidate lncRNAs. The secondary structure of lncRNAs was predicted by RNAfold (Gruber et al., 2008) (http://rna.tbi.univie.ac.at/). LncRNAs with classical stem-loop hairpins were regarded as the putative precursors of miRNA (Zhou et al., 2017). Given that miRNA targets and lncRNAs have highly similar miRNA-binding sites, miRNA can be sequestered by lncRNA (Paschoal et al., 2017). Three kinds of prediction software were used to determine the miRNAs targeted to candidate lncRNAs. The first was TAPIR (Bonnet et al., 2010) (http://bioinformatics.psb.ugent.be/webtools/tapir); this server offers the possibility to search for plant miRNA targets using a fast and precise algorithm (score ≤ 4, free energy ratio ≥ 0.7). The second was psRobot (Wu et al., 2012b) (Version 1.2) with default parameters, which is a widely used online miRNA target prediction tool. The third was psRNATarget (2017 release) with default settings, and psRNATarget was developed to identify plant sRNA targets by (i) analyzing complementary matching between the sRNA sequence and target mRNA sequence using a predefined scoring schema and (ii) by evaluating target site accessibility (Dai, Zhuang & Zhao, 2018), targets with an E value less than 5.0 were retained as potential miRNA targets.

Subcellular localization of lncRNAs

The subcellular localization of S. miltiorrhiza lncRNAs was predicted by lncLocator (Cao et al., 2018). LncLocator is an ensemble classifier-based predictor, which adopts both k-mer features and high-level abstraction features generated by unsupervised deep models and constructs four classifiers by feeding these two types of features to support vector machine and random forest, respectively. The current lncLocator can predict five subcellular localizations of lncRNAs, namely, cytoplasm, nucleus, cytosol, ribosome, and exosome.

Predicted lncRNAs related to the diterpenoid biosynthesis in S. miltiorrhiza

Diterpenoid biosynthetic genes and the genes encoding TFs involved in diterpenoid biosynthesis were screened from the transcriptomic annotation data of S. miltiorrhiza. Our research focused on the downstream of the diterpenoid biosynthetic pathway (Fig. S1).

The functions of lncRNA can be performed on diterpenoid biosynthetic genes/TFs in a cis manner (Kopp & Mendell, 2018), and the lncRNAs and their target genes are considered lncRNA-mRNA/TF pairs. Pearson correlation coefficient (PCC) of the expression of lncRNA and its target gene pair was calculated using the psych (Revelle, 2022) package in R, requiring co-expressed gene pairs with —PCC— ≥ 0.4 and p ≤ 0.05. The basis for predicting cis-regulation was related to the positional relationship of lncRNA genes and coding genes on the genome (Liu et al., 2019). It was determined to be cis-regulatory if the lncRNA gene was within 100 kilobases (kb) upstream or downstream regions of the target genes used (Huang et al., 2018; Wang et al., 2021). used. For genomic location analysis, the co-expressed gene pairs were aligned against two S. miltiorrhiza genomes (Xu et al., 2016; Ma et al., 2021) by BLAST+ (BLAST−2.11.0+) with a cutoff of E value <1e−5. Finally, the candidate lncRNA-mRNA/TF pairs related to the diterpenoid biosynthesis of S. miltiorrhiza were obtained.

Construction of the diterpenoid biosynthesis-related lncRNA-mRNA-TF networks and identification of hub genes

Cytoscape 3.2 (Shannon et al., 2003) was utilized to draw the putative lncRNA-mRNA-TF co-expression network with —PCC— ≥ 0.8 (p ≤ 0.05). Genes with a degree ≥ 10 were considered hub genes.

Candidate lncRNA-mRNA/TF pairs expression profile in MeJA-induced S. miltiorrhiza

MeJA is an effective elicitor for the production of diterpenoid tanshinones in S. miltiorrhiza (Luo et al., 2014). To further explore whether candidate lncRNA-mRNA/TF pairs have a response to the diterpenoid biosynthetic pathway, we performed the quantitative real-time polymerase chain reaction (qRT-PCR) to analyze expression patterns of the detected lncRNA-mRNA/TF pairs under MeJA induction in S. miltiorrhiza. Plant materials treated with MeJA dissolving media were used as a control (0 h), and three biological replications were carried out. The 2(-delta-delta CT) method was used for calculating the relative expression levels of genes. ANOVA was calculated using SPSS (Version 23.0, IBM, USA), and p < 0.05 and p < 0.01 were considered statistically significant. The hub genes were identified by using the cytoHubba (Chin et al., 2014) plug-in of Cytoscape 3.2. To increase the sensitivity and specificity, we proposed Maximal Clique Centrality (MCC) to discover featured nodes. Subsequently, we constructed a network module of these hub genes.

RNA extraction and qRT-PCR

Total RNA was extracted from plant tissues using the RNeasy plant kit (PH-01013-B, Foregene, Chengdu, China). RNA degradation and contamination were monitored using 1% RNase-free agarose gel electrophoresis, and the RNA purity was analyzed using a NanoPhotometer™ -N60 ultra-micro spectrophotometer. The reverse transcription reaction was used RT Easy™ II (With gDNase) kit (Version Number: 1.0-1904, Foregene, Chengdu, China) following the instruction manual. qRT-PCR was performed in a 20 μL reaction volume containing primers by Real-Time PCR EasyTM-SYBR Green I kit (Cat. No. QP-01011/01012/01013/01014, Foregene, Chengdu, China) following the instruction manual. qRT-PCR was carried out in triplicate reactions in the LineGene K Plus real-time PCR detection system (Bioer, Hangzhou, China). Primers were designed using the Primer-BLAST (Ye et al., 2012). Primer sequences were listed in Table S1. β-actin (Yang et al., 2010) gene was selected as a reference since it showed stable expression in the S. miltiorrhiza tissues analyzed compared to others. The reaction program was as follows: 3 min at 95 °C, 10 s at 95 °C, 30 s at 60 °C, and 40 cycles. The temperature was then gradually increased to produce melting curves for amplification specificity verification. The mean value of three replicates was normalized using β-actin as the endogenous control. The 2(-delta-delta CT) method was used for calculating the relative expression of genes.

Identification of candidate lncRNAs

On the basis of the strict lncRNA identification pipeline, a total of 30,347 transcripts with no annotation information were obtained. A total of 21,742 transcripts with length >200 nt and ORF <100 aa were obtained. Subsequently, a total of 21,468 identified lncRNAs were obtained from the intersection of PLEK and CPC2 software prediction results (Fig. 1, Tables S2, S3). Expression transcripts with FPKM ≤ 0.05 were filtered, and a total of 6,651 candidate lncRNAs were selected for further analysis (Table S4).

The 21,468 identified lncRNAs ranged from 201 nt to 4,340 nt in length, the majority of which (approximately 68.95%) were 200–400 nt. The mean length was 392 nt, which was lower than the values observed in S. miltiorrhiza mRNAs (mean length = 1,515 nt). The GC content of lncRNAs was mainly concentrated at 31.28%–37.28% (accounting for 46.32%), whereas mRNAs were mainly concentrated at 39.28%–41.28% (accounting for 36.00%). The mean GC content of lncRNAs was approximately 36.86%, which was slightly lower than that of mRNA sequences (approximately 43.42%) (Fig. 2, Table S5).

Characterization and conservation of lncRNAs in S. miltiorrhiza

Conservation analysis showed that 21.02% and 15.40% of S. miltiorrhiza lncRNAs (6,651 candidate lncRNAs) were conserved compared with Salvia splendens and Salvia hispanica of the same genus, respectively. However, in the same family of the different genera of Mentha longifolia, Scutellaria baicalensis, Pogostemon cablin, and Sesamum indicum, 6.54%, 0.96%, 0.83%, and 0.51% of S. miltiorrhiza lncRNAs were conserved, respectively. Among plants of different families, only 0.05%, 0.02%, and 0.02% of S. miltiorrhiza lncRNAs were conserved in Nicotiana tabacum, Arabidopsis thaliana, and Brassica napus, respectively. There was no match in Selaginella moellendorffii. The results of the conservation analysis are presented in Table S6.

LncRNAs might be used as precursors or target mimics of S. miltiorrhiza miRNAs

In this study, a total of 17 lncRNAs were identified as potential precursors of 41 S. miltiorrhiza miRNAs (Table S7). To improve the prediction accuracy, the intersection of three kinds of software prediction results was selected. A total of 14 lncRNAs were identified as potential targets of 66 S. miltiorrhiza miRNAs from 18 families according to PmiREN2.0 (Table S7).

Subcellular localization predictions indicated that most S. miltiorrhiza lncRNAs were found in the cytoplasm and nucleus

Subcellular localization of lncRNA is closely related to its function. At present, many studies have indicated that lncRNA contains specific RNA motifs with nuclear localization (Zhang et al., 2014). Meanwhile, in the subcellular localization of lncLocator (Cao et al., 2018), a total of 3,381 cytoplasms, 3,316 nucleus, 83 exosomes, 51 cytosols, and 20 ribosomes were obtained from 6,651 candidate lncRNAs. Among our candidate lncRNA-mRNA/TF pairs, 13 lncRNAs were found in the nucleus, 9 lncRNAs were found in the cytoplasm, and 1 lncRNA was found in the exosome (Table S8).

Genes and TFs involved in diterpenoid biosynthesis from transcriptomic data

To investigate the potential lncRNAs involved in diterpenoid biosynthesis, we predicted the potential targets of lncRNAs in cis-regulatory relationships. In transcriptomic annotation, we obtained 46 diterpenoid biosynthetic genes: GPPS (Van (Schie et al., 2007), FPPS, GGPPS, CPS, ent-CPS, KS, ent-KS, CYP76AH1 (Guo et al., 2013; Ma et al., 2016), CYP76AK2, CYP76AK3 (Li et al., 2021), CYP76AK5 (Xiangdong & Lizhi, 2017), ent-kaurene oxidase (KO) (Hedden & Thomas, 2012), ent-kaurenoic acid oxidase (KAO) (Helliwell et al., 2001), GA 2-oxidase (GA2ox), GA 3-oxidase (GA3ox), GA 20-oxidase (GA20ox) (Hedden & Thomas, 2012; Du et al., 2015), and 11 TFs: bHLH148, ERF6, GRAS1, MYB36 and WRKY2 (Li & Lu, 2014; Zhang et al., 2015; Li et al., 2015; Li et al., 2019; Ji et al., 2016), belonging to five TF families: bHLH, ERF, GRAS, MYB, and WRKY. The list of the above genes is shown in Table S9.

LncRNAs related to the diterpenoid biosynthetic pathway

We detected 6455 lncRNAs that were co-expressed with 46 genes and 11 TFs involved in diterpenoid biosynthesis. In the PCC matrix, we obtained 45,198 correlation gene pairs with —PCC—≥ 0.4 (p ≤ 0.05), of which 47,946 pairs were positively correlated and 5,775 pairs were negatively correlated. At the positional level, combining the results from the two S. miltiorrhiza genomes (Xu et al., 2016; Ma et al., 2021), a total of 180 pairs of lncRNA-mRNA/TF were located adjacent to each other within 100 kb (including ±3 kb), of which 23 pairs were co-located and co-expressed. A total of 23 candidate lncRNA-mRNA/TF pairs were obtained (Table S10).

LncRNA-mRNA-TF networks and hub genes

The gene co-expression network was constructed with —PCC— ≥ 0.8 (p ≤ 0.05) as the threshold, the visualization is shown in Fig. S2A. We obtained 1,402 gene pairs with —PCC—≥ 0.8. We obtained 24 mRNAs with a degree > 10, which were considered hub genes (Table S11). In the 23 gene pairs of lncRNA-mRNA/TF, 8 mRNAs were present in these hub genes. Only two gene pairs Smlnc0032870-Sm0012648 (GA2OX) and Smlnc0018769-Sm0037093 (KS2) existed in the network with high expression correlation of 0.85 and 0.80, respectively. The hub genes Sm0012648 and Sm0037093 in these two gene pairs were used to construct the subnetworks (Figs. S2B, S2C).

Tissue specificity of lncRNA expression

Using the qRT-PCR, expression patterns of lncRNAs in the candidate lncRNA-mRNA/TF pairs were analyzed in roots, stems, leaves, and flowers of 2-year-old field greenhouse-grown S. miltiorrhiza plants. Of them, 11 lncRNAs were detected in at least one tissue and showed tissue-specific expression. The other lncRNAs were undetected, suggesting that they could be not expressed or expressed at a low level in the tissues analyzed. Smlnc0012647, Smlnc0032870, Smlnc0042160, Smlnc0063419, and Smlnc0070114 exhibited the highest expression in root tissue. Smlnc0000154, Smlnc0008477, Smlnc0019429, and Smlnc0052170 were more stem-specific. Smlnc0018769 was expressed mainly in flowers and stems. Smlnc0008662 showed high expression in roots, flowers, and leaves and low expression in stems (Fig. 3). Thus, lncRNAs showed obvious differential expression in different tissues, which may be related to their regulatory function. These results suggest that the expression of lncRNAs may be limited to specific tissue types or regulated by development in S. miltiorrhiza.

Time-series expression pattern of candidate lncRNA-mRNA/TF pairs induced by MeJA and module detection in S. miltiorrhiza

To validate whether the candidate lncRNA-mRNA/TF pairs are related to diterpenoid biosynthesis, we analyzed the time-series expression patterns of candidate gene pairs in response to MeJA treatment. MeJA treatment significantly changed the expression of 19 genes in S. miltiorrhiza at least a time-point of MeJA treatment (Fig. S3). As shown in Fig. 4, the expression of most genes showed a downward trend at the 6 h time point, and only the Smlnc0018769-Sm0037093 pair and Smlnc0008662-Sm0063385 pair had an upward trend. The Smlnc0008477-Sm0056000 pair and Smlnc0012647-Smlnc0032870 pair had the same expression trend. The Smlnc0042160-Sm0028870 pair had the same expression trends at 6, 12, and 24 h. The Smlnc0018769-Sm0037093, Smlnc0008662-Sm0063385, and Smlnc0052170-Sm0067296 pairs had opposite expression trends at 48 h, whereas the Smlnc0019429-Sm0026208 pair had opposite trends at 24 and 48 h. Interestingly, the two gene pairs in the co-expression network, Smlnc0032870-Sm0012648 and Smlnc0018769-Sm0037093, which were positively correlated, showed complex expression correlation under the induction of MeJA, even showed opposite expression trends after 24 h. The hub genes: Sm0056000, Sm0037097, Sm0063385, and Sm0067296 in the co-expression network showed significantly differential expression in MeJA-induced S. miltiorrhiza.

The data in Fig. 5A showed that Smlnc0063419, Sm0026208, and Sm0067296 had similar expression patterns, which formed the lncRNA-mRNA-TF module. Sm0026208 and Sm0067296 were annotated with TF WRKY2 and mRNA GA3ox2, respectively. Smlnc0012647, Smlnc0032870, and Sm0009433 had similar expression patterns, as shown in Fig. 5B. Sm0009433 was annotated with TF MYB36. Smlnc0042160 and Sm0037093 had similar expression patterns, as shown in Fig. 5C. Sm0037093 was annotated with mRNA KS2. Through the plug-in cytoHubba, we calculated the top 10 hub genes (Table S12 and Fig. 5D). The modules of Smlnc0063419-Sm0026208-Sm0067296 and Smlnc0012647-Smlnc0032870-Sm0009433, in which all genes were present in the top 10 hub genes.