For DNA extraction, field-grown plant leaves were collected and stored at -80°C. The genomic DNA was extracted using a modified CTAB method (Porebski et al., 1997). For the purpose of analysis, the purified DNA was dissolved in TE buffer. The Nanodrop 2000 was used to measure the DNA’s quantity and quality.

At GeneSeek (Lincoln, Nebraska, USA), the SNP array (Affymetrix 2) was used for the genotyping. No samples were excluded for low call rate (< 0.95) or low quality. After removing SNPs with genotyping error, a call rate of less than 0.95, or a minor allele frequency of less than 0.05, a total of 12,526 SNP markers were retained. STRUCTURE 2.2.3 was used to analyze the population structure of the 102 accessions (Pritchard et al., 2000). In a previous study, the markers’ density and heterozygosity were reported (Zhang H. et al., 2019).

The filtered SNP data was used to calculate the linkage disequilibrium (LD) with the PopLDDecay program (Zhang C. et al., 2019). Standardized disequilibrium coefficients (D’) and squared correlations (r²) were figured. Based on the means of distances between SNP pairs and r² values, LD decay plots were created.

To detect the association between SNPs and flavor-related traits in peanuts, statistical analyses were performed using three different methods: (1) TASSEL 5.0: the PCA (first five PCs) + GLM module (general linear model, GLM) was utilized in the TASSEL software (Bradbury et al., 2007). (2) GAPIT 3.0: association mapping was implemented in R utilizing the MLM module (mixed linear model), Y = Xβ + Zu + e (Lipka et al., 2012). (3) BLUP analysis with mrMLM software (Zhang YW et al., 2020).

The TASSEL module matched the data better than the MLM module, as shown by the Q-Q plots. The TASSEL was used to analyze the phenotype and genotype data. LD pruning was accomplished by PLINK (Purcell et al., 2007) using the indep-pairwise function with a step size of 50, a variant count of 5, and r² threshold of 0.5, which generates 2,363 independent SNPs and LD blocks for this population. Based on the number of independent SNPs and LD blocks, the P-value threshold for significant and suggestive QTLs was calculated (Zhou et al., 2017). The LD block was defined as a collection of adjacent SNPs with a pairwise r² value of at least 0.50 (Gu et al., 2011). A total of 2,363 independent SNPs and LD blocks were retained after LD pruning. Thus, the threshold for significant QTLs was 0.05/2363 = 2.12e−5 with − log10 ^{(P value)} = 4.67. Major or significant QTLs were defined as QTLs that met the significant threshold. In a genome-wide test, “suggestive association” permits one false-positive effect, indicating that the threshold p-value was 1/2363 = 4.23e−4 with − log10 ^{(P value)} = 3.37. Suggestive QTLs were defined as QTLs that were lower than the significant threshold but up to the suggestive threshold. Using qqman, Manhattan plots of the p-value results were created (Turner, 2014).

