RNA-Seq reveals new DELLA targets and regulation in transgenic GA-insensitive grapevines

Background Gibberellins (GAs) and their regulator DELLA are involved in many aspects of plant growth and development and most of our current knowledge in the DELLA-facilitated GA signaling was obtained from the studies of annual species. To understand GA-DELLA signaling in perennial species, we created ten GA-insensitive transgenic grapevines carrying a DELLA mutant allele (Vvgai1) in the background of Vitis vinifera ‘Thompson Seedless’ and conducted comprehensive analysis of their RNA expression profiles in the shoot, leaf and root tissues. Results The transgenic lines showed varying degrees of dwarf stature and other typical DELLA mutant phenotypes tightly correlated with the levels of Vvgai1 expression. A large number of differentially expressed genes (DEGs) were identified in the shoot, leaf and root tissues of the transgenic lines and these DEGs were involved in diverse biological processes; many of the DEGs showed strong tissue specificity and about 30% them carried a DELLA motif. We further discovered unexpected expression patterns of several key flowering induction genes VvCO, VvCOL1 and VvTFL1. Conclusions Our results not only confirmed many previous DELLA study findings in annual species, but also revealed new DELLA targets and responses in grapevine, including the roles of homeodomain transcription factors as potential co-regulators with DELLA in controlling the development of grapevine which uniquely possess both vegetative and reproductive meristems at the same time. The contrasting responses of some key flowering induction pathway genes provides new insights into the divergence of GA-DELLA regulations between annual and perennial species in GA-DELLA signaling. Electronic supplementary material The online version of this article (10.1186/s12870-019-1675-4) contains supplementary material, which is available to authorized users.


Background
Gibberellins arguably has the widest role to play in the whole life cycle of a plant. Plants tightly regulate the availability of GA by coordinating the activation and deactivation of multiple biosynthesis genes of GA in different tissues at various developmental stages [1]. The flow rate and availability level of bioactive GA in various tissues are regulated by a de-repressible system in which GA interacts with its receptor GA INSENSITIVE DWARF1 (GID1) and negative regulator, DELLA, to form a trimeric GA:GID1:-DELLA complex [2]. As a key component in the GA signaling cascade, DELLA genes, such as Rht1 from wheat and GIBBERELLIN INSENSITIVE (GAI) from Arabidopsis, have been extensively studied [3]. These genes all carry a conserved DELLA amino-acid domain in their protein sequences. A critical gain-of-function mutation in the DELLA domain prevents the formation of a stable trimeric complex with GA:GID1 required for its eventual degradation, and thus instead maintains the negative regulation of GA signaling [4]. As a result, GA-deficient plants show many abnormal developmental characteristics such as smaller darker leaves and reduced internode length [5]. In addition to basic characterization of DELLA protein genes in the GA signal transduction [6], many DELLA target genes which control development of various traits in model and annual species, have also been identified [7]. In contrast, progress in understanding DELLA regulation in perennial woody plants lags considerably behind, although some DELLA knowledge has been obtained from transgenic poplar research; metabolic and phenotypic changes in transgenic GA-deficient poplar [8] were found similar to their Arabidopsis mutant counterpart [9]. Similar to Arabidopsis gai mutants, transgenic poplar carrying non-degradable DELLA gene (i.e. gai, a mutant version of GAI) were dwarf but insensitive to GA-mediated restoration to wildtype [10,11].
Grapevine (Vitis) is a perennial woody species with significant economic value. Impact of DELLA-mediated GA on vine growth and development was first found in a natural DELLA mutant derived from a somatic variation in the L1 meristem layer of the grape cultivar V. vinifera 'Pinot Meunier' [12]. The mutant vine carries a copy each of the wildtype GIBBERELLIC ACID-INSENSITIVE1 (VvGAI1) and mutant (Vvgai1) and manifests reduced internode length and short plant stature, which is a hallmark trait of DELLA mutants in many species. The mutant vine also shows continuous flowering and bears inflorescences at most nodes after the first inflorescence appears. Apparently, both GA and DELLA proteins are involved in grape's inflorescence formation, flower induction and berry development [13].
To elucidate DELLA-mediated roles in grapevine development, we created transgenic lines in the background of a table grape cultivar, V. vinifera 'Thompson Seedless' , carrying a grape mutant DELLA gene, noted as Vvgai1 [14]. We analyzed the RNA-Seq profiles of shoot, leaf and root tissues of non-transgenic and representative transgenic lines and affirmed in grapevine the presence of a DELLA-centered feedback mechanism that maintains the GA homeostasis [15] and also the intricate interactions of DELLAs with numerous transcription factors controlling plant development and growth as were found in annual species [16]. We further discovered DELLA's possible roles in the induction of the anlagen, the distinct vegetative meristem for tendril/ inflorescence development in grapevine [17], through coordination with meristem regulators. Moreover, we discovered unexpected expression patterns of several key flowering induction genes, including grapevine CON-STANS (VvCO), CONSTANS 1 (VvCO1), and TER-MINAL FLOWER 1 (VvTFL1 or CENTRORADIALIS), which were in sharp contrast with what were found in annual species. These findings provide insights into how DELLA genes regulate grapevine development and growth, especially in relation to flowering, and fill some critical knowledge gaps in this important research area between annual and perennial species.

