Genome-wide identification and expression analysis of GRAS family transcription factors in tea plant (Camellia sinensis)

GRAS proteins are important transcription factors that play multifarious roles in regulating the growth and development as well as stress responses of plants. Tea plant is an economically important leaf -type beverage crop. Information concerning GRAS family transcription factors in tea plant is insufficient. In this study, 52 CsGRAS genes encoding GRAS proteins were identified from tea plant genome database. Phylogenetic analysis of the identified GRAS proteins from tea plant, Arabidopsis, and rice divided these proteins into at least 13 subgroups. Conserved motif analysis revealed that the gene structure and motif compositions of the proteins were considerably conserved among the same subgroup. Functional divergence analysis indicated that the shifted evolutionary rate might act as a major evolutionary force driving subfamily-specific functional diversification. Transcriptome analysis showed that the transcriptional levels of CsGRAS genes under non-stress conditions varied among different tea plant cultivars. qRT-PCR analysis revealed tissue and development stage-specific expression patterns of CsGRAS genes in tea plant. The expression patterns of CsGRAS genes in response to abiotic stresses and gibberellin treatment suggested the possible multiple functions of these genes. This study provides insights into the potential functions of GRAS genes.

Tea is a non-alcoholic beverage accepted by an increasing number of people because of its beneficial functional components 1 . The active substances, such as tea polyphenols (especially catechins), caffeine, theanine, and polysaccharides, have been reported being a stimulant and an antioxidant, and reduce the frequency of cancers, and neurodegenerative diseases [2][3][4] . Tea plant (Camellia sinensis (L.) O. Kuntze) is an important perennial evergreen leaf-used economic crop that originated from subtropical regions [5][6][7] . Tea plant have a life span of over 100 years 3 . Thriving in the wild, tea plant is seriously affected in terms of growth, yield, and quality by various adverse environmental stresses, including salt, drought, cold, and heat [8][9][10][11] .
With the advent of next-generation sequencing technology, a great number of genomics, transcriptomics, proteomics and metabolomics have been reported, which provide great opportunity for the development of metabolic engineering [12][13][14][15][16][17] . Recently, the high quality and available of genome sequences of tea plant have been published, which greatly promoted the understanding the biological processes and molecular/cellular mechanisms of tea plant 18 .
Transcription factors (TFs) are critical components in physiological processes and regulatory networks of higher plants; TFs serve as a switch of gene regulation, promote or inhibit specific functional gene expression, and maintain the growth and development as well as stress responses of plants [19][20][21] . A large proportion of the plant genome is dedicated to TFs. For example, 7.5% of the Arabidopsis genome encodes for putative TFs 22 . A recent tea plant transcriptome analysis has indicated that at least 2543 unigenes belong to 58 TF families 23 .
The name GRAS family TFs is derived from the first three isolated members, namely, GAI, RGA, and SCR 24,25 . It was initially taken as a plant-specific TF family. However, a recent study has identified the GRAS gene family in bacteria and classified it under the Rossmann fold methyltransferase superfamily 26 . Typically, GRAS TFs share a widely conserved C-terminal GRAS domain and a variable N-terminal region, and total amino acids range from 400 to 770 24,25,27 . The highly conserved C-terminus GRAS domain is constituted by five critical conserved motifs in the following order: LHR I (leucine heptad repeat I), VHIID, LHR II (leucine heptad repeat II), PFYRE,

Results
Identification of CsGRAS TFs in tea plant. A total of 52 CsGRAS proteins were confirmed and used for further analyses (Table S1). Hidden Markov models (HMMs) of these CsGRAS proteins were determined using the Pfam database. All of these CsGRAS proteins contain a conserved GRAS domain (PF03514.13). However, five of the CsGRAS members (CsGRAS11, -12, -16, -24 and -48) possess a DELLA domain (PF12041.7) and thus were considered DELLA proteins (Table S2).
The physical and chemical characteristics of tea plant GRAS proteins were analyzed using ExPasy (Table S3). The length and molecular weight (Mw/Da) of the deduced GRAS proteins ranged from 416 to 1340 amino acids and 46291.14 Da to 146535.66 Da, respectively. The theoretical isoelectric point (pI) of most CsGRASs is slightly acidic (4.88-6.84), and two GRAS proteins are alkalescent. In most CsGRAS proteins, the number of positive amino acids (Arg + Lys) is less than that of negative amino acids (Asp + Glu). The grand average of hydropathy (GRAVY) of all CsGRAS proteins varied from −0.618 to −0.074, indicating that these proteins are hydrophilic. The instability index of most CsGRAS proteins exceeded 40 (42.42-62.37), indicating that these proteins are unstable. No CsGRAS proteins were predicted to contain transmembrane helices by using the TMHHM server. Subcellular location analysis showed that most CsGRAS proteins were located in the nucleus (Table S4).

