# Prediction of genetic value for sweet cherry fruit maturity among environments using a 6K SNP array

Craig M. Hardner, Ben J. Hayes, Satish Kumar, Stijn Vanderzande, Lichun Cai, Julia Piaskowski, José Quero-Garcia, José Antonio Campoy-Corbalan, Teresa Barreneche, Daniela Giovannini, et al.

► To cite this version:

Craig M. Hardner, Ben J. Hayes, Satish Kumar, Stijn Vanderzande, Lichun Cai, et al. Prediction of genetic value for sweet cherry fruit maturity among environments using a 6K SNP array. Horticulture research, 2019, 6 (1), 10.1038/s41438-018-0081-7. hal-02620841

HAL Id: hal-02620841

https://hal.inrae.fr/hal-02620841

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Prediction of genetic value for sweet cherry fruit maturity among environments using a 6K SNP array

Craig M. Hardner<sup>1</sup>, Ben J. Hayes<sup>1</sup>, Satish Kumar<sup>2</sup>, Stijn Vanderzande<sup>3</sup>, Lichun Cai<sup>4</sup>, Julia Piaskowski<sup>3</sup>, José Quero-Garcia<sup>5</sup>, José Antonio Campoy<sup>5</sup>, Teresa Barreneche<sup>5</sup>, Daniela Giovannini<sup>6</sup>, Alessandro Liverani<sup>6</sup>, Gérard Charlot<sup>7</sup>, Miguel Villamil-Castro<sup>1</sup>, Nnadozie Oraguzie<sup>8</sup> and Cameron P. Peace<sup>3</sup>

## Abstract

The timing of fruit maturity is an important trait in sweet cherry production and breeding. Phenotypic variation for phenology of fruit maturity in sweet cherry appears to be under strong genetic control, but that control might be complicated by phenotypic instability across environments. Although such genotype-by-environment interaction (G × E) is a common phenomenon in crop plants, knowledge about it is lacking for fruit maturity timing and other sweet cherry traits. In this study, 1673 genome-wide SNP markers were used to estimate genomic relationships among 597 weakly pedigree-connected individuals evaluated over two seasons at three locations in Europe and one location in the USA, thus sampling eight 'environments'. The combined dataset enabled a single meta-analysis to investigate the environmental stability of genomic predictions. Linkage disequilibrium among marker loci declined rapidly with physical distance, and ordination of the relationship matrix suggested no strong structure among germplasm. The most parsimonious G × E model allowed heterogeneous genetic variance and pairwise covariances among environments. Narrow-sense genomic heritability was very high (0.60–0.83), as was accuracy of predicted breeding values (>0.62). Average correlation of additive effects among environments was high (0.96) and breeding values were highly correlated across locations. Results indicated that genomic models can be used in cherry to accurately predict date of fruit maturity for untested individuals in new environments. Limited G × E for this trait indicated that phenotypes of individuals will be stable across similar environments. Equivalent analyses for other sweet cherry traits, for which multiple years of data are commonly available among breeders and cultivar testers, would be informative for predicting performance of elite selections and cultivars in new environments.

## Introduction

The timing of fruit maturity is important for sweet cherry (Prunus avium L.) production, particularly given the brief shelf life of the fruit. Fruit of an individual tree tends to ripen over a narrow window of 3–7 days and firmness and taste decline rapidly over 4–7 days after

harvest, even under cold storage<sup>1–4</sup>. Similarly, the harvest window for cherry is relatively short<sup>5–7</sup>. Fruit produced in the early or later part of the production season attracts higher prices than in the mid-season<sup>8</sup>.

There is evidence that phenotypic variation in phenology of fruit maturity in sweet cherry is under strong genetic control, suggesting breeding opportunities for this trait. Traditionally, sweet cherry cultivars are divided according to their maturity date into early (e.g., 'Burlat' and 'Chelan'), mid- (e.g., 'Bing') and late season (e.g., 'Sumtare')<sup>9</sup>. Individual broad-sense heritability of fruit maturity timing is high (0.76–0.83)<sup>10,11</sup> and large-effect Quantitative Trait

<figure><img src="image_2.png" /></figure>

<footnote>Correspondence: Craig M. Hardner (craig.hardner@uq.edu.au)</footnote>

<footnote n="1">¹University of Queensland, Queensland Alliance for Agriculture and Food Innovation, Brisbane, QLD 4072, Australia</footnote>

<footnote n="2">²The New Zealand Institute for Plant and Food Research Limited, Hawke's Bay Research Centre, Hastings 4130, New Zealand</footnote>

<footnote>Full list of author information is available at the end of the article.</footnote>
Loci (QTLs) on linkage groups 1, 4, 5, and 6 have been detected within a full-sib family ('Regina' × 'Lapins'), each explaining up to 20% of phenotypic variation¹⁰,¹².

While cherry has been cultivated for more than 2000 years, breeding commenced around the early 1800s¹³,¹⁴. From the early 20th century, modern breeding involving selection of elite parents was undertaken in the USA and Canada, followed by European countries such as the UK, Russia, Ukraine, and others by the 1950s. In the last 30 years, breeding has progressed rapidly, fuelled in part by the adoption of dwarfing rootstocks allowing intensive orchard production¹⁵. Sansavini and Lugli¹⁶ noted that during the period 1991–2004, 230 new cultivars were released, this number being second in stone fruit only to peach. Nevertheless, some very old selections, such as 'Burlat' or 'Bing', and cultivars with unclear origin such as 'Ambrunés', 'Emperor Francis', 'Germersdorfer', 'Hedelfingen', 'Napoleon' or '0900 Ziraat', are still widely planted. Cherry breeding has relied on a narrow genetic base for the main founders¹⁷, and modern cultivars are only a few generations removed from early ancestors¹⁸.

The development of early- and late-maturing cultivars is a major objective in sweet cherry breeding across the globe¹⁰,¹³. Replacing 'Burlat' has been an important target for decades, although most new cultivars of the same early maturity period have undesirable fruit quality in terms of firmness, tolerance to rain-induced fruit cracking, or post-harvest shelf life. The breeding program conducted at Summerland, Canada, has released several commercially successful late-maturing cultivars, such as 'Sumtare' (Sweetheart™), '13S2009' (Staccato™) and '13S2101' (Sovereign™).

Genotype-by-environment interaction (G × E) is a common phenomenon in crop plants that can complicate selection¹⁹,²⁰. Confirming the degree of G × E is an important initial step in developing strategies for commercial deployment of improved germplasm in target production environments¹⁹,²¹. If the degree of G × E is small, germplasm can be selected for overall mean performance across the target region. Alternatively, if the magnitude of G × E is large and repeatable factors can be identified that explain some of the pattern in G × E, higher gains can be achieved by using those factors to subdivide production regions into 'mega-environments' and select germplasm targeted for each. Otherwise, the presence of G × E represents experimental noise, reducing selection accuracy.

Detection of G × E in plants requires replication of genetic effects across multi-environment trials (METs). Commonly, METs in horticulture are established with individuals that are clonally replicated across environments, e.g., refs. ²²⁻²⁵, or are within related families connected through a pedigree structure that is used to estimate an expected relationship matrix²⁶. Thereby, the

performance of an individual in one environment can be genetically correlated with its performance in other environments²⁷,²⁸ and that correlation can be used to predict the performance of genetic treatments (e.g., cultivars, or even specific alleles) in environments in which they have not been directly tested²⁹. However, conducting conventional METs in horticultural tree crops is expensive, particularly due to the large size of experimental units and long juvenility periods. As a consequence, tree fruit breeding and elite selection or cultivar evaluation programs are often localised and have limited replication of germplasm among programs. These features limit the opportunity to evaluate G × E for horticultural tree crops.

Linear mixed models are the preferred approach for analysing data from METs³⁰. Such models support estimation of variance and covariances to describe G × E, incorporation of complex models that might include different experimental designs among trials or spatial effects, and prediction of unbiased genetic values from unbalanced data. A common model used to describe G × E from METs is a simple main effect and homogeneous interaction (G + G × E) linear mixed model. However, this model is restrictive as it unrealistically assumes a common genetic variance across environments and a common covariance (and hence correlation) among all pairs of environments³⁰. More general models allow greater flexibility, although a highly unstructured covariance matrix might be difficult to estimate, particularly when the number of dimensions is large and the data for each dimension are limited³¹⁻³³. The interaction term in the simple G + G × E model is an average of the individual pairwise covariances in the general model³⁴ A parsimonious alternative is the FA parameterisation of the genetic-by-environment covariance matrix, which estimates environmental loadings for a reduced number of hypothetical orthogonal dimensions that maximise differences among environments (λ), and estimates specific variances for each environment for genetic effects not accounted for by the loadings (ψ)³¹.

Genetic relationships among individuals can be described using a DNA marker-based genomic relationship matrix (GRM)³⁵, to approximate identity-by-descent as replication of chromosome segments among individuals deployed across METs. An advantage of the GRM is that it can describe heterogeneity of realised relationships within families that occur as a consequence of Mendelian sampling³⁶. A realised relationship matrix is expected to be more accurate than a pedigree-derived relationship matrix, because the latter assumes a mean expected relationship among relatives. In addition, realised relationships might capture cryptic or undocumented relationships³⁷,³⁸. The GRM could therefore be used to replace the pedigree-derived relationship matrix to study G × E. Considerations in using the GRM with a factor
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Location of field trials included in this study</figcaption></figure>