Results
Transgenic Vvgai1 expression and phenotypes Ten Vvgai1 transgenic vines were generated in this study and five representative ones covering a range of mutant phenotypic variations were chosen for RNA-Seq analysis. The internode length of four transgenic lines were significantly shorter than the NT in various degrees, with approximate three-fold reduction for the most severe transgenic line G02 ( Fig. 1a and b, Table 1). More severe dwarf lines had smaller shoots and leaves that were darker and curled at the edges (Fig. 1c). The distinctive V. vinifera shoot patterning of two consecutive nodes with a tendril followed by node without a tendril [18] was proportionally disrupted with a higher frequency of tendril skips among the shorter dwarf lines. The tendrils also appeared progressively late among the severe dwarf vines, as late as the tenth node in G03 compared to the fifth in the NT (Fig. 1d, Table 1). No tendrils were observed for the most severe G02 line even after 1 year growth in greenhouse conditions (data not shown). More severe dwarf lines had proportionally enlarged and denser roots (Fig. 1e, f ) and showed poorer root establishment.
As expected, Vvgai1 expression was found in all the five transgenic lines but not in the NT using digital qRT-PCR. Among the five transgenic lines, the most severe dwarf G02 showed statistically the highest Vvgai1 expression in all three tissues (Table 1, Fig. 1a). These results suggest that there was a positive correspondence between the relative Vvgai1 expression levels (Vvgai1/ VvGAI1 ratios) and the observed severity of phenotypic changes among the transgenic lines.
Differentially expressed genes associated with transgenic Vvgai1 expression Differential gene expression analyses in the leaf, shoot and root tissues were separately conducted for the five transgenic lines. More severe dwarf lines generally have more DEGs, although such trend was not observed in the root tissue (Fig. 2, Additional file 1: Table S1). To capture maximum numbers of DEGs associated with Vvgai1 with the least false positives, we opted to retain only those DEGs (1.5x-fold changes at FDR ≤ 0.01) in the most severe transgenic line G02 and concurred in at least one of the four remaining transgenic lines. We also identified 25 genes originally filtered out for not having the minimum reads (i.e. 1CPM) in all the 12 libraries, but were in fact expressed only in the transgenic lines and suppressed in the NT, or vice-versa (Additional file 2: Table S2). This recovered set had several cell-wall related genes induced only in the leaf of the transgenic lines and a key flowering gene SHORT VEGETATIVE PHASE (VvSVP) highly suppressed in the root of the transgenic lines. The final numbers of DEGs for subsequent analyses were 153, 719, and 2314 for shoot, leaf and root, respectively, totaling to 2986 genes significantly differentially expressed at least once in the three sampled tissues (Fig. 2).
We validated the RNA-Seq data by analyzing the digital RT-PCR expression of 30 notable genes related to DELLA regulation from the same tissues. The RT-PCR and RNA-Seq datasets were in a good agreement with a general correlation coefficients (r) ranging from 0.77 to 1.0 (Additional file 3: Table S3).

Gene ontology (GO) analysis
The annotated DEGs in the shoot and leaf tissues were comprised of about twice as many suppressed (96 for the shoot and 494 for the leaf ) as induced ones (57 for the shoot and 225 for the leaf ), while the root had about the same (1124 suppressed and 1190 induced) (Fig. 2).  The enriched GO terms have large overlaps of the member genes, and a semantic clustering analysis revealed about ten common interdependent biological themes ( Table 2). This is not a surprise as the arbitrary process of clustering involved a considerable number of master growth regulators with broad biological roles known to coordinate plant's morphogenetic changes and responses to stimuli, such as the transcription factors NAM/ ATAF1,2/CUC2 (NAC), WRKY, SQUAMOSA PRO-MOTER BINDING PROTEIN LIKE (SPL), TB1-CYC-PCF (TCP) and THREE-AMINO-ACID-LOOP-EX-TENSION (TALE) class of homeoproteins. There were also numerous biosynthesis and signaling genes belonging to hormones auxin, cytokinin, brassinosteroids, and ethylene.

