GA3 regulates LR formation in cucumber
To investigate the impact of different GA3 levels during LR development, germinated seeds were initially subjected to exogenous GA3 at concentrations of 10 nM, 50 nM, 100 nM, 1 µM, 10 µM, 50 µM, 100 µM and 200 µM for five days (Additional file 1). We observed that GA3 concentrations of 10 nM, 50 nM, 100 nM, 1 µM, 10 µM, 50 µM and 100 µM had little effect on primary root length, respectively (Fig. 1). A significant reduction of primary root length was only observed for GA3 concentration at 200 µM. However, the number of secondary roots displayed no significant difference among the treatments of the eight GA3 concentrations and Control. Interestingly, the secondary root length was gradually increased by an increase in GA3 concentration in the range of 10 nM-100 µM, whereas at high GA3 concentration (200 µM), the secondary root length did not significantly differ from that of the control seedlings (Fig. 1; Additional file 1). Similarly, the effects of GA3 on the number and length of tertiary roots was also found to be dose dependent (Fig. 1). The tertiary root number and length were both significantly increased under 1 µM, 10 µM, 50 µM and 100 µM GA3 treatment, whereas the tertiary root number and length under higher doses of GA3 (200 µM) treatment didn’t show a significant difference compared with control (Fig. 1). We also investigated the effects of various levels of GA3 on the dry weight of roots and the ratio of root dry weight and shoot dry weight (RDW/SDW). The dry weight of roots was significantly increased by 50 µM GA3 treatment, which also increased the RDW/SDW (Fig. 1).
We also applied uniconazole (Uni, GA biosynthesis inhibitor of ent-Kaurene oxidase) to germinated cucumber seeds for five days. Low concentrations of Uni (10 nM, 50 nM and 100 nM) had no significant effects on the primary root length, secondary root length and number (Additional file 2). However, they were significantly suppressed by Uni and gradually decreased with the increase of Uni concentration from 1 µM to 200 µM (Additional file 2). The promotion and inhibition of LR development by exogenous GA3 and Uni, respectively, indicated that GA plays important roles in regulating LR development in cucumber and we used the 50 µM GA3 concentrations for further studies.
Overview of the transcriptome profiles of the roots in response to GA3 treatment
To better understand how GA regulates LR development and identify genes involved in the pathway in cucumber, we performed RNA-seq analysis using RNA isolated from the LRs at three developmental timepoints (2, 3 and 5 DAG) with or without 50 µM GA3 treatment. All tissues were analyzed in two or three independent biological replicates (17 samples in total). Summary statistics for each of the RNA-seq libraries are shown in Additional file 3. The Spearman correlation coefficient (SCC) analysis showed high reproducibility between the biological replicates, ranging from 0.97 to 0.99 (Fig. 2A). A principal component analysis (PCA) was also performed to obtain a general view of the transcriptomes derived from all 17 tissue samples. The first two principal components (PC1 and PC2) accounted for 70% of the total variance and the samples clustered based on both treatment and developmental timepoints (Fig. 2B). That is, there are different transcript profiles between the GA3-treatment and Control dependent on developmental timepoints.
Differentially expressed genes (DEGs) in response to GA3 treatment
To investigate how transcriptome changes of roots across the three developmental timepoints between GA3-treatment and Control, DEGs were identified by three pairwise comparisons at each developmental timepoint (Fig. 3). A total of 417 unique DEGs exhibited significant lower expression and 447 unique DEGs exhibited significant higher expression at the three timepoints of root development in GA3-treated roots as compared to the Control (Fig. 3A; Additional file 4). Notably, the number of DEGs gradually increased with the development of roots and the largest number of genes (724) exhibited differential expression at 5 DAG followed by 3 DAG (367) and 2 DAG (239) between the Control and GA3-treated roots, indicating more transcriptional changes in gene expression at the later developmental stages.
Various biological process GO terms showing significant enrichment were detected among the genes that were up- and down-regulated in the GA3-treated roots. Notably, GO terms related to cell wall organization/ modification, including ‘Cell wall organization’, ‘Cell wall organization or biogenesis’ and ‘Cell wall modification’, and root cell development/ differentiation, including ‘Root hair cell development’, ‘Root epidermal cell differentiation’ and ‘Lateral root branching’, are significantly enriched in the up-regulated DEGs in the GA3-treated roots (Fig. 3B). This result indicates that GA3 promotes LR development might be associated with the up-regulation cell wall-related genes. In addition, the over-represented GO terms in the down-regulated DEGs in the GA3-treated roots included terpenoid- and gibberellin-related terms, such as ‘Terpenoid biosynthetic process’, ‘Gibberellin mediated signaling pathway’ and ‘Gibberellin metabolic process’, suggesting that exogenous GA3 might inhibits the synthesis and signaling transduction of endogenous GAs.
To systematically explore the RNA-seq data, linear factorial modeling was also applied to identify DEGs significantly affected by GA3 treatment (GT) and the interaction of GT and developmental timepoint (GT*D). Genes that consistently showed higher or lower expression in GA3-treated roots across all developmental timepoints were identified as DEGs significantly affected by GT, while genes responding differently to GT and developmental timepoints were considered as DEGs significantly affected by GT*D. A total of 3397 and 31 genes were identified with significant GT and GT*D effects, respectively (Additional file 5). All DEGs identified with pairwise comparisons and linear factorial modeling were pooled, resulting in 3523 non-redundant DEGs (Additional file 6) that were used for the downstream analyses.
To further determine the expression patterns of the 3523 DEGs, the application of the degPatterns function on the dataset partitioned the DEGs into clusters of co-expressed genes that share similar expression dynamics. Of these, 3321 DEGs were clearly grouped into 20 clusters named C1- C20 (Fig. 4; Additional file 7). The GO enrichment analysis was performed for the genes in each cluster to investigate the biological significance of the distinct dynamic expression patterns. While no significant enriched biological processes were identified for the genes in C1, C3, C4, C5, C6, C8, C9, C11, C13, C14, C15 and C20, the genes in C2, C7, C10, C12, C16, C17, C18 and C19 were enriched for many interesting biological processes GO terms (Additional file 8). Of these, C2, C12, C18 and C19 show genes that are constitutively up-regulated by GA3 treatment. In C12, DEGs that showed increased expression in GA3-treated roots throughout the time course were enriched for lignin-related processes and ‘phenylpropanoid catabolic process’. The DEGs in C18 and C19 were enriched for hormone-related pathways and cell wall related processes. The opposite expression pattern was found in C16 and C17 showing decreased expression in GA3-treated roots at the three timepoints. C16 and C17 additionally showed that these genes are increased and decreased as the development of roots in both GA3-treatment and Control, respectively. The DEGs in C16 enriched for hydrogen peroxide-related processes and four processes related to triterpenoid metabolism, including ‘triterpenoid biosynthetic process’ and ‘pentacyclic triterpenoid biosynthetic process’. The significantly enriched GO terms for C17 included many metabolic processes, such as oxoacids, organic acids and fatty acids. In addition, the genes in C7 and C10 are up-regulated in response to GA3 treatment at 3 and 5 DAG. Top enriched GO terms in C7 and C10 were related with various biosynthetic/ metabolic processes, including phenylpropanoid, ubiquinone, quinone, lignin, aminoglycan, chitin, amino sugar, glucosamine-containing compounds and cell wall macromolecule.
Differentially expressed transcription factors (TFs)
Since TFs are the main regulators of gene expression, we sought out the differentially expressed TFs from the DEG list. There are 259 TFs, accounting for ~ 7.3% of the total DEGs, which were classified into 38 families (Additional file 9). The most abundant TF families were shown in Fig. 5, including bHLH (28), ERF (25), C2H2 (19), GRAS (17), WRKY (17) and MYB (15). Interestingly, 248 out of 259 (95.8%) TFs were significantly affected by GT (Additional file 6; Additional file 9). Of these, 73 TFs were also identified in pairwise comparisons (Additional file 6), including one ETHYLENE RESPONSE FACTOR (ERF) TF, Csa3G895680, significantly affected by GT*D interaction and putatively associated with primary cell wall deposition .
Many orthologs of the TFs were known from other studies to participate in LR development. For example, the putative cucumber CYTOKININ RESPONSE FACTOR 2 (CRF2) (Csa4G051360), encoding an ERF protein, was clustered into C17 and down-regulated in GA3-treated roots (Fig. 4; Additional file 7; Additional file 9). In Arabidopsis, CRF2 has been proposed to control LR initiation and LR formation . Nine out of 17 differentially expressed GRASs were down-regulated in GA3-treated roots (Fig. 4; Additional file 7; Additional file 9), including genes (Csa1G616330 and Csa5G623560) putatively regulating GAs response pathway in root endodermis  and Csa4G061850 putatively involved in endodermal fate in the root meristem . In addition, Csa6G495620 and Csa7G372280 is a closest ortholog of Arabidopsis GRAS TF SHORTROOT (SHR, AT4G37650) and SCARECROW (SCR, AT3G54220), respectively, which regulates both stem cell niche specification and radial patterning in the root [54–56].
Notably, we found that some TFs in one gene family were all up-regulated or down-regulated by GA3-treatment. For example, the seven AP2 genes grouped in C2, C12, C18 and C19 were up-regulated in response to GA3 treatment at 2, 3 and 5 DAG (Fig. 4; Additional file 7; Additional file 9). Recently, AP2 TFs were demonstrated to be involved in LR development in Arabidopsis . Auxin response factors (ARFs) play important roles in auxin signaling pathways that are crucial for root branching and LR formation in plants [6, 58]. The five DEGs belonging to ARF gene family were assigned to C4, C12, C18 and C19 and significantly up-regulated by GA3-treatment (Fig. 4; Additional file 7; Additional file 9). Of these, Csa6G405890 is the closest ortholog of Arabidopsis ARF16 (AT4G30080), which is indispensable for root cap development . Csa2G068680 and Csa2G315390 are one of the orthologs of Arabidopsis ARF6 (AT1G30330) and ARF17 (AT1G77850), respectively. In Arabidopsis, ARF6 positively regulates adventitious roots formation, whereas ARF17 is an inhibitor of the formation of adventitious roots [60, 61]. In addition, GROWTH-REGULATING FACTORS (GRFs) are well known as positive regulators of cell proliferation and expansion and play important roles in regulating root development [62–64]. In our datasets, all five GRF TFs showed consistently decreased expression in response to GA3-treatment in C1, C3 and C15 (Fig. 4; Additional file 7; Additional file 9).
These results suggest that the differentially expressed TFs identified in our datasets may play a key role in regulating GA-mediated LR development in cucumber. However, further study is required to dissect their exact roles in GA-mediated pathways regulating LR development.
DEGs involved in GAs synthesis and signaling pathway
Since GA3 promotes LR development in our study, we further investigated the expression of genes involved in GAs synthesis and signaling pathways (Fig. 6). CPS, KS and KO, which were each encoded by a single gene in Arabidopsis, catalyze the early stages of GAs biosynthesis pathway [65, 66]. Similarly, we also identified one gene encoding CPS, KS and KO, respectively, in cucumber genome. Cucumber CPS (CsCPS) and KO (CsKO) are DEGs with significant GT effects and down-regulated in GA3-treated roots (Fig. 6; Additional file 5). Two KAO genes, CsKAO1 and CsKAO2, were identified in cucumber. The expression of CsKAO2 was repressed by exogenous GA3 treatment, whereas CsKAO1 showed very low expression in our datasets. GA20ox is the key enzyme of bioactive GA synthesis, and cucumber has five GA20ox genes . The transcripts of CsGA20ox1 and CsGA20ox2 were detected in our datasets, and both of them were DEGs with significant GT effects, showing lower expression levels in GA3-treated roots than that of Control. Among the five cucumber GA3ox genes , only CsGA3ox1 expression was detected in both conditions, showing a significant decrease after GA3-treatment throughout the time-course. So far, GA7ox was only identified in pumpkin and cucumber but has not been found in other species, which oxidizes GA12-aldehyde to GA12 and possess mono-oxygenase 7-oxidase activity [68–70]. The expression of the two genes encoding GA7ox also exhibited decreased expression in response to GA3-treatment. GA2ox can inactivate bioactive GAs through catalyzing 2b-hydroxylation . Cucumber GA2ox is putatively encoded by eight genes (CsGA2ox1-CsGA2ox8) , but only CsGA2ox1, CsGA2ox4 and CsGA2ox6 expression were detected in our datasets. Of these, CsGA2ox1 and CsGA2ox4 were DEGs significantly affected by GA3-treatment and showed consistently higher expression in GA3-treated roots compared with Control, whereas the expression of CsGA2ox6 was not significantly affected by GA3-treatment.
We also investigated the expression profiles of genes involved in GA signaling pathways, including genes encoding GA receptor GIBBERELLIN INSENSITIVE DWARF1 (GID1), a positive regulator of GA signaling F-box protein SLEEPY1 (SLY1) as well as GA signaling repressors, including DELLA proteins (DELLAs), and SPINDLY (SPY) [65, 71]. The expression of CsGID1A, CsGID1B and CsSLY1 were significantly repressed by GA3-treatment (Fig. 6). For the GA signaling repressors, two genes (CsGAI2 and CsGAI3) encoding DELLA proteins exhibited up-regulation in response to GA3-treatment, whereas the expression of CsGAI1 and CsSPY were not affected by GA3-treatment. Overall, 20 out of 34 genes in GA synthesis and signaling pathways were expressed in cucumber roots. Of these, 15 genes were identified as DEGs with significant GT effects, suggesting their important roles in GA-mediated LR development.
DEGs involved in auxin synthesis and signaling pathway
Considering the central roles of auxin in almost all steps of LR development and the fact that GAs influence auxin biosynthesis and signaling during root development [6, 72], we further investigated the expression of genes involved in auxin metabolisms, including IAA biosynthesis, inactivation, transport and signal transduction.
The TRYPTOPHAN AMINO-TRANSFERASE OF ARABIDOPSIS/ YUCCA (TAA/YUC) pathway has been proposed as the primary endogenous auxin biosynthesis pathway that involves in various biological processes, including root development [73–75]. There are eight YUC family genes in cucumber (Additional file 10). Of these, Csa3G133910 (CsYUCCA3) and Csa3G619930 (CsYUCCA5) were identified as DEGs with significant GT effects. They were clustered in C19 and significantly up-regulated under GA3-treatment (Fig. 7; Additional file 7). To our surprise, no TAAs were identified in cucumber genome using BLASTP, suggesting that the auxin biosynthesis pathways in cucumber remain elusive. We then checked the expression of the genes involved in IAA inactivation, including genes encoding IAA-amido synthase GRETCHEN HAGEN 3 (GH3), Dioxygenase for Auxin Oxidation (DAO) and IAACARBOXYLMETHYLTRANSFERASE 1 (IAMT1). GH3 is putatively encoded by eight genes in cucumber (Additional file 10), by which IAA can also be conjugated to amino acids . Two out of eight GH3 genes, Csa3G198490 (CsGH3.3) and Csa3G431430 (CsGH3.5), were down-regulated in response to GA3, whereas only the expression of Csa3G088930 (CsGH3.4) was significantly increased in GA3-treated roots (Fig. 7). The putative cucumber IAMT1 (Csa7G081680) showed similar expression levels in the roots of GA3-treatment and Control. In addition, there was not an ortholog of Arabidopsis DAO in cucumber. These results suggested that exogenous GA3 might suppress the processes converting active IAA to inactive forms.
PINFORMED (PIN) family, AUXIN1/LIKE-AUX (AUX/LAX) family and ATP-binding cassette family B (ABCB) family are well known auxin transport proteins in plants [77, 78]. In cucumber, seven, seven and 20 PIN, AUX/LAX and ABCB family members were identified, respectively (Additional file 10). Of these, two PINs (CsPIN1 and CsPIN5), four AUX/LAXs (CsAUX1, CsAUX2, CsLAX1 and CsLAX4) and five ABCBs (CsABCB4, CsABCB5, CsABCB10, CsABCB11 and CsABCB18) were DEGs with significant GT effects (Additional file 5). Notably, only CsLAX1 and CsABCB18 were down-regulated by GA3-treatment. The well-known components of the auxin signal transduction pathway include TRANSPORT INHIBITOR RESPONSE1/AUXIN-RELATED F-BOX (TIR1/AFB) and Aux/IAA proteins and the AUXIN RESPONSE FACTORs (ARFs) [78, 79]. Five, 27 and 16 genes were identified in cucumber TIR1/AFB family, Aux/IAA family and ARF family, respectively (Additional file 10). Notably, a total of 21 genes from the three families were significantly affected by GA3 treatment and 20 of them were constitutively higher expressed in GA3-treated roots compared to Control (Fig. 7).
Taken together, a total of 36 DEGs involved in IAA biosynthesis, inactivation, transport and signal transduction were identified in our datasets. Among them, 31 genes were up-regulated in response to exogenous GA3. To determine if GA3 treatment enhances IAA biosynthesis, we analyzed IAA levels in 5 DAG roots by LC-MS/MS analysis. As shown in Fig. 7C, IAA level was significantly increased in GA3-treated roots. These results suggest that the enhanced auxin metabolisms are associated with GA-mediated root development.
DEGs involved in cell wall biosynthesis
Given that the up-regulated DEGs were enriched for various cell wall-related processes, including metabolic processes of lignin, chitin and cell wall macromolecule which are the primary components of plant cell wall , we also want to know the transcriptional changes of cell wall biosynthesis genes under GA3-treatment in cucumber roots. The Arabidopsis cell wall biosynthesis genes include 10 cellulose synthase (CESA) genes, 30 cellulose synthase- like (CSL) genes (CSLAs, CSLBs, CSLCs, CSLDs, CSLEs and CSLGs) and 12 Glucan Synthase-like (GSL) genes [81, 82] (https://www.arabidopsis.org/browse/genefamily/cellwall.jsp). A total of 47 genes from CESA, CSL and GSL families were identified as putative cell wall biosynthesis genes in cucumber genome. The 47 genes encode eight CESA proteins, 30 CSL proteins and 9 GSL proteins in eight distinct groups (Fig. 8A).
The expression of four CESA, 10 CSL and one GSL genes were significantly affected by GA3-treatment and 11 of them showed higher expression in GA3-treated roots (Fig. 8B), indicating cell wall biosynthesis in roots was likely promoted by GA3-treatment which is consistent with the phenotype of GA3-treated roots (Fig. 1; Additional file 1). Some of them are putatively involved in root development based on previous studies. For example, CsCSLD2 is the closest ortholog of Arabidopsis CSLD3, which is essential for the synthesis of cell wall at the root hair tip [83, 84]. In Arabidopsis, loss-of-function of CSLA9, the closest ortholog of CsCSLA5, had reduced numbers of LR .