GWAS provides new insights into the genetic mechanisms of phytochemicals production and red skin colour in apple

Abstract Understanding the genetic architecture of apple phytochemicals, and their interplay with conventional selection traits, is critical for the development of new apple cultivars with enhanced health benefits. Apple accessions (n = 344) used for this genome-wide association study (GWAS) represented the wide diversity of metabolic profiles in the domesticated and wild Malus genepools. Fruit samples were phenotyped for 34 metabolites, including a stable vitamin C glycoside “ascorbic acid 2-β-glucoside” (AA-2βG), and the accessions were genotyped using the Apple 20 K SNP Array. Several fruit quality traits, including red skin over-colour (OCOL), were also assessed. Wild Malus accessions showed at least 2-fold higher average content of several metabolites (e.g. ascorbic acid, chlorogenic acid, phloridzin, and trilobatin) than Malus domestica accessions. Several new genomic regions and potential candidate genes underpinning the genetic diversity of apple phytochemicals were identified. The percentage of phenotypic variance explained by the best SNP ranged between 3% and 21% for the different metabolites. Novel association signals for OCOL in the syntenic regions on chromosomes 13 and 16 suggested that whole genome duplication has played a role in the evolution of apple red skin colour. Genetic correlations between phytochemicals and sensory traits were moderate. This study will assist in the selection of Malus accessions with specific phytochemical profiles to establish innovative genomics-based breeding strategies for the development of apple cultivars with enhanced nutritional value.


Introduction
Apples are among the most consumed fruit worldwide, and studies have shown a clear trend for higher polyphenol content in whole apples leading to greater overall health benefits [1]. Apple polyphenols can be extracted and used as functional food materials that aid in weight management, satiety and obesity prevention. Consumption of apple polyphenol extract was found to have an inhibitory effect on the accumulation of visceral fat in humans [1]. The apple phytochemicals are widely recognised for their antioxidant activity, suggesting that apple or constituents of apple may contribute to increased total antioxidant capacity in human plasma. Various f lavonoids compounds are found in apple, and studies in humans using f lavonoid-containing extracts showed enhanced cognitive function and memory performance [2].
The visual characteristics (e.g. skin characteristics), eating quality (e.g. texture and f lavour), and storability are among the main fruit quality traits being targeted in apple breeding programmes, but the enhancement of phytochemicals is now gaining traction to select "bio-fortified" apple cultivars. Farneti et al. [3] screened various accessions of Malus domestica and wild Malus species and reported that modern apples have drastically reduced polyphenol content compared with the ancestral heritage and germline cultivars. Davies et al. [4] found that the phenolic content in cultivated apples was 68% lower than in their wild progenitors (Malus sieversii), and suggested that selection against internal f lesh browning and bitterness could have contributed to the lower phenolic content in commercial cultivars. There is an opportunity to utilise the diversity of wild apples to enhance the nutritional content of future cultivars.
QTL mapping studies, using relatively small bi-parental populations (ca. 150) of M. domestica, have identified the upper region of chromosome 16 as a potential hotspot for apple phytochemicals [7,8]. In a genome-wide association study (GWAS) on fruit samples of 136 M. domestica accessions, McClure et al. [9] reported that the concentrations of polyphenols varied by multiple folds across cultivars, and that most polyphenols could be predicted using DNA markers. In order to capture genetic diversity across species, a recent GWAS using 124 accessions of M. domestica and Malus spp. identified marker-trait associations for phytochemicals on each linkage group [10]. Apples with higher ascorbic acid (AsA) content would not only provide additional health benefits, but would also enhance the storage and shelf-life of the fruit. The recent discovery of a more stable glycosylated version of AsA, "ascorbic acid 2β-glucoside" (AA-2βG), in wild Malus species offers opportunities to develop new cultivars with all the associated health benefits of vitamin C [5]. Several studies have been conducted to identify QTL for the AsA biosynthesis [6], but apparently there are no published reports on the genetic control of AA-2βG in fruit crops.
Studies on apple phytochemicals have reported on the genetic variation [3,11] and marker-trait associations [9], but the genetic relationships between apple fruit quality traits (especially sensory traits) and f lavonoids have not been reported. The complex quantitative nature of fruit quality traits and phytochemicals can be affected by genetic and non-genetic factors, hence an understanding of the genetic interplay between these two groups of traits would be vital for parental selection and apple cultivar breeding. In the present GWAS, we used by far the largest sample set (ca. 350) of M. domestica and wild Malus species accessions to assess the diversity and genetic architecture of main bioactive compounds of apple fruit. Marker-trait association analyses were conducted, and the relationships between sensory traits and phytochemicals were also investigated to guide the development of bio-fortified apple cultivars.

