Genome-wide characterization and expression analysis of GRAS gene family in pepper (Capsicum annuum L.)

Identification and annotation of pepper GRAS genes

Whole genome data for pepper cv. CM334 and cv. Zunla-1 were used for this study, and their genomic information were downloaded from http://peppergenome.snu.ac.kr/download.php and http://peppersequence.genomics.cn/, respectively (Kim et al., 2014; Qin et al., 2014). Arabidopsis GRAS genes were obtained from TAIR (https://www.Arabidopsis.org/), whereas rice GRAS genes were downloaded from RGAP (http://rice.plantbiology.msu.edu) (Tian et al., 2004). The tomato GRAS information was obtained from SGN (https://solgenomics.net/) (Niu et al., 2017). The latest Hidden Markov Model (HMM) of GRAS domain (PF03514.11) (http://pfam.sanger.ac.uk/) was used as a BLAST query to search against the entire protein datasets of cv. CM334 and cv. Zunla-1 with an E-value of 1e⁻⁵ using HMMER 3.0 (Huang et al., 2015). Meanwhile, all AtGRAS and OsGRAS proteins were used as queries to search against the two pepper databases using default parameters. The length of the hit out of the range from 350 to 820 aa was rejected. In order to validate their putative accuracy, conserved domains essential for GRAS proteins were evaluated by SMART (http://smart.embl-heidelberg.de/) and PFAM database. Finally, all outputs from two independent databases were aligned and those having similar GRAS core domain were deemed as the same gene. After these stringent criterions, sequences with the presence of GRAS domain were retained for further analysis. In our study, we refer to the variety cv. CM334 as the reference for subsequent whole genome-wide analysis.

Phylogenetic analysis of CaGRAS genes

All screened GRAS proteins from pepper, Arabidopsis, rice and tomato were used for multiple alignments by ClustalW program (Larkin et al., 2007). Gene IDs of GRAS members used in this study were listed in Table S1. Arabidopsis and rice are most common used model plants for researching genetic correlations, and tomato is another model plant of the Solanaceae family, which is closely related to pepper. Then maximum likelihood method was adopted to generate an unrooted phylogenetic tree based on alignment results. Reliability of phylogenetic tree was estimated with 1,000 bootstrapping replicates (Tamura et al., 2013). GRAS members in pepper were further categorized into different subfamilies based on well-established classifications in Arabidopsis (Tian et al., 2004).

Protein property and gene structure analysis