Evolution of GRAS TFs among different plant species.
A comparison analysis of GRAS TFs in different species was performed to understand the evolution of GRAS TFs among different plant species (Fig. 1). GRAS family members were only identified in higher plants; no members were found in lower plants and fungi (Volvox carteri, Chlamydomonas reinhardtii, Ostreococcus lucimarinus and Saccharomyces cerevisiae). Evolutionary analysis indicated that the GRAS proteins may have expanded since the divergence from lower plants to higher plants, and these GRAS factors are closely related to higher plant evolution. The number of GRAS TFs in tea plant (52) was similar with closely related species, Coffea canephora (50), Solanum lycopersicum (53). The number and genome size of the GRAS proteins had no evident regularity between monocots and eudicots as well as woody and herbal plants. Some species, including Populus trichocarpa, Glycine max, and Zea mays, exhibited a relatively large number of GRAS family members (i.e., over 100). The density of the GRAS proteins to genome size was rather low in several plant species, which are particularly evident in two gymnosperms, namely, Picea abies (0.0015 number/Mb) and Pinus taeda (0.0004 number/Mb). This occurrence is mainly due to the large genome size of these two species.
Phylogenetic analysis of GRAS proteins in tea plant, Arabidopsis, and rice. A phylogenetic analysis of GRAS proteins from tea plant, Arabidopsis, and rice was conducted to clarify the evolutionary relationships among tea plant GRAS proteins. We acquired 34 GRAS proteins in Arabidopsis and 60 in rice, which had been proven in a previous study 32 . Meanwhile, some members possess atypical GRAS domains, and some of them are severely truncated. For example, the domain sequences of some GRAS proteins are less than 100 amino acids (e.g., ΨOsGRAS2, ΨOsGRAS3, ΨOsGRAS5, ΨOsGRAS7), whereas the approximate number of typical GRAS domain is 350 amino acids. The short or atypical sequences may be putative pseudogenes and derived from ancient pseudogenization events 32 . Considering that incorporating these fragments results in low reliability in multiple sequence alignment and phylogenetic analysis, we excluded these fragments for later analyses. Finally, 33 GRAS members in Arabidopsis and 50 in rice were identified (Table S5). Similar to those in tea plant, these proteins possess a typical GRAS domain.
Phylogenetic trees were constructed with alignments of 135 full-length GRAS proteins, including 52 from tea plant, 33 from Arabidopsis, and 50 from rice, and the neighbor-joining (NJ), minimum evolution (ME), and Scientific REPORts | (2018) 8:3949 | DOI:10.1038/s41598-018-22275-z maximum likelihood (ML) algorithms were employed ( Fig. 2 and Figs S1 and S2). The NJ and ME trees exhibited similar topologies. Compare with the NJ and ME trees, the ML tree exhibited only minor modifications. In this study, the NJ phylogenetic tree was used for further analysis. On the basis of the phylogenetic analysis and existing research, all GRAS proteins were grouped into 13 subfamilies, namely, HAM, DELLA, AtSCL3, DLT, AtSCR, AtLAS, AtSCL4/7, AtSHR, AtPAT1, Os4, Os19, Os43, and LISCL. The CsGRAS proteins were unevenly distributed. For instance, most of the CsGRAS proteins were distributed in the AtPAT1 (9 members), and LISCL (10 members) subfamilies, whereas no members were found in the AtSCL4/7 subfamily.
Estimation of functional divergence. To investigate the level of functional divergence between gene clusters of the GRAS subfamily, type I and type II functional divergence were estimated using DIVERGE (v2.0) software on the basis of the ML algorithm [47][48][49] , which provides a sequence-based method to identify cluster-specific amino acid sites that are likely responsible for functional changes after gene duplication. Type I functional divergence between two duplicate clusters indicates that some amino acid residues are highly conserved in one duplicate cluster but highly variable in another cluster, suggesting that these sites may have experienced shifted functional constraints 47 . Type II functional divergence between two duplicate clusters indicates that some amino acid residues are highly conserved in both duplicate clusters but exhibits distinct biochemical properties, suggesting that these sites may be responsible for functional specification in different subfamilies 50 . On the basis of the GRAS protein NJ tree, GRAS subfamilies HAM, DELLA, AtSCL3, AtSCR, AtSHR, AtLAS, DLT, AtPAT1, Os4, and LISCL were used in the present study as input for the DIVERGE2 analysis, whereas Os19, Os43, and AtSCL4/7 were excluded because of their small size (Table 1).
Functional distance analysis was conducted and a star-like tree in terms of the functional branch length (b F ) of each subfamily was generated (Fig. 3) to examine the pattern of type I functional divergence for each pair of subfamilies. The functional branch lengths of the subfamilies were measured and ranked as follows . LISCL exhibited the longest functional branch, indicating that the evolutionary conservation may have been changed at many sites for this subfamily and that the derived functional state may be far away from the ancestral state. By contrast, the bF for AtSCR was virtually zero, suggesting that the evolutionary rate of each site has not changed substantially since duplication and that this subfamily is likely to have inherited the ancestral function.
The coefficients of type I functional divergence (θ I ) between any two subfamilies were significantly greater than 0 (p < 0.05, Table 1), suggesting that site-specific rate shift after gene duplication is a common phenomenon in the evolution of the GRAS family. The coefficient of type II functional divergence (θ II ) was smaller than 0 between most pairs but showed relatively high standard errors. To further explore the pattern of functional divergence between GRAS subgroups, posterior analysis was conducted to identify the critical amino acid sites (CAASs) that may be highly relevant to functional divergence. To reduce possible false positives, CAASs between gene clusters were identified using posterior probability (Q k ) > 0.95 as the cut-off. Results showed that the numbers of CAASs were considerably varied among each comparison (Table 1). Some CAASs were supposed to be responsible for type I functional divergence, of which 16,24,24,11,26, and 44 CAASs were identified for pairs AtPAT1/AtSHR, AtPAT1/HAM, AtSHR/LISCL, HAM/AtSCL3, HAM/LISCL, and DELLA/LISCL, respectively. In contrast to type I functional divergence, 52, 94, 54, and 56 type II-related CAASs were identified for pairs AtPAT1 /AtLAS, AtSHR/HAM, HAM/DELLA, and DELLA/LISCL, respectively. No CAASs were identified for type I functional divergence between AtSHR/HAM and DLT/Os4, whereas 94 and 44 respectively predicted sites existed for type II functional divergence. Thus, the rapid change in amino acid physiochemical properties can be mainly attributed to the evolutionary force driving the functional divergence between this pair and secondarily to the shift in evolution rate. This research indicated that site-specific shifts in evolutionary rate and changes in amino acid property are inconsistent acts on the GRAS subfamily over the course of evolution.
Analysis of the sequence alignment and conserved motifs of the GRAS proteins. Multiple alignments indicated that all GRAS proteins contain a widely and conserved C-terminal GRAS domain, which could be divided into five parts (LHRI, VHIID, LHRII, PFYRE, and SAW) (Fig. S3), and a highly variable N-terminal domain (data not shown) as described in previous studies 24,25,27 . The five conserved motifs could be further subdivided into small distinct units. For example, both LHRI and LHRII are composed of two repeat units (A and B); VHIID can be divided into three units (A, B, and C); PFYRE contains three distinct parts (P, FY, and RE); SAW is composed of four units (RVER, W-G, L-W, and SAW). Overall, the results are slightly similar to previous statements.
To further examine the sequence features of GRAS TFs from tea plant, a comparative analysis of the conserved motifs was performed between tea plant and Arabidopsis GRAS proteins. Thirty motifs were predicted to reveal the details of the GRAS protein structure using the MEME program (Table 2 and Fig. S4). In general, GRAS proteins clustered in the same subgroups share similar motif compositions (Fig. 4), indicating functional similarities among members of the same subgroup. Nearly all of the GRAS members contained motifs 1, 3, 5, 7, 9, and 10, suggesting the important roles of these motifs in the GRAS gene family. Motifs 17,20,21,23,25, and 27 were only found in subgroup LISCL; 26 was only found in subgroup HAM; 2 was completely lost in subgroups AtSCR, Os19, and HAM; 15, 18, 22 were only found in subgroup DELLA; 24 and 28 were only found in subgroup AtPAT1; 30 was only found in subgroup AtSCR; 6 was found and 11 was completely lost in subgroups AtPAT1 and LISCL. The differences in motif distribution among the subgroups of GRAS genes revealed that the functions of these genes may have diverged during evolution.