Variation of phytochemicals in apple germplasm
The targeted LC-MS analyses of fruit samples revealed the presence of 34 compounds in the accessions of M. domestica and wild species. Seasonal variation was observed for most phytochemicals for the 15 accessions assessed in both years, with between-year correlations ranging from 0.22 for fatty acids to 0.96 for dihydrochalcones, with an average of 0.65 across all 34 compounds.
While most of the 34 compounds showed substantial variation among the 344 accessions (Supplementary Figure S1), only 14 were significantly (p < 0.01) different between Malus spp. and M. domestica accessions ( Table 1). The concentrations of AsA, chlorogenic acid, epicatechin, phloretin, phloridzin, procyanidins B2, B7 and C1 were significantly higher in Malus spp. accessions by two to three fold ( Table 1). The maximum phloridzin content assessed in wild Malus spp. accessions was 430 μg/g FW compared with a maximum of 115 μg/g FW in M. domestica. Similarly, the highest chlorogenic acid (CGA) content (2746 μg/g FW) in Malus spp. accessions was about 2.5-fold higher than the maximum of 1104 μg/g FW in cultivated apples. The average content of triterpenoids (e.g. maslinic, betulinic, ursolic and corosolic acid) in M. domestica accessions was 11, 8, 60 and 21 μg/g FW, respectively, which did not differ significantly from the wild species (Table 1).
AA-2βG was detected in two-thirds of wild species accessions, with an average of 927 μg/g FW. On the other hand, only about 2% of M. domestica accessions displayed AA-2βG with a very small average content (< 1.0 μg/g FW). The average content of malic acid (MA) in M. domestica accessions was 8682 μg/g FW, which was significantly lower compared with wild accessions (11504.3 μg/g FW) ( Table 1). Trilobatin was detected in the fruit of 9 out of 31 Malus species accessions, and the content in "Malus f loribunda 821" (33.1 μg/g FW) and Malus sieboldii "Snowbright" (29.9 μg/g FW) was about four times higher than in the M. trilobata accession (7.93 μg/g FW). Trilobatin was not considered for GWA analysis due to the low detection frequency. The phytochemicals-based clustering (Fig. 1A) of domesticated and wild species accessions overlapped as many of the compounds were not significantly different between the two genepools. The PCA loadings plot (Fig. 1B) suggested that the accessions of wild species were characterised by a higher content of various compounds, including procyanidin B2, procyanidin B7, procyanidin C1, phloretin, epicatechin, and chlorogenic acid -supporting the results described in Table 1.

Genetic correlations among phytochemicals and fruit quality traits
The degree of the correlations among all pairs of phytochemicals are illustrated in Fig. 2 and Supplementary Figure S3 Although several pairs of compounds were negatively correlated, only few pairs (e.g. corosolic acid and MA; trans-5-p-coumaroyl quinic acid and pomaceic acid) were found to be significant (p < 0.001). Maslinic and corosolic acids showed strong positive correlation with several other compounds (e.g. annurcoic acid and 3-oxo-hydroxy-urs-12-en-28-oic acid_1). The correlation (ca. 0.40) between AsA and AA-2βG was significant, and both of these compounds were also significantly correlated with MA (Fig. 2). The clusters displayed by the dendrogram showed that the f lavanols (procyanidins, catechin, and epicatechin) and flavonols (quercetin compounds) clustered separately. Similarly, several triterpene acids, such as pomaceic, annurcoic, euscaphic, pomolic, maslinic, betulinic and ursolic acids also clustered together (Fig. 2).

Marker-trait associations for phytochemicals
There were multiple genomic regions significantly associated with various phytochemicals on all 17 chromosomes, and the number of significant genomic regions was the highest on chromosome 9 and the lowest on chromosomes 8 and 14 ( Fig. 4; Supplementary Data Table S3). The percentage of phenotypic variance explained by the most significant SNP for each phytochemical ranged between 4% (triterpenes) and 21% (procyanidin) (Supplementary Figure S2).

Fruit quality traits
Three main genomic regions associated with OCOL were identified on chromosomes 13, 16 and 17 (Fig. 5)