Expression profiles of key GA signaling genes
Similar to the currently accepted model in Arabidopsis [6], the grapevine GA signaling are comprised of three DELLA homologues of VvDELLA1, 2 and 3, two GA receptor GID1 of VvGID1a and VvGID1b. Grapevine, however, have two proposed F-box SLY1 genes, VvSLY1a and VvSLY1b [19]. VvDELLA1 (synonymous to VvGAI1) in the NT was more abundant in the shoot and leaf (5 0 CPM), about five folds higher than that in the root (~10 CPM). VvDELLA1 in the transgenic G02 line, however, was differentially expressed with about two folds higher than the NT in the leaf and root, but not in the shoot (Fig. 3a). VvDELLA2 had the highest expression in the root among the three NT tissues (30-240 CPM) and its expression in the root, but not in the other two tissues, was increased more than twice in the transgenic G02 line. VvDELLA3 was consistently expressed at low levels in all three tissues in both NT and transgenic background (~10 CPM). The expression profiles of VvDELLA1 and VvDELLA2 were validated by RT-PCR with high correlation coefficients (r =~0.95, Additional file 3: Table S3). The endogenous VvDELLAs seemed responsive to and inducible by the transgenic Vvgai1.
The two homologues of GA receptor GID1 were abundantly expressed in all tissues in the NT. In the transgenic lines, VvGID1b showed five folds of expression increases in the leaf (~40 CPM to~220 CPM) (Fig. 3b). The two grapevine F-box genes, VvSLY1a and VvSLY1b, consistently showed low levels of expression (< 10 CPM) in both NT and transgenic background (Fig. 3c).
Homologues of the GA biosynthesis genes GA3 oxidase (GA3ox) and GA20 oxidase (GA20ox), another group of recipients of the DELLA-controlled feedback mechanism (like GID1) [20], were induced in all the three sampled tissues but at different levels, seemingly preserving the observed tissue-specificities of this multi-gene family [21]. VvGA20ox2a showed 15-, 2-and a dramatic 50 folds of increases in the shoot, leaf and root, respectively (Fig. 3d). In a similar manner, VvGA3ox1a was induced by more than four folds in the leaf, and 10 folds in the root (Fig. 3e). In consonance, the GA deactivating GA2ox genes were generally suppressed in the transgenic background. The largest suppression was observed for VvGA2ox1a, which decreased by about nine folds (from 28 CPM to 3 CPM) and two folds (76 CPM to 40 CPM) in the shoot and root, respectively. VvGA2ox8a was also suppressed with a  decrease of about 2.5 folds (55 CPM to 20 CPM) in the leaf (Fig. 3f ). These results affirm the positive feedback mechanism facilitated by DELLA to GID1 and the biosynthesis genes GA20ox and GA3ox [15]. Interestingly in this study, we observed that such feedback signal appeared tissue-specific, i.e. different GID1 and GA biosynthesis homologues induced in different tissues.
Expression profiles of some key genes influencing shoot, leaf, and root development A number of highly responsive DEGs observed are regulators of developmental processes. Among the highly suppressed were members of the GROWTH REGULAT-ING FACTOR (GRF) gene family, important regulators mostly found in intercalary meristem and leaf primordia [22]. For example, the homologue GRF5 controlling cell number related to leaf size [23] was suppressed by about 20 folds in the leaf of the G02 line (Fig. 4). The brassinosteroid biosynthesis gene DWARF4, another key gene synergistically associated with auxin [24] and GA regulation [25] and whose loss-of-function mutants were dwarfs with round and dark-green leaves in shortened petioles [26], was suppressed in all three tissues, but especially in the leaf of the G02 line by nine folds (Fig. 4). The TCP transcription factor, a master regulator that influences plant height, leaf curvature [27] and tendril phylloctatic patterning in grapevine [28], was suppressed at four to nine folds in the shoot and leaf of the transgenic G02 line (Fig. 4).
Genes for light-mediated developmental processes were notably induced. ERF105, a gene active during intense light and freezing stress conditions [29,30], was induced by seven and three folds in the leaf and root, respectively, and the PHYTOCHROME INTERACTING FACTOR 1 (PIF1), a major gene involved in chloroplast production through maintenance of cell number in the leaf [31], was increased by three folds in the shoot and leaf and about eight folds in the root (Fig. 4). The cell wall loosening gene EXPANSIN A1 was also induced about two folds in the shoot and six folds in the leaf and root of the transgenic background (Fig. 4).

Auxin and cytokinin
Auxin and cytokinin related genes have extensive cross-talk with GA-DELLA regulation as master regulators for many growth and development processes in plants [32]. Auxin-related genes seemed generally active in all transgenic lines. SHORT HYPOCOTYL 2 (SHY2), a negative regulator of auxin responses [33], was negatively correlated to Vvgai1 expression in the shoot and root with a suppressed response of about two and four folds, respectively, in the transgenic line G02 (Fig. 4). Conversely, the grape homologous gene of TRANSPORT INHIBITOR RESPONSE 1 (TIR1), a positive regulator of auxin [34], was only slightly induced in the shoot but increased by about two folds in the root. In conjunction, two homologues of the auxin biosynthesis gene YUCCA showed induced responses more prominently in the root. On the other hand, related auxin transporters and response factors exhibited various responses in the sampled tissues, reflecting the dynamic auxin flux. The grape homologue of PIN-FORMED3 (PIN3) was suppressed and significantly negatively correlated with Vvgai1 in their expression in the root, which was validated by our qRT-PCR (Additional file 3: Table S3). The grape homologue of AUXIN RESPONSE FACTOR11 (ARF11), one of the response factors antagonized by cytokinin when over-expressed in the root [35] was strongly induced. Interestingly, activities of cytokinin-related genes seemed unusually low in our transgenic lines. The  biosynthesis genes homologous to ISOPENTENYL-TRANSFERASE2 (IPT2) and IPT3 were of low expression levels (10 CPM) in the three sampled tissues. The grape CYTOKININ OXIDASE (CKX) gene family that irreversibly deactivate cytokinins [36] was suppressed in the root, although some were interestingly induced in the leaf (Fig. 4).

