Genetic Dissection of Bloom Time in Low Chilling Sweet Cherry (Prunus avium L.) Using a Multi-Family QTL Approach

Bloom time in sweet cherry (Prunus avium L.) is a highly heritable trait that varies between genotypes and depends on the environmental conditions. Bud-break occurs after chill and heat requirements of each genotype are fulfilled, and dormancy is released. Bloom time is a critical trait for fruit production as matching cultivar adaptation to the growing area is essential for adequate fruit set. Additionally, low chilling cultivars are of interest to extend sweet cherry production to warmer regions, and for the crop adaptation to increasing winter and spring temperatures. The aim of this work is to investigate the genetic control of this trait by analyzing multiple families derived from the low chilling and extra-early flowering local Spanish cultivar ‘Cristobalina’ and other cultivars with higher chilling requirements and medium to late bloom times. Bloom time evaluation in six related sweet cherry populations confirmed a high heritability of this trait, and skewed distribution toward late flowering, revealing possible dominance of the late bloom alleles. SNP genotyping of the six populations (n = 406) resulted in a consensus map of 1269 SNPs. Quantitative trait loci (QTL) analysis using the Bayesian approach implemented by FlexQTL™ software revealed two major QTLs on linkage groups 1 and 2 (qP-BT1.1m and qP-BT2.1m) that explained 47.6% of the phenotypic variation. The QTL on linkage group 1 was mapped to a 0.26 Mbp region that overlaps with the DORMANCY ASSOCIATED MADS-BOX (DAM) genes. This finding is consistent with peach results that indicate that these genes are major determinants of chilling requirement in Prunus. Haplotype analysis of the linkage group 1 and 2 QTL regions showed that ‘Cristobalina’ was the only cultivar tested that contributed early bloom time alleles for these two QTLs. This work contributes to knowledge of the genetic control of chilling requirement and bloom date and will enable marker-assisted selection for low chilling in sweet cherry breeding programs.


INTRODUCTION
Bloom time (BT) is an important horticultural trait in temperate fruit tree species like sweet cherry (Prunus avium L.). Cultivar adaptation to climatic conditions in the growing area is essential for flower production and fruit set. Early blooming cultivars are susceptible to spring frost damage in cold regions (Luedeling, 2012), while late blooming cultivars can exhibit irregular floral development and low fruit set due to warm temperatures during the flowering period (Mahmood et al., 2000;Atkinson et al., 2013). The biological control of BT is complex and is known to depend on environmental signals during the winter and spring seasons Fadón and Rodrigo, 2018). Fruit trees like sweet cherry require a period of chilling temperatures followed by a period of warm temperatures to induce blooming (Lang et al., 1987). In Prunus species, several studies indicate that BT is more dependent on chilling than on heat requirement and that there is large variation in these requirements among individuals of the same species (Campoy et al., 2011;Okie and Blackburn, 2011;Sánchez-Pérez et al., 2012;Castède et al., 2014).
Several genetic studies in Prunus have focused on understanding the genetic control of chilling (CR) and heat requirements (HR) contributing to the differences in BT (reviewed in Abbott et al., 2015). BT in Prunus is a quantitative trait with high heritability (Anderson and Seeley, 1993;Dirlewanger et al., 2012;Castède et al., 2014), and genetic approaches have led to the identification of quantitative trait loci (QTLs) associated with this trait. BT QTLs have been identified on all Prunus linkage groups (LGs) (reviewed in Salazar et al., 2014;Abbott et al., 2015), but major QTLs have been found on LG1 (Fan et al., 2010;Zhebentyayeva et al., 2014;Bielenberg et al., 2015) and LG4 (Sánchez-Pérez et al., 2012;Dirlewanger et al., 2012;Castède et al., 2014) in all the Prunus crop species evaluated to date. In sweet and sour (Prunus cerasus) cherries, several QTLs have been identified for BT and CR. In sweet cherry, Dirlewanger et al. (2012)

