Eremosparton songoricum seeds were collected from the Gurbantunggut Desert Xinjiang province, China (88°24′67′′E, 45°58′11′′N). Seeds were soaked in 98% (v/v) sulfuric acid for 10–15 min to break physical dormancy, after which they were washed with ddH₂O and sown on moist filter paper. Seeds were germinated and grown in covered petri dishes under controlled conditions (25°C, 12 h darkness, 100 mol m^–2 s^–1 light intensity, 60% relative humidity). To each Petri dish (150 mm) containing seeds, 15 mL ddH₂O was added every 2 days.

Two-week-old seedlings with similar primary root lengths were transferred to a culture box containing 400 mL PEG6000 (20% w/v) for the drought treatment. For each time point (0, 6, 12, 24, and 36 h), 10 whole plants were harvested, pooled, and rapidly frozen in liquid nitrogen for RT-qPCR assay, and samples at 0 h served as controls.

A total of 26,442 genes were obtained from E. songoricum genome data. The genome was submitted to CNGB¹ under the project ID: CNP0002419. A Hidden Markov model (HMM) scan with PF00847 (AP2 domain) and PF02362 (B3 domain) was used to mine the AP2/ERF gene sequences in E. songoricum and the other 22 plant genomes (Mistry et al., 2021). The profiles were queried using the HMM search command of TBtools-HMMER software (e-value cutoff 1 × 10^–3) (Chen et al., 2020). The other 22 plant genomes were downloaded from the NCBI genome² and Phytozome 13.0³ (Goodstein et al., 2011). The AP2 domain of AP2/ERFs with a length of approximately 60 amino acids was considered a full-length AP2 domain protein. All predicted peptide sequences contained at least one complete AP2 domain.

The conserved motifs in the AP2/ERF protein sequences were analyzed using the online tool Multiple Expectation Maximization for Motif Elicitation (MEME, version 5.3.3)⁴ with default parameters (Bailey et al., 2015): only motifs with an e-value under 0.1 were retained for further analysis. Conserved motifs were visualized using the Evolview V3 web tool⁵ (Subramanian et al., 2019). Molecular weights and isoelectric points of EsAP2/ERFs were predicted using SwissProt-Expasy.⁶ Subcellular localization analysis was performed using Busca⁷ (Savojardo et al., 2018).

We downloaded 175 Arabidopsis AP2/ERF predicted amino acid sequences from the Plant Transcription Factor Database (PlantTFDB V5.0)⁸ (Tian et al., 2019). The protein sequences of 153 EsAP2/ERFs from E. songoricum and 175 AtAP2/ERFs from A. thaliana were used to construct phylogenetic trees. Multiple sequence alignment was performed using ClustalW (Thompson et al., 1994), phylogenetic trees were constructed by the neighbor-joining method (with 1,000 bootstrap replicates) using MEGA 11, and evolutionary distances were computed using Poisson correction with pairwise deletion (Tamura et al., 2021). The phylogenetic trees were visualized using Evolview web (Subramanian et al., 2019).

Gene duplication type and Ka/Ks were calculated using the MCSCAN module in TBtools (Chen et al., 2020). Information on chromosome distribution (including chromosome length and the starting positions) was retrieved from the local genome annotation and visualized using Circos module in TBtools (Chen et al., 2020). The intron-exon structures of the AP2/ERF genes were assessed based on the annotation file of 21 plant genomes and were visualized using Evolview web (Subramanian et al., 2019). The 3,000-bp region upstream of the transcriptional start site of 153 EsAP2/ERF genes was used for cis-element site analysis using the PlantCARE website.⁹

The STRING database¹⁰ was used to assemble the protein–protein interaction (PPI) networks (Szklarczyk et al., 2021). Based on the annotation of transcriptome and identification in genome, a total of 153 AP2/ERFs and 1,276 other TFs in E. songoricum were used to predict PPI network. The minimum required interaction score was set to medium confidence (0.6). Cytoscape (version 3.8) was used for image processing.

Transcriptome data (CNGB, CNP0003150) was used to analysis the expression abundance of 153 EsAP2/ERF genes under drought stress at 0, 6, 12, 24, and 36 h. Three biological replicates were performed. Differentially expressed genes (DEGs) were determined using FDR < 0.05 and |log2fold change| > 1. TBtools was used to make Heatmap access genes’ relative expression levels.