Phytochemicals variability and between-trait relationships
Similar to earlier studies [3,12], our results revealed large variation within-and between the M. domestica and wild Malus species accessions. Compared with the other metabolite groups, differences in organic acids and DHC content between M. domestica and Malus spp. were more prominent. The fruit from wild accessions also showed significantly higher content of AsA and AA-2βG than M. domestica accessions.
The genomic h 2 of AA-2βG was considerably higher than that of AsA (0.42 vs 0.20), and the between-season correlation for AA-2βG was also higher compared with AsA (0.95 vs 0.70). These findings would suggest that AA-2βG could be a better selection trait than AsA for improving Vitamin C content of future apple cultivars, but further research is needed on the stability of AA-2βG in the human gut and whether this will increase vitamin C bioavailability. Dihydrochalcones (e.g. phloridzin), f lavanols (e.g. catechin, epicatechin) and hydroxycinnamic acids (e.g. chlorogenic acid), which are associated with several health benefits [1], tended to be slightly negatively correlated with the overall eating quality of apple fruits in this study, which would pose challenges for the simultaneous improvement of fruit quality and these metabolites.
Trilobatin has shown strong antioxidant, anti-inf lammatory, and α-glucosidase inhibition activities for the management of diabetes [13]. In our study, trilobatin content was favourably correlated (∼0.25) with sensory sweetness, suggesting that introgression of trilobatin in future cultivars would further enhance both health benefits and taste profiles of apples. We are conducting further studies, using segregating populations, to better understand the genetics of trilobatin content in apple fruit. The finding of a significant correlation between OCOL and ursolic acid in our study is supported by earlier reports of simultaneous production of anthocyanin and ursane/oleanane in plant cell culture [14]. However, the physiological and molecular functions underpinning the co-production of these compounds remain to be determined in apples.

Marker-trait associations and potential candidate genes for phytochemicals
Several ripening and senescence-associated gene functions (e.g. MdEIN3: MD15G1441000, MdERF1: MD01G1177000, ACS12: MD01G1186400, MD09G1239100, MD17G1231600) were located close to the significant SNPs, supporting earlier studies showing linkage between the ripening and phytochemical biosynthesis pathways [15,16]. The EIN3-like protein was also reported to be involved in terpene synthases in kiwifruit [15], and studies on winegrape [17] revealed a potential crosstalk between f lavonol and terpenoid metabolism. This would support the occurrence of a terpenoid cyclase protein (MD11G1284300) about 75 kb from the quercetin-linked SNP.
The occurrence of several UGT proteins in close proximity to the DHC-associated regions in this study is in agreement with earlier reports on the role of UGTs that convert phloretin to phloridzin in Malus species [18]. The DHC-associated regions on Chr1 and Chr5 overlapped with previous studies [7][8][9]. Chalconef lavanone isomerase (CHI) is one of the enzymes involved in the quercetin biosynthesis pathway, and several CHI genes (MD01G1118000, MD01G1118100, etc.) were located in the region associated with quercetin glycosides on chromosome 1, and this region overlaps with a previously published QTL [7,8]. Overexpression of CHI genes was shown to increase f lavonols in "Royal Gala" fruit [19].
This is the first study to identify the upper region (1.80-3.60 Mb) of Chr9 as a major hotspot for f lavanol regulation. This region harbours MYB transcription factors MYB24 (MD09G1030700) and MYB16 (MD09G1054000). Transgenic assays showed that MYB16 affects the regulation of anthocyanin and f lavanol in apples [20,21], and this gene has also been reported as a major regulator of cuticle formation in apples [22]. MYB24 was reported to promote f lavonoid biosynthesis through regulation of FLS1 gene expression in Arabidopsis thaliana [23], and MYB24 also plays a role in f lavonoid regulation in grape berries [24]. The MYB transcription factor MYB101 (MD09G1109900), underpinning a significant region, was shown to play role in chalcone synthase (CHS) leading to f lavonoid biosynthesis in Lotus japonicus [25] and Rosa rugosa [26].
The occurrence of MYB52 (MD06G1235300) and MYB3 (MD15G1411200) within genomic hotspots for phloridzin and CGA is supported by previous reports. For example, overexpression of MdMYB3 in tobacco resulted in upregulation of several f lavonoid pathway genes, including CHS, CHI, UFGT and FLS [27]. The expression of MYB52 was reported higher, along with the content of phloridzin, in russeted apples, suggesting a potential indirect role of this drought-stress response gene in phloridzin biosynthesis [22]. MYB26 (MD09G1191100), underpinning the GWAS signals for CGA on chromosome 9, was shown to be involved in the regulation of phenylpropanoid production in pea [28].
The HCT/HQT and 4CL genes (MD17G1225100 and MD17G1229400) resided in the DHC (3-hydroxyphloridzin) associated genomic regions, and these gene functions were hypothesised to play a role in the regulation of CGA [9] and phloridzin [29]. HCT/HQT are also reportedly linked to chalcone and f lavonoid synthesis in plants [30], and a cluster of HCT/HQT proteins f lanked a f lavanol-associated signal at the distal end of Chr9 which is paralogous to the DHC-associated region on Chr17.
Our GWA results suggest that a CCOAMT gene (MD16G1118200), CAD genes (MD17G1119700; MD17G1119800) and DHQS gene (MD02G1126200) are potential candidates underpinning DHCassociated genomic regions. These gene functions were also proposed as candidate genes for the biosynthesis of hydroxycinnamic acids (CGA) [9], which is supported by the high genetic correlation (ca. 0.65) between phloridzin and CGA in this study. The CCOAMT enzyme was reported to play a key role in the synthesis of quercetin in citrus [31], which gives support to the occurrence of a cluster of CCOAMT genes in the vicinity of a quercetin-associated SNP in the upper region of chromosome 16 (chr16:8480114) in this study.
Overexpression of B-box (BBX) genes in apple callus led to the higher expression of anthocyanin biosynthesis-related genes (MdMYB1, MdANS) and f lavonoid pathway genes (MdCHS, MdF3H and MdDFR) [32], which support the presence of a BBX protein (MD10G10660000) in the f lavanol-associated region on Chr10 in our study. A significant signal for f lavanol compounds in the central region of chromosome 9 was underpinned by 2ODDS gene (MD09G1177300), and this family of gene systems was shown to be involved in the f lavonoid biosynthesis pathway in tobacco [33].
Clusters of 2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase proteins occurred within the f lavonols-linked regions on Chr1 and Chr16. 2OG genes are involved in the synthesis of anthocyanin in crabapples [34], and it was shown that overexpression of anthocyanin biosynthesis genes promoted transcription of some f lavonoid biosynthetic genes and induced quercetin accumulation in crabapples [21]. A study on Arabidopsis thaliana also reported that f lavonol synthase requires 2OG and ferrous iron (Fe) for catalytic function [35]. The concentration of f lavonols was reported negatively correlated with laccase activity in apple cultivars [36], and interestingly the f lavonols-associated hotspot on Chr1 harbours two laccase proteins (MD01G1158900 and MD01G1159400) in our study.

