Genome-wide characterization of GRAS family genes in Medicago truncatula reveals their evolutionary dynamics and functional diversification

The GRAS gene family is a large plant-specific family of transcription factors that are involved in diverse processes during plant development. Medicago truncatula is an ideal model plant for genetic research in legumes, and specifically for studying nodulation, which is crucial for nitrogen fixation. In this study, 59 MtGRAS genes were identified and classified into eight distinct subgroups based on phylogenetic relationships. Motifs located in the C-termini were conserved across the subgroups, while motifs in the N-termini were subfamily specific. Gene duplication was the main evolutionary force for MtGRAS expansion, especially proliferation of the LISCL subgroup. Seventeen duplicated genes showed strong effects of purifying selection and diverse expression patterns, highlighting their functional importance and diversification after duplication. Thirty MtGRAS genes, including NSP1 and NSP2, were preferentially expressed in nodules, indicating possible roles in the process of nodulation. A transcriptome study, combined with gene expression analysis under different stress conditions, suggested potential functions of MtGRAS genes in various biological pathways and stress responses. Taken together, these comprehensive analyses provide basic information for understanding the potential functions of GRAS genes, and will facilitate further discovery of MtGRAS gene functions.


Introduction
Transcription factors play important roles in regulating various plant development and physiological processes. The plant-specific GRAS gene family has been studied in nearly 30 plant species from more than 20 genera [1,2]. Evolutionary analyses have suggested that the GRAS gene family possibly originated from bacteria through horizontal gene transfer [3]. The name PLOS  Genomic analyses of the GRAS gene family have been conducted in various species including Arabidopsis [1], rice [1], Chinese cabbage (Brassica rapa ssp. pekinensis) [9], tomato [12], castor beans [16] and grapevine (Vitis vinifera L.) [10], but have not been explored in legumes. In this study, we systematically and comprehensively analyzed GRAS genes in M. truncatula using comparative genomic strategies and experimental validation. The aims of this study were as follows: (1) identify and classify GRAS genes in M. truncatula; (2) explore the evolutionary dynamics of MtGRAS gene proliferation and uneven distribution; and (3) determine the functional diversity of MtGRAS genes by structure conservation analysis and expression profile analysis in different tissues and stress treatments. These findings provide insights into the molecular functions of MtGRAS genes, and will be helpful for future functional characterization of GRAS genes in legumes.