Homeodomain gene family
The expressions of many homeodomain transcription factors were generally suppressed and a large number of them belonged to the family classes of TALE, WUSCHEL RELATED HOMEOBOX (WOX) and HOMEODOMAIN--LEUCINE ZIPPER (HD-ZIP) (Additional file 1: Table S1). These homeodomain transcription factors are extensively involved in regulation of meristem maintenance and plant form [37], and manifest disruption of shoot phyllotaxy [38]. TALE has two subfamilies, KNOTTED1-like (KNOX) and BELLRINGER (BEL). In the KNOX subfamily, KNOT-TED-LIKE FROM Arabidopsis thaliana1 (KNAT1) and SHOOT MERISTEMLESS (STM) are a homologous pair that synergistically interact with cytokinin in the shoot apical meristem (SAM) [39]. We observed that STM was suppressed in the shoot and leaf of the transgenic line G02; KNAT1 was suppressed only in the shoot (Fig. 4); and KNAT1 was not expressed in the leaf which is in agreement with its negative role in leaf primordia initiation [40]. BEL1 (BELLRINGER 1), a member of the BEL subfamily which forms critical KNOX-BEL heterodimers for proper initiation and maintenance of the apical meristem [41], was induced in all three tissues in the transgenic line G02. WOX family genes are specialized transcription factors found during embryogenesis and lateral organ formation [42]. In transgenic line G02, WOX1 gene was induced in the shoot but suppressed in the leaf, while the WOX13 was suppressed in the root (Fig. 4). Furthermore, we found that the HD-ZIP family member genes were generally suppressed in the shoot and root tissues. HD-ZIP member genes are associated with organ and vascular development [43].

Flower-induction pathways
Four flowering-related genes showed significant responses to Vvgail expression in this study. VvSBP12 is a member of the miR156-targeted grapevine SPL genes [44] and its homologue in Arabidopsis, AtSPL3, was shown to be a target of suppression by DELLA and a target of miRNA156/ − 172 phase-change regulation [45]. The VvSBP12 was suppressed in all three tissues of the G02 line, especially in the shoot (Fig. 5). Interestingly, one of the homologues of a novel gene family associated with miRNA metabolism, HEN1 [46], showed variable response patterns among the three tissues, being induced in the shoot and leaf but strongly suppressed in the root (Additional file 1: Table S1).
On the other hand, the two CONSTANS homologues in grapevine, VvCOL1 and VvCO, were induced in all three tissues (Fig. 5), which is in contrast in Arabidopsis' where they were reported as targets of suppression by DELLA [47]. Interestingly, VvTFL1 gene, a negative regulator of meristem identity [48], was notably suppressed in the shoot but induced in the root (Fig. 5, Table 3, Additional file 1: Table S1). This intriguing expression profile was verified by the qRT-PCR results (r =~0.99) (Additional file 3: Table S3). Furthermore, we observed that the MADS-domain-carrying VvSVP gene was highly suppressed in all three tissues (Fig. 5). As predicted from the suppression status of the key flowering induction genes above, we found that some of the floral meristem identity genes were not induced. Although the floral meristem identity gene FLOWERING LOCUS T (FT) had no detectable transcripts in both NT and transgenic lines, SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) was expressed but not significantly differentially expressed. LFY was found expressed only in the shoot where it was likewise not differentially expressed. In grapevine, these meristem identity genes were expressed in the lateral branches after dormancy [49].

Transcription factor motif enrichment
Transcription factor enrichment analysis among the Vvgai1-correlated DEGs (r ≥ abs (0.70)) revealed a total of 17 significantly over-represented regulatory motifs (Additional file 4: Table S4). Of these the DELLA motif was the most over-represented in all three tissues: 45% of the suppressed DEGs in the shoot (49 out of 110), 37% of the suppressed DEGs in the leaf (269 out of 725) and 34% of the induced DEGs in the root (126 out of the 371). These amounted to 395 DELLA-bearing DEGs which were about 13% of the 2986 Vvgai1-correlated DEGs. To assess the extent of possible co-regulation of DELLA with other transcription factors, we examined these 395 DELLA-bearing DEGs for the abundance of other regulatory motifs. Computational motif prediction based on PlantTFDB 4.0 [50] revealed that aside from the DELLA motif, at least two of the 16 motifs discussed earlier were present in any of these DEGs (Table 3,  Additional file 5: Table S5) Table S5). The remaining motifs were of lower frequency and tend to be tissue-specific. The class I and II TCP motifs, for example, were found in about 4-10% of the DELLA-enriched DEGs, most of which were in the leaf. The R2R3 MYB motif, as represented by three genes (MYB15, MYB61, MYB88) involved in stomatal movement and lignification [51][52][53], was co-located with DELLA in about 4-60% of the DELLA-enriched DEGs in the shoot, 4-15% in the leaf, and 2-15% in the root. The C2H2 zinc finger motif were represented by two IDD genes (IDD5 and IDD11) and co-shared with the DELLA motif in 4-8% of the DEGs in the shoot, 5-8% in the leaf, and 2-17% in the root.