mapped two major BT QTLs on
LGs 1 and 4 and three minor QTLs on LGs 5, 6 and 7. Castède et al. (2014) using three to six years data and two F 1 populations identified BT and CR QTLs in all LGs, with a major and stable QTL for both traits overlapping on LG4. Castède et al. (2014) also detected minor QTLs for both CR and BT on LGs 1 and 7, highlighting the influence of CR in BT in this species. In sour cherry, Wang et al. (2000) investigated BT QTL using an F 1 population and 3-year data, and two major QTLs were identified on LGs 1 and 2. Another QTL study in sour cherry revealed four QTLs for BT on LGs 1, 2, 4, and 5; the most significant allele for LG4 QTL was associated with almost three days bloom delay .
Candidate genes have been reported for the Prunus BT and CR QTL that maps to LG1 (Yamane et al., 2011;Zhebentyayeva et al., 2014;Castède et al., 2015). In peach, a tandem set of six DORMANCY ASSOCIATED MADS-BOX (DAM) genes have been identified in this region (Zhebentyayeva et al., 2014;Romeu et al., 2014;Bielenberg et al., 2015). Studies of these genes was facilitated by the study of a peach non-dormancy mutant termed evergrowing peach mutant (evg) that has a deletion within this QTL region (Bielenberg et al., 2008). Expression analyses in peach reported that DAM5 and DAM6 are not expressed under chilling temperatures (Jimenez et al., 2010) whereas the expression of DAM4 and DAM6 are activated by short photoperiods (Li et al., 2009) suggesting that DAM5 and DAM6 are the main genes underlying this Prunus LG1 CR QTL (Yamane et al., 2011). For the major Prunus BT QTL located on LG4, genes related to temperature sensing (ARP4, EMF2, NUA, and PIE1) and the gibberellin pathway (GA2ox and KS genes) have been proposed as the most promising candidates to control BT (Dirlewanger et al., 2012;Castède et al., 2015).
Most BT QTL studies in Prunus have been done using a single bi-parental population. This strategy is limited because only alleles present and segregating in the two parental cultivars can be detected (Bink et al., 2014). However, knowledge of the effects of these alleles in different genetic backgrounds and other loci not segregating in the bi-parental cross are needed to fully implement marker-assisted selection (MAS). The development of QTL mapping approaches based on multi-parental populations allow the identification a larger number of QTLs and QTL alleles improving the application of these results in MAS for a larger number of genetic backgrounds (Pauly et al., 2012). The Bayesian QTL mapping approach implemented by FlexQTL ™ (Bink et al., 2008 andBink et al., 2014) allows analyzing multiple pedigreelinked progenies at the same time; reducing the limitations derived from QTL analyses using single populations. This approach has been successfully used in recent years for QTL analyses of different traits in Rosaceae species, such as sweet cherry (Rosyara et al., 2013), apple (Bink et al., 2014;Guan et al., 2015;Allard et al., 2016;Di Guardo et al., 2017;Howard et al., 2018), peach (Fresnedo-Ramírez et al., 2015;Fresnedo-Ramírez et al., 2016;Hernández Mora et al., 2017), and strawberry (Roach et al., 2016;Mangandi et al., 2017;Anciro et al., 2018).
Furthermore, previous QTL analyses in sweet cherry used cross-pollinated F 1 populations. Self-pollination is often not possible in sweet cherry due to the gametophytic selfincompatibility system operating in this species (Herrero et al., 2017). However, natural self-compatible sweet cherry mutants like 'Cristobalina' (Wünsch and Hormaza, 2004) or other selfcompatible sweet cherry accessions can be used to generate F 2 populations which can be used for genetic mapping studies. The self-compatible local cultivar 'Cristobalina' and its selfcompatible descendant, the selection 'BC-8', were used to develop two self-pollinated populations for genetic analysis. Genetic maps of these sweet cherry populations were constructed, and their level of homozygosity was previously reported (Calle et al., 2018). These were the first F 2 populations used for genetic map construction in this species and are now available for QTL analysis. The Spanish cultivar 'Cristobalina' comes from the Mediterranean area, and in addition to being self-compatible (Wünsch and Hormaza, 2004;Cachi and Wünsch, 2011;Ono et al., 2018), it has a low CRs (< 550 h; Tabuenca, 1983;Alburquerque et al., 2008) and extra early flowering and maturity dates. These characteristics make 'Cristobalina' an interesting breeding cultivar. Cultivars with low CRs often show early flowering (Alburquerque et al., 2008) and this low chilling requirement is of high interest for extending cherry growing to regions with warmer winters, and in the current context of global warming as a source of adaptation to low chilling.
The objective of this work was to investigate the genetic basis of BT in different genetic backgrounds, including in the low chilling 'Cristobalina'. To achieve this objective, six related sweet cherry populations that descend from 'Cristobalina' and other sweet cherry cultivars with mid and late BTs were used. Four years of BT data from these families was used for QTL analysis using a Bayesian approach implemented in FlexQTL™. Two self-pollination populations were also used to investigate the genetic effects of the QTL alleles. The results obtained broaden the understanding of the genetic control of CR and BT in this species and will allow the implementation for MAS of these traits from this and related plant material.

Bloom Time Phenotyping
BT was evaluated for each progeny and the parental cultivars during 4 years (2015 to 2018). During the flowering season, blooming stage was evaluated 3 days per week (every 2-3 days) and BT was recorded when approximately 75% of the floral buds reached stage F (full bloom; Baggiolini, 1980). BT data were expressed as Julian Days (JD; days from January 1 st ). The best linear unbiased prediction (BLUP) value among years was calculated using the lme4 package of R 3.4.1 software (Bates et al., 2015;R Core Team, 2017).
Mean, minimum and maximum values and the distribution of BTs were estimated in each population per year and for the BLUP values. Mean BTs were then compared between families using analysis of variance (ANOVA) and Tukey post-hoc test (p < 0.05) for BLUP values. Pearson correlation coefficient of BTs between years and BLUP values were also estimated. Statistical analyses were carried out using SPSS statistics v21.0.0 software (IBM, Chicago, IL, USA).

