Gene family and phylogenetic analysis

A total of 16,791 orthologous gene groups were identified among 18 species using OrthoMCL (Li et al., 2003; Supplementary File 4: Data S3), including 959 A. ipaensis-specific families containing 6,443 genes (Supplementary File 1: Table S14). A gene ontology (GO) annotation suggested differentially enriched functional categories in the peanut-specific families, indicating that these gene families may be closely related to the unique Arachis growth characteristics, such as aerial flowers but subterranean fruit, and lipid biosynthesis (Supplementary File 1: Figures S13–S15). Moreover, 1,624 of these orthologous groups were single-copy orthologs (Supplementary File 1: Table S15 and Figure S16). In addition, 6,443 unique paralogs of A. ipaensis genes occurred in species-specific homolog groups, indicating that these groups could arise from genomic structural rearrangements that obscured simple orthology, such as nonallelic recombination or gene conversion, followed by duplication (Supplementary File 1: Table S15 and Figure S16; Varshney et al., 2013). We identified 12,017 orthologous groups common to all five leguminous species (Figure ​(Figure2B),2B), 11,985 groups between A. ipaensis and Ricinus communis (oilseed crop) (Supplementary File 1: Figure S17), 9,099 groups between A. ipaensis and Gramineae/Poaceae crops (Supplementary File 1: Figure S18), and 10,501 orthologous groups are common to A. ipaensis and other distantly related plant species (Supplementary File 1: Figure S19). These results provide an important molecular foundation for comparative biology and for functional mechanistic inferences in A. ipaensis, as well as other species, because simple orthologous family genes often exhibit conserved molecular functions that were maintained during evolution process.

A polygenetic tree based on single-copy orthologous genes showed A. ipaensis and A. duranensis in the same clade, which did not include any other leguminous species, indicating their closer genetic distance and divergence time (Figure ​(Figure2C).2C). Furthermore, a special phylogenetic tree estimated that the divergence of the two species occurred ~2.9 million years ago (Mya) (Figure ​(Figure2D),2D), which was fairly consistent with a previous report (~2.16 Mya) (Bertioli et al., 2016). Syntenic blocks identified between A. ipaensis and other species was found to be extensively conserved (Supplementary File 1: Table S16). The largest number of syntenic blocks was identified between A. ipaensis and G. max (Figure ​(Figure2E).2E). The longest syntenic block (>10 Kb) was observed between A. ipaensis and A. duranensis (Supplementary File 1: Table S16). The numbers of syntenic blocks identified within the respective A. ipaensis and A. duranensis genomes were extremely lower than the number between the two genomes (Supplementary File 1: Figure S20) as well as the number between the A. ipaensis and G. max genomes (Supplementary File 1: Figure S21; Bertioli et al., 2016), indicating that few large-scale genome duplication events occurred in the A. ipaensis genome's evolution or that syntenic blocks were lost after genome duplication events.

The Ks values between paralogous or orthologous genes reveals a mechanism of molecular evolution (Lna, 1996). Distributions of Ks distances between paralogs within A. ipaensis and orthologs among A. ipaensis, leguminous crops and other species were plotted (Figure ​(Figure2F).2F). The A. ipaensis paralogs showed a peak at ~ 0.80, which is similar to those of M. truncatula (~0.80) and L. japonicus (~0.73) (Cannon et al., 2006) but lower than those of A. duranensis (~0.9) and A. ipaensis (~0.95) (Chen et al., 2016). Thus, the whole-genome duplication events of A. duranensis and A. ipaensis occurred around the time that corresponds to a Ks value range of 0.8–0.95. In addition, A. duranensis and A. ipaensis orthologs showed a prominent peak at ~0.04, which is consistent with a previous study (Bertioli et al., 2016). Assuming a synonymous substitution rate per synonymous site of 6.1 × 10⁻⁹ per year for eudicots (Lynch and Conery, 2000), the two species were estimated to have diverged ~3.28 Mya, which is close to the estimation based on the phylogenetic tree (Figure ​(Figure2D).2D). Furthermore, Ks dating suggested the divergence of A. ipaensis and G. max (Ks = ~0.54) at 44.3 Mya and that of A. ipaensis and C. arietinum (Ks = ~0.64) at 52.5 Mya.

