Genome-wide identification, phylogeny and expression analysis of AP2/ERF transcription factors family in Brachypodium distachyon

Background The AP2/ERF transcription factor is one of the most important gene families in plants, which plays the vital role in regulating plant growth and development as well as in response to diverse stresses. Although AP2/ERFs have been thoroughly characterized in many plant species, little is known about this family in the model plant Brachypodium distachyon, especially those involved in the regulatory network of stress processes. Results In this study, a comprehensive genome-wide search was performed to identify AP2/ERF gene family in Brachypodium and a total of 141 BdAP2/ERFs were obtained. Phylogenetic analysis classified them into four subfamilies, of which 112 belonged to ERF, four to RAV and 24 to AP2 as well as one to soloist subfamily respectively, which was in accordance with the number of AP2 domains and gene structure analysis. Chromosomal localization, gene structure, conserved protein motif and cis-regulatory elements as well as gene duplication events analysis were further performed to systematically investigate the evolutionary features of these BdAP2/ERF genes. Furthermore, the regulatory network between BdAP2/ERF and other genes were constructed using the orthology-based method, and 39 BdAP2/ERFs were found to be involved in the regulatory network and 517 network branches were identified. The expression profiles of BdAP2/ERF during development and under diverse stresses were investigated using the available RNA-seq and microarray data and ten tissue-specific and several stress-responsive BdAP2/ERF genes were identified. Finally, 11 AP2/ERF genes were selected to validate their expressions in different tissues and under different stress treatments using RT-PCR method and results verified that these AP2/ERFs were involved in various developmental and physiological processes. Conclusions This study for the first time reported the characteristics of the BdAP2/ERF family, which will provide the invaluable information for further evolutionary and functional studies of AP2/ERF in Brachypodium, and also contribute to better understanding the molecular basis for development and stresses tolerance in this model species and beyond. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-2968-8) contains supplementary material, which is available to authorized users.