Discussion
Global responses to DELLA-mediated regulation were previously reported but most were on annual model species, for a specific tissue or a specific developmental stage [21,54,55]. In this study, we analyzed RNA-Seq profiles in the shoot, leaf, and root of five transgenic GA-insensitive lines whose severity of mutant phenotypes corresponded well with the levels of transgene Vvgai1 expression. By examining differential expression and correlation to Vvgai1 in three tissues, we were able to dissect both global and tissue-specific regulation of DELLA, one of its characteristic regulatory features [56]. Our study has not only confirmed many of the previously reported GA-and DELLA-related genes and GO processes, but also found new DELLA targets and action modes, including the involvement of meristem regulators in co-regulating vine development and identification of several key flowering induction genes which responded to Vvgai1 differently from annual species.

DELLA's impact on the expression of GA signaling genes
In agreement with the previous report [19], we found VvDELLA2 abundantly expressed in all three tissues of the NT, followed by VvDELLA1, while the expression of VvDELLA3 was generally low. Interestingly, the endogenous VvDELLAs were significantly induced, but  Table 1) varied from tissue to tissue, suggesting their regulatory interdependence and tissue-specific actions [19,57]. Equally interesting was the observation that the Vvgai1 (mutant DELLA) and induced endogenous DELLA transcripts in transgenic lines were not degraded as expected in a wildtype background. This was also in line with the observation of very low transcript levels of VvSLY1a and VvSLY1b (Fig. 3b), the two F-box genes that degrade DELLA through ubiquination [58]. Low SLY1 expression was previously found associated with unusually high accumulation of VvDELLA and VvGID1 transcripts [59]. We also noted that the putative DELLA-driven feedback induction of VvGID1 [15] was tissue-specific as only one of the two VvGID1 homologues (VvGID1b) exhibited a positively correlated change to Vvgai1 and prominently in the leaf tissue. On the other hand, different homologues of GA biosynthesis components (GA20 oxidases and GA3 oxidases) were induced in the same manner as in annual [12,21] and perennial species [10].

DELLA and co-regulation networks in grapevine
The 2986 DEGs identified in this study represent a large number of enriched GO terms which overlapped in the three tissues in various degrees and response states, reflecting not only the breadth of GA signaling cascades [55] but also the dynamic coordination of biological functions within and between tissues [1,60]. For example, one of the larger GO biological terms was related to meristem regulation and pattern-specification processes and largely enriched among the suppressed DEGs in the shoot and leaf. Included in the GO term were key meristem regulators of gene families such as GRFs, SPLs, SVPs, and genes such as YABBY, AS1, STM and KNAT1. In the root, however, the same GO term was enriched among the induced DEGs, including the auxin receptor TIR1 and the major root growth and patterning regulators, the GRAS SHORT ROOT (SHR), SCARE-CROW (SCR), and the transcription factor IDD5.
Misregulations of the meristem-related developmental processes in the transgenic vines was apparent all the more by the fact that numerous classes of homeodomain  transcription factors were generally in a suppressed state in the sampled tissues, including the class I KNOX transcription factors which interact with GA pathway genes in the maintenance of SAM [21,37,39]. In Arabidopsis, the class I KNOX genes STM and KNAT1 critically coordinate proper SAM growth by properly localizing the hormones cytokinin and GA [39,61,62]. Transcriptome scans of GA-insensitive Arabidopsis mutants seemed to suggest that both STM and KNAT1 were not particularly expressed at the early stage [15,55,63] as much as during the adult flowering stage [54,64]. Our discovery of a large number of key homeodomain transcription factors among the responsive DEGs at the early growth stages of the GA-insensitive grapevines highlights the unique organ development process in a grapevine shoota separate vegetative and reproductive meristem (i.e. anlagen) in the same apical meristemand the complex regulation of vegetative and reproductive phase transition it entails. Such complex regulation was also reflected by 17 large and diverse transcription factor motifs statistically over-represented among the DEGs. Of which, DELLA was the largest, found in about 30% of the total DEGs, followed by motifs of several meristem regulators such as TALE, TCP, MYB and IDD, and stress-response regulators, WRKY, NAC and bHLH. Many of these motifs were co-located with DELLA motifs in DEGs, providing opportunities for DELLA to interact with other transcription factors to co-regulate diverse biological processes [7,65]. Indeed, in the DELLA-motif-bearing DEGs, motifs of several important regulators were present in varying enrichment in three tissues. In the root, the motifs of AP2/ERF, MYB and IDD were relatively more enriched. The enrichment of IDD motif was notable because the IDD5 homologues were recently shown to act as a transcriptional scaffold in the competitive interaction between SCR and DELLA in mediating root growth [66]. On the other hand, in the suppressed DELLA-bearing DEGs of the shoot and leaf, the motifs of TALE, TCP and MYB were more enriched. TALE was particularly notable because it was the most frequent motif being co-present with DELLA. The putative DELLA-TALE bearing gene families included genes affecting germination (Nuclear Factors Y, IDDs), growth (GATAs), internode elongation (GRFs), shoot and leaf expansion (YABBY, BEL), as well as floral induction phase-change (SPLs, VvTFL1) pathways. The discoveries of these diverse transcription regulators and co-presences of their motifs with DELLA in the DEGs provide a potential explanation for how GA-DELLA signaling could impact so many different aspects of plant growth and development.