Functional identification and interaction network of the GRAS proteins in tea plant.
Annotations of the CsGRAS factors were retrieved from tea plant transcriptome databases and were predicted from integrated function annotation databases (COG, GO, KEGG, Swissprot, TrEMBL, Nr, and Nt) (Table S6). The predicted CsGRAS genes possibly participate in diverse development processes, such as "floral organ morphogenesis", "regulation of seed dormancy process", "post-embryonic development", "seed germination", "post-embryonic morphogenesis", "reproductive structure development", "regulation of post-embryonic development", "organ development", "petal formation", "sepal formation", "radial pattern formation" and "leaf development". CsGRAS genes may be involved in various abiotic and biotic stresses, such as, "response to xenobiotic stimulus", "hyperosmotic response", "response to salt stress", "response to temperature stimulus", "response to osmotic stress", "response to red or far red light" and "defense response to bacterium". As transcriptional regulators, CsGRAS genes may be associated with the biosynthetic and signaling pathways of multi-hormones, such as GA, jasmonic acid, abscisic acid, salicylic acid, auxin, and ethylene. To understand further the interactions of the GRAS genes in tea plant, an interaction network was built using STRING software on the basis of the orthologs in Arabidopsis (Fig. 5). The homologous proteins with the highest bit score were considered STRING proteins. SCL13 acts as a positive regulator of continuous red light signals downstream of phytochrome B (PHYB) and is required for the regulation of hypocotyl elongation during de-etiolation. SCL13 may also be required to modulate phytochrome A (PHYA) signal transduction in a phyB-independent manner. PAT1 (CsGRAS1, 15, 18, 29, and 43) may be involved in PHYA signal transduction. SHR (CsGRAS3, 36, and 40) is involved in radial pattern formation in roots and is required for normal shoot gravitropism, which could directly control the transcription of SCR(CsGRAS23, 30, and 33) and of MGP, RLK, TRI, NUC, and SCL3 (CsGRAS2, 10, 17, 37, and 46) when associated with SCR. GID1A, GID1B, and GID1C interact with specific DELLA proteins that are required for GA signaling that controls root growth, seed germination, and stem elongation. Compared with other DELLA proteins, RGA1 (CsGRAS11, 12 and 48) is the most sensitive to GA application, whereas GAI (CsGRAS16 and CsGRAS24) is less sensitive to GA 51 .
Expression profiles of CsGRAS genes in four tea plant cultivars. Using RNA sequence transcriptome data, we investigated the transcript levels of the putative GRAS genes among different tea plant cultivars under normal condition. A heat map was exhibited to assess the transcription profile of the CsGRAS genes among four tea plant cultivars on the basis of the reads per kilobase per million (RPKM) values ( Fig. 6 and Table S7). All of the 20 CsGRAS genes were widely expressed across a variety of cultivars, apart from CsGRAS20, which was not detected in Tea_T3. This gene may have lost or weakened its functions in Tea_T3 during evolution that its expression became either too low for detection or had spatial and temporal patterns. Overall, most CsGRAS genes (except for CsGRAS19 and CsGRAS20) exhibited a relatively high level (RPKM > 1) in all of the four tea plant cultivars. Some genes, such as CsGRAS7, 11, and 12, were stable and highly expressed. Some obviously exhibited cultivar specificity. For example, CsGRAS3, 16,19, and 20 displayed higher expressions in Tea_T4 than in the three other cultivars; on the contrary, CsGRAS2 and CsGRAS8 showed a lower expression in Tea_T4 than in the others. In addition, CsGRAS1, 4,6,7,9,10,13,14,15, and 18, displayed higher expression levels in Tea_T1 and Tea_T2 than in Tea_T3 and Tea_T4. Among the three genes encoding DELLA proteins, CsGRAS11 and CsGRAS12 were stable and highly expressed in all tea plant cultivars, whereas CsGRAS16 was variable. The diversity of expression indicates that the CsGRAS genes may show functional differentiation among different tea plant cultivars.

