Genome-Wide Analysis of the GRF Family Reveals Their Involvement in Abiotic Stress Response in Cassava

GENERAL REGULATORY FACTOR (GRF) proteins play vital roles in the regulation of plant growth, development, and response to abiotic stress. However, little information is known for this gene family in cassava (Manihot esculenta). In this study, 15 MeGRFs were identified from the cassava genome and were clustered into the ε and the non-ε groups according to phylogenetic, conserved motif, and gene structure analyses. Transcriptomic analyses showed eleven MeGRFs with constitutively high expression in stems, leaves, and storage roots of two cassava genotypes. Expression analyses revealed that the majority of GRFs showed transcriptional changes under cold, osmotic, salt, abscisic acid (ABA), and H2O2 treatments. Six MeGRFs were found to be commonly upregulated by abiotic stress, ABA, and H2O2 treatments, which may be the converging points of multiple signaling pathways. Interaction network analysis identified 18 possible interactors of MeGRFs. Taken together, this study elucidates the transcriptional control of MeGRFs in tissue development and the responses of abiotic stress and related signaling in cassava. Some constitutively expressed, tissue-specific, and abiotic stress-responsive candidate MeGRF genes were identified for the further genetic improvement of crops.


Identification and Evolutionary Analysis
The GRF protein sequences in Arabidopsis and rice were acquired from UniProt and rice genome annotation project (RGAP), respectively [33,34]. The genome sequences of cassava were downloaded from the cassava genome database [35]. The known GRF sequences were used to build hidden markov models (HMM) for searching the cassava dataset using HMMER software [36,37]. The identified cassava GRFs were further validated by basic local alignment search tool (BLAST) analysis using GRFs from rice and Arabidopsis as queries [38]. Based on the conserved domains database (CDD) database, the conserved domain of cassava GRFs were validated [39,40]. To study the evolutionary relationship of cassava GRFs, MEGA 5.0 and Clustal X 2.0 softwares were employed to construct the evolutionary tree with the entire protein sequences of GRFs from cassava, Arabidopsis, and rice [41,42]. The protein sequences of cassava GRFs were also used to construct evolutionary tree and protein motif analysis.

Transcriptome Analysis
Total RNA of stems, leaves, and storage roots in Arg7 and W14 was extracted with plant RNA extraction kit (Tiangen, Beijing, China) for complementary DNA (cDNA) library construction. The sequencing was performed with an Illumina GAII (Illumina, San Diego, CA, USA) following manufacturer's instructions. Adapter sequences were removed with FASTX-toolkit [43]. Clean reads were generated by removing low quality sequences using FastQC [44]. Tophat v.2.0.10 was used to map the clean reads to the cassava genome [45]. Using cufflinks, the transcriptome data was assembled [46]. Fragments per kilobase of transcript per million fragments mapped (FPKM) was employed to calculate gene expression levels. Each sample has three replicates.

Sequence Analysis
Using the ExPASy proteomics server, the molecular weight (MW) and isoelectric points (pI) of cassava GRFs were predicted [47]. The conserved motifs of cassava GRFs were identified with MEME and InterProScan databases [48,49]. The gene structures of cassava GRFs were assessed with the gene structure display server (GSDS) database [50]. Using STRING [51] with option value >0.7, specific interaction network with experimental evidences of GRFs in Arabidopsis were constructed, which identifies 18 high confidence interactive proteins in Arabidopsis. Then, the homologs of these interactive proteins in cassava were identified with reciprocal best BLASTP analysis.

qRT-PCR Analysis
The relative expression levels of cassava GRFs were examined by qRT-PCR analysis using StratageneMx3000P (Stratagene, CA, USA) instrument and SYBR Premix Ex Taq (TaKaRa, Dalian, China). The 2 −∆∆Ct method was used to assess the relative expression of cassava GRF genes [52]. The appropriate primer pairs were selected by melting curve, agarose gel electrophoresis, and sequencing PCR products (Table S1). The amplification efficiency was in the range of 0.91-1.07. The relative expression of cassava GRF genes in each time point was calculated according to the control and treated samples. Each sample has three replicates and three biological experiments were performed.

