Genome-Wide Analysis of the AP2/ERF Superfamily Genes and their Responses to Abiotic Stress in Medicago truncatula

The AP2/ERF superfamily is a large, plant-specific transcription factor family that is involved in many important processes, including plant growth, development, and stress responses. Using Medicago truncatula genome information, we identified and characterized 123 putative AP2/ERF genes, which were named as MtERF1–123. These genes were classified into four families based on phylogenetic analysis, which is consistent with the results of other plant species. MtERF genes are distributed throughout all chromosomes but are clustered on various chromosomes due to genomic tandem and segmental duplication. Using transcriptome, high-throughput sequencing data, and qRT-PCR analysis, we assessed the expression patterns of the MtERF genes in tissues during development and under abiotic stresses. In total, 87 MtERF genes were expressed in plant tissues, most of which were expressed in specific tissues during development or under specific abiotic stress treatments. These results support the notion that MtERF genes are involved in developmental regulation and environmental responses in M. truncatula. Furthermore, a cluster of DREB subfamily members on chromosome 6 was induced by both cold and freezing stress, representing a positive gene regulatory response under low temperature stress, which suggests that these genes might contribute to freezing tolerance to M. truncatula. In summary, our genome-wide characterization, evolutionary analysis, and expression pattern analysis of MtERF genes in M. truncatula provides valuable information for characterizing the molecular functions of these genes and utilizing them to improve stress tolerance in plants.

The AP2/ERF superfamily is a large, plant-specific transcription factor family that is involved in many important processes, including plant growth, development, and stress responses. Using Medicago truncatula genome information, we identified and characterized 123 putative AP2/ERF genes, which were named as MtERF1-123. These genes were classified into four families based on phylogenetic analysis, which is consistent with the results of other plant species. MtERF genes are distributed throughout all chromosomes but are clustered on various chromosomes due to genomic tandem and segmental duplication. Using transcriptome, high-throughput sequencing data, and qRT-PCR analysis, we assessed the expression patterns of the MtERF genes in tissues during development and under abiotic stresses. In total, 87 MtERF genes were expressed in plant tissues, most of which were expressed in specific tissues during development or under specific abiotic stress treatments. These results support the notion that MtERF genes are involved in developmental regulation and environmental responses in M. truncatula. Furthermore, a cluster of DREB subfamily members on chromosome 6 was induced by both cold and freezing stress, representing a positive gene regulatory response under low temperature stress, which suggests that these genes might contribute to freezing tolerance to M. truncatula. In summary, our genome-wide characterization, evolutionary analysis, and expression pattern analysis of MtERF genes in M. truncatula provides valuable information for characterizing the molecular functions of these genes and utilizing them to improve stress tolerance in plants.

