Multi-Locus Genome-Wide Association Studies Reveal Fruit Quality Hotspots in Peach Genome

Peach is one of the most important fruit crops in the world, with the global annual production about 24.6 million tons. The United States is the fourth-largest producer after China, Spain, and Italy. Peach consumption has decreased over the last decade, most likely due to inconsistent quality of the fruit on the market. Thus, marker-assisted selection for fruit quality traits is highly desired in fresh market peach breeding programs and one of the major goals of the RosBREED project. The ability to use DNA information to select for desirable traits would enable peach breeders to efficiently plan crosses and select seedlings with desired quality traits early in the selection process before fruiting. Therefore, we assembled a multi-locus genome wide association study (GWAS) of 620 individuals from three public fresh market peach breeding programs (Arkansas, Texas, and South Carolina). The material was genotyped using 9K SNP array and the traits were phenotyped for three phenological (bloom date, ripening date, and days after bloom) and 11 fruit quality-related traits (blush, fruit diameter, fruit weight, adherence, fruit firmness, redness around pit, fruit texture, pit weight, soluble solid concentration, titratable acidity, and pH) over three seasons (2010, 2011, and 2012). Multi-locus association analyses, carried out using mrMLM 4.0 and FarmCPU R packages, revealed a total of 967 and 180 quantitative trait nucleotides (QTNs), respectively. Among the 88 consistently reliable QTNs detected using multiple multi-locus GWAS methods and/or at least two seasons, 44 were detected for the first time. Fruit quality hotspots were identified on chromosomes 1, 3, 4, 5, 6, and 8. Out of 566 candidate genes detected in the genomic regions harboring the QTN clusters, 435 were functionally annotated. Gene enrichment analyses revealed 68 different gene ontology (GO) terms associated with fruit quality traits. Data reported here advance our understanding of genetic mechanisms underlying important fruit quality traits and further support the development of DNA tools for breeding.


INTRODUCTION
Peach [Prunus persica (L.) Batsch] is a diploid species, with a short juvenile period (2-4 years), relatively simple genome (∼230 Mbp), and one of the best genetically characterized deciduous trees . Peach is the third most cultivated temperate tree fruit in the world, after apple and pear, with a world production of approximately 24.6 million tons (Food and Agricultural Organization of the United Nations (FAOSTAT), 2018). Despite the high production, peach consumption has declined over the past decades. In the United States, peach per capita consumption decreased to 1.3 kg per year compared to ∼3 kg per year in the 1980s (Minas et al., 2018). Inconsistent and low fruit quality is recognized as the major limiting factor for consumer acceptance and, consequently, the low rates of peach consumption (Cirilli et al., 2016).
Peach breeders have always selected for fruit quality with respect to size, color and firmness, as well as tried to expand harvest season (Laurens et al., 2018). Recently, more emphasis is on other traits such as internal quality and postharvest traits (Elsadr et al., 2019).
Recent advances in next-generation high-throughput sequencing and genotyping techniques, such as development of the 9K peach SNP array by the International Peach SNP Consortium (IPSC) (Verde et al., 2012), allow use of DNA information to develop tools for facilitating breeding efforts da Silva Linge et al., 2018). Understanding the genetic mechanisms that control a specific trait would enable peach breeders to efficiently apply marker-assisted breeding (MAB) through the development of DNA diagnostic tools, and consequently select seedlings with desired quality traits early in the selection process before the characters can be evaluated in the field (Abdelghafar et al., 2020).
The link between the genetic markers and a particular trait could be determined using different approaches. Quantitative trait loci analysis (QTL mapping) and genomewide association studies (GWAS) are widely used for dissection of complex genetic traits (Meneses and Orellana, 2013). In peach, several linkage maps have been used in QTL discovery of key fruit quality traits such as fruit size, diameter, firmness, acidity, soluble solid concentration, individual sugars, maturity date, pubescence, blush, fruit texture, and phytochemical compounds (Eduardo et al., 2011;Martínez-García et al., 2013;Pirona et al., 2013;Frett et al., 2014;Vendramin et al., 2014;da Silva Linge et al., 2015;Zeballos et al., 2016;Ciacciulli et al., 2018;Nuñez-Lillo et al., 2019;Abdelghafar et al., 2020). These maps were typically developed for mapping particular traits in a specific parental background with limited recombination events and genetic diversity.
Alternatively, GWAS has the advantage of increasing the recombination events and consequently mapping resolution with a significant reduction of the research time (Zhu et al., 2008). However, false positives due to population structure or kinship among genotypes, or false negatives due to removal of rare alleles that are involved in natural variation are some of the weaknesses of GWAS (Brachi et al., 2011). To deal with this problem, GWAS methods utilizing mixed linear models (MLM), which take into account multiple levels of relatedness, have become standard methodology (Yu et al., 2006). Significant marker-trait association based on the single-locus models, such as the general linear model (GLM) and MLM, were reported for several traits such as fruit pubescence, fruit shape, stone adhesion-flesh texture, fruit flesh color, non-melting/melting flesh, fruit weight, titratable acidity, soluble solid concentration, leaf gland type, flower type, bloom date, fruit development period, maturity date, ripening index, and total sugars Cao et al., 2016Cao et al., , 2019Elsadr et al., 2019;Font I Forcada et al., 2019). Singlelocus models test one locus at a time and fail to match the true genetic model of complex traits that are controlled by numerous loci simultaneously (Xu et al., 2018). Thus, major improvements in GWAS statistical methodology have occurred, and multilocus GWAS methods considering the information of all loci simultaneously have been developed .
Recently, six multi-locus GWAS approaches were integrated into an R package, named mrMLM . The mrMLM 4.0 R package comprises the mrMLM , FASTmrMLM (Tamba and Zhang, 2018), FASTmrEMMA , ISIS EM-BLASSO , pLARmEB , and pKWmEB (Ren et al., 2018) two-step multi-locus GWAS methods. First, various algorithms are used to select all potentially associated markers. Second, these selected markers are put in one model, in which all the effects are obtained by empirical Bayes, and all the non-zero effects are further identified by likelihood ratio test for true Quantitative Trait Nucleotides (QTNs) .
The multi-locus model Fixed and random model Circulating Probability Unification (FarmCPU) uses the associated markers as covariates in a fixed-effect model (FEM) and optimization on the associated covariate markers in a random effect model (REM). FarmCPU adopts REML optimization to replace the criterion that the variance explained by kinship is near zero, which can only be arbitrarily determined. FarmCPU also adopted a binning approach from super to select pseudo QTNs. The whole genome is equally divided into bins, and only one significant marker with the smallest P-value from each bin is selected as the candidate pseudo QTN. These candidate pseudo QTNs are determined by a REM. The candidate pseudo QTNs are first ranked by P-value. Then, the best combinations between the different bins and the number of candidate pseudo QTNs are determined by REM. Finally, the two types of models (FEM and REM) are performed iteratively until no change occurs in the selection of pseudo QTNs (Huang et al., 2018). Thus, FarmCPU decreases the computer time required, provides reliable results by efficiently removing the confounding between the population structure and Kinship, avoiding model over-fitting, and controlling for false positives (Liu et al., 2016).
The objective of this study was to identify significant markertrait association for 14 agronomic traits, using the multi-locus GWAS methods in mrMLM 4.0 and FarmCPU in a U.S. peach diversity germplasm panel of 620 individuals, managed by three public fresh market peach breeding programs at University of Arkansas System Division of Agriculture, Texas A&M University and Clemson University.

