The HMMER3 (Mistry et al., 2013)¹ was used to search FAR1 gene family domains based on the hidden Markov model (HMM) file (PF03101), which was downloaded from the Pfam (Finn et al., 2014) database² by the FAR1 DNA binding domain identity PF03101. Four Arachis genomes, including A. duranensis, A. ipaensis (Bertioli et al., 2016), A. monticola (Yin et al., 2018), and A. hypogaea (Chen et al., 2019), were downloaded from the relevant database exhibited in each previous report. The G. max genome version 2.1 was obtained from the NCBI database (assembly accession: GCF_000004515.5). The key parameters were set as default and the cutoff value was set as 1e-15. The identified FAR1 members were confirmed using NCBI-Conserved Domain Database (CDD),³ SMART,⁴ and Pfam (Finn et al., 2014) (see text footnote 2) databases. All incorrect, repetitive, and non-FAR1 family members were removed.

The protein sequences of FAR1 family members were used to create multiple alignments using ClustalW with default parameter sets. The phylogenetic tree was constructed using MEGA version 7.0 under the neighbor-joining (NJ) method with 1000 bootstrap replicates (Kumar et al., 2016).

According to physical position, all the FAR1 genes were mapped to each reference genome using TBtools version 1.068 (Chen et al., 2020). Gene duplication events were analyzed using MCScanX with default parameters (Wang et al., 2012). The syntonic map was constructed using Python jcvi utility libraries.⁵ Non-synonymous (Ka) and synonymous (Ks) substitution of duplicated FAR1 genes were calculated by using KaKs_Calculator 2.0 (Wang et al., 2010) and plotted using R ggplot2 (Wickham, 2016).

The FAR1 gene structures were analyzed based on reference genome annotation information using our in-house Perl scripts. The conserved domain motifs of FAR1 members were identified by MEME version 4.12.0 (Bailey et al., 2015) with the number of motifs set to 10 and minimum and the maximum length of motifs set to six and 100, respectively (Bailey et al., 2009). The results of the gene structures, conserved domains and phylogenetic tree were plotted by using TBtools version 1.068 (Chen et al., 2020). The in silico elements of the AhFAR1 promoter in 1500 bp regions were searched using the PlantCARE database (Lescot et al., 2002).⁶ The original results were filted and “light responsive” elements were retained. The features of in silico elements were visualized by GSDS 2.0 (Hu et al., 2015).⁷

The transcriptomes of five different typical tissues, including flower, peg, leaf, root, and stem, were obtained from our previous study to analyses the expression levels of the members of the AhFAR1 gene family (Chen et al., 2019). Moreover, the transcriptomes of 11 development stages of shell and seed were downloaded from our previous report to evaluate the AhFAR1 gene expression patterns in underground tissues (Chen et al., 2016b). The expression levels of the AhFAR1 were evaluated using the fragment per kilobase million (FPKM) method. The expression levels were visualized using FPKM standardization data in R pheatmap packages.⁸

Six different typical components of peanut, including root, stem, leaf, peg, and pod, were used to verify AhFAR1 gene expression levels. In addition, five different developmental stages of ovule-carrying peg (peg length = 1–5 cm with 1 cm as step) and pods (including aerial, not swelling subterranean, early swelling subterranean, swelling subterranean, and mature pods) were used to confirm the expression levels of candidate AhFAR1 genes. Total RNA was extracted from each sample using a Plant RNeasy Mini Kit (TIANGEN, Beijing, China). The quality of the RNA extracted from each sample was checked using a Nano Drop (Thermo Scientific, United States). The qRT-PCR was carried out as described by Liu et al. (2019a). DNA-free RNA was used to synthesize the first strand of cDNA. The yellow leaf specific 8 (yls8) gene (Forward Primer: 5′-AACTGCTTAGCTGCTATTACCC-3′, Reverse Primer: 5′-TCGCCAAATAACACGTTGCATT-3′) was used as an internal control. Each measurement was carried out in three experimental replications, and each reaction was performed in biological triplicate replications. The relative expression levels of each target gene was analyzed using 2^{–△ △ CT} method.

A total of 650 FAR1 genes were identified from four Aarchis and its closely related species G. max (Supplementary Table 1). Of them, the most (246) and least (36) FAR1 genes were retrieved from A. hypogaea and G. max, respectively. The total number of FAR1 genes identified from each genome exhibited great consistency with the evolutionary genetic relationships (Supplementary Figure 1). This result suggested that the FAR1 gene family was expanding rapidly in the Arachis species, which was consistent with our previous findings (Chen et al., 2016a; Lu et al., 2018; Yan et al., 2020). All of the identified FAR1 genes were mapped on the chromosomes of their respective genomes based on the physical positions of the genes (Figure 1A and Supplementary Figure 2).