analytic (FA) parameterisation of the genetic-by-environment covariance matrix to model G × E have been described previously³⁹.

The aim of this study was to use a G × E model involving a GRM to combine phenotypic data on fruit maturity timing collected in multiple locations for cherry germplasm having little clonal replication among locations, to thereby estimate patterns in G × E and predict performance of individuals in environments in which they have not been tested. By leveraging the information from existing breeding program datasets in this way, we hypothesised achieving improved accuracy of genetic predictions within and across locations compared to predictions arising from datasets of single breeding programs.

# Methods

## Description of data

Germplasm used consisted of 597 sweet cherry cultivars, selections and unselected offspring, subsets of which were grown at four locations: Prosser, WA, USA (46.291383, -119.746753); Bourran (44.333793, 0.413504) and Balandran (43.757312, 4.461919), France; and Forlì (44.216667, 12.05000), Italy (Fig. 1). Fruit maturity timing was assessed at each location for two seasons (Table 1) to give eight unique location-by-season environments. Most individuals at Balandran were also assessed at Bourran, but there was less overlap among the individuals in France and those at the Italian and Prosser locations, and only five individuals were present at all four locations (Fig. 2).

The 71 individuals assessed at Balandran, France (trial 'Balandrin.FR'), belonged to sweet cherry trials of the 'Chart of evaluation of fruit varieties' program coordinated by Centre Technique Interprofessionnel des Fruits et Légumes (CTIFL) and involving Institut National de la Recherche Agronomique (INRA), regional experimental stations, representatives from growers associations, and nurserymen and private breeders from France and beyond. Balandran is one of the three locations in France where new cultivars are initially evaluated and is 44 m above sea level. The climate is classified as Mediterranean and the trial was planted on Costiéères soil. Each cultivar was grafted onto rootstock 'Colt' and planted as two clonal replicates, adjacent wherever possible. Trees were trained to a vase with four to five permanent branches. The trial was irrigated and managed using standard approaches. Fruit maturity timing was measured as the date (in Julian days) on which 50% of the fruit on both accessions were mature.

The 193 individuals assessed at Bourran, France (trial 'Bourran.FR') were accessions of a germplasm collection maintained by the Prunus Genetic Resources Center of INRA at Bourran (Lot & Garonne), near Bordeaux⁴⁰. The site is 50 m above sea level. The climate is classified as oceanic and the soil type is clay alluvium. Two replicate trees for each individual on the rootstock 'Maxma 14' were planted adjacent in a random design. Approximately half were of French breeding origin and the remainder were from 15 other countries in North America, Asia and Europe. Accessions were historic landrace cultivars (n = 85) and recently bred cultivars
<table><caption>Table 1 Summary of the phenotypic dataset of fruit maturity timing from eight location-by-season environments</caption><thead><tr><th>Trial</th><th>Location</th><th>Country</th><th>Trait</th><th>Season</th><th>No. of individuals</th><th>Mean<sup>a</sup></th><th>Variance<sup>a</sup></th></tr></thead><tbody><tr><td>Balandran.FR</td><td>Balandran</td><td>France</td><td>Fruit maturity date</td><td>1997</td><td>61</td><td>144</td><td>103.8</td></tr><tr><td>Balandran.FR</td><td>Balandran</td><td>France</td><td>Fruit maturity date</td><td>1998</td><td>50</td><td>148</td><td>100.5</td></tr><tr><td>Bourran.FR</td><td>Bourran</td><td>France</td><td>Fruit maturity date</td><td>2014</td><td>187</td><td>150</td><td>81.3</td></tr><tr><td>Bourran.FR</td><td>Bourran</td><td>France</td><td>Fruit maturity date</td><td>2015</td><td>192</td><td>150</td><td>74.6</td></tr><tr><td>Forlí.IT</td><td>Forlí</td><td>Italy</td><td>Date of first harvest</td><td>2014</td><td>50</td><td>149</td><td>141.0</td></tr><tr><td>Forlí.IT</td><td>Forlí</td><td>Italy</td><td>Date of first harvest</td><td>2015</td><td>55</td><td>152</td><td>122.5</td></tr><tr><td>Prosser.US</td><td>Prosser</td><td>USA</td><td>Fruit maturity date</td><td>2011</td><td>231</td><td>186</td><td>60.1</td></tr><tr><td>Prosser.US</td><td>Prosser</td><td>USA</td><td>Fruit maturity date</td><td>2012</td><td>360</td><td>184</td><td>50.6</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Venn diagram of numbers of individuals specific to and shared among four trial locations. FR France, IT Italy, US USA</figcaption></figure>

currently planted in Europe (n = 106). Two accessions of each of the cultivars Noir d'Ecully and Giorgia were also included, each corresponding to a different period of introduction into the INRA collection. Trees were trained to an open vase and were irrigated and managed using a standard commercial regime. Fruit maturity timing was measured as the date (in Julian days) on which at least 70% of the fruit on both replicate trees were mature, using firmness, taste and skin colour indicators.

The pool of 56 individuals assessed at Forlí, Italy (trial 'Forlí.IT'), was composed of traditional germplasm from the Emilia-Romagna territory, modern cultivars from various breeding programs, breeding selections, and a few clones of known cultivars selected in the 1970s after irradiation of bud sticks and are part of the sweet cherry collection at CREA (Council for Agricultural Research and Economics)—Forlí. The site is 34 m above sea level, with Mediterranean climate and silty-clay soil with medium-low organic matter soil. Each individual was represented by two to three adjacent replicate trees graf- ted onto the rootstock 'Colt'. Trees were managed according to integrated production practices, irrigated and trained to a vase. Fruit maturity timing was measured as the date (in Julian days) that the first 10% of fruit on all

replicate trees were ripe for consumption⁴¹ using firmness, taste and skin colour indicators.

The 384 individuals evaluated at Prosser, USA (trial 'Prosser.US') were within the RosBREED sweet cherry Crop Reference Set (n = 268)⁴² and the Pacific Northwest Sweet Cherry Breeding Program at Washington State University. The site is ~200 m above sea level. The climate is classified as semi-arid and the soil is well-drained non-saline from the Aridisols order. The RosBREED Crop Reference Set consists of cultivars, founders, ancestors, advanced selections and unselected offspring representing germplasm used in US breeding programs⁴²,⁴³. A total of 53 Prosser.US individuals were cultivars or breeding selections, with the remainder consisting of unselected offspring from 98 families (53 families with two or more offspring). The germplasm was planted between 2006 and 2008 as single trees of each individual, with offspring often arranged in families. Trees were managed using conventional orchard management practices. All indivi- duals were grown on 'Gisela 6' rootstock except for own-rooted offspring. Fruit maturity timing was measured in 2011 and 2012 within the RosBREED project⁴⁴ as the date (in Julian days) at which 50% of fruit on a tree were mature using firmness, taste and skin colour indicators⁴³.

SNP data generated using the RosBREED cherry 6K SNP array v1⁴⁵ were available from other projects for the 597 individuals from the four locations described above (Fig. 2). Briefly, DNA samples had been extracted using standard methods, e.g.⁴⁰, and run on the array following manufacturer's instructions. Resultant probe intensity data were used as input for GenomeStudio to obtain genotype calls for each SNP of each individual⁴⁶. For SNP data available from RosBREED⁴⁵, CREA and CITA (Centro de Investigación y Tecnología Agroalimentaria de Aragón, Spain), GenomeStudio results were exported and used as input to the software program ASSIsT⁴⁷ to obtain a subset of robust SNPs with reliable genotype calls. A multi-step pipeline using FlexQTL, VisualFlexQTL⁴⁸ and

<footnote>ªUnits are in Julian days</footnote>
<table><caption>**Table 2 Standardised field terms for use in combining phenotypic datasets**</caption><thead><tr><td>Field</td><td>Description</td></tr></thead><tbody><tr><td>Location</td><td>Spatial unit containing one or more trials for which genetic effects are assumed to be homogeneous (i.e., homogeneous genetic variance and genetic correlation of 1)</td></tr><tr><td>Trial</td><td>Spatial unit that contains experimental units grouped by a particular factor (usually planting date)</td></tr><tr><td>Section</td><td>Spatial unit within location for which experimental units are contiguous</td></tr><tr><td>Row</td><td>Spatial coordinate within section that indexes the distribution of experimental unit along the planting row dimension</td></tr><tr><td>Position</td><td>Spatial coordinate within section that indexes experimental unit along dimension perpendicular to planting row</td></tr><tr><td>Block</td><td>Design unit within section that indexes group of experimental units grouped to account for spatial variation</td></tr><tr><td>Plot</td><td>Design unit within block that indexes group of experimental units grouped to account for spatial variation</td></tr><tr><td>Plant</td><td>Individual unit that represents unique propagation event</td></tr><tr><td>Planting date</td><td>Date at which experimental unit was established</td></tr><tr><td>REF_ID</td><td>Standardised genetic treatment</td></tr><tr><td>Experimental unit</td><td>Unit at which genetic treatment is applied, i.e., unit on which sampling may be temporally repeated. Created as the unit combination of design units and genetic treatment</td></tr><tr><td>Year</td><td>Year at which observation is made</td></tr><tr><td>Observation unit</td><td>Created as unit combination of experimental unit and factor defining repeated observations on experimental unit</td></tr><tr><td>y</td><td>Untransformed observation</td></tr></tbody></table>

