Genotype-by-environment and QTL-by-environment interactions in sweet cherry (Prunus avium L.) for flowering date

In sweet cherry (Prunus avium L.), flowering date is strongly dependent on the environment conditions and, therefore, is a trait of major interest for adaptation to climate change. Such trait can be influenced by genotype-by-environment interaction (G×E), that refers to differences in the response of genotypes to different environments. If not taken into account, G×E can reduce selection accuracy and overall genetic gain. However, little is known about G×E in fruit tree species. Flowering date is a highly heritable and polygenic trait for which many quantitative trait loci (QTLs) have been identified. As for the overall genetic performance, differential expression of QTLs in response to environment (QTL-by-environment interaction, QTL×E) can occur. The present study is based on the analysis of a multi-environment trial (MET) suitable for the study of G×E and QTL×E in sweet cherry. It consists of a sweet cherry F1 full-sib family (n = 121) derived from the cross between cultivars ‘Regina’ and ‘Lapins’ and planted in two copies in five locations across four European countries (France, Italy, Slovenia and Spain) covering a large range of climatic conditions. The aim of this work was to study the effect of the environment on flowering date and estimate G×E, to carry QTL detection in different environments in order to study the QTL stability across environments and to estimate QTL×E. A strong effect of the environment on flowering date and its genetic control was highlighted. Two large-effect and environment-specific QTLs with significant QTL×E were identified on linkage groups (LGs) 1 and 4. This work gives new insights into the effect of the environment on a trait of main importance in one of the most economically important fruit crops in temperate regions. Moreover, molecular markers were developed for flowering date and a strategy consisting in using specific markers for warm or cold regions was proposed to optimize marker-assisted selection (MAS) in sweet cherry breeding programs.


Introduction
The phenotype of an individual is the result of a combination of the effect of its genotype, the external environment, and the interaction between the genotype and environmental variations (Genotype-by-environment interactions, G×E). G×E is a common phenomenon referring to differences in the response of genotypes to different environments (Allard and Bradshaw, 1964;Falconer and Mackay, 1996). The presence of G×E can affect the genetic advance obtained from selection (i.e. superior cultivars in one environment are not necessarily superior in another environment) and therefore is a major concern for plant breeders. Two major sources of G×E exist: (1) rank-change interaction (or crossover interaction), when genotypes are ranked in different orders in different environments; and (2) scale-change interaction (or level-of-expression interaction), when genotypic differences vary across environments (heterogeneity of genetic variances across environments). Conventionally, multi-environment trials (METs) are conducted to assess the performance of a common set of genotypes in different environments and evaluate G×E (Smith et al., 2005). Many statistical methods have been developed for a precise description of G×E and nowadays, linear mixed models where genotypes are treated as random effect factors are frequently used (Smith et al., 2005;Malosetti et al., 2013;van Eeuwijk et al., 2016).
In perennial plant species, G×E has been mostly studied in forest tree species for traits related to tree height and trunk diameter (Li et al., 2017;Bird et al., 2022). More recently, horticultural fruit tree species such as apple (Jung et al., 2020), macadamia (Hardner, 2017), sweet cherry (Hardner et al., 2019) and peach  have also been studied. For instance, in sweet cherry, low additive G×E and high additive genomic correlations among environments were estimated for maturity date in a germplasm of 597 cultivars, accessions and unselected offspring planted in three locations in Europe and one location in the USA (Hardner et al., 2019). In apple, genomic G×E explained 18 and 12% of the phenotypic variance of floral emergence and harvest date in a population of 534 genotypes planted across six European countries (Jung et al., 2020). In summary, in these horticultural crops, several traits related to phenology, yield and fruit quality (e.g. sweetness) have been studied. Nevertheless, little is known about the environmental stability of genetic effects for flowering time.
In the current global warming context, flowering date (FD) is a trait of major interest in temperate fruit tree species such as sweet cherry (Prunus avium L.). FD is highly dependent on the climatic conditions, therefore, breeding programs tend to develop early and late blooming cultivars according to their area of production (Quero-Garcıá et al., 2017). Early flowering cultivars are promoted in warm regions to avoid high temperatures during flowering while late cultivars are best suited for cold areas to avoid spring frost damages. Moreover, FD is dependent on the dormancy period in which temperate fruit trees enter during winter to stop meristem activity and prevent frost damages (Lang et al., 1987). The length of this period varies according to climate conditions and individuals, as specific amounts of chill and heat, known as 'chilling requirements' (CRs) and 'heat requirements' (HRs), are required to release dormancy (Alburquerque et al., 2008). However, nowadays, few low-chilling varieties (with CRs varying from 300 to 800 chilling hours) are available (Wenden et al., 2017). This is even more problematic in the context of global warming, with the increase of temperatures in autumn and winter, which has already provoked serious production losses on cultivars with high CR (Quero-Garcia, comm. pers.). Furthermore, increases in spring temperatures with global warming have entailed a significant advance of flowering dates of cherry cultivars in numerous production areas (Luedeling, 2012;Wenden et al., 2017), with a subsequent increase in risk of frost damage. Although some work has been conducted in several breeding programs to investigate tolerance or resistance to cold, and hence to frost damage, very few cultivars carrying these favorable traits have been released so far and none has reached commercial importance (Quero-Garcıá et al., 2017;Wenden et al., 2017). To date, knowledge about G×E and the stability of whole genome effect for FD in sweet cherry is missing.
FD is a quantitative trait with high broad sense heritability and genetic approaches have led to the identification of many FD QTLs, highlighting a complex genetic control. Due to the high genomic synteny within the Prunus genus, FD QTLs have been identified in sweet cherry and other species in similar chromosomal regions (Dirlewanger et al., 2004;Dirlewanger et al., 2012;Castède et al., 2014;Sańchez-Peŕez et al., 2014;Calle et al., 2020). Although QTLs have been detected on all linkage groups (LGs), the two largest-effect loci were located on LGs 1 and 4 (Fan et al., 2010;Dirlewanger et al., 2012;Castède et al., 2014;Sańchez-Peŕez et al., 2014;Cai et al., 2018;Calle et al., 2020;Branchereau et al., 2022). Castède et al. (2014) dissected sweet cherry FD into CRs and HRs and detected QTLs for FD, CRs and HRs co-localizing in the LG4 region. The QTL on LG1 covers the genomic region of the well-known DORMANCY-ASSOCIATED MADS-box (DAM) genes (Bielenberg et al., 2008;Castède et al., 2015;Calle et al., 2021). Recently, candidate genes involved in auxin responses and splicing have been identified in the QTL on LG4 using the 'Regina' sweet cherry genome sequence and transcriptomic analyses (Branchereau et al., 2022).
In sweet cherry, little is known about the stability of QTL effects for FD across environments. As the overall genetic performance, QTLs can be expressed differently in different environments: a QTL can be significant for a given trait in one environment but not in another. This is called QTL-by-environment interaction (QTL×E) and a QTL with large QTL×E interaction is less stable than a QTL with small QTL×E. In Prunus, Asıńs et al. (1994) were the first to use QTL detections as an approach to study G×E. Authors studied the effect of years on QTL detection for several quantitative traits in almond. They highlighted that only few QTLs behaved homogeneously over the years and showed that the presence of G×E was the most likely cause of this phenomenon.
The presence of G×E is a main concern for breeders as it complicates the identification of superior cultivars and reduces gains from selection. The identification of stable genotypes (low G×E) over a wide range of environments is of main importance, especially in locations associated with strong environmental fluctuations. Moreover, information about QTL×E is essential to develop marker-assisted selection (MAS) for FD according to the area of production.
The main objectives of this study were to (i) estimate G×E in a sweet cherry MET, (ii) perform QTL detection in order to study the QTL stability over environments, and (iii) evaluate QTL×E interactions for FD in sweet cherry. This work should contribute to increase the efficiency of sweet cherry breeding programs.