Background
Plant growth, development and productivity are adversely affected by numerous abiotic stresses, such as drought, salt and heat. To survive and flourish under these environmental stresses, plants have developed a complicated response mechanism by repressing or inducing the expression of a series of genes with diverse functions. Transcription factors (TF), as an important group of regulatory proteins, play the central roles in regulation network and signaling pathways of plant development and in response to abiotic stresses. Among them, AP2/ERF (APETALA2/Ethylene Responsive Factor) superfamily is one of the biggest plant TFs, which distinguished by one or two highly conserved ethylene-responsive element-binding factor domains that consisted of 50-60 amino acids [1,2]. Based on sequence similarities and repetitions of AP2 DNA-binding domains, it can be classified into AP2, ERF and RAV families [3]. The members of AP2 family proteins contain two AP2/ERF domains and are further divided into AP2 and AINTEGUMENTA (ANT) monophyletic groups [4,5], while the members of ERF subfamily possesses a AP2/ERF domain with the specific WLG motif and are subdivided into ten group [3], of which Group I to IV belong to the DREB subfamily and group V to X belong to the ERF subfamily. The ERF subfamily is characterized by an additional cisacting element AGCCCGCC of the GCC-box in the promoter regions [6], whereas the DREB subfamily typically binds to dehydration-responsive elementbinding factor, which has a core motif of CCGAC [7]. The RAV family members containing the single AP2/ ERF domain and a specific B3 DNA-binding motif [8]. In addition, other members with an AP2-like domain but lacking additional motifs are often defined as Soloist.
Extensive studies have revealed the crucial role of the AP2/ERF genes playing in plant growth, development and stress responses [4,[9][10][11]. Generally, the AP2 subfamily members were the main factors involving in regulating organ architecture and development, such as leaf epidermal cell determinacy, spikelet meristem differentiation and floral organ patterning [12] as well as seed mass and seed yield [13,14], while the RAV subfamily showed the important functions in plant hormone signal transduction, such as ethylene [15], Brassinosteroid [16], and also involved in response to biotic and aboitic stresses [17,18]. Additionally, the DREB, together with other members in ERF subfamily mainly involved in response to biotic and abiotic stresses, such as water deficit [19], low and high temperature [20,21] and high salinity [22].
B. distachyon, belong to Brachypodium tribe Poaceae family which has a close phylogentic relationships with the major cereal crops, including wheat, barley and rye. It has many favorable features, such as small genome (~300 Mb), diploid accessions, self-fertility, a short lifecycle and easy transformation, which make it an ideal model organism for functional genomic studies of temperate grasses, cereals and biofuel crops [23,24] and now its genome has been completely sequenced [25]. The available genome data facilitated the studies to reveal the gene function and regulation network in this species, and the study of B. distachyon will provide the vital clue for better understand the molecular mechanism of stress response and subsequently improve the abiotic stress tolerance of other cereal crops. So far, the AP2/ERF family has been identified in Arabidopsis [1], Bamboo [26], grapevine [27], maize [28], peach [29] and rice [30]. However, to the best of our knowledge, the systematic identification of AP2/ERF family has not been performed in B. distachyon, limiting the further function analysis of this important gene family.
In this study, a genome-wide bioinformatics analysis was conducted to investigate the genomic organization, phylogenetic relationship and expression profiles of AP2/ERF genes in B. distachyon. The chromosomal localization, gene structures, cis-elements in the promoter region as well as gene duplication and evolutionary mechanisms were subsequently analyzed. By using RNA-seq and microarray expression data, the expression profiling of these identified AP2/ERF genes in different tissue as well as under cold and drought stresses was further investigated. Our study provided a basis for further study on the regulation roles of the AP2/ERF family playing in B. distachyon development and in response to biotic and abiotic stresses, which will not only provide the helpful information on the evolutionary mechanism of this TFs family in plant, but also contribute to revealing the molecular mechanism of development and stresses response in B. distachyon and other cereal crops.

Identification of AP2/ERF gene family in Brachpodium genome
The whole genome data of B. distachyon was available at Ensemble plants database (http://plants.ensembl.org/index.html). The predicted protein sequences were downloaded as the dataset for downstream analysis (v1.0.29). The AP2/ERF domain (PF00847) obtained from PFAM database (http://pfam.xfam.org/) was used as the query for Hidden Markov Model (HMM) search using HMMER 3.0 program with a pre-defined threshold of E <1e −5 . Furhtermore, the AP2/ERF protein sequences ofArabidopsis and rice were obtained from the plant transcription factor database (http://plntfdb.bio.uni-potsdam.de/v3.0/) and then used as query to search against the Brachpodium protein dataset using the BLASTP program with an e-value of 1e-5 and identity of 50 % as the threshold. Furthermore, HMMER and BLAST hits were compared and parsed and then a self-blast of these sequences was performed to remove the redundancy and no any alternative splice variants were considered. After manual correcting, the putative BdAR2/ERF proteins were obtained. Then, the NCBI-CDD web server (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART database (http:// smart.embl-heidelberg.de/webcite) were used to further confirm the predicted BdAR2/ERF genes. The theoretical isoelectric point (PI) and molecular weight (MW) of the obtained proteins were conducted by the compute pI/Mw tool in the ExPASy server (http://www.expasy.org/). The subcellular localization prediction of each gene was predicted using the cello web server (http://cello.life.nctu.edu.tw/).

Multiple sequence alignment and phylogenetic analysis
Multiple sequence alignment was performed using Clustal X v2.0 [31] with the default parameters. An unrooted neighbor joining (NJ) tree with 1000 bootstrap replications was constructed using MEGA 6.0 [32] based on the full-length protein alignment.
Chromosome distribution, gene structure and conserved motif analysis The chromosome distribution of these genes were obtained from the genome annotation information, and then validated by BLASTN search. The exon-intron organizations and splicing phase of these predicted AR2/ERF genes were also investigated based on the annotation file of B. distachyon genome, and then graphically displayed by the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/). Conserved motifs or domains were predicted using the MEME Suite web server (http://meme-suite.org/), with the following parameters: maximum number of motifs set at 25 and optimum with of motifs set from 5 to 200 amino acids.

Promoter analysis and identification of miRNAs targets
The upstream 2 kb genomic DNA sequences of each predicted AR2/ERF genes were extracted from the B. distachyon genome, and then submitted to PLACE database (http://www.dna.affrc.go. jp/PLACE/) to identify the putative cis-regulatory elements in the promoter regions. Furthermore, all the identified AP2/ERF transcripts were searched against the published B. distachyon miRNAs in the miRBase using psRNATarget tool (http://plantgrn.noble.org/psRNATarget/) to predict the AR2/ERF targeted by miRNA.

Gene duplication and synteny analysis
Gene duplication events were identified manually using the method as described by Chen et al. [33]. The segmental duplication events were characterized as copying the whole blocks of genes from one chromosome to another, while contiguous homologous genes with the original duplication on a single chromosome were defined as tandem duplications [34]. For synteny analysis, duplications between B. distachyon AP2/ERF genes, as well as the synteny block of this family among B. distachyon and other 5 grass species (rice, maize, sorghum, foxtail millet and switchgrass) were obtained from the Plant Genome Duplication Database (http://chibba.pgml.uga.edu/duplication/) and the diagrams were visualized using the program Circos v0.67 [35].

Gene expression and network interaction analysis
Microarray data of B. distachyon were obtained from Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/ geo/) and EBI ArrayExpress (https://www.ebi.ac.uk/ arrayexpress/) databases, and then used to detect the expression of the AR2/ERFs in different tissue and in response to abiotic stresses. Additionally, high throughput RNA sequencing data were also retrieved and downloaded from the SRA database (http://www.ncbi.nlm.nih.gov/sra) and then used to detect the differential expression of the AR2/ERF genes by FPKM analysis. A total of 9 RNA data of different tissues at different development stages were used, including anther, pistil, leaves (20 days), seed (5 and 10 days after pollination), endosperm(25 days after pollination), embryo(25 days after pollination), and inflorescence (early and emerging time). Finally, the interaction network which these putative AR2/EFR genes involved in were investigated based on the orthogous genes between B. distachyon and Arabidopsis using the AraNet V2 tool (http:// www.inetbio.org/aranet/) [36].

Plant growth, stress treatment and RT-PCR analysis
Roots, stems, leaves and spikes were collected from twomonths-old Bd21 genotype for RNA extraction and then used for organ-specific expression analysis. The 3 weeks old seedling were subjected to 4°C, 20 % PEG, 150 mM NaCl conditions as cold,drought and salt treatments. After 24 hours treatment, the leaves of plant under these 3 stresses were collected for RNA isolation, respectively. Total RNA was isolated using RNAiso Reagent (TaKaRa, Dalian, China) according to the manufacturer's instructions. Semi-quantitative RT-PCR was employed to determine the transcript levels of 11 randomly selected BdAP2/ERF genes following the method as described by Chen et al. [33]. The primers are listed in Additional file 1: Table S1.

Identification of AP2/ERF family in Brachypodium
Using the method as described above, a total of 141 genes were identified as putative AP2/ERF genes in the Brachypodium genome, accounting for approximately 0.45 % of all annotated Brachypodium genes. Previous study has reported there were 146 AP2/ERF genes in Brachypodium through exploration of genes encoding TF domains to construct TF database [37]. The difference between them were further compared and results found that previous study considered the alternative splices transcripts encoded by the same gene into different AP2/ERF members, which resulted in the increase of the gene number. Since there is no standard nomenclature, the predicted BdAP2/ERF genes were then designated as BdAP2/ERF001 to BdAP2/ERF141 based on their chromosome location and family classification ( Table 1). The detailed sequence information including genomic, transcript, CDS and protein sequence as well as 2 kb upstream has been listed in Additional file 2. Among them, 24 genes containing two repeated AP2/    (Table 2 and Additional file 1: Table S2). Chromosome distribution analysis found that the BdAP2/ERF genes were unevenly distributed on all of the five chromosomes of Brachypodium. In detail, 40 AP2/ERF genes located on the chromosome 3, representing the most abundant regions, followed by the chromosome 1, 2 and 5, with the number of 36, 32 and 18 respectively, while there were only 15 genes on the chromosome 4, which have the minimum number of AP2/ERFs. Interestingly, all the 4 RAV genes located on the chromosome 2, which may be a Brachypodiumspecific feature. The putative proteins of BdAP2/ERFs ranged from 92 to 1338 amino acids in length, with molecular weights (Mw) ranging from 9.8 to 148.6 kDa and theoretical isoelectric points (PI) ranging from 4.33 to 11.63. Subcellular localization analysis indicated that majority of BdAP2/ERFs (108 out of 141, 76.5 %) localized in the nucleus, while 25 genes were predicted to be located in the chloroplast and the remaining 7 genes located in cytoplasmic, mitochondrial, plasma membrane and extra-celluar (Table 1). To further assess the actual existence of these genes identified in this study, all the available Brachypodium expressed sequence tags (EST) were used to search against these genes using the BlastN program. Results showed that most of the AP2/ERFs were supported by EST hits, only 36 genes (25.5 %, 36/111) showed no EST hits. In light of the limit of available ESTs, the not-supported BdAP2/ERF gene might not express under any the used conditions or express with very low level that cannot be detected experimentally.

Phylogenetic relationship, conserved motif and gene structure analysis
To evaluate the evolutionary relationships of BdAP2/ ERF genes, phylogenetic analysis was further conducted based on multiple sequence alignment of all of the BdAP2/ERF together with rice and Arabidopsis AP2/ ERF genes. The phylogenetic tree clustered all the AP2/ ERF genes into three major clades (ERF, AP2 and RAV) depending on their domain composition as described above (Fig. 1). Furthermore, the ERF clades further divided into ten groups. According to the classification criteria in Arabidopsis and rice [3], the ERF superfamily could be further divided into DREB and ERF subfamily. Four groups (group I-IV) of the ERF clades belonged to ERF subfamily, containing 9, 7, 32 and 4 members while the remaining six groups (V-X) were DREB subfamily, having 9, 11, 7, 14, 8 and 6 members, respectively ( Table 2). It's established that DREB subfamily were major factors involved in plant abiotic stress responses and many stress-inducible DREBs have been isolated from numerous plants to date [21-22, 25,]. The identified DREB genes of B. distachyon provided the valuable resource to characterize the stress-responsive genes. Additionally, the bootstrapping values of the nodes in this phylogenetic tree were not very high in every clade, which was consistent with previous studies [3,38]. NJtree reliability was certified by generating another phylogenetic tree by Maximum Parsimony (MP) analysis (Additional file 3: Figure S1), and it was found that nearly all the BdAP2/ERF members were placed within the same topological clusters. Furthermore, the conserved motifs of BdAP2/ERFs were analyzed and compared. A total of 25 conserved motifs were characterized and named as motif1 to motif25 ( Fig. 2 and Additional file 3: Figure S2). Among them, 8 motifs, including motif 1, 2, 3, 4, 6, 7, 16 and 22 were found to be located on the AP2/ERF domain region, while other 17 motifs were corresponded to the regions outside the DNA-binding domain, which was thought to contain either functionally factors, or domains relevant to nuclear localization and transcription regulation [39] (Additional file 1: Table S3). It is noteworthy that proteins within the same group shared one or more motifs that outside the AP2/ERF domain region.  For example, motif 19 and 25 were shared by 9 members in the AP2 subfamily. Motifs 12, 15 and 20 were specifically shared by each member in the RAV subfamily, and the motif 11 was shared by ERF group I as well as motif 8, 9, 10, 14 and 18-23 were specifically presented within the group III members in the ERF subfamily. Finally, the motif 24 was shared by the group V in the DREB subfamily. The proteins within the same subfamilies contained the similar composition of conserved motifs, suggesting the similar function may be shared within each group.
Gene structure analysis of B. distachyon AP2/ERF genes further showed that the member within the subfamily possessed the similar exon-intron structures. As a whole, the number of exon regions ranged from 1 to 12, with an average of 2.65. Most of the ERF subfamily genes (74.33 %) were observed to be intronless, which was consistent with the previous study [1]. In contrast, the AP2 subfamily members contained more intron than ERFs, which had at least four exons (Fig. 3). The highly diverse gene structure suggested that vast differentiation may occur during the B. distachyon genome formation and evolution.

Cis-elements and miRNA targets analysis
In order to understand the possible biological functions and regulation network of these AP2/BdERFs involved in, 2 kb genomic sequences upstream of the 5′-UTR of BdAP2/ERF genes were extracted and used to identify cis-regulatory elements. A total of 276 putative ciselements were found to be presented in at least one BdAP2/ERF gene and only 7 (GT1CONSENSUS, DOFC OREZM, EBOXBNNAPA, MYCCONSENSUSAT, CAA TBOX1, CACTFTPPCA1, WRKY71OS) out of them were presented in the promoter region of all BdAP2/ ERF genes (Additional file 1: Table S4). In addition, 32 cis-elements were detected as gene-specific, such as S2FSORPL21, ABREDISTBBNNAPA and ABREDISTBBN-NAPA were unique to Bradi5g24360, Bradi3g58980 and Bradi5g17620, respectively. The different numbers and types of cis-elements presenting in BdAP2/ERF genes indicated the differential regulatory network which the BdAP2/ERF genes may involve in. Further analysis found that hormones-response (e.g. abscisic acid, gibberellins, auxin, jasmonic acid and ethylene), abiotic stress-related (e.g., drought, extreme temperatures, high salinity, wounding, and disease) and organogenesis-related cis-elements were abundantly presented in the promoter regions of BdAP2/ERF (Additional file 1: Table S5), which indicated that these AP2/ERF genes might have potential functions involving in regulating a variety of stresses response and hormone signaling transduction. Furthermore, the putative microRNAs (miRNAs) targeted BdAP2/ERF genes were also detected in this study and a total of 8 BdAP2/ERFs were predicted to be targeted by seven miRNAs (Additional file 1: Table S6). Although miRNA inhibition mostly involved the transcript cleavage, the BdAP2/ERF006 was predicted to be inhibited to translation. Most predicted microRNA target sites located into CDS region but outside the AP2 domain, whereas for gene BdAP2/ERF051 the cleavage site located in the 3′UTR region. The miRNAs-AP2/ ERF complex identified in this study would be useful in interpreting the post-transcriptional control of gene expression during various stress-induced physiological and cellular processes in B. distachyon as well as other cereal crops.

Gene duplication and synteny analyses of AP2/ERFs between B. distachyon and other three grass species
The tandem and segmental duplication events of BdAP2/ERF genes were investigated through five B. distachyon chromosomes (Fig. 4). Four AP2/ERF gene clusters contained twelve tandem duplicated genes were identified, which located on chromosome 1, 2, 4, respectively. Each cluster had a pair of genes except the cluster located on chromosome 4, which contained six genes belonged to group III of ERF subfamily. Furthermore, 27 pairs of chromosomal segments duplication were also found (Fig. 4). Intriguingly, 3 out of 4 RAV family members showed orthologous relationship, suggesting they may share a common ancestor. To derive the origin and evolutionary relationships of AP2/ERF genes, the comparative syntenic analysis between B. distachyon with other three grass species (rice, sorghum and maize) was performed (Fig. 5a, b, c). Through whole genome-wide syntenic analysis, 44, 49 and 48 % of BdAP2/ERF were identified to be orthologous to rice, sorghum and maize, respectively. Most of BdAP2/ERF genes showed syntenic bias towards particular chromosomes of sorghum, maize, rice, which indicated that the chromosomal rearrangement events like duplication and inversion may predominantly shape the distribution and organization of AP2/ERF genes in these genomes.
The substitution rate of non-synonymous (Ka) versus synonymous (Ks) was an effective measure to examine the positive selection pressure after duplication, wherein Ka/Ks =1 means neutral selection, Ka/Ks <1 stands for purifying selection, and Ka/Ks >1 signifies accelerated evolution with positive selection [40]. Furthermore, the divergence rate of the tandem and segmental duplicated BdAP2/ERF genes was calculated to detect selection influence (Additional file 1: Table S7 and S8). The Ka/Ks ratio for tandem duplicated gene-pairs in B. distachyon AP2/ERF genes ranged from 0.23 to 0.51 with an average of 0.31, whereas Ka/Ks for segmental duplicated gene- These results indicated that the duplicated BdAP2/ERF genes were under strong purifying selection pressure and had gone through substitution elimination and enormous selective constraint by natural selection during the process of evolution since their Ka/Ks ratios were estimated to be lower than one. In addition, the duplication event of these BdAP2/ERF tandem and segmental duplicated genes was estimated to have occurred around~54 and~61 Mya, respectively. Although the BdAP2/ERF gene-pairs of segmental (Ka/Ks = 0.53) and tandem duplication (Ka/Ks = 0.31) events are not under similar evolutionary positive selection pressure, both set of gene pairs revealed that these duplication events may take place simultaneously. Additionally, the Ka/Ks ratios of the orthologous gene-pairs between B. distachyon and other three grass species were also calculated (Additional file 1: Table S9, S10, S11). The average Ka/ Ks value was maximum between B. distachyon and maize (0.47), followed by rice (0.44) and sorghum(0.43), suggesting the genes pairs between B. distachyon and those three grass species appeared to have undergone extensive intense purifying selection.
The divergence time was about 47, 49 and 51Mya for rice, sorghum and maize, respectively. Therefore, it can be concluded that the segmental and tandem duplication events played a major role in evolution and functional divergent of AP2/ERF genes family in B. distachyon as well as other grass species.

Co-expression network between AP2/ERFs and other genes in B. distachyon
To get the preliminary information about the interaction relationship between AP2/ERF and other genes in B. distachyon, we constructed the interaction network of them based on the orthology-based prediction followed the network in Arabidopsis (Fig. 6). A total of 39 AP2/ERFs, with 517 gene pairs of network interactions, were detected. The GO annotations of interacted genes were involved diverse biological process, cellular component and molecular function (Additional file 1: Table S12). For example, symbols BLH6, IAA16, IAA31, ZCW32, LBD41 and HAT3, which play an important role in organ development and response to osmotic stress, were identified as the most closed linked genes with AP2/ERFs. Furthermore, we found AP2/ERF61 and AP2/ERF100 regulated

Expression profiles of BdAP2/ERF genes at different developmental stages and under stresses
The tissue-specific expression profiles of BdAP2/ERF genes at different developmental stages were investigated using RNA-Seq data based on the FPKM analysis. Results found there was high variance in the expression levels among BdAP2/ERF genes (Fig. 7). Several proteins showed relatively high expression in all the tissues, including BdAP2/ERF106, BdAP2/ERF018, BdAP2/ERF113, BdAP2/ ERF108, BdAP2/ERF023, BdAP2/ERF048, BdAP2/ERF0 37, BdAP2/ERF003 and BdAP2/ERF111, suggesting they played the indispensable roles in regulating growth and development. However, three genes, including BdAP2/ ERF119, BdAP2/ERF116 and BdAP2/ERF118 showed very low expression in all the tested organs. Furthermore, the tissue-specific expressed AP2/EFR genes were also identified. BdAP2/ERF083 and BdAP2/ERF064 were found to be predominantly expressed in pistil and leaf, respectively, while BdAP2/ERF005 and BdAP2/ERF006 showed preferential expression in the emerging inflorescence. In addition, six genes namely BdAP2/ERF092, BdAP2/ERF 131, BdAP2/ERF011, BdAP2/ERF012, BdAP2/ERF013 and BdAP2/ERF139 were found to be mainly expressed during pollination, which may contribute to further study of the reproductive growth and seed formation in B. distachyon.
To study the roles of BdAP2/ERF genes in the response to abiotic stresses, the RNA-seq data of B. distachyon under cold treatment (4°C, 24 h) [41] was first used to investigate their expression patterns. Based on the RNA-seq data, a total of 106 BdAP2/ERF genes were detected. Using the fold change method (log2-bias ratio) with more than one fold as criterion, 69 genes were identified as differentially expressed genes (Fig. 8). Among them, 34 genes were up-regulated whereas 35 were down-regulated. Remarkably, BdAP2/ERF122 presented 32 fold up-regulated, while BdAP2/ERF118 showed 122 times down-regulated. Furthermore, the expression profiles of BdAP2/ERF genes under drought stress were also analyzed using the available microarray data [42]. Results found that 16 BdAP2/ERF genes were differentially expressed under drought treatment (Fig. 9). In the expansion zone, five genes were identified as differentially expressed genes, of which one was up-regulated, the remaining four was down-regulated. In the mature zone, we detected eight differentially expressed genes, six genes showed up-regulated while the remaining two showed down-regulated. In the proliferation zone, we characterized three up-regulated genes and three down-regulated genes, respectively. Remarkably, AP2/ERF062 showed down-regulated in all three zones, whereas AP2/ERF022 showed up-regulated expansion in zone and mature zone.
Expression patterns of BdAP2/ERFs in various tissues and under stress treatment by semi-quantitative RT-PCR analysis To further verify the expression of these identified AP2/ ERF genes, 11 BdAP2/ERF genes were randomly selected to detect their expression levels in four tissues and under three stresses treatments through semi-quantitative RT-PCR analysis (Fig. 10). Results showed only one gene (BdAP2/ERF114) was not expressed in these four tissues and the other ten genes were detected to be expressed. Among them, seven genes were found to be expressed in all four tissues with different profiles. In addition, BdAP2/ ERF014 was found to be specifically expressed in stems. BdAP2/ERF076 showed high expression level in stem and leaf, while BdAP2/ERF022 and BdAP2/ERF073 showed high expression level in leaf and spike. Under stress conditions, all of the 11 genes were detected to be expressed. BdAP2/ERF 014, BdAP2/ERF022 and BdAP2/ERF120 were down-regulated under all three stress conditions compared to control, while BdAP2/ERF045, BdAP2/ERF053 and BdAP2/ERF062 was up-regulated under all the treatments. Furthermore, BdAP2/ERF113 showed higher expression under drought treatment, while BdAP2/ERF076 and BdAP2/ERF114 showed high expression under cold and drought treatment respectively, which were consistent with that of RNA-seq and microarray analysis.

Discussion
AP2/ERF superfamily is one of the largest groups of plantspecific transcription factors, which has been widely studied in diverse plant species, such as Arabidopsis, soybean, rice, maize, foxtail millet and switchgrass [1,28,30,43,44]. In this study, we performed a comprehensive search for AP2/ ERF genes throughout Brachypodium genome, and 141 BdAP2/ERF genes were found, accounting for 0.45 % of all the Brachypodium genes, which was similar with the result in rice (0.43 %), maize (0.44 %) and foxtail millet (0.44 %) [43]. While compared to other plants, the number of AP2/ ERF in Brachpodium is much lower than that of rice (174), maize (184) as well as foxtail millet (171), and also slight lower than that of Arabidopsis (148) and grape (149). It has been revealed that the number of AP2/ERF gene family was mainly depending on the number of ERF family members [45]. It's found that there are 112 members in ERF family in Brachypodium, while 122, 132 and 158 in Arabidopsis, rice and maize, respectively. In contrast, the number of AP2 and RAV family members showed no significantly Fig. 7 The expression profiles of BdAP2/ERF genes in different tissue and development stage difference among them, with the value of 28, 24, 34 and 25 in Brachypodium, Arabidopsis, rice and maize respectively. Thus, the lower AP2/ERF gene abundance in Brachpodium may also due to the lower number of ERF and DREB subfamily. Furthermore, the gene density is 0.3972 AP2/ERF genes per Mb in B. distachyon, while the value for rice and Arabidopsis is 0.4047 and 1.1760 respectively. B. distachyon shows closer AP2/ERF density with rice than Arabidopsis, suggesting the specific evolutionary events might occur to regulate the retention and disposition of this gene family between Monocots and Eudicots.
It has been widely revealed that AP2/ERF transcription factors played crucial roles in regulating plant growth, development and response to diverse stresses as well as signal transduction pathway in plants [46]. However, the function of BdAP2/ERFs is not well understood at present. In this study, the expression patterns of these genes in different tissues and under different stresses were systematically investigated to understand their potential function during development and stress response. Results found that a total of 138 BdAP2/ERFs were expressed in at least one tested tissue, indicating they widely involved in growth and development. Compared to AP2 family, the members in ERF family showed higher expression levels in these tissues. We found that the ERF family genes had less intron than AP2 family in Brachypodium, which may cause the quicker response and higher expression of ERF genes during development  [45]. At the same time, BdAP2/ERF genes also showed obvious spatial and temporal expression profiles. For example, BdAP2/ERF064 is specifically expressed in leaf, and BdAP2/ERF005 showed preferential expression in the emerging inflorescence. In addition, six genes having significantly higher expressed during pollination were also identified, which may play the vital role in embryo and endosperm development. AP2/ERF proteins could bind to GCC-box or DRE motifs through the ERF domain, and then regulated the target gene expression under stress conditions [47,48]. Compare to control, 69 BdAP2/ERFs showed differential expression under cold stress and 16 showed differential expression under drought stress, respectively. Among them, BdAP2/ERF 120 (Bradi4g35630) which is a member of DREB subfamily, showed significantly up-regulated under both cold and drought stresses. Previous study have reported Bradi4g35630 encoding a C-repeat binding factor 3-like protein, is a cold-responsive gene and over-expression of this gene could improve the drought, salt and cold tolerance in Brachypodium [49]. Moreover, a total of 5 ACGTATERD1 (element of early responsive to dehydration), 1 DRE1COREZMRAB17 (element of responsive to drought) and 9 MYCCONSENSUSAT (element of responsive to dehydration and cold) cis-elements were identified in the promoter region of BdAP2/ERF120. In addition, BdAP2/ERF053 (Bradi2g27920) is found to be highly expressed in all three stress treatments, which contained 5 LTRECOREATCOR15 (core element of low temperature responsive), 6 EMBP1TAEM (element involving in ABA-mediated stress-signaling pathway) and 1 GT1GMSCAM4 (element required for salt-induced gene expression) cis-elements. We speculated that ciselements were the vital regulators to control the spatial and temporal expression of the BdAP2/ERFs, which integrated other functional proteins with the AP2/ERF transcription factor to form the complex regulatory metabolic network during development and stress response processes [50]. These identified tissue-specific and stress-induced BdAP2/ERF provided the valuable candidates for further functional studies of AP2/ERF genes in B. distachyon as well as in other cereal crops.

Conclusions
Our current study identified and characterized the AP2/ ERF transcription factors in the model grass B. distachyon.  By performing a genome-wide search, a total of 141 BdAP2/ERF genes were obtained. EST hits or full-length cDNA sequences confirmed their actual existence. The chromosome location, exon-intron structure and conserved motif composition as well as phylogenetic relationship of these BdAP2/ERFs were systematically analyzed and compared. BdAP2/ERFs could be classified into four subgroups in accordance with the number of AP2 domains and putative functions. Co-expression network analysis found that 39 BdAP2/ERFs were involved in regulating other B. distachyon genes, and 517 network branches were found. The expression profiles of BdAP2/ ERF genes in various tissues as well as under cold and drought stresses were investigated, and several tissuespecific or stress-induced BdAP2/ERF genes were identified, which could considered as the candidates for further study of their function in plant development and stress response. Our study for the first time reported the organization, structure, evolutionary and expression features of the BdAP2/ERF family, which will facilitate the future functional analysis of BdAP2/ERF genes, and lay the foundation for better understanding the molecular mechanism of plant development and stress physiological processes in B. distachyon and beyond.

Additional files
Additional file 1: Table S1. The primers used in RT-PCR analysis. Table S2. Summary of functional domains presented in the BdAP2/ ERF proteins. Table S3. Detail information of motifs identified in BdAP2/ERF proteins. Table S4. Characteristics of cis-regulatory elements in the promoter region of BdAP2/ERF genes. Table S5. Annotations of 278 putative cis-acting regulatory DNA elements identified in 141 BdERF genes in B. distachyon by PLACE. Table S6. List of putative miRNAs targeted BdAP2/ERF predicted by psRNATarget tool. Table S7. The Ka/Ks ratios and estimated divergence time for tandemly duplicated BdAP2/ERF genes. Table S8. The Ka/Ks ratios and estimated divergence time for segmentally duplicated BdAP2/ERF genes. Table S9. The Ka/Ks ratios and estimated divergence time for orthologous BdAP2/ERF proteins between B. distachyon and rice. Table S10. The Ka/Ks ratios and estimated divergence time for orthologous BdAP2/ERF proteins between B. distachyon and maize. Table S11. The Ka/Ks ratios and estimated divergence time for orthologous BdAP2/ERF proteins between B. distachyon and sorghum. Table S12. Detail information of network of BdAP2/ERF with other genes. Table S13. Detail information of alternative splicing variants. (XLSX 412 kb) Additional file 2: The detailed sequence information of all the BdAP2/ ERF genes, including genomic, transcript, CDS and protein as well as 2 kb upstream sequence. (RAR 245 kb) Additional file 3: Figure S1. A phylogenetic tree of B. distachyon, Arabidopsis and rice AP2/ERF proteins constructed by MP method using MEGA5.0. Then groups are marked I to X. Figure S2. Conserved motifs identified from the AP2/ERF genes in B. distachyon. (ZIP 3757 kb) Abbreviations AP2/ERF, APETALA2/Ethylene responsive factor; DREB, dehydration responsive element binding protein; EST, expressed sequence tag; FPKM, fragments kilobase of exon model per millon mapped reads; Ka, substitution rate of non-synonymous; Ks, substitution rate of synonymous; MP, maximum parsimony; MW, molecular weight; NJ, neighbor joining; PI, isoelectric point; TF, transcription factors