To identify R2R3-MYBs in peanut, the A. hypogaea cv. Tifrunner protein sequences were retrieved from the PeanutBase (https://peanutbase.org). The MYB DNA-binding domain (PF00249) was exploited for the identification of R2R3-MYB genes in the peanut genome by using the HMMER 3.3.2 program at a standard E value <1×10^–5 (http://pfam.xfam.org/search#tabview=tab1). In total, 204 predicted gene models were found with two consecutive repeats of the MYB domain. All but five (199) were also retrieved when performing a BLASTP analysis using the previously identified genes of A. hypogaea against the whole A. thaliana R2R3-MYB gene data set (Dubos et al., 2010) with a cut-off e-value of e^-40. Subsequently, protein sequences were evaluated for the presence of the MYB domain against the repository of the NCBI CDD (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART databases (http://smart.embl-heidelberg.de/). Finally, 196 R2R3-MYB genes were obtained after eliminating incomplete and uncertain sequences. The predicted molecular weights and the theoretical isoelectric point (pI) were obtained by the Expasy proteomics server (https://web.expasy.org/compute_pi/). The prediction of transmembrane helices in AhR2R3-MYB proteins was analyzed by TMHMM - 2.0 (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0).

The 196 R2R3-MYB protein sequences in peanuts were performed by multiple sequence alignment using ClustalX (Larkin et al., 2007). The alignment of 196 peanut R2R3-MYB domains was performed using ClustalX and DNAMAN (Lynnon Biosoft). The amino acid residue distributions of the conserved MYB domains of AhR2R3-MYBs were created using the WebLogo program with default parameters (http://weblogo.berkeley.edu/logo.cgi) (Crooks et al., 2004).

The protein sequences of the 126 A.thaliana R2R3-MYBs were downloaded from the TAIR (http://www.arabidopsis.org/). The protein sequences of R2R3-MYB proteins from A. hypogaea and A.thaliana (protein sequence information is listed in Supplementary Table 1) were aligned by the MAFFT with the FFT-NS-i algorithm (Katoh et al., 2019), and the multiple sequence alignments were used for phylogenetic analysis. The phylogenetic tree was constructed by the neighbor-joining method of MEGA 7.0 with 1000 bootstrap replicates based on the p-distance model and pairwise deletion for gap treatment (Kumar et al., 2016). The phylogenetic tree was retouched by FigTree (http://tree.bio.ed.ac.uk/software/figtree/). For the construction of the phylogenetic trees of R2R3-MYB proteins from Arachis. hypogaea the same method described above was adopted.

Information on the intron-exon position and splicing sites for each gene model in the corresponding chromosome scaffold was downloaded from PeanutBase. Furthermore, the MEME program (https://meme-suite.org/meme/meme_5.3.3/tools/meme) was used for the identification of motifs of 196 R2R3-MYB protein sequences in peanut. The optimized parameters of MEME were employed as follows: the number of motifs that MEME find, 10; and the optimum width of each motif, 6–60 residues (Bailey et al., 2009). The gene exon-intron pattern and MEME results were also visualized by CFVisual_V2.1.5 (Chen et al., 2022).

The information on chromosome length and R2R3-MYB gene locations was acquired from the PeanutBase (https://peanutbase.org) (Supplementary Table 2), and the figure was created by TBtools (Chen et al., 2020). The whole-genome sequences and annotation documents of A. hypogaea were downloaded to PeanutBase (https://peanutbase.org). Then, the One Step MCScanx program of TBtools was executed to analyze the synteny relationships of genomes. We identified gene pairs with physical distance within 100 kb, with no more than 10 genes spaced in between, and in the same subgroup as tandem repeat gene pairs according to Hanada et al. (2008). The duplication pattern of the AhR2R3-MYB genes was visualized by the Amazing Super Circos package of TBtools. The Ka/Ks value was completed by the Simple Ka/Ks Calculator program of TBtools. Duplication time was calculated by the following formula as described by Bertioli et al. (2016): T = Ks/2λ (λ = 8.12 × 10⁻⁹).

The peanut (Arachis hypogaea L.) cultivar ‘Changhua18’, a germplasm resource preserved in our laboratory, was planted in a pot with a 1:1 mixture of nutrient soil and vermiculite, 450 mm long, 335 mm wide, and 170 mm high. When the plants were grown for about 8 weeks, the treatment group was subjected to waterlogging treatment according to Zeng et al. (2021), while the control group was kept under normal growth conditions. The samples (three seedlings per repeat) were collected at 0 h, 6 h, 24 h, 3 days, and 5 days after treatment, respectively. Subsequently, the samples were rapidly frozen using liquid nitrogen and stored at -80°C for RNAseq.

The raw read counts in Arachis hypogaea RNAseq samples were downloaded from PeanutBase (https://peanutbase.org). The data were obtained from 22 tissues at different developmental stages in peanut with three biological repeats (Clevenger et al., 2016) with all raw data deposited as BioSamples SAMN03944933–SAMN03944990. HTSeq was used to generate raw reads that were uniquely mapped on the Arachis hypogaea genome. StringTie and Ballgown were used for the FPKM calculation (Pertea et al., 2016). The transcript profiles for AhR2R3-MYB genes were displayed in TBtools (Chen et al., 2020).

These cDNA libraries generated from the samples were sequenced by Metware Biotechnology Ltd. on the Illumina sequencing platform. (Wuhan, China). Download the reference genome and its annotation files from NCBI (https://ftp.ncbi.nlm.nih.gov/genomes/), use HISAT v2.1.0 to construct the index and compare clean reads to the reference genome. The featureCounts v1.6.2/StringTie v1.3.4d was used to calculate the gene alignment and FPKM. Gene expression patterns were also charted by TBtools. Quantitative RT-PCR (qRT-PCR) was performed to verify the transcriptome data. RNAprep Pure Plant Plus Kit (Tiangen Biotech, Co., Beijing, China) was utilized to extract RNA from control and treated peanut samples. The PrimeScrit^TM RT Kit with gDNA eraser user manual was used to prepare the cDNA (perfect real-time, Takara Biomedical Technology, Ltd., Beijing, China). The primers were designed by TBtools and were shown in Supplementary Table 3. Subsequently, qRT-PCR was performed using the ABI 7500 qRT-PCR detection system (ABI, United States) with SYBR Green Kit (Tiangen, Beijing, China). The ABI 7500 real-time PCR program was 95°C for 15 min, followed by 40 cycles of 95° C for 10 s, and 60°C for 30 s in a 20 µl volume. Three technical repeats of qRT-PCR were carried out, and the relative expression level was determined using 2^-△△Ct technique.The results of differential expression analysis between homeolog pairs of A and B subgenomes across tissues and pod developmental stages in Arachis hypogaea were downloaded from PeanutBase (https://peanutbase.org) (Bertioli et al., 2019). Differential gene expression analysis was performed using the DESeq2 package (v1.14.1) with log₂ fold change >= 1 and Benjamini-Hochberg adjusted P-value < 0.05 as the statistical cutoff for differentially expressed genes. Charts generated by GraphPad Prism 9.

Total branch number, length of root and shoot, pod size and RS phenotypes were assessed using a randomized complete block design and replicated in five environments. SNPs of the R2R3-MYBs were obtained from transcriptome data set of a peanut germplasm population with 146 accessions (alleles in each polymorphism with minor allele frequency >0.05). Association analysis was performed with TASSEL 3.0 using an MLM Q + K model. Figures were generated by GraphPad Prism 9.

A total of 196 R2R3-MYB genes were obtained after genome-wide screening and exclusion in A. hypogaea cv. Tifrunner. The 196 genes were named AhMYB001-AhMYB196 according to their physical location on chromosome. The mRNA length of the AhR2R3-MYBs ranged from 753 to 3712 bps (Supplementary Table 4); the proteins length ranged from 192 to 916 amino acids, with predicted molecular weights from 21.95 to 103.72 kDa; the theoretical isoelectric point (pI) ranged from 4.60 to 10.32. Furthermore, none of the proteins were predicted to contain transmembrane domains (Supplementary Table 4).

To investigate the MYB conserved domains of R2R3-MYB transcription factors in peanut, we performed WebLogo and multiple alignment (ClustalX) analysis using amino acid sequences of the R2R3 repeats (Supplementary Figure 1). The results showed that the R2 and R3 repeat consisted of two repeats of approximately 51 amino acid residues, and all of them contained highly conserved tryptophan residues (Trp, W) (Supplementary Figure 2). The R2 repeat contains three highly conserved W residues at positions 5, 26, and 47, while only two highly conserved W residues were uncovered at positions 24 and 43 in the R3 repeat, and the W residue at position 5 was generally replaced by phenylalanine or isoleucine (Phe, F or Ile, I). Tryptophan residues formed a hydrophobic core maintaining the stability of the helix-turn-helix (HTH) structure in the DNA binding domain, and the basic and polar amino acids adjacent to that tryptophan were thought to be directly involved in DNA binding (Saikumar et al., 1990). Several highly conserved polar amino acid residues were also found around the third W residue in each repeat, for example, asparagine (Asn, N), arginine (Arg, R), and lysine (Lys, K) around the third tryptophan in the R2 repeat, and K and N near the third W residue in the R3 repeat, which might be directly involved in binding to DNA. These polar amino acids, which are conserved around tryptophan, were also highly conserved in different species, such as Arabidopsis and tomato (Stracke et al., 2001; Li et al., 2016). Overall, most of the conserved amino acid residues were mainly distributed between the second and third conserved tryptophan residues, with the first tryptophan residue being relatively less conserved. Therefore, those conserved residues probably maintain the function of the DNA binding domain together with the conserved tryptophan.

To investigate the evolutionary relationship between the R2R3-MYB genes in peanuts and Arabidopsis, the 196 predicted AhR2R3-MYB proteins were subjected to multiple sequence alignment along with 126 Arabidopsis R2R3-MYB proteins, and their evolutionary relationships were inferred by constructing a neighbor-joining phylogenetic tree. The 322 R2R3-MYB proteins were classified into 48 subgroups based on the topology and bootstrap value of the phylogenetic tree (Figure 1 and Supplementary Figure 3). Since the R2R3-MYB genes have been intensively studied in A. thaliana and most subgroups contained at least one AtR2R3-MYB gene, we named these subgroups according to the nomenclature of (Kranz et al., 1998) revised by Stracke et al. (2001) and Dubos et al. (2010). When a subgroup name was not yet determined in A. thaliana, we named the subgroup after the member of A. thaliana with the most distinct functional characteristics. In general, the phylogenetic characteristics of A. thaliana described in this paper were generally consistent with those described previously (Stracke et al., 2001; Dubos et al., 2010). The only exception was the A. thaliana genes of the subgroups 20 and 25, were split in two respectively (S20a and S20b; S25a and S25b), and the genes from 10 and 24 were merged (S10&24). As shown in Figure 1, 40 subgroups contained at least one gene from peanut, and the other eight subgroups contained genes only from peanut, and they were named as new subgroups 1-8 in this study. The distribution of AhR2R3-MYBs in the 40 subgroups was biased, varying from one (NS-8) to 17 members (S14). Notably, the number of R2R3-MYB genes in almost all subgroups was unbiased in both subgenomes, revealing a close association between two subgenomes of cultivated peanut as described by (Bertioli et al., 2016; Bertioli et al., 2019). The topology of the neighbor-joining tree for AhR2R3-MYB genes was in good agreement with the subgroup described above (Supplementary Figure 4).

To investigate the relationship between subgroup classification and function of the peanut R2R3-MYBs, 10 conserved motifs were identified in the AhR2R3-MYBs through MEME program search (Supplementary Figure 5). The DNA binding domain of AhR2R3-MYB was represented by motifs 3, 6, 1, 2. Motifs 1 and 2 contained the amino acid sequence of the third helix forming the MYB domain, which is involved in the recognition and binding of cis-acting elements (Ogata et al., 1996; Jia et al., 2004). Motifs 5 and 8 were only presented in SAtM88, while motif 9 only presented S20a and S20b, suggesting potential specific functions of these subgroups (Figures 2A, B). In general, most of the motif compositions of members in the same subgroup were similar at N-terminal, but differ at C-terminal, and the motif compositions of the members in different subgroups were not identical.

Exon and intron structure analysis showed that all AhR2R3-MYB genes possessed 1 to 21 exons (Figure 2C). In general, the size of introns was variable, but the locus and phase of introns were relatively conservative among subgroups (Jiang et al., 2004). Most of the R2R3-MYB genes (77%) contained two conserved introns, 14% contained one conserved intron, and the rest showed a different number. Genes in the same subgroup have similar gene structures and highly conserved in intron phasing (Figures 2A, C). The multiple introns of a gene provide the opportunity to selectively splice and provide variant proteins that may play different roles in biological processes (Min et al., 2015). Most members of the subgroups are intron-poor (Contain three or fewer) or intron-less genes. However, we found that genes in subgroups NS-3, S18, and SAtM88 possessed an abundant number of introns, and multiple transcripts were present in all three subgroups except for AhMYB177.

Chromosomal localization showed that the A and B subgenomes contained 99 and 97 R2R3-MYB genes, respectively (Figure 3). This suggested that the distributions of R2R3-MYB genes between the two subgenomes were almost not biased. In addition, we found that the distribution of genes on the corresponding chromosomes was similar between A and B subgenomes, except for Chr07 and Chr17, Chr08 and Chr18. This was probably the result of the complex rearrangement event on chromosomes 7, and 8 of two diploid wild ancestors, and subsequent retention of this rearrangement after polyploidization (Bertioli et al., 2016; Bertioli et al., 2019). The R2R3-MYB genes were unevenly distributed among the 20 chromosomes. Chr03 (19), Chr08 (18), Chr13 (18), and Chr18 (15) contained a larger number of R2R3-MYB genes. Chr02 (5), Chr07 (4), Chr10 (4), Chr17 (4), and Chr20 (4) possessed fewer R2R3-MYB genes. The density of R2R3-MYB genes on Chr08 was significantly higher than that on the other 19 chromosomes.

The R2R3-MYB genes were far more abundant in A. hypogaea than in lower terrestrial plants, suggesting that a large-scale gene duplication event occurred during the evolution of the plants (Du et al., 2015). To explore the mechanism of R2R3-MYB gene expansion in the cultivated peanut, we further analyzed the syntenic relationships between the peanut R2R3-MYB genes. Based on synteny analyses, 45 genes out of the 33 syntenic pairs in the A subgenome underwent segmental duplication events, while 39 genes out of the 24 syntenic pairs in the B subgenome were undergoing segmental duplication events (Figure 3 and Supplementary Table 5). 12 and eight genes in the A and B subgenomes, respectively, experienced tandem duplication episodes. In addition, significant numbers of orthologous genes were found between the A, and B subgenomes (Figure 4 and Supplementary Table 5). Based on synteny analyses, 173 of the 196 AhR2R3-MYBs had syntenic relationships between the two subgenomes. We found that some subgroups were expanded mainly by segment duplication, such as subgroups S11, NS-2, NS-4, SAtM85, S20a, S20b, S13, S1, S10&24, SAtM88, SAtM46, SAtM59, and NS-5, while tandem duplication occurred mainly in subgroups S5, NS-6, S20b, and S18 (Figure 3).

To further understand the evolution of R2R3-MYB duplication in peanut, we analyzed the rate of synonymous (Ks) and nonsynonymous substitutions (Ka) in gene duplication pairs (Supplementary Table 5). The Ks values for direct homologous gene pairs between A and B subgenomes range from 0 to 2.47. The frequency distribution peaks at Ks = 0.03, indicating a massive duplication event of R2R3-MYB genes 1.85 million years ago (Mya). Segment duplication and tandem duplication may occur in 147.75-35.91 and 163.87-2.20 Mya, respectively. Ka/Ks ratios for paralogous and orthologous ranged from 0-2.10 with an average of 0.27, whereas the ratios for tandem duplication ranged from 0.13-1.02 with an average of 0.36. The Ka/Ks analysis showed that the orthologous gene pair AhMYB052-AhMYB158 was neutrally selected and the tandem duplication pair AhMYB173-Arahy.QI53CA and the orthologous gene pairs AhMYB010-AhMYB108, Arahy.M1LASL-AhMYB118, and AhMYB031-AhMYB131 were subjected to positive selection, and all other synteny and tandem duplicated genes were subjected to purifying selection.

To further study the expression pattern of the R2R3-MYB genes in different tissues and explore its function in peanut growth and development, the tissue expression profiles of the R2R3-MYB genes were analyzed by using the transcriptome data of 22 peanut tissues (AhMYB015 and AhMYB100 were not detectable in the dataset) (Figure 5 and Supplementary Table 6). The peanut R2R3-MYBs can be clustered into 10 groups according to their tissue expression pattern (Cluster 1-10). Most genes from same phylogenetic subgroup showed similar tissue expression patterns and were clustered into the same cluster, such as AhMYB023, AhMYB024, AhMYB074, AhMYB125, AhMYB126, AhMYB174, and AhMYB175 in NS-6, which highly expressed in fruit and pericarp, were cluster in cluster 3 (Figure 5). All the members of subgroup SAtM35, for example, were clustered in cluster 9 (Figure 5). However, a certain proportion of the AhMYBs from the same phylogenetic subgroup showed differential tissue expression, such as AhMYB028 and AhMYB087 from S5, AhMYB028 was highly expressed in the reproductive shoot, while AhMYB087 was mostly expressed in seed (Figure 5), this might be because genes belonging to the same phylogenetic subgroup exercise similar functions, but in different tissues responses to different developmental processes or different environmental stimuli ( (Dubos et al., 2010).

Interestingly, members of the eight peanut specific subgroups (NS-1 to 8) were distributed in all clusters except 2 and 4.The two NS-1 members (AhMYB032 and AhMYB132) were strong expressed in perianth; the two collinear gene pairs (AhMYB005 and AhMYB104, AhMYB044 and AhMYB146) from NS-2 subgroup were highly expressed in seed and shoot, respectively. AhMYB084 (NS-3) had significant expression both in leaf and shoot; the three members in NS-4 (AhMYB010, AhMYB037 and AhMYB139) expressed mainly in leaf; AhMYB082 and AhMYB111 in NS-5 showed higher expression level in nodule, but AhMYB150 was enriched in peg tip; most of the members in NS-6 were clustered in cluster 3 and were highly expressed in pericarp; members in NS-7 (AhMYB061 and AhMYB163) and NS-8 exhibited strong expression in seed. These subgroups contained totally 30 genes, 21 of which were expressed at high levels in reproductive organs, especially the 10 genes in NS-6 that had high transcript abundance in the early pod development stage after peg tip entry, suggesting an important role of these genes in peanut reproductive development.

The expression of homeologous gene pairs from A and B subgenomes was examined in various tissues and developmental phases (Supplementary Table 7 references from (Bertioli et al., 2019)). The total number of homeologous gene pairs with expression biased towards the A subgenome was similar to the previously reported it did not differ significantly from the number with expression biased towards the B subgenome (P = 0.16, binomial test; n = 47 and 42 for A and B, respectively) (Figure 6) (Chalhoub et al., 2014; Zhang et al., 2015). In nine tissues (lateral leaf, seeding leaf, vegetative shoot tip, reproductive shoot tip, perianth, gynoecium, pattee 6 seed, pattee 7 seed, and pattee 8 seed), of the homologous pairs, there were more A subgenome-highly expressed genes rather than B subgenome, whereas the other 10 tissues exhibited the reverse pattern. These differences were more significant in the four reproductive tissues (P < 0.05, binomial test). 14 homologous gene pairs had biased expression in just one tissue whereas 18 homologous gene pairs exhibited the same bias in several tissues (Supplementary Table 8). Interestingly, we identified four homologous gene pairs (AhMYB006 and AhMYB105, AhMYB011 and AhMYB107, AhMYB053 and AhMYB157, AhMYB074 and AhMYB174) showed opposite biased expression in specific tissues, such as homolog pair AhMYB074 and AhMYB174, AhMYB174 owned a higher expression level in lateral leaf {log2FoldChange (B.vs.A homeolog pair comparison):1.906081616}, while AhMYB074 had a stronger expression level in perianth {log2FoldChange (B.vs.A homeolog pair comparison): -2.204691792}.

We examined the transcript abundance of the R2R3-MYB gene in peanut seedlings after water logging treatment to further investigate the function of the genes in water logging response. 90 genes from 24 subgroups (NS-6, S10&24, S11, S13, S14, S2, S20a, S20b, S21, SAtM85, SAtM88, etc.) were found response to waterlogging. According to the chronological order of the respond genes, they were divided into five categories (Figure 7). 26 genes showed a strong tendency to be down-regulated following the water logging treatment. Nine genes in the transitional response group showed strong elevated after 6 hours of treatment. Three genes (AhMYB121, AhMYB008, and AhMYB089) were considerably up-regulated among the 12 early responded genes, whereas the others showed a tendency of down- and subsequently up-regulation. At the time point of 3 and 5 days after treatment, a total of 27 and 16 genes were up-regulated, respectively. In addition, the results of qRT-PCR showed that the expression patterns of the 10 selected genes under water logging treatment were generally consistent with the transcriptome results (Supplementary Figure 6) which verified the results based on transcriptome data analysis.

To determine the role of R2R3-MYB genes in peanut, we conducted candidate gene association analysis using 59 single-nucleotide polymorphisms in Ad and AiR2R3-MYBs identified from transcriptome data of 146 peanut varieties (Supplementary Table 9) and phenotype of total branch number (TBN), pod length (PL) and root-shoot ratio (RS ratio) variation collected in five environments (Supplementary Table 10). One polymorphic site [A03_122181829 (C/S/G)] was uncovered highly associated with TBN, PL and RS ratio variation in five environments (P < 0.01) (Supplementary Table 11 and Figure 8A), located in the third exon region of AdMYB03-18 (AdMYB03-18) (Figure 8B). A03_122181829 mainly formed three haplotypes (A03_122181829 (C/S/G)) in the associated population (Figure 8C). Analysis results indicated that TBN in haplotype G was notably higher than those in haplotype C, while PL in haplotype C was higher than PL in haplotype G and S, RS ratio in haplotype S was higher than RS ratio in haplotype C and G (Figure 8D).