The phylogenetic analysis revealed that the FAR1 domains could be roughly divided into Arachis and G. max groups. However, these domains identified in Arachis species could not be clearly clustered by different species attributes (Figure 1B). Moreover, the FAR1 identified in different species could be classified into different corresponding groups (Supplementary Figure 3). Deep phylogenetic analysis of the 246 AhFAR1 domains identified in A. hypogaea indicated that these domains could be classified into three main groups. Among them, 99, 11, and 136 members belong to group I, II, and III, respectively (Figure 1C). The protein sequence analysis indicated that the conserved motifs of the domains in each group were different, although they all contained the typical FAR1 conserved motifs, for example the N-terminus contained FYxxY, GFxxR, and CxxxG and the C-terminus had HNHxL (Figure 1D).

The gene structure (exon-intron-UTR), conserved motif organization and phylogenetic tree of all the FAR genes identified in each species were analyzed to gain insight into the difference of the FAR family in different plant species (Figure 2 and Supplementary Figure 4). The results showed that members belonging to the same group had almost similar motifs and gene structures. For example, in A. duranensis, the members of the red group typically contained motif 3 and motif 8. The special motif 10 gene was observed in blue group (Supplementary Figure 4). In A. hypogaea, the distribution of exon-intron-UTR structure and the motif phase corresponded with the clusters of AhFAR1 genes (Figures 1C, 2), especially for group I (yellow) and III (blue), most of which contains motif 1 and motif 2. However, for group II (red) the genes did not contain these two motifs, but contained motif 3 and motif 4 instead (Figure 2). Overall, the conserved motif organizations, similarity of gene structures and phylogenetic tree results could be confidently verified, confirming the reliability of the classifications and evolutions of the FAR1 gene family in different leguminous plant species.

In plants, cis-acting elements of promoter regions often affect the biological functions of genes. Here, in order to identify the cis-acting elements of the 246 AhFAR1 genes obtained from A. hypogaea genome, a 1.5 Kb upstream genomic region of each gene was explored in the PlantCARE database (Lescot et al., 2002). A total of 2,147 cis-acting elements related to “light responsiveness” were identified (Supplementary Table 2 and Supplementary Figure 5). The average number of cis-acting elements identified by each gene was 8.7 (from 2 to 18). In Ahy_A03g015787 and Ahy_B02g058113, 18 cis-acting elements were identified (Supplementary Table 3). Moreover, a total of 31 types of cis-acting elements were obtained. Among them, the Box-4 element was the most widely distributed (601), and the ACA-motif, Pc-CMA2a, CAG-motif and P-box were each only identified once (Supplementary Table 4). Functional annotation indicated that all of the obtained cis-acting elements, such as Box-4, G-box and GT1-motif, were identified to be involved in the photopic response system (Supplementary Table 2). All of these cis-acting elements were functional annotated as light responsiveness through different regulation patterns, such as conserved DNA module, MYB binding site, gapA-CMA1, and gibberellin pathway (Supplementary Table 2). The results regarding of the cis-acting elements suggested that the AhFAR1 gene family could play critical roles in peanut growth and development.

The synonymous (Ks) and non-synonymous (Ka) substitution rates of paralogous and orthologous gene pairs were calculated to study selection pressures of FAR1 genes in Aarchis and its related species G. max in the process of biological evolution. Each of the 13 paralogous gene pairs were identified in subgenome A (At) and subgenome B (Bt), respectively. Moreover, a total of 33 paralogous gene pairs were observed between At and Bt (Supplementary Table 5). There was no significant difference in Ka/Ks value either within or between the two subgenomes (Figure 3A), suggesting that the paralogous AhFAR1 genes had undergone a relaxed selection pressure. Furthermore, a total of 11, 6, 21, 40, and 62 FAR1 paralogous gene pairs were identified in G. max, A. ipaensi, A. duranensiss, A. monticola, and A. hypogaea, respectively (Supplementary Table 6). A large number of paralogous gene pairs were identified in the allotetraploid A. hypogaea, which might be caused by polyploidization or whole genome duplication events during evolution. The average Ka/Ks of paralogous FAR1 in G. max (0.30) was significantly smaller than that calculated for A. duranensiss (0.51), A. monticola (0.42), and A. hypogaea (0.49). However, there was no significant difference in the average Ka/Ks value among the four Aarchis species (Figure 3B). These results indicated that the FAR1 gene family experienced strong evolutionary selection pressure when evolution occurred from Leguminosae, while these genes were under relaxed selection in the subsequent evolution of Arachis species. The average Ka/Ks of orthologs in both G. max vs. A. hypogaea and G. max vs. A. monticola were all smaller than that of the Arachis combination gene pairs, suggesting that the biological functions of the FAR1 gene family underwent a relaxed diversification after the divergence in Arachis species (Figure 3C and Supplementary Table 7).