Putative candidate genes for fruit quality traits
We identified a novel texture-associated signal on Chr9 in this study, and the MYB6 (MD09G1261200), located within 50 kb of the significant SNP, has been shown to play a role in ripening and storability of strawberry [37]. Some other gene functions residing in this region include a senescence-associated protein (MD09G1263800), the S-adenosyl-L-methionine-dependent methyltransferases (SAM-MTase) protein (MD09G1258900), and a dehydroascorbate reductase (DHAR) protein (MD09G1264700). SAM-MTase acts in the biosynthesis pathway of ethylene and polyamines, and plays important roles during the ripening process and may affect fruit firmness [38]. A DHAR protein, known for its role in limiting oxidative stress, was shown to be upregulated in highly firm apple fruit, suggesting a role in maintaining fruit texture [39].
Our study identified novel genomic regions associated with OCOL variation in apple germplasm. Among other gene functions, paralogs of MYB15 and ERF1B underpinned the association signals for OCOL in the syntenic regions on Chr13 and Chr16 -suggesting that whole genome duplication has played a role in the evolution of apple red skin colour. The levels of anthocyanin significantly increased in MdERF1B-overexpressing callus [40,41], while the expression of MYB15 was shown to be upregulated in "Royal Gala" skin under hot conditions, suppressing the regulation of anthocyanin [20]. Further studies are needed to understand the interplay between these two contrasting anthocyanin regulatory genes ERF1B and MYB15.
There were two WRKY transcription factors (WRKY75: MD17G1001500; WRKY28: MD17G1138100) in the OCOL-associated regions on Chr17. Both of these WRKY transcription factors have been reported being involved in the regulation of anthocyanin biosynthesis [42,43]. It has been suggested that the acetyl-CoA pool could be directed to anthocyanin biosynthesis during many different stresses [44], which supports the occurrence of the acetyl-CoA synthase gene (MD17G1003300) within 3 kb of the OCOL-associated SNPs in the upper region of Chr17. NAC42, supporting the OCOL-associated region in the middle of Chr17, regulates anthocyanin accumulation in apple by interacting with MYB10 [45]. Our results suggest that the genetic regulation of apple red skin is complex and potentially controlled by the interplay between several gene families.