Expression profiles of CsGRAS genes in different tissues of tea plant. The expression patterns of
CsGRAS genes in different tissues or different development stages (root, stem, bud, first leaf, second leaf, and third leaf) of two tea plant cultivars ('Huangjinya' and 'Yingshuang') were analyzed using qRT-PCR to explore the potential functions of GRAS genes during the vegetative development of tea plant. Eighteen CsGRAS genes were investigated. The 18 genes were representative and can explain the expression profiles of genes from different subgroups. All of the desigened primers were suitable for qRT-PCR based on the corrected sizes of PCR amplified band and the validated by PCR paired-end sequencing (Figs S5 and S6).
Results indicated that these CsGRAS genes were detected in all of the tissues of the two tea plant cultivars and exhibited variable transcript levels (Fig. 7). CsGRAS1, 3,5,7,9,11,15,16,17, and 18 were expressed at relatively high levels, whereas the other genes were expressed at low levels. Most of the CsGRAS genes presented similar tissue-specific expression patterns between two tea plant cultivars. For example, in two tea plant cultivars, CsGRAS2 and CsGRAS7 showed relatively high expression levels in the roots; CsGRAS16 and CsGRAS19 were abundantly expressed in the buds; CsGRAS3 was highly expressed in the third leaf; CsGRAS8 was highly expressed in the stems and leaves and weakly expressed in the roots and buds. The transcript levels of CsGRAS1 and CsGRAS13 gradually increased in the first, second, and third leaves. Some differences also existed. For