INTRODUCTION
On a global scale, plant growth, and development are threatened by various abiotic stresses, such as extreme temperatures, drought, and high salinity. Due to environmental conditions continuously accumulating, plants are being confronted with increasingly serious challenges to their survival. Plants employ complex regulatory mechanisms to adapt to environmental stresses, undergoing physiological and biochemical changes in response to unfavorable conditions (Zhu, 2002;Chinnusamy et al., 2004;Katagiri, 2004). It is notable that plants contain numerous genes encoding transcription factors (TFs), which regulate the expression of downstream genes by binding to their cis-acting elements Yamaguchi-Shinozaki and Shinozaki, 2006;Le Hir and Bellini, 2013). TFs play important roles in plant growth, development and responses to environmental stress by directly responding to stress or regulating the expression of downstream target genes, confirming their critical roles in plant life cycles. However, only a few TF families have been characterized outside of wellstudied model plant systems, such as rice and Arabidopsis. The APETALA2/ethylene-responsive element binding factor (AP2/ERF) superfamily is one of the largest groups of TFs in plants. These TFs contain at least one AP2 domain. Based on the number of AP2 domains and other DNA binding domains, AP2/ERF TFs are classified into four families, including the AP2, ERF, RAV, and Soloist families (Cao et al., 2001;Sakuma et al., 2002;Mizoi et al., 2012). AP2 family members contain a double, tandemly repeated AP2 domain, while ERF family members contain a single AP2 domain. RAV family members have a single AP2 domain and an additional B3 domain, i.e., a DNA-binding domain commonly found in other TFs (Nakano et al., 2006;Mizoi et al., 2012;Li et al., 2015).
AP2/ERF TFs regulate a number of biological processes, such as plant growth, development, and responses to stress (Mizoi et al., 2012;Matías-Hernández et al., 2014). AP2 family members regulate the development of various plant tissues. For example, CRL5 is expressed in the stem region and regulates the crown root initiation process in rice (Kitomi et al., 2011). Some members of the AP2 family regulate reproductive development in Arabidopsis, including flower, ovule, and sepal development (Kunst et al., 1989;Krizek, 2009). AP2 TFs are also involved in fruit development and the ripening process in tomato and grapevine (Chung et al., 2010;Licausi et al., 2010;Sharma et al., 2010). The ERF family consists of two subfamilies (ERF and DREB), which participate in many developmental and stress response processes. ERF subfamily members, which bind to GCC-boxes, are involved in hormone signaling pathways, such as the ethylene, jasmonic acid, and salicylic acid pathways, which are important for plant development and stress responses (Fujimoto et al., 2000;Oñate-Sánchez and Singh, 2002;Andriankaja et al., 2007;Mantiri et al., 2008a). DREB subfamily members bind to dehydration-responsive element/C-repeat (DRE/CRT) elements, which are present in stress-responsive genes, including RD (responsive to desiccation, RD29) and COR (cold-regulated, COR15) genes (Cao et al., 2001;Xu et al., 2011). Therefore, many members of the DREB subfamily improve the stress tolerance of various plants under different environmental stress, including cold, drought, salinity stress, and so on.
Since the release of the whole-genome sequences of many plant species, the AP2/ERF superfamily has been successfully identified and investigated in plants including Arabidopsis, rice (Nakano et al., 2006), grapevine (Licausi et al., 2010), poplar (Zhuang et al., 2008), and soybean (Zhang et al., 2008). Genome-wide analysis of AP2/ERF TFs has helped elucidate their regulatory functions in plant growth, development and especially stress responses. Medicago truncatula is an annual, diploid legume plant. Due to its features, such as its small genome, self-pollination, high genetic transformation efficiency, nitrogen fixation, and so on, this plant has been chosen as a model legume for molecular genetic and genomic analyses. To date, only 37 AP2/ERF genes were identified based on EST sequences, which is far fewer than the number identified in any other plant species (Zhang et al., 2013b). The functions of AP2/ERF superfamily members in M. truncatula, especially those in responses to abiotic stress, have hardly been reported. ERN and ERN1, with an AP2 domain, were isolated from M. truncatula. Many studies have demonstrated that they have a regulation function in the process of nodulation, as do ERN2 and ERN3 (Middleton et al., 2007;Vernié et al., 2008;Hirsch et al., 2009). MtSERF1, regulated by the plant hormones auxin and cytokinin, plays an important role in somatic embryogenesis (Mantiri et al., 2008a,b). MtERF1-1, a member of the AP2/ERF B3 subgroup, has been shown to mediate resistance to root pathogens in M. truncatula (Anderson et al., 2010). WXP1 and WXP2 can improve drought tolerance in transgenic alfalfa by increasing leaf wax accumulation (Zhang et al., 2005(Zhang et al., , 2007. Pennycooke et al. (2008) isolated MtCBF1-3 from M. truncatula, and demonstrated their action against low temperature stress. MtCBF4 has been shown to respond to abiotic stresses, including cold, drought, salt, and ABA . Interestingly, 12 potential MtCBF genes clustering on chromosome 6 have been shown to play a major role in tolerance to freezing by QTL mapping (Tayeh et al., 2013). However, their mechanism of doing so, has not yet been discovered. Their expression profiles during freezing remain to be determined. Recently, a draft of the M. truncatula genome sequence was completed and released (Young et al., 2011). A number of gene families have been analyzed based on this genome information, such as ARF and CCCH families, which were promoting Medicago Genus and legume genetic research (Zhang et al., 2013a;Gujaria-Verma et al., 2014;Shen et al., 2015).
In this study, we performed a comprehensive analysis of the AP2/ERF superfamily in M. truncatula, including phylogenetic analysis, chromosomal localization, gene duplication analysis, and expression profiling. We also characterized the functions of these TFs in the abiotic stress response via transcriptome analysis. The results of this study will be helpful for future investigations aimed at the functional characterization of these AP2/ERF TFs and their utilization for the genetic improvement of legumes.

MATERIALS AND METHODS
Identification and Classification of the AP2/ERF Genes in Medicago truncatula M. truncatula genome and proteins sequences were downloaded from the JCVI website (M. truncatula Genome Project v4.0, http://www.jcvi.org/medicago/; Young et al., 2011), and Arabidopsis AP2/ERF gene sequences were downloaded from the DATF database (http://datf.cbi. pku.edu.cn; Guo et al., 2005). These Arabidopsis AP2/ERF sequences were utilized for BLAST (Altschul et al., 1990) searches against the M. truncatula genome sequence with the parameters of expected values ≤1E-3 and more than 80% coverage. All BLAST hits were retrieved and searched using the Hidden Markov Model (HMM) profile of the AP2 domain (PF002701), which was downloaded from the Pfam website (pfam.sanger.ac.uk;Finn et al., 2014). The AP2/ERF sequences were confirmed based on the presence of an AP2 domain, and all of the putative AP2/ERF proteins were aligned to Arabidopsis AP2/ERF proteins to classify them into different groups, as described by Nakano et al. (2006). Furthermore, all of the annotation information about putative AP2/ERF genes was retrieved from the M. truncatula genome website, and the number and distribution of introns in AP2/ERF genes were investigated using M. truncatula genome annotation information.
Phylogenetic and Conserved Motif Analysis of the AP2/ERF Genes Multiple alignments of candidate AP2/ERF protein sequences were carried out using ClustalW with default parameters (Thompson et al., 2002). Unrooted phylogenetic trees of all AP2/ERF proteins were generated with MEGA (V4.0) using the neighbor-joining (NJ) method with the following parameters: Poisson correction, pair-wise deletion and 1000 bootstrap replicates (Tamura et al., 2007). Conserved motifs in M. truncatula AP2/ERF TFs were identified using the motif finding tool MEME (Multiple EM for Motif Elicitation, V4.8.1; Bailey et al., 2006). MEME searching was performed across MtERF protein sequencing using the following parameters: (1) optimum motif width was set to ≥10 and ≤200; (2) the maximum number of motifs was set to identify 25 motifs; (3) occurrences of a single motif distributed among the sequences with model: zero or one per sequence (-mod zoops).

Chromosomal Localization and Gene Duplication Analysis of the AP2/ERF Genes
Positional information about all of the AP2/ERF genes was investigated, and diagrams of their chromosome locations in M. truncatula were drawn using the Circos software (http://circos. ca/; Krzywinski et al., 2009), revealing duplications between AP2/ERF genes in M. truncatula. If two genes with similarities of more than 85% were separated by four or fewer gene loci, they were identified as tandem duplications (TD). Others were identified as segmental duplications (SD), separated by more than five genes. In addition, duplications between the AP2/ERF genes were also identified and complemented using the PGDD database (http://chibba.agtec.uga.edu/duplication/; Lee et al., 2013). Duplicated genes between different chromosomes or loci were linked with colored lines in the diagrams using the Circos as previously described.

In silico Expression Analysis of the AP2/ERF Genes during Plant Development
Genome-wide transcriptome data from M. truncatula in different tissues during development were downloaded from the NCBI short read archive database (SRA database; http://www.ncbi. nlm.nih.gov, Accession numbers SRX099057-SRX099062). The transcriptome data were derived from six tissues, including roots, nodules, blades, buds, seedpods, and flowers. All transcriptome data were mapped to the M. truncatula genome using the TopHat (Trapnell et al., 2009), and the expression of MtERF genes was evaluated using the software Cufflinks as previously described (Trapnell et al., 2012). The expression data were analyzed and clustered using hierarchical cluster programs HCLUST of R and CLUSTERGRAM of Matlab (MathWorks, R2012a).

Plant Material and Stress Treatments
Seeds of M. truncatula (cv. Jemalong A17) were germinated and transferred to a mixture of perlite and sand (3:1, V/V). All seedlings were grown in a growth chamber (Conviron E15) at a temperature of 18 (night) and 24 • C (day), a relative humidity of 60-80% and a 14/10 h photoperiod (daytime, 06:00-20:00). The seedlings were irrigated with half-strength Hoagland solution once every 2 days, and after 8 weeks, they were randomly divided into six groups for stress treatments. For cold stress (B group) and freezing stress (C group) treatment, the seedlings were transferred into another chamber with the temperature set at 4 or −8 • C, respectively. For drought stress (D group) and salt stress (E group) treatment, the seedlings were treated with 300 mM mannitol or 200 mM NaCl solution, respectively. For ABA treatment (F group), the seedlings' leaves were sprayed with 100 µM ABA solution. Control (untreated, A group) and treated (B-F groups) seedlings were harvested at 3 h after treatment. For each group, five randomly chosen whole seedlings were pooled to form a biological replicate. All plant samples were frozen in liquid nitrogen and stored at −80 • C until use.

Transcriptome Analysis of the Response of the AP2/ERF Genes to Abiotic Stress
Total RNA was extracted from six samples (one biological replicate sample per group) using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA) following the manufacturer's instructions. The integrity of the RNA was assessed by formaldehyde agarose gel electrophoresis. Total RNA was quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and a Bioanalyzer 2100 (Agilent Technologies, CA). RNA Integrity Number (RIN) values were greater than 8.0 for all samples. Purified RNA samples were sent to BGI-Shenzhen Ltd. (Shenzhen, China) for construction of pair-end cDNA libraries and Illumina sequencing of abiotic stress-treated samples. Processing of raw data, removal of adapter sequences, base-calling, and quality value calculations were performed to produce clean data. Clean reads from six samples were mapped to the M. truncatula genome, and splice junctions were mapped using the TopHat. MtERF gene expression across six treatment samples (groups A-F) were evaluated using the Cufflinks software, and clustered using hierarchical cluster programs HCLUST of R and CLUSTERGRAM of Matlab (Mathworks, R2012a). Compared to controls, MtERF genes with fold changes ≥2 or ≤0.5 were identified as differentially expressed in response to abiotic stresses. The expression of MtERF genes was analyzed and visualized using the package ggplot2 of R platform.   and 72 • C for 1 min, and the experiments were repeated three biological replicates. The C T method was used to calculate relative expression levels of MtERF genes using GAPDH as reference gene. Primers of nine MtERF genes (randomly selected from DREB subfamily) and GAPDH gene used for qRT-PCR detection are listed in Table S1.

Identification of the AP2/ERF TFs in Medicago truncatula
Using homology searches and domain confirmation, we identified 123 putative AP2/ERF TF genes in M. truncatula, designated as MtERF001 to MtERF123, and we determined that these AP2/ERF genes encode putative proteins ranging from 120 to 689 aa in length (see Table 1). Among these TF genes, 98 genes with a single AP2/ERF domain were assigned to the ERF family, and based on the similarity of their encoded amino acid sequences, these genes were further classified into two subfamilies: 50 genes were identified as DREB subfamily members and 48 genes were identified as ERF subfamily members. Of the remaining MtERF genes, 21 genes were grouped into the AP2 family due to their tandemly repeated double AP2/ERF domain, and three genes were classified as RAV family members, as they encode proteins containing a single AP2/ERF domain together with a B3 domain. MtERF123 is homologous to the Arabidopsis Soloist gene (At4g13040) and was therefore designated as Soloist, as shown in Figure 1. As Nakano et al. (2006) previously described, the genes of the ERF family can be subdivided into 10 groups according to their similarity to Arabidopsis ERF sequences, as shown in Figure 2. The DREB subfamily includes group I-IV, containing five, 14, 26, and five members, respectively, while the ERF subfamily consisted of group V-X, with eight, nine, three, six, 21, and two members, respectively. By BLAST search, nine AP2/ERF genes previously identified and totally characterized were also confirmed in the present study, see Table 2. The AP2/ERF family in M. truncatula has relatively few members compared with FIGURE 2 | Phylogenetic tree analysis of the ERF family in Medicago truncatula.  (147), soybean (148), grapevine (149), rice (180), and poplar (202) (Nakano et al., 2006;Zhang et al., 2008;Zhuang et al., 2008;Licausi et al., 2010). The numbers of AP2 and RAV family members varies little among species, ranging from 18 to 29 and from three to six, respectively, which contributes little to the reduction in the total number identified in M. truncatula. This small number may arise from the reduced number of ERF family members, as there are 98 members of the ERF family in M. truncatula and 122, 122, 169, and 145 members in Arabidopsis, grapevine, soybean, rice, and poplar, respectively.

Phylogenetic and Conserved Motif Analysis of the AP2/ERF TFs in Medicago truncatula
To determine the evolutionary relationships between AP2/ERF family proteins in M. truncatula, a phylogenetic tree was constructed based on alignment of full-length sequences of MtERF proteins. The phylogenetic tree confirms that MtERF TFs could be classified into four groups, shown as Figure 1, which is consistent with the classification results obtained using homology searches as described above. The phylogenetic relationships of the 98 ERF family members were assessed in depth. The results show that these genes (except MtERF21) could be divided into 10 groups, as described by Nakano et al. (2006). As shown in Figure 2, groups I-IV were identified as DREB subfamily members, while groups V-X were characterized as ERF subfamily members. The conserved motifs in AP2/ERF family proteins in M. truncatula were investigated using MEME, revealing a total of 25 conserved motifs (designated motifs 1-25), as shown in Figures 3, 4. Motifs 1-6 were found to be similar to the AP2/ERF domain region, while the remaining motifs corresponded to regions outside of the AP2/ERF domain region, which are distributed in specific clades in the phylogenetic tree. Proteins in the same group or subgroup contain similar motifs, while the motifs are divergent among different groups or subgroups. For example, motif 23 is only present in each member of the RAV family, which indicates that this motif is specific to the RAV family. Similarly, motif 19 is shared by members of the AP2 family, motif 16 is present in groups ERF-I and II, and the ERF II group contains motif 14 while the ERF I group does not. Motifs 8-10 are specific to the ERF III group, and motif 8 is shared by each member of the ERF III group. Finally, while motif 9 and 10 are shared by most members (16/26, 62%), motif 11 and 13 are conserved in the ERF V group and motif 12 is shared by groups ERF VI and VII. These results indicate that most motifs are distributed among specific groups, which is correlated with their functional divergence (see Figures S1-S4).

Chromosomal Locations and Duplication of the AP2/ERF TFs in Medicago truncatula
The 123 AP2/ERF TF genes are distributed throughout the eight chromosomes of M. truncatula; their physical locations on chromosomes are shown in Figure 5. Each M. truncatula chromosome contains some AP2/ERF genes in numbers ranging from 10 to 27. Chromosomes 1 and 5 have the highest number of AP2/ERF genes (27 and 21 genes, respectively), while chromosomes 7 and 8 contain the fewest genes (11 and 10, respectively). The AP2/ERF genes are not randomly distributed on each chromosome, as there are some gene  clusters "hot regions" on the chromosomes. For example, chromosome 6 contains 10 AP2/ERF genes (MtERF34-43) in a short chromosome region (∼390 kb), and chromosomes 1, 2, and 5 contain similar gene clusters, as shown in Figure 5. In addition, using gene duplication analysis, we identified 38 pairs of gene duplications, which arose from tandem duplications and segmental duplications. Tandem duplications produced MtERF gene clusters or hot regions, such as the MtERF34-43 cluster on chromosome 6 and the MtERF10-12 cluster on chromosome 2. Segmental duplication produced many homologous AP2/ERF genes on different chromosomes, which expanded the numbers of MtERF genes from different groups. For example, MtERF52, 54, 56, and 57 from group ERF-V are distributed on different chromosomes (MtERF52 on chromosome 2, MtERF54 on chromosome 3, MtERF56, and 57 on chromosome 4), which are products of genome segmental duplication.

Expression Profiles of MtERF Genes in M. truncatula Tissues
We investigated the expression profiles of MtERF genes in various tissues using high-throughput sequencing data from NCBI, including root, nodule, blade, bud, seedpod, and flower tissues, revealing that 75 MtERF genes were expressed in at least one of six tissues. Of these, the expression of 46 genes was detected in root tissue, 46 genes in nodule tissue, 42 genes in blade tissue, 47 genes in bud tissue, 36 genes in seedpod tissue, and 52 genes in flower tissue (see Figure 6A). To further elucidate the transcription patterns of MtERF genes, their expression patterns were clustered across six tissues, as shown in Figure 7. Among these MtERF genes, group A (17 MtERF genes) were highly expressed across all six tissues, while the others exhibited tissuespecific profiles. For example, group B (16 MtERF genes) were highly expressed in root and nodule tissues and not in other tissues. Similarly, group D (six MtERF genes) were specifically expressed in buds, group E (eight MtERF genes) in seedpod and blades and group E (14 MtERF genes) in flowers.
Recently, researchers have demonstrated that microRNA miR172 plays an important role by targeting AP2 TFs in the root nodule symbiosis of legumes (Reynoso et al., 2013;Wang et al., 2014;Nova-Franco et al., 2015). To determine potential regulatory mechanisms of miRNAs and AP2 TFs in the nodulation process of M. truncatula, we submitted M. truncatula miRNA genes (download from miRBase, Van Peer et al., 2014) and MtERF genes to the psRNATarget website (Dai and Zhao, 2011) for the identification of potential target sites. There were four MtERF genes (MtERF101, 106, 108, and 113) of the miR172 family with unambiguously identified cleavage sites (see Figure 8A, and Table S2). Based MtERF genes expression from high-through data, we found that four target MtERF genes had low expression in nodules. Indeed, two of them (MtERF108 and MtERF113) were not expressed at all (Figure 8B).

Expression Responses of MtERF Genes to Abiotic Stress
To investigate the molecular functions of MtERF genes in response to abiotic stress, we performed RNA-seq to detect the expression levels of AP2/ERF TF genes under different stresses, including cold, freezing, drought, salt and ABA. In total, more than 135 M clean reads from control (23,357,742), cold (20,437,484), freezing (23,469,784), drought (22,672,350), salt (21,993,278), and ABA (23,441,046) treated libraries were obtained, respectively; all clean data were submitted to the NCBI SRA database (Accession numbers: SRX1056987-92). After mapping these reads to the M. truncatula genome, we determined that the expression of 75 MtERF genes was detected in at least one library. As shown in Figure 6B, 60 MtERF genes were expressed in the control library, 73 MtERFs under cold stress, 67 MtERFs under freezing stress, 60 MtERFs under drought stress, 62 MtERFs under salt stress, and 58 MtERFs under ABA treatment. Based on expressional profiles of MtERF genes responses to abiotic stresses, they were clustered into three groups, as shown in Figure 9. Among these MtERF genes, group A (14 MtERF genes) were highly expressed by most abiotic stresses, except cold stress. While group B (37 MtERF genes) were highly expressed in response to cold and/or freezing stress, but not in response to other stresses. Group C (24 MtERF genes) were specifically expressed in response to cold stress, fewer in freezing stress and others. Compared to the control library, we identified 48 MtERF genes that were differentially expressed under at least one stress condition, as shown in Figure S5. A total of 35 MtERF genes were differentially expressed under cold stress, while 29 MtERF genes were differentially expressed under freezing stress. A total of 12, 12, and 23 MtERF genes were differentially expressed under drought stress, salt stress and ABA stress, respectively, which was consistent with the MtERF genes clustering results. Notably, four MtERF genes (MtERF022, MtERF023, MtERF043, and MtERF073) were up-regulated under all stress conditions, suggesting that they play important roles in the response of M. truncatula to abiotic stress.
To further confirm the RNA-seq results of these MtERF genes to abiotic stresses, qRT-PCR was performed for nine MtERF genes from the DREB subfamily under abiotic stresses. The expression patterns of most of the MtERF genes in the qRT-PCR analysis were consistent with RNA-Seq analysis, but the magnitude of the fold changes varied between RNA-seq and qRT-PCR experiments (Figure 10). The means of the correlation coefficients of the qRT-PCR validations and the RNA-seq results for the MtERF genes were as high as 0.83, implying that our RNA-seq results were highly reliable.

DISCUSSION
In this study, we performed a comprehensive search for AP2/ERF TF genes throughout the Medicago genome; 123 MtERF genes FIGURE 9 | Expression profile cluster analysis of MtERF genes responses to abiotic stresses.
Frontiers in Plant Science | www.frontiersin.org FIGURE 10 | qRT-PCR validation of MtERF genes in the response to abiotic stress. Fold changes more than 45 were set as 45 for good plot (including two genes, MtERF039, and 043).
were identified and characterized. Previously, Zhang et al. (2013b) identified 37 AP2/ERF genes using EST sequences, whereas we characterized 123 AP2/ERF members from the ERF, AP2, and RAV families. Compared to other plants, the number of AP2/ERF genes in Medicago is slightly lower than that in Arabidopsis (148), soybean (147), cucumber (131), and grapevine (149) and much lower than that in rice (180) and poplar (202). This is most likely due to the lower number of DREB and ERF subfamily members, as there are 98 members in Medicago and 122 in Arabidopsis, 120 in soybean, 122 in grapevine, 145 in rice, and 169 in poplar. By contrast, the number of AP2 and RAV family members is highly conserved among higher plants. For example, 21 members of the AP2 family were identified in this study, while 18, 26, 20, 29, and 26 AP2 genes were identified in Arabidopsis, soybean, grapevine, rice, and poplar, respectively.
In general, TFs always harbor some important conserved domains and motifs for their regulatory function. In MtERFs, 25 motifs were identified based on the MEME results; six motifs (motifs1-6) related to the AP2 domain. Among these motifs, motif1 harbored the largest region of the AP2 domain, including the whole of the β-sheet region and part of the α-helix region (G[VI]R, Gx 4 E, WLG, and AYD elements, as described by Nakano et al., 2006). The other five motifs were short AP2 domains. Motif4 and motif5 contained G[VI]R and Gx 4 E elements, while motif2, motif3, and motif6 contained WLD and AYD elements. All MtERFs contained at least one of the six motifs, indicating that the AP2 domain is highly conserved in MtERF genes. Besides the AP2 domain-related motifs, the other 19 motifs outside the AP2 domain were present in group-specific distributions. For example, motif17 was specifically present in the ERF VI group (also called the ERF B5 group), and contained (L/F) DLN (L/F) xP residues. It was identified as an ERFassociated amphiphilic repression (EAR) motif, essential for repression function (Ohta et al., 2001). Meanwhile, a SP[TV]SVL motif was characterized. It had the potential to be phosphorylated by mitogen activated protein kinase (MAPK) (Nakano et al., 2006). Motif17 may have specific repression functions by virtue of the phosphorylation of the ERF V group. Motif13 contained a unique "EDLL" motif, and was present in the ERF V and IX groups. The "EDLL" motif was previously characterized as a transcriptional activation domain in AP2/ERF TF (Tiwari et al., 2012), implying its participation in the activating function of the ERF V and IX groups (previously called the B6 and B3 groups, respectively). Motif8 was highly conserved in all members of the ERF III group and in four members of the ERF II group (in total, including 30 out of the 45 members of the DREB subfamily). It contained two conserved residue features: LPRP and D[IV]QAA. They have been identified as being essential regulation signatures of the response to various stresses in plants (Albrecht et al., 2001;Qu and Zhu, 2006). Their expression response to abiotic stresses has been shown by RNA-seq analysis. However, except for the three motifs 8, 13 and 17, and the six AP2 domain related motifs (motifs 1-6), the function of the other 16 conserved motifs (65% of the total) identified in present study is uncertain. Their roles need to be further characterized.
AP2/ERF genes in various plant species are differentially expressed in different tissues, indicating that they play important roles in plant tissue development. Genes in the ERF subfamily play many diverse roles, such as functioning in the response to hormonal stimuli and regulating developmental processes in various angiosperms. Meanwhile, members of the AP2 family participate in the regulation of developmental processes, such as flower development and meristem determinacy (Krizek, 2009). In the current study, 75 MtERF genes were found to be expressed in at least one tissue based on high-throughput sequencing data analysis, including DREB (26/50, 52%), ERF (31/48, 64.58%), AP2 (15/21, 71.43%) and RAV (2/3, 66.67%), genes and one Soloist gene (1/1, 100%). Except for two families with fewer members (RAV and Soloist), most members of the ERF and AP2 families participate in the regulation of tissue determinacy. When we compared the expression of MtERF genes in these two families, we found that ERF genes had higher expression levels than AP2 genes, which may be the result of the higher intron content of the AP2 family. Each member of the AP2 family has five to nine introns, while most ERF genes lack introns (74/98) and others contain only one or two introns. Due to the small number of introns, genes in the ERF family respond more quickly and are expressed at higher levels than AP2 family genes. Of the six tissues examined in Medicago, MtERF genes were most highly expressed in flower tissue, which is consistent with previous reports (Krizek, 2009;Matías-Hernández et al., 2014).
Like other leguminous plants, M. truncatula has established a symbiotic relationship with nitrogen fixing rhizobial bacteria, resulting in the formation of specialized lateral organs, called nodules. Nodulation is a complex developmental process involving many molecular signals constituting a genetic regulatory network. TFs, including GRAS, AP2/ERF, NF-Y, and others, have been shown to have important roles in controlling the expression of early nodulation genes (ENODs) and in regulating the later steps of the rhizobial symbiotic interaction (Combier et al., 2006;Andriankaja et al., 2007;Middleton et al., 2007;Vernié et al., 2008;Hirsch et al., 2009;Zanetti et al., 2010;Soyano et al., 2013;Laloum et al., 2014;Baudin et al., 2015). The ERF V group, previously called the B6 group in other studies, contains eight members, four of which (MtERF 52, 53, 55, and 58) were highly expressed in nodule tissue ( Figure S6). MtERF 55 and 58 were only expressed in nodule tissue. This strongly suggests that these genes play unique and important roles in regulating root development and the symbiotic associations with rhizobia. We found that MtERF58 was homologous to ERN and ERN1 ( Table 2). Previous researchers have demonstrated that the ERN gene (similar to ERN1), which is homologous to RAP2.11, interacts with GRAS factors in the regulation of NF-elicited gene transcription rhizobial infection (Andriankaja et al., 2007;Middleton et al., 2007;Vernié et al., 2008). MtERF55 is homologous to PtaERF003, which has been identified as promoting root development in Populus (Trupiano et al., 2013). These results confirm the potential of the ERF V group for the regulation of nodulation. Their functions in the development of other tissues need further investigation.
Nodulation is a complex process in M. truncatula, involving the ERF V group and other AP2/ERF genes, with different regulation patterns. As previously mentioned, many researchers have found that miRNA miR172 is also involved in the nodule development of legumes (Wang et al., 2014;Nova-Franco et al., 2015). It regulates the nodulation process by targeting AP2/ERF TFs, and repressing their expressions. In so doing, it is a negative regulator of root nodulation. In the present study, four AP2/ERF genes (Figure 8A), belonging to the AP2 family, were identified as potential targets of miR172 with unambiguous target sites. Based on expression analysis in silico, we found that these MtERF genes had low levels of expression in nodules, implying that they were potentially down-regulated by miR172 in the nodulation process. These results show the divergent functions of MtERF genes in nodulation. MtERF58, with high levels of expression, is a positive regulator of nodulation in M. truncatula; as confirmed in many studies. However, members of the AP2 family (MtERF 101, 106, 108, 113) were cleaved by miR172, and have low levels of expression in nodules. They may be important repressors of nodulation. This is well-characterized in soybeans and in common beans (Wang et al., 2014;Nova-Franco et al., 2015). However, its function in M. truncatula needs further confirmation.
When adapting to various environmental conditions, plants employ many TF families involved in regulating a wide range of defense responses to environmental clues. AP2/ERF superfamily TFs with ERF domains can bind to GCC-box elements or DRE motifs, thereby regulating gene expression in response to biotic or abiotic stress (Fujimoto et al., 2000;Cao et al., 2001). As previously shown, MtCBF1-3 isolated from M. truncatula, has one AP2 domain, and have key functions in cold acclimatization in M. truncatula and the related species, M. falcate (Pennycooke et al., 2008; Table 2). Similarly, Li et al. (2011) identified the role of MtCBF4 in responding to abiotic stresses, including cold, drought, and salt. Studies in which it was over-expressed conferred improved drought and salt tolerance. Notably, Tayeh et al. (2013) identified a major role for QTL (Mt-FTQTL6) in tolerance against freezing, demonstrated 40% of the phenotypic variation by QTL mapping, and characterized 12 potential MtCBF genes clustering in the Mt-FTQTL6 region. However, the expression profiles of these MtCBF genes in response to freezing stress are still poorly understood, and their function in freezing tolerance in M. truncatula is unknown. In this study, we analyzed the expression of MtERF genes under different stress conditions via transcriptome sequencing. We identified 75 MtERF genes that were expressed under various stress treatments, among which 63 were also expressed during tissue development, while 12 were specifically expressed under stress treatment, as shown in Figure S7. Compared to the control samples, 48 MtERF genes were found to be differentially expressed in response to abiotic stress (see Figure S5), and most MtERF genes were induced by abiotic stress. These MtERF genes include DREB (23/48, 47.92%), ERF (15/48, 31.25%), AP2 (7/48, 14.58%), RAV (2/48, 4.17%), and one Soloist (1/48, 2.08%) gene, the DREB and ERF subfamilies are the largest groups that were responsive to abiotic stress. We found that the DREB and ERF subfamilies were more sensitive to abiotic stress, while ERF and AP2 members were active in tissue development. Among these MtERF genes, 43 MtERF genes were responsive to cold or freezing stress (see Figure S8), while 30 MtERF genes were responsive to salt, drought or ABA-induced stress, as shown in Figure S9. A total of 25 MtERF genes commonly respond to abiotic stress. These include MtERF22, 23, 43, amongst others. MtERF23, which was called MtCBF4 in the report by Li et al. (2011), was previously identified as being induced by cold, drought, salt, and ABA stresses. Our RNA-seq data and qRT-PCR experiments confirmed its high expression under abiotic stresses, (see Figure 10). Of the MtERF genes, 21 were specifically induced or repressed by low temperature stresses, including the previously reported MtCBF 1-3 (Pennycooke et al., 2008). In our study, these were called MtERF 31, 40, and 42, respectively. Most of them are from the DREB subfamily (see Figure S5). Their expression profiles indicate that MtERF genes play more important roles in regulating the response to abiotic stresses M. truncatula, especially DREB subfamily critical function to low temperature stress. Interestingly, 10 MtERF genes (MtERF34-43) are clustered on chromosome 6, and they were tandemly duplicated within an ∼393 kb region, which was previously identified by Tayeh et al. (2013). In the present study, we identified and characterized 10 MtERF genes, belonging to the DREB subfamily. Two previously identified potential MtCBF genes were excluded because the aligning region with Arabidopsis ERF genes was less than 80%. Moreover, our RNA-seq profiles demonstrated that six of these MtERF genes (also called MtCBF genes) were highly up-regulated under both cold and freezing stress (see Figure S5), and qRT-PCR had confirmed their high expression in cold and freezing stresses (MtERF38, 39, and 43, see Figure 10). The expression profiles suggested important and complementary potential roles of the MtERF gene cluster on chromosome 6 under abiotic stresses, especially its major contribution to freezing tolerance in M. truncatula.

CONCLUSIONS
In summary, we identified 123 MtERF genes from the M. truncatula genome sequence. We investigated the classification, evolution, and tissue-specific expression of these MtERF genes, revealing that MtERF genes broadly participate in the regulation of plant tissue development. Meanwhile, we identified 48 candidate MtERF genes that may be involved in abiotic stress responses, their expression profiles were confirmed by qRT-PCR experiment. In particular, a tandem array of MtERF genes on chromosome 6 was identified and found to function in the response to cold and freezing stress, as revealed by RNA-seq analysis. The results of this study will be useful for identifying and characterizing these genes. Further functional analyses of these genes will be performed in the future to enable them to be used for transgenic applications.         Table S1 | List of qRT-PCR validation primers used in the present study.