Total RNA was isolated using the MiniBEST plant RNA kit (Takara, Kyoto, Japan), and first-strand cDNA was synthesized using the PrimeScript™ RT reagent kit (Takara). Twelve EsAP2/ERF genes representative of different subfamilies were selected to explore gene expression patterns. RT-qPCR primers were designed using Primer Premier 5.0, and primer specificities were tested by executing a BLASTn search against local E. songoricum genome data. Primer specificity was further assessed using melting curve analysis of the RT-qPCR. All primers used for RT-qPCR are listed in Supplementary Table 1. RT-qPCR experiments were carried out using the CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, United States) with the SYBR Premix Ex Taq™ kit (Takara). Three biological replicates and three technical replicates of each biological replicate were used for all samples. The RT-qPCR protocol was as follows: initial denaturation at 95°C for 30 s and 40 cycles of 94°C for 5 s and 60–62°C for 30 s. Gene expression levels were calculated relative to the 0 h timepoint using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The EsEF gene was used to normalize RT-qPCR results (Li et al., 2012).

Statistical analyses were performed, and plots were produced using Graph Pad Prism (version 9.0), and all data were analyzed using one-way analysis of variance (ANOVA) at 95% confidence level. Differences were tested using the least significant difference multiple comparison test.

Based on the E. songoricum genome, 153 AP2/ERF genes were identified. A total of 175 AtAP2/ERFs and 153 EsAP2/ERFs were used to construct a phylogenetic tree to classify EsAP2/ERFs. Phylogenetic analyses divided 153 EsAP2/ERF members into four subfamilies: AP2, DREB, ERF, Soloist; RAV subfamily members were not found, based on the genome data (Figure 1). The classification results were confirmed by counting the number of conserved domains (AP2 and B3 domains). Among the confirmed family members, 24 EsAP2/ERFs were assigned to the cluster of the AP2 subfamily based on the phylogenetic tree and presence of the tandemly repeated double AP2 domain (Figure 1). A total of 127 EsAP2/ERFs clustered in the ERF (68 TFs) and DREB (59 TFs) subfamilies, which contained a single AP2 domain. The small Soloist subfamily included two members in E. songoricum (Figure 1). In summary, we identified 24 AP2s, 59 DREBs, 68 ERFs, and two Soloists based on the E. songoricum genome. Considering that some studies reported that transcriptional factors with only one B3 are also considered to belong to the RAV subfamily (Zhao et al., 2017; Chen et al., 2021), we performed a BLASTp search against the E. songoricum genome using the full length of seven Arabidopsis RAVs and 13 soybean RAVs (8 GmRAV only had one B3 domain) as queries, and three hits with high scores were obtained (coverage ≥ 90%; e-value < 1^–10; Supplementary Table 2). The predicted protein lengths of 153 EsAP2/ERF members ranged from 110 to 757 amino acids, molecular weight ranged from 12.7 to 86.5 kD, and the protein isoelectric point ranged from 4.4 to 9.8. Most EsAP2/ERFs (88%) were predicted to be located in the nucleus, and two Soloist members were predicted to be localized in the chloroplast (Supplementary Table 3).

To better understand the evolution of the AP2/ERF family, in addition to E. songoricum and A. thaliana, we also identified the AP2/ERF family in 21 other representative plants ranging from algae and early land plants to angiosperms, including 10 legume species. We found that in three algae (V. carteri, C. braunii, and C. reinhardtii), only AP2, ERF, and Soloist subfamilies were present. The DREB and RAV subfamily began to occur in mosses (Figure 2), and DREB and ERF subfamilies were consistently represented with numerous members among different plant species, whereas the RAV subfamily retained only few members. The Soloist subfamily represented a large proportion in algae, whereas in land plants it comprised few members.

Among ten legume plants, the other nine representative legume species all contained five subfamilies of AP2/ERF, except for E. songoricum, which lacks a canonical RAV, according to the classic classification method. The number of AP2/ERFs in the ten legumes ranged from 140 (L. corniculatus) to 440 (G. max). The proportion of ERF subfamily genes was higher than that of DREBs, which exceeded 40% in the ten legume plants. The DREB subfamily in E. songoricum was represented at the highest proportion (38.56%) among the ten legumes, followed by L. corniculatus (34.29%), and M. truncatula (34.18%). The ratio of the RAV subfamily was <2%, and the numbers of genes ranged from 2 to 6 (Figure 2), Soloist gene numbers varied among legumes, ranging from 1 (C. cajan and A. duranensis) to 16 (A. hypogaea).