Expression analysis of CsGRAS genes under GA treatment in two tea plant cultivars.
Considering that many GRAS proteins are involved in regulating GA response, we examined the response of CsGRAS genes to GA treatment in 'Huangjinya' and 'Yingshuang' through qRT-PCR. Under GA treatment (Fig. 8), three CsGRAS genes (CsGRAS1, 3, and 17) were obviously upregulated in both tea plant cultivars, whereas CsGRAS2, 6, 7, and 15 were obviously downregulated. Three CsGRAS genes (CsGRAS11, 13, and 14) were mainly upregulated in 'Huangjinya' but downregulated in 'Yingshuang' . CsGRAS8 and CsGRAS10 were mainly upregulated in 'Yingshuang' but showed a relatively stable expression level in 'Huangjinya' .

Expression analysis of CsGRAS genes under different abiotic treatments in two tea plant cultivars.
The responses of CsGRAS genes to salt (200 mM NaCl), drought (20% PEG 6000), cold (4 °C), and heat (38 °C) treatments in 'Huangjinya' and 'Yingshuang' were examined through qRT-PCR to elucidate further the role of GRAS genes in tea plant during diverse environmental stresses. Except for CsGRAS3, all of the analyzed genes responded to at least one abiotic stress treatment. Although most CsGRAS genes exhibited a similar tendency in the two tea plant cultivars, some notable exceptions were also observed.

