Genome-wide in silico identification , characterization and transcriptional analysis of the family of growth-regulating factors in common bean ( Phaseolus vulgaris l . ) subjected to polyethylene glycol-induced drought stress

According to most recent findings, growth regulating factors (GRFs) are plant-specific transcription factors (TFs) that play important roles in many processes, including abiotic and biotic stress response mechanisms. Completion of the common bean (Phaseolus vulgaris) genome project has provided researchers with the opportunity to identify all GRF genes in this species. With this aim, a genome-wide in silico study was performed and 10 GRF proteins (called PhvGRFs) were identified in the common bean genome. Conserved and mandatory motifs (QLQ and WRC) were confirmed in all identified PhvGRFs and two segmental duplication events were determined. Most of the PhvGRFs were found to be more similar to Arabidopsis thaliana GRFs than to Zea mays GRFs in a phylogenetic tree. According to the expression analysis of 10 PhvGRFs, inversely related expression patterns were observed in the roots of Yakutiye and Zulbiye cultivars based on their capacity to adopt to drought stress. After drought treatment of the Zulbiye cultivar, a drought-sensitive common bean cultivar, PhvGRF1, PhvGRF2, PhvGRF3, PhvGRF5, PhvGRF6, PhvGRF9 and PhvGRF10 genes were upregulated 2to 4-fold in root tissues, as compared to the untreated control. The trend of PhvGRF1, PhvGRF2, PhvGRF3, PhvGRF5, PhvGRF6, PhvGRF7, PhvGRF9 and PhvGRF10 genes showed a consistent decline of 2to 6-fold in root tissues of the drought-tolerant Yakutiye cultivar subjected to 24 h of drought stress. We demonstrated that the expression patterns of the identified PhvGRFs correlated with the drought-stress response in a cultivar-specific manner in the common bean. We suggest that members of the GRF family can also be used for genetic engineering applications in the common bean.