The DREB subfamily plays key roles in plant abiotic stress responses and can be further divided into A1–A6 groups. In the present study, 59 EsDREBs were classified into 6 groups by constructing a phylogenetic tree with DREBs of A. thaliana. There were 13, 7, 1, 16, 15, and 7 members in the A1–A6 groups of EsDREBs, respectively (Figure 3A). Fifteen conserved motifs were detected in the protein sequence of the EsDREBs (Figures 3B,C). Motifs 1, 2, and 3 were composed of the AP2 domain, of which motif 1 contained the β3 sheet and α-helix, motif 2 contained β1 and β2 sheets of the AP2 domain; motif 3 represented the C-terminal of the AP2 domain. Other motifs were divergent among different groups, such as motif 10 and motif 5 which occurred only in the A1 group; motif 7 was specific to the A2 group, whereas the EAR repressor motif was detected in the C-terminus of four A-5 type EsDREBs (Figure 3D). Multiple sequence alignment analysis of the AP2 domain showed that EsDREBs shared significant amino acid similarity; 9 out of the 59 amino acids in AP2 domains, including the residues 20-P, 26-R, 28-W, 29-L, 30-G, 38-A, 39-A, 41-A, and 43-D, were completely conserved among all EsDREBs (Figure 3D). In addition, consistent with other plant DREBs, the 14th V residue was completely conserved among all EsDREB members except ES05G02560 which lacks some amino acids in the N-terminal region, while the 19th amino acid had multiple patterns, including E/V/L/A/Q (Figure 3E).

EsERF subfamily can be further classified into six groups (Figure 4A). A total of 13, 3, 19, 9, 8, and 16 EsDREBs were identified for B1–B6, respectively (Figures 4A,B). In total, 20 motifs in 68 EsERF proteins were observed, of which motif 2 contained β1 and β2 sheets, motif 1 was composed of β3 and an α-helix, and motif 3 was the C-terminus of the AP2 domain. Furthermore, ERFs in the same group mostly shared a similar motif composition in addition to the AP2 domain; for example, motif 14 was shared among the B2 group, motif 16 was unique to the B4 group, and motif 5 was only present in the B5 group. The transcription repressor and activator motifs EAR (motif 19) and EDLL (motif 6) occurred at the C-end of B1 and B3 EsDREBs, respectively, and four B3 EsDREBs had EDLL motifs, which were distributed in seven B1 EsDREB members (Figures 4B–D). Three out of the 59 amino acid residues, including 27-W, 37-A and 41-Y, which occurred in the β3 sheet and α-helix of all EsERFs, were completely conserved (Figure 4E). The 14th and 19th conserved amino acids were diverse among ERFs, and the patterns of the 14th and 19th position of AP2 domains were A14 and D19 in most EsERFs, while the 14th position also showed a G/V/I/T/C pattern, and the 19th position was changed to H or N in some EsERFs. The AEIR located between the 14th and 19th position was the same in B1–B5, while the four amino acids changed to SEIR or AEIR/K in B6 group EsERFs (Figure 4E).

To gain insights into gene structural diversity, we examined the exon and intron organization of EsAP2/ERF genes. In E. songoricum, the proportions of genes with introns in the DREB, ERF, AP2, and Soloist subfamilies were 11.9, 35.3, 100, and 100%, respectively. All members of EsAP2 (24/24) and EsSoloist (2/2) had five to ten introns in the AP2 subfamily and five to six in the Soloists subfamily (Figure 5). For the ERF subfamily, 24 of 68 EsERFs had introns; among them, 19 EsERFs had only a single intron, four EsERFs had two introns, and one EsERF (ES02G35200) belonging to the B4 group had seven introns, which also had the longest gene length. For the DREB subfamily, only 7 out of 59 EsDREB members had introns (N = 1–2); 6 of 7 EsDREBs only had a single intron, and the remaining one had two introns (Figure 5).