DELLA and root growth in grapevine
One interesting phenotypic change of the transgenic vines was an enlarged root which had poor rooting abilities (data not shown). The root tissue had the largest number of DEGs, compared with the shoot and leaf tissue in the transgenic lines. DELLA motif was statistically over-represented among the induced DEGs. Among these induced DELLA-motif-bearing DEGs were SHR, SCR and SCL3, which are GRAS genes involved in root maintenance and patterning [67]. In the transcriptome profiling of a transgenic GA insensitive poplar, the GRAS SCR was also among the genes found induced [68]. However, we observed numerous cyclin genes, D-type cyclins in particular, in a suppressed state in the leaf and root of our transgenic vines, while they were induced in transgenic poplar [68]. The D-type cyclins are mutual activators of SHR-SCR complex [69] and cytokinin [70] regulatory networks for formative cell divisions in the apical meristem, and formation of lateral organs [71]. Suppression of these cyclins in the present study may compromise their interaction with the SHR-SCR complex and cytokinin, which may explain the observed poor rooting ability in our study, in a contrast to the vigorous roots observed in transgenic poplar.
We deduce that cytokinin was diminished in our transgenic lines based on the relatively low transcript levels of its biosynthesis (i.e. IPT2 and IPT3), signal reception (e.g. AHK4), and response regulators (type B RR) components. The enlarged roots were likely the manifestation of a very low cytokinin level in the transition zone region, where cytokinin is usually kept relatively higher than auxin to favor cell differentiation and establishment of the size homeostasis of root meristems [72]. In contrast, a root tip usually maintains a high level of auxin that favors cell division needed for maintenance of the stem cell niche, which seemed to be the case in our transgenic vines as indicated by the induced state of auxin biosynthesis (i.e. YUCCA), signal receptor (i.e. TIR1), transporter (i.e. PINs) and response factors (i.e. ARFs). Critically, we observed that SHY2, the key negative regulator governing auxin:cytokinin crosstalk and balance [73], was suppressed and negatively correlated to Vvgai1 in the shoot and root. We were unable to find and verify if SHY2 was likewise suppressed in the microarray-based study of the transgenic GA insensitive poplar [68]. However, our transgenic lines seemed chronically deficient in cytokinin and high in auxin from the expression profiles discussed above. The enlarged roots in our transgenic vines were likely due to low levels of cytokinin and thus reduced levels of cell differentiation and proper control of root meristem size. In parallel, the disrupted balance between auxin:cytokinin might result in more accumulation of auxin which can lead to poor rooting ability as previously observed in Arabidopsis [74].

