Characterization of the GRAS gene family reveals their contribution to the high adaptability of wheat

Genome-wide detection and sequence analysis of TaGRAS

The CDS and protein sequences of both high and low confidence bread wheat (Triticum aestivum cv. Chinese Spring) genes were downloaded from https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.0/. High-confidence genes indicated genes with complete gene models (genes contain both start and stop codons) and very good sequence homology to plant proteins from SwissPort and Poaceae proteins in trEMBL but no good sequcence homology in transposon database TREP. Low-confidence genes indicated that genes which have incomplete gene models but high identity to plant proteins from SwissProt; or genes with complete gene models but no significant homology to TREP, UniPoa or Unimagdatabases; or incomplete gene models with very high identity to the annotated Poaceae protein in SwissProt or trEMBL but not in the TREP database. Besides, low confidence genes including genes with incomplete gene models and no significant homology to any genes in SwissProt, trEMBL or TREP database. The complete and detailed parameters can be referred to the previous literature (IWGSC, 2018). The GRAS protein sequences of Arabidopsis (32), rice (Oryza sativa, 56), barley (Hordeum vulgare L., 62) and Brachypodium distachyon (48) were retrieved from previous studies (Niu et al., 2019; Tian et al., 2004; To et al., 2020). These sequences were merged together to query bread wheat proteins by BLASTp. Proteins with Evalue ≤ 1e−5, identity ≥60 and length coverage ≥50% were used to detect the exaistance of GRAS domains with HMMER 3.0 and the Hidden Markov Model (HMM) profile of the GRAS domain (PF03514, http://pfam.xfam.org/, Evalue ≤1e−5) (Finn et al., 2016; Mistry et al., 2013). Then, these candidate GRAS genes were examined to confirm the existence of GRAS domain in SMART (http://smart.embl-heidelberg.de/), PFAM (http://pfam.xfam.org/) and Conserved Domain Database (CDD) in NCBI (https://www.ncbi.nlm.nih.gov/). Proteins with truncated, lacking, or anonymous GRAS domain were eliminated. The confirmed GRAS genes were named based on an abbreviation for the species name Triticum aestivum (Ta) and their order on the chromosomes. The subcellular localization of GRAS proteins were predicted in WoLFPSORT (https://wolfpsort.hgc.jp/). Exon-intron structures were constracted in Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn). ExPASy (https://web.expasy.org/compute_pi/) was used to calculate isoelectric point and molecular weight. Phylogenetic tree was constructed with MEGA 5.0 based on the maximum likelihood method using 1000 bootstraps (Hall, 2013). The phylogenetic tree was subsequently visualized in iTOL (https://itol.embl.de/). Homoeologous genes between the A, B and D subgenome of wheat were detected based on the phylogenetic tree and previous classifications (Ramírez-González et al., 2018; Schilling et al., 2020). The chromosomal location was displayed with circos (Krzywinski et al., 2009). Ka (non-synonymous substitution) and Ks (synonymous substitution) was calculated using ParaAT (Parallel Alignment and back-Translation) software with MA methods (Zhang et al., 2012).

Expression profiles in different tissues and under stress conditions

The expression level (TPM, transcripts per million) of TaGRAS genes were downloaded from http://www.wheat-expression.com/, and the heatmap was generated with R packages “pheatmap” (scale=“none”, cluster_rows=F, cluster_cols=F). TaGRAS genes were classified as 10 clusters with kmeans according to expression level, and then the clusters were divided into groups that preferentially expressed in reproductive, vegetative or ubiquitous tissues. Homoeolog expression patterns were performed as previously described (Krzywinski et al., 2009). Briefly, TPM of triad with 1:1:1 ratio between the three subgenomes was extracted, and the triads with summed TPM lower than 0.5 were discarded. The euclidean distance of normalized expression level of triads was calculated and the triads were classified based on the shortest distance as A∕B∕D dominant, A∕B∕D suppressed or balanced profiles. Ternary plot was generated with R package “vcd”.

ChIP-seq data analysis

The ChIP-seq raw data (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, H3K9ac and H3K9me2) were downloaded from Gene Expression Omnibus Database (accession: GSE121903) (Li et al., 2019b). Trimmomatic was used to remove the reads with low quality. Then the clean reads were mapped to the reference sequence of bread wheat (IWGSC, v1.0) with bowtie2 (Bolger, Lohse & Usadel, 2014; IWGSC, 2018). To ensure the accuracy of results, reads from PCR duplication or with mapping quality below 20 were removed. The reads-enriched regions (peaks) of each mark were detected with MACS (v1.3.7) (Zhang et al., 2008). The peaks with at least 1 base pair overlapped between biological replicates were identified as the credible peaks. Target genes were defined as genes with a peak in the gene body. The TPM data of corresponding stage seedling were downloaded from http://www.wheat-expression.com/, (Table S7). Genes with TPM higher than 0.5 were recorded as expressed genes.

Degradome data analysis

To detected the targets of miR171, the degradome data of 21 and 28 DPA (days post anthesis) grains were downloaded from NCBI Gene Expression Omnibus (accession: GSE65799) (Zhang et al., 2015a). The degradome data of seed embryos, seedling leaves, seedling roots and grains of 8 days after pollination were downloaded from NCBI SRA database (accession: SRP040143) (Sun et al., 2014). First of all, low-quality reads were removed with Trimmomatic (Bolger, Lohse & Usadel, 2014). CleaveLand4 was used to predictthe targets of miR171, the targets with alignment score above 4.5 were removed (Addo-Quaye, Miller & Axtell, 2009).

Detection of TaGRAS in bread wheat

In this study, 188 TaGRAS genes were detected and their molecular weights varied from 26.69 to 164.81 kDa (Table S1). These genes were named from TaGRAS1 to TaGRAS188 according to their chromosomal location. The longest TaGRAS protein was TaGRAS82, which had 819 amino acids, while the shortest was TaGRAS172 with only 245 amino acids. The subcellular localization of TaGRAS proteins were predicted by WoLFPSORT. 42.0% (79/188) of TaGRAS proteins were predicted to be located in nucleus. As previously confirmed, TaGRAS179 (TaMOC1-7A) was predicted to be located in the cell nucleus (Zhang et al., 2015b). The isoelectric point of TaGRAS proteins ranged from 4.73 to 10.21, of which 158 were below 7.00, suggesting that most TaGRAS were slightly acidic, and different TaGRAS proteins played roles in different microenvironment.

Phylogenetic and gene structure analysis of TaGRAS gene families

To analyze the evolutionary relationship of the GRAS protein in Arabidopsis, B. distachyon, rice and wheat, a phylogenetic tree was constructed using the maximum likelihood method with MEGA5.0 (Hall, 2013) (Fig. 1). TaGRAS genes were divided into 12 subfamilies: SHR (containing 12 members), PAT1 (20), DELLA (18), LISCL (70), SCL3 (16), DLT (3), SCR (12), LAS (6), SCL4/7 (3), HAM (18), Os4 (7) and Os19 (3). Consistent with the exon-intron structure in grapevine, rice and tomato, 60.0% (114/188) of TaGRAS genes had no intron (Fig. S1) (Huang et al., 2015). The exception was the subfamily PAT1, which had a significantly higher number of exons than other subfamilies (Fig. 2A). When considering the Ka/Ks, none of the TaGRAS genes were under positive selection. The genes from subfamily DELLA, LAS, SCL4/7, SCR and SHR were under strict purifying selection pressure compared with PAT1and LISCL (Fig. 2B).

Chromosomal location of TaGRAS

A total of 188 TaGRAS genes unevenly distributed on wheat chromosomes with homoeologous group 4 had a higher density of genes, which was significantly higher than expected (56/188, 29.8%; χ²test, P < 0.001; Figs. 3A, 3B). The enrichment in homoeologous group 4 was mainly the result of tandem duplication (two TaGRAS genes were separated by no more than three genes, Table S1). For example, TaGRAS89-94, TaGRAS111-116 and TaGRAS129-135 were tandemly duplicated clusters on chromosomes 4A, 4B and 4D respectively. In addition, the presence of TaGRAS in chromosomal segment C was significantly less frequent than all wheat genes (data from Ramírez-González et al., 2018), as only two subfamilies had genes located in this region, indicating that they were preferentially distributed on the distal end of chromosomes (1.60% vs 10.67%; χ²test, P < 0.001; Table S2). In contrast, in subfamily LAS, no gene was located in distal telomeric segments (segment R1 and R3, Table S2), indicating a subfamily-specific location pattern.

TaGRAS genes exhibit a high rate of homoeolog retention

When considering the total number of GRAS genes, the wheat GRAS gene number was higher than that of B. distachyon, rice or Arabidopsis even when polyploid level was considered (188/3>48, 56, 32, Fig. 3C). The increase was mainly the result of amplification of the LISCL, DELLA and PAT1 subfamilies. Gene duplication events contributed to the expansion of these subfamilies. For example, the BRADI1G15123 gene from the LISCL subfamily of B. distachyon was sister to 23 wheat genes. This was the result of segmental and tandem duplication in the ancestor species of Triticum (before polyploidization of wheat). In some cases, the expansion was supported by low-confidence genes, such as BRADI3G24210, which was sister to four wheat genes, including one low-confidence gene (TaGRAS160, TraesCS5B01G604500LC, Fig. 1).

To further explore the evolution of TaGRAS genes during the polyploidization of wheat, we detected the homoeologous relationship in detail. As shown by the phylogenetic tree, in 71.4% (35/49) of the subclades, one B. distachyon gene and one rice gene were closely related to a triad of wheat homoeologs (genes had a 1:1:1 correspondence across the three A, B and D subgenomes) (Fig. 1). For example, Os03g15680 and BRADI1G67340 were sister to the triads TaGRAS88 (TraesCS4A01G088500), TaGRAS117 (TraesCS4B01G215700) and TaGRAS136 (TraesCS4D01G216200). We noted that TaGRAS genes of subfamilies SCL4/7, LAS, DLT and SCR all presented as triads. Moreover, compared with all of annotated wheat genes (data from Ramírez-González et al., 2018), TaGRAS had a higher proportion of genes belonging to triads (79.8% vs 35.8%) (Table S3), suggesting a high rate of homoeolog retention during hexaploidization.

Gene expression profiles of TaGRAS in various tissues

GRAS proteins are transcriptional regulators that play important roles in plant growth and development (Cheng, 2004; Fleck & Harberd, 2002; Hirsch & Oldroyd, 2014; Tyler et al., 2004). To detect the expression pattern of TaGRAS, RNA-seq data were analyzed. Approximately 69.7% of TaGRAS were expressed in at least one developmental stage (131/188, TPM>0.5, Fig. 4A, Table S4). The majority of genes from the same subfamily displayed similar expression patterns. In general, the gene expression levels of subfamily PAT1, except for TaGRAS149, were much higher than those of other subfamilies (Fig. 4A, Table S4). In addition, genes from subfamily PAT1 were expressed at high lever in all of the investigated tissues, although there was a slight decrease in grains (Fig. 4A, Table S4). Os19, DLT and SHR subfamily genes were preferentially expressed in roots, consistent with their roles in root development. All of the genes from LAS were hardly expressed in different tissues except for the quite low level in roots (Fig. 4A, Table S4). We further hierarchically clustered the genes according to TPM as previously described (Schilling et al., 2020). Fig. S2 showed that the largest subfamily LISCL contained genes from nine clusters, which were preferentially expressed in ubiquitous, reproductive, vegetative tissues or even not expressed. Subfamily DLT preferentially expressed in vegetative tissues (Fig. S2).

As 79.8% of TaGRAS presented as triads, we explored the homoeolog expression bias in leaf/shoot, root, spikelet and grain. We mainly focused on genes with summed expression levels above 0.5 TPM across the traid. According to the expression level and the euclidean distance of A, B and D subgenome homoeologs, the expression bias was classified as balanced, A/B/D-surpressed or A/B/D-dominant patterns. Consistent with all of wheat genes (data from Ramírez-González et al., 2018), 62.4% of the triads of roots showed similar relative abundances from the A, B and D subgenome, namely, the balanced category (Table S5) (Ramírez-González et al., 2018). In reproductive tissue spikelets and grains, approximately 60.7% of the triads were assigned to the balanced category (Figs. 4B, 4C, Fig. S2 and Table S5). In addition, in contrast to previous reports, a higher proportion of D-homoeolog suppression was observed across all the tissues (16.2% vs 5.7%). That mainly results from three triads: PAT1 (TaGRAS148-TaGRAS160-TaGRAS166), DELLA (TaGRAS100-TaGRAS105-TaGRAS123) and LISCL (TaGRAS150-TaGRAS161-TaGRAS167), which displayed D-homoeolog suppression in all detected tissues (Table S5).

TaGRAS genes play roles under biotic and abiotic stresses

To explore the function of TaGRAS genes in response to stresses, the TPM data of bread wheat genes under various treatments were downloaded and analyzed (Table S6). We noted that five LISCL homoeolog genes, TaGRAS178, TaGRAS180, TaGRAS181, TaGRAS185 and TaGRAS186, were induced by various stress treatments (Fig. S3). They were gradually induced after inoculation with Fusarium graminearum (accession number, ERP013829, Fig. 5A). In addition, they showed significant upregulation in leaves and shoots when infected with stripe rust pathogens, powdery mildew pathogens, chit or flg22 (Fig. 5A). Moreover, abiotic stress, such as phosphorous starvation, heat and drought stress or PEG6000 treatment, induced their expression (Fig. 5A). These results imply that these genes may play an important role in various diseases and abiotic resistance. Three genes of subfamily SCL4/7 were mainly expressed in roots and upregulation was detected in leaves after 6 h of heat stress treatment (Fig. 5B). A low-confidence gene from subfamily DELLA, TaGRAS147, was expressed at a very low level in the developmental stages, but showed upregulation after 6 h of heat, drought and both stresses, suggesting that the authenticity of the gene (Fig. 5C).

Epigenetic modification regulated the expression of TaGRAS

According to a previous study, GRAS genes are the targets of miR171 (Llave et al., 2002). Using degradome data of grains, TaGRAS142, TaGRAS155, TaGRAS169 and TaGRAS171 were predicted to be targets of miR171 (Fig. 6A, Fig. S4 and Table S4). In leaves, miR171 was found to target TaGRAS104, TaGRAS175 and TaGRAS183 (Table S4). In addition to miRNAs, other epigenetic modifications play important roles in regulating the expression of genes. Here, we analyzed genes regulated by histone modifications in leaves. Consistent with the function of H3K9me2 in inhibiting the activation of transposons in centromere regions, none of the TaGRAS genes were modified by H3K9me2 (Table S7) (Ayyappan et al., 2015; Kowar et al., 2016). As a repressive mark, H3K27me3 modified 25.0% (32/128) of the nonexpressed TaGRAS (TPM<0.5), such as TaGRAS37 (Fig. 6B, Table S7) (Molitor et al., 2014). In contrast, approximately half of the expressed TaGRAS were targets of H3K36me3, H3K4me3 and H3K9ac (48.3%, 50.0% and 50.0%, respectively, Table S7). This result was comparable with their enrichment at transcript start sites of activated genes (Du et al., 2013; Liu et al., 2010; Yang, Howard & Dean, 2014). For example, Rht-B1 was modified by H3K4me3, H3K36me3, H3K9ac and H3K27ac in leaves (Fig. 6C, Table S7). Several nonexpressed genes were modified by neither miRNA nor histone modifications. These genes may be regulated by DNA methylation (Fig. 6D). Certainly, genetic factors are also crucial for gene expression regulation, such as cis regulatory elements and mutations, which require further study.