Genotyping and QTL Analysis
Genotypes of 417 individuals that include the six populations described above, their parental and ancestor cultivars ( Figure 1) were used for QTL analysis. All these plant material had been previously genotyped using the RosBREED Illumina cherry 6K SNP array v1 (Calle et al., 2018). For QTL mapping using FlexQTL ™ (Bink et al., 2014), the genotyped SNPs of all plant materials were previously filtered. SNPs monomorphic in all populations, that had null alleles, with MAF < 0.05, with more than 5% of missing data, and showing errors in various genotypes were discarded. The selected markers were further checked by analyzing their genetic segregation using FlexQTL ™ and by visual inspection when a double recombination event occurred within an interval smaller than 10 cM. A consensus genetic map from the selected SNPs was constructed. Those SNPs previously mapped in V×C (Calle et al., 2018) were assigned their genetic position. Those SNPs not previously mapped in the V×C map were integrated in the map by using their physical position on the peach genome v2.0.a1 (Verde et al., 2017).
QTL mapping for BT, each year and for BLUP values, using all the plant material, was carried out using a Bayesian multiple QTL model implemented in FlexQTL ™ (Bink et al., 2014). In this work, only additive effects with normal prior distribution were considered. The models were set with a prior distribution of number of QTLs of 1 and 3 in order to assess sensitivity of posterior inference to the prior assumptions. Markov chain Monte Carlo (MCMC) simulation with minimum of 500,000 iterations for each prior number of QTL were performed until at least 100 effective chain samples (ECS) were reached for overall mean (m), the residual variance (s 2 e ), the number of QTLs (N QTL ), and the variance of QTLs ( V QTL) (Gelman and Rubin, 1992;Sorensen and Gianola, 2002). A total of 1,000 samples (500,000 iterations with thinning value of 500) were stored for further posterior inferences. The inference in the number of QTLs was estimated using twice the natural log of Bayes factors (2lnBF) (Kass and Raftery, 1995) obtained after a pairwise comparison of models differing by one QTL from each other. Values of 2lnBF higher than 2, 5, and 10 indicate positive, strong and decisive QTL evidences, respectively. Only QTLs with strong and decisive evidences were considered in this work. The percentage of variation explained (PVE) by each QTL was estimated as [(wVAt1/PV) × 100], where wVAt1 is the weighted variance and PV is the total phenotype variation (Mangandi et al., 2017). The genomic breeding value (GBV) for each individual and parent was predicted using QTL genotype probabilities, intensity and effect size (Bink et al., 2014).

Haplotype Analysis
SNP haplotypes of the two most significant QTLs detected with an average 2lnBF higher than 10, were constructed for the parental cultivars and theirs ancestors. The haplotypes were designed to span the confidence interval with 2lnBF > 10 for these QTLs using BLUP model. The SNP haplotypes were estimated using SNP phase estimation of FlexQTL ™ for all the cultivars, except for 'Bing' and 'Napoleon' that were estimated manually (for QTL on LG2) based on pedigree information and the availability of previously phased haplotypes (Cai et al., 2017). SNP haplotypes were also confirmed based on segregation. Mean phenotypic values of each of these QTL haplotypes and their combined effects were estimated in each segregating class of each population. Individual progenies with recombination events within these QTL regions were excluded from the analysis. For mean comparison of the phenotypic values within each population, Kruskal-Wallis, two-tailed Student's t test and Tukey test (p < 0.01) post-hoc analysis were used. All statistical analysis were done using SPSS statistics v21.0.0 software (IBM, Chicago, IL, USA).