We further selected six model plants, including G. max, M. truncatula, A. thaliana, O. sativa, P. patens, and V. carteri, to perform intron-exon distribution comparisons of AP2/ERFs (Supplementary Figures 1–6). Similar to E. songoricum, almost all genes belonging to the AP2 and Soloist subfamilies in other plant species had introns, while most DREB and ERF genes were without intron (Figure 6 and Supplementary Figures 1–6). Overall, the intron number distribution of five subfamilies in algae (V. carteri) and moss (P. patens) was more dispersed than in angiosperms. Specifically, for the AP2 subfamily, except for P. patens (4.17% of PpAP2s were intronless), the other six species all had introns (mainly seven to nine) (Figure 6A); however, regarding the ERF subfamily, most ERFs in angiosperm plants lack introns, and the intron number centered at 2 (Figure 6B); a similar intron-exon distribution pattern was found in DREB subfamilies (Figure 6D). For Soloist subfamilies, except for V. carteri (14.29%), Soloists in other species all had introns (mainly six to seven; Figure 6C). A large majority of RAVs in angiosperms were intronless, whereas PpRAVs had two to six introns.

Chromosome mapping analysis indicated that 153 AP2/ERF genes were unevenly mapped on eight chromosomes in E. songoricum. Each chromosome contained at least three EsAP2/ERF genes. Chromosome 1 contained the largest number (33 genes) of EsAP2/ERF genes, followed by chromosomes 4 (25 genes), 5 (22 genes), and 6 (22 genes), whereas chromosome 8 contained only three EsAP2/ERFs with two ERFs and one AP2 gene (Figure 7A). Members of the AP2 and ERF subfamilies were widely distributed across all chromosomes. DREB genes did not occur on chromosome 8, whereas two Soloist genes were located on chromosomes 1 and 7. Chromosome 1 had 13 ERFs, 11 DREBs, 8 AP2s, and one Soloist; chromosome 5 had the largest number of DREB genes (15); chromosomes 1 and 4 had the largest number of ERFs (13) (Figure 7B). In order to elucidate the duplication mechanism of the EsAP2/ERF family, the gene replication types were analyzed using MCScanX, and the results are shown in Supplementary Table 4. A total of 150 duplication events were found, including 71 segmental or WGD (Whole Genome Duplication), 46 dispersed, 27 tandem, and six proximal duplication events. Furthermore, 15 link gene pairs were found on chromosomes 1, 2, 4, 5, and 6 (Figure 7A). The Ka/Ks values of the 15 link gene pairs were all <1, indicating that they evolved under the pressure of purifying selection (Supplementary Table 5).

Transcription factors play key roles in regulating gene expression by interacting with cis-acting elements (Nakashima et al., 2009). To understand the putative function of EsAP2/ERF genes in plant stress responses and the regulatory relationship with plant hormones, we analyzed cis-elements of 3,000 bp of promoter regions of 153 EsAP2/ERF genes. A total of 116 cis-elements were found, including six unnamed and 110 reported elements that were involved in distinct biological functions (Supplementary Table 6). We selected 36 classic previously reported cis-elements (Chai and Subudhi, 2016; Shi et al., 2019; Li et al., 2021) belonging to 22 categories and divided them into two major functional categories, i.e., phytohormone-responsive and abiotic or biotic stress-related (Supplementary Figure 7 and Supplementary Tables 7, 8). The 22 categories occurred 4,496 times in the promoter region of 153 EsAP2/ERF genes, among which MYB, MYC, and ABRE elements were the most enriched (Supplementary Figure 7). Specifically, 99% of the EsAP2/ERF genes contained MYB and MYC binding sites that appeared 1,085 and 620 times, representing 37 and 21% of the abiotic and biotic stress-related cis-elements, respectively (Supplementary Table 7). The ABRE and ERE, involved in ABA and ET responses, were found 541 and 260 times in 110 and 107 EsAP2/ERF genes, respectively, and accounted for 52% of the phytohormone-responsive genes (Supplementary Table 7). Moreover, GA-responsive elements (GARE, P-box, and TATC-box), JA-responsive elements (CGTCA motif, TGACG motif), SA-responsive elements (TCA-element), and auxin-responsive elements (TGA, AuxRR) were present in 122, 110, 102, and 76 EsAP2/ERF genes, respectively (Supplementary Figure 7).