With the help of multiple expectation maximization for motif elicitation (MEME, http://meme-suite.org/), conserved motifs of GRAS proteins were scanned with the following parameters: (1) maximum number of motif was 12; (2) optimum motif width was set from 6 to 50 aa (Bailey et al., 2009). These identified motifs were further validated using InterProScan (http://www.ebi.ac.uk/Tools/pfa/iprscan/) (Mulder & Apweiler, 2007). The properties of GRAS proteins were calculated on ExPASy online server (http://web.expasy.org/), such as molecular weight (MW), isoelectric point (pI), instability index and GRAVY value (grand average of hydropathy) (Gasteiger et al., 2003). Based on the relationship of coding sequence and its corresponding genomic DNA sequence, the final exon/intron distribution of each CaGRAS gene was illustrated by GSDS 2.0 (gene structure display server, http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015).

Chromosomal mapping and gene duplication analysis

Physical position of each CaGRAS gene was extracted from pepper genome annotation file, and plotted onto the corresponding chromosome using Mapchart 2.3 (Voorrips, 2002). We renamed each CaGRAS gene according to its ascending chromosomal distribution. Duplicated gene pairs and patterns in pepper were analyzed by using MCScanX and BLASTP (Wang et al., 2012). Tandem duplicated genes were characterized as contiguous homologous genes located in a 100 kb single region or separated by less than five genes, while the whole blocks of genes copying from one chromosome region to another were defined as segmental duplications (Tang et al., 2008). Subsequently, non-synonymous (Ka) and synonymous substitution (Ks) between duplicated CaGRAS gene pairs were calculated by PAL2NAL (http://www.bork.embl.de/pal2nal/) (Suyama, Torrents & Bork, 2006). The microsyntenic map of precise region containing GRAS genes among pepper, tomato and Arabidopsis was created by MCScanX and plotted using Circos (Krzywinski et al., 2009; Wang et al., 2012).

Prediction of CaGRAS protein–protein interaction network

To further clarify the relationships between CaGRASs, a protein–protein interaction network was predicted using their interologs from Arabidopsis. First, specific interolog relationships between Arabidopsis AtGRASs and pepper CaGRASs were mapped from INPARANOID database (http://inparanoid.sbc.su.se/cgi-bin/gene_search.cgi) (Remm, Storm & Sonnhammer, 2001). Then, we retrieved the interaction information among AtGRASs from AraNet database (http://www.functionalnet.org/aranet/) and mapped these attributions to CaGRASs to generate corresponding interaction relationships for pepper (Guo et al., 2015b; Lee et al., 2010). Finally, these interaction networks among CaGRASs were visualized using Cytoscape version 3.4.0 (Shannon et al., 2003).

Expression analysis of CaGRAS genes in different tissues

The public transcriptome data of leaf, stem, root, pericarp and placenta at mature green, breaker, five and 10 days post-breaker, six, 16 and 25 days post-anthesis (PC-MG, PL-MG, PC-B, PL-B, PC-B5, PC-B10, PL-B5, PL-B10, PC-6DPA, PC-16DPA, PC-25DPA, PL-6DPA, PL-16DPA, PL-25DPA) for pepper cv. CM334 have been previously generated (Guo et al., 2015a; Kim et al., 2014). We retrieved the fragments per kilobase per million reads value representing the expression level of each CaGRAS gene and displayed the result using BAR Heatmapper Plus.

Pepper plant preparation and stress treatments

Pepper plants were grown on soil in greenhouse with conditions: 14/10 h photoperiod, 25/20 °C day/night temperature and 60% relative humidity. In this study, pepper seedlings with 6–8 true leaves were randomly divided into five groups, namely control (untreated) and treatment with cold (4±1 °C), salt (300 mM NaCl), drought (400 mM mannitol) and gibberellin solution (20 μM GA). Leaves were sampled at 3 h after the treatment. For each treatment, leaves from five randomly selected seedlings were bulked to form one sample, and six biological replicate samples were immediately frozen in liquid nitrogen and then stored at −80 °C before use.

RNA isolation and qRT-PCR analysis

Total RNA from leaves was extracted using Total RNA kit (BioTeke, Beijing, China) and reversely transcribed into cDNA using M-MLV Reverse Transcriptase (Promega, Madison, WI, USA). Real-time quantitative PCR (qRT-PCR) experiment was done using SYBR GREEN I Master Mix (Applied Biosystems, Waltham, MA, USA) on iCycler iQ™ thermocycle (Bio-Rad, Hercules, CA, USA). Each reaction volume contained 12.5 μl of SYBR GREEN Mix, 1 μl of each primer, 5 μl of 10 × diluted cDNA, and 5.5 μl of nuclease-free water. The reaction program was set as follows: initial polymerase incubation at 95 °C for 10 min, then 40 cycles of 95 °C for 15 s, 60 °C for 45 s. Melting curve analysis was conducted with heating the PCR product from 60 °C to 95 °C for verifying the specificity of the primers. The relative expression levels of CaGRAS genes were calculated based on the comparative Ct method using the 2^−ΔΔCt method with the actin1 as an internal reference gene. Primer pairs were designed by Primer Premier 5.0 and checked by NCBI Primer BLAST (Table S3).

Genome-wide identification of GRAS gene family in pepper

We employed two different approaches to identify GRAS genes in pepper genome. Totally, 50 non-redundant CaGRAS genes were found from variety cv. CM334, concurrent with the corresponding genes from cv. Zunla-1 (Table 1). Nearly all these proteins contained one representative GRAS domain (PF03514.11), with the exception of three CaGRASs (CA00g84110, CA01g26680 and CA00g84090) that had more than one such domain. The molecular mass and length of CaGRAS proteins varied greatly, with MWs ranging from 48 to 87 KDa and length from 419 to 801 aa. The average theoretical pI was 6.1, implying that most CaGRAS proteins were weakly acidic. Only CaGRAS21 was stable because of its instability index less than 40, whereas the rest were considered as unstable. All CaGRASs were predicted to be hydrophilic due to the less GRAVY value (<0) of each protein. Most of CaGRAS proteins contained large percentage of aliphatic amino acids, with predicted aliphatic index ranging from 65.74 to 95.76. Interestingly, most of CaGRAS genes (84%) were intronless, while seven members had just one intron. Only one CaGRAS gene had two introns (Fig. 1).

Chromosomal localization and gene duplication analysis of CaGRAS genes

Except for six members (CaGRAS45-50) mapped to the scaffolds, the remaining 44 CaGRAS genes were unevenly distributed across 11 out of 12 pepper chromosomes. Among those anchored members, Chr7 occupied the largest number of GRAS genes (n = 7; 15.22%), followed by Chr1 (n = 6; 13.04%) and three chromosomes (Chr2, Chr5 and Chr12) each having five members. Additionally, four GRAS genes were located on Chr4, while three genes were detected on Chr3, Chr4 and Chr9, respectively. Two GRAS genes were found on Chr8, and only one was on Chr10. Notably, most of CaGRAS genes were gathered at both ends of chromosomes.

Furthermore, we analyzed the duplication events of CaGRAS gene in pepper genome since gene duplication acts importantly on the occurrence of novel functions and gene family expansion. As shown in Fig. 2, two pairs of tandem duplicated genes (CaGRAS4/5 and CaGRAS20/21) located on Chr1 and Chr5, respectively. Additionally, 10 pairs of CaGRAS genes were identified as segmental duplications (Fig. 3). We found that all duplicated gene pairs had Ka/Ks ratios less than 0.5, suggesting these genes experienced strong purifying selection pressure during evolution processes (Table 2). Clearly, segmental duplication played a more prominent role in the expansion of pepper GRAS genes than tandem duplication.

In order to understand the phylogenesis of GRAS gene family, microsynteny analysis was employed for the precise region containing GRAS genes in pepper, tomato and Arabidopsis. (Fig. 3). A total of 37, 15 and seven orthologous gene pairs were identified in the cross of pepper and tomato, pepper and Arabidopsis, tomato and Arabidopsis, respectively. Several such regions were also found between different chromosomes in pepper. These data indicate that GRAS gene family is highly conserved, and pepper CaGRAS genes are more closely to those of tomato than that in Arabidopsis. The GRAS genes with microsynteny may evolve from the same ancestor.

Phylogenetic analysis, classification and functional characterization of CaGRAS family

To uncover the evolutionary relationships among CaGRAS proteins and their classifications, we performed a phylogenetic analysis using 189 full-length GRAS proteins (32 from Arabidopsis, 56 from rice, 51 from tomato and 50 from pepper). An unrooted phylogenetic tree was constructed (Fig. 4), demonstrating that all of these GRAS proteins could be classified into eleven distinct subfamilies based on clade support values and classification from Arabidopsis and rice. They were termed as DELLA, PAT1, SCL3, SHR, SCR, LISCL, HAM, LAS, DLT, Ca_GRAS and Os4, respectively. Of these, nine were named specifically according to the findings of previous study (Tian et al., 2004). The remaining two were named by the members of species origin. For example, the subfamily Ca_GRAS consisted of six GRAS members from tomato and six GRAS proteins from pepper, indicating that it might be a Solanaceae-specific group. Os4, a rice-specific subfamily, only contained eleven GRAS members from rice.

Generally, the genes clustered into a group tend to possess similar function and structure. Therefore, we could predict the potential function of CaGRAS members based on Arabidopsis or rice homologues in the same branch. For example, within the PAT1 subfamily, AtPAT1 and AtSCL13 were previously shown to be involved in phyA and phyB signaling pathway, respectively. Hence, we inferred that CaGRAS28 and CaGRAS18 in the same subfamily may also play a significant role in phytochrome signal transduction. DELLA subfamily included two CaGRAS members (CsGRAS14 and CaGRAS41), five AtGRAS and six OsGRAS. All these members contain the complete DELLA and TVHYNPS motifs (Fig. 5). Previous studies reported that DELLA proteins mainly regulate GA signal transduction pathway (Zhang et al., 2011) implying that CaGRAS14 and CaGRAS41 may have the similar role. The SCL3 subfamily consisted of two CaGRAS members (CaGRAS2 and CaGRAS44) and one AtGRAS (AtSCL3). This subfamily may mediate GA homeostasis through integrating other signals, because AtSCL3 was found to regulate root cell elongation by integrating multiple signals in Arabidopsis (Zhang et al., 2011). For subfamilies SHR and SCR, AtSHR and AtSCR were detected to function importantly in maintaining stem cell and root meristem. It is reasonable to predict that those pepper GRAS homologs in these two subfamilies may possess the similar functions (Di Laurenzio et al., 1996). LISCL subfamily consisted of nine CaGRAS, six AtGRAS and three OsGRAS members, and their biological roles are mostly unknown although a homolog member (LISCL) from Lilium longiflorum was proven to play an important regulatory role during microsporogenesis (Morohashi et al., 2003). The first HAM gene in the HAM subfamily was isolated from petunia and proved to promote shoot indeterminacy (Stuurman, Jaggi & Kuhlemeier, 2002). CaGRAS3 in the HAM subfamily was also demonstrated to be involved in shoot apical meristem organization and axillary meristem development (David-Schwartz et al., 2013). The LAS subfamily comprised two members from pepper, three from rice and three from Arabidopsis. AtLAS proteins in this subfamily mainly function to regulate and promote the initiation of axillary meristems (Liang et al., 2014). The DLT subfamily, the smallest group, contained six members (one from pepper, one from Arabidopsis, two from rice and two from tomato). The members of this group have been previously shown to participate in brassinosteroid signal pathway responsible for the plant height (Tong et al., 2009). For the Ca_GRAS subfamily having six CaGRAS and six SlGRAS members, no Arabidopsis and rice GRAS homolog was grouped into this subfamily, indicating that these genes may be Solanaceae-specific. The function of this subfamily awaits further exploration.

To investigate the common feature of pepper GRAS proteins in more detail, we used MEME suite to identify their conserved motifs and sequence logos. A total of 11 conserved motifs (named Motif 1–11) were identified, with more motifs locating at C-terminus than at N-terminus. Moreover, the motifs from the same subfamily nearly hold the similar patterns (Fig. 5). We then matched up the motifs with corresponding GRAS domain. It was found that Motif 10 and 4 is in LHRI domain at N-terminus, followed by Motif 7 and 1 in VHIID domain, Motif 6 and 8 in LHRII domain, Motif 9, 3 and 11 in PFYRE domain, and Motif 2 and 5 in SAW domain at C-terminus (Fig. 5). Of the 10 subfamilies of CaGRAS, members from PAT1 and LISCL subfamilies all contained the 11 conserved motifs identified.

Prediction of CaGRAS protein–protein interaction network

Due to unavailable reference for pepper interactome data, we predicted the protein–protein interaction relationships of CaGRAS members based on the interologs from Arabidopsis. We only obtained the interaction information for 19 CaGRAS proteins, and generated a complex interaction network using these proteins (Fig. 6). In general, the members from the SCL3 subfamily (CaGRAS2 and CaGRAS44) owned more interaction partners than others. These were consistent with their working mechanisms, considering the fact that AtSCL3 protein could regulate GA homeostasis by integrating other signal pathway, although such a relationship needs to be confirmed (Zhang et al., 2011). CaGRAS33, a member of LAS subfamily, directly interacted with nine CaGRAS members, while CaGRAS7 from the same subfamily only had three interaction partners. Surprisingly, no interaction partner was detected for CaGRAS proteins from DELLA and DLT subfamily. Our interaction networks may provide important clues for understanding the functions of unknown proteins.

Expression analysis of CaGRAS genes in various tissues and fruit developmental stages

We used online available transcriptome data of three tissues (leaf, stem and root) and seven developmental stages of pericarp and placenta (mature green, breaker, five and 10 days post-breaker, six, 16, 25 days post-anthesis) to investigate the expression patterns of pepper GRAS genes (Fig. 7). The RPKM value for each of those CaGRAS genes was listed in Table S2. The transcripts for the other 12 CaGRAS genes were not detected in any tissues (RPKM < 0.001), which may be the result of pseudogenes. Generally, 25 CaGRAS genes were detected to express in all tissues, with only five members (CaGRAS8, CaGRAS16, CaGRAS29, CaGRAS38 and CaGRAS48) showing high expression levels (PPKM > 10). A number of CaGRAS genes exhibited a certain degree of tissue specificity. For example, CaGRAS18 and CaGRAS27 were only expressed in pepper pericarp. CaGRAS35 and CaGRAS43 were highly expressed in leaf while the transcripts of CaGRAS30 and CaGRAS34 largely accumulated in stem rather than in other tissues. Tissue-specific expression of these genes showed that they may highly participate in the corresponding tissue development. CaGRAS28 homologous with AtPAT1 showed high expression level in leaves, which is in line with AtPAT1 function as a positive regulator in phyA signal pathway (Bolle, Koncz & Chua, 2000). Several CaGRAS genes exhibited constitutive expression levels at most stages of pericarp development. For example, CaGRAS7 and CaGRAS42 displayed a relatively higher expression at green fruit stage (PC_6DPA and PC_16DPA), and then decreased gradually towards fruit ripening. This expression pattern implied that CaGRAS7 and CaGRAS42 may function importantly in the early fruit development. In addition, the similar expression patterns were often detected for gene pairs from duplication event, but not for all such genes. For instance, in the CaGRAS18/40 duplicated region, CaGRAS40 was highly expressed, whereas the other showed the opposite expression pattern. These differences implied that duplicated GRAS gene pairs may have diverged evolutionary outcomes.

Response of CaGRAS genes to different stress treatments

In order to elucidate the functions of CaGRAS genes responsive to GA stimuli, qRT-PCR was performed to examine the expression of such genes in seedling leaves after treatment with GA. In this study, 14 CaGRAS genes showed obvious changes in response to GA treatment (Fig. 8). Of them, the most upregulated gene was CaGRAS37, while the most downregulated gene was CaGRAS10. To broaden our knowledge regarding how these genes are affected by GA, we conducted a comprehensive analysis on cis-elements in the promoter regions of such 14 CaGRAS genes using PlantCARE (Lescot et al., 2002). Additionally, 12 CaGRAS genes were detected to contain at least one GA responsive element (GARE) in their promoter sequences, again confirming the function of these genes in mediating GA signal pathway in pepper (Table S4).

We further examined the expression levels of CaGRAS genes under abiotic stresses, including salt, drought and cold treatments. Compared to the control group, the expression of 12 CaGRAS genes were highly affected by these treatments, indicating that those genes may have diverse functions involving in plant responses to abiotic stresses. The downregulated expression was detected for six, two and four CaGRAS genes, respectively, under cold, drought and salt stresses. Furthermore, we found the upregulated genes exhibit a group-specific expression. For example, the expression of CaGRAS genes from the DELLA subfamily was significantly induced under cold stress. The genes in the SCL3 subfamily was highly upregulated under drought stress, and the genes in the PAT1 subfamily were highly induced by GA and other four stress treatments. Therefore, it is possible that different CaGRAS members function in different stress responses.