DELLA and grapevine flowering-induction
The induced responses of VvCO and VvCOL1 to Vvgail discovered in this study was in a sharp contrast to the findings that AtCO was suppressed by DELLA in Arabidopsis [47]. Recent study showed that VvCO and VvCOL1 positively regulated lateral bud induction and dormancy of grapevine; but in contrast to the CO and COL1 in Arabidopsis, the diurnal oscillation of VvCO and VvCOL1 was less affected by continuous light and had a higher amplitude that instead peaked at dawn [75]. Since it was shown that lower temperature favors tendril over inflorescence in grapevine [76] and GA signaling exhibited diurnal periodicity in Arabidopsis [77], GA and the photoperiod pathway genes in grapevine might have been adaptively modified to perceive warmer temperatures and light intensity instead, and less of daylength as dominant stimuli for flower induction [76].
Similar to GA-insensitive Arabidopsis [45], expression of several VvSPL genes were suppressed in the shoot and leaf of our transgenic lines. The GA-insensitive Arabidopsis was late-flowering because DELLA sequestered the SPL regulators and essentially prolonged the juvenile phase. However, the distinction between juvenile and adult phases is less defined in grapevine and obscured by the simultaneous presence of vegetative and reproductive meristems in the same shoot [78]. Recent transcriptome survey of tendril and inflorescence development revealed that VvSPLs were particularly expressed in the early stages of both tendril and inflorescence development [79], suggesting that the equivalent of phase-transition in grapevine takes place when the anlagen differentiates towards tendril or inflorescence formation. At least one of the VvSPLs, VvSBP12, was reported to carry the regulation motif of miRNA156/172 [44], which raises the possibility of the role of miRNAs in grapevine shoot architecture and anlagen developmental processes. VvSBP12's homologues in annual model plants mediate numerous phase-related developmental processes, including initiation of the first true leaves in Arabidopsis (AtSPL13 [80]), grain enlargement (OsSPL16 [81]), and ear glume development in maize (TGA1 [82]). The suppression of VvSPL genes in this study might explain the observed progressive severity of tendril-less nodes correlated with high levels of Vvgai1 expression. We noted that HEN1, a dicer-like gene involved in miRNA metabolism [46] was induced in the shoot and leaf. HEN1, however, was suppressed in the root, while another dicer homologue, CARPEL FAC-TORY (CAF) was not correlated to Vvgai1, suggesting the existence of complex interactions among DELLA, SPLs and miRNA.
We found that the gene VvSVP, a MADS box gene also involved in phase-transition dependent flowering in Arabidopsis [45,83], was significantly suppressed in leaf and shoot tissues. The suppressed expression of VvSVP and VvSBP12 in the transgenic plants was in line with the observed delayed tendril emergence and a progressive disruption in shoot phyllotaxy. This is also supported by the fact that VvSVP was considered as a positive regulator of floral transition in grapevine [13,79,84], in contrast to being a negative regulator in Arabidopsis [83].
Interestingly, we also discovered an antagonistic flowering gene, VvTFL1, induced in the root, suppressed in the shoot and not expressed in the leaf. Since VvTFL1 is cell mobile [85], the induced VvTFL1 in the root could migrate and function as an antagonist of the meristem identity gene LFY and floral integrator VvFT in grapevine shoot. This finding may suggest that roots are involved in the over-all regulation of flowering in grapevine via GA-DELLA signaling cascade.

Conclusion
Grapevine as a perennial species shared a similar DELLA-centered feedback mechanism as annual species for maintaining GA homeostasis and controlling plant development and growth through intricate interactions of DELLAs with numerous and specialized transcription factors. However, some of these interaction outcomes could be very different between perennial and annual species, as illustrated in this study that the expression behaviors of certain flowering induction and development pathway genes in the grapevine were in sharp contrast with that of annual species. It was interesting to observe that homeodomain transcription factors seemed to play important roles in coordinating vegetative and reproductive transition in grapevine in which both vegetative and reproductive meristems simultaneously co-exist in the same shoot.

