Transcriptome analysis reveals the GRAS family in Salvia miltiorrhiza hairy roots in response to gibberellin

Background: Salvia miltiorrhiza is one of the most important traditional Chinese medicinal plants with high medicinal value. Gibberellins are growth-promoting phytohormones that regulate numerous growth and developmental processes. However, their role in regulating secondary metabolism has not been investigated. Results: In this study, we found that GA can promote hairy roots growth and increase the content of tanshinones and phenolic acids. Transcriptome sequencing revealed that secondary metabolism pathway genes were enriched in the GA-responded. After further analysis of the changes of GA signaling pathway genes, it was found that the GRAS transcription factors family had a significant response to GA. We identified 35 SmGRAS genes in S. miltiorrhiza . SmGRAS genes can be divided into ten subfamilies in which members of the same subfamily showed similar conserved motifs and gene structures, suggesting the possible conserved functions. Conclusions: Most SmGRAS genes are significantly responsive to GA, indicating that they may play an important role in the GA signaling pathway and participate in GA regulation of root growth and secondary metabolism in S. miltiorrhiza .

do not meet the ever-increasing market demand. Methods of improving the secondary metabolites have been tried, such as adversity stress, addition elicitor, overexpressing or suppressing key enzyme genes or transcription factors involved in the biosynthetic pathways of these secondary metabolites.
However, the regulation of gibberellin on secondary metabolites biosynthesis remains unknown.
Gibberellins (GAs) are growth-promoting phytohormones that regulate numerous growth and developmental processes throughout the whole life cycle of the plant, including seed germination, root and stem elongation, and flower development [10]. Since the 1950s, more than 130 GAs have been identified in various plants (http://www.plant-hormones.info/gibberellin_nomenclature.htm) [11].
Only a few of them, such as G1, G3, G4 and G7, are bioactive [10]. GAs biosynthesis and catabolism pathways in plants have been well characterized. GAs are biosynthesized from the common precursor trans-geranylgeranyl diphosphate (GGPP), formed via the methylerythritol phosphate pathway, through the sequential action of two terpene cyclases CPS and KS, followed by oxidation by cytochrome P450 monooxygenases and then by 2-oxoglutarate-dependent dioxygenases [12] Subsequently, GA functions in plants through its signaling pathways [13]. Binding of GA to the nuclear receptor, GID1, causes a conformational change in the N-Ex of protein that promotes its association with GRAS domain of DELLA protein. This stable complex enables efficient SCF SLY1 recognition and subsequent degradation of DELLA by the proteasome [14].
The plant-specific GAI-RGA-and-SCR (GRAS) family proteins function as transcriptional regulators and play critical roles in GA signaling [15]. Most GRAS proteins comprise an N-terminal less-conserved variable region and a C-terminal conserved GRAS domain. Typical GRAS domains comprise 5 conserved sequence motifs: leucine heptad repeats I (LHRI), VHIID, leucine heptad repeats II (LHRII), PFYRE and SAW. Flanked by the two leucine-rich regions, the VHIID motif is present in all GRAS family members [16]. Based on amino acid sequences, the GRAS family is divided into 10 distinct subfamilies: DELLA, SCL3, LAS, SCL4/7, SCR, SHR, SCL9 (LISCL), HAM, PAT1, and DLT [16]. Different subfamily protein sequences have different characteristics and perform different functions. For example, DELLA proteins possess a conserved DELLA sequence motif in the N-terminal region. They function as repressors of GA and act as key regulatory targets in the GA signaling pathway in regulating plant growth [17,18]. SCL3 functions as a repressor of DELLA, which can positively regulate the GA signaling pathway and control GA homeostasis in Arabidopsis thaliana root development [19,20]. The SCL subfamily participates in root cell elongation, GA/DELLA signaling and stress response [15]. The VHIID and PFYRE motifs in the GRAS domain of SHR are essential for the interaction between SCR and SHR [16]. They could form a complex in order to participate in regulating root-related developmental processes in Arabidopsis [21,22]. The PAT1 subfamily has been shown to mediate phytochrome and defense signaling pathways [23]. And LISCL has two conserved subfamily-restricted acidic motifs in the N-domain and has been reported to be involved in stress response and adventitious root formation in response to auxin [16].
Although GA plays an important role in many aspects of plant growth and development, little is known about its role in regulating secondary metabolism. As diterpenoids, GA and tanshinones have common precursors GGPP [9]. There might be some correlation between GA and the biosynthesis of tanshinones. In addition, GRAS has crucial roles in the GA signaling pathway. Therefore, we speculated that SmGRASs might be involved in the GA regulating secondary metabolism in S. miltiorrhiza. To address these questions, we treated the hairy roots of wild-type S. miltiorrhiza with GA and determined the changes in biomass, tanshinones, phenolic acids, total phenolics, and total flavonoids. At the same time, we also measured the transcriptome changes and specifically analyzed the transcriptional level changes of the secondary metabolic pathway and GA signaling pathway genes. Finally, the bioinformatics of all SmGRAS family genes in S. miltiorrhiza and their responses to GA were analyzed. Our results revealed the possible pathways of GA regulating secondary metabolism and the response of SmGRASs to GA in regulating the secondary metabolism in S. miltiorrhiza, which provided a reference for the GA signaling pathway to regulate secondary metabolism.

Results
GA treatment effects root growth and secondary metabolism After 6-day cultivation, the GA-treated S. miltiorrhiza hairy roots showed bigger and redder than the control (Fig. 1A). The fresh and dry weights of the GA-treated hairy roots were all significantly higher than the controls (Fig. 1B, C). These results showed that the GA application promoted the growth of hairy roots. In order to investigate the changes of secondary metabolites affected by GA application, the contents of phenolic acids and tanshinones were determined after GA treatment. The two phenolic acids and four tanshinones contents in the hairy roots were all significantly increased after treated with GA (Fig. 1D, E). In addition, the contents of total phenolics and total flavonoids were also determined. However, the contents of total phenolics and total flavonoids were all decreased after GA treatment (Fig. 1F, G). Collectively, the data indicated that the GA treatment promoted the root growth, increased the accumulation of phenolic acids and tanshinones, while decreased the total phenolics and total flavonoids contents in the S. miltiorrhiza hairy roots.

Transcriptome-scale Analysis Of Ga-responsive Genes
To gain a comprehensive overview of the GA-responsive genes, we performed a transcriptome analysis of CK and GA-treated hairy roots. There were 10321 differentially expressed genes (DEGs) that were annotated in the volcano plot. The comparison of CK and GA-treated hairy roots revealed 4945 were GA-induced genes and 5376 were GA-repressed genes ( Fig. 2A). To verify the results from RNA-seq, 10 genes were selected for Quantitative Reverse-Transcription PCR (qRT-PCR) analysis. The results showed the expression patterns consistent with the RNA-seq, indicating that the RNA-seq data was reliable (Fig. S1). The global functional analysis of the DEGs was revealed that the "biological processes", "metabolic processes" and "cellular processes" were the top three categories in the most enriched gene ontology (GO) terms (Fig. 2B). And these DEGs identified were next assessed using the Kyoto Encyclopedia of Genes and Genomes (KEGG). The most significantly enriched term was biosynthesis of secondary metabolites, following by ribosome, plant-pathogen interaction and some primary and secondary metabolic pathways (Fig. 2C).
Secondary metabolism pathway genes in response to GA treatment In order to further explore the effect of GA on secondary metabolism, we used the Mapman program to analyze the transcriptome data (Fig. 3). The result showed that most DEGs of secondary metabolism were GA induced, especially the shikimate pathway, MVA pathway, simple phenols, betaines, wax, and anthocyanins. Diterpenoids such as GA and tanshinones are biosynthesized by MVA and MEP pathways. Most of the DEGs in the MVA pathway were GA induced at the transcriptional level, which might be related to the fact that GA induced tanshinones biosynthesis. Phenolic acids are mainly biosynthesized by shikimate and phenylpropane pathways. The results showed that most of the shikimate pathway DEGs and some of the phenylpropane pathway DEGs were induced by GA. In addition, GA also affected the biosynthesis pathways of flavonoids and alkaloids-like, wax and glucosinolates. In summary, GA regulated the transcription of many secondary metabolite synthesis pathway genes.

GA biosynthetic and signaling pathways genes in response to GA treatment
To investigate the effect of GA treatment on GA biosynthesis signaling pathway, we further analyzed and summarized the DEGs of GA biosynthesis and signaling pathway (Fig. 4). The results showed that the changes in transcriptional levels of GA biosynthesis pathway genes were diverse, the expressions of some DEGs were increased and others were decreased. In the downstream GA signaling pathway, the expression of GA receptor GID1 was increased. And the expression levels of most GRAS family genes, which were the key regulators of the GA signaling pathway, were increased. The expressions of F-box proteins SCF in the GA signaling pathway were different under GA treatment. In conclusion, GA could regulate the expressions of biosynthetic pathway genes and downstream signaling pathway genes, and further regulating many downstream physiological processes, such as cell growth, secondary metabolism, and plant resistance.

Identification and phylogenetic analysis of GRAS proteins in S. miltiorrhiza
In order to study the roles of SmGRAS in GA regulation of root growth and secondary metabolism of S. miltiorrhiza, we conducted a comprehensive analysis of the GRAS family genes in S. miltiorrhiza. We used HMMER to screen the protein sequences based on the HMM profiles from the S. miltiorrhiza genome database to identify putative GRAS proteins. We identified 35 SmGRAS proteins from S. miltiorrhiza, named GRAS1 ~ 35. The putative amino acid sequences of SmGRAS1 ~ 35 contained the conserved GRAS domain for the GRAS protein [15]. To study the evolutionary relationships of SmGRAS genes, we constructed a phylogenetic tree analysis with dicotyledons of Arabidopsis and monocotyledons of rice. We revealed that the SmGRAS proteins were divided into 10 subfamilies ( Fig. 5). Among these GRAS family proteins, there were 9 proteins from the PAT1 subfamily, 6 proteins from the LISCL subfamily, 5 proteins from the SHR subfamily, 4 proteins from the SCL and DELLA subfamilies, 2 proteins from the DLT and SCR subfamilies, and the other 3 subfamilies only have 1 protein. Therefore, we speculated that these SmGRAS genes might be similar to Arabidopsis GRAS genes of the same subfamily. They might be involved in the GA signaling pathway, root and stem cortex development, light morphogenesis and resistance to stress.

SmGRAS proteins sequence alignments and conserved motifs
In order to further confirm that these 35 genes belong to the GRAS family, we used DNAMAN and online MEME to perform multiple sequences alignment and conservative domain analysis on them.
The multiple sequences alignment result showed that the amino acid sequences of all these 35 proteins have high identities (Fig. 6A). Almost all SmGRAS proteins contained conserved GRAS domain consisting of LHRI, VHIID, LHRII, PFYRE and SAW. We also used online MEME to identify the conserved motifs of full-length SmGRAS proteins (Fig. 6B, C). We showed the 10 most conserved motifs located in the GRAS domain in  Expression analysis of SmGRAS genes in response to GA treatment As the key regulator of the GA signaling pathway, the transcriptional levels of GRASs were significantly affected by GA. We comprehensively analyzed the transcription level changes of these 35 SmGRAS family genes after 100 µm GA treatment for 2 hours. The expression levels of these SmGRAS genes changed a lot. The heatmap showed that 15 SmGRAS genes were GA induced, and 11 SmGRAS genes were GA repressed, while the other 9 genes did not change significantly under GA treatment (Fig. 7). The most significantly increased of these was SmGRAS5, which increased 4-fold expression levels and was followed with SmGRAS20 (2.3-fold) and SmGRAS14 (2.1-fold). The most significantly reduced of these genes was SmGRAS28, which reduced about 90% of the control. The expressions of SmGRAS31, SmGRAS8, SmGRAS11, and SmGRAS12 were all fell by more than half.

Discussion
It is well known that bioactive GAs are diterpene phytohormones that regulate growth and development throughout the whole life cycle [14]. There are many reports on the important roles of GA in plant growth, development and stress [10][11][12][13], but few reports on the relationship between GA and secondary metabolism. In addition, GA shares the same biosynthetic pathway and precursor substances with diterpenoid metabolites tanshinones, which were the main secondary metabolites of S. miltiorrhiza [9]. We speculated that there might have some correlation between GA and tanshinones. Therefore, we treated the hairy roots of S. miltiorrhiza with GA and found that it not only promoted the root growth but also increased the accumulation of tanshinones and phenolic acids. These results make us begin to pay attention to the specific mechanism of GA regulating secondary metabolism.
GA regulates many growth and development processes such as seed germination, root and shoot elongation, metabolism, stress, flowering and fruit patterning [10,12]. The GA signaling pathway includes the biosynthesis of the active GA, perception, signal transduction and inactivation [12]. It is now clear that GAs accumulate in the elongating endodermal cells of Arabidopsis root, and play central roles in growth regulation through key transcription factor DELLA (GRAS family member) [10,13]. Several genes of the GRAS transcription factor family in the GA signaling pathway have been reported to regulate the development of the cortex in the root. For instance, GRAS family genes form a complex SHR-SCR-SCL3 to regulate middle cortex formation in the Arabidopsis root [24]. SCL28 also plays a role in the root growth response to stress-induced microtubule organization in Arabidopsis [25]. Moreover, there is an interaction between energy metabolism and the GAmediated control of growth that coordinates cell wall extension, lipid metabolism and secondary metabolism in Arabidopsis [26]. Similarly, the overexpression of HaGRASL reduces the metabolic flow of GAs in Arabidopsis, and this modification could be relevant in axillary meristem development [27]. After treating the wild-type hairy roots of S. miltiorrhiza with GA, we found that GA promoted the root growth, promoted the accumulation of tanshinones and phenolic acids, but inhibited the contents of total phenolics and total flavonoids. In order to explore the possible molecular mechanism, we performed transcriptome analysis. The results showed that the gene response of the secondary metabolic pathway was significant after GA treatment. We speculated that GA regulated the accumulation of secondary metabolites by regulating genes in the secondary metabolites biosynthetic pathway.
After further analysis of the DEGs in secondary metabolic pathways, we found that most of the DEGs in MVA, shikimate and phenylpropane pathways are GA induced at the transcriptional level. This result also explains a series of composition content changes in the hairy root of S. miltiorrhiza at the transcriptional level after GA treatment, such as increased tanshinones and phenolic acids contents. However, there are too many components in the content of total phenolics and total flavonoids, so it is difficult to find the corresponding law through the changes of these DEGs. In addition, we also deeply analyzed the expressions of DEGs in GA signaling pathway.
The results showed that the gene expression levels of most key regulatory factors GRAS family genes and GA receptor GID1 in the pathway were significantly increased. Therefore, it is speculated that GA induction affects the GA signaling pathway and causes the response of key regulatory factors GRASs in the GA signaling pathway, and further regulating the downstream physiological processes such as growth, secondary metabolism, and stress response.
In order to study the roles of GRAS in GA regulation of growth and secondary metabolism, we used HMMER to search for all SmGRAS protein sequences based on the HMM profiles from the S. miltiorrhiza genome database.
SHR, SCR and SCL3 formed a complex to regulate middle cortex formation. DELLA has a central role in suppressing GA signaling [28]. DLT was involved in the brassinosteroid signaling, SCL4/7 was involved in the response to environmental stresses and LAS participated in the formation of lateral shoots [15]. Therefore, we speculated that SmGRASs might also be involved in the regulation of the above processes in S. miltiorrhiza. In

Plant materials and treatment
The leaves of S. miltiorrhiza were from the sterile seedling line preserved in our laboratory at Northwest A&F University in Yangling, Shaanxi Province, China. The establishment of S. miltiorrhiza hairy roots was derived from aseptic leaves of S. miltiorrhiza infected with Agrobacterium rhizogenes (ATCC15834), as previously reported [29]. Samples of the hairy roots (0.3 g fresh weight) were cultured in 100 ml beaker flasks contains 50 ml of the 6,7-V liquid medium on an orbital shaker 120 rpm·min − 1 , 25℃ in the dark and sub-cultured every 30 days.
For GA treatment, GA 3 (Sigma, CA, USA) stock solution was filter sterilized through 0.22 µm filters and added to cultures to the 21-day-old hairy roots to a desired final concentration of 100 µM. After 2 h and 6 days of treatment with GA, hairy roots were collected for the determination of qRT-PCR and HPLC analysis. Hairy roots without GA treatment were used as controls. All treatments were performed in three independent biological replicates.

HPLC Analysis
The 21-day-old hairy roots were treated with 100 µM GA 3 for 6 days and collected in three biological replicates.
0.04 g of the dried hairy roots was powdered and extracted by 8 mL of 70% methanol overnight and then sonicating the sample for 45 min. The mixture was centrifuged at 8000 g for 10 min, and the supernatant was filtered through a 0.2-µm filter and analyzed by HPLC, according to the general method in our laboratory that was described previously [30]. The GA 3 concentration analysis was measured by HPLC as previously described [31].

Determination of total phenols and total flavonoids
Total phenolics and total flavonoids were detected as previously described with some modifications [5,32]. The absorbance of the samples for the total phenolics and total flavonoids analyses was measured at 765 nm and 506 nm, respectively. Gallic acid and rutin (Solarbio, China) were used to construct a calibration curve to determine the total phenolics and total flavonoids contents, respectively.

RNA-seq Library Construction And Sequencing
At 2 h after GA treatment, GA-treated hairy roots (GA) and controls (CK) were collected from three biological replicates and analyzed by transcriptome technology. Total RNA was extracted using the RNA extraction kit following the manufacturer's instruction (Tiangen Biotech, Beijing, China). After measuring the quality, the strandspecific RNA-Seq libraries were constructed and sequenced on the Illumina PE150 platform (Novogene, Tianjin, China). The high-quality clean data were calculated and used for downstream analysis.

Identification Of DEGs And Functional Enrichment
The reference genome and gene model annotation files were downloaded from genome websites [33]. The fragments per kilobase of transcript per million mapped reads (FPKM) were used to determine the relative expression levels of each gene. DEGs were determined using the DESeq R package [34]. The expression levels of DEGs were considered significantly differentially expressed genes with an adjusted P-value < 0.05 and |fold-change|≥2. The GO term enrichment of DEGs was evaluated using the GOseq R package [35]. The statistical enrichment of DEGs in the KEGG was identified with KOBAS [36]. Mapman program was used to analyze the transcriptome data of metabolic and signal pathways.

qRT-PCR Validation
A total of 10 DEGs were randomly chosen to verify the RNA-seq data (Fig. S1). The 10 DEGs specific primers were designed by primer 5 (Table S1). qRT-PCR was performed on a real-  [39].

Analysis of conserved motifs and gene structures
The Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) was used to obtain the gene structure of introns and exons based on the CDS and corresponding to genic sequences. The GRAS domains were aligned and the conserved sites were checked manually for their corresponding amino acid residues, which were shaded using DNAMAN software. The conserved motif analysis of SmGRASs was performed by online MEME tools with 20 motif numbers. The theoretical isoelectric point (pI) and molecular weight (Mw) were predicted by the ExPASy server

Ethics approcal and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The data setssupporting the conclusions of this article are included with in the article and its additional files.

Competing interests
The authors declare that they have no competing interests.

Funding
This work was supported by the National Natural Science Foundation of China (No. 81773835) for the design and analysis the data of the study.

Authers' contributions
WL and ZL conceived and designed the experiment; WL, CL, JL, and ZB performed the experiments; WL analyzed the data and wrote the article with input from other authors. All authors read and approved the final manuscript.  DEGs of secondary metabolism pathway between CK and GA treatment hairy roots. DEGs marked in red indicated they were GA induced, and blue ones indicated they were GA repressed.

Figure 4
A model for a possible mechanism of the regulation of GA to root growth and secondary metabolism.
DEGs marked in red indicated they were GA induced, and blue ones indicated they were GA repressed.  Gene expression profiles of SmGRAS members under control and GA treatment. Genes marked in red indicated they were GA induced, and blue ones indicated they were GA repressed.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. Table S3.xls Figure S1.tif Table S2.xlsx Table S4.xls Figure S2.tif Table S1.xlsx Table S5.docx