Genome-wide identification and characterization of GRAS transcription factors in tomato (Solanum lycopersicum)

Identification of tomato GRAS genes

We retrieved the tomato GRAS genome sequences, protein sequences, and annotation information from SGN (https://solgenomics.net/) (The Tomato Genome Consortium, 2012). The Arabidopsis GRAS genes’ family sequences and annotation information was downloaded from TAIR (http://www.arabidopsis.org/) (Swarbreck et al., 2008), whereas the rice GRAS transcription factor sequences and annotation information was obtained from RGAP (http://rice.plantbiology.msu.edu/index.shtml) (Ouyang et al., 2007). The grape GRAS gene family sequences were downloaded from plantTFDB (http://planttfdb.cbi.pku.edu.cn/). The HMM model of GRAS transcription factor was downloaded from PFAM (http://pfam.xfam.org/) (Finn et al., 2010), known as PF03514. The HMM model was used as a query to identify the tomato GRAS genes containing the GRAS domain with a cut-off E-value of 1e⁻⁵ in the HMMER software. Consequently, we identified a total of 54 GRAS genes in tomato, 34 from Arabidopsis, 56 from rice, and 37 from the grapevine. Then, we conducted a quality check using the Simple Modular Architecture Research Tool SMART (http://smart.embl-heidelberg.de/) (Letunic, Doerks & Bork, 2012) to confirm the presence of GRAS domain on the candidate GRAS genes. These data were used for subsequent analysis.

Phylogenetic analysis for tomato GRAS genes

MEGA 6.0 was used to analyze the phylogenetics of genome-wide GRAS gene family based on the whole set of GRAS protein sequences from tomato (solanum lycopersicum), Arabidopsis (Arabidopsis thaliana), rice (Sativa oryza), and grapevine (Vitis vinifera). The Arabidopsis is one of the most commonly used plants in Cruciferae for studying the genetic correlations; Arabidopsis to Cruciferae, grapevine to Vitaceae, and rice to Gramineae (Tamura et al., 2011). In this building process, several shorter amino acids (Solyc01g090950.1.1, Solyc04g011630.1.1, Solyc12g049320.1.1) were excluded, the domain length of these sequences were shorter than half of the typical GRAS domain (350 amino acids) with low similarity among the tomato GRAS family (Huang et al., 2015). According to a previous study, the classification was made on the phylogenetic tree using the Evolview software (http://www.evolgenius.info/).

Structure analysis of tomato GRAS genes

Using the NCBI platform (https://www.ncbi.nlm.nih.gov/), the conserved domains of 54 GRAS genes were visualized. In order to search for the potentially conserved motifs in the complete amino acid sequence of tomato GRAS proteins, the Multiple EM for Motif Elicitation (http://meme-suite.org/) (Bailey et al., 2006) was used with default parameters, except that the number of motifs was set to 10. In order to present the characteristic of every subfamily, the 54 GRAS genes were used to build a Maximum Likelihood tree based on the JTT matrix-based model. Other parameters were same as above. The secondary structures of the tomato GRAS genes were generated using the Gene Structure Display Sever (http://gsds.cbi.pku.edu.cn/) (Guo et al., 2007).

Evolution analysis of tomato GRAS genes

The sequences were compared using the GRAS genes in tomato, Arabidopsis, and rice; the entire protein sequences were used to identify the orthologous and paralogous genes using the software OrthoMCL (http://orthomcl.org/orthomcl/) with an E-value of 1e⁻⁵ and a match cut-off value 50 for against all BLASTp alignments (Li, Stoeckert Jr & Roos, 2003; Li et al., 2003). The MCscanX software (Tang et al., 2008) was used to identify the collinear block based on the tomato genomes; if a gene had more than one transcript, only the first transcript in the annotation was used. The underlying mechanism of tandem duplication showed that the two genes were physically close to each other with the genes residing within 20 kb (Liu et al., 2014), segmental duplication resulted from the whole genome duplication accompanied by a comprehensive gene loss (Tang et al., 2008), and those large duplication events can be deduced by anchor genes in collinear blocks (Cannon et al., 2004). The chromosomal localization and homologous collinear relationship of GRAS genes were visualized using the Circos program (Krzywinski et al., 2009).

Expression pattern analysis of GRAS genes

We utilized the Illumina RNA-seq data of tomato download from SGN (https://solgenomics.net/), which were reported previously. To confirm the expression patterns of the GRAS genes, the FPKM was used to represent the expression level of each tomato gene. The data on GRAS genes’ expression was retrieved to display the consequences using the HemI program (http://hemi.biocuckoo.org/); the expression data were amplified 100-fold (Deng et al., 2014). The transcriptomic data were extracted from twelve tissues of S. pimpinellifolium (LA1589), and S. pimpinellifolium was abundant in genetic variation, which can better assess the evolution of the GRAS genes. The organism comprised of the following: A, newly developed leaves approximately 5 mm long; B, mature green leaflets; C, flower buds 10 days before anthesis or younger; D, flowers at anthesis (0DPA); E, 10 days post-anthesis fruit (10DPA); F, 20 days post-anthesis fruit (20 DPA); G, breaker stage ripening fruit (33DPA); H, another set of 10DPA fruit that was collected in a separate greenhouse for comparison. The following tissues were collected from seeds that were germinated and grown for 7 days in a Petri dish under glowing lights; I, whole root; J, hypocotyl from below the cotyledons to above the root zone; K, cotyledons; L, vegetative meristems.

In the case of difficulty, other public data were retrieved to reveal the relationship between GRAS genes and biotic stress. The transcriptome of leaves of resistant (RG-PtoR) and susceptible (RG-prf3 and RG-prf19) tomato plants treated with pstDC3000 in 4 and 6 h were sequenced.

Genome-wide identification and annotation of the GRAS genes in tomato

In order to identify the GRAS proteins in tomato, we downloaded the raw data sequence of tomato GRAS transcription factors from the SGN database. The bioinformatics approach retrieved 54 tomato GRAS genes, 34 Arabidopsis GRAS genes, 56 rice GRAS genes, and 37 grape GRAS genes. In addition, we obtained basic information in connection with the tomato GRAS proteins (Table S1). The length of GRAS proteins ranged from 125–864 aa and the tomato GRAS genes were almost distributed across the 12 chromosomes uniformly; the highest content was on chromosome 1, containing eight GRAS genes. However, at least three GRAS genes were located on the other chromosomes. The annotation information represented that the GRAS domain always assembled on the C-terminal, suggesting that the conserved terminal was supported by the stable structure of the GRAS domain; this phenomenon could be used for the subsequent comparison analysis. The smaller HMM E-value provided reliable screening results. Notably, a number of amino acid residues in some GRAS proteins were less than that in the typical GRAS domain, such as Solyc04g011630.1.1 or Solyc06g076290.1.1 genes.

Genome-wide sequence alignment and phylogenetic analysis

The topology tree contains 178 proteins: 56 from rice, 34 from Arabidopsis, 37 from grapevine, and 51 from tomato. However, the remaining three proteins from tomato were neglected because of the shorter domain sequence. Basing on the topology structure, clade support value, and previous classification from rice and Arabidopsis, as shown in Fig. 1 (Tian et al., 2004), all of these proteins were divided into ten groups, named specifically after a common feature or one of the members: PAT1, SHR, LISCL, Gv6, SCR, LAS, SCL3, DELLA, HAM, and BolA, respectively. Of these, eight were named according to a previous report (Tian et al., 2004). The remaining two groups were named by the characteristic of the members. The tree showed that the four species GRAS proteins were randomly distributed in the ten groups, and hence, the members of GRAS transcription factors in these four species are not represented equally. For example, the subfamily SCL3 contains one or two proteins from Arabidopsis, tomato, and grapevine, respectively, and seven members from rice. Strikingly, some clades do not contain the members of GRAS proteins from Arabidopsis or rice; for instance, Gv6 and BolA groups. This might be attributed to the loss of the corresponding member following the separation of the last common ancestor if not the problem of assembling or annotation of the Arabidopsis or rice genomes. We found that the Gv6 group consisted of six tomato proteins and two grapevine proteins, indicating that it is a fruit tree species-specific clade.

The genes with orthologs frequently tend to be clustered together, and the subfamily members in a major group share similar gene structure and function. Therefore, we observed that the distribution of specific genes provides information on the role of the other genes in the same clade. For example, within the PAT1 group, AT5G48150.1 of Arabidopsis was previously demonstrated to participate in the process of phytochrome signal transduction. In addition, the Solyc07g047950.1.1 was highly similar to the AT5G48150.1 protein, and hence, we inferred that Solyc07g047950.1.1 protein also played a crucial role in phytochrome signal transduction. PAT1 and LISCL subfamily comprised of more proteins than the other subfamilies as well as that reported in previous studies.

Members in the same sub-branch were marked by the same color and surrounded with an appropriate color line that was commanded by the subfamily name. The red circle corresponds to the tomato GRAS proteins, the green circle corresponds to the Arabidopsis GRAS proteins, the blue circle corresponds to the rice GRAS proteins, and the pink circle corresponds to the grapevine GRAS proteins.

Structural analysis of GRAS genes (motif analysis and exon/intron analysis)

To achieve a general overview of the conserved features of the tomato GRAS proteins, we performed a multiple sequence alignment of the 54 tomato GRAS genes using MEME software; the conserved motifs of all proteins are shown in Fig. 2B. Ten default motifs were identified by MEME analysis. We found that 80% GRAS proteins contained >7 motifs, and 24% GRAS proteins encompassed 10 motifs. Moreover, members of the same clade of phylogenetic tree shared similar motif organization with respect to either gene length or motif number (Liu et al., 2014). On the other hand, several proteins only contained one or two motifs, which might be attributed to the short duration, resulting in a shorter domain. Also, the motifs were more likely to be located in the C-terminal than the N-terminal. The motif logo is shown in Fig. 2A. According to the previous studies on GRAS motifs (Tian et al., 2004), the five putative GRAS domains were exhibited in the logo of motifs; we speculated that each LHRI, LHRII, VHIID, PFYRE, and SAW motifs could be divided into different units. This phenomenon demonstrated that maximum tomato GRAS genes are conserved, and the conserved region is localized on the C-terminal. Furthermore, we conducted a domain analysis (Fig. 2C) and exon/intron analysis (Fig. 2D) of tomato GRAS genes in order to explore the diversity and functionality. The visualization of the domains revealed that nearly all the GRAS proteins domains were primarily distributed in the C-terminal, further supporting the theory that the C-terminal is highly conserved. Nevertheless, there were some specific genes, such as the whole Solyc06g076290.1.1 gene has one GRAS domain, Solyc10g086380.1.1, Solyc11g011260.1.1, Solyc07g063940.1.1, Solyc01g090950.2.1, and Solyc09g090830.2.2 genes contained two domains. In addition, the exon/intron analysis displayed similar distribution characteristics. The Solyc06g076290.1.1 and Solyc090830.2.2 genes are comprised of multiple exons; however, the majority of GRAS proteins possessed only one exon and one intron, which might result from intron gain or loss event during evolution.

Comparison and expansion analysis of tomato GRAS genes

To evaluate the evolutionary relationships among tomato GRAS genes, we made a comparative analysis to identify the co-orthologous, orthologous, and paralogous GRAS genes among tomato, Arabidopsis, and rice using the OrthoMCL software, as shown in Fig. 3. We identified 22 co-orthologous and 46 orthologous gene pairs between tomato and Arabidopsis, 45 co-orthologous and 48 orthologous gene pairs between tomato and rice, and 19 co-orthologous and 29 orthologous gene pairs between Arabidopsis and rice. Moreover, 38 tomato GRAS genes (65%) have paralogous genes, which was higher than that in rice (63%) and Arabidopsis (19–56%). The orthologous genes commonly share a similar structure and biological function. The number of orthologous gene pairs between tomato and rice was more than that between tomato and Arabidopsis, which indicated that tomato is similar to rice.

Moreover, the duplication events were discovered in the evolution of tomato GRAS transcription factors (Cannon et al., 2004). Tandem and segmental duplication were vital for the expansion of the GRAS family. To reveal the relationship between gene duplication and amplification, the syntenic regions were analyzed by MCscanX software. As shown in Fig. 4, finally, we obtained five tandem duplication gene pairs (18.5%) (the details of the duplication events were summarized in Table S2), and the result suggested that the origination of GRAS genes applied to the tandem duplication events. We also assessed the contribution of segmental duplications; 18 genes (33.3%) with duplications were harbored in collinear blocks, indicating their robust participation in the expansion of GRAS genes. The same observations were also made in other protein families (Wang et al., 2016).

Analysis of expression pattern tomato GRAS genes in different tissues

In this study, we analyzed different expression levels of GRAS proteins in various tissues regarding the published RNA-seq data. The heat map was constructed to show the expression profile and the genes listed according to their branch (Fig. 5). We found that only the Solyc01g059960.1.1 gene was not detected in the RNA-seq data that might be due to the absence of temporal expression (Song et al., 2014). The other tomato GRAS genes were obtained in at least one tissue. As a whole, we found that the same group of GRAS genes shared a similar expression pattern; for instance, nearly all the PAT1 subfamily members showed a higher expression level than that of other groups, and the Gv6 subfamily showed a low expression level in all tissues. The GRAS genes on some branches exhibited a tissue-specific expression pattern. For example, the Solyc11g005610.1.1, Solyc03g123400.1.1, and Solyc09g018460.1.1 genes belonged to the HAM subfamily that showed a low expression in those tissues except the hypocotyls, thereby suggesting that the three genes contributed considerably to the development of the hypocotyl. The Solyc04g014830.1.1 gene only expresses in the anthesis stage, and thus, would be closely related to the flower opening. Solyc11g011260.1.1, Solyc12g005340.1.1, and Solyc11g012510.1.1 genes are always expressed at a high level in all the tissues, indicating their crucial role in the whole growth process of tomato. We also determined that some genes present a time-specific expression, such as Solyc01g100200.2.1, Solyc01g008910.2.1, and Solyc07g052960.1.1 which show high expression level in the stage of fruit ripening; the Solyc10g074680.1.1 gene was expressed in the later stage. During the mature period of leaf and flower, the expression of the GRAS genes remains stable, and only two or three genes show a difference; for example, the Solyc11g013150.1.1 gene merely did not express in newly developed leaves, the Solyc01g.008910.2.1 only expressed in flowers, and Solyc02g092570.1.1 only expressed in the newly developed leaves.

According to the public database, using the visualized tool HemI to present the relationship between pstDC3000 and GRAS genes from tomato (Fig. 6), we deduced that the number of higher expression genes in susceptible tomato (54%) were more than that in the resistant tomato (41%). The group of PAT1 continued to show the highest expression level, and the members of the PAT1 subfamily might regulate some critical physiological processes in tomato growth that are related to resistance. During the same infection time, the two varieties of susceptible plant shared the same tendency of the expression of GRAS genes. As the infection time continued, the SISCL group genes’ expression was on a downward trend in the resistance tomato, such as Solyc06g076280.1.1, Solyc01g100200.2.2, and Solyc06g076290.1.1.; interestingly, the Solyc05g053090.1.1 gene was on an increasing trend. In the susceptible tomato, the Solyc02g085600.1.1 and Solyc07g063940.1.1 genes’ expressions were increasing, whereas the expression of the other members of the HAM, LISCL, and PAT1 groups were declining. In 4 h, the expression of Solyc01g100200.2.1, Solyc07g063940.1.1, and Solyc12g005340.1.1 genes was different between the resistance tomato and the susceptible tomato. In 6 h, the Solyc12g005340.1.1 gene expression differed largely, and the results suggested that the Solyc12g005340.1.1 gene regulated the resistance of the plant disease.