Bloom Time Phenotyping, Segregation, and Heritability
BT was evaluated in the parental cultivars and populations during 4 consecutive years (Figure 1; Supplementary Table S1; Supplementary Figure S1). The parental cultivar 'BC-8' was phenotyped only the first year as the tree was in poor health in subsequent years. 'Cristobalina' was the earliest parental cultivar to bloom in all of the years [BLUP value: 69 JDs; Figure 1 and Figure 2]. 'Ambrunés' 'BC-8' and 'Brooks' showed midseason flowering ( Figure 1 and 2), while 'Vic' and 'Lambert' exhibited late blooming with BLUP values of 95 and 97 JDs, respectively ( Figure 1 and 2). CR were fulfilled for 'Cristobalina' (550 h; Tabuenca, 1983) between mid-December to the first week of January during the 4 years of analysis (Supplementary Figure  S2). The rest of the parental cultivars, 'Ambrunés' 'Brooks,' and 'Lambert' all had higher CR (900 to 1100 h; Tabuenca, 1983;Alburquerque et al., 2008). During the 4 years evaluated, the CR of these three cultivars were not fulfilled until mid-January to late February (Supplementary Figure S2).
The same blooming order (extra-early, mid, and late blooming) was observed for the parental cultivars each year, but differences in the BTs were observed between years (Supplementary Table S1; Supplementary Figure S1). BT was earliest in 2017 for the mid and late cultivars, which bloomed 16 to 17 days earlier than the average date of the rest of the years. However, for 'Cristobalina' the earliest bloom period occurred in 2016 (Supplementary Table S1; Supplementary Figure S1), while the latest bloom period for all parental cultivars occurred in 2018 (Supplementary Table S1; Supplementary Figure S1). In 2016, fulfilment of the CR and HR for 'Cristobalina' occurred early resulting in an early bloom. However, this was followed by a period of cold temperatures that delayed the flowering of the rest of the cultivars extending the blooming season (Supplementary Figure S2). In 2017, a high accumulation of chill hours during the early winter followed by a period of warmer temperatures in the beginning of February resulted in an earlier bloom for all the cultivars and a shorter blooming period (Supplementary Figure  S2). In 2015 and 2018, although large amounts of chill were accumulated during the early winter, cold temperatures in February delayed bloom. Years 2017 and 2018 were colder and CR were fulfilled earlier in the year, but BT was earlier in 2016.
In the populations, different numbers of offspring were phenotyped each year, ranging from 258 (64%) in 2015 to 367 (90%) in 2018 (Supplementary Table S1). Only individuals (N = 360) with phenotypic data from two or more years were used to estimate BLUP values. The bloom period varied between years and populations, from 7 to 24 days. On average A×C, B×C, and L×C showed shorter bloom periods (10 to 12 days) than B×C F2, C×C and V×C (16 to 18 days) (Supplementary Table S1; Supplementary Figure S1). 'Cristobalina' self-pollination (C×C) was the earliest population to bloom (Figure 1 and   Table S1; Supplementary Figure S1). It was especially noticeable in 2017 when warm temperatures resulted in earlier flowering for all genotypes and the shortest bloom period (18 days) of all years. In 2016, warm temperatures in mid-February, resulted in very early bloom of some individuals from C×C population, but a cold period later on delayed the flowering of the rest of the population, resulting in the largest bloom period in all years evaluated.
BT distribution varied between populations and years. Only the smallest populations (B×C and L×C) fitted a normal distribution for all evaluated years and for BLUP values (Shapiro-Wilk test; Prob < W: 0.083-0.263). Populations B×C F2 and A×C fitted a normal distribution only two of the years (Shapiro-Wilk test; Prob < W: 0.085-0.664), whereas the remaining two populations (C×C and V×C), which are the largest, did not fit normal distributions in any of the years. BT of all the progenies together also did not show normal distribution for any of the years or BLUP values ( Figure 2). In all cases, skewed distributions towards medium and late BT were observed. In fact, only some C×C offspring were extra-early blooming, and only some B×C F2 offspring were early blooming. The rest of the plant materials were medium to late blooming, even though all populations (except B×C F2) had 'Cristobalina' (extra-early blooming) as one of the parental cultivar. Transgressive segregation toward late blooming was also observed in the 4 years and for BLUP values for all populations except L×C. On the other side, transgressive segregation towards early blooming was only observed in the self-pollination populations (C×C and B×C F2; Figure 2; Supplementary Table S1; Supplementary Figure S1).
Highly significant (p < 0.01) and positive correlations were observed for BTs between years and BLUP value (r = 0.897 to 0.966; Supplementary Table S2). BT broad-sense heritability (H 2 ) estimated for all populations together was 0.97, and for each population individually, the H 2 was also high, with values ranging from 0.85 to 0.96 (Supplementary Table S1).

QTL Analysis
Quality filtering of the SNP markers resulted in 1,269 (22.3%) SNPs selected for map construction (Supplementary Tables S3 and S4) and QTL analyses. These selected SNPs covered a total genetic length of 721.98 cM with an average marker density of 0.57 cM (1 SNP per 176 kb) and a maximum gap between markers of 10.95 cM (1.43 Mbps) located on LG 7 (Supplementary Table S4).
QTL analysis using BLUP values of the 4 years in all the populations revealed three significant BT QTLs located on LGs 1, 2 and 4, of which qP-BT1.1 m and qP-BT2.1 m were found with decisive evidences (2lnBF > 10) (Table 1 and Figure 3). QTL on LG1, qP-BT1.1 m , explained the largest percentage of phenotypic variation (PVE; 32.4%) and was associated with a mean additive effect of 7.4 days (Table 1 and Figure 3). The other decisive QTL, qP-BT2.1 m , had a PVE of 15.2% and an additive effect of 5.9 days ( Table 1). Also a QTL on LG4 (qP-BT4.1 m ) with strong evidence was identified, although it had a lower additive effect (3.6 days) and PVE (6.0%) than the other two major QTLs (Table 1 and Figure 3).
In the QTL analyses for individual years, the two major QTLs detected on LGs 1 and 2 in the 4-year analysis were also detected every single year with strong or decisive evidence ( Table 1). In these years, the PVE of qP-BT1.1 m and qP-BT2.1 m ranged from 14.2% to 60.9%, and from 15.0% to 23.0%, respectively. The lowest and highest PVE for qP-BT1.1 m were observed in 2016 and 2018, respectively, while the mean additive effect ranged from 7 to 10 days between 2016 and 2015, respectively. For qP-BT2.1 m , the lowest and highest PVE were observed in 2015 and 2016, respectively, and the mean additive effect ranged from 3 to 11 days in 2017 and 2016, respectively. In 2016, the PVE and mean additive effect exhibited by qP-BT2.1 m was larger than that observed by qP-BT1.1 m ( Table 1). In other years, like 2017 and 2018, in which the effects of these two QTLs were lower, additional QTLs with minor effects were detected.