Identification and Phylogenetic Analysis of GRF Gene Family in Cassava
Both BLAST and Hidden Markov Model searches were employed to extensively identify cassava GRF genes with GRF sequences from Arabidopsis and rice as query. As a result, 15 predicted full-length MeGRFs (MeGRF1-MeGRF15) were identified from cassava and their sequences were shown in Table S2. The cassava MeGRF proteins ranged from 232 to 264 amino acid residues in length, and their relative molecular mass varied from 26.25 kDa to 30.06 kDa, with the pIs in the range of 4.72 to 5.03 (Table S2).

Conserved Motif and Gene Structure Analyses of GRF Gene Family in Cassava
To better understand the structural divergence and prediction the function of MeGRF proteins, a total of ten conserved motifs of cassava GRFs were found using MEME and further annotated by InterPro Scan 5 [48,49] (Figure 2). Seven motifs (motifs 1, 2, 3, 4, 6, 7, 8) were annotated as 14-3-3 protein motifs, which are the basic characteristics of the GRF family. Based on the motif assay, all the identified cassava GRFs showed at least six 14-3-3 protein motifs, indicating their typical characteristic of GRF family. Interestingly, some close homologous GRFs exhibited the same motif

Conserved Motif and Gene Structure Analyses of GRF Gene Family in Cassava
To better understand the structural divergence and prediction the function of MeGRF proteins, a total of ten conserved motifs of cassava GRFs were found using MEME and further annotated Genes 2018, 9,  by InterPro Scan 5 [48,49] (Figure 2). Seven motifs (motifs 1, 2, 3, 4, 6, 7, 8) were annotated as 14-3-3 protein motifs, which are the basic characteristics of the GRF family. Based on the motif assay, all the identified cassava GRFs showed at least six 14-3-3 protein motifs, indicating their typical characteristic of GRF family. Interestingly, some close homologous GRFs exhibited the same motif organization, including MeGRF1/3 specially having motif 9, MeGRF11/12 uniquely showing motif 5, and MeGRF8/10 uniquely containing motif 10, which further supports the phylogenetic classification of GRF family.

Figure 2.
Phylogenetic and motif analyses of cassava GRFs. All cassava GRF protein sequences were determined with MEME software and were classed into two groups based on phylogenetic relationship [48].
Additionally, exon-intron organization among the coding sequence was examined to understand the structural diversity and evolution of cassava GRF family. As shown in Figure 3, it is clear that the ε group contained 4-7 exons, whereas the non-ε group showed 4 exons. The distinct exon-intron organization between the ε group and the non-ε group indicated their diversity during evolution, further supporting the phylogenetic classification of cassava GRFs.

Expression Profiles of MeGRF Genes in Different Tissues
To study the expression profiles of MeGRF genes in different tissues, transcriptomic analyses were carried out from samples of leaves, stems, and storage roots in a wild subspecies (W14) and cultivated variety Arg7 ( Figure 4; Table S4). The results revealed that 11 GRFs (MeGRF-1, -4, -5, -6, -7, -8, 10, -11, -12, -13, -15) showed constitutively high expression (FPKM > 29) in all the tested tissues of the two genotypes. Interestingly, MeGRF2 showed tissue specific expression pattern, with transcripts lacking in leaves and low expression (FPKM < 10) in stems and storage roots of the two genotypes. MeGRF9 showed high expression in stems and storage roots of Arg7, but low expression in leaves of Arg7 and in all the tested tissues of W14. MeGRF14 had high expression in storage roots of Arg7, whereas low expression in other tissues tested. These results implied the same or differential roles of cassava GRF genes in tissue development.

