Jung et al. Horticulture Research (2020)7:189 https://doi.org/10.1038/s41438-020-00408-8 Horticulture Research www.nature.com/hortres ARTICLE Open Access The apple REFPOP—a reference population for genomics-assisted breeding in apple Michaela Jung1,2, Morgane Roth2,16, Maria José Aranzana 3,4, Annemarie Auwerkerken5, Marco Bink6,7, Caroline Denancé8, Christian Dujak4, Charles-Eric Durel8, Carolina Font i Forcada3, Celia M. Cantin3,9, Walter Guerra10, Nicholas P. Howard 11,12, Beat Keller1,2, Mariusz Lewandowski13, Matthew Ordidge14, Marijn Rymenants5,6,15, Nadia Sanin10, Bruno Studer1, Edward Zurawicz13, François Laurens8, Andrea Patocchi2 and Hélène Muranty8 Abstract Breeding of apple is a long-term and costly process due to the time and space requirements for screening selection candidates. Genomics-assisted breeding utilizes genomic and phenotypic information to increase the selection efficiency in breeding programs, and measurements of phenotypes in different environments can facilitate the application of the approach under various climatic conditions. Here we present an apple reference population: the apple REFPOP, a large collection formed of 534 genotypes planted in six European countries, as a unique tool to accelerate apple breeding. The population consisted of 269 accessions and 265 progeny from 27 parental combinations, representing the diversity in cultivated apple and current European breeding material, respectively. A high-density genome-wide dataset of 303,239 SNPs was produced as a combined output of two SNP arrays of different densities using marker imputation with an imputation accuracy of 0.95. Based on the genotypic data, linkage 1234567890():,; 1234567890():,; 1234567890():,; 1234567890():,; disequilibrium was low and population structure was weak. Two well-studied phenological traits of horticultural importance were measured. We found marker–trait associations in several previously identified genomic regions and maximum predictive abilities of 0.57 and 0.75 for floral emergence and harvest date, respectively. With decreasing SNP density, the detection of significant marker–trait associations varied depending on trait architecture. Regardless of the trait, 10,000 SNPs sufficed to maximize genomic prediction ability. We confirm the suitability of the apple REFPOP design for genomics-assisted breeding, especially for breeding programs using related germplasm, and emphasize the advantages of a coordinated and multinational effort for customizing apple breeding methods in the genomics era. Introduction aim to create new cultivars adapted to consumer demands Apple (Malus × domestica Borkh.) is one of the most and changing climate, these goals could be difficult to economically valuable fruit crops in temperate regions1. reach within the narrow elite genetic pool of modern Thousands of cultivars are grown in national and private breeding material4. repositories around the globe. Extensive genetic variation Since the advent of genomics, genotyping tools have described in the European apple germplasm illustrates the begun to produce affordable genome-wide marker data. available genetic diversity among cultivars2,3. Only a Large datasets are being analyzed to explore fraction of the existing apple cultivars is commercially genotype–phenotype relationships in genome-wide asso- used. Although ongoing breeding programs worldwide ciation studies (GWAS) and to allow genomic prediction. Particularly, genomic prediction5 has revolutionized breeding and more than doubled genetic progress of Correspondence: Hélène Muranty (helene.muranty@inrae.fr) 1 major livestock such as cattle6. The method relies upon Molecular Plant Breeding, Institute of Agricultural Sciences, ETH Zurich, 8092 models fitted to broad datasets of genotypes and pheno- Zurich, Switzerland 2 Breeding Research group, Agroscope 8820 Wädenswil, Switzerland types from a training population. The aim is to predict the Full list of author information is available at the end of the article © The Author(s) 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Jung et al. Horticulture Research (2020)7:189 Page 2 of 16 agronomic performance of breeding material related to analysis. Our further aim was to discuss the use of the the training set based on the marker information alone. apple REFPOP for genomic prediction in multiple envir- For the application of genomic prediction in fruit trees, onments, across multiple traits and multi-management apple became a model species due to its economic practices and to facilitate the improvement of apple importance and the range of available research resources7. breeding using the established apple REFPOP. Genomic prediction in apple was tested for the first time by Kumar et al. in 20128. However, this and further stu- Results dies have been based on a limited number of genetic Composition of the apple REFPOP markers and/or been carried out at a local scale9–14. The apple REFPOP was ultimately built with (i) 269 The prediction accuracy of genomic prediction models diploid accessions representing a wide range of genetic generally increases with the number of markers used to diversity in apple, originating from various geographic genotype the training population, reaching a plateau regions around the globe and (ii) 265 diploid progeny from depending on the architecture of the trait, the number of 27 parental combinations from several European breeding individuals in the training population, the size of the programs (Supplementary Tables 1 and 2). Similar to a genome and linkage disequilibrium15. Low population previous study2, 194 European accessions were classified structure2 and rapid linkage disequilibrium decay in highly according to their origin into (i) Northern and Eastern diverse apple germplasm3 underlie the need for dense SNP Europe incl. Russia and Baltic countries (NEE, number of marker datasets for GWAS. Considerable progress has accessions n = 28), (ii) Western and Central Europe been made in the development of genomic resources in (WCE, n = 134), (iii) Southern Europe with accessions apple16. Creation of low and medium density SNP arrays, from Spain, Italy, and Portugal (SE, n = 22), and (iv) such as the Illumina Infinium® 20K SNP genotyping array Southeastern Europe containing accessions from Romania, (20K array)17 was followed by the establishment of the Northern Macedonia, Moldova, Bulgaria, and Turkey Affymetrix Axiom® Apple 480K SNP genotyping array (SEE, n = 10). The 69 non-European accessions originated (480K array) with more than 480K SNPs18. Using this SNP from (i) Australia and New Zealand (ANZ, n = 8), (ii) array for apple, markers significantly associated with Canada (CAN, n = 16), (iii) Japan (JPN, n = 9), (iv) United phenological traits in several germplasm collections have States of America (USA, n = 34), and (v) South Africa been successfully discovered by GWAS3. Although the (ZAF, n = 2). Six accessions previously analyzed with the high marker density of the 480K array may lead to higher 480K array but of unknown origin were also included in prediction accuracies than achieved before, no genomic the apple REFPOP. prediction study was conducted using this array so far. Commercial apple breeding programs still cannot afford Marker imputation and validation of the imputation large-scale genotyping of their germplasm with expensive output tools. Therefore, a balance between genotyping density, Missing marker values in the reference set obtained from costs, and predictive ability should be found for the the 480K array (see “Material and methods” section, Table 1) application of genomic prediction in breeding. were inferred through a minor imputation step to obtain a Here, we present an apple reference population: the full dataset of 303,239 SNPs covering the 17 apple chro- apple REFPOP. The population has been replicated across mosomes. This minor imputation step was performed for a six environments in Europe and designed for comparing collection of 1356 accessions (including 269 apple REFPOP two management practices, which will allow for a thor- accessions) and six progeny genotyped in separate stu- ough and unique study of the effects of genotype, envir- dies18,19. Subsequently, a major marker imputation applied onment, and management as well as their interactions on to the 259 of the 265 progeny that were genotyped with the apple phenotypes. The main objectives of this study were 20K array (the remaining 6 progeny with 480K array data to consolidate the high-density SNP marker dataset for all available were included in the reference set) increased the apple REFPOP genotypes, apply the SNP marker dataset marker density from 7060 to 303,239 SNPs (97.7% of marker when describing population characteristics of the apple values imputed). Additionally, imputation accuracy of 0.96 REFPOP and prove suitability of the apple REPOP design and 0.94 when calculating Pearson correlations between for genomics-assisted breeding. Success of the genomics- imputed and original values across individuals and across assisted breeding may depend on characteristics such as markers, respectively, was estimated for 86 genotypes from marker density, trait architecture or size of the training two parental combinations, ‘Fuji’ × ‘Pinova’ and ‘Golden population10,15. For the first time in apple, these aspects Delicious’ × ‘Renetta Grigia di Torriana’. could be tested with (i) GWAS and genomic prediction using the high-density 480K array marker dataset, (ii) a Linkage disequilibrium and population structure comparison of the effects of SNP density on GWAS and From the complete set of genome-wide SNP data for genomic prediction, and (iii) the prediction precision both apple REFPOP accessions and progeny, rapid decay Jung et al. Horticulture Research (2020)7:189 Page 3 of 16 Table 1 Overview of the unique genotypes used in the genotype imputation Number of unique genotypes Set of data Population type SNP array resolution 269 Reference Accessions, apple REFPOP 480K 1,087 Reference Accessions, additional material 480K 6a Reference 1 parental combination, apple REFPOP 480K 259 Imputation 27 parental combinations, apple REFPOP 20K 86b Validation 2 parental combinations, additional material 480K 6 progeny of ‘Fuji’ × ‘Pinova’ a 40 progeny of ‘Fuji’ × ‘Pinova’ and 46 progeny of ‘Golden Delicious’ × ‘Renetta Grigia di Torriana’ b Fig. 1 Linkage disequilibrium decay in the apple reference population. Linkage disequilibrium with a loess smoother for a distances betweenSNPs across the span of chromosomes, and b for SNPs within a 5 kb distance of linkage disequilibrium was found (Fig. 1a). The loess ADMIXTURE analysis revealed that two local minima of smoother fitted to r2 values of SNPs within a 5 kb distance the cross-validation error were reached at the number of dropped below a threshold of 0.2 at distance of 2.52 kb, clusters K = 14 and K = 17 (Supplementary Fig. 1). For the the curve being very flat (Fig. 1b). Average r2 calculated at first minimum (K = 14), genotypes sorted by cluster 100, 5, 1 kb and 100 bp was 0.14, 0.21, 0.21, and 0.24, membership within groups defined by geographic region of respectively. Separate analysis of the accession and pro- origin appeared highly admixed (Fig. 2c). A PCA of the geny group showed a similar pattern of linkage dis- progeny group with 13 parents of the crosses included as equilibrium in both groups (Supplementary Figs. 1 and 2). supplementary individuals showed that the first two prin- In an unrooted neighbor-joining tree of apple REFPOP cipal components described 10.4% of the total variance in accessions (Fig. 2a), the non-European accessions clustered the genomic data and the parents fell among the many at the upper side of the tree with a transition towards the small clusters formed by members of distinct parental European accessions clustering at the lower branches. In combinations (Supplementary Fig. 4). From all population the principal component analysis (PCA) of the accession structure analyses, we concluded that the apple REFPOP group (Fig. 2b), the first two principal components was composed of diverse germplasm with very weak explained only 7.6% of the total variance in genetic markers. population structure together with high levels of admixture. The first principal component showed a slight differentia- tion between European and non-European accessions with Phenotypic analyses the majority of the non-European accessions positioned on Using the spatial analysis of field trials with splines, the right side of the plot. The second principal component spatial patterns were captured appropriately for all envir- displayed a weak latitudinal cline in the European acces- onments and traits as suggested by the homogeneous sions with the southern and northern European accessions residuals (Eq. (1)). Estimated effective dimensions of the placed towards the opposite extremes of the second com- spatial model, which are helpful for characterizing the ponent. Progeny, added to the PCA as supplementary importance of model components20, showed a generally individuals after PCA loadings were estimated, were larger field variation due to spatial effects for floral grouped tightly together among the accessions although emergence than for harvest date (Supplementary Table 3). they did not form any separate cluster (Fig. 2b). The Differences between environments for both floral Jung et al. Horticulture Research (2020)7:189 Page 4 of 16 Fig. 2 Structure of the apple reference population. a Unrooted neighbor-joining tree of the accession group, colors correspond to the legend in “b”. b Principal component analysis of the accession group with progeny group as supplementary individuals encircled with a normal confidence ellipse (constructed using a multivariate normal distribution, level 0.95). Plot of the first two principal components with their respective proportion of the total variance shown within brackets. c ADMIXTURE bar plot of the accession group for K = 14. Labels in plots “a” to “c” refer to the geographic origin of genotypes: ZAF (South Africa), JPN (Japan), ANZ (Australia and New Zealand), CAN (Canada), USA (United States of America), WCE (Western and Central Europe), NEE (Northern and Eastern Europe), SE (Southern Europe), SEE (Southeastern Europe), U (accessions of unknown geographic origin), and P representing the progeny group in plot “b”. In plot “c”, each group of genotypes with a common geographic origin is labeled at its right side emergence and harvest date were visualized using adjusted various environments, with an apparent broader variation phenotypic values of each tree in every environment among genotypes for harvest date than for floral emer- (Fig. 3a, b). Based on this adjusted phenotypic data from gence (Fig. 4a, b). When modeling the effects of envir- 2018, trees began flowering during a narrow period ran- onment, genotype and their interaction (see “Material and ging between 16 days in Switzerland and 47 days in France. methods” section, Eq. (5)), the proportion of variance in On average, the flowering began earlier at the southern floral emergence explained by the environment, genotype, sites (Spain and Italy). By contrast, harvest dates were and genotype by environment interaction was 43%, 22%, distributed across a much longer time span of 99–150 days and 18%, respectively (Fig. 4c). For harvest date, the in Switzerland and Spain, respectively. As with the floral environment explained only 5%, whilst genotype emergence, the mean harvest dates tended to be the ear- explained 74% and genotype by environment interaction liest in the southern European sites. When comparing the explained 12% of the overall variance (Fig. 4c). Confidence adjusted data to the raw data, we found a small systematic intervals for the estimated variances of random effects of increase in individual-location clonal mean heritability for genotype, genotype by environment interaction, and all environments and both traits (Fig. 3c). Values of the residuals were distinct from zero and therefore, the gen- individual-location clonal mean heritability after the otypes differed in the evaluated trait and interacted with adjustment were generally high for both traits and larger the environment (Supplementary Figs. 5 and 6). for harvest date than for floral emergence. Floral emer- gence data from Poland with individual-location clonal Genome-wide association and prediction analyses mean heritability below 0.1 were excluded. Using the The application of GWAS to apple REFPOP dataset adjusted phenotypic values of each tree from the remain- identified three markers associated with floral emergence ing environments to estimate the phenotypic least-square (Fig. 5a, see Supplementary Table 4 for a list of p-values). means of genotypes across environments, the obtained Reported p-values were closely below the log-transformed values were distributed across a narrow period of 33 days Bonferroni-corrected significance threshold for GWAS for floral emergence whereas the values of harvest date performed with the full set of 480K array SNPs; all three were distributed over 97 days (Fig. 3d). associations were significant at density of 150,000 SNPs. Predicted values of genotypes adjusted for spatial het- Two SNPs were located at the top of chromosome 9 erogeneity within each environment showed that indivi- (proportion of explained phenotypic variance r2 of 0.07 dual genotypes appeared to respond differently to the and 0.03) and a third one on chromosome 11 (r2 = 0.10). Jung et al. Horticulture Research (2020)7:189 Page 5 of 16 Fig. 3 Distributions and heritability of the phenotypic data. a, b Violin plots of floral emergence and harvest date for individual environments using the adjusted phenotypic values of each tree. Gray and black circles denote mean and median values, respectively. c Individual-location clonal mean heritability for two analyzed traits with values before and after the correction of spatial heterogeneity. d Density plot of phenotypic least- square means of genotypes across environments with environmental effects removed, calculated from the adjusted phenotypic values of each tree corrected for spatial heterogeneity within environments. The environments were labeled with codes: Belgium (BEL), Switzerland (CHE), Spain (ESP), France (FRA), and Italy (ITA) Fig. 4 Visualization of variability between traits (floral emergence and harvest date), genotypes (both apple REFPOP groups) and environments. a, b Predicted values of genotypes adjusted for spatial heterogeneity within each environment, 30 randomly chosen genotypes were highlighted with colors. Order of the environments corresponds to their latitude. The environments were labeled with codes: Spain (ESP), Italy (ITA), Switzerland (CHE), France (FRA), and Belgium (BEL). c Stacked bar plots with the variance of the fixed effect of environment and the random effects of genotype, genotype by environment interaction, and residuals; calculated from the model following Eq. (5) (see “Material and methods” section) Four SNPs were identified to be significantly associated 3 (r2 = 0.39) with another significantly associated marker with harvest date using the full set of 480K array SNPs at a distance of 5.3 Mb (r2 ≈ 0). Two SNPs on chromo- (Fig. 5b, see Supplementary Table 5 for a list of p-values). some 10 (r2 = 0.15) and on chromosome 16 (r2 = 0.11) The strongest association was found on the chromosome were also significantly associated with harvest date. The Jung et al. Horticulture Research (2020)7:189 Page 6 of 16 Fig. 5 Results of the genome-wide association study (GWAS) and genomic prediction analysis. Manhattan plots for a floral emergence and b harvest date with log-transformed p-values obtained with GWAS and Bonferroni-corrected significance threshold indicated with dashed line. c Genomic predictive ability measured with Pearson correlation coefficient, multi-location clonal mean heritability, and average predictive ability. Comparison of the number of significant associations in d GWAS and e genomic predictive ability measured with the respective mean values and their 95% confidence intervals under various SNP densities obtained through three feature selection strategies (see “Materials and methods” section). Plots “d” and “e” share a common legend placed in “e”. f, g Expected precision of genomic estimated breeding values (GEBVs) with different training population sizes N, number of effective markers M in linkage disequilibrium with the genes underlying the trait and two heritability values h2 of f 0.5 and g 0.8. The N values correspond to a minimum of 10 individuals, and to the sizes of the accession group (N = 269), the whole apple REFPOP (N = 534) and a population of an approximately double size of the apple REFPOP (N = 1000) QQ-plots for both traits indicated a good model fit associations was higher for harvest date than for floral (Supplementary Figs. 7 and 8). emergence at all densities (Fig. 5d). For floral emergence, With the full set of 480K array SNPs, floral emergence the number of significant associations increased up to was predicted using the RR-BLUP model with an average density of 150,000 SNPs. Beyond this density, the number predictive ability of 0.57, whereas the multi-location clo- of significant associations decreased. The number nal mean heritability was equal to 0.86 for this trait remained at zero for the 12K marker set (12,374 SNPs) (Fig. 5c). Average predictive ability of harvest date was thinned for linkage disequilibrium. For harvest date, the 0.75 with the multi-location clonal mean heritability of number of significant associations was low (~1) for SNP 0.97 for the trait. densities of 0.5–1K. Across all densities larger than 1K, Comparison of GWAS performance under various SNP the number of significant associations with harvest date densities showed that the number of significant remained high (≥2). For the 12K marker set thinned for Jung et al. Horticulture Research (2020)7:189 Page 7 of 16 linkage disequilibrium, two significant associations on one training and breeding populations is a crucial factor for chromosome were found. Similar to GWAS, predictive successful genomic prediction13, the European breeding ability of genomic prediction was higher for harvest date programs can directly profit from the outcomes of the than for floral emergence at all densities (Fig. 5e). For project because their breeding material has been included both traits, the plateau in predictive ability was reached at in the apple REFPOP. Furthermore, due to the high a density of 10K SNPs and the difference between pre- admixture of the accessions, which stem from across the dictive ability at 500 and 10,000 SNPs was ~0.1. globe, the levels of diversity in the apple REFPOP appear The precision analysis of genomic prediction (Fig. 5f, g) to be adequate to predict the performance of a broad allowed the assessment of the trait architecture and her- spectrum of novel breeding material. itability ranges that may lead to satisfying precision of genomic estimated breeding values. For a training popu- Expanding genomic information via imputation lation with a size of N = 10, very high precision values Lately, resources for high-density marker genotyping in (precision ≥ 0.8) may be expected only for simple, oligo- apple became available with the development of the 480K genic traits (number of markers underlying the trait M ≤ array18, but using this array remained costly. Here, we 10). With a training population of the size of the apple combined the high-density marker information of the REFPOP accession group (N = 269), more complex traits accessions with the recently inferred pedigrees for with M = 100 may be predicted with very high precision if numerous cultivars19 to perform a large-scale marker the heritability was high (h2 = 0.8), but with a decrease in imputation of the progeny group genotyped with the 20K precision for a moderate heritability (h2 = 0.5). A training array17. Although marker imputation has been routinely population size equal to that of the whole apple REFPOP applied at large scales in well-studied organisms such as (N = 534) may be sufficient to predict traits with M = 100 cattle before23, here we report for the first time a large for both high and moderate heritabilities with very high scale, high-density imputation in apple. The imputation precision. The precision analysis indicated that highly has multiplied the amount of marker information by ~43× precise predictions for traits with very complex archi- for the progeny group at relatively low cost. A high tecture (100 < M ≤ 1000) may be possible if heritability imputation accuracy was achieved, which is comparable was high and the training population had approximately to a similar study in poplar24. To evaluate the imputation twice the size of the apple REFPOP (N = 1000). However, accuracy, two parental combinations ‘Fuji’ × ‘Pinova’ and the increase in estimated precision from N = 534 to N = ‘Golden Delicious’ × ‘Renetta Grigia di Torriana’ consist- 1000 was comparably smaller than the increase from N = ing together of 86 progeny with known high-density 269 to N = 534. marker data were used. The parents ‘Fuji’ and ‘Pinova’ are present in the pedigree of five and three parental com- Discussion binations of the imputation set, respectively. The parent Composition of the apple REFPOP ‘Golden Delicious’ can be found in all of the pedigrees This study is the first report on the apple REFPOP, a from the imputation set, and therefore, relates to the reference population created to advance apple breeding whole progeny group of the apple REFPOP. Although the with genomics-assisted methods. The population is divi- pedigree of ‘Renetta Grigia di Torriana’ is unknown, both ded into an accession and a progeny group. Genotypes of parental combinations ‘Fuji’ × ‘Pinova’ and ‘Golden Deli- its accession group (all diploid) were chosen to represent cious’ × ‘Renetta Grigia di Torriana’ are strongly related to a wide range of geographic origins and capitalize on the the imputation set and therefore may provide a useful previously available high-density SNP array data. The estimate of the imputation accuracy in other imputed progeny group stemmed from eight different European parental combinations. An attempt to improve the breeding programs to represent current tendencies in imputation accuracy might include expanding the refer- European apple breeding. For the progeny group, the ence set with all 19 parents of the parental combinations choice of genotypes was mainly based on genetic distance as well as some of their ancestors for which high-density sampling within progeny groups derived from a number marker genotype data are not yet available. of parental combinations. On the contrary, the choice of genotypes for the accession group did not correspond to a Insights into apple population genomics core collection design which generally aims at maximizing Earlier studies with lower density of genome-wide SNPs genetic diversity21. The incorporation of both cultivated found that their number of markers was insufficient for accessions and modern breeding material is expected to GWAS because of a rapid decay in linkage disequilibrium promote the applicability of the apple REFPOP for in apple10,25. In concordance with previous findings, we genomics-assisted breeding beyond the scope of most found a strong linkage disequilibrium decay over short current breeding programs, whose genetic bases are distances in the apple REFPOP dataset, the pattern of generally narrow22. Also, since relatedness between decay being very similar to that described by Urrestarazu Jung et al. Horticulture Research (2020)7:189 Page 8 of 16 et al.3 The linkage disequilibrium decay was observed in factor NAC18.1 listed as gene MD03G1222600 on the the apple REFPOP from a distance as short as ~2.5 kb. GDDH13 v1.1 genome26. Other studies have identified Given the genome size of a doubled haploid derivative of associations between NAC18.1 and harvest date, and the ‘Golden Delicious’, estimated to 651 Mb26, and the num- gene is a known member of a family of conserved tran- ber of markers in this study, our marker resolution cor- scriptional regulators involved in ripening3,10. Another responded to one marker per ~2.1 kb on average. Based SNP on chromosome 3 associated with harvest date was on this value, each group of loci in linkage disequilibrium located ~1 Mb upstream from a marker reported before should be on average represented by one marker. Thus, by Urrestarazu et al.3 The remaining two SNP associa- the marker coverage we report seems appropriate for tions were found on chromosomes 10 and 16 where QTL genome-wide analyses, which is supported by our GWAS for harvest date have also been discovered before3,31. The and prediction analyses. Due to the heterogeneity in the identification of SNPs associated with both phenology marker density over the genome, the statistical power traits in well-characterized genomic regions and for one could be reduced in specific regions with low SNP density. possible minor QTL indicates the suitability of the apple Partial differentiation between European and non- REFPOP for discovering other alleles with large and small European accessions in our material is consistent with effects on trait variability in apple. Although novel significant differences between Old World and New marker–trait associations could not be identified and their World varieties found by Migicovsky et al.10 The weak number was lower than in Urrestarazu et al.3, who also latitudinal cline in the European accessions may reflect used genotypic data of the 480K array, the number of the population structure reported in this germplasm by reported marker–trait associations in this work was Urrestarazu et al.2 The overall weak population structure higher than in GWAS studies using lower SNP densities in the apple REFPOP accession group, which was similar obtained by genotyping by sequencing10,32. to a previously described weak genetic structure in apples Applying a genomic prediction model with cross-vali- of European origin, presumably reflected the generally dation, we were able to predict both phenology traits highly admixed apple germplasm that is a result of the with moderate to high predictive ability when compared prominent gene flow characteristic of the cultivated apple to the predictions for different apple traits reported gene pool2,27. previously8–11,13, although cross-validation may have inflated predictive ability compared to a potential inde- An efficient design for genomics-assisted breeding pendent validation with a test set33. To our knowledge, We applied GWAS using the multi-locus mixed-model floral emergence has not been predicted in apple before; (MLMM) method28, which accounts for potential con- in this work, an average predictive ability of 0.57 was founding effects of kinship and population structure. reached for this trait. The predictions for harvest date Thus, this method allows to combine highly related plant had an average predictive ability of 0.75, which was material, like in a pedigree-based-analysis, and more higher than any previously reached accuracies of this, or diversified, unrelated individuals to integrate different equivalent, traits in apple10,11. The presented genomic levels of linkage disequilibrium in the same analysis, prediction methodology may be directly applied for the hence maximizing the chances of finding regions asso- breeding of floral emergence and harvest date. In parti- ciated with the target trait. By performing GWAS on the cular, breeding programs using germplasm related to the apple REFPOP dataset, SNPs associated with phenotypic apple REFPOP may capitalize on this work. variation could be identified and assigned to previously The high SNP density in this study allowed for powerful known genomic regions. A genomic region located on GWAS and genomic prediction analyses, with overall chromosome 9 was significantly associated with the trait lower performance for floral emergence than harvest date. floral emergence. Quantitative trait loci (QTL) at close Floral emergence appeared under weak genetic control locations were previously identified for floral bud break, with the majority of phenotypic variance explained by the floral emergence, and flowering period3,29,30. The second effects of environment and genotype by environment SNP association with floral emergence was located on interaction. Markers significantly associated with the trait chromosome 11 approximately 1 Mb downstream from a explained a low proportion of the phenotypic variance in SNP identified by Urrestarazu et al.3 This association may our study, pointing to a complex genetic architecture of be related to a minor QTL on chromosome 11 discovered floral emergence with many influential genomic regions using best linear unbiased predictors for genotype by year yet uncovered. With SNP density increasing up to 150,000 interaction effect of one season for the mapping of QTLs SNPs, increasing amount of the phenotypic variance of associated to bud break date (a trait highly correlated to floral emergence can be explained with GWAS. floral emergence) in the study of Allard et al.30 For harvest Environment-specific GWAS or GWAS with phenotypes date, one significantly associated SNP on chromosome 3 from across several seasons and locations may allow for was located 14,610 bp upstream from the transcription improved GWAS performance in floral emergence. For Jung et al. Horticulture Research (2020)7:189 Page 9 of 16 harvest date, a trait under strong genetic control with a harvest date reported in sweet cherry34. The replication of large proportion of the phenotypic variance explained by a the apple REFPOP across six environments will enable the few major genomic regions, our results suggest that SNP inclusion of these interactions into genomic prediction density as low as 10,000 markers may be required to models. Furthermore, GWAS across separate environ- discover genomic associations with the trait. The 10,000 ments can be performed in the future to identify marker subset was also sufficient to reach a plateau of the environment-specific associations and evaluate the stabi- predictive ability in genomic prediction of both traits, with lity of associations across environments. the difference in predictive ability between traits possibly The apple REFPOP was also designed for comparing attributable to the trait architecture. As previously shown different management practices: one part of every orchard in American cranberry15, the SNP density at which the was grown under the conventional practice of each plateau of predictive ability is reached may be specific to region, to evaluate the response of the germplasm to the apple REFPOP independently of the traits, but may environmental effects; the second part can be managed in differ for other populations of the same species. These order to evaluate response to managements such as results can have a practical impact on apple breeding, reduction of irrigation or pesticide application. However, where cost-effective genotyping of 10K SNPs may be the second management practice has not been applied so sufficient for precise genomic prediction. In breeding far to allow the trees to mature. Incorporating genotype research, the same genotyping coverage may be adequate by management and genotype by environment by man- to perform GWAS of oligogenic traits. For complex traits agement interactions in prediction models may help select under strong environmental control such as floral emer- new material adapted to drier climates or with stronger gence, high-density SNP marker datasets remain desirable resistance to apple pests and diseases. in GWAS. Future modeling including the effects of In addition to the multi-environment and multi- environment and genotype by environment interaction management design, protocols for phenotyping of var- may contribute to a higher precision of genomic predic- ious traits have been applied since 2018 to evaluate the tion in complex traits. apple REFPOP. More than 10 different traits including In our study, the size of the SNP marker dataset seemed yield, fruit quality and phenology are being simulta- to affect the genomic prediction ability stronger than the neously phenotyped with the same method at the six feature selection method. On the contrary, a higher power environments. Genotypes are replicated at least twice at of GWAS to detect marker–trait associations can be each of the environments; each tree is evaluated for all obtained with a smaller SNP marker dataset than with a traits separately to allow for variance decomposition up to larger set of SNPs depending on the representation of the individual level (i.e., the tree). Using multi-trait genes underlying the traits in the SNP marker dataset. genomic prediction models, prediction of traits with low This way, more marker–trait associations could be heritability or labor-intensive phenotyping can be sup- revealed with the 7060 SNPs of the 20K array than using ported by genetically correlated traits with higher herit- SNPs pruned for linkage disequilibrium and such pruning ability and available phenotypes35,36. should be avoided prior to GWAS. In the light of prediction precision analysis, the lower Conclusion predictive ability for floral emergence may be explained by This study benefits from a collaborative European a higher complexity of this trait, together with lower approach dedicated to the improvement of apple breeding heritability. Importantly, the results were very likely via genomics-assisted methods. A reference population, impacted by the low phenotypic variability of floral which sampled diverse apple germplasm and current emergence in the season of 2018. Additional phenotyping European breeding material at larger scales than most seasons will likely contribute to a better representation of current breeding programs, has been established. An the flowering variability in the population. extensive set of high-density SNP marker data has been assembled via cost-effective validated marker imputation Prospect of the apple REFPOP for multi-environment, while making use of the recently available SNP arrays and multi-management and multi-trait testing pedigrees. The imputation method of localized haplotype The performance of breeding material in tested as well clustering together with the consolidated high-density as untested but similar environments can be predicted SNP marker dataset can be implemented as a standard for accurately using genomic prediction models taking into cost-effective genomics-assisted breeding. Our diversity account genotype by environment interactions34. We and quantitative genetics analyses showed that the refer- found a moderate but noticeable effect of the genotype by ence population is representative of the current apple environment interaction on both phenology traits eval- diversity and breeding material, and that the associated uated in 2018, which contrasts with the limited effect of genotypic resources and experimental design allow for the genotype by environment interaction on a trait similar to development and application of genomics-assisted Jung et al. Horticulture Research (2020)7:189 Page 10 of 16 breeding methods in apple. This work emphasized the Wageningen University and Research (WUR), Wagenin- positive effects of high marker density on GWAS and the gen, The Netherlands. Availability of budwood from trees role of trait architecture in both GWAS and genomic affected the numbers of chosen accessions. Therefore, prediction. The apple REFPOP with its unique multi- accessions often had to be either (i) replaced with a dif- environment, multi-management and multi-trait design ferent accession of the same MUNQ or (ii) excluded in represents a rich source of data for future environment- case no other accession corresponding to the unique specific GWAS and genomic predictions. Particularly, the genotype was available or could provide enough budwood. predictions produced with models accounting for the Additionally, the triploid accession ‘Biesterfelder Renette’ interaction effects between genotype, environment, and (MUNQ 1106.1) and the accession ‘Karinable’ (MUNQ management, as well as using multi-trait modeling can 7828) with no available SNP data were excluded from the help untangle the effects underlying the traits and ulti- analysis, although planted in the orchards. mately improve the efficiency and success of apple The progeny group of the apple REFPOP included 27 breeding. The apple REFPOP will be the cornerstone of full-sib parental combinations previously used in the Eur- many future projects including the application of genomic opean project FruitBreedomics7. These full-sib parental selection in apple and the work on apple in the EU- combinations originated from eight different breeding H2020-INVITE project (2019–2024). programs and they were obtained from 32 parents while 13 of the parents were included in the accession group. For Materials and methods most parental combinations, genotypic information was Composition of apple REFPOP plant material available ahead of this study, with different genotyping The apple REFPOP collection composed of accessions density depending on the parental combination. Twenty- and progeny was designed and established by the colla- two parental combinations had been genotyped with the borators of the FruitBreedomics project7. The accession 20K array7,17,38. In addition, in the frame of a pilot study of group consisted of old and modern diploid accessions genomic selection9, three other parental combinations representing a wide range of genetic diversity in apple. were genotyped with a custom 512 SNPs array using the Simple sequence repeat data obtained by Fernández-Fer- TaqMan OpenArray technology covering the whole gen- nández37, Lassois et al.21 and Urrestarazu et al.2 were used ome at a very low density. Finally, for two parental com- to allocate unique genotype code (so-called MUNQ, for binations, whole genome data were not available. Malus UNiQue genotype code as described by Muranty A subset of 10 individuals per parental combination was et al.19) to the accessions included in these studies, which chosen to form the progeny group of the REFPOP. For the resulted in 1292 unique genotypes available for further 25 parental combinations for whom the whole genome choice. When possible, passport data of the accessions marker data was available, 10 individuals were chosen belonging to each unique genotype were used to identify using a genetic distance sampling strategy39. In all cases, its country or region of origin. To form the accession the 10 individuals formed the center of each of the clus- group, a subset of the unique genotypes was created. All ters defined by genetic distance. Where budwood was not possible unique genotypes were chosen when <15 geno- available in sufficient numbers, they were either (i) types were available per each country or region of origin. replaced by individuals closest to the center of the same For origins represented by a larger number of accessions, cluster or (ii) excluded from the progeny group leading to priority was given to the genotypes already analyzed with fewer than 10 genotypes per parental combination being the Axiom®Apple480K array18. Additional genotypes of chosen (as in the case of parental combinations X6679 × each origin were chosen randomly so that the overall X6808 and X6679 × X6417). From the parental combi- number of selected genotypes per origin was proportional nation ‘Jonathan’ × ‘Prima’, 11 progeny were included in to the number of genotypes of this origin in the whole the REFPOP. One progeny of the ‘Dalinette’ × X6681 collection of unique genotypes. Unique genotypes pre- parental combination (NOVADI/0830) was found as tri- viously analyzed with the 480K array but of unknown ploid using 20K genotypic data and thus excluded from origin were also included. Additionally, five accessions the analysis, although planted in the orchards. For the two (‘Red Winter’, ‘O53T136’, ‘Priscilla-NL’, ‘P7 R4A4’, and parental combinations without whole genome marker X6398) considered as founders in the progeny group data available, individuals were chosen randomly. pedigree were added to the accession group. Most acces- sions were chosen from the apple germplasm collection of Multiplication of plant material and planting design the National Fruit Collection, Brogdale, United Kingdom In 2015, budwood from each apple genotype was col- in order to simplify collection of the plant material. The lected and grafted onto ‘M9’ rootstocks. The grafting was budwood of the founders stemmed from two sources: performed in three different nurseries, i.e., at (i) INRAE France’s National Research Institute for Agriculture, Food Angers, France, (ii) Better3Fruit, Rillaar, Belgium, and (iii) and Environment (INRAE), Angers, France and Consorzio Italiano Vivaisti, San Giuseppe, Italy. The Jung et al. Horticulture Research (2020)7:189 Page 11 of 16 following year, grafted trees were planted across six reconstructed in a recent study using all genotyped contrasting environments, each located in (i) Rillaar, diploid unique genotypes19. For more details about SNP Belgium, (ii) Angers, France, (iii) Laimburg, Italy, (iv) filtering, see Supplementary Methods 1. The SNP posi- Skierniewice, Poland, (v) Lleida, Spain, and (vi) Wädens- tions consistent with the apple reference genome based wil, Switzerland. The environments ranged across several on the doubled haploid GDDH13 v1.126 were used. biogeographical regions in Europe, i.e., the Mediterranean Markers unassigned to the 17 apple chromosomes were in Spain, Atlantic in France and Belgium, Alpine in Italy excluded, resulting in a dataset of 303,239 biallelic SNPs. and Continental biogeographical region in Switzerland Genomic data for the apple REFPOP progeny group and Poland40. Depending on the environment, planting were generated using the 20K array. For 210 progeny from distance within and between tree rows ranged from 0.9 to 22 parental combinations, the data were already avail- 1.3 m and 3.2–3.6 m, respectively. All orchards were able7,17,38. The 49 remaining progeny from parental divided into two parts: the first part was to be managed combinations X338 × ‘Braeburn’, ACW 11303 × ACW with the common agricultural practice of the country, the 18522, ACW 13652 × ACW 11567, ‘Dalinette’ × X6681, second was to receive alternative management conditions and X6398 × ‘Pinova’ were genotyped with the 20K array (e.g., low pesticide or water input). So far, both parts have within the framework of this study. All 20K array SNPs been managed in the same way. The first part consisted of were filtered and the allele data for the remaining SNPs two randomized complete blocks, totaling together 1068 were curated to ensure the data made logical marker trees, each block containing one replicate per genotype. inheritance and co-segregation patterns following the The second part also consisted of two randomized com- methods and principles described by Vanderzande et al.41 plete blocks, each containing one representative of The set of 20K-array-generated markers was aligned to approx. one-third of the genotypes (on average 184 gen- the 480K array marker set and 7060 of the 20K array SNPs otypes). Every environment shared about half of these were retained for further analysis. genotypes (92 on average) with one of the other countries (country pairs: Belgium–Italy, France–Switzerland, Marker imputation Poland–Spain). Additionally, each environment shared on Due to the difference in marker resolution of the SNP average 34 genotypes with three further countries. The arrays that were used to generate the SNP marker datasets, blocks within parts were used to ensure that two repli- marker imputation was performed to provide high-density cates of the same genotype were not planted at a close SNP marker information across the whole apple REFPOP. proximity to each other, but were not intended to block Genotypes included in the apple REFPOP and additional for environmental effects. The two parts were physically genotypes involved only in the imputation were used either separated by a row of trees comprised of additional as (i) reference for the imputation or (ii) imputation set or representatives of genotypes included in the apple (iii) validation set (Table 1). First, a reference set of 480K REFPOP as well as other material. Since the alternative array data was formed from 1356 accessions and six pro- management regimes were not applied in the initial years, geny of the cross ‘Fuji’ × ‘Pinova’ that had been previously the trees in the barrier, as well as trees in what would later used for validation of the 480K array and the analysis of be the alternative management scenario may have con- apple pedigrees18,19. Second, the imputation set was tributed phenotypic values for traits in 2018. Three con- formed of the remaining 259 progeny from 27 parental trol genotypes ‘Gala’, ‘Golden Delicious’, and Modì® combinations (i.e., including the remaining four progeny ‘CIVG198’ were each replicated 48 and 18 times in the from the parental combination ‘Fuji’ × ‘Pinova’ that were first and the second part of the orchard, respectively. not included in the validation set). Third, 40 additional progeny of the parental combination ‘Fuji’ × ‘Pinova’ and 46 Molecular marker genotyping progeny of ‘Golden Delicious’ × ‘Renetta Grigia di Torriana’ A SNP marker dataset for the apple REFPOP was pro- which had all been genotyped using the 480K array18, but duced using two overlapping SNP arrays of different none of which had been chosen for inclusion in the apple resolution. The genotypic data generated with the 480K REFPOP, were designated as a validation set. array for a total of 1356 unique genotypes were retrieved The imputation was performed with the localized hap- from previous studies3,19. The SNP marker dataset lotype clustering implemented in the software Beagle 4.0 included all apple REFPOP accessions (including 13 par- using pedigree information42. In the first step, reference ents of the apple REFPOP progeny group) and additional genotypes (Table 1, reference set of data) along with the unique genotypes later used as a reference set for marker recently inferred pedigrees19 were supplied to the pro- imputation. The applied filtering strategy differed from gram for the inference of haplotype phase and minor the original one described by Bianco et al.18 in discarding marker imputation in the reference set. To prepare data the quality prediction based on metrics of the SNP clus- for the second step of imputation, SNP density of the ters and making use of the pedigree, which was validation set (480K array data, see also Table 1, validation Jung et al. Horticulture Research (2020)7:189 Page 12 of 16 set of data) was decreased to the density of the imputation Phenotype scoring set (20K array data, see also Table 1, imputation set of Two phenotypic traits were evaluated at the six plan- data) in order to spike the imputation set with known tation sites in 2018. Floral emergence of each tree was samples of the validation set. In the second step of the recorded as the date when 10% of flowers opened47. imputation, phased reference genotypes along with the Harvest date was measured as the date when more than pedigrees were used to impute missing marker values in 50% of the fruits reached full physiological maturity, as the extended imputation set (i.e., both imputation set and determined by iodine coloration or expert knowledge47. reduced-density validation set samples). Imputation Both traits were evaluated for each of the replicate trees accuracy in the validation set was then evaluated by individually. Deviation from the phenotyping protocol for computing the Pearson correlation between the imputed harvest date led to exclusion of harvest date recorded in and original high-density SNP genotypic values in the Poland. After measurements, dates were converted to validation set controls. counts of days starting at the beginning of the year in which they were measured. Genomic data analyses Linkage disequilibrium Phenotypic data analyses Linkage disequilibrium statistics were calculated as a Raw phenotypic values were corrected for spatial het- square of the correlation coefficient (r2) between pairs of erogeneity individually within environments to obtain the SNPs on each chromosome with the R package snpStats43. adjusted phenotypic values of each tree. The corrected as To reduce computational time, SNPs of each chromo- well as uncorrected (raw) phenotypic values were used to some were sampled randomly to include one tenth of the estimate the individual-location clonal mean heritability. markers per chromosome in the r2 calculation, which The adjusted phenotypic values of each tree were further resulted in 4.6 × 108 marker combinations. A loess used to fit a mixed model including the effects of geno- smoother (α = 0.5) was fitted to 100,000 randomly chosen type, environment and their interaction. The variance of r2 values across the whole span of chromosomes. Addi- each effect was calculated from the fitted mixed model to tionally, the loess smoother (α = 0.5) was calculated for all estimate multi-location clonal mean heritability and the obtained r2 values for pairs of SNPs within a 5 kb distance. fraction of phenotypic variation associated with the To determine the distance between SNPs at which the r2 effects. Finally, the adjusted phenotypic values of each tree dropped below 0.2, average r2 was calculated at 100, 5, 1, were used to obtain phenotypic least-square means of and 0.1 kb (100 bp) as a mean of all r2 within a window of genotypes across environments, i.e., a single mean phe- 100 bp around each of the values. notypic value for each genotype. Population structure Correction of spatial heterogeneity The neighbor-joining method as implemented in the R To account for spatial variation of the complete block package ape44 was used to estimate and visualize an design, e.g., due to different soil composition or water unrooted neighbor-joining tree. Principal component availability in the orchards, and to predict adjusted phe- analysis (PCA) with supplementary individuals was per- notypic values of each tree, spatial heterogeneity in the formed with the R package FactoMineR45. Tenfold cross- phenotypic data was modeled separately for each envir- validation for the number of populations K = {1, 2, …, 20} onment and trait using the spatial analysis of field trials was performed with ADMIXTURE 1.3, a program for with splines (SpATS)20. To specify the smooth compo- estimating ancestry in unrelated individuals46. ADMIX- nent, a two-dimensional penalized tensor-product of TURE was used with default settings and a subset of marginal B-spline basis functions based on the P-spline markers filtered according to the program’s manual. To ANOVA approach (PS-ANOVA) was defined with the avoid spurious effects of high linkage disequilibrium default settings as a function of covariates plantation row between adjacent markers, SNPs were removed in sliding and column, further denoted as f(u, v). The following windows of 50 SNPs advanced by 10 SNPs when squared linear mixed model was fitted correlation was >0.1 for pairs of variants, leading to a subset of 12,374 SNPs. The correct K value was identified y ¼ f ðu; vÞ þ Z g cg þ Z r cr þ Z c cc þ ε; ð1Þ at the lowest cross-validation error. As structure is best estimated among unrelated (or weakly related) indivi- with y being the vector of phenotypic values measured for duals, the neighbor-joining tree and ADMIXTURE were each tree, u and v denoting the numeric positions, i.e., estimated for the accession group only. For the population rows and columns, the vectors cg, cr and cc being the structure analyses, the European accessions were divided random effect coefficients for the genotypes, rows and into several broad regions of origin to compensate for the columns (as factors) associated with the design matrices uncertainty around the exact origin of old varieties2. Zg, Zr and Zc, respectively, the ε denoting the random Jung et al. Horticulture Research (2020)7:189 Page 13 of 16 error vector. Spatial independence assumption of the For the remaining environments, multi-location clonal error vector was visually assessed using a residuals’ spatial mean heritability was estimated for each trait with a pooled plot for each model fit as described by Rodríguez-Álvarez analysis across environments using mixed-effects models et al.20 To characterize the importance of model fitted via restricted maximum likelihood (R package lmer48) components, effective dimensions associated with each random factor (cg, cr, and cc), the PS-ANOVA spatial yijk ¼ μ þ gi þ lj þ glij þ εijk ; ð5Þ trend (fv(v), fu(u),uhv(v),vhu(u) and fu,v(u,v)) and the total effective dimension EDs (sum of partial effective dimen- where yijk was the kth adjusted phenotypic value of each sions associated with each component of the PS-ANOVA tree from genotype i in environment j, µ was the grand spatial trend) were assessed20. Predicted values of mean, gi was the random effect of the ith genotype, lj was genotypes adjusted for spatial heterogeneity within each the fixed effect of the jth environment, glij the interaction environment were produced with random model terms of effect between the ith genotype and jth environment the smooth component and the random row and column (random) and εijk was the error term. Then, multi-location effects excluded from the predictions. These values were clonal mean heritability was calculated using Eq. (3) with used to visualize variability between traits and sites. the phenotypic variance estimated as Residuals for each tree (i.e., each replicate of a genotype) were extracted from the model fit and summed with the σ 2gl σ 2ε σ 2p ¼ σ 2g þ þ ; ð6Þ corresponding predicted values of genotypes adjusted for nl nl nr spatial heterogeneity within each environment to obtain adjusted phenotypic values of each tree. where σ 2gl was the genotype by environment interaction variance and nl the number of environments. Broad-sense heritability and phenotypic variation From the model fit according to Eq. (5), variance of each To evaluate the efficiency of the spatial correction random effect was calculated. The fraction of phenotypic method when comparing raw data with data adjusted for variation associated with the fixed effect was estimated as spatial heterogeneity, individual-location clonal mean the variance of the vector of values predicted from the heritability H2 was estimated for each of the traits and model fit when all random effects were set to zero. To environments before and after the spatial correction. A assess the precision of the random effects, confidence random-effects model was fitted for each environment via intervals for the variance components were estimated restricted maximum likelihood (R package lmer48) using profiling likelihood method. yik ¼ μ þ gi þ εik ; ð2Þ Phenotypic least-square means Phenotypic least-square mean of each genotype across where yik was the kth phenotypic value from genotype i environments was estimated from the adjusted pheno- (adjusted and non-adjusted phenotypic values of each typic values of each tree corrected for spatial hetero- tree), µ was the grand mean, gi was the random effect of geneity within the environments. First, a multiple linear the ith genotype and εik was the error term. Individual- regression model was fitted for each trait location clonal mean heritability was calculated from variance components of the model as total genotypic yijk ¼ μ þ gi þ lj þ εijk ð7Þ variance σ 2g over the phenotypic variance σ 2p where yijk was the kth adjusted phenotypic value of each tree from genotype i in environment j, µ was the grand σ 2g H2 ¼ ; ð3Þ mean, gi was the effect of the ith genotype, lj was the effect σ 2p of the jth environment and εijk was the error term. Second, phenotypic least-square means (or LS-means) of genotypes where the phenotypic variance was calculated from the across environments were calculated with the R package genotypic variance, error variance σ 2ε and the mean doBy49. The phenotypic least square means of genotypes number of replications nr . across environments were used later for genome-wide association analyses and genomic prediction. σ 2ε σ 2p ¼ σ 2g þ ; ð4Þ nr Genome-wide association and prediction analysis Genome-wide association study The individual-location clonal mean heritability was The multi-locus mixed model (MLMM) method28, a used to eliminate one trial with the heritability value stepwise mixed-model regression for mapping complex below 0.1. traits under population structure, was applied to perform Jung et al. Horticulture Research (2020)7:189 Page 14 of 16 GWAS for both apple REFPOP groups together. As a subsets were built for densities of 500, 1000, 5000, 10,000, response variable, the phenotypic least-square means of 50,000, 100,000, 150,000, 200,000, and 250,000 markers. As genotypes across environments were used. MLMM was suggested by Bermingham et al.52, the markers were chosen used as implemented in the R package GAPIT 3.050 with evenly spaced from a random starting point. Since the 480K default settings and SNPs with low minor allele frequency array was designed to cover not only the genic regions, but (0.05) removed. Marker–trait associations were found as also to reach a uniform coverage of non-genic regions18, the significant for p-values falling below a Bonferroni- order of the markers on chromosomes was used as a proxy corrected significance threshold α* = α/m with α = 0.05 for their physical distance when sampling. For each of the and m representing the number of tested markers. The nine SNP densities, the SNP choice was repeated with 10 proportion of phenotypic variance explained by each SNP different seeds resulting in 90 subsets. For these subsets, marker significantly associated to the phenotypic least- GWAS and genomic prediction were performed as descri- square means of genotypes across environments was bed above, with the five-fold cross-validation of the genomic estimated as a coefficient of determination (r2). The r2 was prediction being repeated 10 times for each subset. The estimated from a simple linear regression model fitted number of significant associations detected by GWAS was using the numeric marker values of a single SNP as pre- determined as the number of regions containing SNPs with dictor and phenotypic least-square means of genotypes p-values falling below a Bonferroni-corrected significance across environments as response. threshold. Mean number of the significant associations in GWAS and mean prediction ability of the genomic pre- Genomic prediction diction as well as their 95% confidence interval were cal- The RR-BLUP model51 was used for genomic prediction culated for each number of chosen markers. of breeding values. The model was defined as Analysis of prediction precision y ¼ WGu þ ε ð8Þ Expectations of the precision of genomic estimated breeding values were approximated based on Eq. (9), i.e., with the y being the vector of phenotypic least-square Eq. (5) in the original article by Elsen53. The approxima- means of genotypes across environments, W the design tion of the precision (^r 2 ) of genomic estimated breeding matrix relating genotypes to y, G the SNP matrix, u  values was determined for different values of parameters Nð0; Iσ 2u Þ the vector of SNP marker effects, and ε the vector influenced by the experimental design (i.e., species, of errors. Five-fold cross-validation was applied with the population size and composition, SNP density, environ- model and data of the whole apple REFPOP, each run ments, or trait architecture) as masking 20% of the genotypes as validation set. Marker effects (BLUPs) were generated from the remaining 80% of Nh2 ^r 2 ffi ð9Þ the genotypes to make predictions for the validation set. Nh2 þ Mð1  h2 Þ Predictive ability was estimated as an average value of Pearson correlation coefficient calculated between observed with M the number of loci in linkage disequilibrium with phenotypic least-square means of genotypes across envir- genes underlying the trait, N the population size, and h2 the onments from each validation set and predictions for the heritability. Parameter values were chosen to encompass same genotypes. The cross-validation was repeated 100 and extrapolate beyond the apple REFPOP design: N taking times, the five folds being chosen randomly without values of 10, 269, 534, and 1000, M between 1 and 106 and replacement before each of the repetitions. h2 being equal to 0.5 (moderate) or 0.8 (high). To interpret the output, trait architecture was classified using M into Comparison of GWAS and genomic prediction performance oligogenic (M ≤ 10), complex (10 < M ≤ 100), and very under various SNP densities complex (100 < M ≤ 1000). Precision of genomic estimated Additionally to the analyses with the full set of 303,239 breeding values was considered very high when equal or SNPs produced with the 480K array, different subsets of larger than 0.8. All statistical analyses and data formatting SNPs were used to perform the GWAS and genomic pre- in this article were performed with R54 and visualized with diction to investigate the effect of feature selection approa- the R package ggplot255, unless stated otherwise. ches on model performance. SNPs in subsets were chosen according to three main feature selection strategies, (i) the Acknowledgements The authors wish to thank Dr. Eric van de Weg (Wageningen University and SNP set from the 20K array (7060 SNPs) available in the full Research, The Netherlands (WUR)) and Dr. Mario Di Guardo (WUR and SNP marker dataset, (ii) the SNP set thinned according to Fondazione Edmund Mach, Italy; presently Dipartimento di Agricoltura, linkage disequilibrium (12,374 SNPs, for details see popu- Alimentazione e Ambiente, University of Catania, Catania, Italy) for making raw genotypic data available. Additional genotyping was supported by an INRAE lation structure analysis described above), and (iii) an SelGen grant (project named GDivSelGen: “Efficient use of genetic diversity in unsupervised SNP choice. For this latter strategy, data genomic selection”). The authors thank the field technicians and staff at INRAE Jung et al. Horticulture Research (2020)7:189 Page 15 of 16 Experimental Unit (UE Horti), Angers, France and technical staff at other apple 4. Kumar, S., Volz, R. K., Chagné, D. & Gardiner, S. in Genomics of Plant Genetic REFPOP sites for the maintenance of the orchards. We thank (i) INRAE Resources: Vol. 2. Crop Productivity, Food Security and Nutritional Quality (eds Experimental Unit (UE Horti), Angers, France, (ii) Better3Fruit, Rillaar, Belgium, Tuberosa, R., Graner, A. & Frison, E.) 387–416 (Springer Netherlands, Dordrecht, and (iii) Consorzio Italiano Vivaisti, San Giuseppe, Italy for grafting of the plant 2014). material. This work was partially supported by the project RIS3CAT (COTPA- 5. Meuwissen, T. H. E., Hayes, B. J. & Goddard, M. E. Prediction of total genetic FRUIT3CAT) financed by the European Regional Development Fund through value using genome-wide dense marker maps. Genetics 157, 1819 (2001). the FEDER frame of Catalonia 2014–2020 and by the CERCA Program from 6. Georges, M., Charlier, C. & Hayes, B. Harnessing genomic information for Generalitat de Catalunya. We acknowledge financial support from the Spanish livestock improvement. Nat. Rev. Genet. 20, 135–156 (2019). Ministry of Economy and Competitiveness through the “Severo Ochoa 7. Laurens, F. et al. An integrated approach for increasing breeding efficiency in Program for Centres of Excellence in R&D” 2016–2019 (SEV‐20150533). This apple and peach in Europe. Hortic. Res. 5, 11 (2018). work has been partially funded by the EU seventh Framework Program, the 8. Kumar, S. et al. Genomic selection for fruit quality traits in apple (Mal- FruitBreedomics project No. 265582: Integrated Approach for Increasing us×domestica Borkh.). PLoS ONE 7, e36674 (2012). Breeding Efficiency in Fruit Tree Crops. The views expressed in this work are 9. Muranty, H. et al. Accuracy and responses of genomic selection on key traits in the sole responsibility of the authors and do not necessarily reflect the views of apple breeding. Hortic. Res. 2, 15060 (2015). the European Commission. Christian Dujak was supported by “DON CARLOS 10. Migicovsky, Z. et al. Genome to phenome mapping in apple using historical ANTONIO LOPEZ” Abroad Postgraduate Scholarship Program, BECAL-Paraguay. data. Plant Genome 9, https://doi.org/10.3835/plantgenome2015.11.0113 (2016). 11. McClure, K. A. et al. A genome-wide association study of apple quality and Author details scab resistance. Plant Genome 11, 170075 (2018). 1 12. McClure, K. A. et al. Genome-wide association studies in apple reveal loci of Molecular Plant Breeding, Institute of Agricultural Sciences, ETH Zurich, 8092 Zurich, Switzerland. 2Breeding Research group, Agroscope 8820 Wädenswil, large effect controlling apple polyphenols. Hortic. Res. 6, 107 (2019). Switzerland. 3IRTA (Institut de Recerca i Tecnologia Agroalimentàries), 08140 13. Roth, M. et al. Genomic prediction of fruit texture and training population Caldes de Montbui, Barcelona, Spain. 4Centre for Research in Agricultural optimization towards the application of genomic selection in apple. Hortic. Genomics (CRAG) CSIC-IRTA-UAB-UB, Campus UAB, 08193 Bellaterra, Barcelona, Res. 7, 148 (2020). Spain. 5Better3fruit N.V., 3202 Rillaar, Belgium. 6Biometris, Wageningen 14. Kumar, S., Hilario, E., Deng, C. H. & Molloy, C. Turbocharging introgression University and Research, 6708 PB Wageningen, The Netherlands. 7Hendrix breeding of perennial fruit crops: a case study on apple. Hortic. Res. 7, 47 (2020). Genetics Research, Technology and Services B.V., PO Box 1145830AC Boxmeer, 15. Covarrubias-Pazaran, G. et al. Multivariate GBLUP improves accuracy of The Netherlands. 8IRHS, Université d’Angers, INRAE, Institut Agro, SFR 4207 genomic selection for yield and fruit weight in biparental populations of QuaSaV, 49071 Beaucouzé, France. 9ARAID (Fundación Aragonesa para la Vaccinium macrocarpon Ait. Front. Plant Sci. 9, 1310 (2018). Investigación y el Desarrollo), 50018 Zaragoza, Spain. 10Research Centre 16. Peace, C. P. et al. Apple whole genome sequences: recent advances and new Laimburg, 39040 Auer, Italy. 11Department of Horticultural Science, University prospects. Hortic. Res. 6, 59 (2019). of Minnesota, St. Paul, MN 55108, USA. 12Institute of Biology and Environmental 17. Bianco, L. et al. Development and validation of a 20K single nucleotide Sciences, University of Oldenburg, 26129 Oldenburg, Germany. 13Research polymorphism (SNP) whole genome genotyping array for apple (Malus × Institute of Horticulture, 96-100 Skierniewice, Poland. 14School of Agriculture, domestica Borkh). PLoS ONE 9, e110377 (2014). Policy and Development, University of Reading, Whiteknights RG6 6AR 18. Bianco, L. et al. Development and validation of the Axiom®Apple480K SNP Reading, UK. 15Laboratory for Plant Genetics and Crop Improvement, KU genotyping array. Plant J. 86, 62–74 (2016). Leuven B-3001 Leuven, Belgium. 16Present address: GAFL, INRAE, 84140 19. Muranty, H. et al. Using whole-genome SNP data to reconstruct a large multi- Montfavet, France generation pedigree in apple germplasm. BMC Plant Biol. 20, 2 (2020). 20. Rodríguez-Álvarez, M. X., Boer, M. P., van Eeuwijk, F. A. & Eilers, P. H. Correcting Data availability for spatial heterogeneity in plant breeding experiments with P-splines. Spat. All SNP genotypic data generated with the 480K array used in this study have Stat. 23, 52–71 (2018). been deposited in the INRAe dataset archive (https://data.inrae.fr/) at https:// 21. Lassois, L. et al. Genetic diversity, population structure, parentage analysis, and doi.org/10.15454/IOPGYF. All SNP genotypic data generated using the 20K construction of core collections in the French apple germplasm based on SSR array used in this study have been deposited in the INRAe dataset archive markers. Plant Mol. Biol. Report. 34, 827–844 (2016). (https://data.inrae.fr/) at https://doi.org/10.15454/1ERHGX. The imputed SNP 22. Kumar, S., Volz, R. K., Alspach, P. A. & Bus, V. G. M. Development of a recurrent marker datasets generated during the current study are available from the apple breeding programme in New Zealand: a synthesis of results, and a corresponding author on reasonable request. The phenotypic data generated proposed revised breeding strategy. Euphytica 173, 207–222 (2010). during the study (phenotypic least-square means of genotypes across 23. McCarthy, S. et al. A reference panel of 64,976 haplotypes for genotype environments) are available in the Supplementary Table 6. imputation. Nat. Genet. 48, 1279–1283 (2016). 24. Pégard, M. et al. Sequence imputation from low density single nucleotide polymorphism panel in a black poplar breeding population. BMC Genom. 20, Conflict of interest 302 (2019). The authors declare that they have no conflict of interest. 25. Vanderzande, S., Micheletti, D., Troggio, M., Davey, M. W. & Keulemans, J. Genetic diversity, population structure, and linkage disequilibrium of elite and Supplementary Information accompanies this paper at (https://doi.org/ local apple accessions from Belgium using the IRSC array. Tree Genet. Genom. 10.1038/s41438-020-00408-8). 13, 125 (2017). 26. Daccord, N. et al. High-quality de novo assembly of the apple genome and Received: 23 April 2020 Revised: 25 August 2020 Accepted: 6 September methylome dynamics of early fruit development. Nat. Genet. 49, 1099–1106 2020 (2017). 27. Cornille, A. et al. New insight into the history of domesticated apple: sec- ondary contribution of the European wild apple to the genome of cultivated varieties. PLoS Genet. 8, e1002703 (2012). 28. Segura, V. et al. An efficient multi-locus mixed-model approach for genome- References wide association studies in structured populations. Nat. Genet. 44, 825–830 1. Food and Agriculture Organization of the United Nations. FAOSTAT (2012). (FAO, 2017). Access 19 Mar 2019. http://www.fao.org/faostat/en/#home. 29. Celton, J.-M. et al. Deciphering the genetic determinism of bud phenology 2. Urrestarazu, J. et al. Analysis of the genetic diversity and structure across a in apple progenies: a new insight into chilling and heat requirement effects wide range of germplasm reveals prominent gene flow in apple at the on flowering dates and positional candidate. New Phytol. 192, 378–392 European level. BMC Plant Biol. 16, 130 (2016). (2011). 3. Urrestarazu, J. et al. Genome-wide association mapping of 30. Allard, A. et al. Detecting QTLs and putative candidate genes involved in flowering and ripening periods in apple. Front. Plant Sci. 8, 1923 budbreak and flowering time in an apple multiparental population. J. Exp. Bot. (2017). 67, 2875–2888 (2016). Jung et al. Horticulture Research (2020)7:189 Page 16 of 16 31. Chagné, D. et al. Genetic and environmental control of fruit maturation, dry 42. Browning, S. R. & Browning, B. L. Rapid and accurate haplotype phasing and matter and firmness in apple (Malus × domestica Borkh.). Hortic. Res. 1, 14046 missing-data inference for whole-genome association studies by use of (2014). localized haplotype clustering. Am. J. Hum. Genet. 81, 1084–1097 (2007). 32. Larsen, B. et al. Genome-wide association studies in apple reveal loci for aroma 43. Clayton, D. snpStats: SnpMatrix and XSnpMatrix Classes and Methods (2017). volatiles, sugar composition, and harvest date. Plant Genome 12, 180104 (2019). 44. Paradis, E., Claude, J. & Strimmer, K. APE: analyses of phylogenetics and evo- 33. Gezan, S. A., Osorio, L. F., Verma, S. & Whitaker, V. M. An experimental validation lution in R language. Bioinformatics 20, 289–290 (2004). of genomic selection in octoploid strawberry. Hortic. Res. 4, 16070, https://doi. 45. Lê, S., Josse, J. & Husson, F. FactoMineR: an R package for multivariate analysis. org/10.1038/hortres.2016.70 (2017). J. Stat. Softw. 25, 1–18 (2008). 34. Hardner, C. M. et al. Prediction of genetic value for sweet cherry fruit maturity 46. Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of among environments using a 6K SNP array. Hortic. Res. 6, 6 (2019). ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009). 35. Schulthess, A. W. et al. Multiple-trait- and selection indices-genomic predic- 47. International Union for the Protection of New Varieties of Plants (UPOV). tions for grain yield and protein content in rye for feeding purposes. Theor. Guidelines for the conduct of tests for distinctness, uniformity and stability: TG/14/ Appl. Genet. 129, 273–287 (2016). 9 Apple (fruit varieties) (International Union for the Protection of New Varieties 36. Lado, B. et al. Resource allocation optimization with multi-trait genomic pre- of Plants (UPOV), Geneva, 2005). diction for bread wheat (Triticum aestivum L.) baking quality. Theor. Appl. Genet. 48. Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models 131, 2719–2731 (2018). using lme4. J. Stat. Softw. 67, 1–48 (2015). 37. Fernández-Fernández, F. SID 5 (Research Project Final Report) (Fingerprinting 49. Højsgaard, S. & Halekoh, U. doBy: Groupwise statistics, LSmeans, Linear Contrasts, the National Apple & Pear Collections, Defra, 2010). Utilities (2018). 38. Bink, M. C. A. M. et al. Bayesian QTL analyses using pedigreed families of an 50. Tang, Y. et al. GAPIT version 2: an enhanced integrated tool for genomic outcrossing species, with application to fruit firmness in apple. Theor. Appl. association and prediction. Plant Genome 9, https://doi.org/10.3835/ Genet. 127, 1073–1090 (2014). plantgenome2015.11.0120 (2016). 39. Jansen, J. & van Hintum, T. Genetic distance sampling: a novel sampling 51. Endelman, J. B. Ridge regression and other kernels for genomic selection with method for obtaining core collections using genetic distances with an R package rrBLUP. Plant Genome 4, 250–255 (2011). application to cultivated lettuce. Theor. Appl. Genet. 114, 421–428 (2007). 52. Bermingham, M. L. et al. Application of high-dimensional feature selection: 40. European Environment Agency. Biogeographical Regions in Europe https:// evaluation for genomic prediction in man. Sci. Rep. 5, 10312 (2015). www.eea.europa.eu/data-and-maps/figures/biogeographical-regions-in- 53. Elsen, J.-M. An analytical framework to derive the expected precision of europe-2 (2017). genomic selection. Genet. Sel. Evol. 49, 95 (2017). 41. Vanderzande, S. et al. High-quality, genome-wide SNP genotypic data for 54. R Core Team. R: A Language And Environment For Statistical Computing pedigreed germplasm of the diploid outbreeding species apple, peach, and (R Foundation for Statistical Computing, Vienna, Austria, 2014). sweet cherry through a common workflow. PLoS ONE 14, e0210928 (2019). 55. Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).