Horticulture Research, 2022, 9: uhac028 https://doi.org/10.1093/hr/uhac028 Article Genetic architecture and genomic predictive ability of apple quantitative traits across environments Michaela Jung1 , 2 , *, Beat Keller1 , 2 , Morgane Roth1 , 3 , Maria José Aranzana4 , 5 , Annemarie Auwerkerken6 , Walter Guerra7 , Mehdi Al-Rifaï8 , Mariusz Lewandowski9 , Nadia Sanin7 , Marijn Rymenants6 , 10 , Frédérique Didelot11 , Christian Dujak5 , Carolina Font i Forcada4 , Andrea Knauf1 , 2 , François Laurens8 , Bruno Studer2 , Hélène Muranty8 and Andrea Patocchi1 1 Agroscope, Breeding Research Group, 8820 Wädenswil, Switzerland 2 Molecular Plant Breeding, Institute of Agricultural Sciences, ETH Zurich, 8092 Zurich, Switzerland 3 GAFL, INRAE, 84140 Montfavet, France 4 IRTA (Institut de Recerca i Tecnologia Agroalimentàries), 08140 Caldes de Montbui, Barcelona, Spain 5 Centre for Research in Agricultural Genomics (CRAG) CSIC-IRTA-UAB-UB, Campus UAB, 08193 Bellaterra, Barcelona, Spain 6 Better3fruit N.V., 3202 Rillaar, Belgium 7 Research Centre Laimburg, 39040 Auer, Italy 8 Univ Angers, Institut Agro, INRAE, IRHS, SFR QuaSaV, F-49000 Angers, France 9 The National Institute of Horticultural Research, Konstytucji 3 Maja 1/3, 96-100 Skierniewice, Poland 10 Laboratory for Plant Genetics and Crop Improvement, KU Leuven, B-3001 Leuven, Belgium 11 Unité expérimentale Horticole, INRAE, F-49000 Angers, France *Corresponding author. E-mail: michaela.jung@usys.ethz.ch Abstract Implementation of genomic tools is desirable to increase the efficiency of apple breeding. Recently, the multi-environment apple reference population (apple REFPOP) proved useful for rediscovering loci, estimating genomic predictive ability, and studying genotype by environment interactions (G × E). So far, only two phenological traits were investigated using the apple REFPOP, although the population may be valuable when dissecting genetic architecture and reporting predictive abilities for additional key traits in apple breeding. Here we show contrasting genetic architecture and genomic predictive abilities for 30 quantitative traits across up to six European locations using the apple REFPOP. A total of 59 stable and 277 location-specific associations were found using GWAS, 69.2% of which are novel when compared with 41 reviewed publications. Average genomic predictive abilities of 0.18–0.88 were estimated using main-effect univariate, main-effect multivariate, multi-environment univariate, and multi-environment multivariate models. The G × E accounted for up to 24% of the phenotypic variability. This most comprehensive genomic study in apple in terms of trait- environment combinations provided knowledge of trait biology and prediction models that can be readily applied for marker-assisted or genomic selection, thus facilitating increased breeding efficiency. Introduction Apple breeding is labor- and time-intensive, but Apple (Malus × domestica Borkh.) is the third most pro- selection efficiency can be improved by integrating duced fruit crop worldwide [1]. Since its domestication DNA-informed techniques into the breeding process in the Tian Shan mountains of Central Asia, the culti- [7]. Marker-assisted selection allows breeders to predict vated apple developed into a separated near-panmictic the value of a target trait based on its association species [2]. Over the centuries, thousands of apple cul- with a genetic marker. The method leads to removal of tivars have been propagated and conserved thanks to inferior seedlings without phenotyping, thus increasing grafting [3]. Extensive relatedness among cultivars with selection intensity and/or reducing the labor costs when a strong inf luence of a few founders through the history decreasing the number of individuals passing to the of apple breeding has been reported despite their high next selection step [7]. Quantitative trait locus (QTL) genetic diversity [4–6]. Only a fraction of the existing mapping has been traditionally used to investigate the cultivars are grown commercialy [3] and they require an genetic basis of variation in traits such as pathogen intensive use of pesticides for crop protection. To diver- resistance, phenology, and some fruit quality traits [8–11]. sify apple production, it is desirable to produce new cul- To bridge the gap between the discovery of marker-trait tivars for sustainable intensive agriculture and adapted associations and their application in breeding, protocols to future climates, while remaining attractive to con- that transfer the knowledge obtained by QTL analyses sumers. into DNA tests were established [12, 13]. However, Received: 25 June 2021; Accepted: 11 January 2022; Published: 19 February 2022; Corrected and Typeset: 1 April 2022 © The Author(s) 2022. Published by Oxford University Press on behalf of Nanjing Agricultural University. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 2 | Horticulture Research, 2022, 9: uhac028 marker-assisted selection in apple remains restricted measurements from multi-environment trials, the effect to a limited number of traits associated with single of environment on individuals remains constant at genes or a handful of large-effect QTL, such as pathogen single locations, but the G × E leads to changes in the resistance and fruit firmness, acidity, or color [14]. DNA- ranking of genotypes across locations. With an increasing informed selection is rarely deployed in apple when proportion of G × E effect relative to genotypic effect, breeding for quantitative traits with complex genetic both heritability of average effect across environments architecture, though this task became feasible with the and response to selection decrease [40]. A noticeable recent technological developments in apple genomics. effect of contrasting European environments and G × E In the genomics era, advancements in genotyping on two apple phenology traits – floral emergence and and sequencing technologies led to a broad range of harvest date – has been reported, which demands testing new tools for genetic analyses. In the case of apple, the multi-environment modelling approaches in apple several reference genomes have been produced [15– [36]. A location-specific GWAS may be used to identify 19], single nucleotide polymorphism (SNP) genotyping loci with stable effects across environments and loci arrays of different densities such as 20 K or 480 K specific to individual locations [41]. Multi-environment SNPs have been developed [20, 21], and genotyping- prediction models can account for G × E by explicitly by-sequencing methods have been adopted [22, 23]. modeling interactions between all available markers Genome-wide association study (GWAS) emerged as a and environments [42]. These models can outperform method for exploring the genetic basis of quantitative more simple modelling approaches that ignore G × E [42– traits [24]. GWAS in apple has been used to identify 44]. Additionally, taking advantage of information that associations between markers and various traits such traits provide about one another, a multivariate (also as fruit quality and phenology traits [22, 23, 25–29]. called multi-trait) genomic prediction can be applied. The associations found with GWAS can be translated This method may be useful in case the assessment of one into DNA tests for marker-assisted selection. Besides trait remains costly, but another correlated trait with less GWAS, genomic selection was developed to exploit the expensive measurement is available or can be assessed effects of genome-wide variation at loci of both large more easily [45]. The multivariate prediction can also and small effects on quantitative traits using a single be extended to a multi-environment approach when model [30] and is sometimes called marker-assisted treating measurements from different environments as selection on a genome-wide scale [31]. For genomic distinct traits [46]. selection, prediction models are first trained with A population of 269 diverse apple accessions from phenotypic and genomic data of a training population. across the globe and 265 progeny from 27 parental In a second step, the models predict the performance combinations originating in recent European breeding of breeding material based on the genomic data alone. programs constituted our apple reference population These genomic estimated breeding values are then used (apple REFPOP) [36]. The apple REFPOP had a high-density to make selections among the breeding material, thus genomic dataset of 303 K SNPs and was deemed suitable increasing the breeding efficiency and genetic gain. for the application of genomics-assisted breeding [36]. Several studies have assessed genomic predictive ability Combined with extensive phenotypic information, the for apple quantitative traits related to fruit quality and apple REFPOP provided the groundwork for marker- phenology [22, 23, 29, 32–36]. Genomic selection can assisted and genomic selection across contrasting Euro- double genetic gain, as demonstrated by yield traits in pean environments. Hence, 30 traits related to productiv- dairy cattle [37], but the accuracy of genomic prediction ity, tree vigor, phenology, and fruit quality were measured for yield traits in apple has not been studied. Analyses of in the apple REFPOP during up to three years and at up to genomic datasets beyond 100 K SNPs have been limited six locations with various climatic conditions of Europe to flowering and harvest time (GWAS and genomic (Belgium, France, Italy, Poland, Spain, and Switzerland). prediction) [26, 36], fruit firmness and skin color (GWAS) GWAS was performed to dissect the genetic architecture [28, 38]. Marker density, trait architecture, and heritability of the studied traits, identify associated loci stable have been shown to differentially affect prediction across locations and location-specific loci, and to observe performance in simulated data and for apple [34, 36, 39] signs of selection on loci of large effect. Integrating and their impact on genomic analyses should therefore genetic diversity of the cultivated apple accessions and be further empirically tested. Moreover, GWAS for the progeny with phenotypic data measured in multiple same traits measured at different locations, the effect of environments, the goal was to estimate predictive ability genotype by environment interaction (G × E) on genomic and patterns of G × E for key traits in apple breeding. predictive ability, and predictions with multivariate Compared to previous genomic prediction studies genomic prediction models have not been evaluated yet carried out for a smaller number of genotypes and/or in apple. environments, this study aimed to obtain improved Plants are known for their strong phenotypic response predictive ability for productivity, tree vigor, phenology, to environmental factors, a phenomenon regularly and fruit quality traits using main-effect univariate, tested in plant breeding using multi-environment trials. main-effect multivariate, multi-environment univariate, In general, when statistical models are applied to and multi-environment multivariate genomic prediction Jung et al. | 3 models. Finally, a critical analysis of our results provided recommendations for future implementation of genomic prediction tools in apple breeding. Results Phenotypic data analysis The accession and progeny groups of the apple REFPOP were evaluated for 30 quantitative traits at up to six loca- tions. The measurements for ten traits were collected at one location, while the remaining 20 traits were available from at least two locations (three traits were measured in two locations, three traits in four locations, eleven traits in five locations and three traits in six locations, Figure 1, Supplementary Table 1). Most traits [25] were assessed during three seasons while five traits were measured during two seasons. The traits showed unimodal as well as multimodal distributions (Supplementary Figure 1). Differences of various extent between the accession and progeny groups were observed (Supplementary Figure 2). Removing environmental effects from the phenotypic data, best linear unbiased prediction of random effects of genotypes, hereafter called clonal values, were produced Figure 1. Locations and the respective number of phenotyping seasons for each trait. Locations of the measurements are labeled as: BEL – across all locations and separately for each location. Belgium, CHE – Switzerland, ESP – Spain, FRA – France, ITA – Italy, POL – As expected, high phenotypic clonal values correlations Poland. Traits measured at a single location are labeled with an asterisk. and genomic breeding values correlations (>0.7) between traits were observed within trait categories, namely the phenology, productivity, fruit size, outer fruit, inner fruit, Supplementary Table 3). No significant associations and vigor category (Figure 2a). A few moderate positive were observed for trunk diameter and russet cover in the phenotypic clonal values correlations (0.3–0.7) were across-location GWAS. In the location-specific GWAS, 309 found between trait categories such as harvest date and significant marker-trait associations for all 30 traits were fruit firmness (0.51), yellow color and russet cover (0.55), discovered (Figure 3b, Supplementary Table 3). Of these soluble solids content and russet cover (0.36), or between 309 marker-trait associations, 32 associations for twelve yield (weight and number of fruits) and vigor trait traits were shared between the location-specific GWAS category (0.36–0.51, Figure 2a). High average phenotypic and the across-location GWAS (Supplementary Table 3). clonal values correlations were observed between the The coefficient of determination (R2 ) of significant environments (combinations of location and year) for associations was the largest for red over color (0.71), harvest date (0.82 [0.73, 0.95]) or red over color (0.80 [0.62, green color (0.55) and harvest date (0.42, Figure 3c, 0.92]) whereas low average phenotypic clonal values Supplementary Table 3). correlations (<0.3) were present between environments Significant associations with different traits co- for f lowering intensity (0.18 [−0.49, 0.68]) and trunk localized at identical positions or occurred very close increment (0.16 [−0.31, 0.55], Supplementary Table 2, in some genomic regions (distance between marker posi- Supplementary Figure 3). A shift of the progeny group tions below 100 kb, Figure 3c, Supplementary Table 3). compared to the accession group towards smaller, In the across-location GWAS, a marker significantly more numerous and less russeted fruits was observed associated with harvest date on chromosome 3 (posi- (Figure 2b). tion 30 681 581 bp) was located next to two markers associated with fruit firmness (positions 30 587 378 Genome-wide association studies and 30 590 166 bp). The same marker on the position Across-location GWAS for 20 traits measured at more 30 681 581 bp was also associated with harvest date, than one location (Figure 1, Supplementary Table 1) ground color, overall russet frequency and soluble and location-specific GWAS for all 30 traits were solids content measured at several different locations used to explore the genetic basis of the assessed (location-specific GWAS). Similarly, the association with traits. The quantile-quantile plots showed that the harvest date on chromosome 16 (position 9 023 861 bp) observed and expected distributions of p-values cor- was closely located to a marker associated with fruit responded well and no apparent inf lation of p-values firmness (position 8 985 888 bp) in the across-location was found (Supplementary Figure 4 and 5). Across- GWAS. The traits related to bitter pit analyzed in the location GWAS revealed 59 significant (−log10 (p) > across-location GWAS, i.e. bitter pit frequency and grade, 6.74) marker-trait associations in 18 traits (Figure 3a, showed significant associations on chromosome 16, 4 | Horticulture Research, 2022, 9: uhac028 Figure 2. Exploratory phenotypic data analysis of the studied quantitative apple traits. a Pairwise correlations between traits with the phenotypic clonal values and genomic breeding values correlations in the lower and upper triangular part, respectively. Phenotypic clonal values correlation was assessed as Pearson correlation between pairs of global clonal values (across-location clonal values with the addition of location-specific clonal values for traits measured at a single location), the genomic breeding values correlation as Pearson correlation between pairs of genomic breeding values estimated from a G-BLUP model. Trait categories are outlined along the vertical axis. Traits measured at a single location are labeled with an asterisk. b Principal component analysis biplot based on global clonal values. position 7 681 416 bp. Several associations with traits 2 005 502 bp on chromosome 8 was shared between measuring fruit skin russet in the across-location GWAS fruit diameter and fruit volume from Switzerland and co-localized on chromosome 12 (position 23 013 281 bp, single fruit weight from Belgium. On chromosome 11, russet frequency on cheek and in the eye) and 17 (posi- fruit diameter, fruit volume and single fruit weight tion 27 249 890 bp, overall russet frequency and russet from Switzerland, as well as single fruit weight from frequency at stalk). A marker at position 18 679 105 bp Belgium, shared the association at position 18 521 895 bp. on chromosome 1 was associated with both single Additionally, position 3 622 193 bp on chromosome 11 fruit weight from the across-location GWAS and fruit was shared between the associations of fruit length diameter from Switzerland (found with the location- and single fruit weight from Switzerland. For red over specific GWAS). The association with marker at position color and green color, the association with a marker Jung et al. | 5 Figure 3. Significant marker-trait associations found by GWAS. a Distribution of the significant associations and corresponding p-values from across-location GWAS over the 17 apple chromosomes. b Distribution of the significant associations and corresponding p-values from location-specific GWAS over the 17 apple chromosomes. Locations are labeled as BEL (Belgium), CHE (Switzerland), ESP (Spain), FRA (France) and ITA (Italy). a-b Size of the symbols indicate the −log10 (p). The x-axis shows chromosome numbers. c Physical positions (in bp) of the significant associations on chromosomes with their respective coefficients of determination (R2 ) from the across-location GWAS complemented with the location-specific GWAS for traits measured at a single location. Size of the symbols indicate the R2 . The x-axis shows chromosome numbers. 6 | Horticulture Research, 2022, 9: uhac028 Figure 4. Comparison of the significant marker-trait associations with previously published associations. a Venn diagram comparing the unique associations, which were either previously published (former), reported in the across-location GWAS (present) or the location-specific GWAS (present per location). Color intensity and the values ref lect the number of associations per diagram area. b Scatterplot of unique associations comparing published associations (former) with the merged across-location and location-specific GWAS (present). The traits were assembled into trait groups based on similarity between the approaches to the trait measurement. Symbol size ref lects the number of markers used in the studies. In case more than one publication reported an association in the same chromosome segment, only the report with the largest number of markers is shown (see Supplementary Table 4 for the complete list of previously published associations). a-b Positions of associations were assigned to three chromosome segments: top, center and bottom. Only the unique combinations of trait groups with segments and type of study (former or present) are shown. on chromosome 9 (position 33 799 120 bp) occurred in (Supplementary Figure 7). Among these major associ- across-location and four location-specific GWAS, while a ations, changes in the frequency of alleles with an close marker (position 33 801 013 bp, less than 2 kb away) increasing effect on trait phenotypes were quantified was associated in the two other location-specific GWAS. in 30 ancestral accessions (five ancestor generations of Additional significant marker-trait associations occurred the progeny group, Supplementary Table 5) and all 265 in the same genomic regions among the location-specific progenies included in the apple REFPOP (Figure 5a). Com- GWAS and between the across-location and location- pared to the ancestral accessions, the frequency of the specific GWAS (Supplementary Table 3). allele with an increasing effect on phenotype (Figure 5c) Previous reports on QTL mapping and GWAS in apple was higher in the progeny for the alleles associated with were extensively reviewed and 41 publications reporting later harvest date and increased flowering intensity, on traits measured similarly to our own were found and titratable acidity, fruit firmness and trunk increment taken for comparison (Supplementary Table 4). In the lit- (Figure 5a). For the marker associated with green color erature, in the across-location GWAS and in the location- and red over color, the allele frequencies were equivalent specific GWAS, 166, 52 and 172 unique combinations for ancestors and progeny, which reflected the minor of chromosome segments with traits were discovered, allele frequency of nearly 0.5 for both traits (Figure 5b,d). respectively (Figure 4a). Out of all segment-trait com- On a closer look at the allele frequencies across the binations across our GWAS, 30.8% overlapped with the accession and progeny groups for the markers closely previously published results of QTL mapping or GWAS. associated with harvest date and fruit firmness on chro- All previously published segment-trait combinations for mosome 3 (Figure 5b), the allele associated with later the trait groups bitter pit and trunk were also detected harvest date and firmer fruits was fixed in all progeny, in our study, whereas no overlap between the former while the allele with a decreasing effect on the phenotype and present associations was found for ground color and was present with a frequency below 0.1 in the whole sugar trait groups (Figure 4b, Supplementary Figure 6). apple REFPOP (Figure 5a-d). The allele associated with larger trunk increment on chromosome 1 was found in progeny known to segre- Allele frequency dynamics over generations gate for Rvi6, and it was present in only two accessions Eleven major significant marker-trait associations (‘Prima’ and X6398) that are also known to carry the (R2 >0.1) were identified in the global GWAS results apple scab resistance gene Rvi6, which is located about (across-location GWAS with the addition of location- 1.8 Mb from the SNP associated with trunk increment specific GWAS for traits measured at a single location (Figure 5b-c). only, Supplementary Table 3). Boxplots of the across- Squared Pearson’s correlations in a window of ∼3000 location clonal values against the dosage of the reference markers surrounding each of the major significant allele (0, 1, 2) for the eleven associations showed marker-trait associations showed that markers in linkage mostly additive effects of the alleles on phenotypes disequilibrium extended over larger distances around Jung et al. | 7 Figure 5. Allele frequency dynamics of the major significant marker-trait associations. a-d The associations were chosen based on the coefficient of determination (R2 >0.1) from the global GWAS. a For each association, frequency of the allele with increasing effect on trait phenotypes in the apple REFPOP is shown. For the progeny group (progeny) and its five ancestor generations (ancestors), the allele frequencies are shown as points connected with a line. Out of all known ancestors, the allele frequency was estimated for 30 accessions included in the apple REFPOP. Colors of the points and lines correspond to chromosome locations of the associated SNPs. b Allelic combinations carried by the apple REFPOP genotypes, sorted according to geographic origin of accessions and affiliation of progeny to parental combinations (the x-axis was labeled according to Supplementary Table 1 and 2 in Jung et al. [36]). c Global clonal values of traits and their standard error for each allelic combination, centered to mean 0 and scaled to standard deviation of 1. d Frequency of the minor allele in the whole apple REFPOP. b-d The legend and y-axis are shared between plots. In d, the color of an allele corresponds to the color of the homozygous allelic combination of the same allele in b and c. some marker-trait associations (Supplementary Figure 8). (G-BLUP), random forest (RF), BayesCπ and reproducing When visually compared with other loci, the associations kernel Hilbert spaces regression (RKHS)). with harvest date and fruit firmness on chromosome 3 as Similar performance of five main-effect models well as red over color and green color on chromosome 9 was observed (Figure 6a). The rt estimated for the were found in genomic regions of the highest linkage dis- baseline model G-BLUP varied between 0.28 and 0.78 equilibrium between markers (Supplementary Figure 8). (Supplementary Table 6). When the predictive ability of The markers associated with trunk increment and Rvi6 the G-BLUP model was averaged over all traits (r), the also showed signs of linkage disequilibrium among them obtained r was equal to 0.50. RF showed an rt higher (Supplementary Figure 8). than G-BLUP for 9 out of 30 traits and an r of 0.49. The remaining associations (R2 ≤0.1) reported by the BayesCπ , RKHS and the main-effect multivariate model global GWAS showed various trends in allele frequencies with an unstructured covariance matrix of the random across generations such as increased frequency of alleles effect (MTM.UN) showed an r of 0.50, 0.51 and 0.50 and associated with increased weight of fruits in the progeny exceeded rt of G-BLUP in one, twelve and ten traits, (Supplementary Figure 9). The individual parental com- respectively. Among the main-effect univariate genomic binations of the progeny group were often fixed for single prediction models, the rt ranged between 0.18 for russet alleles of the remaining associations (R2 ≤0.1) from the cover and 0.88 for red over color, both extreme values global GWAS (Supplementary Figure 10). observed with RF (Supplementary Table 6). Allele frequencies equivalent to those observed for the When compared with G-BLUP, the main-effect multi- whole progeny group were also found for the mean allele variate model MTM.UN showed an improved predictive frequency for 10-times repeated resampling of 30 geno- ability for several traits when they were modelled in com- types from the progeny group (Supplementary Figure 11). bination with a correlated trait (genomic breeding val- ues correlation larger than 0.3, Figure 6a, Supplementary Table 6). The inclusion of floral emergence as correlated Genomic prediction trait improved rt of full flowering from 0.43 to 0.46 and The best predictive ability across all eight models from 0.43 to 0.47 for end of flowering. Similar response compared in this study was found for the traits har- was observed for a combination of weight of fruits with vest date, green color and red over color (Figure 6, flowering intensity that improved rt of flowering inten- Supplementary Table 6). The lowest predictive ability sity from 0.58 to 0.59. Fitting the model using fruit length was found for traits related to bitter pit and russet as showed an increased rt of single fruit weight (difference well as yellow color. Additionally, the predictive ability in rt of 0.01) and using single fruit weight led to an for f lowering intensity and trunk increment with the increase in rt for fruit diameter, fruit length, maximum multi-environment models remained strongly below fruit size and fruit volume (difference in rt of 0.01, 0.03, the average predictive abilities per trait (rt ) of the 0.02 and 0.03, respectively). Using soluble solids content corresponding main-effect models (i.e. genomic-BLUP resulted in an increase of rt for russet cover (difference 8 | Horticulture Research, 2022, 9: uhac028 Figure 6. Genomic predictive ability in apple quantitative traits using eight genomic prediction models and two cross-validation scenarios. a Predictive ability of four main-effect univariate models, i.e. random forest (RF), BayesCπ , Bayesian reproducing kernel Hilbert spaces regression (RKHS) and genomic-BLUP (G-BLUP), and one main-effect multivariate model with an unstructured covariance matrix of the random effect (MTM.UN). The models were applied with a five-fold cross-validation where 20% of the genotypes were masked in each of the five runs. The MTM.UN was used in case a trait showed genomic breeding values correlation larger than 0.3 with at least one other trait. b Predictive ability of two multi-environment univariate models, i.e. across-environment G-BLUP (G-BLUP.E) and marker by environment interaction G-BLUP (G-BLUP.E.G × E), and the multi-environment multivariate factor-analytic model (MTM.FA). The models were applied under two five-fold cross-validation scenarios CV1 and CV2. The CV1 was applied for all traits using G-BLUP.E and G-BLUP.E.G × E and for traits measured in at least three environments using MTM.FA. The CV2 was applied for traits measured in Switzerland and in at least a one other location. a-b Predictive ability was estimated as a Pearson correlation coefficient between the observed and the predicted values of genotypes whose phenotypes were masked in a five-fold cross-validation. For the multi-environment models, the correlation coefficients were estimated for each environment separately. In the box plot, the bottom and top line of the boxes indicate the 25th percentile and 75th percentile quartiles (the interquartile range), the center line indicates the median (50th percentile). The whiskers extend from the bottom and top line up to 1.5-times the interquartile range. The points beyond the 1.5-times the interquartile range from the bottom and top line are labeled as dots. Jung et al. | 9 in rt of 0.01), while using russet frequency at cheek led to an improved rt of russet frequency at stalk (difference in predictive ability of 0.01). Predictive abilities for all possible combinations of correlated traits can be found in Supplementary Table 7. When comparing two multi-environment univari- ate models – across-environment G-BLUP (G-BLUP.E) and marker by environment interaction G-BLUP (G- BLUP.E.G × E) – and the multi-environment multivariate factor-analytic model (MTM.FA), the prediction perfor- mance of G-BLUP.E, G-BLUP.E.G × E and MTM.FA was generally lower under the first cross-validation scenario (CV1) than under the second cross-validation scenario (CV2, Figure 6b, Supplementary Table 6). For all traits, the G-BLUP.E.CV1, G-BLUP.E.G × E.CV1 and MTM.FA.CV1 showed lower rt than the main-effect G-BLUP, the r being equal to 0.40, 0.40 and 0.36, respectively. The G- BLUP.E.G × E.CV1 performed better than G-BLUP.E.CV1 for 14 out of 30 traits. The G-BLUP.E.CV2 and G- BLUP.E.G × E.CV2 outperformed G-BLUP for 13 out of 20 traits. The G-BLUP.E.CV2 and G-BLUP.E.G × E.CV2 both showed r equal to 0.57. The increase in rt from G- BLUP to G-BLUP.E.CV2 (0.35) as well as from G-BLUP to G-BLUP.E.G × E.CV2 (0.36) was the most pronounced for russet cover. The performance of G-BLUP.E.CV2 and G-BLUP.E.G × E.CV2 remained below the level of G-BLUP Figure 7. Synthesis of phenotypic and genomic analyses. predictions for productivity traits (f lowering intensity, Across-environment clonal mean heritability, genomic heritability, weight and number of fruits), ground color, soluble average predictive ability (rt ) for the main-effect G-BLUP and the proportion of phenotypic variance explained by the effect of each solids content, fruit firmness and trunk increment. The significantly associated marker (SNP 1–8), genotype (G), environment (E) G-BLUP.E.G × E.CV2 performed better than G-BLUP.E.CV2 and genotype by environment interaction (G × E). The significantly for 8 out of 20 traits. The r of MTM.FA.CV2 was equal to associated markers corresponded to results of the global GWAS. SNPs associated with each trait were sorted according to the proportion of 0.52 and therefore similar to G-BLUP, however, the model phenotypic variance they explained, i.e. SNP 1 represented the outperformed G-BLUP for nine out of 20 predicted traits association explaining the most variance within a trait. Proportions of (Supplementary Table 6). The MTM.FA showed higher phenotypic variance components were used to estimate clusters of traits outlined along the vertical axis. Within each cluster, the traits predictive ability than both G-BLUP.E and G-BLUP.E.G × E were sorted according to rt . for two traits under CV1 and five traits under CV2 (Supplementary Table 6). Synthesis of phenotypic and genomic analyses very low proportion of the total variance for floral emer- The across-environment clonal mean heritability was gence (13.1%), flowering intensity (11.4%), trunk diame- generally very high in the evaluated traits, the value ter (10.9%) and trunk increment (8.7%). The major pro- being close to one for harvest date and red over portion of the phenotypic variance was explained by the color and not lower than 0.80 for all the other traits effect of environment for floral emergence (73.9%) and with the exception of full f lowering (0.74), end of trunk diameter (66.3%). The lowest effect of environment flowering (0.79) and water core grade (0.79, Figure 7, was found for traits measured at only one location over Supplementary Table 6). The genomic heritability, which two or three years such as fruit diameter or water core is the proportion of phenotypic variance explained by frequency, both showing an effect of environment (i.e. the markers, was larger than 0.80 for harvest date, floral year) below 1%. The effect of G × E was the most pro- emergence, green color and red over color, the value nounced for productivity traits, i.e. flowering intensity was not lower than 0.40 for all the other traits with the (23.7%), weight of fruits (20.8%) and number of fruits exception of bitter bit frequency (0.33) and grade (0.39, (21.6%). The proportion of the G × E effect was the low- Figure 7, Supplementary Table 6). est for harvest date (4.7%), floral emergence (5.2%), red The effects of genotype and significantly associated over color (5.9%) and trunk diameter (4.2%) among the markers together explained a substantial part of the phe- traits measured at more than one location and for end notypic variance of traits, the largest sums of these geno- of flowering (5.7%), fruit volume (5.9%) and green color typic effects were observed for harvest date (82.8%) and (3.9%) among the traits measured at one location. A high red over color (74.6%, Figure 7, Supplementary Table 6). proportion of the phenotypic variance remained unex- Altogether, the sum of the genotypic effects explained a plained by the model parameters for flowering intensity 10 | Horticulture Research, 2022, 9: uhac028 (47.5%), bitter pit grade (53.4%) and trunk increment and fruit firmness, which can be also observed empir- (55.1%). ically for early cultivars that soften more, may be an Three clusters of traits were determined from the hier- example of pleiotropy or linkage disequilibrium between archical clustering of the table of traits by proportion loci. Harvest date and fruit firmness are known to be of phenotypic variance explained by different sources regulated by ethylene production [48] and associated (Figure 7, Supplementary Figure 12). A strong genotypic with loci present on chromosomes 3 (NAC18.1), 10 (Md- effect and a comparably low effect of environment and ACO1, Md-PG1), 15 (Md-ACS1) and 16 [22, 49–52]. G × E was observed for 13 traits assigned to the cluster In this work, closely located (distance <100 kb) asso- one. Most of the phenotypic variance was explained by ciations with both harvest date and fruit firmness were the effect of environment in f loral emergence and trunk found on chromosome 3. Migicovsky et al. [22] reported diameter, which were grouped in cluster two. Finally, 15 an overlap between associations with harvest time and traits with a pronounced effect of environment and/or fruit firmness on chromosome 3 falling within the coding G × E were grouped in cluster three. region of NAC18.1. The authors hypothesized that the lack of associations on other chromosomes was likely due to low SNP density around the causal loci (the study Discussion used a GBS-derived 8 K SNP dataset for 689 genotypes). The larger number of associations reported here might Both study aims were successfully achieved. Dissection be a result of the high SNP density (303 K SNPs) deployed of the genetic architecture of 30 key traits in apple in GWAS, however, not all previously reported loci were breeding and identification of associated loci stable re-discovered. across locations and location-specific loci using GWAS The SNPs associated with harvest date and fruit firm- provided 336 stable and location-specific marker-trait ness on chromosome 10 were further apart (∼6 Mb). For associations discovered across all traits under study. harvest date, one of the associations on chromosome 10 Of these loci, 69.2% were novel when compared with was stable across locations and several associations were previously published associations. Additionally, signs of location specific. However, the association on chromo- selection were found for associations of large effect. The some 10 with fruit firmness was found for the Italian multi-environment and multi-trait experimental design location only. It has been shown that chromosome 10 of the apple REFPOP allowed to estimate patterns of contains more than one QTL controlling fruit firmness G × E for the studied traits and assess their predictive [49–51], but stable across-location association with fruit ability. Our results showed that G × E accounted for firmness on chromosome 10 was missing in our study. up to 24% of the phenotypic variability and genomic One of the known loci on chromosome 10, the Md-PG1 predictive abilities of up to 0.88 were observed. Based on gene, is responsible for the loss of fruit firmness after the outputs of this research, recommendations for future storage [51, 53]. In apple REFPOP, fruit firmness was implementation of the prediction tools in apple breeding measured within one week after the harvest date and were derived. this very short storage period might have contributed to the less pronounced effect of the locus Md-PG1 in Discovered loci overlap between association our GWAS. studies and traits Two associations with harvest date measured in Italy Our GWAS enlightened the architecture of analyzed but no association with fruit firmness were found on traits as well as the identification of numerous marker- chromosome 15. Although a marker for Md-ACS1 related trait associations stable across, and specific to, the to ethylene production was previously mapped on chro- locations of the apple REFPOP. The particular design mosome 15 [50], and QTL for fruit firmness was dis- of the experiment, including the diversity of the plant covered on the same chromosome [49], these markers material used (accessions and small progeny groups), did not co-locate, but rather, mapped at the opposite multiple locations, and multiple years of evaluation, extremes of chromosome 15 [49, 50]. Likewise, the con- are likely to have resulted in about two thirds of the nection between harvest date and fruit firmness on chro- discovered associations being novel when compared with mosome 15 could not be confirmed here. the loci published in studies spanning more than two Our GWAS showed associations with harvest date and decades. Our study design also allowed us to replicate the fruit firmness on chromosome 16, which were located identification of many previously known loci associated 38 kb apart. In the past, loci associated with harvest date with the studied traits. and fruit firmness have been reported in the same region The association of one locus with two or more seem- on chromosome 16 [26, 49]. The role of this locus in the ingly independent traits (i.e. caused by pleiotropy) and regulation of harvest date and fruit firmness remains linkage disequilibrium between loci associated with dif- unknown and requires further research. ferent traits are frequent for complex traits [47]. The In practice, ripeness of fruit (harvest date) is decided GWAS performed in this study showed several marker- based on ground color and starch content. The GWAS trait associations at identical or close positions for dif- results showed that the association on chromosome 3 ferent traits. The interdependency between harvest date was not only found for harvest date and nearby markers Jung et al. | 11 associated with fruit firmness, but also corresponded to (Better3fruit, Belgium) than in parents used in the associations with ground color and soluble solids con- US breeding programs [11, 56]. The desired acidity tent. This might be explained by the fact that these level might depend on local climate of the breeding traits are used to define ripeness and thus harvest date. program and market preferences [56]. The increase in Further, the association of the NAC18.1 locus on chromo- frequency of the allele contributing to higher acidity in some 3 with overall russet frequency would support the the progeny may indicate a current preference towards known enhanced expression of NAC transcription factors more acidic apples in European breeding, but further in russet skin [54]. investigation is needed to clarify the trend. The last Co-localizations between associations found for dif- locus of large effect showing allele frequency dynamics ferent measures of bitter pit on chromosome 16, russet between generations was found for trunk increment. on chromosomes 12 and 17, fruit size on chromosomes The increased frequency of allele associated with an 1, 8 and 11, and skin color on chromosome 9 are likely increase in trunk increment may have occurred in the the result of relatedness among trait measurements. The progeny due to its potential impact on productivity measures that are easiest to score can be used in future suggested by moderate positive phenotypic clonal values to phenotype these traits. correlations between tree vigor (trunk diameter and increment) and yield-related traits. Additionally, the Signs of selection in marker-trait associations of marker associated with trunk increment was 1.8 Mb large effect apart from a SNP marker associated with Rvi6 gene The design of apple REFPOP allowed for the discovery responsible for resistance against apple scab [10]. These of major marker-trait associations and for the analysis two markers (AX-115183752 for trunk increment and AX- of changes in allele frequency between 30 ancestral 115182989 (also called Rvi6_42M10SP6_R193) for apple accessions and 265 progeny included in the apple scab) showed a correlation of 0.15 and occurred within REFPOP. Although the progeny group did not undergo a region of increased linkage disequilibrium between phenotypic selection before the establishment of the markers (Supplementary Figure 8). All accessions were apple REFPOP, the parents of the progeny were a matter homozygous for the reference allele of AX-115183752 of choice within the European breeding programs. This associated with decreased trunk increment (Figure 7c) allowed us to assess the impact of the past pheno- except for ‘Prima’ and X6398, which were heterozygous. typic selection on the studied germplasm. Comparing The scab-resistant accessions ‘Prima’ and X6398 (which ancestors with the progeny, higher frequencies of the is a second-generation offspring of ‘Prima’ [57]) but also alleles associated with later harvest date and increased ‘Priscilla-NL’ (known to be heterozygous for Rvi6 [58]), flowering intensity, titratable acidity, fruit firmness and were also heterozygous for AX-115182989. All other trunk increment were found for the progeny. Of these accessions were homozygous for the reference allele traits, harvest date and fruit firmness are correlated, not associated with Rvi6. The allele on chromosome 1 probably due to pleiotropy or linkage disequilibrium of associated with increased trunk increment may have causal loci, as it was shown in this and previous studies been co-selected with the Rvi6 locus responsible for [22]. Consequently, the consistently higher frequency of resistance against apple scab. alleles contributing to later harvest and firmer apples Signs of intense selection for red skin were recently in the progeny is because the softening of harvested detected in the USDA germplasm collection when com- apples is undesirable and likely selected against [55]. pared with progenitor species of the cultivated apple [5]. Signs of selection for increased firmness were also Our results show that the associations with red over recently found in USDA germplasm collection [5]. Our color and green color, which phenotypically mirrored red study also showed fixation of the late-harvest and high- over color and was associated with the same marker, did firmness alleles on chromosome 3 in the whole progeny not show changes in allele frequency between ances- group, which suggests a loss of genetic diversity in the tors and progeny included in the apple REFPOP. Some modern breeding material at this locus. For f lowering parental combinations showed almost exclusively the intensity, a trait positively correlated with apple yield allele increasing red skin color, other parental combi- (i.e. weight and number of fruits, Figure 2a), a new nations exhibited a lack of the allele. This uneven dis- locus was discovered on chromosome 14. The increased tribution of the alleles in the progeny group pointed to frequency of the allele contributing to higher f lowering different directions of selection for fruit skin color in the intensity in the progeny, its presence in all parental European breeding programs (Figure 5b). genotypes, and fixation in some parental combinations may be the result of breeding for high yield. The Performance of the main-effect univariate major locus found for acidity on chromosome 16 was genomic prediction models consistent with the Ma locus frequently detected in Main-effect univariate genomic prediction models were various germplasm [8, 11]. The total number of the applied to the global clonal values. The observed small high-acidity alleles for Ma and Ma3, which is another differences between genomic predictive abilities of regularly detected acidity locus, was shown to be various models (Figure 6a) were in accordance with higher in parents of a European breeding program previous model comparisons where distinctions among 12 | Horticulture Research, 2022, 9: uhac028 models were negligible [39, 59]. The extremes in pre- Role of genotype by environment interactions in dictive ability between traits were found with random multi-environment univariate genomic forest, which allowed for the overall highest predictive prediction ability among all compared models for red over color. The multi-environment univariate genomic prediction The explanation for the striking performance of random models either assumed that effects of markers were forest for red over color might be found in the results the same across environments (across-environment G- of our GWAS. This trait of oligogenic architecture was BLUP, called here G-BLUP.E) or additionally estimated associated with a few small-effect loci and one locus marker effects separately for each environment and thus of large effect explaining 61% of the red over color considered the G × E (marker by environment interaction phenotypic variance measured in the apple REFPOP. G-BLUP, called here G-BLUP.E.G × E) [42]. The average High correlations between many markers, i.e. linkage accuracy of the G-BLUP.E.G × E model across traits was disequilibrium, were found in the vicinity of the large- only slightly higher than the accuracy of the G-BLUP.E. effect locus (Supplementary Figure 8). Random forest In contrast, the G-BLUP.E.G × E model had substantially is known to assign higher importance to correlated greater predictive ability than the G-BLUP.E model when predictor variables (based on the amount of accu- applied in wheat [42]. In the latter study, a productivity racy decrease when a variable – here marker – is trait was measured under simulated conditions of mega- excluded during the decision tree construction) [60], environments and the effect of G × E explained up to which may have contributed to the particularly high ∼60% of the phenotypic variance [42]. Our work only predictive ability found for red over color with random focused on European environments and the largest forest. proportion of phenotypic variance assigned to G × E The predictive ability for red over color reached ∼0.4 was 24% for a productivity trait (flowering intensity). in several former prediction studies [22, 23, 29, 34] and Furthermore, the average proportion of G × E across traits was approximately doubled in our work, which demon- was approximately 12%, which may explain the mostly strated the potential of the current study design for negligible differences between the G-BLUP.E and G- accurate genomic predictions. For harvest date, the cur- BLUP.E.G × E models. Our results were in line with the low rently reported predictive ability of 0.78 was only slightly interaction of additive genetic effects with location of up higher than the accuracy of 0.75 obtained with the initial to ∼6% obtained for apple fruit quality traits measured apple REFPOP dataset measured during one year [36], but at two locations in New Zealand [33], the generally these accuracies showed a considerable improvement stable genetic effects found for apple fruit quality traits over other accuracies of approximately 0.5–0.6 reported assessed across the commercial production region of the elsewhere [22, 23, 29]. As shown before [36], these results State of Washington [61], and the limited G × E reported underline the suitability of apple REFPOP design for the for fruit maturity timing in sweet cherry across conti- application of genomic prediction. nents [62]. For approximately half of the tested traits, the Opposite to harvest date and red over color, the G-BLUP.E.G × E did not outperform G-BLUP.E. For these predictive ability of yellow color and russet cover was traits, the G-BLUP.E ignoring G × E may be sufficient to low, although the genotypic effects explained 45% and account for the environmental effects across European 47% of the phenotypic variance, respectively. The across- sites because it is computationally simpler and therefore environment clonal-mean heritability of russet cover less demanding in terms of computational resources. was high (0.97), while the heritability for yellow color Traits such as flowering intensity, soluble solids content, was slightly lower (0.81, Figure 7). Yellow color showed trunk increment or traits related to fruit size and russet a moderate phenotypic clonal values correlation of 0.55 showed an improved performance under G-BLUP.E.G × E with russet cover, suggesting that the phenotyping device when compared to G-BLUP.E. For traits where the might have classified some russet skin as yellow color. predictive ability was greater for the G-BLUP.E.G × E Symptoms of powdery mildew could have been misin- model compared to the G-BLUP.E model, the G × E should terpreted as russet skin. The decreased performance of be considered in marker effects estimated separately genomic prediction models might stem from inaccurate for each environment when making predictions across phenotyping methods, insufficient SNP density in the environments. The highest improvement of predictive associated regions, or other factors, all of which were ability with G-BLUP.E.G × E when compared to G-BLUP.E outside the scope of this work. was found for flowering intensity, the difference between All main-effect univariate genomic prediction models the models amounting to 0.07 (Figure 6b). This result as well as other genomic prediction models compared might be explained by the highest contribution of G × E to in this study depended on predictions of clonal values the phenotypic variance of flowering intensity among all obtained during the phenotypic data analysis. Further traits (Figure 7). A comparably high contribution of G × E extension of the models should consider if adding a term was also found for weight of fruits and number of fruits, to account for the permanent non-genetic effect of the though no improvement with G-BLUP.E.G × E model was tree over years into the mixed-effects models (Equation observed for these traits. When comparing the relative 1) improves accuracy of the clonal values predictions and contributions of variance components to the phenotypic so increases the genomic predictive ability. variance of flowering intensity, weight of fruits and Jung et al. | 13 number of fruits, the proportions of G × E were similar the multi-environment multivariate genomic prediction in the three compared traits, but the effects of genotype model (MTM.FA) showed the potential to perform equally and environment explained a higher proportion of the well for six (CV1) and three traits (CV2) and was able to variance for weight of fruits and number of fruits than outperform both models for two (CV1) and five traits for f lowering intensity. This may have contributed to the (CV2). Except for the noticeable increase in predictive surprisingly lower accuracy of the G-BLUP.E.G × E model ability for trunk increment under CV2 that could not be when compared with G-BLUP.E for weight of fruits and explained by our analyses, the performance of MTM.FA number of fruits, but additional investigations may be was similar to G-BLUP.E and G-BLUP.E.G × E, which needed to clarify this result in the future. establishes the multivariate model as a useful tool for The G-BLUP.E.G × E model assumes positive correla- multi-environment genomic prediction in apple. tions between environments and is therefore mostly suitable for the joint analysis of correlated environments Two approaches to genomic prediction addressed [42, 63]. As shown by Lopez-Cruz et al. [42] and in our with cross-validation scenarios study, this assumption of G-BLUP.E.G × E resulted in the The cross-validation scenarios CV1 and CV2 were applied best model performance for traits showing high positive with multi-environment genomic prediction models to correlations of the adjusted phenotypic values of each test two genomic prediction approaches typically faced genotype between environments (here harvest date and in breeding. The CV1 imitated evaluation of breeding red over color) and the worst performance for traits material that was yet untested in field trials. The CV2 exhibiting low correlations of the adjusted phenotypic was implemented to simulate incomplete field trials values of each genotype between environments (here where breeding material was evaluated in some but flowering intensity and trunk increment, Figure 6b, not all target environments. More specifically, the CV2 Supplementary Table 2, Supplementary Figure 3). For investigated a situation where the breeding material has flowering intensity and trunk increment, multivariate been evaluated at one location (the breeding site, in this prediction of the environments or prediction with a case Switzerland) and the material’s potential over other different G × E model not assuming positive correlations European sites was predicted without its assessment in between environments might be more appropriate than a multi-environment trial, which may increase selection the currently applied approach [42, 64]. efficiency at latter stages of evaluation. As CV2 provided more phenotypic information to the models than CV1, a Multivariate models as a useful element in the higher genomic predictive ability was found under CV2 genomic prediction toolbox when compared with CV1, which was anticipated [33, Multivariate (also called multi-trait) models were shown 42]. The CV2 was tested by calibrating the model with to improve predictive ability for traits that are costly Swiss observations only. The application of CV2 could to phenotype when a correlated trait less expensive to be extended to other apple REFPOP locations to provide phenotype was available [45]. In our study, when the useful information for the breeding programs located at predictive ability of the main-effect multivariate model these sites. The choice of cross-validation scenario did MTM.UN was compared with the baseline model G-BLUP, not affect the general ranking of the average genomic several combinations of related and unrelated traits predictive abilities estimated for the evaluated traits. led to increased accuracy. For the related traits with a high genomic breeding values correlation (Figure 2a), Implications for apple breeding prediction of traits measured at one location were Phenotypic variance decomposition into genetic, envi- often improved when a related trait measured across ronmental, G × E and residual effects was compared with different locations was included. This was the case for the results of GWAS and genomic prediction as well as the combination of f loral emergence with full f lowering heritability estimates. The comprehensive comparison and end of f lowering and for single fruit weight combined indicated three classes of traits with contrasting genetic with fruit diameter, fruit length, maximum fruit size architecture and prediction performance. Characteristics and fruit volume. Inclusion of soluble solids content of these trait classes and proposals for their efficient in MTM.UN resulted in increased predictive ability for prediction strategies are described in the following para- russet cover, although the traits showed only a moderate graphs. genomic breeding values correlation and no obvious The first class included harvest date and red over explanation for this result could be found. Our study color that showed a few loci of large effect and some supports the potential of multivariate models to borrow additional loci of low effect, the highest predictive abil- information that correlated traits provide about one ities, and the highest across-environment clonal-mean another and identified trait combinations that can be heritability among all traits. Both traits showed a very successful under the multivariate setup. high proportion of the genotypic effect explaining ∼75% In place of the correlated traits, environments of a of the phenotypic variance. For harvest date and red single trait can be implemented in a multivariate model over color, the marker with the largest effect explained [46]. Compared to the multi-environment univariate 52% and 59% of the phenotypic variance and all marker genomic prediction models G-BLUP.E and G-BLUP.E.G × E, effects in genomic prediction captured together 88% and 14 | Horticulture Research, 2022, 9: uhac028 85% of the phenotypic variance (i.e. genomic heritability measured at one location, the combinations of locations of 0.88 and 0.85), respectively. Selection for these traits and years were effectively equal to years only. Although exhibiting a strong genetic effect of one locus could measurements from additional locations would improve be done using marker-assisted selection, although only estimation of the effect of G × E in traits with unavail- a part of the variance would be explained by a single able multi-location data, our results were generally able marker. Better results can be achieved using genomic to support that multi-environment genomic prediction prediction, as this was able to explain a substantially models can be efficient when applying genomic selection larger amount of the phenotypic variance. Other traits to various traits in apple. Decomposition of the effects of such as fruit firmness, titratable acidity, end of f lowering environment and G × E into parts associated with loca- or traits related to fruit size and water core were grouped tions, years and their interaction could indicate whether in the same cluster as harvest date and red over color defining breeding zones would be useful for apple, but it (Figure 7). These traits showed a strong genotypic effect was out of the scope of this study. and a comparably low effect of environment and G × E, The decision to apply either marker-assisted or suggesting that selection for the traits would be efficient genomic selection can be based on genetic architec- when performed using main-effect genomic prediction ture of traits of interest and resources available in models rather than multi-environment prediction. a breeding program. For breeding of yet genetically The second class of traits was represented by flo- unexplored traits, variance decomposition of historical ral emergence and trunk diameter displaying a high phenotypic data prior to genomic analyses may help proportion of the environmental effect (∼70%) and a describe trait architecture, assign traits to one of the similar ratio of variance explained by genotypic effects three classes described in the previous paragraphs, compared to variance explained by G × E effects (∼2.5). and finally determine the most appropriate method of The genomic predictive ability did not considerably devi- genomics-assisted breeding. From all traits explored in ate from the average accuracy over all traits. Several this study, the marker-trait associations with large and marker associations with these traits were identified stable effects across environments found for harvest using location-specific GWAS. However, in the across- date, flowering intensity, green color, red over color, location GWAS, only one association explaining a very titratable acidity, fruit firmness and trunk increment small part of phenotypic variance (f loral emergence) or could be implemented into DNA tests for marker-assisted no association (trunk diameter) were discovered. Conse- selection. These tests would allow for a reduction of labor quently, such traits predominantly driven by the effect costs in a breeding program when removing inferior of environment can be successfully selected based on seedlings without phenotyping [7]. Although generally genomic prediction, but the lack of associations stable requiring more statistical competences than marker- across environments limits the applicability of marker- assisted selection, genomic selection can make use assisted selection to this class of traits. of both large- and small-effect associations between In the third class, the productivity traits (f lowering markers and traits when accommodating thousands of intensity, weight of fruits and number of fruits) showed marker effects in a single genomic prediction model. For the largest proportion of variance explained by G × E all studied traits, our results showed that marker effects (∼20%), with similar amounts of variance explained by estimated in genomic prediction were able to capture genotypic effects for weight of fruits and number of a larger proportion of the phenotypic variance than fruits, but half as much variance explained by genotypic individual markers associated with the traits. Therefore, effects for f lowering intensity (Figure 7). As a conse- genomic selection should become the preferred method quence, only f lowering intensity showed higher predic- of genomics-assisted breeding for all quantitative traits tive ability with G-BLUP.E.G × E than G-BLUP.E model. As explored in this study to ultimately increase their shown above, the G × E should be considered when mak- breeding efficiency and genetic gain. ing predictions across environments for traits respond- ing positively to the G-BLUP.E.G × E model, but G-BLUP.E may be sufficient for other traits to account for the Conclusion environmental effects. To our knowledge, this is the first This study laid the groundwork for marker-assisted and report of genomic prediction for apple yield components genomic selection across European environments for and our results can aid the establishment of produc- 30 quantitative apple traits. The apple REFPOP experi- tivity predictions in apple breeding. Other traits falling mental design facilitated identification of a multitude within the same cluster as the productivity traits, namely of novel and known marker-trait associations. Our full f lowering, ground color, yellow color, soluble solids multi-environment trial provided accurate genomics- content, trunk increment, and traits related to bitter pit estimated breeding values for apple genotypes under and russet, showed a pronounced effect of environment various environmental conditions. Limited G × E detected and/or G × E (Figure 7). For the majority of these and in this work suggested consistent performance of other traits in our study, the effect of G × E was estimated genotypes across different European environments for based on environments that represented combinations most studied traits. Utilizing our dataset, genomic of locations and years. However, for the part of traits selection of traits related to yield may lead to higher Jung et al. | 15 productivity and increased genetic gain in the future iQS4 v.1.0, the instrument Pimprenelle (Setop, France)). [37]. The genomic prediction models developed here Two phenology traits measured in 2018, i.e. floral emer- can be readily used for selecting germplasm in breeding gence and harvest date, were previously analyzed by Jung programs, thus providing breeders with tools increasing et al. [36]. selection efficiency. Application of our study design to other horticultural crops such as peach [65] can promote Phenotypic data analyses broader use of genomics-assisted breeding in the future. Spatial heterogeneity was modeled separately for each trait and environment (combined factor of location and year) using the spatial analysis of field trials with splines Methods (SpATS) to correct for the replicate effects and differences Plant material due to soil characteristics [68]. Phenotypic values of traits Plant material in this study was comprised of the apple adjusted for spatial heterogeneity within each environ- REFPOP, which was designed and established by the ment were estimated at the level of trees (adjusted phe- collaborators of the FruitBreedomics project [66] as notypic values of each tree) and genotypes (adjusted described by Jung et al. [36]. The apple REFPOP consisted phenotypic values of each genotype) as described before of 534 genotypes from two groups of diploid germplasm. [36]. The accession group consisted of 269 accessions of Further analyses were performed to estimate trait her- European and non-European origin representing the itability and remove the effects of location and year diversity in cultivated apple. The progeny group of 265 from the collected phenotypes. The general statistical genotypes stemmed from 27 parental combinations model for the following phenotypic data analyses fitted produced in the current European breeding programs. via restricted maximum likelihood (R package lme4 [69]) In 2016, the apple REFPOP was planted in six locations was: representing several biogeographical regions in Europe, in (i) Rillaar, Belgium, (ii) Angers, France, (iii) Laimburg, y = Xβ + Zb + ε (Equation 1) Italy, (iv) Skierniewice, Poland, (v) Lleida, Spain and (vi) where y was a vector of the response variable, X the Wädenswil, Switzerland (one location per country). Every design matrix for the fixed effects, β the vector of fixed genotype was replicated at least twice per location. effects, Z the design matrix for the random effects, b All plants included in this study were treated with the vector of random effects and ε the vector of random   agricultural practice common to each location. Calcium errors. The b was a q × 1 vector assuming b ∼ N 0,  spraying was avoided due to its inf luence on bitter pit. where  was a variance–covariance matrix of the ran- Flowers were not thinned, but the fruits were hand- thinned after the June fruit drop and up to two apples   for the N × 1 vector of ran- dom effects. The assumptions dom errors were ε ∼ N 0, Iσε2 with N × N identity matrix per fruit cluster were retained. I and the variance σε2 , the N being the number of trees. Genotyping To assess the reliability of environment-specific data, The plant material was accompanied by a high-density a random-effects model was first fitted separately for genome-wide SNP marker dataset, which was produced each trait and environment to estimate an environment- as reported by Jung et al. [36]. Brief ly, SNPs from two specific clonal mean heritability. Applying the Equation overlapping SNP arrays of different resolution, (i) the 1, the response y was a vector of the adjusted phenotypic Illumina Infinium® 20 K SNP genotyping array [20] and (ii) values of each tree. On the place of X, a vector of ones the Affymetrix Axiom® Apple 480 K SNP genotyping array was used to model the intercept β. The design matrix [21], were curated and then joined applying imputation Z reflected that the genotype was the grouping factor with Beagle 4.0 [67] using the recently inferred pedigrees defining the random effects. The environment-specific [4]. Non-polymorphic markers were removed to obtain clonal mean heritability was calculated from the vari- a set of 303 148 biallelic SNPs. Positions of SNPs were ance components of the random-effects model as: based on the apple reference genome obtained from the doubled haploid GDDH13 (v1.1) [16]. σg2 H2 = (Equation 2) σp2 Phenotyping Thirty phenotypic traits related to phenology, produc- where the phenotypic variance σp2 = σg2 + σε2 /nr was tivity, fruit size, outer fruit, inner fruit, and vigor were obtained from the genotypic variance σg2 , error variance evaluated at up to six locations of the apple REFPOP σε2 and the mean number of genotype replications nr . during up to three seasons (2018–2020). Trunk diame- The environment-specific clonal mean heritability was ter was measured in 2017 in some locations, enabling used to eliminate location-year-trait combinations with for a trunk increment calculation for 2018. The traits a heritability value below 0.1. were recorded as described in the Supplementary Meth- For the remaining location-year combinations, a ods, the measurements being performed either visu- single-trait mixed-effects model following the Equation ally or using automatic devices (sorting machine Greefa 1 was fitted to the vector of the adjusted phenotypic 16 | Horticulture Research, 2022, 9: uhac028 values of each tree as response (y) to estimate the across- clonal values were scaled and centered and their missing environment clonal mean heritability. The effects of values for each trait were replaced with the mean of environments, i.e. combination of location and years, the global clonal values of the same trait. A multivariate were used as fixed effects and the effects of genotypes normal distribution was assumed for the ellipses in the and genotype by environment interactions as random biplot. effects. Estimated variances of the model components were used to evaluate the across-environment clonal Genome-wide association studies mean heritability calculated using the Equation 2 with As one of the currently most powerful procedures the phenotypic variance estimated as: for identification of loci associated with complex traits in terms of computational speed and statistical 2 σge σε2 power [71, 72], the Bayesian-information and linkage- σp2 = σg2 + + (Equation 3) ne ne nr disequilibrium iteratively nested keyway (BLINK) [72] was chosen to perform the GWAS. BLINK applies two fixed 2 where σge was the genotype by environment interac- effect models and one filtering process for the choice of tion variance and ne represented the number of envi- associated markers that are not in linkage disequilibrium ronments. As the minimum number of genotype repli- with each other as covariates. The process is repeated cates at each location was only two, the variation among until all markers are tested and the choice of associated genotype replicates that stemmed from differences inde- markers is optimized using Bayesian information criteria pendent of the environmental inf luence may not be [72]. BLINK implemented in the R package GAPIT 3.0 [73] accurately captured. This may contribute to uncertainty was applied using an n × m matrix for a population of in the estimation of phenotypic variance components. size n = 534 genotypes (i.e. accessions and progeny) The mixed-effects model (Equation 1) did not account with m = 303 148 markers, with across-location clonal for the permanent non-genetic effect of the tree over values (across-location GWAS) or location-specific clonal years, which could lead to inf lated values of clonal-mean values (location-specific GWAS) as the response. BLINK heritability. was used with two principal components and the minor To predict across-location clonal values (and location- allele frequency threshold was set to 0.05. Marker- specific clonal values when only single location data was trait associations were identified as significant for p- available), an additional mixed-effects model following values falling below a Bonferroni-corrected significance the Equation 1 was fitted to the adjusted phenotypic threshold α ∗ = α/m with α = 0.05 (−log10 (p) > 6.74). values of each tree (y) using the effects of location, year The proportion of phenotypic variance explained by and their interaction as fixed effects and the effects each significantly associated SNP was assessed with a of genotypes as random effects. Due to the skewness coefficient of determination (R2 ). The R2 was estimated of their distributions, y-values of the traits weight of from a linear regression model, which was fitted with fruits, number of fruits and trunk diameter were log- a vector of SNP marker values (coded as 1, 2, 3) as transformed. BLUPs (b̂) extracted from the model were predictor and either the across-location clonal values or further denoted as across-location clonal values. To location-specific clonal values as response. GWAS based estimate the location-specific clonal values, a model on the across-location clonal values with the addition of according to the Equation 1 was fitted with a subset of location-specific clonal values, in cases where traits were the adjusted phenotypic values of each tree from single measured at a single location only, was further denoted locations (y) using the effects of years as fixed effects and as the global GWAS. The position of the last SNP on a the effects of genotypes as random effects. The across- chromosome was used to estimate chromosome length, location clonal values with the addition of location- which was used to divide each chromosome into three specific clonal values for traits measured at a single equal segments, i.e. top, center and bottom. The marker- location were further denoted as the global clonal values. trait associations were assigned to these chromosome It should be noted that due to the different estimation segments based on their positions to allow for a of the two elements of global clonal values (the use of subsequent comparison with published associations. fixed effects of location, year and their interaction to Previous reports on QTL mapping and GWAS in apple estimate the across-location clonal values versus the were reviewed to perform an extensive comparison with fixed effect of year used to estimate the location-specific our GWAS results (Supplementary Table 4). Published clonal values), the elements of the global clonal values results for traits measured similarly to the traits studied are not strictly comparable. in the present work were considered, with the traits The global clonal values were used to assess pheno- being assembled into trait groups: harvest time (harvest typic clonal values correlation as the Pearson correlation date and similar), flowering time (floral emergence, full between pairs of traits. The correlation between pairs of flowering, end of flowering and similar), productivity environments within traits was calculated as the Pear- (flowering intensity, weight of fruits, number of fruits and son correlation between the adjusted phenotypic val- similar), fruit size (single fruit weight, fruit diameter, fruit ues of each genotype within environments. To estimate length, maximum fruit size, fruit volume and similar), the principal component analysis biplot [70], the global ground color (ground color, yellow color and similar), Jung et al. | 17 over color (red over color, green color and similar), bitter variance that was explained by the SNPs (phenotypic pit (bitter pit frequency, bitter pit grade and similar), variance associated with the random effects of SNPs) and russet (russet cover, russet frequency overall, at stalk, that remained unexplained by the SNPs (phenotypic vari- on cheek and in the eye and similar), acidity (titratable ance associated with the random effects of genotypes). acidity and similar), sugar (soluble solids content and Due to the low number of levels of the random effect similar), firmness (fruit firmness and similar), water terms of SNPs, uncertainty in the estimation of their core (water core frequency, water core grade and variance may have been introduced. similar) and trunk (trunk diameter, trunk increment A centered and scaled table (mean 0, standard devia- and similar). The positions of published associations tion of 1) of trait by proportion of phenotype variances within respective chromosomes were visually assigned explained by different sources (genotypic, environmen- to the three chromosome segments, i.e. top, center tal, genotype by environment interaction, and residual and bottom. The total number of markers used was effects) was constructed and hierarchical clustering fol- recorded (Supplementary Table 4). Where the number of lowing Ward [74] was applied to the distance matrix of overlapping markers between the maternal and paternal this table. The number of clusters was estimated from a linkage maps was not provided in a publication, the dendrogram, which was cut where the distance between marker numbers for both maps were summed. splits was the largest. In the global GWAS results, the allele frequency was studied over generations. The ancestors of genotypes Genomic prediction were identified making use of the apple pedigrees of Main-effect genomic prediction Muranty et al. [4]. For all significant marker-trait asso- Four univariate and one multivariate main-effect genomic ciations from the global GWAS, frequency of the allele prediction models were used to evaluate predictive associated with increased phenotypic value was esti- ability for 30 phenotypic traits. These models differed mated for the progeny group and for its five ancestor in the way of estimating the marker effects. For all generations. For comparison with the allele frequency models, the n × m matrix for a population of size n = 534 in the progeny group, allele frequency was estimated genotypes with m = 303 148 markers was centered by separately for the 30 accessions of the progeny group that subtracting the column means from their corresponding were included in the apple REFPOP. Additionally, mean columns, then scaled by dividing the columns by their allele frequencies and standard errors were estimated for standard deviations, and further denoted as the additive 10-times repeated resampling of 30 genotypes from the genomic matrix M. The models were fitted with the progeny group. global clonal values predicted during the phenotypic For major significant marker-trait associations with data analysis. R2 > 0.1 reported in the global GWAS, linkage disequilib- The first univariate model examined was regression rium was estimated as squared Pearson’s correlations in with random forest (RF) [75], which is a non-parametric a window of 3000 markers surrounding each of the asso- model that may be able to capture non-additive effects ciations. A smaller window size was used for associations in addition to additive effects. For RF, the response y was located towards the end of a chromosome. defined as a n × 1 vector of the global clonal values and A mixed-effects model was used for every trait to esti- the columns of the matrix M were used as predictors. The mate proportions of phenotypic variance explained by number of decision trees in the RF was set to 500 and the the model components as described in Jung et al. [36]. The number of variables randomly sampled as candidates at mixed-effects model following the Equation 1 was fitted each split was (rounded down) mtry = m/3. to the vector of the adjusted phenotypic values of each For the remaining three univariate and one multivari- tree as response (y) using the effects of environments ate main-effect genomic prediction models, the general as fixed effects and the random effects of each SNP random effects model was defined as: significantly associated with the trait (a factor of the respective SNP values), the remaining random effects of y = 1μ + u + ε (Equation 4) genotypes and genotype by environment interactions. In cases where traits with no marker-trait associations were where y was a response vector of the global clonal values, found in the global GWAS, the additional random effects μ was an intercept, u was a term used to specify random of significantly associated SNPs were omitted from the effects and ε was a vector of residuals. model. The proportion of phenotypic variance associated The second univariate model BayesCπ is a parametric with the fixed effect of environment was estimated as the model, which estimates prior probability π that a genetic variance of the vector of values predicted from the model marker has zero effect [76]. Following the Equation 4, the fit when all random effects were set to zero. The propor- response y was a n × 1 vector of the global clonal values m tions of phenotypic variance explained by the random for one trait, the term u = k=1 zk ak with zk being an effects of significantly associated SNPs and genotypes n × 1 vector of the number of copies of one allele at the were summed to obtain the genotypic variance. The SNPs marker k and ak being the additive effect of the marker significantly associated with traits and the genotypes k. The prior for ak depended on the variance σa2 and the k were treated as random effects to approximate genotypic prior probability π that a marker k had zero effect, the 18 | Horticulture Research, 2022, 9: uhac028 priors of all marker effects having a common variance genomic breeding values correlations larger than 0.3 σa2 = σa2 . The π parameter was treated as an unknown were retained for the MTM.UN analysis. In case a trait k with uniform [0, 1] prior. The random vector of residual was included in more than one pair of traits, the result effects followed a normal distribution ε ∼ N(0, Iσε2 ) with for the pair with the highest average predictive ability for n × n identity matrix I and the variance σε2 . this trait was reported. The third examined univariate model, the semi- With all models, a five-fold cross-validation repeated parametric Bayesian reproducing kernel Hilbert spaces five times was performed, generating 25 estimates of pre- regression (RKHS), is able to capture additive as well as dictive ability. The folds were chosen randomly without non-additive effects and was implemented here using replacement to mask phenotypes of 20% of the geno- a multi-kernel approach [77]. The multi-kernel RKHS types in each run. Predictive ability was estimated as model was fitted according to the Equation 4 with the a Pearson correlation coefficient between phenotypes response y being a n × 1 vector of the global clonal values of the masked genotypes (global clonal values) and the for one trait, and L = 3 n-dimensional vectors of the predicted values for the same genotypes (genomic breed- L random effects u. The vectors u = l=1 ul followed a ing values, i.e. the average predictions from all individ- distribution u ∼ N(0, Kl σul ), with Kl being the reproducing 2 ual regression trees for RF and the estimated posterior kernel evaluated at the lth value of the bandwidth means of random effects excluding the residuals for parameter h = {h1 , . . . , hL } = {0.1, 0.5, 2.5} and the BayesCπ , RKHS, G-BLUP and MTM.UN). variance σul2 . For each random effect, the kernel matrix BayesCπ, RKHS, G-BLUP and MTM.UN were applied K = {K(xi , xi )} was an n×n matrix K(xi , xi ) = exp{−h×Dii }, with 12 000 iterations of the Gibbs sampler, a thinning  m (x −x )2  where D = Dii = k=1 ik ik was the average of 5, and a burn-in of 2000 discarded samples. The RF m model was implemented in the R package ranger [79], squared-Euclidean distance matrix between genotypes, the models BayesCπ , RKHS and G-BLUP in the R package and xik the element on line i (genotype i) and column k BGLR80 and the MTM.UN model in the R package MTM (kth marker) of the matrix M. The residual effect assumed (http://quantgen.github.io/MTM/vignette.html). ε ∼ N(0, Iσε2 ). The fourth univariate model genomic-BLUP (G-BLUP) Multi-environment genomic prediction was fitted using a semi-parametric RKHS algorithm. To facilitate efficient incorporation of a large number of To explore the effects of genotypes, environments and markers, the additive genomic relationship matrix tradi- their interaction in genomic prediction, the predictive tionally used for genomic prediction was replaced in G- ability for 30 traits was estimated using an across- BLUP by a genomic relationship matrix [78]. The genomic environment and a marker by environment interaction relationship matrix G was computed as G = MM /m and univariate genomic prediction algorithms that assumed used to fit the G-BLUP model following the Equation 4. constant or changing random marker effects across The response y was a n × 1 vector of the global clonal environments, respectively. The random effects model values for one trait. The n-dimensional vector of random for the examined univariate multi-environment models effects followed u ∼ N(0, Gσu2 ) with variance σu2 and the that were reported by Lopez-Cruz et al. [42] was: model residuals assuming ε ∼ N(0, Iσε2 ). The fifth model applied was a multivariate model with y = 1μ + u + ε (Equation 5) an unstructured covariance matrix of the random effect (here abbreviated as MTM.UN). The model was fitted for chosen pairs of traits using the Bayesian multivariate where the response y was a vector of the adjusted phe- Gaussian model environment MTM (http://quantgen.gi notypic values of each genotype of length n × r (with r thub.io/MTM/vignette.html). The MTM.UN followed the equal to the number of environments, the environments Equation 4 with the response vector y, which was a vector being represented as the combined factor of location of the global clonal values for t traits with the length and year), μ was the vector with an intercept for each n · t and t = 2 being the number of traits used in the environment, u represented the vector of random effects model. The vector of random effects with dimension of length n × r and ε was a vector of residuals. Of the n · t followed u ∼ N(0, U ⊗ G) where U was an unstruc- two univariate models, the across-environment G-BLUP tured (within-genotype) covariance matrix of the random model (G-BLUP.E) assumed that marker effects were con- effects with dimension t × t. Model residuals assumed stant across environments with u ∼ N(0, G0 σu2 ) where ε ∼ N(0, R ⊗ I) with R being a t × t (within-genotype) G0 = J ⊗ G, the J being an r × r matrix of ones. The unstructured covariance matrix of the residual effect. To model residuals assumed ε ∼ N(0, Iσε2 ). Additionally to choose the pairs of traits for MTM.UN, the G-BLUP model the constant effects of markers across environments as was applied as described above using all genotypes to assumed in the previous model, the marker by environ- estimate genomic breeding values (estimated posterior ment interaction G-BLUP model (G-BLUP.E.G × E) allowed means of random effects excluding the residuals), which the marker effects to change across environments. The were then used to obtain pairwise genomic breeding random effects were defined as u = u0 + u1 where u0 ∼ values correlations between traits. The pairs with the 2 N(0, G0 σu0 ) represented the random effects common to Jung et al. | 19 all environments and u1 ∼ N(0, G1 ) the random devia- All three multi-environment genomic prediction mod- tions of the effects for specific environments with: els were applied with 12 000 iterations of the Gibbs sam- pler, a thinning of 5 and a burn-in of 2000 discarded ⎡ ⎤ 2 σu1 G 0 0 samples. The models G-BLUP.E and G-BLUP.E.G × E were G1 = ⎣ 0 2 σu2 G 0 ⎦ implemented in the R package BGLR [80], the model 2 0 0 σu3 G MTM.FA in the R package MTM (http://quantgen.githu b.io/MTM/vignette.html). assuming r = 3 here for easier notation. The model residuals assumed ε ∼ N(0, Iσε2 ). Genomic heritability To test multivariate multi-environment genomic pre- The BayesCπ model was applied for each trait as diction for the traits under study, a multivariate multi- described before but trained with a full set of the global environment factor-analytic model (here abbreviated as clonal values as response. The genomic heritability MTM.FA) was fitted to the genomic and phenotypic data h2 = Vg /(Vg + Ve ) was estimated as the proportion of using the Bayesian multivariate Gaussian model envi- phenotypic variance explained by the markers, where Vg ronment MTM (http://quantgen.github.io/MTM/vignette. and Ve represented the amount of phenotypic variance html). The traits measured at only one location dur- explained and unexplained by the markers, respectively ing two seasons (full f lowering, end of f lowering, fruit [81, 82]. The genomic heritability was calculated from volume, water core frequency and water core grade) the marker effects saved in each iteration of the Gibbs were not modeled using MTM.FA because the analy- sampler and averaged over iterations to obtain the mean sis required at least three environments. Following the genomic heritability per trait. Equation 5, the vector of random effects assumed u ∼ N(0, C ⊗ G) where C was an r × r genetic-by-environment covariance matrix. For the factor analysis, the C = BB + Acknowledgements where B was a matrix of loadings (regressions of the The authors thank the field technicians and staff, original random effects into common factors) and  especially Sylvain Hanteville, at INRAe IRHS and Exper- was a diagonal matrix whose entries gave the variances imental Unit (UE Horti), Angers, France, and technical of environment-specific factors. The loadings were esti- staff at other apple REFPOP sites for the maintenance mated for all environments and the variance of the Gaus- of the orchards and phenotypic data collection. We sian prior assigned to the unknown loadings was set to thank Dr. Graham Dow for English language editing. 100. The model residuals assumed ε ∼ N(0, R ⊗ I) with R Phenotypic data collection was partially supported by being an r × r (within-genotype) unstructured covariance the Horizon 2020 Framework Program of the European matrix of the residual effect and I the n-dimensional Union under grant agreement No 817970 (project INVITE: identity matrix. “Innovations in plant variety testing in Europe to foster The folds of a five-fold cross-validation were chosen the introduction of new varieties better adapted to randomly without replacement. The cross-validation varying biotic and abiotic conditions and to more was repeated under two scenarios. In the first cross- sustainable crop management practices”). This work validation scenario (CV1), the phenotypes of 20% of was partially supported by the project RIS3CAT (COTPA- the genotypes were masked across all environments. FRUIT3CAT) financed by the European Regional Devel- For the second cross-validation scenario (CV2), the opment Fund through the FEDER frame of Catalonia phenotypes of 20% of the genotypes were masked across 2014-2020 and by the CERCA Program from Generalitat all environments except for three Swiss environments, de Catalunya. We acknowledge financial support from i.e. phenotypes of all genotypes from the environments the Spanish Ministry of Economy and Competitiveness “CHE.2018”, “CHE.2019” and “CHE.2020” were used for through the “Severo Ochoa Programme for Centres model training. Ten traits were measured in only one of Excellence in R&D” 2016-2019 (SEV-20150533) and location and therefore excluded from CV2 (i.e. full 2020-2023 (CEX2019-000902-S). C.D. was supported by flowering, end of f lowering, fruit diameter, fruit length, “DON CARLOS ANTONIO LOPEZ” Abroad Postgraduate maximum fruit size, fruit volume, yellow color, green Scholarship Program, BECAL-Paraguay. We dedicate this color, water core frequency and water core grade). paper to Prof. Edward Zurawicz of the National Institute Predictive ability was estimated as a Pearson correlation of Horticultural Research in Skierniewice, Poland who coefficient between the phenotypes of the masked co-promoted this study, but sadly recently passed away. genotypes (adjusted phenotypic values of each genotype in each environment) and the predicted values for these genotypes in each environment (estimated posterior Author contributions means of random effects excluding the residuals). The M.J.A., W.G., F.L., H.M. and A.P. conceived the research correlations were estimated for each predicted environ- plans; B.S., H.M. and A.P. supervised the project; M.J., B.K., ment separately. The models following the Equation 5 M. Roth, M.J.A., A.A., M.A., M.L., N.S., M. Rymenants, F.D., did not account for the permanent non-genetic effect of C.D., C.F.F. and A.K. contributed to data collection; M.J. the tree over years. carried out the data analysis and wrote the article in 20 | Horticulture Research, 2022, 9: uhac028 consultation with B.K., M. Roth, H.M. and A.P.; M.J.A., A.A., breeding insights from a multi-family study. Tree Genet Genomes. W.G., N.S., M. Rymenants, F.L. and B.S. provided critical 2019;15:18. https://doi.org/10.1007/s11295-019-1324-y. feedback to the article. 12. Baumgartner IO, Kellerhals M, Costa F et al. Development of SNP-based assays for disease resistance and fruit qual- ity traits in apple (Malus × domestica Borkh.) and valida- Data availability tion in breeding pilot studies. Tree Genet Genomes. 2016;12:35. https://doi.org/10.1007/s11295-016-0994-y. All SNP genotypic data generated with the 480 K array 13. Iezzoni AF, McFerson J, Luby J et al. RosBREED: bridging the used in this study have been deposited in the INRAe chasm between discovery and application to enable DNA- dataset archive (https://data.inrae.fr/) at https://doi.o informed breeding in rosaceous crops. Hortic Res. 2020;7:177. rg/10.15454/IOPGYF. All SNP genotypic data generated https://doi.org/10.1038/s41438-020-00398-7. using the 20 K array used in this study have been 14. Chagné D, Vanderzande S, Kirk C et al. Validation of SNP markers deposited in the INRAe dataset archive at https:// for fruit quality and disease resistance loci in apple (Malus × doi.org/10.15454/1ERHGX. The raw phenotypic data are domestica Borkh.) using the OpenArray® platform. Hortic Res. available in the INRAe dataset archive at https://doi.o 2019;6:30. https://doi.org/10.1038/s41438-018-0114-2. rg/10.15454/VARJYJ. 15. Velasco R, Zharkikh A, Affourtit J et al. The genome of the domes- ticated apple (Malus × domestica Borkh.). Nat Genet. 2010;42: 833–9. https://doi.org/10.1038/ng.654. Competing interests 16. Daccord N, Celton J-M, Linsmith G et al. High-quality de novo assembly of the apple genome and methylome dynam- The authors declare no competing interests. ics of early fruit development. Nat Genet. 2017;49:1099–106. https://doi.org/10.1038/ng.3886. Supplementary data 17. Zhang L, Hu J, Han X et al. A high-quality apple genome assembly reveals the association of a retrotransposon Supplementary data is available at Horticulture Research and red fruit colour. Nat Commun. 2019;10:1494. Journal online. https://doi.org/10.1038/s41467-019-09518-x. 18. Sun X, Jiao C, Schwaninger H et al. Phased diploid genome assemblies and pan-genomes provide insights into the genetic References history of apple domestication. Nat Genet. 2020;52:1423–32. 1. FAOSTAT (Food and Agriculture Organization of the United https://doi.org/10.1038/s41588-020-00723-9. Nations, 2019). 19. Broggini GAL, Schlathölter I, Russo G et al. Chromosome- 2. Cornille A, Giraud T, Smulders MJM et al. The domestication scale de novo diploid assembly of the apple cultivar and evolutionary ecology of apples. Trends Genet. 2014;30:57–65. ‘Gala Galaxy’. bioRxiv. 2020); 2020.2004.2025.058891. https://doi.org/10.1016/j.tig.2013.10.002. https://doi.org/10.1101/2020.04.25.058891. 3. Way RD, Aldwinckle HS, Lamb RC et al. Apples (Malus). Acta Hortic. 20. Bianco L, Cestaro A, Sargent DJ et al. Development and validation 1991;3–46. https://doi.org/10.17660/ActaHortic.1991.290.1. of a 20K single nucleotide polymorphism (SNP) whole genome 4. Muranty H, Denancé C, Feugey L et al. Using whole- genotyping array for apple (Malus × domestica Borkh). PLoS One. genome SNP data to reconstruct a large multi-generation 2014;9:e110377. https://doi.org/10.1371/journal.pone.0110377. pedigree in apple germplasm. BMC Plant Biol. 2020;20:2. 21. Bianco L, Cestaro A, Linsmith G et al. Development and vali- https://doi.org/10.1186/s12870-019-2171-6. dation of the Axiom® Apple480K SNP genotyping array. Plant J. 5. Migicovsky Z, Gardner KM, Richards C et al. Genomic 2016;86:62–74. https://doi.org/10.1111/tpj.13145. consequences of apple improvement. Hortic Res. 2021;8:9. 22. Migicovsky Z, Gardner KM, Money D et al. Genome to phenome https://doi.org/10.1038/s41438-020-00441-7. mapping in apple using historical data. The Plant Genome. 2016;9. 6. Urrestarazu J, Denancé C, Ravon E et al. Analysis of https://doi.org/10.3835/plantgenome2015.11.0113. the genetic diversity and structure across a wide range 23. McClure KA, Gardner KM, Douglas GM et al. A genome-wide of germplasm reveals prominent gene flow in apple association study of apple quality and scab resistance. The at the European level. BMC Plant Biol. 2016;16:130. Plant Genome. 2018;11:170075. https://doi.org/10.3835/plantgeno https://doi.org/10.1186/s12870-016-0818-0. me2017.08.0075. 7. Wannemuehler SD, Luby JJ, Yue C et al. A cost–benefit analy- 24. Hirschhorn JN, Daly MJ. Genome-wide association studies for sis of DNA informed apple breeding. HortScience. 2019;54:1998. common diseases and complex traits. Nat Rev Genet. 2005;6: https://doi.org/10.21273/hortsci14173-19. 95–108. https://doi.org/10.1038/nrg1521. 8. Maliepaard C, Alston FH, van Arkel G et al. Aligning male 25. Kumar S, Garrick DJ, Bink MCAM et al. Novel and female linkage maps of apple (Malus pumila Mill.) genomic approaches unravel genetic architecture of using multi-allelic markers. Theor Appl Genet. 1998;97:60–73. complex traits in apple. BMC Genomics. 2013;14:393. https://doi.org/10.1007/s001220050867. https://doi.org/10.1186/1471-2164-14-393. 9. Kenis K, Keulemans J, Davey MW. Identification and stability of 26. Urrestarazu J, Muranty H, Denancé C et al. Genome-wide QTLs for fruit quality traits in apple. Tree Genet Genomes. 2008;4: association mapping of flowering and ripening periods in apple. 647–61. https://doi.org/10.1007/s11295-008-0140-6. Front Plant Sci. 2017;8:1923. https://doi.org/10.3389/ 10. Jänsch M, GAL B, Weger J et al. Identification of SNPs linked fpls.2017.01923. to eight apple disease resistance loci. Mol Breed. 2015;35:45. 27. Larsen B, Migicovsky Z, Jeppesen AA et al. Genome-wide asso- https://doi.org/10.1007/s11032-015-0242-4. ciation studies in apple reveal loci for aroma volatiles, sugar 11. Verma S, Evans K, Guan Y et al. Two large-effect QTLs, Ma composition, and harvest date. The Plant Genome. 2019;12:180104. and Ma3, determine genetic potential for acidity in apple fruit: https://doi.org/10.3835/plantgenome2018.12.0104. Jung et al. | 21 28. Hu Y, Han Z, Sun Y et al. ERF4 affects fruit firmness through environmental data. Theor Appl Genet. 2014;127:595–607. TPL4 by reducing ethylene production. Plant J. 2020;103:937–50. https://doi.org/10.1007/s00122-013-2243-1. https://doi.org/10.1111/tpj.14884. 44. Tsai H-Y, Cericola F, Edriss V et al. Use of multiple 29. Minamikawa MF, Kunihisa M, Noshita K et al. Tracing traits genomic prediction, genotype by environment founder haplotypes of Japanese apple varieties: application interactions and spatial effect to improve prediction in genomic prediction and genome-wide association study. accuracy in yield data. PLoS One. 2020;15:e0232665. Hortic Res. 2021;8:49. https://doi.org/10.1038/s41438-021- https://doi.org/10.1371/journal.pone.0232665. 00485-3. 45. Lado B, Vásquez D, Quincke M et al. Resource allocation opti- 30. Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total mization with multi-trait genomic prediction for bread wheat genetic value using genome-wide dense marker maps. Genetics. (Triticum aestivum L.) baking quality. Theor Appl Genet. 2018;131: 2001;157:1819. 2719–31. https://doi.org/10.1007/s00122-018-3186-3. 31. Meuwissen T. Genomic selection: marker assisted selection 46. Gianola D, Fernando RL. A multiple-trait Bayesian LASSO for on a genome wide scale. J Anim Breed Genet. 2007;124:321–2. genome-enabled analysis and prediction of complex traits. https://doi.org/10.1111/j.1439-0388.2007.00708.x. Genetics. 2020;214:305–31. https://doi.org/10.1534/genetics. 32. Kumar S, Chagné D, Bink MCAM et al. Genomic selection 119.302934. for fruit quality traits in apple (Malus × domestica Borkh.). 47. Watanabe K, Stringer S, Frei O et al. A global overview of PLoS One. 2012;7:e36674. https://doi.org/10.1371/journal.pone. pleiotropy and genetic architecture in complex traits. Nat Genet. 0036674. 2019;51:1339–48. https://doi.org/10.1038/s41588-019-0481-0. 33. Kumar S, Molloy C, Muñoz P et al. Genome-enabled 48. Johnston JW, Gunaseelan K, Pidakala P et al. Co-ordination of estimates of additive and nonadditive genetic vari- early and late ripening events in apples is regulated through ances and prediction of apple phenotypes across differential sensitivities to ethylene. J Exp Bot. 2009;60:2689–99. environments. G3 Genes|Genomes|Genetics. 2015;5:2711–8. https://doi.org/10.1093/jxb/erp122. https://doi.org/10.1534/g3.115.021105. 49. Chagné D, Dayatilake D, Diack R et al. Genetic and environ- 34. Muranty H, Troggio M, Sadok IB et al. Accuracy and responses mental control of fruit maturation, dry matter and firmness of genomic selection on key traits in apple breeding. Hortic Res. in apple (Malus × domestica Borkh.). Hortic Res. 2014;1:14046. 2015;2:15060. https://doi.org/10.1038/hortres.2015.60. https://doi.org/10.1038/hortres.2014.46. 35. Roth M, Muranty H, Guardo MD et al. Genomic prediction of 50. Costa F, Stella S, Van de Weg WE et al. Role of the genes fruit texture and training population optimization towards the Md-ACO1 and Md-ACS1 in ethylene production and shelf life application of genomic selection in apple. Hortic Res. 2020;7:148. of apple (Malus domestica Borkh). Euphytica. 2005;141:181–90. https://doi.org/10.1038/s41438-020-00370-5. https://doi.org/10.1007/s10681-005-6805-4. 36. Jung M, Roth M, Aranzana MJ et al. The apple REFPOP— 51. Costa F, Peace CP, Stella S et al. QTL dynamics for fruit firm- a reference population for genomics-assisted breeding in ness and softening around an ethylene-dependent polygalactur- apple. Hortic Res. 2020;7:189. https://doi.org/10.1038/s41438-020- onase gene in apple (Malus × domestica Borkh.). J Exp Bot. 2010;61: 00408-8. 3029–39. https://doi.org/10.1093/jxb/erq130. 37. García-Ruiz A, Cole JB, VanRaden PM et al. Changes in genetic 52. Longhi S, Moretto M, Viola R et al. Comprehensive QTL map- selection differentials and generation intervals in US Holstein ping survey dissects the complex fruit texture physiology in dairy cattle as a result of genomic selection. Proc Natl Acad apple (Malus × domestica Borkh.). J Exp Bot. 2012;63:1107–21. Sci. 2016;113:E3995–4004. https://doi.org/10.1073/pnas.15190 https://doi.org/10.1093/jxb/err326. 61113. 53. Longhi S, Hamblin MT, Trainotti L et al. A candidate gene based 38. Duan N, Bai Y, Chen X et al. Genome re-sequencing approach validates Md-PG1 as the main responsible for a QTL reveals the history of apple and supports a two-stage impacting fruit texture in apple (Malus × domestica Borkh.). model for fruit enlargement. Nat Commun. 2017;8:249. BMC Plant Biol. 2013;13:37. https://doi.org/10.1186/1471-2229- https://doi.org/10.1038/s41467-017-00336-7. 13-37. 39. Howard R, Carriquiry AL, Beavis WD. Parametric 54. Legay S, Guerriero G, Deleruelle A et al. Apple russeting and nonparametric statistical methods for genomic as seen through the RNA-seq lens: strong alterations selection of traits with additive and epistatic genetic in the exocarp cell wall. Plant Mol Biol. 2015;88:21–40. architectures. G3: Genes|Genomes|Genetics. 2014;4:1027–46. https://doi.org/10.1007/s11103-015-0303-4. https://doi.org/10.1534/g3.114.010298. 55. Johnston JW, Hewett EW, Hertog MLATM. Postharvest 40. Cooper M, DeLacy IH. Relationships among analytical softening of apple (Malus domestica) fruit: a review. N Z J Crop methods used to study genotypic variation and genotype- Hortic Sci. 2002;30:145–60. https://doi.org/10.1080/01140671.2002. by-environment interaction in plant breeding multi- 9514210. environment experiments. Theor Appl Genet. 1994;88:561–72. 56. Rymenants M, van de Weg E, Auwerkerken A et al. Detec- https://doi.org/10.1007/BF01240919. tion of QTL for apple fruit acidity and sweetness using sen- 41. Snape JW, Foulkes MJ, Simmonds J et al. Dissecting gene × sorial evaluation in multiple pedigreed full-sib families. Tree environmental effects on wheat yields via QTL and phys- Genet Genomes. 2020;16:71. https://doi.org/10.1007/s11295-020- iological analysis. Euphytica. 2007;154:401–8. https://doi.org/ 01466-8. 10.1007/s10681-006-9208-2. 57. van de Weg E, Guardo MD, Jänsch M et al. Epistatic 42. Lopez-Cruz M, Crossa J, Bonnett D et al. Increased prediction fire blight resistance QTL alleles in the apple cultivar accuracy in wheat breeding trials using a marker × environment ‘Enterprise’ and selection X-6398 discovered and characterized interaction genomic selection model. G3: Genes|Genomes|Genetics. through pedigree-informed analysis. Mol Breed. 2017;38:5. 2015;5:569–82. https://doi.org/10.1534/g3.114.016097. https://doi.org/10.1007/s11032-017-0755-0. 43. Jarquín D, Crossa J, Lacaze X et al. A reaction norm model 58. Evans KM, Patocchi A, Rezzonico F et al. Genotyping of pedigreed for genomic selection using high-dimensional genomic and apple breeding material with a genome-covering set of SSRs: 22 | Horticulture Research, 2022, 9: uhac028 trueness-to-type of cultivars and their parentages. Mol Breed. 70. Gabriel KR. The biplot graphic display of matrices with applica- 2011;28:535–47. https://doi.org/10.1007/s11032-010-9502-5. tion to principal component analysis. Biometrika. 1971;58:453–67. 59. Heslot N, Yang H-P, Sorrells ME et al. Genomic selection in plant https://doi.org/10.1093/biomet/58.3.453. breeding: a comparison of models. Crop Sci. 2012;52:146–60. 71. Tibbs Cortes L, Zhang Z, Yu J. Status and prospects of https://doi.org/10.2135/cropsci2011.06.0297. genome-wide association studies in plants. The Plant Genome. 60. Strobl C, Boulesteix A-L, Kneib T et al. Conditional variable 2021;14:e20077. https://doi.org/10.1002/tpg2.20077. importance for random forests. BMC Bioinformatics. 2008;9:307. 72. Huang M, Liu X, Zhou Y et al. BLINK: a package for the https://doi.org/10.1186/1471-2105-9-307. next level of genome-wide association studies with both 61. Hardner CM, Evans K, Brien C et al. Genetic architecture individuals and markers in the millions. GigaScience. 2018;8. of apple fruit quality traits following storage and implica- https://doi.org/10.1093/gigascience/giy154. tions for genetic improvement. Tree Genet Genomes. 2016;12:20. 73. Tang Y, Liu X, Wang J et al. GAPIT version 2: an https://doi.org/10.1007/s11295-016-0977-z. enhanced integrated tool for genomic association and 62. Hardner CM, Hayes BJ, Kumar S et al. Prediction of prediction. The Plant Genome. 2016;9. https://doi.org/10.3835/ genetic value for sweet cherry fruit maturity among plantgenome2015.11.0120. environments using a 6K SNP array. Hortic Res. 2019;6:6. 74. Ward JH. Hierarchical grouping to optimize an https://doi.org/10.1038/s41438-018-0081-7. objective function. J Am Stat Assoc. 1963;58:236–44. 63. Crossa J, Pérez-Rodríguez P, Cuevas J et al. Genomic https://doi.org/10.1080/01621459.1963.10500845. selection in plant breeding: methods, models, and 75. Breiman L. Random forests. Mach Learn. 2001;45:5–32. perspectives. Trends Plant Sci. 2017;22:961–75. https://doi.org/ https://doi.org/10.1023/A:1010933404324. 10.1016/j.tplants.2017.08.011. 76. Habier D, Fernando RL, Kizilkaya K et al. Extension of the 64. Cuevas J, Crossa J, Montesinos-López OA et al. Bayesian bayesian alphabet for genomic selection. BMC Bioinformatics. genomic prediction with genotype × environment interac- 2011;12:186. https://doi.org/10.1186/1471-2105-12-186. tion kernel models. G3 Genes|Genomes|Genetics. 2017;7:41–53. 77. de los Campos G, Gianola D, Rosa GJM et al. Semi-parametric https://doi.org/10.1534/g3.116.035584. genomic-enabled prediction of genetic values using reproduc- 65. Cirilli M, Micali S, Aranzana MJ et al. The multisite PeachRef- ing kernel Hilbert spaces methods. Genet Res. 2010;92:295–308. Pop collection: a true cultural heritage and international https://doi.org/10.1017/S0016672310000285. scientific tool for fruit trees. Plant Physiol. 2020;184:632–46. 78. VanRaden PM. Efficient methods to compute genomic https://doi.org/10.1104/pp.19.01412. predictions. J Dairy Sci. 2008;91:4414–23. https://doi.org/ 66. Laurens F, Aranzana MJ, Arus P et al. An integrated approach 10.3168/jds.2007-0980. for increasing breeding efficiency in apple and peach in Europe. 79. Wright MN, Ziegler A. Ranger: a fast implementation of random Hortic Res. 2018;5:11. https://doi.org/10.1038/s41438-018- forests for high dimensional data in C++ and R. J Stat Softw. 0016-3. 2017;77:1–17. https://doi.org/10.18637/jss.v077.i01. 67. Browning SR, Browning BL. Rapid and accurate haplotype phas- 80. Pérez P, de los Campos G. Genome-wide regression and predic- ing and missing-data inference for whole-genome association tion with the BGLR statistical package. Genetics. 2014;198:483–95. studies by use of localized haplotype clustering. Am J Hum Genet. https://doi.org/10.1534/genetics.114.164442. 2007;81:1084–97. https://doi.org/10.1086/521987. 81. de los Campos G, Sorensen D, Gianola D. Genomic 68. Rodríguez-Álvarez MX, Boer MP, van Eeuwijk FA et al. heritability: what is it. PLoS Genet. 2015;11:e1005048. Correcting for spatial heterogeneity in plant breeding https://doi.org/10.1371/journal.pgen.1005048. experiments with P-splines. Spatial Statistics. 2018;23:52–71. 82. Lehermeier C, de Los Campos G, Wimmer V et al. Genomic https://doi.org/10.1016/j.spasta.2017.10.003. variance estimates: with or without disequilibrium covari- 69. Bates D, Mächler M, Bolker B et al. Fitting linear mixed-effects ances? J Anim Breed Genet. 2017;134:232–41. https://doi.org/ models using lme4. J Stat Softw. 2015;67:1–48. 10.1111/jbg.12268.