Pedihaplotyper⁴⁹ to detect genotyping errors was applied to obtain a high-quality genotypic dataset of 1617 SNPs. Data for 1215 high-quality SNPs loci were available from INRA, and for 1636 SNPs from CREA.

A data-curation pipeline was developed to manage germplasm, SNP and phenotypic data contributed from the various sources. The first step in the process was standardisation of names of individuals common across datasets. Initially, the format of an individual's name was simplified by removing special characters, changing all upper case letters to lower case, and inserting an underscore between characters separated by a space. Formatted names were then compared to an existing dictionary of original names, formatted names and final names (which were the formatted final names used to link data in cases of multiple sources of data for the same individual) created from previous data sources. Full and partial matching were used to identify potential matches in the existing dictionary. New names that could not be reconciled were added to the dictionary.

Curation of SNP data involved standardisation of SNP names across the datasets against NCBI standard identities. SNP data from individuals common to multiple locations were compared to ensure consistency and fill in missing genotypic calls. The very few inconsistencies in genotypes among common individuals were set to missing, except in one case: 'Noire de Meched' genotyped by INRA and found to be substantially different to an individual with the same name in the CREA collection (the

former believed to be the original, true-to-type cultivar introduced from Iran). In this case, cultivar names were modified to distinguish the difference. SNP loci with missing data of >30% or a minor allele frequency <0.05 were removed, leaving 1273 loci for 550 individuals with unique genome-wide DNA profiles (i.e., 47 of the original 597 individuals were determined to be duplicate genotypes). The missing 33,704 (4% of the total) individual SNP locus genotypes were imputed with a hidden Markov Model implemented in Beagle⁵⁰,⁵¹ in the R package synbreed⁵² using only flanking markers, i.e., without considering family information.

Phenotypic data curation began with incorporation of standardised individual names. Next, field names in original data files that identified experimental and sampling design factors were identified and renamed using standardised terms (Table 2). A field representing the experimental unit in each trial was created to identify the smallest unit to which a unique random genetic treatment was applied⁵³. Fields not in the list of standard design and sampling effects were retained. Phenotypic data was reshaped to long format (i.e., each record represented an individual observation).

## Statistical methods

A multivariate genomic prediction linear mixed model was implemented to estimate model parameters and predict fixed and random effects for fruit maturity data collected from individuals planted across the four
<table><caption>**Table 3** Estimates of additive genomic variation (vG), residual variation (vR), phenotypic variation (vP) and narrow-sense heritability ($h^2$) for eight location-by-season environments</caption><thead><tr><th>Location</th><th>Season</th><th>vG</th><th>vR</th><th>vP</th><th>h<sup>2</sup></th></tr></thead><tbody><tr><td>Balandran.FR</td><td>1997</td><td>111.3</td><td>13.8</td><td>125.1</td><td>0.89</td></tr><tr><td>Balandran.FR</td><td>1999</td><td>101.4</td><td>10.5</td><td>111.9</td><td>0.91</td></tr><tr><td>Bourran.FR</td><td>2014</td><td>60.3</td><td>12.7</td><td>73.0</td><td>0.83</td></tr><tr><td>Bourran.FR</td><td>2015</td><td>57.1</td><td>8.1</td><td>65.2</td><td>0.88</td></tr><tr><td>Forlí.IT</td><td>2014</td><td>110.4</td><td>19.3</td><td>129.7</td><td>0.85</td></tr><tr><td>Forlí.IT</td><td>2015</td><td>98.4</td><td>20.0</td><td>118.4</td><td>0.83</td></tr><tr><td>Prosser.US</td><td>2011</td><td>27.6</td><td>18.1</td><td>45.7</td><td>0.60</td></tr><tr><td>Prosser.US</td><td>2012</td><td>26.5</td><td>17.6</td><td>44.1</td><td>0.60</td></tr></tbody></table>

locations and assessed over Table 3 two seasons, i.e., in eight location-by-season environments. The general mixed model for predicting genomic breeding values of $n_G$ individuals, evaluated at $n_L$ locations and $n_{Sj}$ seasons within the $j^{th}$ location, was:

$$ \mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{Z}\mathbf{u} + \mathbf{r}, $$

where $\mathbf{y}$ was a vector of observations, $\mathbf{b}$ was vector of unknown fixed effects (mean effect of location-by-season environment, and year of planting at the Balandran location—the only location where this factor was relevant), $\mathbf{X}$ was the design matrix for fixed effects, $\mathbf{u}$ was a vector of unknown random effects (genomic breeding value for each of the $n_G$ individuals, at each of the $n_E$ location-by-season environments), $\mathbf{Z}$ was the design matrix for random effects and $\mathbf{r}$ was a vector of unknown residual effects, with variation of observations assumed to follow a multivariate normal distribution:

$$ \operatorname{var}(\mathbf{y}) = \mathbf{Z}\mathbf{G}\mathbf{Z}^T + \mathbf{R}, $$

where $\mathbf{G}$ was the variance-covariance of additive genomic effects for the $n_G$ individuals among the $n_E$ location-by-season environments and $\mathbf{R}$ was the variance-covariance among residual effects.

The $(n_G \times n_E) \times (n_G \times n_E)$ matrix $\mathbf{G}$ was separated into the $n_G \times n_G$ additive genomic relationship matrix ($\mathbf{G}_A$), estimated following VanRaden³⁵ method 1 using the complete imputed individual-by-SNP matrix and the $n_E \times n_E$ additive genetic-by-environment covariance matrix ($\mathbf{G}_{A \times E}$), i.e.:

$$ \mathbf{G} = \mathbf{G}_A \otimes \mathbf{G}_{A \times E}, $$

where $\otimes$ was the Kronecker product matrix operation. The $n_E \times n_E$ additive genomic-by-environment covariance

matrix, $\mathbf{G}_{A \times E}$, can be modelled using various structures³⁰. Assuming homogeneous genomic variances across all environments and perfect genomic correlations (i.e., $r_A = 1$) among all pairwise environmental combinations leads to $\mathbf{G}_{A \times E} = \sigma_A^2 \otimes \mathbf{1}_{n_E}$, where $\mathbf{1}_{n_E}$ was a $n_E \times n_E$ matrix of 1 s. A compound structure, i.e., homogeneous variances on the diagonal and homogenous covariances on the off-diagonals, is equivalent to a main effect interaction model (i.e., $A + A \times E)^{30}$. Alternatively, a FA parameterisation models $\mathbf{G}_{A \times E}$ as a reduced set of parameters that are environmental loadings ($\Lambda_{A \times E}$) on hypothetical orthogonal environmental dimensions that explain the maximum variation at each observed environmental dimension and environment-specific variation not explained by the loadings ($\Psi_{A \times E}$). In this study, only a first-order FA model was examined. The variance-covariance matrix of residual effects ($\mathbf{R}$) was modelled as independent blocks with each block representing the variance-covariance among observations across the seasons at the $j^{th}$ location ($\mathbf{R}_j$), which was separated into $\mathbf{R}_{Sj}$ (residual-by-season covariance matrix for the $j^{th}$ location) and $\mathbf{I}_{nj}$ (the correlation matrix of residuals, an identity matrix of dimensions $n_j$, the number of experimental units evaluated at the $j^{th}$ location).

Parameters of the random and residual models were estimated using average information approaches implemented in the statistical software ASReml⁵⁴. Independent analyses for each specific location-by-season environment were undertaken to review the distribution of residuals. Best linear unbiased estimates (E-BLUEs) of fixed effects and best linear unbiased predictions (E-BLUPs) of random effects (breeding value) were produced following Henderson's⁵⁵ mixed model equations using estimated $\mathbf{G}$ and $\mathbf{R}$ matrices.

The pattern of linkage disequilibrium across the experimental population was examined by plotting individual pairwise squared correlations of SNP locus genotypes among individuals ($\rho^2$) estimated using the software PLINK⁵⁶ against pairwise physical distances. Average linkage disequilibrium (LD) decay in physical distance was converted to genetic distance using the map presented in Guajardo et al.⁵⁷

An eigen-value decomposition of the individual-by-location GRM ($\mathbf{C}$) was undertaken to examine the germplasm genetic structure across the four locations. The matrix $\mathbf{C}$ was estimated as:

$$ \mathbf{M}\mathbf{G}_A\mathbf{M}^T $$

where $\mathbf{M}$ was an $m \times n_G$ incidence matrix (where $m$ was the number of unique individual-by-location realisations and $n_G$ was previously defined) that mapped the incidence of an individual onto specific individual-by-environment realisations. A pedigree-derived additive genetic
relationship matrix was estimated from known pedigree records to examine the similarity between expected (pedigree-derived) and estimated realised GRMs.

The average additive genomic correlation among the eight location-by-season environments was estimated as:

$$ \bar{r}_{A \times E} = \frac{\sigma_A^2}{\sigma_A^2 + \sigma_{A \times E}^2}, $$