Estimation of QTL Genotypes and Genomic Breeding Value
QTL genotype estimation was carried out for QTL regions with either strong or decisive evidence using BLUP value (Figure 4) for the parental cultivars and the ancestors in the collection (Figure 4). For the major QTLs on LG1 (qP-BT1.1 m ) and LG2 (qP-BT2.1 m ), 'Cristobalina' was predicted to be homozygous for alleles associated with early bloom (qq = low phenotype value) for the LG1 QTL and predicted to be heterozygous (Qq) for the LG2 QTL (Figure 4). 'BC-8' an offspring from the cross of 'Brooks' and 'Cristobalina' was heterozygous for the early bloom allele for the LG1 QTL (qP-BT1.1 m ). No other parental cultivar was predicted to have early BT alleles for these two QTL, instead the remaining parental cultivars were predicted to be homozygous (QQ) for LG1 and LG2 QTL alleles for later BT (Figure 4). This indicates that 'Cristobalina' is the only cultivar that contributed early BT alleles for the major QTLs qP-BT1.1 m and qP-BT2.1 m of all the plant materials. For the QTL on LG4 (qP-BT4.1 m ), only 'Rainier' and its offspring 'Brooks' were predicted to be heterozygous for early BT alleles, while the rest, including 'Cristobalina' were predicted to be homozygous (QQ) for late bloom alleles.
Differences in the predicted genotypes were used to estimate genome breeding value (GBV) of parents and ancestors ( Figure 4). Differences of as much as 27 JDs (almost 1 month) were observed between the GBV of the earliest and latest blooming cultivars (Figure 4). All but 'Cristobalina' had GBV associated with delayed flowering. 'Cristobalina' had the lowest GVB, due to the relative abundance of alleles predicted to result in earlier flowering of 8.9 JDs. In contrast, 'Vic' and 'Lambert' cultivars had the GBV most associated with late flowering (17.0 and 17.9 JDs, respectively; Figure 4).

QTL Haplotype and Genotype Analysis
The haplotypes (alleles) of the two most significant QTLs detected, those on LG1 (qP-BT1.1 m ) and LG2 (qP-BT2.1 m ), were constructed for the parental cultivars and their ancestors (Supplementary Table S5) to investigate their phenotypic effects ( Figure 5). Haplotypes of qP-BT1.1 m and qP-BT2.1 m were constructed with 4 and 14 SNPs, respectively (Supplementary  Table S5). H1-c was only found in 'Cristobalina' (homozygous for this allele) and in the selection 'BC-8' (heterozygous). The remaining two haplotypes, H1-a and -b, were found in the rest of cultivars (Supplementary Table S5).
For the QTL qP-BT2.1 m , 10 haplotypes were identified (H2-a to H2-j). All the parental cultivars had unique haplotypes except 'Brooks' and 'Lambert' that shared H2-d, and 'BC-8' that shares haplotypes with its parental cultivars, 'Brooks' and 'Cristobalina' (Supplementary Table S5). The estimation of the mean phenotypic values of these QTL haplotypes in the F 2 populations (C×C and B×C F2) revealed that for the QTL qP-BT1.1 m , those individuals that were homozygous for H1-c, like 'Cristobalina' (cc), showed the earliest BTs ( Figure 5). As in C×C, this QTL was not segregating and all the progeny were also 'cc,' with a mean BT of 75 JDs. In B×C F2, this QTL was segregating but segregation distortion was observed as no 'aa' individuals were identified. In the two remaining segregating classes, individuals with the 'cc' genotype showed a mean difference of almost 7 JDs earlier blooming that those with 'ac' genotype ( Figure 5).
For the QTL qP-BT2.1 m , H2-f was associated with early flowering. In C×C ('Cristobalina'; 'ef'), this QTL segregated in three classes, with offspring that were 'ef' and 'ff' flowering on average 7 JDs earlier that those that were 'ee' ( Figure 5). As no significant differences were observed between these two classes, H2-f appeared to be dominant to H2-e ( Figure 5). B×C F2 also segregated in three classes for this QTL, but no significant differences were observed among them ( Figure 5). Since H2-f was not inherited in 'BC-8' (ce), the effect of this haplotype could not be investigated in this population.
The interaction of these two major QTLs showed that those individuals homozygous for H1-c were the earliest to bloom for both populations. Within these, those that also had H2-f showed the earliest BT ( Figure 5).
In the F 1 populations (Supplementary Figure S3), for qP-BT1.1 m , "ac" genotypes were always earlier blooming (approx. 2 to 3 days that those that are "bc"; Supplementary Figure S3), indicating that H1-a was associated with earlier BT compared to H1-b. For qP-BT2.1 m , genotypes with H2-f showed earlier BT (1 to 7 JDs) compared to individuals without it. In addition, segregation distortion against H2-f was observed in all the populations, being most evident in B×C, as none of the progeny had this haplotype (Supplementary Figure S3). A genotype interaction of both QTLs in the F 1 populations was identified as individuals heterozygous for H1-a/H2-f, showed earlier BT than those individuals with other genotype combinations (Supplementary Figure S3).