Figure 2.
Phylogenetic and motif analyses of cassava GRFs. All cassava GRF protein sequences were determined with MEME software and were classed into two groups based on phylogenetic relationship [48].
Additionally, exon-intron organization among the coding sequence was examined to understand the structural diversity and evolution of cassava GRF family. As shown in Figure 3, it is clear that the ε group contained 4-7 exons, whereas the non-ε group showed 4 exons. The distinct exon-intron organization between the ε group and the non-ε group indicated their diversity during evolution, further supporting the phylogenetic classification of cassava GRFs.

Expression Profiles of MeGRF Genes in Different Tissues
To study the expression profiles of MeGRF genes in different tissues, transcriptomic analyses were carried out from samples of leaves, stems, and storage roots in a wild subspecies (W14) and cultivated variety Arg7 ( Figure 4; Table S4). The results revealed that 11 GRFs (MeGRF-1, -4, -5, -6, -7, -8, 10, -11, -12, -13, -15) showed constitutively high expression (FPKM > 29) in all the tested tissues of the two genotypes. Interestingly, MeGRF2 showed tissue specific expression pattern, with transcripts lacking in leaves and low expression (FPKM < 10) in stems and storage roots of the two genotypes. MeGRF9 showed high expression in stems and storage roots of Arg7, but low expression in leaves of Arg7 and in all the tested tissues of W14. MeGRF14 had high expression in storage roots of Arg7, whereas low expression in other tissues tested. These results implied the same or differential roles of cassava GRF genes in tissue development.