The graphic trend of the Ka/Ks (ω) and Ks between the orthologs of A. duranensis and A. ipaensis formed three clusters, such as Ks = 0–0.3, 0.5–1.5, and >1.5, and the ω values decreased as the Ks values increased (Supplementary File 1: Figure S22). The genes with Ks ≥ 1.5 are attributed to pan-eudicot palaeoploidization, and the genes with lower ω ratios are considered to be under neutral selection. Here, the 45 A. ipaensis genes with ω ratios > 1 may be under positive selection pressure (Supplementary File 1: Figure S23).

Peanut is an allotetraploid species that may have originated from a single resent hybridization event between two diploid species, followed by polyploidization. Cytogenetic, phylogeographic and molecular evidence indicates that A. duranensis and A. ipaensis are the most likely donors of the A and B subgenomes, respectively (Kochert et al., 1996; Seijo et al., 2007; Robledo et al., 2009; Robledo and Seijo, 2010; Moretzsohn et al., 2013). A previous study estimated the divergence of the two species at ~2.88 Mya (Moretzsohn et al., 2013). The estimation using a comparative genomics analyses between them was ~2.9 Mya, which was fairly consistent with our report. Moreover, sequence comparisons with tetraploid cultivated peanut estimated the divergence times of A. duranensis and A. ipaensis from the A and B subgenomes of A. hypogaea as ~247,000 and ~9,400 years, respectively (Bertioli et al., 2016).

Comparative genomics analyses of chromosomal structure and synteny between A. duranensis and A. ipaensis indicated that some chromosomes shared a conservative linear structure that was partially in accordance with our other analyses (Supplementary File 1: Figure S20). Other analyses showed a large inversion in one or both arms of a chromosome (Bertioli et al., 2016). In contrast, chromosomes 07 and 08 have undergone complex rearrangements that were consistent with cytogenetic observations (Seijo et al., 2007; Nielen et al., 2010). Importantly, a genomic comparison showed a fundamentally one-to-one correspondence between the diploid chromosomes and cultivated peanut linkage groups. However, the A. duranensis chromosomes were less similar to A. hypogaea sequences than those of A. ipaensis (Bertioli et al., 2016). These results may help to uncover potential mechanisms of hybridization events in the future.

Disease resistances and nucleotide-binding site (NBS)-leucine-rich repeat (LRR) encoding genes

Plant NBS-LRR proteins encoded by resistance genes (R genes) play key roles in the responses of plants to various pathogens. The R genes can be classified into various subfamilies based on the present of different domain, such as CC-NB-LRR, TIR-NB-LRR, ser/thr-LRR, Kin-LRR, and others (e.g., Mol and Asc-1; Sanseverino et al., 2010). The A. ipaensis genomic assembly contains 1,437 putative disease R genes as assessed by a screening of the PRG database (Supplementary File 1: Table S17; Supplementary File 5: Data S4; Sanseverino et al., 2010). Compared with other legumes, the A. ipaensis genome possesses more R genes than the G. max and M. truncatula genomes but less than the A. duranensis and C. cajan genomes. Moreover, these R genes tend to cluster on different scaffolds. For example, several large clusters containing 6–10 R genes are located on six different scaffolds (Supplementary File 1: Figure S24). The NL subfamily of genes, which confers resistance against pests and diseases, is the second largest R gene-containing family, and these genes can be clustered reasonably into different individual clades in A. ipaensis, A. duranensis, and A. thaliana, indicating that gene divergence occurred during the evolution of the three species (Supplementary File 1: Figure S25). In addition, we analyzed protein motifs in the most homology of the top 20 R genes found in PRG database using MEME (Bailey et al., 2009), and the results showed highly conserved LRR motifs (Supplementary File 1: Figure S26). However, further investigation is required to determine the biological functions of these R genes.