DISCUSSION
'Cristobalina' the earliest blooming parental cultivar, has a low CR (176-550 h under 7°C;Tabuenca, 1983;Alburquerque et al., 2008); which is consistent with its native origin close to the Mediterranean coast in eastern Spain (Herrero, 1964). Plant's CR are typically correlated to the climate in the area of origin . In our experimental location, which FIGURE 4 | Posterior estimates of parental quantitative trait loci (QTL) genotype probabilities in QTL regions with strong and decisive QTL evidences (2lnBF > 5) for best linear unbiased prediction (BLUP) values. Red, green and blue colors represent positive evidence for QTL genotypes QQ, Qq, and qq genotypes, respectevely. "Q" and "q" denote alleles with high and low phenotype values, respectively. Gray colors indicate unclear genotype estimation. Genome breeding value (GBV) for each cultivar is indicated at right. experiences a higher chilling accumulation, it is likely that the early blooming of 'Cristobalina' is due to its low CR, as earlier blooming has been observed in cultivars with lower CR (Alburquerque et al., 2008;Castède et al., 2014). 'Cristobalina' is also self-compatible (Wünsch and Hormaza, 2004) due to a single mutation affecting pollen tube growth (Cachi and Wünsch, 2011;Ono et al., 2018). Natural self-compatible mutations are rare in sweet cherry; however, this mutation may be especially beneficial in this low chill cultivar because mating partners with overlapping flowering times would be scarce (Cachi et al., 2014).
The same BT (JD) order, early, mid-and late bloom, of the parental cultivars in the 4 years, independently of the year temperatures, confirms the genetic determination of this trait. As previously demonstrated, BT in cherries is a quantitative trait with high heritability (Dirlewanger et al., 2012;Castède et al., 2014;Cai et al., 2018). High heritability values for this trait were also observed in this work for the 4 years (0.85 to 0.96), and these values are in the same range as those estimated previously for sweet and sour cherry (0.88 to 0.96; Wang et al., 2000;Dirlewanger et al., 2012;Castède et al., 2014;Cai et al., 2018;Piaskowski et al., 2018). However, BT differences between years are highly dependent on environmental conditions and how these conditions impact CR fulfillment. For example, the coldest winters did not result in the earliest BT.
Within the populations, only individuals from F 2 populations (C×C and B×C F2) showed transgressive segregation towards early blooming, whereas F 1 populations showed transgressive segregation and skewed distribution towards late flowering, revealing possible dominance of the late bloom alleles in this plant material. In sweet cherry, skewed segregation in F 1 populations towards high CR, but not late bloom, was also observed (Castède et al., 2014). This was also the case in almond, where Late bloom (Lb) is dominant (Ballester et al., 2001), and in Japanese plum and apricot F 1 populations from low CR cultivars (Campoy et al., 2011;Salazar et al., 2016;Kitamura et al., 2018). However, transgressive segregation toward both early and late blooming was observed in peach F 2 populations (Fan et al., 2010;Bielenberg et al., 2015). The effect of (recessive) alleles in the homozygous state will be possible to detect in F 2 populations (like C×C and B×C F2), which may explain why transgressive segregation towards early bloom in our plant material is only observed in the F 2 populations. The extended bloom period observed in the larger populations (C×C, V×C) may also have resulted from additional climatic variation experienced during this longer BT duration. This effect was also observed by Castède et al. (2014).
Two Major BT QTLs on LG1 (qP-BT1.1 m ) and LG2 (qP-BT2.1 m ) The major QTL identified on LG1 (qP-BT1.1 m ) has been previously detected in sweet cherry (Dirlewanger et al., 2012;Castède et al., 2014). However, in these works, the variation explained and additive effect of this QTL were lower than observed herein. Dirlewanger et al. (2012) first identified this LG1 QTL in 'Lapins' with a PVE ranging from 9.3 to 17.5% and an additive effect of 2.5 days. Castède et al. (2014) reported the same QTL region for two of three years for a 'Regina' × 'Garnet' population, with similar additive effect (1.4 days) and mean PVE (8%). In our work, this QTL represented 32.4 of PVE and has an additive effect of 7.4 JDs. These results indicate that BT of our plant material, in our environmental conditions, was determined by this QTL in a larger proportion than in earlier works in sweet cherry.
A CR QTL overlapping with this LG1 BT QTL was also identified in sweet cherry (Castède et al., 2014), confirming the correlation between both traits and the relevance of CR for BT in the species. This BT and CR QTL on LG1 has also been described in other Prunus species like peach, apricot and almond, as a main QTL controlling these traits (Olukolu et al., 2009;Fan et al., 2010;Dirlewanger et al., 2012;Sánchez-Pérez et al., 2012;Salazar et al., 2013;Zhebentyayeva et al., 2014;Romeu et al., 2014;Bielenberg et al., 2015). In our work, the high significance and effect of this QTL indicates that, in our conditions, the BT of the plant material analyzed is probably more dependent on CR than other materials analyzed in different environments.
Candidate genes of the LG1 QTL have been described in peach and sweet cherry (Bielenberg et al., 2008;Fan et al., 2010;Castède et al., 2015). The position of this QTL overlaps with the region where six DORMANCY ASSOCIATED MADS-box (DAM1-6) genes have been identified as major genes controlling CR and BT in peach (Bielenberg et al., 2008;Fan et al., 2010). In the evergrowing peach mutant, that lacks response to winter cold, four of these genes are deleted and the other two are not expressed (Bielenberg et al., 2008). Castède et al. (2015) mapped two of these DAM genes (DAM5 and DAM6) within the interval of this QTL in 'Lapins' sweet cherry. It is likely that the 'Cristobalina' alleles of these genes are contributing to low CR and early blooming, and the large effect of this QTL in the plant material analyzed. 'Cristobalina' contributed H1-c for this QTL, which was associated with earlier flowering, as BT is earliest (7 days) when H1-c is homozygous, as is the case for 'Cristobalina' and in all the C×C population. Previously, a large amount of homozygosity was observed in 'Cristobalina' and therefore also in the self-pollinated population (Calle et al., 2018). More specifically, a large homozygous region at the bottom of LG1 (26.23 to 47.81 Mbp), overlapping with this BT QTL was observed (Calle et al., 2018). A smaller difference (approx. 2 to 3 days) observed between the two remaining haplotypes (H1-a, -b) is in agreement with the finding that this QTL was detected at lower PVE in other works (Dirlewanger et al., 2012;Castède et al., 2014), where the allele H1-c was probably not present.
The second major QTL was identified on LG2 (qP-BT2.1 m ). This QTL also overlaps with a CR and BT QTL previously described in sweet cherry (Castède et al., 2014). The PVE and additive effect of this QTL in previous work (3.6%-6.5% PVE; 0.8-2.8 day; Castède et al., 2014) was also lower than observed in our study 5.5 JDs). This QTL has also been identified in apricot and in the interspecific cross of peach and P. davidiana, but explained a lower PVE than herein (Quilot et al., 2004;Olukolu et al., 2009;Dirlewanger et al., 2012). SOC1, a MADS-box gene, has been identified as a strong candidate gene for CR and BT underlying this QTL in sweet cherry and apricot (Trainin et al., 2013;Castède et al., 2015). However, the physical position of this gene (Castède et al., 2015) is not within the interval of the QTL detected in this work. Among other candidate genes identified in this region in sweet cherry (Castède et al., 2015), only the candidate gene, FAR-RED IMPAIRED RESPONSE 1 (FAR1) is within the interval of this QTL in our work. FAR1 has been described as a negative regulator of seasonal growth and flowering time in Arabidopsis and the loss of function of this gene resulted in plants with early flowering (Ritter et al., 2018). Therefore, this gene seems a good candidate gene for BT regulation in the genus and further work to characterize this gene in this plant material is ongoing. A larger number of haplotypes (10) were detected for this QTL, maybe due to the haplotypes being constructed across a larger genomic region, and only H2-f from 'Cristobalina' was shown to associate with earlier bloom (7 days). As observed for the QTL detected on LG1, 'Cristobalina' alleles for the underlying genes are likely responsible for the higher effect of this QTL in this plant material.
Segregation distortion was observed for some populations in both major QTLs on LGs 1 and 2. Segregation distortion in these genomic regions was previously detected in these populations (Calle et al., 2018) and in other Prunus species (Fan et al., 2010;Bielenberg et al., 2015). This distortion may be associated with segregation of lethal recessives alleles. However, since a relationship between seed and bud dormancy control has been reported (Leida et al., 2012;Abbot et al., 2015;Wang et al., 2016), it is possible that differences in seed dormancy may have affected seed germination and survival resulting in segregation distortion.