where $\hat{\sigma}_A^2$ was the off-diagonals of the compound symmetry $\mathbf{G}_{A \times E}$ structure, or the variance of the additive genetic main effect across environments for the $A + A \times E$ model, and $\hat{\sigma}_{A \times E}^2$ was the difference between the diagonal and the off-diagonals, or the average additive genomic-by-environment variance.

The additive genomic-by-environment covariance matrix of the FA model was estimated as:

$$ \hat{\mathbf{G}}_{A \times E} = \hat{\mathbf{\Lambda}}_{A \times E} \hat{\mathbf{\Lambda}}_{A \times E}^{T} + \hat{\mathbf{\Psi}}_{A \times E}, $$

where $\hat{\mathbf{\Lambda}}_{A \times E}$ was the $n_E \times k$ matrix of estimated loadings for the hypothetical orthogonal factor for the additive genomic variance at the $n_E$ environments, and $\hat{\mathbf{\Psi}}_A$ is a $n_E \times n_E$ diagonal matrix of environment-specific variances.

Individual narrow-sense genomic estimated heritability for the $k^{th}$ season at the $j^{th}$ location was estimated as:

$$ \hat{h}_{jk}^2 = \frac{\hat{\sigma}_{A_{jk}}^2}{\hat{\sigma}_{A_{jk}}^2 + \hat{\sigma}_{R_{jk}}^2}, $$

where $\hat{\sigma}_{A_{jk}}^2$ and $\hat{\sigma}_{R_{jk}}^2$ were the estimated additive genomic and residual variances for each of the $n_E$ location-by-season environments. Genomic breeding values for fruit maturity timing were predicted using estimated parameters from the FA model.

Two cross-validation scenarios were used to estimate prediction accuracy (PACC). Firstly, a within-location validation population was created by subsampling individuals within locations, as the main interest of this study was performance predictions of individuals in environments in which they have not been tested, but phenotypes of other genotyped individuals have been used to characterise the environment of interest. Five-fold validation within locations was used. Adjusted phenotypes were obtained by removing from the phenotypic observations the effects of location, season, location-by-season, and year of planting at the Balandran location estimated from the fit of the FA model to the full dataset. Individuals at each location were then randomly assigned to one of the five sets to define location-by-set validation populations. For each validation analysis, adjusted phenotypes for the individuals in the particular location-by-set validation population were set to missing and parameters for the FA $A \times E$ model were re-estimated using the training

population and then used to predict additive genomic values for the validation population. Predictive ability (PA) was estimated as the correlation between the genomic predicted breeding values for the particular location-by-set validation population and the original adjusted phenotypes. Prediction accuracy was estimated for the particular location-by-set validation population by dividing PA by the square-root of the individual narrow-sense genomic estimated heritability of genomic values from the fit of the FA $A \times E$ model to the full dataset, as independent heritability estimates were unavailable. Standard error of PACC was estimated as $1/\sqrt{n_{G \times Lk}}$, where $n_{G \times Lk}$ was the number of individuals in the location-by-set validation population$^{58}$.

A fourfold cross-validation was undertaken for a second scenario of dropping out the complete phenotypic data from one entire location for each validation iteration. Predictive ability was estimated as the correlation among the adjusted phenotypic performance of individuals in the removed location and the mean of the predicted breeding values for these individuals across the three preserved locations. Prediction accuracy was again estimated using genomic narrow-sense heritability of the removed environments estimated from the fit to the full dataset.

## Results

### Structure of germplasm

The average LD between SNPs within a 100-kb window was 0.36 with the value for 20% of pairs being >0.80 (Fig. 3). There was a rapid decay in LD with distance within 250–500 kb with observed $\rho^2$-value <0.2 between 140 and 240 kb, which corresponds to ~0.30 and 0.52 cM. At 1 cM/465 kb, LD declined to <0.14.

The first two vectors of an eigen decomposition of the estimated individual-by-location GRM explained 12% of the variation in the matrix (Fig. 4). No obvious structure was apparent in the overall relationships among and within germplasm assessed at the four locations. Contributions of elements of the relationship matrix estimated from pedigree information were quite different

<figure><img src="image_5.png" /><figcaption>Fig. 3 Linkage disequilibrium (LD) by physical distance (kb) for individual loci (grey dots) and averaged for each set of 1000 SNP pairs ranked according to the distance between SNPs (black dots)</figcaption></figure>
from the distributions for realised relationships estimated using SNP data (Fig. 5). While the mode for the distribution of diagonal elements for both the pedigree- and SNP-derived relationship coefficients was one, and 0 for the off-diagonals, the distributions of pedigree-derived coefficients were highly skewed towards these values, in

contrast to the more balanced distribution of GRM coefficients (both diagonal and off-diagonal) around the mode. The highest value for the pedigree-derived diagonals was 1.5, compared to 1.8 for the estimated realised relationship coefficients. The variance of diagonal elements was greater for the GRM matrix (0.0228) compared to the pedigree-derived relationship matrix (0.0024), as was the case for the off-diagonals (0.0160 compared to 0.0103).

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Plot of the first two vectors from an eigen decomposition of the individual-by-location genomic relationship matrix by testing location</figcaption></figure>

## Model fit

The compound symmetry $A + A \times E$ model (log-likelihood = -2554.2) was a significantly ($p < 0.001$) better fit to the data than a model that did not account for $A \times E$ (logl = -2560.7). Estimated average $\bar{r}_{A \times E}$ was 0.96. The FA $A \times E$ model (logl = -2519.2) was a significantly better fit ($p < 0.001$) than the compound symmetry ($A + A \times E$) model.

Best linear unbiased estimation of the location-by-season environment effects for fruit maturity timing using the FA $A \times E$ model (Fig. 6) indicated that on average, and after accounting for unbalanced genetic treatments across

GRM - diagonals

<figure><img src="image_7.png" /></figure>

GRM - off-diagonals

<figure><img src="image_8.png" /></figure>

A - diagonals

<figure><img src="image_9.png" /></figure>

A - off-diagonals

<figure><img src="image_10.png" /><figcaption>**Fig. 5** Diagonal and off-diagonal estimated realised (SNP-derived, GRM) and expected (pedigree-derived, A) additive relationship matrices elements for germplasm evaluated in this study</figcaption></figure>
environments, fruit maturity was significantly earlier at Balandran, France, in 1997 (144 Julian days) than all other environments. However, there was little absolute difference among the European environments, averaging <9 days. In contrast, average timing of fruit maturity at Prosser USA differed only by 2 days between year, but was around 36 days later than in Europe.

Estimated phenotypic variance differed among locations, but was more consistent between seasons within locations (Table 3). The most variable locations were Balandran and Forlí, with Prosser the least. The heterogeneity in phenotypic variance was mostly a consequence of heterogeneity in estimated additive genomic variance, particularly the lower additive genetic variance at Prosser. Estimates of individual narrow-sense heritability were >0.83 for all European environments, but lower (0.60) for the Prosser environments.

<figure><img src="image_11.png" /><figcaption>**Fig. 6** Least square means of location-by-season effects from most parsimonious G × E model. Environments with the same letter were not significantly different (p = 0.05)</figcaption></figure>

Estimated additive genomic correlations among the European locations and between years within these locations from a FA A × E model with the GRM were close to 1 (Table 4). In contrast, the lowest estimated additive genomic correlation was between seasons at the Prosser location (r<sub>A</sub> = 0.8). Estimates of the correlations between the Prosser and European environments were between 0.87 and 0.92.

## Accuracy

The overall mean genomic predictive ability of using genome-wide SNP information from phenotyped individuals at a particular location to predict performance of individuals untested at that location was 0.81, with mean PACC 0.91 ± 0.24 (Table 5, within location). Predictive ability of individual location-by-set-by-season validation populations ranged from 0.60 to 0.98, with individual PACC between 0.65 ± 0.30 to 1.06 ± 0.15 (full results not presented). Predictive ability averaged across sets within seasons was positively correlated (r = 0.80, n = 8) with estimated narrow-sense heritability.

Predictive ability and accuracy for the second examined scenario of predicting into an environment for which no phenotypic information is available (Table 5, across location) was generally lower than for the first scenario, except for prediction into the Forlí location. The lowest prediction accuracy was observed for prediction of genomic breeding value at Prosser.

## Predicted G × E effects

The high estimates of additive genomic correlation among environments supported accurate prediction of genomic breeding values for fruit maturity timing into the various environments for all individuals, even if they had not been phenotypically assessed in those environments (Fig. 7). As suggested by the high estimated additive