To identify potential interacting proteins with the AP2/ERF transcription factors, a protein - protein interaction (PPI) network was generated with the STRING database (Szklarczyk et al., 2021). A total of 49 EsAP2/ERFs including 23 EsERFs, 19 EsDREBs, 6 EsAP2s, one EsSoloist interacted with 160 other TFs in the network (Supplementary Figure 8 and Supplementary Table 9), in which EsDREBs had the largest number of interacted TFs, then followed by EsERFs. EsAP2/ERFs mainly interacted with WRKYs (17), HSFs (14), MYBs (11), STZs (9), and ABFs (8). The A2 type of DREB (ES02G19290) directly interacted with other 20 TFs, such as ABI5 and ABFs. AP2 (ES04G19400) and Soloist (ES01G17860) only interacted with members of EsAP2/ERFs.

To further determine the involvement of EsAP2/ERF genes in response to drought stress, we analyzed the expression pattern of EsAP2/ERF based on transcriptome data (CNP0003150). As a result, a total of 120 EsAP2/ERF genes were quantified by mapping reads among 153 EsAP2/ERFs (120/153, 78.43%), which contained 53 ERFs, 48 DREBs, 17 AP2s, and 2 Soloists. Eighty-one differentially expressed EsAP2/ERFs genes included 32 ERFs (32/53, 60.4%), 39 DREBs (39/48, 81.2%) and 10 AP2s (10/17, 58.8%). Differentially expressed EsAP2/ERFs were mainly distributed at 6 and 36 h, included 37 and 34 genes, respectively, and notably more EsAP2/ERF genes were up-regulated than down-regulated (75/32) (Supplementary Figure 9). Differentially expressed EsDREBs included 28 up-regulated and 11 down-regulated genes, respectively (Figure 8A). Up-regulated DREB genes consisted of A1, A4, and A5 group, most of genes significantly up-regulated at 6 and 36 h, among them, A4 (ES01G27140) and A5 (ES03G05960) had the highest expression with more than 16-fold (converted by |Log2 FC|) increase compared with that at 0 h. Differentially expressed EsERFs included 23 up-regulated and nine down-regulated genes, respectively (Figure 8B), of which ES02G22910 (B4 type of ERF) had the highest expression with more than 16-fold (converted by |Log2 FC|) increase compared with that at 0 h. EsSoloist and EsAP2 subfamily genes were slight induced after drought stress (Figure 8C).

RT-qPCR results showed that expression profiles of 12 genes, representative of different family members of EsAP2/ERFs including five DREBs, five ERFs, one AP2 and one Soloist demonstrated diverse gene expression patterns during drought stress treatment. The expression patterns of most selected genes were consistent with the transcriptome results (Figure 9). Eight out of twelve genes were strongly up-regulated in drought stress-treated E. songoricum, which highly elevated and reached the peak at 24 or 36 h (Figures 9A–H). Similar like RNA-Seq results, ES03G05530 and ES02G22910 also were strongly induced with drought stress and increased almost 40-fold compared with that at 0 h (Figures 9D,G). Two genes including B2 type of ERF (ES01G00530) and Soloist (ES01G17860) were slightly changed (Figures 9I,J), and two genes (ES01G02260 and ES07G11560) exhibited a decline pattern after drought stress (Figures 9K,L).

The AP2/ERF gene family is one of the largest groups of TFs in plants, and it is characterized by at least one AP2 domain, which plays important roles in plant development and stress responses (Sakuma et al., 2002; Nakano et al., 2006; Agarwal et al., 2017). In the current study, genome-wide analysis of the AP2/ERF gene family from the desert legume E. songoricum was conducted, and 153 putative AP2/ERF genes were identified. These 153 genes covered the AP2, ERF, Soloist, and DREB subfamilies based on domain numbers and sequence similarities. The classic RAV gene nomenclature is characterized by the presence of a C-terminal B3 domain and N-terminal AP2 domain (Kagaya et al., 1999). In the present study, we were unable to identify a canonical RAV protein containing both B3 and AP2 domains; however, the other nine legume plants all had canonical RAV genes (N = 2–6; Figure 2). Previous studies on soybean and rice introduced a classification method with BLASTp searching using AtRAVs protein sequence as query, and genes with only one B3 domain were also considered to be the RAV genes (Zhao et al., 2017; Chen et al., 2021), it is worth noting that these GmRAVs showed different expression patterns under abiotic stresses compared with canonical GmRAVs, implying that they may have different functions. Using this method, we obtained three hits that were probably RAV proteins (Supplementary Table 4); however, further evidence is needed to confirm whether those candidates were EsRAVs. The classification results obtained from the phylogenetic tree and conserved domain numbers were in accordance with the motif and amino acid composition patterns, suggesting that the conserved motifs and amino acid sites can be helpful in classifying the specific subfamily/group of AP2/ERF (Shu et al., 2016; Li X. et al., 2017, Li et al., 2018). Moreover, we found that genes belonging to the same subfamily or specific group shared similar intron-exon organization, which has also been reported in other species (Song et al., 2016; Zhang et al., 2022), indicating that intron-exon structure can also be a useful basis for AP2/ERF family classification.