INTRODUCTION
Growth regulation factors (GRFs), plant-specific proteins classified in the transcription factor (TF) family of proteins, have important roles in the control of leaf, floral organ and root development and plant longevity [1][2][3][4].In addition, recent studies have established that GRFs play an important role in abiotic and biotic stress response mechanisms [5,6].The first member of the GRFs was identified in rice (OsGRF1) 15 years ago [1].To date, many genome-wide in silico studies have been performed to identify GRF members in eudicots, monocots and other plant species [2,5,[7][8][9][10][11][12][13].According to previous studies, all members of GRF family contain conserved QLQ (glutamine, leucine, glutamine) and WRC (tryptophan, arginine, cysteine) regions, which are found in the N-terminal of GRF members.QLQ motif-bearing proteins are found in all eukaryotes, while WRC is a plant-specific protein motif.TQL is an another motif that is a semi-conserved region among all GRF members [1].
AtGRF7, a member of the GRF family, acts as a repressor of stress defense genes under normal conditions, and its expression was decreased when Arabidopsis thaliana was subjected to osmotic stress [5].AtGRF1 and AtGRF3 genes play an important role in the response to biotic stress caused by the Heterodera schachtii nematode [6].However, whether specific functions of plant GRFs are associated with stress response is currently unknown for many plant species.tions of GRF family members have been performed for different plant species to date [2,4,7,14].Herein we identified the GRF members in common bean (Phaseolus vulgaris) using a computational method for genome-wide assessment for the first time.The expression levels of the identified GRFs were analyzed using qRT-PCR in root and leaf tissues of two different common bean cultivars -drought-tolerant and a drought-sensitive plants -that were subjected to moderate polyethylene glycol (PEG) stress.

Determination of conserved motifs and phylogenetic analysis
MEME software (http://meme-suite.org/,Multiple Em for Motif Elicitation) was used to identify conserved motifs among all confirmed PhvGRF proteins [18].Subsequently, all protein sequences were aligned with ClustalW using MEGA-6 software [19]; phylogenetic analysis was performed with a neighbor-joining tree (used parameters: Poisson correction, bootstrap analysis with 1000 replicates and pairwise deletion) [11,20].
Zea mays and Arabidopsis thaliana GRF proteins were retrieved from the Phytozome database and aligned with ClustalW using MEGA-6 software.Comparative phylogenetic analysis between ZmGRFs (14 protein sequences), AtGRFs (9 protein sequences) and PhvGRFs (10 protein sequences) was performed by neighbor-joining tree (used parameters: Poisson correction, bootstrap analysis with 1000 replicates and pairwise deletion).

Chromosomal distribution and gene structure of PhvGRFs
Gene duplications were checked according to previously described criteria [21], and chromosomal distributions of PhvGRF genes based on their physical position (bp) obtained from the Phytozome database, were plotted using MapGene2Chrom web v2 (http://mg2c.iask.in/mg2c_v2.0/)[21].The Gene Structure Display Server (http://gsds.cbi.pku.edu.cn) and NCBI ORF finder tool (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) were used to perform structural analyses (exon intron numbers, locations and conserved domain locations) and to determine the open reading frames (ORFs) of the PhvGRF genes, respectively.The Prot-Param tool (http://web.expasy.org/protparam/) was used to compute the physiochemical characteristics (amino acid number, molecular weight and theoretical isoelectric point − pI) of identified PhvGRF genes.

Functional analysis of PhvGRF genes
Ten identified PhvGRF genes were assessed according to their molecular functions, biological process and cellular localizations using Blast2GO functional annotation and Genomics software [22].

Growth of plants and PEG application
Two common bean cultivars, the drought-tolerant Yakutiye and drought-sensitive Zulbiye, were selected for comparative expression analysis of identified Ph-vGRF genes.Seeds of the bean cultivars were germinated and grown hydroponically in pots containing 0. 3-, 100 μM Fe and 1 μM B 3+ .Common bean seedlings were incubated in a controlled environmental growth chamber (25°C day/20°C night, 16-h light/8-h dark photoperiod with 300 μmol m −2 s −1 light intensity) with relative humidity from 55 to 70%.Drought stress was applied using Hoagland's solution containing 100 mM PEG (for moderate PEG-induced drought stress) for 24 h after the seedlings reached the first trifoliate stage in the growth chamber.Following stress application, the root and leaf tissues of the two different common bean cultivars were collected to be used in qRT-PCR analysis.Three biological replicates were used for the qRT-PCR reactions.

RNA extraction and cDNA synthesis
Total RNA extraction was performed using a Nu-cleoSpin RNA kit (Macherey-Nagel, Germany) according to the kit protocol.RNA quantity/quality was measured with a Nanodrop ND-Spectrophotometer Lite (Thermo Scientific, USA) and was also confirmed by gel electrophoresis in 1.5% agarose.cDNA synthesis was performed with 2 μg of RNA, 2.5 μM anchored-oligo(dT)18, 1X Transcriptor High Fidelity Reverse Transcriptase Reaction Buffer, 20 U Protector RNase Inhibitor, 1 mM deoxynucleotide Mix, 5 mM DTT, and 10 U Transcriptor High Fidelity Reverse Transcriptase using the High Fidelity cDNA Synthesis Kit (Roche).The following incubation conditions were applied: 10 min at 65 o C, 30 min at 55 o C and 5 min at 85 o C. cDNA quantity/quality was also measured with a Nanodrop ND-Spectrophotometer Lite.

Quantitative Real-Time PCR
Quantitative real-time PCR (qRT-PCR) was performed with a Light Cycler® Nano System (Roche) thermal cycler.The primer sequences of the target genes and the housekeeping gene, which is used for normalization, were designed with the Primer3 program based on the sequences of predicted PhvGRFs (Supplementary Table S1).All qRT-PCR reactions were performed in three independent biological and technical triplicates with a template-free control to check for contaminations.Amplifications of PCR products were monitored using SYBR Green I dye.After predenaturation for 10 min at 95 o C, 45 cycles of 15 s at 95 o C, 20 s at 60 o C and 20 s at 72 o C were applied.Melting curve analysis was performed to confirm the presence of a single product and absence of primerdimers.Data collection for quantification was done during the annealing period.

Statistical analysis
The abundance of target gene transcripts was normalized to the Actin gene (ACT) and set relative to the control plants according to the 2 -∆∆CT method [23].Changes in relative expression levels (REL) of the gene were checked for statistical significance according to one-way analysis of variance (ANOVA).Fisher's least significant difference test (at 0.05 significant level) was performed.

RESULTS
GRF proteins of Arabidopsis thaliana, Carica papaya, Physcomitrella patens, Populus trichocarpa, Selaginella moellendorffii, Vitis vinifera and Zea mays were used as query sequences to identify the GRF genes in the Phaseolus vulgaris genome.Ten identified predicted non-redundant GRF proteins were subjected to the Pfam and SMART domain searches, and the presence of the mandatory characteristic motifs (QLQ and WRC) of the GRF protein family was confirmed (Fig. 1).In addition, a zinc finger motif (1 His and 3 Cys residues) was also identified in WRC of PhvGRFs (Fig. 1).
The average intron and exon numbers in PhvGRF genes were 2.4 and 3.4, respectively, and all genes had 2 or more introns (Table 1, Fig. 3).Five common mo-tifs were observed according to the motif distribution analysis (Table 2, Fig. 4).Motifs 1 and 2 contain WRC and QLQ domains, respectively.These domains are conserved and mandatory GRF protein domains.FFD and TQL domains, which were detected in motif 4 and 5, are other GRF domains.The most common motif types, motif 1, motif 2 and motif 3, were observed in PhvGRF1, PhvGRF2, PhvGRF4, PhvGRF5, PhvGRF6, PhvGRF7, PhvGRF8 and PhvGRF10, whereas Ph-vGRF9 had only motif 1 and motif 2 (Fig. 4).
A phylogenetic tree was constructed using MEGA6 software based on the neighbor-joining method and the evolutional relationship of 10 PhvGRF proteins was determined (Fig. 5).PhvGRF proteins were classified into 2 main groups, named main group I and main group II, according to the phylogenetic tree.Main group I, which contained all genes except for PhvGRF9, was divided into two subgroups.PhvGRF6 and PhvGRF7 were found to be closest to each other, with a bootstrap value of 100%; PhvGRF8 was on a   branch close to them, with a 100% bootstrap value.PhvGRF2 and PhvGRF4, as well as PhvGRF5 and Ph-vGRF10, revealed 100% bootstrap values.PhvGRF9 was revealed as main group II by itself (Fig. 5).
Construction of phylogeny is very important for the functional evaluation of gene/protein families.For this reason, 10 predicted Phaseolus vulgaris GRFs (PhvGRF), 16 Zea mays GRFs (ZmGRF) and 9 Arabidopsis thaliana GRFs (AtGRF) were used for the construction of a phylogenetic tree with a bootstrapneighbor joining method (Fig. 6).According to the phylogenetic tree, most of the PhvGRF proteins were found closer to the Arabidopsis thaliana GRF proteins than to Zea mays GRF proteins (Fig. 6).
To understand the significance of the expression patterns of the PhvGRF genes during drought stress, the expression levels of 10 identified PhvGRF genes were determined in the root and leaf tissues of two different common bean varieties with different adaptation capacities against drought stress.According to the qRT-PCR results, changes in expression levels were observed in all common bean GRF genes except for PhvGRF4 and PhvGRF8 in response to drought stress (Fig. 7).Most of the PhvGRF genes displayed more consistent changes in expression in the root tissues of the Yakutiye and Zulbiye common bean cultivars than in leaf tissues.A decreasing trend in the expression of PhvGRF1, PhvGRF2, PhvGRF3, PhvGRF5, PhvGRF6, PhvGRF7 and PhvGRF10 genes was observed in the root tissues of Yakutiye after a 24-h drought treatment.In contrast to Yakutiye, an increasing expression trend was observed in the root tissues of Zulbiye for PhvGRF1, PhvGRF2, PhvGRF3, PhvGRF5, PhvGRF6, PhvGRF9 and PhvGRF10 genes.The expression pattern for the PhvGRF7 gene was only detected in the root and leaf tissues of the Yakutiye cultivar under drought stress, whereas there was no expression of the Ph-vGRF7 gene in the Zulbiye cultivar.At the same time, no significant expression was observed for PhvGRF4 and PhvGRF8 genes in either of the two cultivars.Hence the expression data of PhvGRF4 and PhvGRF8 genes were not included in Fig. 7.

DISCUSSION
In previous studies, GRF proteins of several plant species were investigated.A total of 24 GRF proteins were characterized in cucumber, melon and watermelon plant species [10].In the current study, we identified and characterized 10 PhvGRF genes in the Phaseolus vulgaris genome.When comparing the highly conserved domains (WRC, QLQ, TQL and FFD) of Phaseolus vulgaris with other previously studied plant species, we found similar domains in Brachypodium distachyon, Zea mays and the Cucurbitaceae family [7,10,11] (Table 2; Fig. 4).
It is well known from previous studies that GRF proteins are not well conserved because of the different exon-intron organization in the genomes of Brachypodium distachyon, Arabidopsis thaliana, Zea mays and Cucurbitaceae family [5,7,10,11].Similarly, different numbers of introns and different exon-intron organizations were observed for Phaseolus vulgaris GRFs in this study.
Gene duplications are important events that have significant roles in genomic and organismal evolution and that can originate from unequal crossing-over, retrotransposition or chromosomal duplication [24].The plasticity of an organism for adaptation to different environmental conditions would be limited without the occurrence of gene duplication events [24].Two segmental duplication events among PhvGRF6-PhvGRF7 and PhvGRF2-PhvGRF4 genes were observed in the current study.The schematic display of the conserved MEME motifs for the GRF proteins revealed that the PhvGRFs (PhvGRF6, PhvGRF7, PhvGRF2, PhvGRF4), which are encoded by these segmental duplicated genes, contain all the 5 motifs identified.At the same time, these genes were clustered together with the highest bootstrap value (100%) in subgroup I according to the phylogenetic analysis.In a previous study on the identification of GRF proteins in Brachypodium, one segmental duplication event was observed, while tandem duplication events, as described in the current study, were not previously observed [11].No segmental or tandem duplication events were observed in Cucurbitaceae family members in contrast to our findings [10].
When we evaluated the phylogenetic tree of Phaseolus vulgaris, Arabidopsis thaliana and Zea mays GRFs, it was observed that PhvGRFs are phylogenetically closer to AtGRFs with very high bootstrap values than to ZmGRFs.This could be due to the fact that Arabidopsis thaliana and Phaseolus vulgaris are dicot plant species, whereas Zea mays is a monocot species.This correlation was also found and explained in previous studies on the genome-wide identification of GRFs, which can be attributed to the orthology of GRF genes in monocot and dicot plant groups [10,11,14].
For the putative functional analysis of the PhvGRF proteins, 10 PhvGRFs were assigned the term "GO" using Blast2GO software.All PhvGRFs were annotated with the terms in the same GO groups as: regulation of transcription-DNA templated, developmental process under the biological process category, ATP binding under the molecular function category and nucleus under the cellular component category (Table 3).All data showed that PhvGRF proteins have similar roles in biological processes and molecular functions, and the same localization in plant cells.
Previous studies of Arabidopsis thaliana have shown that some GRF genes play a role in the regulation of the stress response mechanism [5].However, it is not known which common bean GRF proteins regulate abiotic stress responses such as drought stress.Hence, we analyzed the expression levels of 10 identified PhvGRF genes in the drought-tolerant Yakutiye and drought-sensitive Zulbiye common bean cultivars subjected to moderate drought stress for 24 h.Drought-tolerant and drought-sensitive cultivars were selected to make possible a comparative analysis of PhvGRFs genes in response to drought stress conditions in common bean.
According to the expression analysis of 10 Ph-vGRFs, inversely related expression patterns were observed in the roots of Yakutiye and Zulbiye cultivars based on their adaptation capacity against drought stress.PhvGRF1, PhvGRF2, PhvGRF3, PhvGRF5, Ph-vGRF6, PhvGRF9 and PhvGRF10 genes were upregulated 2-to 4-fold compared to the untreated controls after drought treatment in the root tissues of Zulbiye, which is a drought-sensitive common bean cultivar.Interestingly, the trend of PhvGRF1, PhvGRF2, Ph-vGRF3, PhvGRF5, PhvGRF6, PhvGRF7 and PhvGRF10 genes showed a consistent decline of 2-to 6-fold in the root tissues of the drought-tolerant Yakutiye cultivar subjected to 24 h drought stress.This significant correlation might be due to the involvement of PhvGRF genes in the stress defense mechanism against moderate drought stress in common bean.Some members of AtGRF genes were shown to be upregulated under normal conditions to suppress and regulate the expression of stress responsive genes, while some members were downregulated [4].Recent studies showed that GRFs act as a link between plant growth and defense signaling and stress responses [2,25].
The expression trend of PhvGRF genes in two common bean cultivars in response to drought stress has revealed a higher correlation in root tissues than in leaves.This difference could be explained by the higher sensitivity to drought stress in root tips, which are an actively growing tissue in plants [26].

CONCLUSION
Genome-wide in silico identification, characterization and expression analyses of Phaseolus vulgaris GRF proteins and genes offer an insight into their potential role in stress-related mechanisms.In addition, we have demonstrated that the expression of PhvGRFs was correlated with drought stress response in a cultivar-specific manner in common bean.The members of the GRF gene family may also be used for genetic engineering applications in common bean, which is an economically important legume crop.

Fig. 2 .
Fig. 2. Distribution of 10 PhvGRF genes in common bean chromosomes.(segmental duplications between the chromosomes are shown by the same symbol).

Fig. 4 .
Fig. 4. Schematic display of the conserved MEME motifs for the PhvGRF proteins.

Fig. 5 .
Fig. 5. Phylogenetic tree of Phaseolus vulgaris GRF proteins.The phylogenetic tree was constructed with MEGA6 using the neighbor-joining method.

Fig. 7 .
Fig. 7. Expression patterns of PhvGRF genes in the droughttolerant Yakutiye and drought-sensitive Zulbiye bean cultivars under moderate drought stress (C − control; ZR − Zulbiye cultivar roots treated with drought stress; YR − Yakutiye cultivar roots treated with drought stress; ZL − Zulbiye cultivar leaves treated with drought stress; YL − Yakutiye cultivar leaves treated with drought stress).

Table 1 .
Physiochemical, structural and sequence characteristics of GRF genes in Phaseolus vulgaris.

Table 2 .
The most conserved protein motifs in GRF protein sequences of Phaseolus vulgaris.Residues in bold show WRC, QLQ, FFD, and TQL domains, respectively.

Table 3 .
Blast2GO annotation details of Phaseolus vulgaris GRF protein sequences.