Breeding for phytochemicals
The content of several key phytochemicals was significantly higher in wild apple accessions, suggesting that these accessions could be used as parents for introgression breeding to develop phytonutrient-rich apple cultivars. In this germplasm study, the correlations between eating quality traits and phytochemicals were moderate, which suggested that it would be feasible to develop tailored breeding strategies for the simultaneous improvement of these two groups of traits. Studies using segregating populations would be helpful to further understand the genetic correlations between eating quality and phytochemicals, which would guide future cultivar development strategies. The complex polygenic architecture of most metabolites in this study suggested that genomic selection would be a better option to reduce the generation interval for introgression breeding [46].

Plant material
The ex situ Malus germplasm collection maintained at the Hawke's Bay site (39 • 40'S 176 • 53 E) of The New Zealand Institute for Plant and Food Research Limited was established as grafted trees in late 1950s when budwood from heritage cultivars and crab-apples was imported mainly from the United Kingdom. In 2013, each accession was re-propagated onto "M9" rootstock and maintained in the orchard following standard commercial practices. This study was performed on 344 accessions, comprising 313 M. domestica accessions and 31 accessions of wild Malus species (Supplementary Data Table S1).

Metabolite phenotyping
Metabolite phenotyping was conducted over a period of two fruiting seasons. Of the 344 accessions, 165 and 179 were assessed in 2019 and 2020, respectively. There were 15 accessions common between the two years to estimate seasonal effect and obtain seasonally adjusted phenotypes of all accessions.

Sample preparation
We followed fruit harvesting and storage protocols as reported earlier by Kumar et al. [47]. Brief ly, the fruit were stored for four weeks at 0.5 • C, then a further 1 day at 20 • C before evaluation. For each accession, a total of 16 cores were cut from four apples at four positions around the equatorial circumference of each apple using a #5 corer (9 mm internal diameter). Each core contained apple skin and cortex (but no apple core or seeds) and was trimmed to approximately 1 cm if required. The total weight of the 16 cores was approximately 10 g and the actual weight was recorded (two decimal places). The 16 apple cores were added to a 250 mL Schott bottle containing 100 mL solvent (ethanol:water:formic acid 80:20:2 with the internal standard naringenin at 1.05 mg/L) and homogenized using the Omni General Laboratory Homogenizer (GLH850, Kennesaw, GA, USA). Samples were left overnight at 4 • C to extract and then an aliquot was transferred to a centrifuge tube for analysis.
Analysis of ascorbic acid (AsA) and ascorbic acid 2-β-glucoside (AA-2βG) by UHPLC Immediately after extraction, a portion of the extract was diluted 1:1 with 2 mM tris(2-carboxyethyl) phosphine hydrochloride (TCEP.HCl) (Sigma-Aldrich, Sydney, Australia) in MeCN/formic acid/water (90:5:5) to preserve the AsA. The actual volume used was 200 μL extract +200 μL 2 mM TCEP in MeCN/formic acid/water (90:5:5) directly into LC vials. The UHPLC system Dionex Ultimate ® 3000 Rapid Separation LC was composed of a pump, autosampler, column oven, and Photodiode Array (PDA) detector. The LC column was a Luna HILIC 150 mm x 4.6 mm, 3 μm (Phenomenex) and was maintained at 40 • C. The f low was 1000 μL min −1 . The solvents were A = 20 mM HCOONH 4 , pH = 5.0 and B = 100% acetonitrile. The solvent composition was isocratic at 20% A and 80% B. The run time for each analysis was 10 min. The injection volume for samples and standards was 3 μL. Data were processed using Chromeleon (ThermoScientific). The response was calibrated using AsA and AA-2αG, therefore AA-2βG was quantified as AA-2αG equivalents. It is noteworthy that AsA is unstable and can degrade during the analytical process. The sample preparation methods were not optimised specifically for AsA, so the measurements of AsA in this study are estimates only. In contrast, AA-2βG is a very stable form of Vitamin C [5], and reliable quantitative results were obtained.