Expression cassette and transformation
A binary vector construct, named as pVv::VvGAI L38H , was used in this study, containing a modified 5.1 kb-long fragment of the grape VvGAI genomic sequence cloned from V. vinifera 'Pinot Noir'. The modified VvGAI genomic sequence has a 2171-bp sequence upstream of the ATG start codon, a 1773-bp coding sequence with an introduced point mutation, which results in replacing a leucine amino acid (L) by a histidine (H) at the amino acid position 38 in the DELLA domain of the encoded GAI protein, and a 1157-bp sequence downstream of the stop codon. Other detail information about the construct was previously reported [14].
The Ralph M. Parsons Plant Transformation Facility at the University of California, Davis (https://ptf.ucdavis .edu/services) generated the transgenic grapevines through Agrobacterium-based transformation in the V. vinifera 'Thompson Seedless' background for this study. Ten transgenic Vvgai1 vines were generated and five of the representative ones were used in this study after preliminary evaluation (Table 1). A non-transgenic V. vinifera 'Thompson Seedless' was provided as a control.

Propagation of transgenic grapevines and phenotyping
Freshly received transgenic grapevines were potted in 4-in. pots using Promix Biofungicide media (Premiere Horticulture, Canada). Green stem cuttings from individual transgenic lines were used to generate multiple clones for subsequent experiments. Vines were maintained in 1-gal pots in a greenhouse by applying routine vine management practices.
Phenotypic traits, such as plant height, internode length and tendril distribution pattern were observed from five biological replicates of the transgenic lines and NT. In the case of leaf weight, the measurements were from 5 to 10 pooled leaves in each biological replicate. Similarly, total chlorophyll content was assessed from pooled leaf discs of five expanded leaves of each vine following an acetone extraction method [86]. Root characteristics were assessed from depotted pots after 5-month of growth in the greenhouse.

RNA-Seq library preparation and sequence reads processing
All tissues for RNA-Seq profiling were collected from 5-month-old vine clones of individual transgenic lines. Duplicated samples from at least two clonal vines were collected for each transgenic line and the non-transgenic control. Three types of tissues were taken. The shoot tissue sample contained 1-3 cm shoot tips of main branches, comprised of the apical region and two young, unfolding leaves. The leaf tissue sample was pooled 5th and 6th expanded leaves from the top and the root was the pool of several 3 cm-clippings of roots from depotted plants. All samples were taken fresh, flash frozen and stored at − 80°C until further processing.
RNA extraction and RNA-Seq library preparation were performed as previously described [87]. RNA-Seq libraries were multiplexed for paired-end 2 × 100-bp (root and shoot tissue samples) or single-end 100-bp (leaf tissue samples) sequencing using Illumina HiSeq 2000 at the Cornell University Biotechnology Resource Center, Ithaca, NY.

Digital RT-PCR validation
Thirty genes were selected for digital RT-PCR verification of the RNA-Seq gene expression. cDNAs were synthesized from the same mRNAs processed for RNA-Seq libraries. The Qiagen multiplex PCR plus kit was used for the 1st round of amplification of 150-200 bp fragments in individual genes. For each tissue sample, four multiplex reactions were performed to cover 40 different regions in 30 genes. For the 2nd round PCR amplification, the Illumina Nextera dual-indexing primer system was used to barcode each individual sample and the PCR reaction was purified by AmpureXP and eluted in 10 μl water. Equal amount of DNA from individual libraries was pooled together for Illumina HighSeq sequencing (2 × 150 paired end).
Three sequence fragments each with an allele-specific SNP were amplified to confidently distinguish the transgenic GAI allele (Vvgai1, [14]) from non-transgenic one (VvGAI1, GSVIVG01011710001 [19]) in the 'Thompson Seedless' background. The ratios of Vvgai1/VvGAI1 in the three regions were then averaged and the mean ratio was used to reflect the relative qRT-PCR expression levels of the transgene Vvgai1 in a tissue. The mean ratios derived from the qRT-PCR data were further used for a correlation analysis of the RNA-Seq expression levels of individual endogenous genes with that of Vvgai1 across different transgenic vines.

Differential expression analysis
For each tissue, we constructed and analyzed 12 RNA-Seq library samples: five transgenic lines and one non-transgenic Thompson Seedless, each with two biological replicates. The shoot library samples were single-end sequenced while the other two tissues were paired-end sequenced.
The abundance of each read was determined by HTseq [88] and differential expression analysis was done using EdgeR [89], following the analysis workflow described in our previous work [28]. To focus on the genes which were most likely affected by Vvgai1, we included in the further analysis only those differentially expressed genes (1.5x-fold changes at FDR ≤ 0.01) identified in the strongest transgenic line G02 and also in at least one other transgenic line when compared with the non-transgenic control (NT). We used count per million mapped reads (CPM) as the normalized expression unit in assessing the expression level of a gene.

Functional analysis
Gene Ontology (GO) analyses were done for the DEGs correlated to Vvgai1 (r > = 0.50) in individual tissues using Plant MetGenMAP [90]. To simplify pattern interpretation, the resulting lists of enriched GO terms (FDR ≤ 0.10) were clustered into similar functions using REVIGO [91]. To assess the relative responses of endogenous genes to the transgenic Vvgai1 expression, pairwise correlation was computed between the expression variation of each DEG among the transgenic lines with the relative expression level of Vvgai1 (ratio of Vvgail1/VvGAI1) in each tissue. The average correlation indices of the genes in a particular enriched biological GO term was then taken as the relative gauge of association to Vvgai1. Transcription factor enrichment and putative pairwise regulator-target interactions were facilitated by using PlantTFDB 4.0 and PlantRegMap [50]. Gene annotation, cross-referencing, and sequence homology analyses were conducted using TAIR [92], Phytozome 12 [93], PlantTFDB 4.0 and PlantMetGenMap [90]. Whenever possible, we used the gene names reported in relevant grapevine research. Otherwise we followed the nomenclature indicated in the Genoscope's 12X V. vinifera genome database (http://www.genoscope.cns.fr/externe/GenomeBrowser/Vitis/) in conjunction with the TAIR database (https://www.arabidopsis.org/). For providing further clarity and avoiding potential confusion, whenever appropriate we added a prefix "Vv" to a gene name referred in the main text as presented in Additional file 1: Table S1.

Additional files
Additional file 1: Table S1. A list of DEGs and their expression correlation to Vvgai1 in the shoot, leaf and root tissue. (XLSX 530 kb) Additional file 2: Table S2. A summary of the uniquely induced and suppressed genes (XLSX 11 kb) Additional file 3: Table S3. Pearson correlation (r) of standardized RNA-Seq and digital qRT-PCR expressions for a selection of 30 genes in the transgenic lines. (XLSX 11 kb) Additional file 4: Table S4. Transcription factor motifs significantly over-represented among the DEGs correlated to Vvgai1 (r ≥ 0.70) in the shoot, leaf and root tissues. (XLSX 13 kb) Additional file 5: Table S5