Discussion
In recent years, a large number of TFs have been identified and analyzed in various species through bioinformatics technology. Compared with that of other TFs families, the research progress of the GRAS family is relatively slower. To date, insufficient information is available about GRAS TFs in tea plant. In the present study, 52 CsGRAS genes were identified from tea plant genome database.
Gene loss and duplication events are the dominant driving forces in species evolution, supplying the raw materials for the generation of novel gene functions and promoting the formation of gene families. No GRAS family members were found in lower plants. However, a recent study has identified the GRAS gene family in bacteria and classified it under the Rossmann fold methyltransferase superfamily 26 . Although methyltransferase activity is possibly lacking in most GRAS proteins in plants, the GRAS proteins in plants still bind a similar substrate to those in bacteria 32 . Exon-intron organization analysis on several species, such as Prunus mume, Arabidopsis, rice, populous, and grapevine, reported that most GRAS genes either lack introns or possess only a single intron 27,30,32,35 . Thus, horizontal gene transfer from ancient prokaryote genomes of bacteria may explain the origin of plant GRAS genes.
Gene duplication and differentiation are the major pathways of origin for new genes and for the differentiation of gene function 52 . Members clustered in similar subgroups commonly exhibit functional similarities. Therefore, phylogenetic analysis could facilitate functional genomics. In the present study, we identified 135 GRAS proteins with full-length domain sequences from tea plant, Arabidopsis, and rice, and classified them into 13 subfamilies on the basis of sequence structures and phylogenetic relationships. All of the five predicted tea plant DELLA (CsGRAS11, -12, -16, -24, and 48) proteins together with other DELLA proteins (AtGAI, AtRGA, AtRGL1, AtRGL2, AtRGL3, and OsSLR1) were grouped into the DELLA subfamily. Some rice GRAS members (OsGRAS1 and OsGRAS27) that were similar to the other DELLA proteins but lacked the DELLA domain were also grouped into the DELLA subfamily. Previously, Liu et al. grouped OsGRAS1 and OsGRAS27 together with the other DELLA proteins into the DELLA subfamily 32 . The results are largely consistent with previous reports. Considerable bootstrapping value supports existed for many of the defined subfamilies in the unrooted tree, although some poor supporting values remained for several clusters. This result was expected because the GRAS protein sequences used in this study had an average of approximately 600 amino acid-lengths, a constraint imposed by a large number of substitutable residues.
Phylogenetic analysis revealed the following several pairs of closely related orthologous GRAS proteins between tea plant and Arabidopsis/rice: CsGRAS8 and AtSCL15, CsGRAS11 and AtRGA and AtGAI, CsGRAS2 and AtSCL3, CsGRAS19 and AtSCL28, CsGRAS20 and AtSCR, CsGRAS3 and AtSHR, CsGRAS1 and OsCIGR1, CsGRAS5 and AtSCL1, CsGRAS4 and AtSCL9, and CsGRAS10 and CsGRAS17 with AtSCL13. Orthologs generally retained similar functions. AtSCL3, which shows a high similarity to CsGRAS2, functions as a repressor of DELLA and controls root development through the GA pathway 53,54 . CsGRAS19, which shares a high similarity to OsGRAS29 (also known as DLT), is involved in brassinosteroid signaling 55 . CsGRAS3 and CsGRAS20 show high degrees of similarity to AtSHR and AtSCR, respectively, which are associated with root and shoot radial pattern formation [41][42][43] . CsGRAS10 and CsGRAS17 show high degrees of similarity to AtSCL13, which is a positive regulator of phytochrome B signaling 56 . The results suggest the possible functions of GRAS genes in tea plant.
The structural divergence of gene sequences has played vital roles in the evolution of gene families, which is an adaptive process for speciation and leads to the efficient use of natural resources or adaptation to unfavorable conditions. Comparisons of GRAS proteins between tea plant and Arabidopsis were conducted to characterize the GRAS sequences. Most GRAS proteins clustered into similar subgroups sharing similar motifs, suggesting their group-specific functions. Significant differences were observed among the different groups. For example, six motifs (motifs 17, 20, 21, 23, 25, and 27) were only present in subgroup LISCL; two motifs (motifs 24, and 28) were only found in subgroup AtPAT1, indicating their specific functions to other subgroup members. Among them, the DELLA (motif 15) and TVHVNP (motif 18) domains function in GA perception 29,51,57 . The distributions of conserved motifs also reflect the relations of different subgroups. For example, motif 6 was found and 11 was completely lost in subgroups AtPAT1 and LISCL. This result indicates the close evolutionary relationships between subgroups AtPAT1 and LISCL. Therefore, structural analysis also provides a clue to locate which subgroup of GRAS genes is the ancient origin 58 .
Type I and type II functional divergence between subfamilies of GRAS genes were estimated. Significant type I functional divergence was detected between any two subfamilies, suggesting that site-specific rate shift after gene duplication is the main force driving the functional divergence in the evolution of the GRAS family. Type II functional divergence between most subfamily pairs was smaller than 0, whereas the standard errors were relatively high, suggesting that the GRAS gene family may adopt type II in different degrees. In the present study, the functional divergence of GRAS genes was not closely correlated with expression profiling in different organs and samples under different abiotic stresses. It may contribute potentially to gene expression profiling accumulation.
The functional annotations of CsGRAS proteins were obtained from the well-characterized biological function of similar genes. It provided a functional context for the CsGRAS genes in tea plant with modules. CsGRAS genes were predicted to have diverse functions and involved in various physiological processes, such as responses to abiotic and biotic stresses, hormones, tissue development, and secondary metabolites. The gene expression levels could provide critical clues to assess their possible functions. On the basis of previous transcriptome sequencing data, the transcript levels of the CsGRAS genes among different tea plant cultivars under normal conditions were investigated. In consideration these cultivars were obtained from different regions of China 7 , the transcription levels were varied among the four tea plant cultivars. CsGRAS7 and two DELLA genes (CsGRAS11 and CsGRAS12) were stable and highly expressed in the four cultivars, indicating their critical roles in tea leaf development and metabolite biosynthesis.
Tissue-specific genes promote the development of particular organs or tissues, and tissue-specific expressions are helpful in defining the precise nature and function of genes. Eighteen tested CsGRAS genes widely exist in the tissues (the root, stem, bud, and 1st, 2nd and 3rd leaf) between the two tea plant cultivars ('Huangjinya' and 'Yingshuang'). Some results showed great error bars in expression levels, which may mainly because the differences among biological individuals. The expression patterns of different genes varied greatly; some genes showed similar tissue-specific expression patterns between two tea plant cultivars. The transcript levels of CsGRAS1 and CsGRAS13 gradually increased in the 1st, 2nd and 3rd leaf, indicating those two genes may function in the later development of different tissues of tea plant. It is observed that the expression levels of most of the CsGRAS genes in 1st, 2nd and 3rd leaf were higher for 'Yingshuang' compared to 'Huangjinya' . That may mainly because CsGRAS genes play greater roles in 'Yingshuang' compared to 'Huangjinya' during leaf development. The DELLA gene CsGRAS11 showed high expression levels in all tested tissues, which is consistent with DELLA genes that serve as the main signaling hub in the regulation of multifarious growth and developmental processes 59,60 . The DELLA gene CsGRAS16 and DLT gene CsGRAS19 showed similar expression patterns and relatively higher expression levels in the buds than in the other tissues, indicating that these genes play a role in bud development. CsGRAS2 and CsGRAS7 showed high transcript levels in the roots. Previous studies reported that the CsGRAS2 homolog AtSCL3 regulates root development through the GA pathway 53,54 .
GA is involved in multiple aspects of plant growth and development, such as seed development and germination, stem elongation, and flower development. The prediction of the functional interaction network indicated that the GRAS proteins act as important components of the GA signaling pathway. Expression analysis indicated that most GRAS genes were induced or inhibited in at least one tea plant cultivar. In Arabidopsis, DELLA proteins CsGRAS2 is an orthologous gene in tea plant, and it is evidently downregulated under GA treatment.
Harmful environmental conditions, such as salt, drought, cold, and heat, result in irreversible damages to the growth and development of tea plant. GRAS genes reportedly play potential regulatory roles in plant response to abiotic stresses ( Fig. 13) 33,[61][62][63][64] . However, information about this family in tea plant is currently lacking. In the present study, expression analysis indicated that most GRAS genes in tea plant could be affected by one or more various stress treatments, suggesting that CsGRAS family genes play important roles in the response to adversity stresses. Some expression trend differences exist between the cultivars 'Huangjinya' and 'Yingshuang' . Most of CsGRAS genes studied have shown higher levels of expression in 'Huangjinya' compared to 'Yingshuang' under salt treatment, indicating CsGRAS genes may have more important roles in 'Huangjinya' than that in 'Yingshuang' under salt stress. The expression pattern of CsGRAS6 under heat treatment has shown exactly opposite pattern in the two tea plant cultivars. Those different trends of expression of CsGRAS genes maybe involved in the response to abiotic stress in the two tea plant cultivars. Several highly induced genes were identified in subfamilies AtPAT1 (CsGRAS1, 10, and 17) and LISCL (CsGRAS14). Two members of the AtPAT1 subfamily genes from rice, CIGR1 and CIGR2, are involved in GA and stress response 65 . AtSCL13, an AtPAT1 subfamily gene from Arabidopsis, responds to various abiotic stresses and phytohormones 56 . The rice gene CIGR1 is orthologous to CsGRAS1, and the Arabidopsis gene AtSCL13 is orthologous to CsGRAS10 and CsGRAS17. Likewise, CsGRAS14 shares a strong sequence similarity to AtSCL14, which plays vital roles in the activation of stress-inducible promoters 66 . A previous expression analysis on Castor Beans indicated that LISCL and AtPAT1 subfamily genes respond to all of the four abiotic stresses (drought, salt, cold, and heat) 34 . These genes could be preferentially utilized for further plant functional studies. Several other genes respond to abiotic stresses. OsGRAS23, an AtSCL3 subfamily gene from rice, is involved in drought stress response through regulating the expression of stress responsive genes 67 . CsGRAS2, also an AtSCL3 subfamily gene in tea plant, is evidently upregulated under drought, salt, heat, and cold treatments, especially under heat stress. DELLA proteins participate in multiple abiotic stresses, such as low temperature and phosphate starvation 68,69 . The DELLA member CsGRAS16 is upregulated under all of the four abiotic treatments in both tea plant cultivars, whereas the expression of the CsGRAS11 genes dramatically changes under the four abiotic stress treatments in 'Huangjinya' but is relatively stable in 'Yingshuang' .
Till now, tea plant transgenic system was not efficient. It was difficult to obtained transgenic tea plants. For further verify the function of those CsGRAS genes of tea plant, we shall carry out the CsGRAS genes overexpression and gene knockout using tea plant system, Arabidopsis or other model plants by transgenic method. In the future, using new and modified Sanger and next-generation sequencing techniques, more quality tea plant genome data shall be performed. We could integrate various genomics, transcriptomics, proteomics, and metabolomics to understand the CsGRAS genes function of tea plant.