Plant material
An F 1 sweet cherry full-sib family derived from the cross between 'Regina' and 'Lapins' cultivars was used for this study. 'Regina' is a late blooming German cultivar, whereas 'Lapins' is an early-intermediate blooming cultivar from Canada. This family, hereafter called R×L, consists of clones of 121 hybrids planted in a MET in five locations across four European countries: Forli (northeastern Italy), Maribor (north-eastern Slovenia), Murcia (southeastern Spain), Nimes (south-eastern France) and Toulenne (southwestern France) (Table S1; Figure 1A). Daylenght is similar in the different locations ( Figure 1B) while temperatures and precipitations are highly contrasted (Figures 1C,D). In Maribor, the climate is continental with cold winters and quite warm summers. Important rain precipitations are observed all year Location and environment characterization of the five experimental sites across Europe. (A), location of the five sites in Europe. (B), average day length in the five sites. (C), temperature deviation to the overall monthly mean. (D), average monthly precipitations in the five sites. In (B-D), the climatic data used is from 2010 to 2021. A color scale from blue to red was chosen to represent the sites according to the temperature data.
long. It is highly contrasted to Murcia, which is dry year-round and has mild winters and very hot summers. Climates in Forli, Toulenne and Nimes are intermediate ( Figures 1B-D). Orchards were irrigated in all locations except Maribor and Toulenne. Trees were planted every 2.5 to three meters in rows separated by five meters. G×E studies require replication of genetic effects across MET, therefore, the 121 R×L genotypes were grafted (clonally replicated) in two copies in all environments (rootstocks: Maxma Delbard ® 14 or Colt) and planted in a random design. However, genotypes are not all present in the five locations and the number of replicates per genotype varies between sites (Table S2).

Flowering date phenotyping
Two FD stages were scored in all locations: beginning of flowering (BF), when approximately 10% of the floral buds reached full bloom, and full flowering (FF), when 75% of the floral buds reached full bloom. Trees were observed from three to four times a week during the season to score the different flowering stages in Julian days (JDs). FD was assessed at each location for several seasons: FD was scored in 2018, and 2021in 2017in and 2021in 2017in , 2018in 2017in , 2018in and 2021, 2017, 2018and 2021 in Toulenne. 'Environments' were considered as the combination of the location and the year (season) of the trial. Therefore, the MET consisted of twenty unique location × season environments.

Environment characterization
As temperature is the most important climatic factor for flowering in perennials like sweet cherry, the twenty environments were characterized using temperature data from October to May, a period covering dormancy (endodormancy and ecodormancy) and flowering. In each location separately, from October to May, a daily mean temperature was calculated using temperature data from 2010 to 2021. Then, the temperature data from each season we studied was represented as a deviation to the overall mean. Plots are available in Figures S1-S5, for Forli, Maribor, Murcia, Nimes and Toulenne, respectively. This type of representation allows to visualize in each location when the temperatures were either lower or higher than the mean.
2.4 Flowering date distribution, correlations and heritabilities 2.4.1 Phenotypic data review Distribution, mean, minimum and maximum values of BF and FF were estimated for each location-by-season environment. Additionally, Spearman correlation coefficients between years within each location and between locations for each year were calculated.