<table><caption>**Table 4** Estimates of additive genomic correlation of fruit maturity timing (Julian days) among eight location-by-season environments from the fit of a FA A × E model using a genomic relationship matrix</caption><thead><tr><th rowspan="2"></th><th colspan="2">Balandran.FR</th><th colspan="2">Bourran.FR</th><th colspan="2">Forlí.IT</th><th colspan="2">Prosser.US</th></tr><tr><th>1997</th><th>1999</th><th>2014</th><th>2015</th><th>2014</th><th>2015</th><th>2011</th><th>2012</th></tr></thead><tbody><tr><td rowspan="2">Balandran.FR</td><td>1997</td><td>1</td><td>1</td><td>0.98</td><td>1</td><td>1</td><td>0.89</td><td>0.92</td></tr><tr><td>1999</td><td></td><td>1</td><td>0.98</td><td>1</td><td>1</td><td>0.89</td><td>0.92</td></tr><tr><td rowspan="2">Bourran.FR</td><td>2014</td><td></td><td></td><td>0.98</td><td>1</td><td>1</td><td>0.89</td><td>0.92</td></tr><tr><td>2015</td><td></td><td></td><td></td><td>0.98</td><td>0.98</td><td>0.87</td><td>0.90</td></tr><tr><td rowspan="2">Forlí.IT</td><td>2014</td><td></td><td></td><td></td><td></td><td>1</td><td>0.89</td><td>0.92</td></tr><tr><td>2015</td><td></td><td></td><td></td><td></td><td></td><td>0.89</td><td>0.92</td></tr><tr><td rowspan="2">Prosser.US</td><td>2011</td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.82</td></tr><tr><td>2012</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>**Table 5** Predictive ability, prediction accuracy and standard error of prediction accuracy by location and season estimated from within-location and across-location validation</caption><thead><tr><th rowspan="2">Location</th><th rowspan="2">Season</th><th colspan="5">Within location</th><th colspan="5">Across location</th></tr><tr><th>nG × L</th><th>nG × LS</th><th>PA</th><th>PACC</th><th>se.PACC</th><th>nG × L</th><th>nG × LS</th><th>PA</th><th>PACC</th><th>se.PACC</th></tr></thead><tbody><tr><td rowspan="2">Balandran.FR</td><td>1997</td><td>14</td><td>12</td><td>0.86</td><td>0.90</td><td>0.29</td><td>63</td><td>61</td><td>0.67</td><td>0.70</td><td>0.13</td></tr><tr><td>1999</td><td></td><td>10</td><td>0.88</td><td>0.91</td><td>0.32</td><td></td><td>50</td><td>0.76</td><td>0.79</td><td>0.14</td></tr><tr><td rowspan="2">Bourran.FR</td><td>2014</td><td>38</td><td>37</td><td>0.81</td><td>0.89</td><td>0.16</td><td>193</td><td>187</td><td>0.59</td><td>0.65</td><td>0.07</td></tr><tr><td>2015</td><td></td><td>38</td><td>0.84</td><td>0.90</td><td>0.16</td><td></td><td>192</td><td>0.64</td><td>0.68</td><td>0.07</td></tr><tr><td rowspan="2">Forli.IT</td><td>2014</td><td>11</td><td>10</td><td>0.81</td><td>0.88</td><td>0.32</td><td>56</td><td>50</td><td>0.84</td><td>0.91</td><td>0.14</td></tr><tr><td>2015</td><td></td><td>11</td><td>0.79</td><td>0.86</td><td>0.30</td><td></td><td>55</td><td>0.86</td><td>0.94</td><td>0.13</td></tr><tr><td rowspan="2">Prosser.US</td><td>2011</td><td>76</td><td>46</td><td>0.79</td><td>1.02</td><td>0.15</td><td>384</td><td>231</td><td>0.43</td><td>0.55</td><td>0.07</td></tr><tr><td>2012</td><td></td><td>71</td><td>0.73</td><td>0.95</td><td>0.12</td><td></td><td>360</td><td>0.44</td><td>0.57</td><td>0.05</td></tr></tbody></table>

Also shown is the average number of individuals in each validation population per location (nG × L) and per location-by-season (nG × LS) for each validation approach PA predictive ability, PACC prediction accuracy, se.PACC standard error of prediction accuracy

<figure><img src="image_12.png" /><figcaption>**Fig. 7** Scatter plot of genomic predicted breeding value of 550 individuals for fruit maturity timing (Julian days) averaged across three locations in Europe compared to genomic predicted breeding value in Prosser, USA</figcaption></figure>

genomic correlations between the European environments and Prosser, the relative ranking of predictions for individuals of average fruit maturity timing across the three European locations was largely consistent with the ranking at Prosser. The earliest predicted average fruit maturity timings across the three European locations were for 'Precoce d'Isigny', 'Belle de St Denis', 'Abouriou', 'Guigne Précoce du Marché' and 'Scwecja', which were 22–25 days earlier than the mean (30 May). Individuals predicted to the latest (19–23 days later than the mean) across the European locations were 'Ambrunés' and four selections from the Pacific Northwest Sweet Cherry

Breeding Program. At Prosser, the five predicted earliest individuals were four of the five earliest across the European locations (all except 'Scwecja') and one selection from the US breeding program, which compared to the mean fruit maturity date for Prosser (2 July) were 12–13 days earlier. The five individuals predicted to be the latest in the European locations were also those predicted to be the latest at Prosser (13–14 days later than the mean date).

## Discussion

### Genetic structure of germplasm

The degree of decay in LD among loci with increasing physical distance across the sweet cherry germplasm in this study is comparable, but slightly weaker, than that reported for only the Bourran germplasm⁴⁰. In this current study, the threshold for which LD was <0.2 occurred between 140 and 240 kb, suggesting slightly stronger LD than observed in Campoy et al.⁴⁰ (0.2 threshold around 100 kb). The stronger LD in our study might be a consequence of the inclusion of full-sib families from the Prosser location. The degree of decay in LD with increasing distance observed in this study and that of Campoy et al.⁴⁰ is similar to that recently reported in a Belgian apple collection⁵⁹, but higher than that reported for peach (<0.2 between 800 and 1800 kb⁶⁰) and an earlier study in apple 1000 kb⁶¹. Campoy et al.⁴⁰ suggested that a lower LD decay reported for peach is a consequence of the presence of self-compatibility in that crop leading to longer regions of homozygosity from inbreeding through selfing. Higher LD in the apple study⁶¹ was likely a consequence of small effective population size in the seven pedigree-connected full-sib families used.

The lack of obvious genetic structure among germplasm displayed as the first two eigenvalue loadings for the SNP-
derived relationship matrix suggests that there is con-
siderable identity of chromosome segments among
germplasm evaluated in the four locations, with indivi-
duals being as related across locations as within the same
location. These results are consistent with the hypothesis
of a bottleneck during domestication of this crop⁶². Using
the same genetic data for the Bourran population as used
here, Campoy et al.⁴⁰ identified two main groups of
individuals representing landraces and bred cultivars.
Further analyses decomposed the groups into nine clus-
ters, with North American-derived cultivars in the Bour-
ran collection allocated to only two of these clusters⁴⁰.
The observed wider diversity of the RosBREED germ-
plasm is not surprising as individuals were deliberately
chosen to maximise diversity by including germplasm that
had not been commonly used previously in North
American breeding programs along with standard North
American cultivars and their available ancestors⁴².

The contrast in distribution of diagonals and off-
diagonals between the SNP and pedigree-derived rela-
tionship matrices observed here has been discussed by
others⁶³. The SNP-derived relationship matrix (GRM) is
able to capture unknown pedigree relationships, and
heterogeneity in relationships due to Mendelian sampling,
that is not possible using only pedigree information to
estimate relationship coefficients ³⁶. Differences between
the pedigree-derived and SNP-derived relationship
matrices might also be a consequence of incorrect or
incomplete pedigree information or improper scaling.
The high frequency of pedigree-derived pairwise rela-
tionship coefficients (off-diagonals) with a value of 0,
compared to the more balanced distribution of realised
relationships estimated from genotypic data, suggests that
pedigree records are incomplete. Germplasm individuals
are generally assumed to be unrelated in the base gen-
eration of the pedigree-derived relationship matrix, but
they might instead be related and those relationships can
be captured in the GRM. Estimation of the GRM is also
dependent on the estimated allele frequency for each
SNP³⁵. Ideally, allele frequencies for the base population
are needed, but they are difficult to obtain. Allele fre-
quencies in our study were computed from available
genotyped accessions and not from a base population,
which could have resulted in improper scaling. None-
theless, Wang et al.⁶⁴ suggests with strong quality control
(e.g., call rate ≥0.9, addressing parent-progeny conflicts,
minor allele frequency >0.05, good pedigree depth), dif-
ferences between the elements of matrices can be
reduced.

The higher variance of the GRM diagonal and off-
diagonal elements, compared to those for *A*, suggests
predictions using a realised relationship matrix will be
more accurate³⁶. In this study in particular, the GRM
might also more accurately estimate additive effects

compared to a relationship matrix based on expected
relationship coefficients derived from pedigree informa-
tion, as there were only a small number of offspring per
family and so there might have been insufficient replica-
tion of Mendelian sampling events within each family to
accurately obtain a family average, which is an important
source of information for prediction of genetic effects⁶⁵.

Estimates >1 for the expected relationship coefficients of an individual with itself derived from pedigree information, and large realised relationship coefficients estimated from genotypic data, suggest that some individuals included in this study are inbred, as was also observed in western larch by, e.g., Klápste et al.⁶⁶ This conclusion agrees with general knowledge of sweet cherry cultivar pedigrees, which involved a strong bottleneck, repeated use of a small number of parents, and crossing of individuals that shared recent ancestors⁴⁰,⁶². For example, breeding programs have continuously used the cultivar 'Stella' and its descendants in breeding for self-compatibility, presumably dramatically reducing the genetic diversity of released cultivars¹⁴.

Model fitting