Identification of MtGRAS genes
The current genome sequence and annotation files (Mt4.0v1) of M. truncatula were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html). The most updated Hidden Markov Model (HMM) for the GRAS gene family, PF03514.11, was downloaded from the Pfam database (http://pfam.xfam.org). Using PF03514.11 as a query, we conducted a BLAST search against the entire protein dataset of M. truncatula with a cut-off E-value of 1e-10 using the blastall v2.2.18 package. Subsequently, all hit protein sequences were extracted using custom Perl scripts. Then, the integrity of the GRAS domain was evaluated using SMART tools [45], and candidate MtGRAS proteins composed of a truncated GRAS domain were identified. To obtain the gene structure, a GFF3 annotation file involving precise position information of introns and exons in each MtGRAS was retrieved from the genomic dataset, and uploaded to the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) [46]. Peptide length, molecular weight, and isoelectric point of each MtGRAS protein were calculated using the online ExPasy program (http://www.expasy.org/) [47].

Classification and conservation analysis of MtGRAS genes
The identified MtGRAS proteins were combined with the well-classified Arabidopsis and rice GRAS proteins and aligned using ClustalW [48]. Then, a phylogenetic tree was constructed in MEGA5 software using the neighbor-joining method with 1000 bootstrap replicates [49]. The phylogenetic tree was visualized using Evolview (http://www.evolgenius.info/evolview/) [50]. The MtGRAS genes were further categorized into different subgroups according to homology with GRAS genes in Arabidopsis and rice. The conserved motif analysis of MtGRAS was conducted using the motif finding tool, MEME (Multiple EM for Motif Elicitation, v4.10.0) with 20 motif numbers, and the order of motifs in each MtGRAS was evaluated by MAST [51]. The targets of miR171 were predicted in silico using the website (http://plantgrn.noble.org/ psRNATarget/). The Pv-miR171 genes were identified based on the homology searching stemloop sequence of osa-miR171, which obtained from the website of miRbase (http://www. mirbase.org/index.shtml). similarity of more than 70% and coverage greater than 75%. In addition, gene duplication information was also identified based on public data in the Plant Genome Duplication Database (PGDD, http://chibba.agtec.uga.edu/duplication/) [53]. If two homologous genes were separated by five or fewer genes, they were identified as tandem duplications (TD), while if two genes were separated by more than five genes or distributed in different chromosomes, they were referred to as segmental duplications (SD). To determine the evolutionary pressure acting on duplicated genes, Ka (non-synonymous substitution) and Ks (synonymous substitution) values were calculated using KaKs_Calculator 2.0 [54].
Expression profile analysis of MtGRAS genes using RNA-seq Illumina single read sequencing data for the transcriptome of M. truncatula were obtained from the NCBI Short Read Archive (accession numbers SRX099059-SRX099062). This dataset contained six different tissues including root, flower, bud, seedpod, blade, and nodule. We aligned all reads from each tissue to the reference genome of M. truncatula (Mt4.0v1) using tophat v2.1.0 [55]. Subsequently, the expression level for each MtGRAS was measured using Cufflink v2.1.1 [56], and the FPKM (fragment per kilobase per million mapped reads) representing the gene expression level of each MtGRAS was extracted with custom Perl scripts. A heatmap of the MtGRAS expression profile was created using Mev v4.8.1 [57].

Biotic and abiotic stress treatments
For hormone treatment, three-week-old M. truncatula (cv. Jemalong A17) seedlings were soaked in liquid MS medium with 30 μM gibberellin (GA3). For cold treatment, seedlings were grown in a greenhouse (12/12h photoperiod, 18-24˚C) and transferred to a cold chamber maintained at 4±1˚C. For salt treatment, 200 mM NaCl was sprayed on the leaves. Seedlings soaked in liquid MS medium without any treatment were used as control. Whole plants were sampled at 3h and 6h after treatment. For each treatment, six randomly chosen seedlings were pooled together to form a biological replicate. All plant samples were frozen in liquid nitrogen and stored at -80˚C until use.

Expression levels of MtGRAS genes under stress treatments
Total RNA was extracted from control and stress-treated samples using Trizol reagent (Invitrogen) based on the manufacturer's instructions. cDNA was synthesized using approximately 2 μg of RNA according to the manufacturer's protocol. Real-time quantitative PCR (qRT-PCR) was performed using SYBR Green mix (TaKaRa) on a LightCycler480 Real-Time PCR Detection System (Roche). The fold-change of expression was calculated with ACTIN as the internal reference gene. All the primers for qRT-PCR are listed in S1 Table.

Results
Identification and structural analysis of MtGRAS genes Using PF03154.11 as a query, we identified 59 GRAS genes in M. truncatula. Most of these genes contained a complete GRAS domain except five (MtGARS6, 23, 25, 26, and 49). To further elucidate the cause and consequences of MtGRAS gene expansion, we selected all the genes identified for further analysis, and named them based on their distribution and linear order on the respective chromosomes ( Table 1). The peptide length of MtGRAS varied greatly ranging from 69 amino acids (MtGRAS25) to 1,155 amino acids (MtGRAS55). Nearly 52 (88%) MtGRAS genes were intronless, which is consistent with most previous studies [1,12,58], while three members (MtGRAS25, 37, and 58) contained just one intron, and four members

Phylogenetic categories and conserved motif analysis of MtGRAS genes
To fully classify the MtGRAS gene family, 59 MtGRAS genes were analyzed with 32 and 53 GRAS genes in Arabidopsis and rice, respectively, to construct an unrooted phylogenetic tree using the neighbor-joining (NJ) method in MEGA5. Eight subfamilies were identified based on clade support values, the topology of the phylogenetic tree, and the previous classification of GRAS families in Arabidopsis and rice. We identified 9, 13,8,8,7,8,4, and 2 MtGRAS genes in the PAT1, LISCL, SHR, SCR, SCL3, HAM, DELLA, and LAS sub-branches, respectively (Fig 1).
To explore the potential biological functions of genes in each sub-branch, a detailed comparative analysis was performed (Fig 1). The PAT1 subfamily included nine MtGRAS genes including one member (MtGRAS52) and two members (MtGRAS7 and MtGRAS45) that have high sequence similarity to AtPAT1 and AtSCL13; MtGRAS11 was also in the PAT1 subfamily and was the closest homolog of CIGR1 in rice. The LISCL subfamily consists of 13 MtGRAS members. Four MtGRAS genes (MtGRAS36, 20, 19, and 18) share high homology with AtSCL9. In addition, another five members (MtGRAS34, 15, 14, 13, and 16) had high similarity with AtSCL14, demonstrating that they may function in stress-related processes [59,60]. The SCR and SHR subfamilies are crucial for stem cell maintenance that occurs during root and shoot development [27]. In our study, two homologous genes (MtGRAS47 and MtGRAS41) of SHR were identified, and one gene (MtGRAS58) shared high similarity with SCR. SCL3 regulates root cell elongation by integrating multiple signals in Arabidopsis [28]. Seven MtGRAS genes (MtGRAS26, 24, 23, 25, 38, 6 and 46) belonged to the SCL3 sub-branch, implying a similar function in root development. MtGRAS29, which belonged to the HAM subgroup, has been suggested to participate in nodule morphogenesis [2]. MtGRAS50, the closest paralog of MtGRAS29, may function in the same pathway as well. In the DELLA subfamily, four members (MtGRAS54, 63, 43, and 28) had the highest similarity to RGA and GAI. Two MtGRAS members (MtGRAS5 and 39) were identified as part of the LAS subgroup, which has several genes that have been found to regulate meristem formation [29][30][31].
Using multiple sequence alignment, the characteristic conserved domains located in the C-termini were identified including VHIID, LHRI, LHRII, PFYRE, and SAW (Fig 2 and S2-S5 Figs). We further explored conserved motifs in MtGRAS using MEME tools [51]. In total, 20 conserved motifs were found, and most of them had a similar distribution within the same subgroup (Fig 3). The logo of these motifs is listed in S6 Fig. The motifs located in the GRAS domain, including LHRI (motif6, motif9), VHIID (motif5, motf1), LHRII (motif13, motif7, motif10), PFYRE (motif4, motif11), and SAW (motif2, motif14, motif3), were shared across almost all MtGRAS members. In addition, motif8 was situated between LHRII and PFYRE and conserved among most MtGRAS subfamilies, suggesting its functional importance. Other motifs located outside the GRAS domain regions showed subgroup specific patterns. Motif16 was located between LHR1 and VHIID, and was specific to the PAT1 and LISCL subgroups, The GRAS gene family analysis in Medicago truncatula reveal their evolution and functions while motif17 was nested within LHRII and was LISCL specific. Motif12 and motif18 were located in the N-termini and were also only present in the LISCL subgroup (Fig 3).

Chromosomal distributions and duplication analysis of MtGRAS genes
Fifty-nine MtGRAS genes were mapped to the chromosomes of M. truncatula; MtGRAS1 was excluded because it was positioned on a scaffold ( Table 1). The distribution of MtGRAS genes among the chromosomes was uneven. Chr2 and chr4 are the "hot regions", and contained 15 (25.9%) and 13 (22.4%) MtGRAS genes, respectively; chr6 is the "cold region", and contained only one (1.7%) MtGRAS gene. Moreover, 5, 8, 7, 5, and 4 MtGRAS genes were found on chr1, chr3, chr5, chr7, and chr8, respectively (Fig 4). Based on these distributions, we explored duplication events of MtGRAS genes. Seventeen duplicated MtGRAS gene pairs were identified. Most duplication events occurred in chr2 (n = 7), chr3 (n = 6), and chr4 (n = 5); one duplication event occurred in chr1, chr5, and chr8; and no duplication events were identified in chr6 and chr7 (Fig 4). We also verified the types of duplication. The results suggested that The GRAS gene family analysis in Medicago truncatula reveal their evolution and functions 11 duplicated gene pairs arose from tandem duplications, while six pairs were segmental duplications (Table 2).
To understand the evolutionary process of gene duplications, we evaluated the positions of duplicated genes. Seven duplicated MtGRAS genes (MtGRAS9, 14, 15, 13, 19, 24, 2, and 49) were positioned near the telomeres of each chromosome, and three genes (MtGRAS24, 29, and 38) were located around the centromeres, implying that the highly repetitive nature of these regions may lead to gene duplication (Fig 4). To reveal the evolutionary dynamics of duplicated MtGRAS genes, we calculated nonsynonymous substitution rates (Ka) and synonymous substitution rates (Ks) between duplicated genes. All of the duplicated gene pairs have Ka/Ks values <1, suggesting that purifying selection acted on them ( Table 2). In our data, expression of most MtGRAS genes was identified in at least in one tissue. Expression of nine MtGRAS genes (MtGRAS4, 10, 23, 25, 26, 27, 48, 49 and 54) was not detected in our transcriptome data, which may be the result of spatial and temporal expression patterns or uncharacterized pseudogenes. The FPKM values of each MtGRAS are shown in S2 Table and the expression profiles were clustered across six tissues (Fig 5). Among the 50 expressed MtGRAS genes, 38 were highly expressed (FPKM >1), while the other 12 had low  Table). The high proportion of MtGRAS genes expressed in nodules indicates that additional MtGRAS genes besides NSP1 and NSP2 participate in the process of nodulation. Furthermore, three Arabidopsis GRAS genes (SCL6, 22, and 27) in the HAM subfamily are post-transcriptionally regulated by miR171. Here, the two closest homologs of SCL6, MtGRAS1 and MtGRAS31, were identified as having a putative binding site for miR171 (Fig 6A). Two Pv-miR171 genes were detected in the Medicago genome. The expression pattern of Pv-miR171 genes were negatively correlated with their targets. Both MtGRAS1 and 31 exhibited highest expression in nodules, while Pv-miR171 genes (Pv-miR171-1 and Pv-miR171-2) showed the lowest expression in this tissue. In buds, MtGRAS1 and 31 accumulated the least transcript, but the transcript of Pv-miR171 genes, especially Pv-miR171-1, was the highest among different tissues (Fig 6B).
Since duplicated genes can exhibit significant variation in gene expression, we next explored the expression divergence of 17 sets of MtGRAS duplicated genes. Detailed expression information of these genes is listed in S3 Table. Eleven duplicated gene pairs shared similar expression patterns with different expression levels. For example, both MtGRAS47 and its duplicated gene MtGRAS41, had higher expression in bud and nodule, but showed lower expression in blade and flower (Fig 7A). This pattern was also observed in gene pairs MtGRAS2/ 30, 19/18, 19/20, 29/50, 9/8, 9/ (Fig 7B and S7 Fig). For example, MtGRAS61 was expressed in six tissues, and had the highest expression in flower, but the transcript level of the duplicated gene MtGRAS49 was not detected (Fig 7B). Intriguingly, two duplicate gene pairs, MtGRAS32/33 and 38/46 -showed different expression profiles (Fig 7C  and 7D). In the gene pair, MtGRAS32/33, MtGRAS32 had higher expression in nodule and blade, whereas MtGRAS33 had higher expression in nodule and flower. MtGRAS38 was highly expressed in bud, whereas the transcripts of its homologous gene MtGRAS46 were enriched in nodules. In all, the expression differences of duplicated genes presented here implies that genes may have various evolutionary outcomes after duplication.

Responses of MtGRAS genes to different stress treatments
We further examined changes in transcript abundance in response to stress treatments including GA3, salt, and cold using qPCR. Twelve MtGRAS genes were used to assess responses to The GRAS gene family analysis in Medicago truncatula reveal their evolution and functions treatments, including: MtGRAS32 and 35 (LISCL subfamily); 60 and 47 (SHR subfamily); 45 and 51 (PAT1 subfamily); 50 and 61 (HAM subfamily); 37 (SCR subfamily); 39 (LAS subfamily); and 38 and 46 (SCL3 subfamily). All of these genes exhibited differential expression in response to at least one stress treatment (Fig 8). After GA3 treatment, the expression levels of most MtGRAS genes (10) (MtGRAS32 and 35). The expression levels of six genes (MtGRAS60, 47, 50, 61, 39, and 38) decreased at 3h but increased at 6h, while three genes (MtGRAS51, MtGRAS37, and MtGRAS46) increased in linear order and reached the highest expression level at 6h. Finally, the expression of MtGRAS45 was highest at 3h but decreased at 6h (Fig 8). The GRAS gene family analysis in Medicago truncatula reveal their evolution and functions

Discussion
GRAS transcription factors play essential roles in regulating plant growth and development. However, the prevalence and functional diversity of the GRAS family in M. truncatula have not been thoroughly investigated. In this study, we performed a comprehensive analysis of the GRAS gene family in M. truncatula. The features of MtGRAS genes, including their chromosomal distribution, phylogenetic classification, expression profiles, and responses to various stresses were explored. Results of these analyses allowed us to study the evolution of the GRAS family and draw hypotheses about the potential functions of unidentified genes. Our results demonstrated that duplication was the evolutionary force behind GRAS gene family expansion. First, the number of GRAS was different among species. In this study, we identified 59 GRAS genes in M. truncatula, which is lower than the number in Populus (106) [15], but higher than in other species such as Arabidopsis (33) [1], rice (57) [1], Chinese cabbage (46) [9], tomato (53) [12] and grapevine (52) [10]. Gene duplication might cause these differences in numbers of GRAS family members. Second, MtGRAS genes were unevenly distributed among the chromosomes, with the "hot regions" on chr2 (15 members) and chr4 (13 members), and the "cold region" on chr6 (1 member). Interestingly, more duplication events were found in the hot regions (7/chr2, 4/chr4) (Fig 4). We further related the duplication events to chromosome positions, and found that duplicated MtGRAS gene pairs tend to be located in regions with low levels of conservation (10 members), such as peritelomeres and pericentromeres; tandem and segmental duplication were enriched in these regions because of the composition of repetitive elements (Fig 4). Furthermore, among the 17 duplicated MtGRAS gene pairs, seven belonged to the LISCL subfamily and clustered in chr2 and chr4, suggesting that gene duplication might cause proliferation of the LISCL subfamily. Previous studies have found that duplication is common in the GRAS gene family. For example, 2/ 34, 15/53, 17/60, and 40/106 GRAS genes were identified as duplicated genes in Arabidopsis [1], tomato [12], rice [1], and Populus [15], respectively, which further validates the contribution of duplication to expansion of the GRAS gene family. Third, nearly 88% of MtGRAS genes were intronless, consistent with other species such as tomato (77.4%) [12], Prunus mume (82.2%) [58], and Arabidopsis (67.6%) [1]. Intronless genes have also been discovered in other large gene families, such as DEAD box RNA helicase [61] and F-box transcription factors [62]. Rapid duplication after horizontal gene transfer from bacteria is the main contributor to the high proportion of intronless genes.
To our knowledge, transcription factors belonging to the same taxonomic clade exhibit recent evolutionary origins and specific conserved motifs associated with functional specification. Because of this, a comparison of homologous genes in the MtGRAS family, including The GRAS gene family analysis in Medicago truncatula reveal their evolution and functions protein sequences and expression profiles, would be an effective method to predict the function of uncharacterized genes. In this study, 50 MtGRAS genes were expressed in at least one tissue according to our transcription analysis, and the expression patterns varied across a variety of tissues, as previously reported in Populus, Prunus mume, tomato, and grapevine. Noticeably, more than half (30/59) of the MtGRAS genes had the highest expression in the nodule, and 5, 7, 4, and 4 members were preferentially expressed in blade, bud, flower, and root, respectively. This result suggested that the functions of MtGRAS have dramatically diverged. Genes belonging to the LAS subfamily have been found to participate in regulating axillary meristem development. For example, the mutation of MOC1 resulted in the phenotype of no tillers except for a main culm in rice [29]. In an Arabidopsis knockout, the homolog of MOC1, named AtLAS, led to an inability to form lateral shoots during vegetative development [31]. In M. truncatula, two MtGRAS genes (MtGRAS5 and MtGRAS39) belonged to the LAS subgroup. Interestingly, both of them were expressed highest in bud, which indicated that they might also play a vital role in axillary meristem formation (Fig 5 and S2 Table). The DLT gene in rice and its orthologs, AtSCL28 in Arabidopsis and VviGRAS8a in Vinus vinfera, modulate the expression of a brassinosteroid-related gene [59,63]. In tomato, the ortholog of DLT was validated to be involved in the flower-fruit transition by mediating brassinosteroid signaling [12]. In our study, MtGRAS12 was the homolog of DLT, and was preferentially highly expressed in bud, but had low expression in flower (Fig 5 and S2 Table). The result suggested that MtGRAS12 might function in response to brassinosteroid signaling during bud development. DELLA genes participate in various developmental processes including flower development, stem elongation, and seed germination [28]. In addition, the DELLA proteins also participate in hormone signaling pathways, such as the gibberellin, jasmonate, auxin, brassinosteroid, and ethylene pathways [64]. In our data, the closest homolog of RGA and GAI, MtGRAS28, was highly expressed in different tissues (FPKM>40) including root, seedpod, and blade, supporting a role in diverse developmental processes (Fig 5 and S2 Table).
Generally, the evolutionary fate of duplicated genes includes nonfunctionalization, neofunctionalization, or subfunctionalization [65,66]. We further evaluated the evolutionary dynamics and consequences of duplicated MtGRAS genes. All of the duplicated genes were under purifying selection (Ka/Ks <1), implying that these genes were still strongly controlled after duplication (Table 2). We next examined the divergence of expression in 17 sets of duplicated MtGRAS gene pairs. Eleven duplicated genes showed similar expression patterns to the original gene, but with different expression levels; four duplicated genes were not expressed in our transcription dataset (Fig 7 and S3 Fig). Furthermore, the duplicated gene pairs MtGRAS38/46 and MtGRAS32/33 exhibited different expression patterns, suggesting that novel functions might evolve after duplication (Fig 7). Further efforts need to be made to elucidate the functional diversity of duplicated genes.
Previous studies have demonstrated that the GRAS protein could interact with ERN to regulate gene expression during rhizobial infection [67]. Two MtGRAS genes, NSP1 and NSP2, belonging to SHR and HAM, respectively, were associated with enhancing Nod factor elicitation [2]. In our study, more than half (30/50) of the MtGRAS genes had the highest expression in nodules, which implied that MtGRAS genes other than NSP1 and 2 might participate in the process of nodulation. In Arabidopsis, three GRAS members in the HAM subfamily were posttranscriptionally regulated by miR171 (AtSCL6, 22, and 27). Interestingly, in the present study, the two closest homologs of AtSCL6 (MtGRAS1 and 31) were found to have a putative binding site for miR171. Both of these genes were highly expressed in nodules, especially the MtGRAS1 gene, which had an FPKM value higher than 40 (Fig 6). These results indicated that mi172 might also be involved in nodule development of legumes [37][38][39].
Numerous studies have found that transcription factors in the GRAS family could be influenced by various biotic and abiotic stresses. Gibberellin, auxin, brassinosteroid, abscisic acid, ethylene, and salicylic acid also play important roles in a diverse array of developmental processes including germination, flowering time, and stem elongation [64]. GAI, RGA, and RGL in the DELLA subfamily were repressors of gibberellin signaling [4][5][6]. Loss of function in the Arabidopsis mutants scr and shr, resulted in hypersensitization to abscisic acid [6]. BnSCL1 in Brassica napus showed differential dose responses to auxin in shoots and roots [68]. Recently, a study in tomato demonstrated that the expression level of GRAS genes could be modulated by signaling of multiple phytohormones, including gibberellin, auxin, brassinosteroid, ethylene, and salicylic acid [12]. Additionally, several studies have revealed the participation of GRAS genes in response to abiotic stresses, such as cold, drought, salt, and heat. In Arabidopsis, over-expression of a poplar GRAS gene, PtSCL7, enhanced tolerance to salt and drought stress [69]. SCL14 in Arabidopsis was involved in the activation of a broad-spectrum detoxification network, and its ortholog in rice, OsGRAS23, was involved in regulating the drought stress response [23,59,60]. The gene BoGRAS was significantly upregulated during heat stress in Brassica oleracea [70]. In our study, 12 MtGRAS genes from different subgroups were randomly selected to explore their responses to biotic (GA3) and abiotic (NaCl and 4˚C) stresses. We found that nearly all MtGRAS genes could be affected by different stress treatments (Fig 8).
Most MtGRAS genes (10/12) were downregulated after treatment with GA3, while only two genes were upregulated, implying that most MtGRAS genes had negative roles in response to this hormone. Under salt treatment, the expression levels of seven genes decreased, and the expression levels of three genes increased, suggesting that MtGRAS genes modulate the signaling of response to salt through complicated networks. The majority of genes (10/12) increased their expression level under the 4˚C treatment, indicating they might positively regulate the response gene in the cold condition. In addition, most MtGRAS genes could be influenced by both hormone and abiotic stress treatments, indicating the coordinated response of these two environmental determinants.

Conclusions
In this study, we performed a genome-wide analysis of the GRAS gene family in M. truncatula based on publicly available genome data. Fifty-nine MtGRAS genes were identified and categorized into eight subfamilies by phylogenetic analysis. Conserved motif analysis combined with expression profile measurement in different tissues and environmental treatments demonstrated the functional conservation and diversity of MtGRAS genes. The evolutionary dynamics of MtGRAS family members was further inferred by analyzing the cause and consequence of duplicated MtGRAS gene pairs. We foresee that these results will be of great value for further functional characterization of the MtGRAS gene family and for genetic improvements in agronomic traits or stress tolerance in legumes.