Analyses of G×E
In order to describe the phenotype profile of the 121 R×L individuals across environments, norms of reaction were obtained with 'ggplot2' R package (Wickham, 2016). For each individual in each environment, the average phenotypic value was calculated from the replicates. Norms of reaction were obtained across locations in each year of study and across years in each location of study.
A G×E model was fitted to multi-trial data to estimate the genetic architecture of traits. As the MET is unbalanced (e.g. different number of replicates per genotype, some genotypes lacking in several environments, different years of measurement available in different locations), the mixed model approach was used. To reduce the heterogeneity of variance between environments and hence potential influence of heterogeneity of variance on G×E, FD observations within each environment were scaled by the raw phenotypic standard deviation of their respective environment (Hardner, 2017). Analyses were performed for BF and conducted in R using ASReml-R package version 4 (Butler et al., 2017).
The general model of the phenotype of an i th individual in the k th block at the l th trial for the j th season was y ljki = m + e lj + eb ljk + g i + ge lji + r ljki where m was the general mean across all trials, blocks within trials, seasons at trials, and individuals e lj was the fixed effect of the lj th environment (i.e. j th season at the l th trial) eb ljk was the fixed effect of the k th block at lj th environment g i was the average total (i.e., additive + non-additive) genetic effect of the i th individual across environments with distribution N (0,Gg) where Gg was the variance-covariance matrix among average total genetic effects across environments given as Gg = I*v where I was the identity matrix of relationships among individuals and vg was the unknown total genetic variance across environments (where * was the Kronecker product) ge lji was the random environment specific total genetic effect of the i th individual at the lj th environment with distribution N(0,Gge) where Gge was the variance-covariance matrix among environment specific total genetic effects at each of the lj th environments given as Gge = I*Vge where I was the identity matrix of relationships among individuals, and Vge was the covariance matrix of environment specific total genetic effects among environments r ljki was the residual effect for each phenotype observation with distribution N(0,R) where R was a block diagonal matrix of residual covariance among seasons for individuals at each trial, i.e. Vr l where Vr l was given as I l *Vr l where I l was an identity matrix among individuals at the l th trial and Vr l was the residual covariance matrix among seasons at the l th trial. Parameters of the model (i.e. covariance components) were estimated using Restricted Maximum Likelihood implemented in the R package ASREML-R v4 (Butler et al., 2017). Tests of significance for fixed effects were done with Wald test (Kenward and Roger, 1997). The diagonal of the variance covariance matrix was referred to the interaction variance in each environment lj, v G×E (lj), in other words, the variance in each environment that was not explained by the variance of the main effect of the genotype across environments (v G ).
The total genetic variance in environment lj was var GEI (lj)=v G +v G×E (lj) . The total genetic covariance between environments lj and l'j' was cov GEI (lj,l'j')= v G +cov G×E (lj,l'j') , where cov G×E (lj,l'j') was the covariance from the variance-covariance matrix between environments lj and l'j'. Therefore, genetic correlation between environments lj and l'j' was estimated as This correlation was used to estimate the magnitude of G×E due to ranking changes (when cor(lj,l'j') <1). Heatmaps of pair-wise genetic correlations were generated using the 'ggplot2' R package (Wickham, 2016).

Heritabilities
Broad-sense heritability was estimated in each location from the analysis of variance based on the following mixed model: where y ijk is the phenotypic value of the k th replicate of the i th individual in the j th season, μ is the mean value of the trait, g i is the random genotypic effect of individual i, s j is the fixed effect of season j, b k is the fixed effect of the block k (or replication), and e is the residual of the model. This linear mixed-effects model was fitted in R using the lme4 package (Bates et al., 2015). Broad-sense heritability of individual location clonal means (H 2 ) was then estimated as: where s 2 g is the genetic variance, s 2 e the residual variance, n is the number of seasons and r is the number of replicates per genotype.
Broad-sense heritability of clonal means across the complete MET (called 'MET broad-sense heritability') was estimated with a pool analysis across environments (i.e. location × season) using the following mixed model: where y ijk is the phenotypic value of the k th replicate of the i th individual in the j th environment, μ is the mean value of the trait, G i is the random genotypic effect of individual i, E j is the fixed effect of the environment j, GE ij is the random effect of the interaction between the i th genotype and the j th environment, and e is the error term. MET broad-sense heritability (H 2 MET ) was then estimated using the following equation: where s 2 g , s 2 ge , and s 2 e are the genotypic, genotype-byenvironment interaction, and error variance components, respectively, and e and r are the number of environments and of replicates within each environment, respectively.
Moreover, the fraction of phenotypic variation explained by the genotype, the environment and their interaction (G×E) was estimated from the mixed model fitted for the calculation of the MET heritability using the insight R package (Lüdecke et al., 2019).