This study employed a FA parameterisation to estimate
parameters of the genetic-by-environment covariance
matrix. While FA models have similarities with additive
main effects multiplicative interactions (AMMI)³⁰,⁶⁷
approaches, AMMI is a fixed effects approach, whereas
FA mixed models offer more flexibility such as fitting
genotypes and genotype-by-environment interaction as
random effects³⁰. A mixed model approach more easily
deals with unbalanced datasets and supports incorpora-
tion of heterogeneous variances and covariances among
environments—all of which are essential features of our
dataset. Superiority of FA models over AMMI models
have been demonstrated using empirical MET eucalypt
and teak datasets⁶⁸.

Despite the high average additive genomic correlation
among environments ($\bar{r}_{A \times E}$) estimated from a simple
compound symmetry model, the FA $G \times E$ model was a
significantly better fit to the data as it allowed hetero-
geneous variances among locations and heterogeneous
pairwise covariances (and hence correlations) among
locations. The compound symmetry structure constrains
all variances and all covariances to be the same, which is
an unrealistic assumption particularly when the environ-
mental sampling is complex as demonstrated here. In
contrast, the FA model captured the greater stability of
genetic effects among and within the European locations
compared to genetic effects at the Prosser location.

This study only applied an additive genetic effects
model. Theoretically, inclusion of significant non-additive
effects, if present, is expected to increase PACC and
support selection for clonal value (i.e., total genetic effects
= additive + non-additive genetic effects) where individuals can be vegetatively propagated⁶⁹⁻⁷¹ as sweet cherry can be and is commercially. However, separating additive and non-additive genetic effects requires powerful designs or clonal replication of individuals. Moreover, fitting these effects separately in genomic prediction models have not shown meaningful improvement in accuracy, e.g. ref. ⁷², unless dominance was the predominant genetic effect and broad-sense heritability was >0.6⁷¹. Dominance effects for days to maturity have been estimated to compose a minor component of genetic variance in sweet cherry¹¹.

The FA G × E model supported the combination of adjustment of phenotypes for fixed effects and estimation of genetic parameters into a multivariate single-step analysis. The realised relationship matrix was used to describe correlated genetic effects across environments so that data from various sources could be combined. Prediction accuracy from multivariate approaches is generally higher than that from univariate approaches and supports the prediction of genetic effects for individuals in environments in which they have not been tested by leveraging of correlated information²⁹,⁷³,⁷⁴. In contrast to the single-step approach taken here, other studies have applied a two-stage analysis which adjusts phenotypes for fixed effects prior to fitting of the genetic model, e.g., refs. ³⁶,⁷². A single-step analysis as reported here was possible due to the small size of the population. Single-step analyses allow incorporation of complex models in the analysis⁷⁵. While Welham et al.⁷⁶ suggest that PACC from a single-step analysis is similar to that from a two-stage analysis, de los Campos et al.⁷⁷ have reported that two-stage approaches for genomic prediction might produce biased genetic marker effects and correlation among residuals.

The estimates of PACC in all environments reported here are considerably higher than those reported previously for other traits in other horticultural tree crops⁷²,⁷⁸,⁷⁹. The within-location validation strategy adopted here is similar to the second validation strategy adopted by Kumar et al.⁷² Here, training populations were constructed by excluding a subset of individuals within a location to simulate the strategy of predicting the performance of an individual within a target environment for which phenotypic data was available from other genotyped individuals. The high accuracy of this strategy (>0.8) might be explained in part by the high apparent heritability of fruit maturity timing (>0.60). The high prediction accuracies of the model for individuals at the Balandran and Bourran locations (>0.85) are also likely to be a consequence of the high genomic correlation (0.98), and large number of individuals replicated (68), between these locations. The high PACC for individuals at the Prosser location (>0.95) is likely to be a consequence of the random sampling strategy used to form the validation population as we did not account for the known complex

family (and ancestral) structure. It is well known that close relationships between the training and validation populations increases PACC⁸⁰. The lower accuracy of the across-location validation strategy of predicting across environments (0.50–0.95) is probably a consequence of the confounding of G × E with within-environment PACC.

## Implications of results for germplasm management

This study appears to be the first to quantify differences among locations in timing of fruit maturity in sweet cherry, particularly on a global scale. The advantage of undertaking a combined multivariate analysis like presented here is that location means are unbiased because they are adjusted for differences in germplasm composition among environments. In this way, comparisons are not confounded by the use of different germplasm in each environment. Differences among seasons within a location were not as great as differences between the European locations and Prosser, suggesting that fruit maturity across the European locations is generally earlier than at the Prosser location, at least for the seasons sampled here. More locations and seasons evaluated could be used to develop a more general predictive equation. A joint analysis of flowering timing and fruit maturity timing could help determine the extent of this phenological phenomenon.

High heritability estimates for timing of fruit maturity in sweet cherry have been reported previously¹⁰, suggesting that an individual's breeding value for timing of fruit maturity is reliably predicted from the phenotype of that individual. However, our model only included additive genetic effects, and others⁷²,⁸¹ have demonstrated that heritability estimates can be inflated if large non-additive effects are present but omitted from the training model. Deflated estimates of additive genetic variances in complex models might also be a consequence of correlation among additive and non-additive genetic effects⁸².

We are not aware of any previously published estimates of G × E for any traits, including fruit maturity timing, in cherry. Despite the detection of significant G × E in this study, G × E is somewhat low compared to the degree reported for most other traits in horticultural tree crops, e.g., ref. ⁸³. Differences in ranking between seasons at Prosser was the major source of G × E detected in this study. The limited number of locations in this study prohibits the development and testing of a hypothesis to explain this pattern. The generally high correlation among all locations suggests that the timing of fruit maturity of an individual relative to other individuals is expected to be stable across production environments similar to these test locations. This conclusion is supported by the high correlation of predicted genetic values at Prosser with the average of predicted values across the three European
locations (Fig. 7). Results also fit with the general obser-
vation that traits with high heritability exhibit low levels of
G × E, which might be due in part to numerical
constraints²².

The high narrow-sense heritability for fruit maturity
timing reported here provides clear evidence of the large
genetic gain that can be achieved by selecting for parti-
cular fruit maturity timing. The high genetic correlation
among environments for this trait (>0.82) suggests that
individuals can be selected for elite average performance
across environments similar to those examined here, and
elite germplasm in one environment for fruit maturity is
expected to be elite in other environments. However, it is
difficult to generalise this recommendation to new
environments. In addition, other traits might exhibit dif-
ferent G × E patterns. However, the genetic potential of
genotyped but untested individuals can be predicted using
the methods outlined here if genomic and phenotypic
data are available for other individuals and are used to
estimate a prediction model.

Extension of method

This study is based on a model of individual SNP loci in
linkage disequilibrium with chromosome segments con-
taining numerous small-effect QTLs such that replication
of SNP alleles across environments models the replication
of causal alleles at the small-effect QTLs. However, as
previous studies¹⁰,¹² have suggested the presence of large-
effect QTLs contributing to genetic variation of fruit
maturity timing, the model of genetic architecture used in
this study might not be accurate. Incorporation of known
QTL effects in the linear model should improve accuracy,
as it is expected that such a model would more closely
approximate the true distribution of QTL effects⁷⁷. One of
the challenges of building these models is that char-
acterisation of the variation around QTL regions as
individual biallelic SNPs might not capture the complexity
of variation in the region. This variation might be better
modelled using the same genotypic data by describing
SNP haplotypes that contain QTL regions. The level of
LD in cherry is sufficiently high that this information
could be used to phase the SNP genotypes (for example,
with Beagle⁵¹). Otherwise, the model could use separately
obtained genotypic data that describes the allelic com-
position of each individual at large-effect QTLs with
multiple effective alleles.

Incorporation of data from additional environments
would be useful to develop a more general description of
G × E patterns for the trait, enable development and
testing of hypotheses to explain the G × E patterns, and
increase PACC. Only a minor degree of G × E was
detected in this study and thus there is insufficient var-
iation to develop and test hypotheses of environmental
drivers of G × E for fruit maturity timing. Increasing the

environmental range might induce greater G × E. In this
case, it might be that repeatable prediction of G × E could
be achieved through incorporation of environmental
variables or through use of crop growth models that
associate variability in underlying physiological process
with variation in environmental drivers<sup>84,85</sup>.

This study has demonstrated that genome-wide SNP-
derived relationship matrices can be used to quantitatively
describe cryptic genetic relationships. In addition, this
study has demonstrated that such realised relationships
can be used to combine phenotypic data collected on
individuals evaluated across disparate locations into a
single analysis without the need for clonal replication or
knowledge of pedigree relationships. These approaches
now allow the latent value of existing, otherwise uncon-
nected, historical phenotypic datasets collected by
numerous local breeding and cultivar testing programs to
be exploited. Such analyses can increase knowledge of the
genetic architecture of important traits to support ongo-
ing genetic improvement, as well as having the immediate
practical value of predicting genetic effects of individuals
(i.e., breeding or clonal value) in untested environments.
This type of analysis would be informative for other
important horticultural traits of sweet cherry, such as
bloom timing, productivity, fruit size, fruit firmness, fruit
sweetness and fruit acidity, for which phenotypic data are
largely available among breeders and cultivar testers.

Acknowledgements

