To accurately collect all members of AhAP2/ERF genes and avoid nonspecific amplification, multiple-database searches were performed. The A. hypogaea cv. Tifrunner genome sequences were downloaded from peanutbase (https://www.peanutbase.org/). The AP2 domain (PF00847) profile was obtained from the pfam database (http://pfam.janelia.org), which was used to match each member of AP2/ERF protein in genomes using HMMER 3.1 software (E-value<1e⁻⁵) (Finn et al., 2011). To avoid the omission of AP2/ERF members, we also performed searches in the Transcription Factor database (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., 2014). All protein sequences acquired were then verified for the AP2 domain by using the SMART (Simple Modular Architecture Research Tool: http://smart.embl-heidelberg.de/) (Letunic et al., 2012) and Pfam (http://pfam.xfam.org) databases (Finn et al., 2008) and NCBI (https://www.ncbi.nlm.nih.gov/cdd). Proteins lacking compete AP2 domains were identified by manual examination. Physichochemical profiling of AP2/ERF genes was performed by using online ExPASy (Gasteiger, 2005; Panu et al., 2012). The subcellular localization analysis of curated AP2/ERF superfamily genes was conducted on the Plant-Ploc server (http://www.csbio.sjtu.edu.cn/bioinf/plant/) (Chou and Shen, 2008).

The AP2 domains were extracted based on results of SMART (Simple Modular Architecture Research Tool: http://smart.embl-heidelberg.de/) (Letunic et al., 2012). Multiple-sequence alignment was executed by DNAMAN and CLUSTAL program (Thompson et al., 1994; Burner and Legendre, 2000). To construct phylogenetic trees, MEGA 7.0 software was used with the neighbor-joining model (1,000 replicates) (Tamura et al., 2013). AP2/ERF family gene names in A. hypogaea cv. Tifrunner were given according to the ascending order of location on chromosomes.

Information of the intron–exon structure was obtained from the reference peanut genome (A. hypogaea cv. Tifrunner, https://www.peanutbase.org/). The Multiple Expectation Maximization for Motif Elicitation program (MEME, http://meme-suite.org/tools/meme) was used to identify potential conserved motifs shared by 185 AhAP2/ERF genes (Bailey et al., 2009). Basic information extraction and preliminary drawing of sequence structure were conducted using TBtools (South China Agricultural University, Guangdong, China) (Chen et al., 2020).

MapChart 2.3 software developed by Wageningen University and Research in Wageningen, Netherlands, is used to locate genes on chromosomes (Voorrips, 2002). Tandem and segmental genes, Ka/Ks values, and circos figures for chromosome locations with AP2/ERF duplication links were completed by TBtools software (South China Agricultural University, Guangdong, China) (Chen et al., 2020). Duplication and divergence time were calculated by the following formula as described by Bertioli et al. (2016): T = Ks / 2 λ ( λ = 8 . 12 × 10 - 9 )

Approximately 1,500-bp upstream sequences of the AhAP2/ERF genes were used to get a better knowledge of the potential function of promoter. PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002) was used to identify the cis-regulatory elements exited in the gene promoters related to stress responses and hormone effects, and these results were visualized by TBtools.

Prediction of the protein interaction network was conducted on the basis of the STRING database (https://string-db.org/. accessed on January 28, 2021) (Szklarczyk et al., 2019). Arabidopsis thaliana L., the well-characterized model plant, was the subject organism (combined score≥ 0.4). PPI networks were constructed by Cytoscape software v 3.8.0 (Shannon et al., 2003).

Seeds of Arachis hypogaea L., ‘YUANZA 9102’, laboratory homozygous material, were sown in 392 cm³ (7 cm long, 7 cm wide, and 8 cm high) pots which were filled with a mixture of vermiculite and perlite (3:1 v/v). Plants were put in a fully controlled growth room (relative humidity: 70%, 16 h/8 h light/dark; 30°C/28°C day/night; light intensity 17,000 lx). Watering was stopped in one part of pots (drought treatment) when seedlings were 5 weeks old, whereas the watering regime remained unchanged in the control plants (every 4 days). Roots were collected from the control and treatment groups every 4 days from seedlings aged 5–9 weeks. All samples were frozen in liquid nitrogen immediately and then stored at -80°C for RNA isolation.

RNAprep Pure Plant Plus Kit (Tiangen Biotech, Co., Beijing, China) was utilized to extract RNA from control and treated peanut samples. The cDNA was prepared by following the user manual of the PrimeScrit™ RT Kit with gDNA eraser (perfect real-time, Takara Biomedical Technology, Ltd., Beijing, China). Thirty-five genes from each family were designed by Primer Premier 5.0 and are shown in Supplementary Table S1. The alcohol dehydrogenase class III (AhADH3, Arahy. VYWU26.2) (forward primer: 5′-GAC​GCT​TGG​CGA​GAT​CAA​CA-3′, reverse primer: 5′-AAC​CGG​ACA​ACC​ACC​ACA​TG-3′) was selected as the internal reference control (Brand and Hovav, 2010). Subsequently, quantitative RT-PCR was performed using the ABI QS5 qRT-PCR detection system (ABI, United States) and SYBR Green Kit (Tiangen, Beijing, China). An ABI QS5 real-time PCR system was used under the following procedure: 95°C for 15 min, followed by 40 cycles of 95°C for 10 s, and 60°C for 32 s in a 20 µl volume. Each PCR assay was carried out in three biological replicates, of which each replicate corresponded to three technical repeats. Relative expression levels of the genes were calculated using the 2^-△△Ct method (Livak and Schmittgen, 2001).

A total of 185 unigenes with the AP2 domain were characterized as AP2/ERF TFs in A. hypogaea cv. Tifrunner (Supplementary Table S2). Depending on the sequence characteristics, the AP2 domains, phylogenetic tree analysis, and the classification system established by the group of Yamaguchi-Shinozaki (Tamura et al., 2011) and Shinshi (Nakano et al., 2006), the AP2/ERF superfamily genes are mainly classified into AP2 (APETALA2), ERF (VI-X, ethylene-responsive-element-binding), DREB (I-V, dehydration-responsive-element-binding), RAV (related to ABI3/VP), and Soloist (few unclassified factors) subfamilies (Figure 1 and Supplementary Figure S1). Among these, the 27 to encode two AP2 domains and the 4 to one AP2 domain together with one B3 domain were thus assigned to the AP2 and RAV families, respectively. Based on the similarity of amino acid sequences with the AP2 domain, the 117 genes were further assigned to the ERF 76) and DREB 41) subfamilies, respectively. Thirty-two members with a single AP2 domain but were distinct from the ERF or DREB subfamily were classified into the AP2 subfamily (Supplementary Table S2). The remaining five genes with independent clades from others were identified as Soloist genes. Subsequently, all superfamily members were named according to the order on the chromosomes of each family member to distinguish from each other for the study (Supplementary Table S3).

To evaluate the phylogenic relationship and classification of the ERF and DREB subfamily, multiple-alignment analyses was performed on the protein sequences with the AP2 domain acquired from peanut (117), Arabidopsis (139), and rice (139), as suggested by Nakano et al. (2006) (Supplementary Figures S2, S3). The NJ phylogenetic tree divided the ERF and DREB subfamilies of peanut and Arabidopsis into 10 subgroups (DREB-I-V and ERF-VI to X) following the classification as described by Nakano et al. (2006) (Supplementary Figure S4). The phylogenetic tree of ERF and DREB subfamily proteins of peanut and rice also exhibited similar results (Supplementary Figure S5). Overall, current findings of the phylogenetic tree demonstrated that classification of the peanut ERF and DREB subfamily proteins is similar to the Arabidopsis and rice ERF family (Supplementary Table S4).

Molecular property analysis showed that MW of AP2/ERF superfamily members varied from 14.46 to 82.54 kDa. Most of AP2/ERF superfamily genes showed their localization in the nucleus. Moreover, the negative GRAVITY values suggested the globular hydrophilic nature of the AP2/ERF proteins. Interestingly, the members of the same families or clades shared similar physical properties, indicating the functions conservatively in the same clades and differentially among subfamilies.

Structural analysis of AP2/ERF genes is helpful for us to fully understand the conservative characteristics of peanut AP2/ERF protein and analyze its evolutionary differences. Numbers of introns varied among AP2/ERF subfamilies (Figure 2). All the AP2 family genes contain 3–10 introns, whereas most members of ERF, RAV, and Soloist subfamilies have only 1–2 introns or do not possess introns. A similar phenomenon was discovered in Arabidopsis thaliana L. and Cucumis sativus L., where most Arabidopsis genes of the ERF family do not possess introns (Sakuma et al., 2002), and 83% of CsERF genes do not have intron (Ito et al., 2014). What is more, members of the ERF family clustered into one branch have a similar gene structure (Figure 2). There was some subgroup specificity: members of groups DREB-III, DREB-II, DREB-I, and ERF-VI did not possess introns, and the genes owning two introns were in the groups ERF-VII and ERF-X, possibly attributed to number changes of introns during evolution, whereas the number and position of the introns were relatively conserved in the same group of plant species.

To provide evidence for the classification and the functional conversation of AP2/ERF superfamily genes among different groups, 10 conserved motifs (motifs 1–10) were analyzed by using MEME software (Figure 2). Among these, motifs 1 and 2 were dominantly present in the AP2 domain regions of all family members. The proteins of the same group showed identical numbers and arrangements of motifs, which are different among the various clades. For example, nine motifs were detected in AP2 and ERF families, and four motifs in RAV and Soloist families. Furthermore, the number and arrangement of motifs in the RAV (motifs 2, 1, 4, 10), ERF, and Soloist (motifs 3, 2, 1, 4) families showed high similarity. Motif 8 was only detected in few AP2 family members, signifying the meticulousness of the motif in the AP2 family. Remarkably, motifs in the same group showed great similarity, indicating the functional conservation in different groups. Comparing the intron–exon structure and conserved motif analysis, it is clear that the members of the same group showed great similarity of characteristics, indicating that most of AhAP2/ERF genes were highly conserved among groups.

AP2/ERF genes showed random distribution on the 20 chromosomes of the peanut genome (Figure 3 and Supplementary Table S3). Maximum numbers of AhAP2/ERF genes (16 genes) were located on Chromosomes 6 and 15, while Chr17 had the least number of genes (4 genes). Other chromosomes had a random number of allocated AhAP2/ERF genes (5–13 genes). Interestingly, most members are distributed at both ends of the chromosome. A similar study on Arabidopsis and other species showed consistent findings (Nakano et al., 2006; Guo et al., 2016). This location similarity of genes on chromosomes indicates functional consistency.

Gene duplication events (segmental or tandem) play a significant role in the expansion and evolution of gene families in plant species (Baloglu, 2014). In total, 99 duplicated gene pairs, which were also named as homoeologous genes, were identified: 14 tandem and 85 segmental duplications (Figure 4). Segmental gene duplication mainly occurred in the A. hypogaea cv. Tifrunner genome rather than tandem duplication event (Figure 4 and Supplementary Table S5). No duplication events occurred in group I of the ERF family, whereas more segmental duplication events occurred in other groups of the ERF family, implying that the biggest members of the families might arise from a higher frequency of segmental duplication, when adapting to various environmental shifts. In contrast, tandem duplication has a confined benefaction to the gene family expansion as compared to the segmental duplication. Similar studies on Arabidopsis, rice, sorghum (Wang et al., 2011), common bean (Kavas et al., 2015; 2016), and cucumber (Baloglu, 2014) showed compatible findings. In general, segmental duplication might be the main driving force for the AP2/ERF gene family expansion in peanut genome.

Moreover, comparative orthologous analysis was conducted among A. hypogaea cv. Tifrunner, A. duranensis, A. ipaensis, Medicago truncatula L., and Glycine max L. to characterize the evolutionary patterns of AhAP2/ERF genes with Leguminosae species (Figure 5). In total, 140, 133, 314, and 145 orthologous gene pairs were found with A. duranensis, A. ipaensis, Medicago truncatula L., and Glycine max L., respectively (Supplementary Tables S6–S9). The Ka/Ks for segmental duplication was 0.02–0.92 with an average of 0.28, while the ratio of tandem duplication ranged from 0.29 to 0.68 with an average of 0.42 (Supplementary Table S5). These segmental and tandem duplications may occur in ∼3.82–43.68 Mya, respectively. In addition, the Ka/Ks ratio of ortholog gene pairs between A. hypogaea cv. Tifrunner and Medicago truncatula L. (0.25) and A. hypogaea cv. Tifrunner and Glycine max L. (0.24) were strongly subjected to pure selection (Lynch and Conery, 2000). The divergence times were 64.07 and 66.44 Mya for Medicago truncatula L. and Glycine max L., respectively.

Specific cis-element motifs can be recognized by TFs and participate in gene expression regulation. In order to further study the potential regulatory mechanism of the AP2/ERF gene in diversified biological processes, especially in plant drought stress response, the 1.5 kb upstream sequence of the AP2/ERF gene translation start site was submitted to the PLANTCARE database to detect cis elements (Cui et al., 2018; Meng et al., 2020). A total of 56 known cis-elements (30 light-related elements, 11 hormone-related elements, 8 tissue-specific elements, and 7 stress-related elements) were detected (Figure 6 and Supplementary Table S10). ABRE, AuxRR core, CGTCA motifs, GARE motifs, O2 site, P-box, TATC-box, TCA element, TGA element, and TGACG motifs involved in hormonal responses are found in 63.8, 5.9, 49.2, 9.2%,18.9, 18.4, 10.8, 32.4, 25.9, and 49.2% of AhAP2/ERF promoters, respectively. Meanwhile, there are a large number of stress-related elements, including MBS (drought inducibility), TC-rich repeats (defense and stress responsiveness), WUN motif (wound responsiveness), LTR (low-temperature responsiveness), ARE (essential for the anaerobic induction), and GC motif (anoxic specific inducibility) (Supplementary Table S10). Moreover, there was a divergence in the percentage of cis-acting elements in promoter regions of various families (Figure 6D). For example, all the RAV family members contained ARE and O2-site elements, whereas 74.6, 80.0, and 60.7% family members of AP2, Soloist, and ERF families possessed ARE elements, and 10.2, 20, and 20.5% of those families exhibited O2-site elements in the promoter region, respectively. Notably, MBS, an important cis-element related to the plant drought-inducibility process, was detected in the promoters of 25.4% of AP2, 50% of RAV, 60% of Soloist, and 22.2% of ERF family members. As a major hormone in plant response to drought stress, ABRE possessed 57.6% of AP2, 75% of RAV, 40% of Soloist, and 67.5% of ERF family members. TC-rich repeats, a cis-acting element involved in defense and stress responsiveness, were discovered only in the promoters of AP2 and ERF family members. WUN motif, a wound-responsive element, was only detected in the ERF family. The variants in the characterization of cis-acting elements implied the functional discrepancy in different families.

What is more, certain ERF members in groups III (AhDREB5/29/41), IV (AhDREB15/27/39), V (AhDREB21/26), VI (AhERF25/52), VII (AhERF56), VIII (AhERF3/19/30/34/43/53/62/64/65/70/74), and IX (AhERF5/31/71/72) possess a relatively large number of MBS elements in the promoter regions, implying the members’ main role in the regulation of plant drought responsiveness and the function variations among groups (d). In other words, cis-acting elements in the same group showed great similarity, indicating the functional conservation in the same groups or clades (Figure 6D).

To understand the synergy among peanut AP2/ERF TFs during their regulatory process, an interaction network was drawn using Arabidopsis ortholog genes (Figure 7). A sum of 31 gene pairs with a combined score over zero value was deliberated to have the interaction with others (Supplementary Table S11). AhAP2-29, AhAP2-37, AhERF-47, and AhRAV-1 had more than three nodes and protein pairs and were involved in more powerful crossing networks, suggesting their core role in peanut. However, other members were only regulated by a few numbers of genes, indicating its less important role in transcriptional level regulation (Figure 7). Interestingly, proteins from different families showed complex interaction. For example, core members AP2-39 or RAV-3 interact with genes from DREB and ERF subfamilies, implying the exitance of functional interaction among subfamilies. These interrelationships will provide a reference for studying the regulatory functions of AhAP2/ERF genes in peanut.

For better knowing the possible regulatory roles of AP2/ERF family genes in peanut response to drought stress, 35 representative genes from AP2 (5), RAV (3), and ERF (27, members from each group) were selected to verify whether their expression levels would be induced under drought stress conditions by qRT-PCR (Figure 8), especially for the ones that possess MBS elements (drought-inducibility), ABRE, TC-rich elements, ARE, and WUN-motif elements in the promoter regions.

The treatment group with no watering was used to stimulate the expression of plant defense genes. The expression analysis of AhAP2/ERF genes responsive to drought stress could present useful information to identify their implied role as candidate genes to mitigate drought stress severity. As shown in Figure 8, under the no watering condition, except AhDREB-3, AhDREB8, AhERF-4, AhERF-13, AhERF-66, AhERF-72, AhERF-75, AhAP2-1, AhAP2-6, AhAP2-13, AhRAV-1, AhRAV-2, and AhRAV-3 exhibited downregulation, whereas the remaining genes were significantly upregulated and subsequently downregulated, at all time points compared with those at 0 h. Furthermore, the peak expression of AhDREB-1, AhDREB-4, AhDREB-5, AhDREB-33, AhERF-7, AhERF-8, AhERF-69, and AhAP2-19 was discovered on the 4th day, but that of AhDREB-14, AhDREB-18, AhDREB-20, AhDREB-23, AhERF-47, AhAP2-6, and AhAP2-37 was detected at the 8th day. Interestingly, the members in the same group showed semblable expression trends, which indicates the function consistency. Notably, collinear genes AhDREB-9 and AhDREB-23, AhAP-26 and AhAP2-37, AhRAV-1, and AhRAV-3 showed highly similar expression patterns under drought stress treatment, indicating that their biological functions also have a certain similarity.