Genes within the ±80 kb region of the significant and suggestive SNPs associated with peanut flavor were detected from the Tifrunner genome sequences (https://peanutbase.org/). The Tifrunner genome annotation file was used to annotate the identified genes. RNA-Seq data under BioProject PRJNA573070 were downloaded from the NCBI Sequence Read Archive (SRA) database and used for differential expressed genes (DEGs) analysis. Four seed developmental stages (yellow, orange, brown, and black) of four genotypes (AABB, aaBB, AAbb, and aabb) of F8 breeding lines were used for RNA-Seq. For details, see the paper by Zhang et al. (2021).

A total of 102 U.S. peanut mini-core accessions planted in 2013 and 2014 were used in this research. Thirty-three traits associated with flavor, which includes moisture, oil content, tocopherol (total tocopherols, alpha-tocopherol, beta-tocopherol, gamma-tocopherol, and delta-tocopherol), sugars (inositol, glucose, fructose, sucrose, raffinose, stachyose, and total sugars), and descriptive quality (roast peanut, sweet aromatic, raw beany, dark roast, woody hulls skins, earthy musty, cardboardy, plastic chemical, painty, fruity fermented, metallic, spice, sour, sweet, salty, bitter, tongue-throat burn, astringent, and ashy) were measured. Significant variability was identified for these traits among the accessions (Table 1). The total sugars varied significantly (ranging from 22502.56 to 40462.04 mcg/g on average 29297.07 mcg/g). Among the total sugars, variation in sucrose ranged from 19999.40 to 37097.90 mcg/g with an average of 25280.23 mcg/g. In addition, the amounts of total tocopherols (293.49-493.25 mcg/g) and gamma-tocopherol (99.38-276.68 mcg/g) also showed significant variation. There is no significant difference in the traits of descriptive quality.

To investigate the relationships and the main factors underlying trait variation, Pearson correlation and PCA were conducted. For the results of the Pearson correlation, salty has no significant relationship with other traits (Dataset 1). PCA results revealed that each PC can only explain a small variance, such as PC1 explaining 18.11% and PC2 explaining 15.32% of total variance (Table 2A; Dataset 2). This indicated that all of the PCs contributed to total variance without significant difference, so we performed condensed PCA selecting one or two traits from each of the 10 PCs according to the first PCA results. Condensed PCA included 13 traits resulting in five PCs and they explained 68.89% of total variance (Table 2B; Dataset 3). Considering the normal distribution and component percentage of the phenotype data, three traits (total sugars, sucrose, and total tocopherols) were selected for the following GWAS analysis.

After removing SNPs with genotyping error, a call rate of less than 0.95, or a minor allele frequency of less than 0.05, 12,526 SNP markers were retained. No samples were excluded because of low call rates (< 0.95) or poor quality. Population structure results indicated that the 102 mini-core accessions were grouped into primarily two subpopulations according to the delta K (Figure 1A), which was in line with our previous study (Zhang H. et al., 2019). Under K = 2, Group 1 contains 52 mini-core genotypes, most of which belong to hypogaea, while Group 2 mainly includes fastigiata (Figure 1B). Linkage disequilibrium (LD) was estimated from r² (r²< 0.2 was considered unlinked) against distances between each marker in the 102 U.S. mini-core collections. In this population, the LD decay was approximately 0.16 Mb and the r² was 0.2 (Figure 1C).

The selected traits (total sugars, sucrose, and total tocopherols) displayed near-normal distribution for both years (Figure 2). In general, the distributions of each trait varied little from year to year. For total sugars and sucrose, the ratings ranged from 20000 to 45000 and 15000 to 40000, respectively, in both years 2013 and 2014. For total tocopherols, the ratings ranged from 250 to 550 in 2013 and 250 to 500 in 2014. The median of rating scores for every trait in various years is displayed in Figure S1. The estimated broad-sense heritability of total sugars, sucrose, and total tocopherols was 87.47%, 88.76%, and 87.15%, respectively (Table 1).

To choose the most suitable model for GWAS analysis of the selected traits in the population, the PCA + GLM and MLM models were tested. Quantile-quantile (Q-Q) plots suggested that the PCA + GLM model provided a better fit with the expected distribution (Figures 3, S2, S3). Therefore, the PCA + GLM model was chosen for subsequent analysis. In total, 67, 91, and 20 QTLs were identified for total sugars, sucrose, and total tocopherols, respectively (Table 3). Among them, there were seven significant QTLs for total sugars and 22 for sucrose; however, significant QTLs were not identified for total tocopherols. A total of 178 QTLs were found to be distributed in over 15 chromosomes; 81 QTLs were found to be distributed in over six chromosomes of the A sub-genome, and 97 QTLs were found to be mapped in over nine chromosomes of the B sub-genome (Table S2). Of the 29 significant QTLs, 12 were situated on the A sub-genome and 17 were on the B sub-genome (Table S3). ChrB04 had a maximum of 28 QTLs, followed by ChrA04 with 26 QTLs, and ChrA01, ChrA07, ChrA08, ChrA10, and ChrB07 had no QTLs (Table S2). ChrB04 and ChrB05 had the same number of significant QTLs, which is seven (Table S3). In addition, two stable significant QTLs (AX-176797514 and AX-176812379) associated with total sugars were identified in both years and both of them explained 16.13% of the phenotypic variance (Table 4). Similarly, three stable significant QTLs (AX-147248761, AX-176796778, and AX-147248574) located on ChrB04, and two stable significant QTLs (AX-176797514 and AX-176812379) located on ChrB05 were detected to be related to sucrose. On ChrB05, two significant QTL regions (AX-176797514 and AX-176812379) were found to be related to sucrose and total sugars, respectively. More importantly, QTL AX-147221247 on ChrA04 was significantly associated with sucrose and total sugars. For this QTL, our analysis indicated that the C/C genotype had higher sucrose and total sugars contents, while the T/C or T/T genotype had lower sucrose and total sugars contents (Figure 4).

Figure 3 depicts the Manhattan plots of the GWAS results. The genes near significant and suggestive QTLs within ± 80 kb regions are shown in Table 5, Dataset 4, and Dataset 5. A total of four potential candidate genes associated with total sugars near significant QTLs were detected, including one stable gene arahy.BXRJ0F, which spans a genomic region between 19,236,308 and 19,239,067 bp on ChrB05 encoding nucleotide-diphospho-sugar transferase superfamily protein (Table 5). The functional annotation of the other three genes in major QTL regions incorporates two genes (arahy.V3ZQ6H and arahy.94J2Q8) that encode UDP-glycosyltransferase and one gene (arahy.0B2YB3) that encodes acyl-CoA thioesterase. For sucrose, within the significant associated region, eight genes were found, including nucleotide-diphospho-sugar transferase superfamily protein (arahy.BXRJ0F), Gag-pol polyprotein n (arahy.0KX8XR), probable polygalacturonase (arahy.ETW0X5), acyl-CoA thioesterase (arahy.0B2YB3), UDP-glycosyltransferase (arahy.V3ZQ6H), UDP-Glycosyltransferase superfamily protein (arahy.94J2Q8), ATP-binding/protein serine/threonine kinase (arahy.8E8ZZ7), and lysosomal alpha-mannosidase-like (arahy.X2JSTS). The stable gene arahy.BXRJ0F related to sucrose was also detected in both years. In addition, four genes (arahy.BXRJ0F, arahy.V3ZQ6H, arahy.94J2Q8, and arahy.0B2YB3) near significant QTL regions were associated with both total sugars and sucrose, which includes one stable gene (arahy.BXRJ0F).

PCA is a useful method to extract main factors from complex traits that are highly correlated, which covers eigenvectors, standard deviation, and covariance (Karamizadeh et al., 2013). Correlation coefficients were calculated from all 33 investigated traits (Dataset 1). Salty was proved to have no significant correlation with other traits. Total sugars and total tocopherols are not significantly correlated with each other; however, both of them have a significantly positive relationship with moisture and a negative relationship with oil content. Moisture and oil contents have a significantly negative relationship with each other. After Pearson correlation calculation, 32 correlated traits were selected to conduct PCA. Each of the traits contributed to a small variance and the first three PCs can only explain 43.91% of phenotype variance during the first PCA. Then, the condensed PCA was performed to reduce the PCs, and finally, three traits (total sugars, sucrose, and total tocopherols) were selected which can explain the most variance of peanut flavor. In addition, for descriptive analysis, astringent, sour, metallic, and roast_peanutty also need to be considered when detecting peanut flavor. From the Pearson correlation and PCA study, we know which traits have bigger effects on flavor.

GWAS is an effective method to identify QTLs that are linked to the interest traits. However, many variables can affect the consequences of association mapping, which includes genetic diversity, linkage disequilibrium (LD), population structure and size, density of markers, and errors in phenotyping and genotyping data (Gordon and Finch, 2005). At least two methodologies can be used to validate the reliability of GWAS results including approving the QTLs related to the characteristic in various populations and mutual validation through combining association mapping with genetic mapping in Recombinant Inbred Line (RIL) or F2 groups. Abundant genetic diversity, rapid LD decay within a diverse population, and higher density of markers could increase the mapping resolution. In our study, a total of 102 U.S. peanut mini-core accessions were used, and it covered all peanut botanical varieties (var. fastigiata, var. vulgaris, var. hypogaea, var.aequatoriana, var. peruviana, and var.hirsuta) (Table S1). Population structure analysis revealed that although the principal component analysis plot did not show any distinct clusters, the mini-core collections were primarily divided into two subpopulations (Figure S4). Population structure, allele frequency, recombination rate, and selection could also influence the precision of GWAS by affecting the LD decay (Liu et al., 2020). In this panel, the LD decay was around 160 kb for the entire genome, which is consistent with a former report (Zhang et al., 2020). GWAS has been extensively utilized in peanut research including different traits, such as ionomic (Zhang H. et al., 2019), leaf spots (Zhang et al., 2020), domestication (Zhang et al., 2017), and yield-related traits (Wang et al., 2019); however, no reports concerning GWAS of peanut flavor-related traits have been reported. The findings from this study will have the potential to be used in the genetic improvement of peanuts, as well as provide important information for flavor research.

Peanut flavor is a complex agronomic trait influenced by many factors. Peanut seed ingredients including free amino acids, reducing sugars, and amino acids released during protein denaturation are important precursors for flavor development through the Maillard reaction pathway. The scoring test and the hedonic scale method are the two most widely used sensory evaluation methods (Oupadissakoon and Young, 1984). An evaluation sheet for peanuts’ off-flavor and texture was developed by Holaday et al. (Holaday et al., 1964), and Thomas et al. (Thomas et al., 1968) utilized the rank preference test to evaluate peanuts’ quality. Arthur D. Little Co. (Cairncross and Sjöström, 2004) first developed a method for descriptive flavor profile analysis and applied it to foods. Professionally trained panelists were hired to evaluate the peanut flavor by eating peanuts; however, this approach is subjective, time and energy-wasting, and not very accurate. In addition, it is impractical to screen for a variety of flavor chemicals and to select for good flavor by tasting a large number of peanuts. With the rapid advances in molecular breeding tools and QTLs/genes discovery, marker-assisted selection could be used to select peanut with desired flavor.

By using PCA and GWAS, three important flavor chemicals were detected to have major positive effects on peanut flavor, and sugars are the major contributors to peanut flavor. The related QTLs and genes were obtained, and the function of the potential candidate genes detected from this research is related to the covalent addition and transfer of sugars. Within the genomic region of 160 kb, gene arahy.V3ZQ6H (120289043-120292411 bp) and arahy.94J2Q8 (120315438-120323907 bp), coding for UDP-Glycosyltransferase superfamily, were found to be located on chromosome A04, which catalyzes the covalent addition of sugars from nucleotide UDP-sugar donors to functional groups (amine, hydroxyl, or carboxyl) on an expansive scope of lipophilic molecules (Meech et al., 2019). Glycosyltransferases include multiple gene families that play a role in a process called glycosylation in plant secondary metabolites, which is necessary for the maintenance of cellular homeostasis (Yu et al., 2017). The potential polygalacturonase is encoded by gene arahy.ETW0X5 plays a role in the biosynthesis of carbohydrates. Some other candidate genes including arahy.BXRJ0F (nucleotide-diphospho-sugar transferase superfamily), arahy.0KX8XR (Gag-pol polyprotein), arahy.0B2YB3 (acyl-CoA thioesterase), arahy.8E8ZZ7 (ATP-binding/protein serine/threonine kinase), and arahy.X2JSTS (lysosomal alpha-mannosidase) were found participating in carbohydrate metabolism near significant QTLs regions. In addition, gene expression analysis indicated that most of these genes were involved in the carbohydrate metabolic process and were mainly expressed at the early seed development stages, such as genes arahy.X2JSTS, arahy.BXRJ0F, arahy.ETW0X5, and arahy.94J2Q8. There was no significant difference for gene arahy.0B2YB3 at different seed development stages in four peanut genotypes (Figure 5).

Haplotype analysis was conducted for the identified genes; however, no SNPs were detected in the gene regions due to the limitation of the SNP numbers. In this research, the SNPs were from the SNP array, not genome resequencing, which limited the analysis. In addition, the number of accessions also had an influence on the results. To address these problems, whole genome resequencing and multi-environment joint GWAS via 3VmrMLM should be conducted in the future, and the candidate genes should be further confirmed by haplotype analysis and multi-omics analysis.