To explore the origin and evolution of the AP2/ERF gene family in plants, we identified AP2/ERF family genes in 22 other representative plant species ranging from algae to land plants and compared the gene numbers and structures across these 21 plants. Consistent with previous studies, AP2 and ERF already existed in algae, whereas RAV and DREB began to appear in mosses (Mizoi et al., 2012; Song et al., 2016). Our results also supported the phenomenon that DREB and ERF occurred at a large proportion, whereas Soloist and RAV were small subfamilies with very few genes in land plants (Agarwal et al., 2016; Shu et al., 2016; Li X. et al., 2017, Li et al., 2018). In the 10 legume species, the number of AP2/ERFs varied from 140 to 440, and G. max had the largest number of AP2/ERF genes, followed by G. soja and A. hypogaea, which can be explained by independent whole-genome duplications at the species level (Liu et al., 2020; Zhao et al., 2021). Notably, the number of DREB genes in E. songoricum accounted for the highest proportion among the ten legume species, which can be partially explained by the adaptation of E. songoricum to harsh desert environments.

It is known that ancestral species have intron-rich genes, and most plant species experienced extensive loss or insertion of introns due to selective pressure (Roy and Penny, 2007; Rogozin et al., 2012). Low intron gain rates and intron number reduction are common in eukaryotic evolution (Roy and Penny, 2006, 2007). In the current study, we selected seven representative species ranging from algae and land plants to compare AP2/ERF gene structures. Our results showed that almost all AP2 and Soloist subfamily genes of E. songoricum and other plant species had introns, whereas most DREB, ERF, and RAV genes were intronless (Figure 6). Angiosperms such as Arabidopsis, E. songoricum, and O. sativa shared the same intron-exon structures, whereas the number of introns in algae and moss was more variable than in angiosperms. These results are consistent with the intro-exon evolution patterns of AP2/ERF genes of various higher species (Dossa et al., 2016; Zeng et al., 2016; Zhang et al., 2020; Ahmed et al., 2021; Cao et al., 2021).

AP2/ERF genes play important roles in plant abiotic and biotic stress responses, especially in the DREB subfamily which exerts important effects under different abiotic stress conditions, including drought, heat, and salt stress (Sakuma et al., 2002; Shi et al., 2019). Eremosparton songoricum is an extremely drought-tolerant legume species; hence, it is a good material for stress-tolerance gene isolation (Li et al., 2014). Cis element analysis showed that all EsAP2/ERF genes were enriched in ABRE, MYB, and MYC cis-elements, suggesting that EsAP2/ERF genes may participate in the regulation of abiotic and biotic stress responses and plant hormone transduction pathways. Classification results of the EsAP2/ERF gene family showed that DREB genes had the largest proportion, compared to the other nine legume plants, indicating the evolutionary expansion of DREBs for adaptation to harsh desert environments. We previously cloned a DREB2 gene from E. songoricum, which was associated with strong resistance to drought, salinity, cold, and heat (Li et al., 2014, 2016). Moreover, EsDREBs had the largest number of differentially expressed genes and were highly induced in response to drought (39/48, 81.3%), in which 28 EsDREBs genes were significantly upregulated in response to drought stress, ES01G27140 and ES03G05960 were strongly increased by more than 16-fold compared to normal condition (Figure 8A). Meanwhile, EsERFs also remarkably increased after drought stress (Figures 8B,C, 9), indicating that EsDREBs/EsERFs play important roles in E. songoricum response to drought stress and they are good stress tolerance genes for further functional verification.