Conclusion
The GRAS TF family has been characterized in several plant species and associated with various critical development and physiological processes. However, information about this gene family in tea plant is lacking. In this work, 52 genes encoding CsGRAS family TFs were identified from tea plant. The phylogenetic relationship, conserved motif composition, and functional divergence of these CsGRAS genes were systematically analyzed and compared, which provided insights into the functional diversity of the GRAS gene family. The expression analyses of CsGRAS genes among different tissues (root, stem, bud, first leaf, second leaf, and third leaf), GA treatment, and abiotic stresses (salt, drought, heat, and cold) provide valuable information about the gene functions of this TF family in tea plant. This study could serve as a reference for future function investigations and molecular breeding of tea plant.

Materials and Methods
Retrieval of putative GRAS genes of tea plant. All candidate CsGRAS genes were derived from the tea plant genome database 18 . The acquired amino acid sequences were analyzed using the BLASTP search at NCBI (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Only the sequences with full-length GRAS domain were considered as CsGRAS proteins and used for further analyses. HMMs of the selected CsGRASs were analyzed using the Pfam database (http://pfam.sanger.ac.uk) 70 . The chemical and physical characteristics of the CsGRAS proteins were predicted by ProtParam (http://web.expasy.org/protparam/). To confirm the subcellular localization of the identified CsGRAS proteins, protein sequences were predicted using WoLF PSORT (http://wolfpsort.org). TMHHM server v2.0 (http://www.cbs.dtu.dk/services/TMHMM/) was used to predict the membrane-bound CsGRAS members. The GRAS family TF databases of other species were downloaded from the plant TFDB database (http://planttfdb.cbi.edu.cn/) 71 . Phylogenetic and evolutionary analyses of the GRAS TFs. Multiple sequence alignment of the identified GRAS protein sequences was performed using the ClustalX1.83 program with default parameters and adjusted manually 72 . The sequence alignments were also used for the subsequent phylogenetic analysis. Phylogenetic trees were generated with MEGA5.02 software by using the NJ, ME, and ML algorithms, and the reliability of the obtained trees was assessed with a bootstrap value of 1000 73,74 . Functional divergence analysis was performed using the DIVERGE (v2.0) software with the algorithm based on ML procedures 47,50 . MEME v4.11.1 (http://meme-suite.org/tools/meme) was employed to identify the conserved motifs with default parameters, except for the maximum number of motifs (30). The functional interacting networks of proteins were constructed using STRING software with the confidence parameter set at 0.15 threshold 75 .
RNA-seq data analysis. The expression profiles of CsGRAS genes in the four tea plant cultivars were assessed using transcriptome sequencing data, which was previously generated and analyzed by Wu et al. 7 . The tea plant samples used included mid-leaf 'Yunnanshilixiang' (Tea_T1) from Yunnan Province, small-leaf 'Chawansanhao' (Tea_T2) from Jiangsu Province, large-leaf 'Ruchengmaoyecha' (Tea_T3) from Hunan Province, and small-leaf ' Anjibaicha' (Tea_T4) from Zhejiang Province. These four tea plant samples exhibited different environmental adaptation modes and leaf sizes 7 . RPKM values were retrieved and normalized to estimate gene expression values. Heat maps were generated, and hierarchical clustering was conducted using HemI1.0 software (http://hemi.biocuckoo.org/faq.php). Tea plant materials and stress treatments. In China, 'Huangjinya' and 'Yingshuang' are important tea plant cultivars that share several identical traits, such as high yield and good quality. 'Huangjinya' is a light-sensitive albino cultivar that shows potential for processing high-quality green tea even in summer but exhibits weak stress resistance 76,77 . 'Yingshuang' possesses appropriate phenol ammonia content and good resistance to abiotic stresses, especially cold stress. Two-year-old vegetative propagated cuttings of tea plants, such as 'Huangjinya' and 'Yingshuang' , were grown in a chamber at the Tea Science Research Institute of Nanjing Agricultural University (Nanjing, China). The seedlings were grown in pots containing a mixture of perlite, vermiculite, and sphagnum (ratio, 1:2:3). The chamber conditions were 23 °C temperature, 14/10 h day/night period, and 70% relative humidity for subsequent stress treatments. To determine the tissue-specific expression patterns of CsGRAS genes, the roots, stems, buds, and first, second, and third leaves were collected from 'Huangjinya' and 'Yingshuang' under normal conditions, rapidly frozen in liquid nitrogen, and then stored at −80 °C for RNA extraction. Those two tea plant cultivars were subject to different treatments, such as cold, heat, high salinity, drought, and GA treatments. The tea plants were treated according to the procedure basis on the previous reports with a slight modification 77,78 . In brief, for extreme temperature treatments, the seedlings were transferred to another chamber and maintained at 4 °C and 38 °C. For drought and salt treatments, the seedlings were irrigated with PEG6000 (20%) and NaCl (200 mM) solutions, respectively. For GA treatment, the seedlings were sprayed with GA3 (1 mM). After 0, 2 (short time), and 24 h (long time) treatments, the young leaves of tea plant were collected, rapidly frozen in liquid nitrogen, and then stored at −80 °C.
RNA isolation and qRT-PCR analysis. The total RNA was isolated from plant samples using the Quick RNA Isolation Kit (Huayueyang, Beijing, China) in accordance with the manufacturer's protocol. The purified RNA (1.0 µg) was reverse-transcribed into cDNA in 20 µL of reaction volume using a PrimeScript RT reagent kit (TaKaRa, Dalian, China). The cDNA reaction mixture was diluted 20-fold using deionized water for fluorescence detection.
The specific primers for tea plant CsGRAS genes were designed from the non-conserved N-terminal region using Primer Premier 5.0 software. Primer specificity was confirmed by the melting curve with a single sharp peak    (Figs S5 and S6).The Csactin gene served as an internal control to normalize the expression levels 78,79 . All primer sequences are presented in Table 3. qRT-PCR assays were completed on the Bio-Rad CFX96 fluorescence quantitative PCR platform with the following program: 95 °C for 3 min; 40 cycles of 95 °C for 5 s for denaturation and 58 °C for 30 s for annealing and extension; and 61 cycles 65 °C for 10 s for melting curve analysis. Each reaction was conducted in a 20 µL reaction mixture containing 2 µL of diluted template cDNA, 0.4 µL of each specific primer, 10 µL of SYBR Premix Ex-Taq (TaKaRa, Dalian, China), and 7.2 µL of ddH 2 O. The experiments were independently repeated thrice, and relative expression levels were measured using the 2 −ΔΔCt method 80 . Significant differences were detected using Duncan's multiple-range test at the 0.05 level with SPSS 17.0 software.