Identification of genes related to biological nitrogen fixation

Nitrogen is one of the most important plants require nutrients, and in agriculture nitrogen availability influences both crop yield and quality. Leguminous plants, such as peanut, soybean, and Medicago, can transform molecular nitrogen into available ammonia nitrogen through the leguminous-root-nodule bacteria nitrogen-fixing system that results from the symbiotic interactions between leguminous plants and rhizobia (Figure ​(Figure3A).3A). In the A. ipaensis and A. duranensis genomic assemblies, 16 and 38 root-nodule developmental genes respectively, have been identified (Supplementary File 1: Table S18; Supplementary File 6: Data S5). As expected, there are greater numbers of nodulation-related genes present in leguminous plants than in non-leguminous plants, such as A. thaliana, O. sativa, and Z. mays (Supplementary File 1: Figure S27).

Open in a separate window

Figure 3

Biological nitrogen fixation in leguminous plants. (A) Genes involved in nodule initiation, development and signal recognition pathway. (B) Protein sequence alignment of Nod related genes identified in A. ipaensis and A. duranensis. (C) Phylogenetic tree of nodule development genes and their homologs from A. ipaensis and A. duranensis. (D) Identification of high conserved domains of leucine-rich repeat (LRR) receptor kinases. Red dashed boxes represent LRR conserved motif. (A) The rhizobium (blue) attach to the surface of root hair cell. After swelling, deformation, curling and infection thread, the bacteria are released into cells via endocytosis then a vacuole-like structures (symbiosomes), in which the bacteria convert N₂ to NH₃, formed. But how is the Nod signal transmitted? Initially, the rhizobia-derived signal is perceived by LysM-type protein receptor kinases, such as NRF1 and 5 (Radutoiu et al., 2003) and SYM10 (Schneider et al., 2002) identified in L. japonicus and P. sativuml, followed by a downstream leucine-rich receptor kinase, for example SYMRK (Stracke et al., 2002 and Capoen et al., 2005), NORK (Endre et al., 2002), DMI2 (Catoira et al., 2000), and SYM19 (Stracke et al., 2004) from L. japonicus, Sesbania rostrata, M. sativa, M. truncatula, and P. sativuml, respectively. Then, the Nod factor (NF) signal is processed through a signal transduction cascade involving proteins including ion channels [MDI1(Ané et al., 2004), CASTOR (Imaizumi-Anraku et al., 2005), POLLUX (Imaizumi-Anraku et al., 2005), and SYM8 (Edwards et al., 2007)], calcium-calmodulin-dependent kinase (CaCaDK) (MDI3 and SYM9) (Lévy et al., 2004) and transcription factors [NSP1 (Smit et al., 2005), NSP2 (Kaló et al., 2005), SYM7 (Kaló et al., 2005), NIN (Schauser et al., 1999), and SYM35 (Borisov et al., 2003)]. Finally, rhizobia infection occurred primarily through uncharacterized target genes that may be activated by these TFs.

