A F₆ RIL population including 242 lines was developed from a cross between homozygous cultivars “ZH16” and “sd-H1.” Field experiments including the parents and segregating populations in years 2015 (Wuchang, E 114° 34′/N 30° 59′) and 2016 (Wuchang, E 114° 34′/N 30° 59′ and Yangluo, E 114° 52′/N 30° 59′) were conducted at Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences, Hubei Province, China. The seeds were sown at the end of April of each experimental year. Each accession was planted in a single row, with 8 ± 10 plants in each row and 10-cm intervals between plants within each row, and 30-cm intervals between rows. The parental genotypes were also sown after every 50 rows as the controls. Randomized blocks were arranged, and average 10 individuals for each line in each block were selected for phenotype investigationfor PH, TBN, PL, PW, LWRP, SL, SW, LWRS, HPW, and HSW. The length and width of the two-seeded pods and seeds were measured using a parallel rule. The pods and seeds were also weighed on a digital scale. The PL/PW and SL/SW were calculated using PL and PW and SL and SW.

Broad-sense heritability was calculated using the equation h² = σ_g ²/(σ_g² +σ_ge²/n +σ_e²/nr). In the equation, σ_g², σ_ge², and σ_e ² were defined as the genetic variance, the interaction variance between genotype and environment, and the residual error variance, respectively. The n and r in the equation were defined as the number of environments and the number of replications, respectively. The general linear model (GLM) was used to estimate each variance component by SAS software. The Pearson's correlation coefficient was calculated for pair-wise comparison matrix of traits across the environments.