Plant Material, DNA Isolation, Quantification, and Genotyping
The material used in this study represents the U.S. peach breeding germplasm assembled under the RosBREED project (Iezzoni et al., 2010(Iezzoni et al., , 2020Peace et al., 2014). A total of 72 cultivars/advanced selections and 548 individuals from three public fresh market peach breeding programs: University of Arkansas System Division of Agriculture (AR), Clemson University (SC), and Texas A&M University (TX), were chosen to effectively represent alleles currently found within North American fresh market peach breeding germplasm (Supplementary Table 1).
Peach DNA was extracted from young leaves using the E-Z 96 Tissue DNA Kit (Omega Bio-Tek, Inc., Norcross, GA, United States). DNA was quantitated with the QuantiT PicoGreen Assay (Invitrogen, Carlsbad, CA, United States), using the Victor multi-plate reader (Perkin Elmer Inc., San Jose, CA, United States). The final DNA concentrations were adjusted to a minimum of 50 ng/µL and submitted to the Research Technology Support Facility at Michigan State University (East Lansing, MI, United States).
Samples were genotyped with the IPSC peach 9K SNP array v1 (Verde et al., 2012). The SNP data curation was performed using the workflow for high-resolution genetic marker data described in Vanderzande et al. (2019). After the SNP data curation, a total of 4005 SNPs distributed over the eight peach chromosomes remained and were used in the multi-locus GWAS (Supplementary Table 2).

Phenotypic Data
Phenotypic data were recorded over three seasons (2010)(2011)(2012) at each fresh market peach program. Bloom data (BD; Julian days) were visually assessed in the field and recorded for each tree when 60-80% of the blossoms were open. Ripening date (RD; Julian days) was determined when 20% of fruits were at commercial harvest by visually inspecting the presence of a few soft fruits in the field for maturity two times per week. Days after bloom (DAB; Julian days) was calculated as the number of days between the date of full bloom and ripening date.
Approximately 20 fruits were harvested for phenotyping. A five firm fruit sample was selected for the following traits evaluations: Blush (0-5 scale, 0 = none, and 3 = 40-60%, 5 > 90% red blush on fruit surface) subjective scales were used as described by Frett et al. (2014). Fruit diameter (FDIA; mm) was evaluated with a millimeter caliper, while fruit weight (FW; g) was measured as the average weight of the five selected peaches. Flesh adherence (ADH) was evaluated using 1 -4 scale where 1 = Freestone; 2 = Semi-freestone; 3 = Semi-clingstone; and 4 = Clingstone. Fruit firmness (FF; N) was measured using an electronic fruit texture analyzer (FTA) fitted with an 8-mm diameter tip (GÜSS Fruit Texture Analyzer; GÜSS Manufacturing (Pty) Ltd., Strand, South Africa). All readings were recorded as kilogram-force (kgf) and then converted to Newton (N) by multiplying the reading by 9.807. Redness around Pit (RP) was measured following the scale 1 = red; 0 = no red. The fruit texture (FT) was evaluated using the scale 1 = melting; 2 = non-melting. Pit weight (PW; g) was measured as the average weight of the five selected pits.
For biochemical traits, a composite sample of one approximately 2 cm wide longitudinal slice from each of the five fruits was used to extract juice with a juicer for the measurement of soluble solid concentration (SSC) using a digital refractometer, pH with a pH meter and titratable acidity (TA) using an automatic titrator (DL 22 Food and Beverage analyzer,Mettler Toledo,Columbus,OH,United States). TA was obtained by the titration of solution of 6 g of the peach juice diluted with 50 mL of distilled water to pH 8.2 with 0.1N NaOH and expressed as milliequivalents of malic acid. The following equation was used to calculate titratable acidity (the milliequivalent factor used corresponded to malic acid, 0.067): Titratable acidity (%) = [NaOH titrated ml × 0.1N (NaOH) × milliequivalent factor × 100] 6 g of juice Descriptive Analysis, Genetic Diversity, and Population Structure The descriptive analysis and the correlations between the traits were performed using the software Past (Hammer et al., 2001). The genetic diversity analysis was performed using the GenAlEx software (Peakall and Smouse, 2012). The narrow sense heritability was calculated using the R package Sommer (Covarrubias-Pazaran, 2016) using the h2.fun: h2.fun (object, data, gTerm, eTerm) where: object represents a model fitted with the mmer function; data represents the dataset used to fit the model provided in the object argument; gTerm is a character vector specifying the genetic terms fitted in the model; and eTerm is a character vector specifying the environment term fitted in the model. For the level from the eTerm (environment) the heritability is calculated as: "PEV" is the predicted error variance for the genotype, "md" is the mean value from the diagonal of the relationship (genomic) matrix "G" and where "Vg" refers to the genotype variance. The model included in the h2.fun was: where "K" refers to the genomic relationship matrix. Population structure, multidimensional scaling (MDS) and Bayesian clustering were performed with fastSTRUCTURE (Raj et al., 2014). The MDS was performed using TASSEL (Bradbury et al., 2007). The MDS results were plotted with the R package "scatterplot3D" (Ligges and Mächler, 2002). The fastSTRUCTURE was run with a "simple prior" option and remaining default parameters. The number of populations (K), ranging from 1 to 20, and the most probable number of populations was chosen for running the built-in script for multiple choices of K. The admixture proportions of each genotype, estimated by fastSTRUCTURE, were visualized using DISTRUCT plots (Rosenberg, 2004). Accessions were assigned to a specific subpopulation when the estimated membership coefficients (Q) were above 0.80.
Linkage disequilibrium (LD) was measured by correlation coefficients (r 2 ) for all pairs of SNPs. The LD decay were calculated using PopLDdecay (Zhang et al., 2018a) with the following parameters: -MaxDist 3000 kb -MAF 0.05.

Genome-Wide Association Study
To validate and increase the accuracy of the multi-locus GWAS results, we used mrMLM 4.0  and FarmCPU (Liu et al., 2016). The six multi-locus GWAS methods (mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS EM-BLASSO) from mrMLM 4.0 R package were used. The SNP data were converted to character, as described in the user manual, the population structure was the Q matrix obtained from fastSTRUCTURE and the kinship matrix was calculated by mrMLM 4.0. All parameters in GWAS were set at default values. The significantly associated SNPs were determined by the critical threshold of LOD score ≥ 3 as described in previous studies . Concerning FarmCPU, the SNP data were converted to numerical using the R package GAPIT (Lipka et al., 2012). Principle component analysis (PCA) was conducted using TASSEL 5.0, and the first three components were incorporated as covariates in the GWAS model. Bonferroni-corrected P-value threshold was set at p < 0.01.
We considered a QTN reliable when: QTNs repeatedly detected in at least four methods and/or two seasons using the mrMLM 4.0; QTN consistently detected in two seasons using FarmCPU; QTNs detected in at least three methods in mrMLM 4.0 and also identified in the FarmCPU approach. These QTNs were named as "qtn" + trait name abbreviation + scaffold + detected QTL order on chromosome.

Candidate Genes
The candidate gene analysis was performed using two strategies. First, the candidate gene analyses were performed within the haploblock regions in which a QTN was detected with at least three methods in mrMLM 4.0 and with FarmCPU. Haploblock regions encompassing the associated SNPs were determined in PLINK 1.9 (Chang et al., 2015) using the flag "blocks" restricted to 500 kb. From the Prunus persica Whole Genome v2.0 Assembly & Annotation v2.1 (Verde et al., 2017) in Genome Database for Rosaceae 1 (Jung et al., 2018), a systematic search was conducted to compile the predicted candidate genes associated with the quality traits. The candidate genes were further analyzed for GO (gene ontology) enrichment using GOseq 1.42.0 R package (Young et al., 2010). The GO terms were considered significantly enriched or depleted at FDR < 0.05. The enriched GO terms were visualized using REVIGO semantic similarities (Supek et al.,

Phenotypic Data
Six hundred twenty individuals from the three fresh market public peach breeding programs were evaluated for 14 different fruit quality traits over 3 years (2010-2012) ( Table 1) Highly significant (P < 0.01) correlations were observed between the traits (Supplementary Table 3). The highest correlation was observed between the FDIA_2011 and FW_2011 (0.92). As expected, a significant negative correlation was detected between the pH and TA (−0.65 and −0.64 in 2011 and 2012, respectively). The traits DAB and RD revealed a significant positive correlation in the years analyzed (0.82 and 0.87). Concerning the correlation between years, BD showed the highest correlation (0.99), followed by RD (0.92).
The narrow sense heritability (h 2 ) was estimated for all 14 traits (Supplementary Table 4
The population structure was analyzed with Multidimensional scaling (MDS) and fastSTRUCTURE. The MDS revealed two main groups, in which the second group could be divided in two clusters. The first group comprised the individuals from the TX breeding program, while the second group grouped the individuals from the AR and SC breeding programs 1 | Descriptive analysis of 14 phenotypic traits observed in 620 individuals from three U.S. public fresh market peach breeding programs (Univ. of Arkansas, Texas A&M and Clemson Univ.) over three seasons (2010)(2011)(2012).
(Supplementary Figure 1). Population structure analysis with fastSTRUCTURE suggested a number of K between 2 and 19. However, the population stratification for K = 3 showed clear differences between groups based mainly on the pedigree information of the individuals belonging to each group (Supplementary Figure 1). The first group accounted for the individuals related to 'Tropic Beauty' , 'TX2293_3' , 'TX2B136' , 'TXW1293_1' , and 'TXW1490_1'. The second group comprised individuals linked to ' A_663' , ' A_760' , and 'Bolinha' , while the third group contained individuals linked either to 'Clayton' and/or 'O'Henry'. The LD decayed with increase of physical distance between SNPs in all groups (Supplementary Figure 2). Considering the admixed individuals, the average of r 2 was 0.16. The physical distance over which LD decayed to half of its maximum value was around 540 kb. Different patterns of LD decays were observed in the three different groups. Group 3 revealed the highest average of r 2 (0.32) and the longest physical distances in which LD decayed to half of its maximum value (1,620 kb), while group 2 showed shortest distance (480 kb). In the group 1, the LD decayed of its maximum value of r 2 in ∼540 kb.

Multi-Locus Genome-Wide Association Study
GWAS using the six multi-locus methods in the R package mrMLM 4.0 revealed a total of 967 QTNs associated with 14 traits (Supplementary Table 5). The highest number of associated SNPs was observed on chromosome 4 (99) and the lowest in chromosome 7 (23). Significant QTNs detected in at least four methods in the same season, were detected for almost all traits except TA. In addition, consistently associated QTNs identified in at least two seasons were detected for BD, RD, DAB, ADH, RP, SSC, Blush, FF, FT, TA, and pH. Furthermore, SNPs associated with more than one trait were identified on chromosome 1 (BD The multi-locus model FarmCPU revealed a total of 180 QTNs (Supplementary Table 6). The highest number of QTNs were detected on chromosome 4 (33), while the smallest was observed on chromosome 7 (6). Consistently associated SNPs over at least two seasons were identified for BD, RD, ADH, RP, and SSC. In addition, SNPs associated with more than one trait were detected on chromosomes 1 (DAB and RD), 3 (Blush and RP; FW and RD), 4 (Blush and RD; RD and SSC; FF and SSC), 5 (TA and pH), and 8 (BD and RD; FT and SSC).
To ensure reliable results, further analyses included only QTNs that met the following conditions: QTNs detected in at least four methods in mrMLM and/or detected in at least two seasons using mrMLM (Table 2 and
The qtnFDIA_7.1 associated with FDIA was located on chromosome 7, revealed LOD scores ranging from 3.0 to 5.1 and explained 4.8-7.9% of the phenotypic variance observed.
The qtnTA_5.1 on chromosome 5 was associated with TA and detected in two seasons. The LOD score varied from 3.1 to 20.4 and explained 1.0-38.0% of the phenotypic variation.

QTNs Detected in Two Seasons Using FarmCPU
The qtnBD_4.2 was consistently associated with BD in the two seasons and revealed a p values of 1.37E-07 and 7.74E-07, respectively (Table 3 and Figure 1).
The qtnRD_4.7 located on chromosome 4 (10.9 Mbp) was significantly associated with RD in all three seasons. In addition, the qtnRD_4.8 and qtnRD_8.1 were associated with RD in 2011 and 2012 and were located at 12.5 Mbp (chromosome 4) and 2.5 Mbp (chromosome 8), respectively.
For RP, the qtnRP_3.1 located at 18.17 Mbp on chromosome 3 was detected in two seasons with p values of 6.64E-11 and 2.93E-16, respectively.
A total of 566 candidate genes (CG) were detected within the haploblock regions for the significantly associated QTNs (Supplementary Table 7), from which 93 CG were detected in the regions for BD, 89 for RD, 29 for DAB, 22 for Blush, 39 for FDIA, 24 for FW, 26 for ADH, 31 for FF, 148 for RP, 12 for FT and 90 for SSC. The gene ontology (GO) annotations were retrieved for 435 CG. The GO enrichment analysis revealed 68 GO terms in all three GO aspects, biological process, molecular function, and cellular component. Twenty-six GO terms (78 genes) were described as biological processes, 32 GO terms (108 genes) with the molecular function, and 10 GO terms (36 genes) with the cellular component (Supplementary Table 8). The GO term cluster representatives were joined into "superclusters" of terms loosely related to cellulose microfibril organization, THO complex part of transcription export complex and sulfotransferase activity in the biological process, cellular component and molecular function, respectively (Figure 2).

Hotspots in Peach Genome
The reliable QTNs revealed fruit quality hotspots in the peach genome (Figure 1). On chromosome 1, three reliable QTNs associated with BD, Blush and FW in the interval of 25.5-27.2 Mbp were identified. In addition, at the bottom of the same chromosome (43.6-46.4 Mbp), we also observed QTNs associated with BD and FF. On chromosome 3, a hotspot involving the quality traits Blush and RP was observed in the region located at 18.2-20.2 Mbp. The majority of the reliable QTNs detected were located on chromosome 4 (0.6-19.0 Mbp), especially concentrated in the genomic region located at 9.0-12.5 Mbp with QTNs associated with DAB, FW, RP, ADH, RD, FF, and SSC. A hotspot was also observed on the top of chromosome 5 (0.3-3.7 Mbp) with significant signals associated with FT, pH, SSC, and TA. In the genomic region on chromosome 6 spanning 28.8-29.5 Mbp, QTNs involving PW, RP, and FW were detected. Furthermore, on top of chromosome 8 (2.5-5.1 Mbp), a hotspot with reliable QTNs associated with RD, RP, and FT was observed.

DISCUSSION
We have analyzed peach germplasm containing 620 individuals from three U.S. public fresh market breeding programs [University of Arkansas System Division of Agriculture (AR), Clemson University (SC) and Texas A&M University (TX)] for 14 traits over three seasons (2010, 2011, and 2012). Phenotypic variation was observed between individuals and seasons, and the mean values for BD, RD, FW, and SSC were lower than those reported in the Spanish and European germplasm (Hernández Mora et al., 2017;Font I Forcada et al., 2019).However, average values for RD and DAB observed in our study were in agreement with the values reported in the University of Guelph's peach germplasm, comprised of accessions originating from different regions across North America (Elsadr et al., 2019). A high and significant correlation between FW and FDIA (0.92) was previously observed in peach Abdelghafar et al., 2020), as well as the positive correlation between RD and DAB (Elsadr et al., 2019) and the negative correlation between TA and pH (Abidi et al., 2011). In addition, the high estimated narrow sense heritability coefficients observed in this study ranging from 0.68 to 0.87, suggesting that the phenotypic variations of all traits are mainly affected by genetic factors, and therefore this dataset can be used for further genetic analyses. The mean observed heterozygosity (Ho = 0.36) in the U.S. peach germplasm was greater than that observed in the germplasm from four European, one Chinese and one Brazilian peach collections reported in previous studies Thurow et al., 2019). In addition, the mean inbreeding coefficient of 0.05 indicated a low level of inbreeding. The low mean of the inbreeding coefficient observed in this study could be attributed to the diverse material, including F1 and F2 populations with different genetic backgrounds.
The multidimensional scaling (MDS) clustered material into two groups, a group of individuals from TX fresh market breeding program related to 'Tropic Beauty' , 'TX2293_3' , 'TX2B136' , 'TXW1293_1' and 'TXW1490_1' , and the second group, comprised of individuals from the AR and SC fresh market breeding programs. Breeding material from AR and SC clustered in one main group, due to the common founders or pedigreelinkages of some breeding populations. Nevertheless, the second group could further be separated in two clusters in which the first cluster grouped the individuals linked to ' A_663' , ' A_760' , and 'Bolinha' , and the second cluster contained individuals linked to 'Clayton' and/or 'O'Henry'. The population structure indicated by fastSTRUCTURE, between K = 2 and 19, supported MDS clustering, as the K = 3 reflected grouping based on the pedigree background, and number of the breeding programs in the panel.
The population structure influences LD patterns within the genome (Thurow et al., 2019). The LD detected in this study decayed much slower in comparison with the observed by Thurow et al. (2019) and faster than the observed by Micheletti et al. (2015). The difference could be explained by the genetic material analyzed and the methods used for analyzing the LD decay. In addition, a slower decay of LD is expected in selfing materials. The LD decay over distance also determines the number of markers required to cover the genome. Considering the LD decay of our dataset (540 Kb), approximately 421 SNPs covering the total peach genome (227.4 Mb) should be sufficient to perform the GWAS. However, domestication regions containing key genes require more SNPs due to the faster LD decay (Cao et al., 2016).

Multi-Locus GWAS
In order to control the false positive rate in GWAS analysis, conservative correction methods such as false discovery rate (FDR) and Bonferroni correction are frequently adopted in association studies. However, these corrections are often too conservative for detecting many important loci. Thus, multilocus GWAS methods have been recommended to overcome the problem of stringent correction . In this study, we have successfully performed a genome-wide association study using six multi-locus GWAS methods (mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS EM-BLASSO) comprised in mrMLM 4.0 and FarmCPU R packages.The mrMLM 4.0 adopts the critical probability value or log of odds (LOD), a less stringent significance threshold while FarmCPU requires Bonferroni correction to detect QTNs. The multi-locus methods detected 967 and 180 QTNs using mrMLM 4.0 and FarmCPU, respectively, allowing the identification of important regions in the peach genome that control fruit quality traits. Furthermore, consistently reliable QTNs (88) for all traits were detected using different multi-locus GWAS methods and/or at least two seasons (Tables 2-4). Half of the reliable QTNs (44) detected have already been reported using different progenies, germplasm and approaches. However, to our knowledge, the other 44 reliable QTNs controlling fruit quality traits in peach have not been previously described.
One of the main goals of breeding programs is the development of commercial varieties with predictable bloom time to adapt to various target environments. Therefore, understanding the genetic architecture of phenology-related traits represents a key prerequisite to enable the development of varieties adapted to different climates (Gogorcena et al., 2020). Reliable QTNs associated with bloom date (BD), the qtnBD_1.1, qtnBD_1.2, qtnBD_1.3, qtnBD_1.4, qtnBD_4.1, qtnBD_4.3, qtnBD_7.1 collocate near or in the same regions previously reported using QTL mapping and pedigree-based analysis (PBA) (Fan et al., 2010;Romeu et al., 2014;Bielenberg et al., 2015;Hernández Mora et al., 2017;Rawandoozi et al., 2020a). The fact that these regions were identified following different approaches (linkage analysis, PBA analysis and GWAS) in diverse genetic material, makes them an interesting source of allelic variation for BD in peach.
QTL clusters for RD and DAB were commonly detected in peach (Fresnedo-Ramírez et al., 2015;Hernández Mora et al., 2017;Elsadr et al., 2019;Rawandoozi et al., 2020a). In this study, two reliable QTNs (qtnDAB_4.1 and qtnDAB_4.2) associated with DAB overlapped with the QTNs associated with RD. The position of the qtnDAB_4.1 (chr 4: 10.7 Mbp) matched the associated SNPs identified in a panel of 132 peach accessions genotypically characterized via genotyping by sequencing (GBS) approach (Elsadr et al., 2019). In addition, the SNP_IGA_410398 was emphasized as a predictive SNP for RD and DAB in haplotype analysis in a DAB QTL detected using PBA approach (Rawandoozi et al., 2020a). The qtnDAB_4.2 (10.9 Mbp) was close to QTL for DAB detected using PBA in the European germplasm (Hernández Mora et al., 2017). Thus, our results confirmed the location of RD and DAB associated regions in the peach genome, and due to their importance for breeding, could be useful in selection of various phenology patterns in future studies.
Previous QTL analyses in peach have identified a major QTL for blush on chromosome 3 accounting, in average, for 63.7% of observed phenotypic variation (Frett et al., 2014). The qtnBlush_3.1, qtnBlush_3.2 and qtnBlush_3.3 were located within the genetic interval of the major QTL for blush on chromosome 3. QTL regions for blush on chromosome 4 were mapped approximately at 10.5 -11.5 Mbp (Rawandoozi et al., 2020b);11.2 -14.1 Mbp (Hernández Mora et al., 2017); and at 19.8 and 28.6 Mb (Shi et al., 2020). The qtnBlush_4.1 was detected at 6.6 Mbp and accounted for the highest phenotypic variation observed. Moreover, we identified two QTNs on chromosome 1 (qtnBlush_1.1: 26.9 Mbp; qtnBlush_1.2: 2.5 Mbp). Analyzing an F1 peach population derived from the cross between "Shahong" and "Hongfurong, " Shi et al. (2020) also observed a QTL associated with blush on chromosome 1. However, the genetic interval was approximately at 21.5 Mbp. Although the percentage of the phenotypic variation explained was low, a QTN associated with blush (qtnBlush_5.1) on chromosome 5 was also detected.
Understanding the genetic control of fruit diameter and weight is an important goal of breeding programs due to the importance of these traits for the fresh market (Yue et al., 2014). QTL regions associated with fruit diameter and weight have been detected in all chromosomes Fresnedo-Ramírez et al., 2016;Zeballos et al., 2016;Hernández Mora et al., 2017;Cao et al., 2019;Abdelghafar et al., 2020;Shi et al., 2020). In this study, we identified two reliable QTNs associated with FDIA (qtnFDIA_7.1 and qtnFDIA_1.1). The qtnFDIA_7.1, on chromosome 7, is in the vicinity to fruit width and fruit depth QTLs (qP-Fwd7.2 and qP-Fd7.2) reported by da Silva  using F 2 progeny resulting from a cross between an ornamental peach PI91459 ("NJ Weeping") × "Bounty." The qtnFDIA_1.1 (36.9 Mbp) was identified in a different region of chromosome 1 when compared with previous linkage analyses, where QTLs were located approximately at 11 Mbp ; 27-28 Mbp (Hernández Mora et al., 2017); 41 Mbp (Zeballos et al., 2016) and43 Mbp (da Silva Linge et al., 2015;Abdelghafar et al., 2020). Concerning position of qtnFW_1.2 and qtnFW_6.2 matched the QTL interval identified in European peach germplasm (Hernández Mora et al., 2017). The FW QTN qtnFW_2.1 and qtnFW_3.1 could be the FW QTLs mapped in an interspecific cross between peach and a wild relative Prunus davidiana (Quilot et al., 2004;Desnoues et al., 2016). On chromosome 4 qtnFW_4.1 was close to the FW QTL reported by da Silva Linge et al. (2015) and qtnFW_4.2 was in the same genetic interval of the QTL identified by Shi et al. (2020).
Fruit firmness (FF) represents an essential indicator of fruit quality for peach consumers. For this reason, several authors have investigated the genetic mechanisms controlling this trait in peach (Peace et al., 2005;Eduardo et al., 2011Eduardo et al., , 2015Martínez-García et al., 2013;Nuñez-Lillo et al., 2015;Zeballos et al., 2016;Serra et al., 2017;Carrasco-Valenzuela et al., 2019). The qtnFF_4.1 and qtnFF_4.2, reported in this study, were located in chromosome 4 with the position matching the QTL interval associated with firmness loss mapped by Serra et al. (2017). Moreover, the qtnFF_4.3 was in the same genetic region in which Carrasco-Valenzuela et al. (2019)  Flesh adherence to the pit (ADH) is another factor determining overall peach fruit quality, with consumers preferring freestone or semi freestone characteristics (Olmstead et al., 2015). Previous studies have shown that ADH is inherited and controlled by the Freestone-Melting (F-M) locus on chromosome 4, with genes encoding endopolygalacturonase (endoPG) associated with this trait (Peace et al., 2005;Gu et al., 2016). We detected qtnADH_4.1 and qtnADH_4.2 explaining the majority of the phenotypic variation close to the genetic region where the endoPG gene is located. Similar to Shi et al. (2020), we detected significant genetic regions associated to ADH in different regions of chromosome 4 and also in other chromosomes.
The QTNs associated with redness around the pit (RP) (qtnRP_3.1 and qtnRP_3.2) located approximately at 18.2 and 18.7 Mb, on chromosome 3, matched the position of the associated signals to flesh color around the stone detected in the recent GWAS using genome structural variations (SVs) (Guo et al., 2020). In addition, the Cs locus associated with red color around the pit was previously mapped in the middle of chromosome 3 (Yamamoto et al., 2001). Interestingly, the SNP_IGA_341962 (qtnRP_3.1) was also associated with blush (qtnBlush_3.2). Therefore, the QTNs associated with RP identified on chromosome 4 (qtnRP_4.1, qtnRP_4.2 and qtnRP_4.3) were close to the associated signals detected by Guo et al. (2020) and in a different region of the associated SNPs reported by Cao et al. (2016), while the QTNs detected on chromosome 6 (qtnRP_6.1) and 8 (qtnRP_8.1) were located in different regions when compared with previous studies (Cao et al., 2016;Guo et al., 2020).
Concerning pit weight (PW), the qtnPW_6.1, identified on chromosome 6, is close to the QTL (qSW6; 24.6 Mb) mapped in the interspecific cross between almond × peach population (Donoso et al., 2016). Cao et al. (2016) also detected a region significantly associated with PW on chromosome 6; however, the location was approximately at 26.9 Mbp.
Soluble solid concentration is one of the most important quality traits in peach, with consumers expecting an enhanced sugar content or sweetness perception for the low acid types (Cirilli et al., 2016). Therefore, SSC has been a target trait in several studies involving intra-and interspecific progenies and germplasm to access the genetic potential and consequently improve the sugar content in new cultivars. We detected QTNs associated with SSC on chromosome 1, 4, 5, and 6. The qtnSSC_5.1 on chromosome 5 that explained the majority of the phenotypic variation was in agreement with the QTL interval reported by Hernández Mora et al. (2017) using a PBA analysis in European peach germplasm. In the same chromosome, we also identified qtnSSC_5.2 (0.7 Mb) and qtnSSC_5.3 (13.0 Mb) whose positions matched QTLs mapped in previous studies using different germplasm and approaches (Nuñez-Lillo et al., 2019;Abdelghafar et al., 2020;Rawandoozi et al., 2020b). Furthermore, qtnSSC_4.1 on chromosome 4 (6.6 Mb) was close to the QTL (qSSC.V-Ch4-2007a) detected by Zeballos et al. (2016), while qtnSSC_4.3 was near MD locus reported by Eduardo et al. (2011). On the other hand, the qtnSSC_1.1 (chromosome 1: 17.5 Mb) and qtnSSC_6.1 (chromosome 6: 21.2 Mbp) were located in a different region in comparison to the QTLs or associated markers previously detected on those chromosomes (Fresnedo-Ramírez et al., 2015;Cao et al., 2016;Hernández Mora et al., 2017;Li et al., 2019;Shi et al., 2020). Concerning the traits pH and TA, the qtnTA_5.1 and qtnpH_5.1 collocated with the major locus for low-acid fruit (D-locus) previously reported in peach (Boudehri et al., 2009).

Fruit Quality Hotspots in Peach Genome
Hotspot regions detected on chromosomes 1, 3, 4, 5, 6, and 8 controlled several fruit quality traits. The detection of hotspots in the genome indicates that genes related to certain traits are more densely concentrated in certain genomic regions . The main hotspot on chromosome 4 (9.0-12.5 Mbp) included reliable QTNs for DAB, FW, RP, ADH, RD, FF, and SSC detected in different seasons and/or approaches and represents a target region for future breeding studies in peach. A QTL hotspot associated with quality traits was previously reported in peach on chromosome 4 (Cantín et al., 2010;Eduardo et al., 2011). However, the study was performed using SSR markers and the QTLs were detected in low-density linkage maps. Using highdensity SNP maps, Rawandoozi et al. (2020a) reported QTLs for DAB and RD within the genetic region detected. Likewise, Hernández Mora et al. (2017) detected a hotspot for blush, SSC, RD and DAB in a wider genetic interval located at 11.2-14.1 Mbp in European germplasm. Moreover, Desnoues et al. (2016) identified a QTL hotspot in the same chromosome related to individual sugars and FW, although in a different location. In addition, the hotspot on chromosome 5 (0.3 to 3.7 Mbp) matched with the QTL hotspot for SSC and TA reported by Hernández Mora et al. (2017). Therefore, this study reinforces the importance of breeding programs targeting the improvement of fruit quality traits in peach focusing on the chromosome 4 and also demonstrated the necessity to promote further studies for the hotspot regions in chromosome 1, 3, 5, 6, and 8.

Candidate Genes
Candidate genes (566) were identified within the haploblock regions encompassing the QTNs detected using the mrMLM 4.0 and FarmCPU, and the GO enrichment approach narrowed down the initial CG list (222) and revealed over-representation of certain GO terms (68). RNA binding proteins and serinetype endopeptidase inhibitor-related genes were identified, and previous studies revealed involvement in the regulation of flowering time (Steffen et al., 2019;Zhang F. et al., 2019). In addition, genes functionally annotated as 2-oxoglutaratedependent dioxygenase, drug transmembrane transport, antiporter activity, pyridine nucleotide biosynthetic process, and chromatin assembly or disassembly were associated with fruit ripening in tomato, apricot, grape, peach, apple, and strawberry (Hanana et al., 2007;Farinati et al., 2017;Decros et al., 2019;Ding et al., 2020;García-Gómez et al., 2020). Furthermore, previous studies have shown that molybdopterin cofactor plays an important role in the metabolic control of avocado fruit growth and final fruit size (Cowan et al., 2001) and ion/H+ exchanger genes (GO: ion transmembrane transport) were critical for providing pH regulation (Pittman, 2012). Lastly, among the CG, Prupe.4G262200 and Prupe.4G261900 coding for endopolygalacturonases (GO: polygalacturonase activity) were previously involved in the inheritance of fruit texture and flesh adherence to the stone in peach (Peace et al., 2005;Gu et al., 2016).

CONCLUSION
We successfully performed a multi-locus GWAS using mrMLM 4.0 and FarmCPU in 620 individuals from three public fresh market peach breeding programs. A total of 88 reliable QTNs were consistently detected in at least two seasons and/or in different methods. Hotspots for quality traits were identified on chromosomes 1, 3, 4, 5, 6, and 8. Candidate genes for quality traits were identified in the vicinity of the reliable QTNs detected using mrMLM 4.0 and FarmCPU. Furthermore, we observed that the position of the previously reported candidate genes for fruit-related traits (BD, Blush, DAB, ADH, RP, pH, and TA) matched with the position of the hotspots detected on chromosomes 1, 3, 4, and 5. Therefore, the information reported in this study supports the development of DNA tools for MAS in peach. Moreover, the importance of chromosome 4 hotspot in breeding for improvement of fruit quality is reinforced, and also emphasized the necessity to further study the hotspot regions on chromosomes 1, 3, 5, 6, and 8.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: www.rosaceae.org, tfGDR1048/b.

AUTHOR CONTRIBUTIONS
CS: formal analysis and writing -original draft. LC: SNP data curation and review. WF: candidate genes analysis. MW, JC, and DB: resources and writing -review & editing. ZR: phenotypic analysis and review. KG: conceptualization, funding acquisition, resources, supervision, and writing -review & editing. All authors have read and approved the final manuscript.