Differential Expression of MeGRF Genes in Response to Cold, Osmotic, Salt, ABA, and H2O2 treatments
To investigate the transcriptional responses of MeGRF genes to abiotic stress and related signaling, cassava seedlings were subjected to cold, osmotic, salt, ABA, and H2O2 treatments (           expression levels. NTC indicates no treatment controls (mean value = 1). Data are means ± SD of n = 3 biological experiments. Means denoted by the same letter do not significantly differ at p < 0.05 as determined by Duncan's multiple range test.

Analysis of GRF Family Interaction Network
To characterize the possible interaction networks of cassava GRFs, detailed analysis of GRFs from Arabidopsis and cassava were performed. Firstly, we applied STRING database [51] to construct the interaction networks of GRFs in Arabidopsis. We found that 13 GRFs were involved in the interaction network and these GRFs interacted with 18 target proteins, including Omethyltransferase, H + -ATPases, RING/FYVE/PHD zinc finger-containing protein, protein kinases, protein phosphatase 2C, and ascorbate peroxidase ( Figure 10; Table S5). Secondly, the homologs of these proteins involved in the interaction network were identified from cassava with reciprocal best BLASTP analysis ( Figure  10; Table S6). Thus, the potential interaction networks of GRFs were constructed in cassava. These results would lay a foundation for further investigation of GRF functions in cassava.

Analysis of GRF Family Interaction Network
To characterize the possible interaction networks of cassava GRFs, detailed analysis of GRFs from Arabidopsis and cassava were performed. Firstly, we applied STRING database [51] to construct the interaction networks of GRFs in Arabidopsis. We found that 13 GRFs were involved in the interaction network and these GRFs interacted with 18 target proteins, including O-methyltransferase, H + -ATPases, RING/FYVE/PHD zinc finger-containing protein, protein kinases, protein phosphatase 2C, and ascorbate peroxidase ( Figure 10; Table S5). Secondly, the homologs of these proteins involved in the interaction network were identified from cassava with reciprocal best BLASTP analysis ( Figure 10; Table S6). Thus, the potential interaction networks of GRFs were constructed in cassava. These results would lay a foundation for further investigation of GRF functions in cassava.

Discussion
Due to the regulatory role of GRFs in plant growth, development and response to abiotic stress, and the limited information for this gene family in cassava, it is essential to study the possible role of GRF genes in cassava. Here, we identified 15 GRFs from the cassava genome, which is expanded in comparison to GRFs from rice, Arabidopsis, common bean, tomato, P. trichocarpa, and B. distachyon [2,[9][10][11][12][15][16][17]19]. Phylogenetic analysis showed that cassava GRFs were classified into the ε group and the non-ε group (Figure 1). This is consistent with previous classification of GRFs from Arabidopsis, rice, common bean, etc. [2,9,15]. Exon-intron organization analysis suggested that the ε group MeGRFs had more exons and introns than the non-ε group (Figure 3). This phenomenon is also observed in Arabidopsis, rice, common bean, M. truncatula, and B. rapa [2,9,15,18,53]. Besides, the exon number of GRFs is conserved among cassava (4-7), Arabidopsis (3)(4)(5)(6), and rice (4-7). Conserved motif analysis showed that at least six 14-3-3 protein motifs existed in both ε group and non-ε group of cassava GRFs, indicating their conservation of protein sequences (Figure 2). Together, these evidences support the classification and conservation of cassava GRF family.

Discussion
Due to the regulatory role of GRFs in plant growth, development and response to abiotic stress, and the limited information for this gene family in cassava, it is essential to study the possible role of GRF genes in cassava. Here, we identified 15 GRFs from the cassava genome, which is expanded in comparison to GRFs from rice, Arabidopsis, common bean, tomato, P. trichocarpa, and B. distachyon [2,[9][10][11][12][15][16][17]19]. Phylogenetic analysis showed that cassava GRFs were classified into the ε group and the non-ε group (Figure 1). This is consistent with previous classification of GRFs from Arabidopsis, rice, common bean, etc. [2,9,15]. Exon-intron organization analysis suggested that the ε group MeGRFs had more exons and introns than the non-ε group (Figure 3). This phenomenon is also observed in Arabidopsis, rice, common bean, M. truncatula, and B. rapa [2,9,15,18,53]. Besides, the exon number of GRFs is conserved among cassava (4-7), Arabidopsis (3-6), and rice (4)(5)(6)(7). Conserved motif analysis showed that at least six 14-3-3 protein motifs existed in both ε group and non-ε group of cassava GRFs, indicating their conservation of protein sequences ( Figure 2). Together, these evidences support the classification and conservation of cassava GRF family.
Accumulated evidences have revealed the important functions of 14-3-3 proteins in plant growth and development [20][21][22][23][24][25]. Investigation of the tissue expression patterns of 14-3-3 proteins would provide some clues on tissue development. In this study, we found that MeGRF-6, -8, -10, and -15 in the ε group and MeGRF-1, -4, -5, -7, -11, -12, and -13 in the non-ε group exhibited constitutive high expression levels (FPKM > 29) in leaves, stems, and storage roots of W14 and Arg7 (Figure 4). In Arabidopsis, AtGF10 from the ε group, and AtGF4 and AtGF6 from the non-ε group showed abundant expression in shoots and roots [3]. Their homologous genes MeGRF15 (homologs of AtGF10), MeGRF5 (homologs of AtGF4), and MeGRF-11, -12, -13 (homologs of AtGF6) in cassava also had high expression (Figure 1; Figure 4). In B. distachyon, seven out of 8 eightGRF genes showed low expression in roots, whereas high expression in stems, leaves, and spikelets [19]. In banana, five out of 25 GRF genes had constitutive high expression in roots, leaves, and fruits of two varieties [6]. In M. truncatula, eight GRF genes exhibited lower expression in leaves compared with in other tissues of roots, shoots, and flowers [53]. In common bean, most of the GRFs displayed high expression in flowers and stems, while low expression in pods and leaves [15]. In mesohexaploid B. rapa, most of the GRF genes showed low expression abundance in various tissues of roots, stems, leaves, and siliques [18]. Collectively, these studies revealed the tissue expression diversity of GRFs in various plant species. Compared with GRFs in other plants, the great number of GRFs with constitutive high expression in cassava indicates their important function in cassava development.
Biochemical and genetic evidences also confirmed the regulatory role of GRFs in plants response to abiotic stress and hormones [19,[26][27][28][29]54,55]. To better understand cassava GRFs mediated transcriptional responses under abiotic stress and related signaling, MeGRFs expression were examined under various treatments. The results showed that MeGRFs could widely respond to cold, osmotic, salt, ABA, and H 2 O 2 treatments at transcriptional levels (almost half members showing induction and half members showing repression under each treatment), suggesting their potential role in abiotic stress response (Figures 5-9). The significant changes of GRFs at the transcriptional level under abiotic stress were also observed in other plants. In banana, nine, 13, and twelve MaGRF genes showed induction after cold, salt, and osmotic treatments, respectively, whereas 10, six, and seven MaGRF genes showed repression under the corresponding treatments [6]. In B. distachyon, GRF genes could be transcriptionally induced or repressed after osmotic, salt, ABA, and H 2 O 2 treatments [19]. In common bean, all the identified GRFs were upregulated after cold treatment, and were induced or repressed upon drought and salt treatments [15]. In mesohexaploid B. rapa, most of the GRFs showed upregulation after salt, ABA, or SA treatments, but downregulation after dehydration or heat treatments [18]. These evidences are in accord with our expression analysis of cassava GRFs, further supporting the possible role of MeGRFs in abiotic stress responses.
Notably, MeGRF3, MeGRF5, MeGRF6, and MeGRF11 were commonly upregulated by salt, osmotic, and ABA treatments; MeGRF4 was commonly upregulated by cold, salt, ABA, and H 2 O 2 treatments; and MeGRF10 was commonly upregulated by cold, osmotic, and H 2 O 2 treatments (Figures 5-9). Numerous evidences indicated the positive role of GRFs in plants response to abiotic stress through affecting ABA pathway, stomatal behavior, ROS balance, and ions transports [19,[26][27][28][29]. Thus, these cassava GRFs may be the converging points when cassava responds to abiotic stress, ABA signaling and H 2 O 2 , and can serve as candidates for genetic improvement of crop tolerance to abiotic stress.
As an important regulatory factor, GRFs elaborate their function through interacting with different clients. There is a need to investigate GRF mediated interaction network. In this study, we predicted 18 possible targets of 13 GRFs in cassava, including enzymes, transporters, protein kinases, and transcription factors ( Figure 10; Table S6). These interactions have been confirmed in Arabidopsis [2,[4][5][6][7]. The interaction relationship between GRFs and their targets, and the expression of cassava GRFs in each interaction group were shown in Table S7. This provided some clues for investigating the expression and function of GRFs and their targets. Further experimental validations would deepen the understanding of GRF functions in cassava.
In conclusion, this study identified 15 GRFs from cassava and investigated their phylogenetic classification, protein motif, and gene structure. Transcriptomic analysis showed the constitutively expressed or tissue specifically expressed MeGRFs. Expression analysis revealed the involvement of MeGRFs under abiotic stress and related signaling and identified some important candidates for improving crop resistances to multiple stresses. Furthermore, the GRF mediated interaction network was characterized, which would facilitate further study of their function in cassava. This systematic study will advance the understanding of GRF-mediated signal cascades in regulating cassava development and abiotic stress response, thereby supplying candidates for crop breeding.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/9/2/110/s1. Table S1: The primer sequences used for qRT-PCR. Table S2: Characteristics of GRFs in cassava. Table S3: The accession number of GRFs in Arabidopsis and rice. Table S4: Expression data of cassava GRFs in different tissues of W14 and Arg7. Table S5: The protein interaction relationship in the GRF-mediated interaction network in Arabidopsis. Table S6: The homologous genes in cassava from the interaction network. Table S7. The protein interaction relationship in the GRF-mediated interaction network and the expression of GRFs in cassava.