The genomic DNA was extracted using young leaves from each accession by a modified CTAB method. Through scanning the restriction endounclease sites on reference genome, the HaeIII endonuclease was chosen to digest the genomic DNA. The details of the SLAF-seq strategy and library construction were previously described (Sun et al., 2013; Zhang et al., 2015b). Finally, DNA fragments belonging to 400–530 bp in size were excised and diluted for paired-end sequencing which was performed on an Illumina HiSeq 2500 system (Illumina, Inc., San Diego, CA, USA). All SLAF-seq raw data from this study have been submitted to the NCBI SRA database (http://www.ncbi.nlm.nih.gov/sra/) under accession no. SRP108694.

The same procedure for read clustering and SNP genotyping was applied according to Sun et al. (2013) and Zhang et al. (2015b). After discarding the low-quality reads, the remaining reads were mapped to the reference genome by BWA software (Li and Durbin, 2009). The two progenitor genome sequences of tetraploid peanut were concatenated as the reference genome for read mapping (Bertioli et al., 2016). Only uniquely mapped reads were considered to identify SLAF markers (SLAFs). The SLAFs with less than four alleles were considered as polymorphic and potential markers. The Bayesian approach was used for genotype scoring, and the markers with >30% of missing data were further filtered out before genetic map construction, as described previously (Zhang et al., 2015b).

The procedure used for linkage map construction was the same as Zhang et al. (2015b), including ordering of markers by HighMap strategy (Zhang et al., 2015b), error correction by SMOOTH strategy (van Os et al., 2005), and imputation of missing genotypes by k-nearest neighbor algorithm (Huang et al., 2012). The multipoint method of maximum likelihood was applied for adding the skewed markers to the genetic map (Xu and Hu, 2009). The genetic distance between markers was calculated using the Kosambi mapping function, and also compared with the marker orders of each group obtained by MSTMap (Wu et al., 2008). For inconsistent regions of marker orders, we used more rigorous parameters (distance < 20 cM between adjacent markers and P-value < 10⁻⁵) to re-calculate and adjust them. A bin is defined as a position on the genetic map containing a clustering of markers showing a genetic distance of zero between each other.

The QTLs were detected using the composite interval mapping method in Windows QTL Cartographer 2.5 (Silva Lda et al., 2012). The permutation test was repeated 1000 times using 1.0 cM step and 5 control markers, and the LOD threshold (P < 0.05) were obtained for declaring significant QTLs. QTLs were considered as effective when the LOD scores were between 3.0 and the LOD threshold from the permutation test. Positive and negative additive effects mean that the favorable alleles were derived from parent “ZH16” and “sd-H1,” respectively. The QTL meta-analysis was performed to integrate the QTLs from different environments and traits using BioMeractor 4.2 software (Sosnowski et al., 2012).

According to the results of pilot experiment, the HaeIII enzyme was selected for SLAF library construction. Considering the low rate of DNA polymorphisms in peanut, we sequenced a high amount of data of the parents and the RIL population to obtain a sufficient number of markers. A total of 524.83 Gb of data containing 2,624.07 M paired-end reads of 100 bp in length were obtained. The Q30 ratio was 90.40% and GC content was 43.25% in average. Of these high-quality data (Supplementary Table S1), 12.60 Gb were from the male parent “sd-H1” with 62,997,078 reads, and 11.59 Gb were from the female parent “ZH16” with 57,928,805 reads. Read numbers for the F6 population ranged from 6,800,354 to 18,885,615 with an average of 10,343,589. The number of SLAFs in male and female parents was 737,780 and 732,877, corresponding to 49,660,470 and 45,483,812 reads, respectively. The average sequencing depth for each SLAF marker was 67.31- and 62.06-fold in male and female parents, respectively. In the F6 population (Supplementary Table S1), the number of SLAFs ranged from 489,875 to 680,491 with an average of 589,904, and the coverage ranged from 9.04 to 27.92-fold with an average of 13.75-fold.

After filtering repetitive SLAFs, 824,866 SLAFs were detected, and 7.6% of these were polymorphic (Table 1). According to the genotype encoding rule, we obtained the number of SLAFs corresponding to eight segregation patterns (ab × cd, ef × eg, hk × hk, lm × ll, nn × np, aa × bb, ab × cc, and cc × ab) (Supplementary Figure S1). 28,720 polymorphic SLAFs containing 31,526 SNPs belonging to aa × bb segregation pattern were used in the following construction of linkage map because the two parents used for the cross are homozygous lines with a genotype of aa or bb. Finally, the low-quality SLAFs were discarded when they are with a parental sequence depth of < 10 × , completeness < 70%, and significant segregation distortion (P < 0.001).

The final map included 3,630 markers belonging to 2,636 bins on the 20 linkage groups (LGs) (Figure 1 and Supplementary Figure S2) and was 2,098.14 cM in length with an average inter-marker distance of 0.58 cM. As shown in Table 2, B07 was the largest LG which has 429 markers, with the length of 144.64 cM and an average distance of 0.34 cM. On the contrary, A02 was the smallest LG which has only 43 markers, with the length of 77.58 cM and an average distance of 1.80 cM. The linkage degree between the markers was reflected by “Gap ≤ 5,” which ranged between 93.15 and 100%, and showed an average value of 98.34%. The largest gap on this map was 15.44 cM, which was located on A09, followed by 14.10 cM on A07.

There are three different types of markers in genetic map including 3,519 “SNP_only,” 101 “InDel_only,” and 10 “SNP&InDel” markers, accounting for 96.94, 2.78, and 0.28%, respectively (Supplementary Figure S3). Among the 3,519 markers of the “SNP_only” type, 95.28% had a single SNP locus and the others had two or three SNP loci. Two transition types of SNPs, including Y (T/C) and R (A/G), accounted for 36.05 and 35.05% of all markers, respectively. The other four transversion types of SNPs, including S (G/C), M (A/C), K (G/T), and W (A/T), ranged from 3.48 to 11.15% of all of the SNPs (Supplementary Figure S3). To confirm the authenticity of the identified SNPs, 19 out of 20 randomly selected SNPs were confirmed by Sanger sequencing (Supplementary Table S2).

Among the 3,630 markers, chi-square testing revealed that 688 (19%) showed significant segregation distortion, with 337 (49%) favoring the elite “ZH16” allele and 351 (51%) favoring the germplasm sd-H1 allele (Table 2). Segregation distortion regions (SDRs) were defined when at least four skewed markers were clustered. We found 20 SDRs including 688 distorted segregation markers on 7 chromosomes as follows: A02 (37), A06 (116), A08 (41), B02 (132), B06 (22), B09 (48), and B10 (248) (Supplementary Figure S4).

To assess the quality of this genetic map, a haplotype map was generated for each individual of the F₆ population (Supplementary Figure S5), which reflects the double crossover and the recombination events (West et al., 2006). We further conducted a comparison between the genetic and physical positions of the markers based on the reference genome. A high collinearity between them indicated that the markers were placed accurately within each LG (Figure 2). A large inversed segment on chromosome A05 (15–45 cM) was found, which was in opposite orientation between genetic and physical positions (Figure 2). Subsequently, the recombination rates along chromosomes were calculated by comparing the genetic distance to the physical distance (Mb) (Supplementary Figure S6). It varied among different chromosomes, ranging from 1.04 centimorgans per megabase (cM/Mb) for chromosome A03 to 4.75 cM/Mb for chromosome A08, and a genome-wide average recombination rate of 2.42 cM/Mb. In addition, we anchored 62 scaffolds accounting for 5.8 Mb of the total length to specific chromosomal positions using genetic markers on scaffolds, which allowed updating of the genome assembly (Supplementary Table S3).

Parent “ZH16” is a high-yield Spanish type cultivar with larger pods and seed size and two seeds in each pod, whereas parent “sd-H1” is a low-yield Valencia type cultivar with smaller pods and seeds and three seeds in each pod. The two parents showed significant differences (P < 0.001) in 14 yield-related traits, which were evaluated in 3 different environments in RILs (Table 3). The traits examined showed approximately the same phenotypic data distribution for all 3 environments (Table 3). A continuous distribution were shown for all traits in RIL population (Figure 3), indicating that these traits were controlled by multiple genes and belonged to a quantitative inheritance pattern.

We conducted ANOVA analysis for all 14 traits across the environments, and calculated the effects of genotype (G), environment (E), and genotype-environment interactions (G × E) on the traits (Table 4). They showed high broad-sense heritability, ranging from 65.7 to 86.2%, suggesting a major role of genetic factors in the expression of these traits as well as a considerable proportion of environmental variation affecting these traits. To elucidate their relationships, correlation coefficients (r) between the traits were calculated across the various environments (Figure 4). Plant height (PH) and lateral branch length (LBL) were strongly positively correlated, with a correlation coefficient of 0.97. Significant positive relationships also have been found for internode number (IN) and lateral branch length (LBL) (r = 0.70), Pod length (PL) and pod width (PW) (r = 0.76), PW and seed width (SW) (r = 0.80), hundred-pod weight (HPW) and hundred-seed weight (HSW) (r = 0.71). The significant phenotypic correlations among these traits coincided with the results of QTL co-localization (see next).

A total of 62 QTLs for all of the 14 traits were detected on 12 chromosomes across the 3 environments (“WuH15,” “WuH16,” and “YangL16”), and explained 4.03 to 18.9% of the observed phenotypic variation, respectively (Supplementary Table S4). Furthermore, 33 consensus QTLs were obtained by the trait-by-trait meta-analysis. For the five plant architecture traits, We have detected five QTLs for PH, three QTLs for LBL, four QTLs for total branch number (TBN), five QTLs for fruiting branch number (FBN) and four QTLs for INN, most QTLs for these traits were located on A01, A06, A10, B01, B06, B07, and B10. It is showed that several QTLs for different traits were co-localized, such as one region (41.1–46.0 cM) on A06 for traits PH (qPHA06.1 from “WuH16”) and LBL (qLBLA06 from “WuH16”), one region on A06 (10.2–14.4 cM) for traits FBN (qFBNA06.1 from “WuH15”) and TBN (qTBNA06 from “WuH15”), and one region on B01 (67.5–74.0 cM) for traits FBN (qFBNB01 from “WuH15”) and TBN (qTBNB01 from “WuH15”), respectively.For the pod- and seed-related traits, we have identified three QTLs for PL, six QTLs for PW, four QTLs for length-width ratio of pod (LWRP), six QTLs for HSW, two QTLs for HPW and thirteen QTLs for SL, their identified QTLs were localized to A04, B02, B06, B07, and B08, and co-localized to two genomic regions: B06 (121.5–141.0 cM) and B07 (131.2–142.9 cM) (Figure 5), including qPLB06.1–2, qPWB06.1–3, qSLB06.1–3, qHPWB06, qHSWB06.1-2 from “WuH15”, qPWB07.1–3 from “WuH16” and “YangL16,” qHSWB07.1–2 from “WuH15” and “WuH16”, qSLB07.1–4 and qSWB07.1–4 from “WuH15,” “WuH16” and “YangL16.” These results showed that the pleiotropic QTLs contained multiple tightly linked genes for different traits or the gene that affects multiple traits (Hall et al., 2006).

All these annotated genes were identified underlying predicted QTLs and then blasted to find homologous genes with functional descriptions in A. thaliana (Supplementary Table S5). We conducted GO analysis of these genes, and found that many of these distributed in GO categories, such as “GO:0032502 developmental process,” “GO:0010926 anatomical structure formation,” “GO:0040007 growth,” “GO:0065007 biological regulation,” “GO:0030528 transcription regulator activity” (Supplementary Figure S7). We further checked the functional annotation and description of genes within two obviously co-located QTLs of seed- and pod-related traits in B06 and B07 since these two QTL regions were also confirmed by our following Bulk segregant analysis (BSA) of the trait of interest (Unpublished data). There are 63 and 76 candidate genes in the two regions, respectively, and some candidate genes were found to be associated with yield-related traits based on previous studies. In the pleiotropic QTL on B06, two genes, Araip.10014506 and Araip.10014509, were homologous to AT1G08840 and AT1G80260, respectively, which were reported to be related to embryo development (McElver et al., 2001; Jia et al., 2016). In the pleiotropic QTL on B07, the Araip.10036332 and Araip.10036333, a pair of tandemly duplicated genes, were homologous to aap8 (AMINO ACID PERMEASE 8, AT1G10010). This gene in A. thaliana plays an important role in the regulation of carbon metabolism and transport by sink strength (Santiago and Tegeder, 2016, 2017). Carbon fixation and assimilation as well as sucrose partitioning to siliques were strongly decreased when app8 mutant plants transitioned to the reproductive phase. Based on these results, these candidate genes within QTLs should be further investigated, including gene cloning and functional analysis.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.