We would like to thank Dr. Ana Wünsch (Centro de Investigación y Tecnología Agroalimentaria de Aragón, Spain) for contribution of genotypic data for the cultivar 'Hedelfingen'. This work was supported by the USDA National Institute of Food and Agriculture (NIFA)—Specialty Crop Research Initiative project, 'RosBREED: Combining disease resistance with horticultural quality in new rosaceous cultivars' (grant number 2014-51181-22378).

Author details

¹University of Queensland, Queensland Alliance for Agriculture and Food Innovation, Brisbane, QLD 4072, Australia. ²The New Zealand Institute for Plant and Food Research Limited, Hawke's Bay Research Centre, Hastings 4130, New Zealand. ³Department of Horticulture, Washington State University, Pullman, WA 99164, USA. ⁴Department of Horticulture, Michigan State University, East Lansing, MI 48824, USA. ⁵UMR 1332 BFP, INRA, University of Bordeaux, 33140 Nouvelle-Aquitaine, France. ⁶Council for Agricultural Research and Economics (CREA), Fruit Unit of Forlì, Via la Canapona, 1 bis, 47121 Emilia-Romagna, Italy. ⁷Centre Technique Interprofessionnel des Fruits et Légumes (CTIFL), 751 Chemin de Balandran, 30127 Bellegarde, France. ⁸Department of Horticulture, Washington State University, Irrigated Agriculture Research and Extension Center, 24106N Bunn Road, Prosser, WA 99350, USA

Conflict of interest

The authors declare that they have no conflict of interest.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.
## References

1. Mattheis, J. P. & Roberts, R. G. Fumigation of sweet cherry (*Prunus avium* Bing) fruit with low-molecular-weight aldehydes for postharvest decay control. *Plant Dis.* **77**, 810–814 (1993).

2. Alique, R., Martinez, M. A. & Alonso, J. Influence of the modified atmosphere packaging on shelf life and quality of Navalinda sweet cherry. *Eur. Food Res. Technol.* **217**, 416–420 (2003).

3. Wani, A. A., Singh, P., Guld, K., Wani, M. H. & Langowskib, H. C. *Food Packaging and Shelf-life* 86–99 (2014).

4. Kupferman, G. & Sanderson, P. Temperature management and modified atmosphere packing to preserve sweet cherry fruit quality. *Acta Hortic.* **667**, 523–528 (2001).

5. Webster, A. D. in *Cherries: Crop Physiology, Production and Uses* (eds Webster, D. A. & Looney, N. E.) (CAB International, Oxford, 1996).

6. Iglesias, I. et al. El cultivo del cerezo en Espana: produccion, consumo e intercambios comerciales. *Fruticoltura* **48**, 6–39 (2016).

7. Palasciano, M. Ciliegio: tutte le varietà per i nuovi impianti. *L'informatore Agrario* **32**, 21–24 (2015).

8. San Martino, L., Manavella, F. A., Garcia, D. A. & Salato, G. Phenology and fruit quality of nine sweet cherry cultivars in south Patagonia. *Acta Hortic.* **795**, 841–848 (2008).

9. Bargioni, G. in *Cherries: Crop Physiology, Production and Uses* 73–112 (eds Webster, D. A. & Looney, N. E.) (CAB International, Oxford, 1996).

10. Dirlewanger, E. et al. Comparison of the genetic determinism of two key phenological traits, flowering and maturity dates, in three *Prunus* species: peach, apricot and sweet cherry. *Heredity* **109**, 280–292 (2012).

11. Piaskowski, J. et al. Genomic heritability estimates in sweet cherry reveal non-additive genetic variance is relevant for industry-prioritized traits. *BMC Genet.* **19**, 16 (2018).

12. Quero-Garcia, J. et al. QTL detection of important agronomic traits for sweet cherry breeding. *Acta Hortic.* **1020**, 57–64 (2014).

13. Quero-Garcia, J., Schuster, M., Lopez-Ortega, G. & Charlot, G. in *Cherries: Botany, Production and Uses* 60–94 (eds Quero-Garcia, J., lezzoni, A., Pulawska, J. & Lang, G.) (CAB International, Wallingford, 2017).

14. lezzoni, A. et al. in *Cherries: Botany, Production and Uses* 35–59 (eds Quero-Garcia, J., lezzoni, A., Pulawska, J. & Lang, G.) (CAB International, Wallingford, 2017).

15. Hrotko, K. & Rozpara, E. in *Cherries: Botany, Production and Uses* (eds Quero-Garcia, J., lezzoni, A., Pulawska, J. & Lang, G.) (CAB International, Wallingford, 2017).

16. Sansavini, S. & Lugli, S. Sweet cherry breeding programs in Europe and Asia. *Acta Hortic.* **795**, 41–58 (2008).

17. Choi, C. & Kappel, F. Inbreeding, coancestry, and founding clones of sweet cherries from North America. *J. Am. Soc. Hortic. Sci.* **129**, 535–543 (2004).

18. lezzoni, A., Schmidt, H. & Albertini, A. in *Genetic Resources of Temperate Fruit and Nut Crops* Vol. 1, 111–173 (eds Moore, J. N. & Ballington Jr, J. R.) (ISHS, Wageningen, 1990).

19. Allard, R. W. & Bradshaw, A. D. Implications of genotype-environment interaction in applied plant breeding. *Crop Sci.* **4**, 503–508 (1964).

20. Cooper, M. & Delacy, I. H. Relationships among analytical methods used to study genotypic variation and genotype-by-environment interaction in plant-breeding multi-environment experiments. *Theor. Appl. Genet.* **88**, 561–572 (1994).

21. Matheson, A. C. & Cotterill, P. P. Utility of genotype × environment interactions. *For. Ecol. Manag.* **30**, 159–174 (1990).

22. Hardner, C. M., Evans, K., Brien, C., Bliss, F. & Peace, C. Genetic architecture of apple fruit quality traits following storage and implications for genetic improvement. *Tree Genet Genomes* **12** (2016) https://doi.org/10.1007/s11295-016-0977-z.

23. Miller, S. et al. Performance of apple cultivars in the 1995 NE-183 regional project planting: III. Fruit sensory characteristics. *J. Am. Pomol. Soc.* **59**, 28–43 (2005).

24. Miller, S. et al. Performance of apple cultivars in the 1995 NE-183 regional project planting: II. Fruit quality characteristics. *J. Am. Pomol. Soc.* **58**, 65–77 (2004).

25. Miller, S. S. et al. Effect of cultivar and site on fruit quality as demonstrated by the NE-183 regional project on apple cultivars. *Hortotechnology* **15**, 886–895 (2005).

26. Hardner, C. M., Bally, I. S. E. & Wright, C. L. Prediction of breeding values for average fruit weight in mango using a multivariate individual mixed model. *Euphytica* **186**, 463–477 (2012).

27. Falconer, D. S. The problem of environment and selection. *Am. Nat.* **86**, 293–298 (1952).

28. Burdon, R. D. Genetic correlation as a concept for studying genotype-environment interaction in forest tree breeding. *Silvae Genet.* **26**, 168–175 (1977).

29. Henderson, C. R. Best linear unbiased prediction of breeding values not in model for records. *J. Dairy Sci.* **60**, 783–787 (1977).

30. Smith, A. B., Cullis, B. R. & Thompson, R. The analysis of crop cultivar breeding and evaluation trials: an overview of current mixed model approaches. *J. Agric. Sci.* **143**, 449–462 (2005).

31. Smith, A., Cullis, B. & Thompson, R. Analyzing variety by environment data using multiplicative mixed models and adjustments for spatial field trend. *Biometrics* **57**, 1138–1147 (2001).

32. Kelly, A. M., Smith, A. B., Eccleston, J. A. & Cullis, B. R. The accuracy of varietal selection using factor analytic models for multi-environment plant breeding trials. *Crop Sci.* **47**, 1063–1070 (2007).

33. Hardner, C. M., Dieters, M., Dale, G., DeLacy, I. & Basford, K. E. Patterns of genotype-by-environment interaction in diameter at breast height at age 3 for eucalypt hybrid clones grown for reafforestation of lands affected by salinity. *Tree Genet. Genomes* **6**, 833–851 (2010).

34. Arief, V. N. et al. Evaluating testing strategies for plant breeding field trials: redesigning a CIMMYT International Wheat Nursery. *Crop Sci.* **55**, 164–177 (2015).

35. VanRaden, P. M. Efficient methods to compute genomic predictions. *J. Dairy Sci.* **91**, 4414–4423 (2008).

36. Hayes, B. J., Visscher, P. M. & Goddard, M. E. Increased accuracy of artificial selection by using the realized relationship matrix. *Genet. Res.* **91**, 47–60 (2009).

37. Visscher, P. M., Hill, W. G. & Wray, N. R. Heritability in the genomics era - concepts and misconceptions. *Nat. Rev. Genet.* **9**, 255–266 (2008).

38. Speed, D. & Balding, D. J. Relatedness in the post-genomic era: is it still useful? *Nat. Rev. Genet.* **16**, 33–44 (2015).

39. Burguono, J., de los Campos, G., Weigel, K. & Crossa, J. Genomic prediction of breeding values when modeling genotype × environment interaction using pedigree and dense molecular markers. *Crop Sci.* **52**, 707–719 (2012).