Other Minor QTLs
A minor effect QTL identified by BLUP values in this work was located on LG4. This QTL (qP-BT4.1 m ) has also been previously detected in cherries (Dirlewanger et al., 2012;Castède et al., 2014;Cai et al., 2018) and other Prunus species (Fan et al., 2010;Sánchez-Pérez et al., 2012;Zhebentyayeva et al., 2014;Bielenberg et al., 2015;Kitamura et al., 2018). This QTL has been reported as the major QTL controlling CR and BT (17.5% to 47.2% PVE) in sweet (Dirlewanger et al., 2012;Castède et al., 2014) and sour cherry , almond (Sanchez-Pérez et al., 2012) and Japanese apricot (Kitamura et al., 2018). However in our work, this QTL explained a smaller part of the variation (6.0%) ( Table  1) and was not detected all years. Several works indicated that the LG4 QTL had a larger effect on BT of high chill cultivars (Castède et al., 2014;Kitamura et al., 2018), while in low chill cultivars, as in this work, the variation in BT is more dependent in the QTL on LG1 (Fan et al., 2010;Sánchez-Pérez et al., 2012;Zhebentyayeva et al., 2014;Salazar et al., 2016). LG1 candidate genes (DAM1-6) are related to CR, and therefore these genes may have a larger contribution to BT of low chilling cultivars. In contrast, BT for high CR cultivars may be less dependent on CR, and the underlying gene(s) for the LG4 QTL has yet to be determined (Zhebentyayeva et al., 2014;Castède et al., 2015).
Three additional minor QTLs were detected in single year analyses, qP-BT1.2, qP-BT5.1, and qP-BT7.1. The QTLs on LGs 2 and 7, qP-BT1.2 and qP-BT7.1, were previously detected in sweet cherry also with small effects (Dirlewanger et al., 2012;Castède et al., 2014), but in this work it is shown that in certain environmental conditions, like those in 2016, qP-BT7.1, may have a large effect. The QTL on LG5, qP-BT5.1, has not been previously reported in sweet or sour cherry, but it has been described in the same genomic region in peach (Bielenberg et al., 2015). 'Cristobalina' was the only cultivar in this work which is heterozygous for this region (data not shown), and thus the identification of this QTL was probably due to the presence of this cultivar, and is probably associated with a rare allele found in 'Cristobalina.'