Nitrogen-fixing root nodules are important symbiotic organs that provide an epiphytic site for rhizobia to convert atmospheric nitrogen to ammonia, and supply its host plant with fixed nitrogen. In return, the rhizobia gain photosynthates from the plant (Figure ​(Figure3A).3A). In leguminous plants, multiple genes are involved in the formation and development of root nodules, as well as in the autoregulation of the nodulation (AON) process, which is a systemic feedback loop used to avoid an excessive energy expenditure from “unwanted” nodulation (Figure ​(Figure3A;3A; Supplementary File 7: Data S6). Here, four homologous LRR receptor kinase genes were identified in A. ipaensis (XP_004512550.1-D2 and XP_007158329.1-D2) and A. duranensis (XP_015956675.1 and XP_015963325.1) (Figure ​(Figure3B;3B; Supplementary File 1: Figure S28; Supplementary File 7: Data S6). A phylogenetic tree showed that the four homologous genes were clustered into an independent clade, together with other LRR receptor kinase genes (Figure ​(Figure3C).3C). Interestingly, these four genes contain multiple common motifs, including a conserved LRR motif, indicating a similar biological function (Figure ​(Figure3D).3D). The GO analyses suggested that the four homologous genes are involved in ion binding and signal transducer activity (Supplementary File 1: Figures S29–S32). More importantly, the proteins encoded by the four genes showed similar three-dimensional structures and localized on the cell membrane (Supplementary File 1: Figures S29–S32).

We also identified two other nodule development-related genes (XP_015934647.1 and XP_015939255.1) that are homologous to the TF genes of the GRAS family in A. duranensis. One gene is homologous with MtNSP2 and PsSYM7 (Kaló et al., 2005), while the other is an ortholog of MtNSP1 (Imaizumi-Anraku et al., 2005) (Figures 3A,B; Supplementary File 1: Figure S33; Supplementary File 7: Data S6). The phylogenetic tree indicated that the two homologs were classified into the TF category but appeared in different branches (Figure ​(Figure3C).3C). In addition, the GO enrichment indicated that the two genes participate in the regulation of multiple biological processes, such as nucleic acid-binding TF and signal transducer activities (Supplementary File 1: Figures S34, S35). The three-dimensional structures of the two proteins were completely dissimilarity, and the two proteins localized in the nucleus (Supplementary File 1: Figures S34, S35). These results could provide candidate genes and basic bioinformation for further functional studies of nodule formation in leguminous crops.

Genetic mechanism of drought adaptation

Peanut (A. hypogaea L.) is a typical upland crop in tropical, subtropical, and warm temperate climates. Drought adaptation plays a central role in their growth and development. During drought stress, TFs, such as MYB, MYB-related, NAC, WRKY, bZIP, and ERF, are involved in numerous physiological responses (Shinozaki and Yamaguchi-Shinozaki, 2000) (Supplementary File 1: Figure S36). Here, the total number of TF genes identified in upland crops was greater than that found in hygrophilous plants (Supplementary File 1: Table S19 and Figure S37). Notably, in A. ipaensis we identified 185 MYB and 129 MYB-related TFs (Supplementary File 1: Table S19), most of which contain a highly conserved DNA-binding domain, and they are key factors in regulatory networks controlling development, metabolism and responses to biotic and abiotic stresses (Dubos et al., 2010). The second large number of drought tolerance-related TFs, with 170 members, is the ERF family (Supplementary File 1: Table S19). The ERF proteins, sharing a conserved 58–59 amino-acid domain, are key regulators linked to responses to plant stresses, such as cold, drought and pathogen attack (Supplementary File 1: Figure S38; Singh et al., 2002). In A. duranensis, A. ipaensis and A. hypogaea species, sets of 51, 57 and 53 ERF TF family proteins, respectively, were obtained from the Plant Transcription Factor Database (Jin et al., 2015, 2017; Supplementary File 1: Figures S39–S41). These TF proteins contained different DNA-binding domains and can be categorized into different branches based on different motif permutation structures, indicating the distinct functional and evolutionary features of ERF TFs in different Arachis species (Supplementary File 1: Figures S39–S41).

Heat-shock proteins (Hsps)/chaperones are important defense mechanism members against abiotic stresses, such as drought, salinity and extreme temperatures (Wang et al., 2004; Supplementary File 1: Figure S42). Drought stress is a common factors that induces Hsp expression (Kimpel et al., 1990; Sun et al., 2002). To elucidate the cause of drought tolerance, five major families of Hsps/chaperones were identified in upland crops and hygrophilous plants (Supplementary File 1: Table S20 and Figure S43). As expected, the total number of Hsps/chaperones obtained in upland crops was much great than in hygrophilous plants (Supplementary File 1: Figure S43).In particular, A. ipaensis and G. soja had 118 and 34 Hsp70 subfamilies, respectively, compared with only 1 in rice (Supplementary File 1: Table S20). The great number of Hsps/chaperones detected in A. ipaensis and G. soja indicates the nature of drought adaptation in upland crops.

The subtilisin-like protease (SDD1) gene family is involved in the regulation of stomatal density and distribution to adjust for drought stress by modulating the apertures of these pores flanked by two guard cells (Berger and Altmann, 2000). In the expanded gene families, 39 and 40 SDD1 genes were identified in A. ipaensis and A. duranensis, respectively (Supplementary File 8: Data S7). These gene families were divided into different clusters according to their related functions but showed a pattern of cross-distribution in each cluster based on their different genetic relationships (Supplementary File 1: Figure S44).

Whole-genome shotgun sequencing and de novo assembly

Whole-genome shotgun sequencing was performed under the HiSeq2500 Sequencing System with 11 paired-end sequencing libraries, including 250, 500, 800 bp, 2, 5, 10, and 20 Kb using the standard protocol provided by Illumina (San Diego, USA).

SOAPdenovo2 (version 2.04.4) (Luo et al., 2012) was employed with optimized parameters to construct contigs and original scaffolds as previous described (Chen et al., 2016). Subsequently, SSPACE (version 2.0) (Boetzer et al., 2011) was used to link the scaffolds constructed by the SOAPdenovo2 as previous described (Chen et al., 2016).

The genome size was estimated based on the 17 K-mer distribution using the total length of sequence reads divided by sequencing depth, and the frequency of each 17-mer were calculated from the whole genome sequenced reads to evaluate the sequencing depth. Subsequently, the A. ipaensis genome size was calculated by following the algorithm: Genome size = K-mer number/Peak depth (Bertioli et al., 2016).

The gene coverage of the assembled genome was comprehensively evaluated using available public transcript sequence tags or expressed sequence tags. Core eukaryotic genes identified by CEGMA v.2.3 (Parra et al., 2007) were remapped to the A. ipaensis genome assembly by BLAT (Kent, 2002) to evaluate the quality of the assembly. CEGMA data were downloaded from the Korf Lab research group at the Genome Center, UC Davis (http://korflab.ucdavis.edu/datasets/cegma/#SCT6).

Gene prediction and function annotation

To annotate the A. ipaensis genome, an automated genome annotation pipeline MAKER was performed to produce de novo gene prediction, infer 5′ and 3′ UTR, and integrate these data to generate final downstream gene models with quality control statistics (Cantarel et al., 2008). All predicted genes were functionally annotated as previous described (Chen et al., 2016). The annotation was conducted using the BLAST+ (version 2.2.27) with 1e-5 as the E-value threshold to against the SwissProt and TrEMBL databases (Bairoch and Apweiler, 2000). To infer functions for the predicted genes, InterProScan (version 4.7) (Zdobnov and Apweiler, 2001) was used to search the predicted genes against the protein signature from InterPro with default parameters. All genes were also aligned against to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Kanehisa et al., 2004).

In order to evaluate the conservation of A. ipaensis ICG_8206 gene model, the BLASTP was used to query the A. ipaensis ICG_8206 proteome against the proteomes of other plant species (Supplementary File 1: Table S7) with an E value of 1e-10 as cut-off (Supplementary File 1: Table S8).

Gene family analysis

All the predicted gene models were analyzed using OrthoMCL (Li et al., 2003) to identify shared and specific gene families among 17 species (Supplementary File 1: Table S7). In the first step, inter and intra species BLASTP with an E-value cutoff of 1e-10 was performed to detect reciprocal best hit pairs between species (putative orthologs), as well as sets of genes within species (putative co-orthologs or in-paralogs). The reciprocal best hit matrix served as the basis for ortholog definition using OrthoMCL. Subsequently, orthologous groups were organized into species-specific and higher taxonomic level groups by requiring that at least one sequence from each clade under comparison be present in the intersecting set. Finally, based on fourfold degenerate sites of single-copy ortholog genes in all species, a phylogenetic tree was constructed using MEGA v6.0 (Tamura et al., 2013) and PhyML v3.0 (Guindon et al., 2010).

To identify TFs in A. ipaensis, the PlantTFDB database was used to search TFs in other plant species (http://planttfdb.cbi.pku.edu.cn/). The predicted genes were used to BLAST search against the PlantTFDB (E-value: 1e-10). The FAR1 motif was predicted using the Multiple Expectation Maximization for Motif Elicitation (MEME)/Motif Alignment and Search Tool (MAST) system (http://meme-suite.org/) (Bailey et al., 2009) and visualized using the TBtools (version 0.4999) (https://github.com/CJ-Chen/TBtools).

Non-coding RNAs and repetitive sequence annotation

Non-coding RNAs were predicted by aligned A. ipaensis genome assembly to against the Rfam databese (version 12.1) (Nawrocki et al., 2015). The pre-tRNAs were identified using tRNAscan-SE (Lowe and Eddy, 1997), pre-rRNAs were predicted using RNAmmer (Lagesen et al., 2007), pre-snRNAs were annotated using INFERNAL (Nawrocki et al., 2009) and others were also obtained by BLAST search against the Rfam database.

The RepeatMasker (Chen, 2004), RepeatProteinMask (http://repeatmasker.org/), Tendem Repeats Finder (TRF) (Benson, 1999) and RepeatModeler (Smith and Hubley, 2014) were performed to identify repetitive sequences through homolog and de novo prediction. The RepeatMasker and RepeatProteinMask were used to screen the A. ipaensis genome against the RepBase database (http://www.girinst.org/). The transposable elements (TEs) were classified as described without consideration of the gaps in the genome assembly (Wicker et al., 2007).

Identification of SSRs and SNPs

MIcroSAtellite (http://pgrc.ipk-gatersleben.de/misa/) was used to mine SSRs in A. ipaensis genome, and primer 3 v3.0 was used for primer design (Thiel et al., 2003; Untergasser et al., 2012). A SSR was defined with at least 6 repeats for di-nucleotide motifs or 4 repeats for tri-, tetra-, penta-, and hexa-nucleotide motifs. The maximum number of interrupting nucleotides in a compound SSR was set as 100.

Reads from five re-sequenced genotypes including two A-genome genotypes (ICG_8123 and ICG_8138) and three B-genome genotypes (ICG_8960, ICG_8209, and ICG_13160) were used to identify genome SNP and InDel variations (Chen et al., 2016). Total of these sequenced reads were aligned to the reference genome (ICG_8026) using the Burrows Wheeler Aligner program (BWA) (Li and Durbin, 2009). Subsequently, SNPs and InDels were identified using GATK v3.5 (http://www.broadinstitute.org/gatk) with default parameters, respectively.

Evolutionary and syntenic block analyses

The phylogenetic tree was constructed based on single-copy orthologous genes shared in A. ipaensis and other 17 plants (Supplementary File 1: Table S7) using MEGA v6.0 with the maximum-liklihood algorithm (Tamura et al., 2013).

Syntenic blocks between the genomes of A. ipaensis and other plants were identified using the MCScanX with default parameters (Wang et al., 2012) and visualized on the genome using Circos (Krzywinski et al., 2009). Genomic sequences were first aligned annotated genes based on amino acid sequence using Promer package of Mummer (version 3.22) (Delcher et al., 2002). Whole genome dot plots were generated using Mummerplot (Delcher et al., 2002) and Gunplot 5.0 (www.gnuplot.info/). Ks values of the homologs within collinearity blocks were calculated using the perl script, add_ka_and_ks_to_collinearity.pl included in MCScanX package, and the median of Ks values was considered to be a representative of the collinearity blocks.

We sincerely thank all the participants. In particular, we thank Jianan Zhang and Haofa Lan for helping us to analysis and upload the sequencing data. We thank Lesley Benyon, Ph.D., from Liwen Bianji, Edanz Group China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.