Analysis of polyphenolic metabolites by LC-HRAM-MS
Each sample extract was diluted 1:1 with methanol (actual volumes: 200 μL extract +200 μL methanol) and then analysed for polyphenols by Liquid Chromatography -High Resolution Accurate Mass -Mass Spectrometry (LC-HRAM-MS) as described earlier [5]. LC-MS grade acetonitrile (Supelco, Bellefonte, PA, USA), ethanol, methanol and formic acid were purchased from Thermo Fisher Scientific (Waltham, MA, USA). Ultrapure water was obtained from a Milli-Q Intergral3 system (Millipore, Burlington, MA, USA). The LC column was a Luna Omega C18 100x2.1 mm, 1.6 μm (Phenomenex, Torrance, CA, USA) and was maintained at 40 • C. The f low was 400 μL min −1 . The solvents were A = 0.2% formic acid and B = 100% acetonitrile. The solvent gradient was: initial composition 90% A 10% B, 0-0.5 min; linear gradient to 60% A 40% B, 0.5-7 min; linear gradient to 5% A 95% B, 5-12 min; composition held at 5% A 95% B, 12-15 min; linear gradient to 90% A 10% B, 15-15.2 min; to return to the initial conditions before another sample injection at 18 min. The injection volume for samples and standards was 1 μL. The micrOTOF QII parameters were same as those described earlier [5]. Data were processed using TASQ (Bruker Daltonics, Bremen, Germany). The chemical data for the metabolites measured in this study are presented in Supplementary Table S2.

Fruit quality traits
Thin wedges of cortical tissue were cut from each of the two apples from each accession for sensory evaluation by one trained assessor. Following Kumar et al. [47], several sensory traits including, firmness (FIRM), crispness (CRISP), juiciness (JUIC), sweetness (SWET), sourness (SOUR), astringency (ASTR), f lavour intensity (FINT) and overall eating quality (EQUA) were assessed on a scale from 0 (= lowest) to 9 (= highest). Red skin over-colour (OCOL) coverage was scored visually on a scale 0 (no red colour) to 9 (completely red). Average fruit weight (AVFW) was measured as the average weight (g) of the six fruits per accession. The f lesh browning disorder (FBD) was also recorded as reported earlier by Kumar et al. [51]. As fruit of some Malus species accessions were very small and/or had an unpleasant taste, such accessions were assessed only for AVFW and OCOL.

SNP array genotyping and quality control
Leaf DNA from all 344 accessions were genotyped using the 20 K Infinium SNP array [48], and SNP calling was conducted using Genome Studio (Illumina Inc., San Diego, CA, USA). SNPs with minor allele frequency (MAF) <0.02, and missing data frequency > 10% were discarded, and the missing genotypes at the retained SNPs were imputed using k-nearest neighbour method. The coordinates of the filtered SNPs were identified on the GDDH13 genome assembly [49] by blasting their f lanking sequences to the GDDH13 genome and detaining the best hits with the aligned length matching the f lanking sequence length and having ≥90% identity.

Data analysis
Exploratory analyses, including trait distribution, correlation between traits, and principal component analysis (PCA) were conducted using the R software (R Core Team, 2019). The twosample t-test was used to test whether the population means of two groups (M. domestica and Malus spp.) are equal or not. The genetic structure of the 344 accessions was estimated using PCA, Neighbour Joining (NJ) tree, and model-based Bayesian clustering using ADMIXTURE software. Estimates of linkage disequilibrium (LD: r 2 ) were calculated between SNPs located on the same linkage group (LG) using R package "LDcorSV". The estimates of pair-wise LD between SNPs were averaged in the increments of one kilobase (kb) distance, and the LD decay curve was fitted using a standard logarithmic function.
The best linear unbiased estimates (BLUEs), adjusted for year effect, were used as phenotypes for GWAS. A mixed linear model (MLM) approach, accounting for cryptic relatedness (K-matrix) and population stratification (Q-matrix), implemented in R package GAPIT [50], was used for genotype-phenotype association analyses. A genome-wide significance threshold of α = 0.01 was used to identify significant genotype-phenotype associations. The putative candidate genes within 500 kb upstream and downstream of the most significant SNPs were retrieved using GDDH13 genome v.1.1 Jbrowser (https://www.rosaceae.org/ jbrowse/). Genomic best linear unbiased prediction (GBLUP) of the additive genetic value of all accessions were obtained for each trait, and then the product-moment correlations between GBLUP for different traits were used as estimates of genetic correlations. Following Kumar et al. [51], the estimates of genomic heritability (h 2 ) were also calculated for each trait using R package GAPIT.