Breeding and Genome Breeding Value
The predicted genotypes for the QTL identified were used to calculate breeding value. This estimation for the parental and ancestor cultivars studied offers powerful information for breeding with these cultivars. 'Cristobalina' can be used for breeding for low CR cultivars as this work shows it is the only evaluated cultivar that exhibited early flowering due to the presence of early bloom and low chill requirements alleles in the two major QTLs affecting these traits. A similar situation was observed in peach (Hernández Mora et al., 2017), where the lowest breeding values correlated with early flowering were identified in peach landraces. This highlights the benefits of introducing exotic germplasm in breeding programs to widen the range of trait variation. Specifically for breeding for low CR cultivars with 'Cristobalina' selecting for H1-c and H2-f from QTLs qP-BT1.1 m and qP-BT2.1 m , respectively, is predicted to result in earlier blooming offspring. However, recovery of both haplotypes (H1-c/H2-f) together, may require a large number of progeny, as segregation distortion against the earlier haplotype H2-f was observed. At the same time, embryo rescue and in vitro embryo culture may be required to obtain low chilling descendants from crosses with 'Cristobalina' as the maternal parent.
If breeding for late blooming, allele H1-b rather than H1-a, should be selected for the QTL on LG1. As this QTL interval has been much narrowed in this work, and a good representation of sweet cherry breeding founders and parental cultivars is included herein, this information will also be useful for sweet cherry breeding of other plant material that do not include 'Cristobalina'. For the QTL on LG2, no conclusive evidence of late blooming haplotypes that were sufficiently predictive to be used in breeding recommendations were observed for the haplotypes detected in the parental and ancestor cultivars. For the QTL on LG4 that had a minor effect in this work, but high effect in other plant material with higher CR, selecting offspring from cultivars such as 'Rainier' and 'Brooks,' which are heterozygous for early and late bloom alleles in this QTL, would be useful for introducing an earlier allele.

CONCLUSIONS
Multi-year analysis of six pedigree-linked populations from different genetic backgrounds that descended from the low chilling 'Cristobalina,' resulted in the identification of robust BT QTLs for this highly heritable trait. Two major QTLs located on LGs 1 and 2 were identified that explained 47.6% of total phenotype variation. This work, representing the first genetic analysis of F 2 populations in sweet cherry and possible with the self-compatible 'Cristobalina,' was instrumental to characterizing the haplotype effects of these QTLs. Since the B×C F2 population resulted from self-pollination, it was possible to compare the effect of homozygous alleles in 'Cristobalina' not segregating in the other F 1 populations. BT is an essential component of cultivar adaptation to low-chill growing conditions and this trait is currently of high interest to breeders to extend sweet cherry growing to warmer areas. The identification and characterization of the haplotypes of these QTLs will enable marker-assisted breeding for this trait. The discovery of the major QTL on LG1 is consistent with the DAM gene(s) as the CR determinant in Prunus, and further suggests that 'Cristobalina' is homozygous for a unique early mutant of one or more of the DAM genes. Further work is ongoing to characterize these genes in this plant material.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the Genome database for Rosaceae, https://www.rosaceae.org/accession https://www.rosaceae.org/accession number tfGDR1040.

AUTHOR CONTRIBUTIONS
AC carried out phenotyping, data analysis and interpretation, and manuscript writing. LC advised on QTL analysis and interpretation. LC, AI, and AW participated in experimental and data analysis supervision and design and in manuscript revision. AW was also a major contributor in manuscript writing. All authors read and approved the final manuscript.