QTL analyses
The R×L family was genotyped using single nucleotide polymorphism (SNP) markers from the RosBREED cherry 6K Illumina Infinium ® SNP array (Peace et al., 2012) and genetic maps have already been published (Klagges et al., 2013;Castède et al., 2014). QTL detection analyses were performed for BF and FF using the Multiple Interval Mapping (MIM) method implemented in MultiQTL V2.6 software (http://www.multiqtl.com). In each environment, the genotype means were used to perform QTL mapping. Analyses were carried out separately for 'Regina' and 'Lapins' parental maps by using the 'single QTL model' (i.e. one QTL per LG). For each location, both single-year and multi-year models were utilized. Moreover, a multi-location-multi-year analysis was performed. Both multi-year and multi-locationmulti-year analyses were carried through the multi-environment model available in MultiQTL. In all analyses, QTL significance thresholds were determined by chromosome-wide permutation tests (1000 iterations) as described in Dirlewanger et al. (2012). A wide-genome type I error of 5% was chosen and used to calculate the type I error at the chromosome level as explained in Saintagne et al. (2004). When performing multi-environment analyses, a single QTL position (in cM) and a single LOD (logarithm of the odds ratio) value are given while values of percentage of variation explained (PVE) are estimated for each environment. For ease of reading, the mean PVE value across environments is presented.

Analysis of QTL×E interactions
QTL×E analyses were performed on a selection of QTLs that, in multi-location-multi-year analyses, explained the largest part of the phenotypic variation, had the highest LOD values, and were consistently significant for the two FD stages. For each QTL, we selected the two closest flanking markers and created, for each R×L hybrid, a variable containing the genotypes of the two markers (with the code AB for heterozygous and AA for homozygous).
In order to estimate the strength of the interaction for each QTL, a step-wise approach was undertaken. Firstly, we studied each QTL independently, in single-QTL models, in order to test the significance of the QTL main effect and the significance of the interaction.
Single QTL models are an extension of the G×E (or non-QTL) model where the total genetic effect of the i th individual, g i , is decomposed into q i and x i where q i is the QTL main effect and x i is the effect of the background genotype. The general model was: y ljki = m + e lj + eb ljk + q i + qe lji + x i + xe lji + r ljki where m was the general mean across all trials, blocks within trials, seasons at trials, and individuals e lj was the fixed effect of the lj th environment (i.e. j th season at the l th trial) eb ljk was the fixed effect of the k th block at lj th environment q i was the fixed QTL main effect in the i th individual qe lji was the fixed environment specific QTL effect of the i th individual at the lj th environment x i was the background genetic effect of the i th individual with distribution N(0,Gg) where Gg was the variance-covariance matrix among total genetic effects given as Gg = I*vg where I was the identity matrix of relationships among individuals and vg was the u n k n o w n t o t a l g e n e t i c v a r i a n c e ( w h e r e * w a s t h e Kronecker product) xe lji was the random environment specific background genetic effect of the i th individual at the lj th environment with distribution N (0,Gge) where Gge was the variance-covariance matrix among environment specific total genetic effects at each of the lj th environments given as Gge = I*Vge where I was the identity matrix of relationships among individuals, and Vge was the covariance matrix of environment specific total genetic effects among environments r ljki was the residual effect for each phenotype observation with distribution N(0,R) where R was a block diagonal matrix of residual covariance among seasons for individuals at each trial, i.e. Vr l where Vr l was given as I l *Vr l where I l was an identity matrix among individuals at the l th trial and Vr l was the residual covariance matrix among seasons at the l th trial.
In single-QTL models, the variance due to other QTLs is accounted for by the background genotype.
Then, we selected the QTLs that showed either a significant main effect (p-value< 0.05), a significant interaction effect, or both, and grouped them in a "complete" model to test whether these effects remained significant when other QTLs were taken into account. Therefore, the multiple-QTL model is an extension of the single QTL model: where q qi was the fixed main effect of the q th QTL in the i th individual qe qlji was the fixed environment specific q th QTL effect of the i th individual at the lj th environment x' i was the background genetic effect of the i th individual with distribution N(0,Gg) where Gg was the variance-covariance matrix among total genetic effects given as Gg = I*vg where I was the identity matrix of relationships among individuals and vg is the unknown total genetic variance (where * was the Kronecker product) x'e lji was the random environment specific background genetic effect of the i th individual at the lj th environment with distribution N (0,Gge) where Gge was the variance-covariance matrix among environment specific total genetic effects at each of the lj th environments given as Gge = I*Vge where I was the identity matrix of relationships among individuals, and Vge was the covariance matrix of environment specific total genetic effects among environments and other terms are identical to the single-QTL model.

Development and analysis of KASP markers
In this section, we aimed to develop and/or analyze KASP markers within QTLs exhibiting the largest QTL×E interactions: QTLs on LGs 1 and 4. SNPs located in the confidence interval (CI) of the QTL on LG1 (47.4-52.3 Mb) were identified through the mapping of 'Regina' and 'Lapins' RNA-sequencing data (Maldonado et al., 2019;Vimont et al., 2019) and re-sequencing 'Lapins' data (Pinosio et al., 2020) on the 'Regina' genome sequence (Le Dantec et al., 2020) using the Integrative Genomics Viewers (IGV) software (Thorvaldsdottir et al., 2013) (http:// software.broadinstitute.org/software/igv/). Those SNPs, KASP_LG1_50.880 and KASP_LG1_52.362 (named accordingly to their physical position on the 'Regina' sweet cherry LG1 in kb, 50.880.246 and 52.362.301 in bp respectively) were then used to develop Kompetitive Allele Specific PCR (KASP) markers, as described in Branchereau et al. (2022). Moreover, markers KASP_9.936 and KASP_9.958 located on LG4 (named accordingly to their physical position on the 'Regina' sweet cherry LG4 in kb, 9.935.681 and 9.957.746 in bp, respectively), developed and validated in Branchereau et al. (2022), were used to genotype the population. Allele effect was studied in the five locations separately through analyses of variances in R software, as described in Branchereau et al. (2022). Details concerning the four KASP markers are given in the supplementary Table S3 (position, sequence, primers, Tm).  (Table S5)

G×E in the R×L population
The genotype, environment, and G×E effects explained 7.6%, 83.8% and 2.3% of the variance, respectively (Figure 3). The high proportion of variation among environmental means in Figure 3 is supported by highly significant (p< 2.2e-16) differences among the mean effect of each environment. Significant G×E interactions were observed. Estimates of pair-wise genetic correlations between environments for BF ranged from 0.50 to 0.99 and averaged 0.80 (Figure 4 and Table S8). Genetic correlations between seasons within each location were high ( For every year of study, reaction norms ( Figure 5) were not parallel (i.e. different slopes), meaning that genotypes were ranked in different orders in different locations. Genotypes responded differently to various locations, especially in 2018 and 2019. To a lesser extent, changes in individuals ranking were also observed in the five locations across different years ( Figure S6), most importantly in Maribor and Murcia. In summary, these nonparallel reaction norms showed that rank-change G×E interactions occurred in the MET, and that genotype × location interactions were stronger than genotype × year interactions.

QTL analyses for flowering date
This section will be divided into two parts. In the first one, we will present QTL analyses conducted for BF and FF in each location Variance of the fixed effect of environment and the random effects of the genotype, the genotype-by-environment interaction and residuals in flowering date (beginning of flowering).  Branchereau et al. 10.3389/fpls.2023.1142974 Frontiers in Plant Science frontiersin.org separately, with single-year and multi-year approaches. In the second part, results from a multi-locationmulti-year QTL analysis will be detailed.

QTL detection in each location
QTLs detected for BF with multi-year and single year analyses in each location are presented in Table 1 (only mean values of PVE and d values across years for each QTL are presented) and Table S9, respectively. QTLs detected for FF are presented in Table S10.
For multi-year QTL analysis at Forli, four QTLs were detected for BF on LGs R4 and R5 of 'Regina' and LGs L1 and L6 of 'Lapins' (Table 1). Across single-year analyses for this location, the QTL on LG R4 was the only locus to be significant every season and explained the largest part of the phenotypic variation (Table S9). The highest PVE value for the QTL on LG R4 (29.8%) was found in 2018. The QTL on LG L1 was only significant in 2019, where it explained 19.4% of the phenotypic variation. In 2020, an additional QTL was detected on LG L5, explaining 10.8% of the phenotypic variation.
In Murcia, the QTL on LG R4 was not significant. With the multi-year analysis, QTLs were detected on LGs R1, R7, L1, L6 and L7 (Table 1). With PVE values ranging from 16.2 to 29.9% in individual years, the QTL on LG L1 showed the largest effect in Murcia (Table S9). However, it was not significant in 2019. In 2019, a QTL on LG R3 was detected, explaining 10.9% of the phenotypic variation.
In Nimes, five QTLs were detected with the multi-year analysis, on LGs R4, R5, R8, L5 and L6 (Table 1). The QTL on LG R4 was the largest effect QTL, with an average PVE value equal to 21.5% in multi-year analysis, and single-year PVE values ranging from 14.2 to 28.5% (Table S9). The QTL on LG L6 was also significant in 2017 and 2021, when it explained 15.2 and 16.6% of the phenotypic variation, respectively.
Finally, in Toulenne, QTLs on LGs R2, R4, R5, L1, L5 and L6 were detected with the multi-year approach (Table 1). All these QTLs were detected at least for one year (Table S9). Just like in Forli,

FIGURE 5
Reaction norms for beginning of flowering (BF) across locations in each year of study. Norm of reaction of each R×L hybrid is shown by a colored line.  In summary, the QTL on LG R4 was significant in every singleyear or multi-year analysis in all locations except Murcia. In Murcia, the QTL on LG L1 was the major locus in both single-year and multiyear analyses. With multi-year analyses, the QTL on LG L1 was also significant in Forli and Toulenne. QTLs on LGs R1, R7 and L7 were only significant in Murcia, while QTLs on LGs R2, R8 and L4 were only significant in Toulenne, Nimes and Maribor, respectively.

Temperatures in each location and QTL detection
Temperatures in each site for every year of evaluation are described in supplementary figures S1 for Forli, S2 for Maribor, S3 for Murcia, S4 for Nimes, S5 for Toulenne.
Winter temperatures in 2017-2018 at Forli, (i.e. from December 2017 to the end of March 2018) were the lowest, with a long period of cold occurred in February ( Figure S1). This year corresponded to the highest estimated PVE value for the QTL on LG R4 (29.8%).
In Maribor, winter temperatures between December and February were lower in 2017 than in 2019 and 2021 ( Figure S2). The QTL on LG R4 with the highest effect was observed in 2017 with a PVE value reaching 41.7% (Tables 1, S9). In Murcia, the temperatures during the 2018/2019 winter season were not different from the other years; however, a long period of cold was observed in October/November 2018, as well as in January 2019 ( Figure S3), the only year in which a QTL on LG R3 was identified.
Temperatures during the month of January in both 2017 and 2021, were particularly low in Nimes ( Figure S4). For these two years, a QTL on LG L6 was significant.
In Toulenne, years 2016, 2017 and 2018 were relatively different from each other and year 2021 did not show any particular specificity which could be related to the detection of QTL on LG L1 ( Figure S5). LG, linkage group; L, distance from the beginning of the chromosome to the point of maximum LOD in the interval; CI, confidence interval; Physical position of flanking markers on 'Regina' v1 genome sequence in mega base pairs (Mb); LOD, logarithm of the odds ratio; PVE mean, mean value of PVE (phenotypic variance explained by the QTL in percentage of the total variation) within separate locations over several years in the multi-environment analysis; d mean, mean value of d (difference X(A) -X(B) according to the environment of evaluation, where A and B are the two homozygotes at the marker loci) in the multi-environment analysis; (+/-), the sign varies according to the environment of evaluation).

Multi-location -multi-year QTL analysis
QTLs detected with the twenty environments of the MET with the multi-environment approach are presented in Table 2 for BF and in Table S11 for FF. For BF, 12 QTLs were detected, on all LGs of 'Regina' and LGs L1, L4, L5 and L6 of 'Lapins'. Only three of them showed an overall mean PVE higher than 5%: QTLs on LGs R4 (PVE: 20.9%, LOD: 149.6), L1 (PVE: 7.2%, LOD: 39.2) and L6 (PVE: 7.3%, LOD: 42.2). QTL on LG R4 explained in average from 20.1 to 36.1% of the phenotypic variation in Forli, Maribor, Nimes and Toulenne, while it explained only 3.6% in Murcia. A significant negative correlation (-0.75, p = 0.00012) was found between the PVE value of the QTL for BF on LG R4 and the temperature in the 20 environments between October and March ( Figure S7). The opposite situation was found for the QTL on LG L1. The QTL on LG L1 explained in average 14.8% of the variation in Murcia, and from 3.1 to 7.7% in the other four locations. A correlation coefficient equal to 0.54 (p = 0.013) was found between the PVE value of the QTL on LG L1 and the temperature in the 20 environments ( Figure  S8). Finally, the QTL on LG L6 showed PVE values higher in Nimes (10.2%) and Toulenne (12.4%), compared to the other three locations (from 2.9 to 5.5%). For this QTL, no correlation with temperature data was found (correlation coefficient: -0.056, p = 0.82).
For most QTLs, genetic and physical CIs were reduced with multi-location-multi-year analysis. For instance, QTLs on LGs R4 and L5 were detected within a CI of less than 0.5 cM (8. respectively). The CI of the QTL on LG L1 (136.9-152.2 cM, 47.43-52.25 Mb) was much reduced compared to the one obtained in multi-year analyses in  and  Mb), however, it was close to the one found in .
In contrast to the results from the QTL analyses, QTLs on LGs R5, R8 and L5 were not significant (neither QTL main effect nor QTL×E interaction) for BF in the multi-environment QTL + background linear model (Table 3). The most significant QTLs for both main and interaction BF effects were QTLs on LGs R4 and L1. QTLs on LGs L6 and R3 showed significant interactions with environment, while QTL on LG R7 showed a significant main effect.
Therefore, we selected QTLs on LGs R3, R4, R7, L1 and L6 and built a complete model combining all these loci. Wald test results are presented in Table 3 (b). In the complete model, the QTL on LG R3 was not anymore significant. No differences were observed for the QTL on LG R7. QTLs on LGs R4, L1 and L6 remained the most significant loci. For the QTL on LG R4, both QTL main effect and QTL×E interaction effect increased in the complete model. LG, linkage group; L, distance from the beginning of the chromosome to the point of maximum LOD in the interval; CI, confidence interval; Physical position of flanking markers on 'Regina' v1 genome sequence in mega base pairs (Mb); LOD, logarithm of the odds ratio; PVE (phenotypic variance explained by the QTL in percentage of the total variation) mean in each location, mean value of PVE within separate location over several years in the multi-environment analysis; PVE overall mean, mean value of PVE in the multi-environment analysis; d mean, mean value of d (difference X(A) -X(B) according to the environment of evaluation, where A and B are the two homozygotes at the marker loci) in the multi-environment analysis; (+/), the sign varies according to the environment of evaluation.
Concerning the QTL on LG L1, QTL main effect increased, while the interaction with the environment slightly decreased. Finally, an important increase of the main effect of the QTL on LG L6 was observed in the complete model.

KASP markers in the QTL for flowering date on LG1
Two KASP markers (KASP_LG1_50.880 and KASP_LG1_52.362) have been developed in the CI of the QTL on LG1 and used to genotype the R×L population. Both parental cultivars 'Regina' and 'Lapins' were heterozygous for these markers; therefore, three genotypes were found in the progeny (Table 4).
For KASP_LG1_50.880, significant phenotypic differences between genotypes were found in all locations expect Maribor. In Forli, Murcia, Nimes and Toulenne, the largest BF differences were observed between hybrids with both homozygous genotypes (A:A and T:T) and between hybrids with A:T and T:T genotypes. No significant differences were observed between hybrids with A:A and A:T genotypes, except in Forli. Individuals with either A:A or A:T genotypes were flowering later than individuals with T:T genotype. The allelic effect at this marker was the highest in Murcia (up to 1.7 days of difference between A:A and T:T).
For KASP_ LG1_52.362, significant BF differences between hybrids with the three types of genotypes were found in Forli, Murcia and Toulenne. Murcia was the only location where significant differences were found between all three allelic classes. Homozygous A:A individuals were flowering 1.2, 2.7 and 1.2 days later than homozygous G:G individuals in Forli, Murcia and Toulenne, respectively.
The population was also genotyped with two KASP markers located within the QTL on LG4 (Branchereau et al., 2022). For both markers, only two genotypic classes were found in the population. In Forli, Maribor, Nimes and Toulenne, heterozygous individuals were flowering much later than homozygous individuals (p-values< 2.2e-16). The largest differences were found in Maribor: 3.5 days for KASP_9.936 and 3.7 days for KASP_9.958. In Forli, Nimes and Toulenne, differences were from 1.7 to 2.4 days for KASP_9.936 and from 1.8 to 2.5 days for KASP_9.958. On the other hand, in Murcia, very small significant differences were found (0.7 days for both markers).

Flowering date in the MET
To our knowledge, this study is the first in sweet cherry to report a MET with twenty unique location × year environments. This study shows that FD is highly dependent on the environment, and therefore is not stable across years and locations. Indeed, FD varies between years in a same location as well as between locations for a same year. An interesting observation was made in Nimes in 2017, when FD was very early. Low temperatures scored between December 2016 and the end of January 2017 may have led to a full satisfaction of the CRs, and an important increase of the temperature in February 2017 may have induced an early satisfaction of the HRs. This is a good example of climatic conditions that can induce important advances in FD, which can have dramatic consequences if the temperatures decrease again (spring frost damages) (Wenden et al., 2017). Correlations between years within a location were higher than correlations between locations. Nevertheless, estimates of broad-sense heritability were high, even at the whole MET level, in the same range as those estimated in prior studies in sweet cherry and other Prunus species (Dirlewanger et al., 2012;Castède et al., 2014;Cai et al., 2018;Calle et al., 2020;Branchereau et al., 2022). Castède et al. (2014) estimated similar heritabilities using the same R×L population planted on own roots and evaluated at Toulenne, suggesting that grafting did not have any impact on heritability in our study. Heritability values for both FD stages were higher in Toulenne than in other locations. This can be explained by a larger number of years of measurements available in this location. The G×E linear mixed model revealed strong environment and interaction effects, and the estimation of pair-wise genetic correlations between environments suggested genotypes ranking modifications across environments, especially in Murcia. This was confirmed by non-parallel reaction norms, highlighting significant rank-change G×E interactions in our MET. Jung et al. (2020) found a significant effect of G×E on floral emergence (equivalent to beginning of flowering) in apple, and highlighted as well the strong effect of the environment on this trait. The results we obtained for FD contrast with the limited effect of G×E on maturity date previously reported in sweet cherry (Hardner et al., 2019). In Hardner et al. (2019), genetic correlations for maturity date were much higher (from 0.82 to 1.0, averaged 0.95) than the one we calculated for FD (from 0.50 to 0.99, averaged 0.80). However, different environments were studied, and the environmental and climatic range has a strong impact on the estimation of G×E. Moreover, FD might be more dependent on the environment than maturity date. Indeed, FD depends on dormancy release which is closely related to temperature to fulfill CRs and HRs, and therefore requires chill accumulation followed by heat accumulation, while maturity date may be primarily dependent on heat accumulation after flowering (Gucci et al., 1991).

Major FD QTL on LGs 1 of 'Lapins' and 4 of 'Regina'
In all locations, few QTLs were detected with the single year analyses. Multi-year analyses improved detections: more QTLs were detected when combining several years of phenotypic data, and CIs were reduced. Using the 20 environments altogether further increased the power and the accuracy of the QTL detection. Many loci accounting for a very small proportion of the phenotypic variation were significant in the multi-locationmulti-year analysis. Multi-environment analysis allowed to detect a much higher number of QTLs than multi-year analyses in separate locations. For some loci (e.g. QTLs on LGs R4 and L5), the genetic position of the CI (in cM) was reduced. However, due to the low marker density of the genetic maps, it did not improve the physical position. Overall, the large number of QTLs detected confirmed the complex polygenic control of FD (Dirlewanger et al., 2012;Castède et al., 2014;Branchereau et al., 2022). Only three loci explained more than 5% of the phenotypic variation. Castède et al. (2014) studied the same R×L population planted on own roots in Toulenne and detected QTLs for FF on LGs R4, R5 and L1 in single year analyses. In the study reported here, the QTL on LG R4 was the only one to be significant every year. In the multiyear analysis combining altogether six years of measurements (period 2006-2012), QTLs were found on LGs R4, R5, R8, L1 and L2 (Castède et al., 2014). In our study, we confirmed that the QTL on LG R4 was the most stable QTL in Toulenne, but also in Forli, Maribor and Nimes. The multi-year analysis we performed in Toulenne (five years) for FF led to the detection of QTLs on LGs R2, R4, R5, L1, L5 and L6. Therefore, QTLs in common in both studies were QTLs on LG R4 and L1 (the QTL on LG R5 mapped in different chromosomal regions in both studies). These differences may be due to year, tree age, rootstock or micro-environmental effects, as well as any combination of these factors and confirm the complexity of the genetic determinism of this trait. In Branchereau et al. (2022), QTLs on LGs 1 and 4 were also detected in the 'Regina' × 'Garnet' population planted in Toulenne. QTL on LG1 was found in 'Garnet' cultivar but was not stable across years, likewise the QTL on LG1 in 'Lapins' cultivar in Toulenne. On the other hand, the QTL on LG4 of 'Regina' was very stable, detected every year of study (i.e. over ten years) (Branchereau et al., 2022), as observed in this study in Forli, Maribor, Nimes and Toulenne. This QTL was detected in a physical CI of the same range in both studies.
Conducting QTL mapping in each location separately and then comparing the results allowed us to discover that the QTL on LG R4 was found in all locations except Murcia, and that, in this location, the major QTL was located on LG L1. The QTL on LG L1 was also detected across some years in Forli and Toulenne. This LG1 QTL has been identified in sweet cherry cultivars 'Cristobalina', 'Garnet' and 'Lapins' in the chromosomal region of the DAM genes (Dirlewanger et al., 2012;Castède et al., 2014;Calle et al., 2020;Branchereau et al., 2022). 'Cristobalina' is an extra-early blooming cultivar with very low CRs, and a recent study revealed that it carries structural mutations in the DAM genes region that might be responsible of this phenotype (Calle et al., 2021). Therefore, the QTL on LG1, covering DAM genes, seems to be related to the CRs.
Nevertheless, the situation is probably highly complex. Indeed, Castède et al. (2014) reported CR QTLs on LGs R1 and G1 by working with a population derived from the cross between 'Regina' and 'Garnet' but the corresponding peaks mapped in a clearly different genomic position as the one carrying the DAM genes. The same result was observed a few years later by conducting the same type of analysis on population 'Regina' × 'Lapins', in which CR QTLs were detected on LG L1 but again, in an upstream chromosomal region (Quero-Garcia, comm. pers.). On the other hand, the LG4 QTL might be more related to HRs but also to CRs since Castède et al. (2014) found clear co-localizations of bloom date, CR and HR QTLs on LG R4. We calculated significant positive and negative correlations between the temperature and the proportion of the phenotypic variance explained by QTLs on LGs L1 and R4, respectively. This suggests that the QTL on LG L1 plays a major role in warm region environments (where HRs are easily fulfilled but not CRs), while QTL on LG R4 is more significant in colder regions (where CRs are easily fulfilled but not HRs). This result might contribute to future experimental designs aimed at elucidating the complex regulation of genes involved in FD, underlying these QTL regions, in particular its interaction with temperature.
In our MET, QTLs on LGs 1 and 4 may be described as 'conditionally neutral' QTLs, because they are detected in only specific environments (El-Soda et al., 2014). If a QTL is detected in some environments but not in others, it implies QTL×E interactions.

QTL×E interactions and MAS
In both single-QTL and multiple-QTL models, loci with most significant and largest QTL×E interactions were QTLs on LGs R4 and L1. This is in accordance with the QTL detections that revealed that both loci were large-effect environment-specific QTLs.
In a context of MAS, both QTL main effect and QTL×E effect should help in the selection of hybrids particularly adapted to specific environments. For instance, a breeder aiming at the release of new cultivars for cold production areas, will put more weight on the QTL of LG4, if using cultivar 'Regina' as a parent, by selecting 'late flowering' alleles in order to avoid the risk of frost damages. On the opposite, a cultivar adapted to regions characterized by warm winters with a lack of chill, which will derive from 'Lapins', will need to inherit 'early flowering' alleles of LG1 QTL, most likely associated to low CRs. Finally, for intermediate environments such as the ones represented in our study by the sites of Toulenne, Forli or Nimes, it might be advisable to combine 'early flowering' alleles from the LG1 QTL of 'Lapins' with 'late flowering' alleles of the LG4 QTL of 'Regina', by trying to combine in a single hybrid 'sufficiently' low CRs with 'sufficiently' high HRs.
These hypotheses/strategies were supported by the analysis of four KASP markers located within the QTLs on LGs 1 and 4. Indeed, the effect of the LG1 KASP markers was larger in Murcia (the warmest environment of the MET) than in other locations, while LG4 KASP markers (Branchereau et al., 2022) played a major role in Maribor, Toulenne, Nimes and Forli. Additionally, two markers developed by Calle et al. (2021) and located within the DAM genes on the LG1 might be useful. Both markers were developed from the extra-early cultivar 'Cristobalina' for the detection of structural mutations (within the DAM genes region) associated to early FD and low CRs in sweet cherry (Calle et al., 2021). More recently, in peach, KASP markers were developed in the region spanning from 43.58 to 43.78 Mb on the chromosome 1 (Pp01) and validated to predict CRs (Demirel et al., 2023). These markers, located near the DAM genes, are located 1 Mb upstream from those that we developed within the LG1 QTL (orthologous positions of KASP_LG1_50.880 and KASP_LG1_52.362 on the peach chromosome 1 are 44.60 Mb and 45.78 Mb, respectively). However, the strong association we found for KASP_LG1_50.880 and KASP_LG1_52.362 with FD in Murcia shows that these markers could be useful for selection in warm environments.
Therefore, all these markers should contribute to establish a complete MAS strategy for FD, their choice depending of the climatic conditions of the place where cherry trees will be planted. This study demonstrates that MAS should be performed with different markers if the climate is warm or cold to select well adapted genotypes to a specific region. In climates with warm winters, genotypes with alleles responsible of early flowering at markers on LG1 should be selected, while in climates with cold winter, genotypes with alleles responsible of late flowering on LG4 should be selected to avoid frost damage. In intermediate environments, selection should be done on both loci by screening genotypes with 'early flowering' alleles on LG1 and 'late flowering' alleles on LG4.

Conclusion
To our knowledge, this study is the first in sweet cherry to perform QTL analyses in a complex MET and estimating G×E and QTL×E interactions for FD. We showed that FD is highly dependent on the environment with important inter-annual and inter-location variations. Differences of individuals ranking between environments were the major source of G×E detected in this study. QTL×E plays a major role in adaptation to environment changes. Our study revealed that two major FD loci in sweet cherry, located on LGs 1 and 4, exhibited strong QTL×E. Therefore, this study provides relevant information for the choice of stable QTLs in specific environments in order to target them in MAS. Molecular markers have been developed in both loci, and therefore could be used simultaneously to start a complete MAS strategy for FD and develop new cultivars well adapted to their cultivation area. Molecular breeding based on these markers could be undertaken to select genotypes for specific climatic conditions. This study focused on FD, a key trait in sweet cherry breeding, but other traits related to fruit quality could be studied as well. Therefore, this unique sweet cherry MET paves the way for molecular breeding strategies.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.