According to the collinearity analyzes among G. max and Arachis species, only 13 FAR1 gene pairs showed a syntenic relationship in both G. max vs. A. duranensis (8) and vs. A. ipaenesis (5), which was significantly less than that identified among Arachis plants (50) (Figure 3D and Supplementary Table 8). This result suggested that the FAR1 gene family was substantially expanded during polyloidization, evolution, and domestication. In addition, many more FAR1 homologous gene pairs were identified among A. monticola vs. A. duranensis (8) and vs. A. ipaenesis (25) than that of among A. hypogaea vs. A. duranensis (2) and vs. A. ipaenesis (8) (Supplementary Table 8 and Supplementary Figure 6). This finding might support the conclusion that A. hypogaea was domesticated from an intermediate species, the wild tetraploid A. monticola, which was formed from its two diploid progenitors, A. duranensis and A. ipaenesis (Yin et al., 2019). Interestingly, some FAR1 genes were associated with at least three or more syntenic gene pairs, such as Araip.M4D5S displaying collinearity with GLYMA_13G341600, EVM0033853, and Aradu.8P3KF in different plants (Figure 3D and Supplementary Table 8), indicating that these genes play key roles in the FAR1 gene family and were highly conserved during evolution.

In order to investigate the expression patterns of AhFAR1 in different tissues (flower, peg, leaf, root, and stem) and in the different development stages of the pod (11 stages of shell and seed), the transcriptome data obtained from our previous study (Chen et al., 2016b,2019) was re-interrogated. The results revealed that a total of 157 AhFAR1 genes were not expressed in any detected tissue, which was speculated to be because these genes had other special expression patterns and could not be detected in our previous study (Supplementary Table 9). However, as showed in Figure 4A, a total of 89 AhFAR1 genes showed tissue-specific expression in peanut. For example, 12 genes in flower, 18 genes in pod, 15 genes in leaf, 31 genes in leaf and root, and 13 genes in stem showed the highest expression levels of the different tissues. In depth analysis of the expression patterns of the 246 AhFAR1 genes in 11 developmental stages of shell and seed showed that 84 of them were expressed in at least one development stage (Supplementary Table 10), and some of the 84 genes exhibited stage-specific expression in shell or seed (Figure 4B). Moreover, a lot of genes were highly expressed in whole-seed developmental stages but expressed at low levels in the shell (Figure 4B). These findings indicated that the AhFAR1 genes might be involved in shell or seed development in peanut. Unequal expression of the AhFAR1 genes between the two subgenomes (At > Bt) was observed (Figure 4C) indicating a divergence of gene biological function compared to in other polyploid plants (Wu et al., 2018). GO term enrichment analysis revealed that the AhFAR1 genes are mainly involved ion, nucleic acid, and compound binding (Figure 4D).

Previous reports indicated that the FAR1 gene family is involved in the light response in plants (Casal, 2000; Wang and Deng, 2003; Huang et al., 2012), and expanded in A. hypogaea, A. duranensis, and A. ipaensis (Chen et al., 2016a; Lu et al., 2018; Yan et al., 2020). From these findings, we hypothesized that the FAR1 genes might be related to the peanut’s unique pod development pattern, such as its aerial flower and subterranean fruit. After pollination, the ovule-carrying peg (P0 stage) forms and elongates with peg-bending growth to bury the fertilized ovule into the soil (P1 stage), and then the pod develops continuously until maturation (P2–P10 stages) (Figure 4E). In this study, seven AhFAR1 members were specifically expressed in peg and pod (shell or seed) as identified through the transcriptome analyses (Supplementary Figure 7 and Figures 4A,B). The quantitative real-time PCR (qRT-PCR) in A. hypogaea cv. Fuhuasheng confirmed that three of the seven AhFAR1 genes, including Ahy_A10g049543, Ahy_A06g026579, and Ahy_A10g048401, had tissue-specific expression, especially in the pod and peg (Figure 4F). Moreover, examination of the five different development stages of ovule-carrying peg revealed that Ahy_A10g049543 was highly expressed in middle and later stages but the other two genes possessed the opposite pattern (Figure 4G). In addition, the relative expression in the different development stages of the pod indicated that Ahy_A10g049543 was highly expressed in the early stage of the pod, Ahy_A06g026579 was mainly expressed in mature stage of the pod and Ahy_A10g048401 was expressed in middle and late stages of the pod, such as early swelling, swelling, and mature stages (Figure 4H). The above results suggested that the AhFAR1 genes showed specific expression responses in different tissues or development stages and might play an important role in regulating pod development in peanut. The biological function verification of the three candidate AhFAR1 genes should not be ignored in future.