40. Campoy, U. A. et al. Genetic diversity, linkage disequilibrium, population structure and construction of a core collection of *Prunus avium* L. landraces and bred cultivars. *BMC Plant Biol.* **16**, 15 (2016).

41. Meier, U. et al. Phänologische entwicklungsstadien des kernobstes (*Malus domestica* Borkh. und *Pyrus communis* L.), des Steinobstes (*Prunus-Arten*), der Johnnisbeere (*Ribes-Arten*) und der Erdbbeere (*Fragaria × ananassa* Duch.). *Nachr. Dtsch. Pflanzenschutzd.* **46**, 141–153 (1994).

42. Peace, C. P., Luby, J. J., van de Weg, W. E., Bink, M. C. A. M. & Lezzoni, A. F. A strategy for developing representative germplasm sets for systematic QTL validation, demonstrated for apple, peach, and sweet cherry. *Tree Genet. Genomes* 1–16 (2014) *Tree Genet. Genomes* 10:1–16.

43. Chavoshi, M. et al. Phenotyping protocol for sweet cherry (*Prunus avium* L.) to facilitate an understanding of trait inheritance. *J. Am. Pomol. Soc.* **68**, 125–134 (2014).

44. Lezzoni, A. et al. RosBREED: enabling marker-assisted breeding in Rosaceae. *Acta Hortic.* **859**, 389–394 (2010).

45. Peace, C. et al. Development and evaluation of a genome-wide 6K SNP array for diploid sweet cherry and tetraploid sour cherry. *PLoS ONE* **7**, 48305 (2012).

46. GenomeStudio Genotyping Module v1.0, User Guide (Illumina Inc., San Diego, 2010).

47. Di Guardo, M. et al. ASSiST: an automatic SNP scoring tool for in- and out-breeding species. *Bioinformatics* **31**, 3873–3874 (2015).

48. Bink, M. et al. Bayesian QTL analyses using pedigreed families of an out-crossing species, with application to fruit firmness in apple. *Theor. Appl. Genet.* **127**, 1073–1090 (2014).

49. Voorrips, R. E., Bink, M., Kruisselbrink, J. W., Koehorst-van Putten, H. J. J. & van de Weg, W. E. PediHaplotyper: software for consistent assignment of marker haplotypes in pedigrees. *Mol. Breed.* **36**, 10 (2016).

50. Browning, S. R. & Browning, B. L. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. *Am. J. Hum. Genet.* **81**, 1084–1097 (2007).

51. Browning, B. L. & Browning, S. R. Improving the accuracy and efficiency of identity-by-descent detection in population data. *Genetics* **194**, 459 (2013).

52. Wimmer, V., Albrecht, T., Auinger, H. J. & Schon, C. C. synbreed: a framework for the analysis of genomic prediction data using *R. Bioinformatics* **28**, 2086–2087 (2012).
53. Brien, C. J. & Demetrio, C. G. B. Using the randomisation in specifying the ANOVA model and table for properly and improperly replicated grazing trials. *Aust. J. Exp. Agric.* **38**, 325–334 (1998).

54. Gilmour, A. R., Gogel, B., Cullis, B. R. & Thompson, R. ASReml User Guide Release 3.0 (VSN International Ltd, Hemel Hempstead, 2009).

55. Henderson, C. R. Best linear unbiased estimation and prediction under a selection model. *Biometrics* **31**, 423–447 (1975).

56. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. *Gigascience* **4**, 16 (2015).

57. Guajardo, V. et al. Construction of high density sweet cherry (*Prunus avium* L.) linkage maps using microsatellite markers and SNPs detected by genotyping-by-sequencing (GBS). *PLoS ONE* **10**, 17 (2015).

58. Fisher, R. A. *Statistical Methods for Research Workers* (Oliver and Boyd, Edinburgh 1954).

59. Vanderzande, S., Micheletti, D., Troggio, M., Davey, M. W. & Keulemans, J. Genetic diversity, population structure, and linkage disequilibrium of elite and local apple accessions from Belgium using the IRSC array. *Tree Genet. Genomes* **13** (2017) https://doi.org/10.1007/s11295-017-1206-0.

60. Micheletti, D. et al. Whole-genome analysis of diversity and SNP-major gene association in peach germplasm. *PLoS ONE* **10**, 19 (2015).

61. Kumar, S. et al. Genomic selection for fruit quality traits in apple (*Malus* × domestica Borkh). *PLoS ONE* **7** (2012).

62. Mariette, S. et al. Population structure and genetic bottleneck in sweet cherry estimated with SSRs and the gametophytic self-incompatibility locus. *BMC Genet.* **11**, 13 (2010).

63. Bérénos, C., Ellis, P. A., Pilkington, J. G. & Pemberton, J. M. Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches. *Mol. Ecol.* **23**, 3434–3451 (2014).

64. Wang, H., Misztal, I. & Legarra, A. Differences between genomic-based and pedigree-based relationships in a chicken population, as a function of quality control and pedigree links among individuals. *J. Anim. Breed. Genet.* **131**, 445–451 (2014).

65. Sales, J. & Hill, W. G. Effect of sampling errors on efficiency of selection indexes 1. Use of information from relatives for single trait improvement. *Animal Production. Anim. Sci.* **2**, 1–17 (1976).

66. Klápste, J., Lstibúrek, M. & El-Kassaby, Y. A. Estimates of genetic parameters and breeding values from western larch open-pollinated families using marker-based relationship. *Tree Genet. Genomes* **10**, 241–249 (2014).

67. Gauch, H. G. Statistical analysis of yield trials by AMMI and GGE. *Crop Sci.* **46**, 1488–1500 (2006).

68. Resende, M. V. D. & Thompson, R. Factor analytic multiplicative mixed models in the analysis of multiple experiments. *Rev. Mat. Estat.* **22**, 31–52 (2004).

69. Henderson, C. R. Best linear unbiased prediction of nonadditive genetic merits on noninbred populations. *J. Anim. Sci.* **60**, 111–117 (1985) https://doi.org/10.1371/journal.pone.0036674.

70. de Boer, I. J. M. & van Arendonk, J. A. M. Prediction of additive and dominance effects in selected or unselected populations with inbreeding. *Theor. Appl. Genet.* **84**, 451–459 (1992).

71. Denis, M. & Bouvet, J. M. Efficiency of genomic selection with models including dominance effect in the context of eucalyptus breeding. *Tree Genet. Genomes* **9**, 37–51 (2013).

72. Kumar, S. et al. Genome-enabled estimates of additive and nonadditive genetic variances and prediction of apple phenotypes across environments. *G3 Genes Genomes Genet.* **5**, 2711–2718 (2015).

73. Hardner, C. M., Healey, A. L., Downes, G., Herberling, M. & Gore, P. L. Improving prediction accuracy and selection of open-pollinated seed-lots in *Eucalyptus dunnii* maiden using a multivariate mixed model approach. *Ann. For. Sci.* **73**, 1035–1046 (2016).

74. Thompson, R. & Meyer, K. A review of theoretical aspects in the estimation of breeding values for multi-trait selection. *Livest. Prod. Sci.* **15**, 299–313 (1986).

75. Oakey, H. et al. Genomic selection in multi-environment crop trials. *G3 Genes Genomes Genet.* **6**, 1313–1326 (2016).

76. Welham, S. J., Gogel, B. J., Smith, A. B., Thompson, R. & Cullis, B. R. A comparison of analysis methods for late stage variety evaluation trials. *Aust. N. Z. J. Stat.* **52**, 125–149 (2010).

77. de los Campos, G., Hickey, J. M., Pong-Wong, R., Daetwyler, H. D. & Calus, M. P. L. Whole-genome regression and prediction methods applied to plant and animal breeding. *Genetics* **193**, 327–345 (2013).

78. Biscarini, F. et al. Genome-enabled predictions for fruit weight and quality from repeated records in European peach progenies. *BMC Genomics.* **18**, 15 (2017).

79. Muranty, H. et al. Accuracy and responses of genomic selection on key traits in apple breeding. *Hortic. Res.* **2**, 12 (2015).

80. Habier, D., Fernando, R. L. & Dekkers, J. C. M. The impact of genetic relationship information on genome-assisted breeding values. *Genetics* **177**, 2389–2397 (2007).

81. Munoz, P. R. et al. Unraveling additive from nonadditive effects using genomic relationship matrices. *Genetics* **198**, 1759–1768 (2014).

82. Melchinger, A. E. Means, variance, and covariances between relatives in hybrid populations with disequilibrium in parent populations. In *Proc. 2nd International Conference on Quantitative Genetics* (eds. Weir, B. S., Eisen, E. J. & Goodman, M. M.) 400–415 (Sinauer, Sunderland, 1988).

83. Li, Y. J., Suontama, M., Burdon, R. D. & Dungey, H. S. Genotype by environment interactions in forest tree breeding: review of methodology and perspectives on research and application. *Tree Genet. Genomes* **13**, 18 (2017).

84. Technow, F., Messina, C. D., Totir, L. R. & Cooper, M. Integrating crop growth models with whole genome prediction through approximate Bayesian computation. *PLoS ONE* **10**, 20 (2015).

85. van Eeuwijk, F. A., Bustos-Korts, D. V. & Malosetti, M. What should students in plant breeding know about the statistical aspects of genotype × environment interactions? *Crop Sci.* **56**, 2119–2140 (2016).