Inheritance and QTL analysis of chilling and heat requirements for flowering in an interspecific almond x peach (Texas x Earlygold) F2 population

Blooming in temperate fruit species is triggered by chilling and heat requirements (CR and HR), with a wide range of requirements within the same species. CR for flower bud dormancy release has become a limiting factor for geographical adaptation of fruit trees in warmer regions. The present study investigated the genetic basis of CR and HR to break dormancy and flowering time (FT) in an almond x peach F2 progeny. FT, HR and CR were evaluated over two consecutive years (2015/2016 and 2016/2017). Seven out of the eight identified quantitative trait loci (QTLs) were found in both periods of analysis. They affected eight traits, and included a consistent QTL for breaking dormancy, CR and HR. Two of them, affecting FT and HR for FT (GDHF), colocalized in G1, and the remaining QTLs, affecting chilling and heat requirements, both influenced by dormancy breaking (DB), were located in G6. These results indicate that factors not related to DB affect flowering time in this population. Implications of the results in peach breeding are discussed.


Introduction
Flowering, an essential and complex developmental process in plants (Castède et al. 2015), is regulated by a number of external signals and internal elements Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10681-020-02588-9) contains supplementary material, which is available to authorized users. (Hanke et al. 2007). Its correct completion is fundamental for the commercial production of seeds and fruits (Zhang and Taylor 2011).
In fruit tree orchards, there must be synchronization between flowering phenology and climatic conditions (Castède et al. 2015). As flowering is crucial for sexual reproduction, buds of perennial species in temperate regions become dormant (cease growth) during the winter months to survive. Endodormancy requires a certain amount of chilling for the transition to ecodormancy, whereas ecodormancy, requires a certain amount of heat to start the flowering process (Castède et al. 2015). Consequently, dormancy and flowering are linked, and breeders must select cultivars whose CR and flowering time match local climatic conditions (Bielenberg et al. 2015). However, it has been demonstrated that global warming can advance or delay flowering and/or fruiting of temperate fruit trees (Heide 1993;Ramirez and Kallarackal 2015;Rivero et al. 2016;Woznicki, et al. 2019), having an unknown and undesired effect on the productivity of fruit crop species. Globally, the temperature has increased by approximately 0.6°C over the past 100 years (Walther 2002). The forecast is for this trend to continue, with studies already having observed the decrease in winter chill and the resulting changes in phenological events (Menzel et al. 2005(Menzel et al. , 2006. The Prunus genus, within the Rosaceae family, is characterized by species that grow in areas with wellmarked seasons and are adapted to survive cold winters and dry summers (Dirlewanger et al. 2012). Various models have been proposed to measure the accumulation of CR in deciduous fruit-growing areas (Alburquerque et al. 2008). The effect and genetic basis of HR on flowering is not yet well understood and studies on this topic are scarce. Previous studies in Prunus have suggested that CR has a stronger effect on flowering time than HR (Couvillon and Erez 1985;Campoy et al. 2012;Alburquerque et al. 2008;Okie and Blackburn 2011;Sánchez-Pérez et al. 2012). In this genus, as in most woody perennials, the physiology and biochemistry of the flowering process is poorly understood (Dirlewanger et al. 2012;Woznicki et al. 2019). Recent reports suggest that intrinsic and environmental signaling interact and dynamically affect the extent of bud dormancy (Castède et al. 2015;Woznicki et al. 2019). Peach (Prunus persica (L.) Batsch) is an economically important species that provides an excellent system for the genetic analysis of CR and FT due to the ample variation for both traits among peach cultivars (CR between 50 and 1050 h) (Zhebentyayeva et al. 2014). Both CR and FT are inherited as quantitative traits (Fan et al. 2010;Hauagge and Cummins 1991), but their molecular regulation is not yet fully understood. Some quantitative trait loci (QTLs) associated with FT have been found in Prunus species (reviewed in Salazar et al. 2014): almond (Sánchez-Pérez et al. 2012Silva et al. 2005); peach (Dirlewanger et al. 2012;Bielenberg et al. 2015); almond x peach (Donoso et al. 2015); apricot (Campoy et al. 2013;Dirlewanger et al. 2012;Kitamura et al. 2018); sweet cherry (Dirlewanger et al. 2012;Calle et al. 2020) and sour cherry (Wang et al. 2000). The separate effect of CR and HR on FT has only been reported in a few studies in peach (Fan et al. 2010), apricot (Olukolu et al. 2009) and cherry (Castède et al. 2014). Since these traits can only be evaluated 2-3 years after seed germination in most Prunus, the identification of genetic markers linked with CR, HR and FT would be a valuable tool to select genotypes at the seedling stage and make the breeding process much more efficient (Bielenberg et al. 2015).
In this work, we studied an F 2 almond x peach progeny (TxE, almond 'Texas' x peach 'Earlygold') developed at IRTA and used as the Prunus reference map (Dirlewanger et al. 2004), for which a highdensity linkage map was available (Donoso et al. 2015). The main goals were to study the inheritance of CR, HR and FT, and to test the existing chill unit (CU) models in Gimenells, Lleida (Spain) (latitude 0°23 0 E / longitude 41°39 0 N), with a temperate semi-arid climate, to identify which of them fits best with the climatic conditions in one of the major areas of peach and almond production in the world.

Plant material
The Prunus reference interspecific almond x peach F 2 progeny (T x E), obtained by selfing a hybrid individual ('MB 1.37 0 ) from a cross between almond 'Texas' and peach 'Earlygold' (Donoso et al. 2015), and its parents was used for this study (Tables S1 and  S2). From the original progeny of 111 hybrids, we phenotyped the 72 trees that were still alive. Trees of T x E are at the IRTA Experimental Station of Lleida in Gimenells (Spain) grafted on 'Garnem' (Felipe 2009) rootstocks. Standard agricultural practices were applied.

Phenotyping
The parents, hybrid and TxE offspring were evaluated over two seasons (2015/2016 and 2016/2017, that we refer to as 2016 and 2017, respectively) following a forcing protocol widely used in temperate fruit trees (Campoy et al. 2011). Traits phenotyped were chilling requirement (CR), flowering time (FT), and heat requirement (HR). Three one-year-old fruiting branches for each individual were randomly collected once a week from November 1st until chilling requirements were reached. At least 30 flower buds were collected from the three branches for each sampling date. The bases of the branches were placed in water in a growth chamber at 25°C, under white fluorescent tubes with a 16 h:8 h, light:dark photoperiod to force floral bud break (Ruiz et al. 2007;Sánchez-Pérez et al. 2012). After 7 d, the phenological stage of the flower buds was observed. The date of dormancy breaking (DB) was established when 50% of flower buds were at phenological growth stage 53 according to the international Biologische Bundesanstalt, Bundessortenamt et CHemische Industrie (BBCH) scale (Meier et al. 1994;Alburquerque et al. 2008). Three chill models were then used to calculate chilling accumulation from October 1st until dormancy release, corresponding to the CR. The flowering time (FT) was scored as the number of Julian days when 50% of flowers were open. Also, the length of the period between DB and FT was calculated as the number of days between both events (DJD, increment of Julian Days). For each genotype, the whole tree was observed by the same person every 1 or 2 days during the flowering period. A scheme of the traits evaluated can be found in Fig. 1.

Weather data
Hourly temperatures from October 1st to flowering time for both years were obtained from the Gimenells weather station of the Generalitat de Catalunya, in the same area as the studied population (https://ruralcat. gencat.cat).

Chilling and forcing models
For this study, we used the Chill Hours (CH), Utah (CU) and Dynamic (CP) Models. The Chill Hours Model (CH) (Weinberger 1950) is the oldest and simplest model, which considers all hours with temperatures between 0 and 7C as effective for chill accumulation. The Utah Model (Richardson et al. 1974), which measures chill in Chill Units (CU), contains a weighted function attributing chilling efficiencies to different temperature ranges, including negative contributions by high temperatures, and it is particularly used in cooler areas of temperate zones (Dennis 2003). The Dynamic Model (Erez and Couvillon 1987) was developed for warmer areas. It considers that dormancy cessation occurs in two steps, the first being reversible and the second irreversible, and CR are calculated as chill portions (CP). It adopts a process-based concept of chill accumulation: an intermediate chill product is first formed through bud exposure to low temperatures, and once a critical amount of this intermediate has accumulated, it is transformed into a Chill Portion (CP). The CP is then retained until the end of the chilling period (Erez and Fishman 1998). As with the Utah model, temperatures have different effects on dormancy, but the temperature ranges differ in the two models (Alburquerque et al. 2008;Byrne 2003).
To describe heat accumulation during the later stages of tree dormancy, we used the model proposed by Richardson et al. (1974), which calculates Growing Degree Hours (GDH) between dormancy release and flowering date. According to this model, heat builds up when hourly temperatures are between 4.5 and 36°C (at different rates depending on the maximum temperature), with maximum accumulation at an optimal temperature (25°C). Additionally, we also calculated GDH between 1st October and flowering date (GDHF).

Chilling and heat requirements
Chilling and heat accumulation were calculated for the two consecutive dormancy seasons (2016 and 2017) with the hourly temperatures measured in the field. Chill was computed according to the Chill Hours (CH), Utah (CU) and Dynamic (CP) models and heat according to the GDH Model. Correlations between chilling and heating requirements and blooming were determined using Partial Least Squares (PLS) regression (Luedeling and Gassner 2012). Since heat cannot have an effect after bloom, the flowering time was considered as the end of the forcing period. CR and HR were estimated as the sum of all daily chill and heat accumulated during the chilling and forcing periods. Heat accumulation was calculated as the number of GDH from DB to FT, and the length of this period of heat accumulation was calculated (DJD) and used as a complement to GDH. Annual heat accumulation up to the date of flowering was also recorded (GDHF).

QTL analysis
For QTL analysis, we used the TxE genetic map described by Donoso et al. (2015), which was constructed using 1948 molecular markers (SNPs and SSRs), covering a total genetic distance of 472.1 cM. The interval mapping method with the MapQTL 6.0 software package (Van Ooijen et al. 2009) was used for QTL analysis of the phenotyped traits. QTLs were considered consistent when the LOD C 3.0 in both seasons, or with a LOD C 3.0 one year and LOD C 2.0 the other year. QTLs were considered as major QTLs when they explained more than 20% of phenotypic variation in both years of study (Tanksley 1993). QTL positions were drawn using the MapChart 2.1 software (Voorrips 2002).
Gene action was estimated following the guidelines of Tanksley (1993) with the ratio, d/a, between the

Temperature and chilling accumulation
Maximum and minimum daily temperatures during the consecutive years studied in Gimenells (Lleida) are shown in Supplementary Fig. S1. Higher maximum and minimum temperatures were registered during the winter in 2016, meaning a warmer winter compared to the following year. However, warmer maximum temperatures were registered at the end of the winter in 2017, reaching 25°C at the beginning of March.
Chilling accumulation was very similar over the winter in both seasons (Fig. 3) using the CU model and the CP model, but the accumulated CH was higher in 2017. Spearman correlations between models were very high for both seasons (Fig. 4).

Chill requirements to break dormancy
The range of CR for dormancy breaking (DB) for the parental lines ['MB 1.37 0 (H), 'Texas' (T) and 'Earlygold' (E)] was between 42 and 64 CP. Of these, 'Earlygold' had the lowest CR (42 and 45 CP for consecutive years), whereas the requirements for H (60 and 59 CP) and 'Texas' (54 and 64 CP) were similar. There was transgressive segregation for CR, with values from 33 to 71 CP (579-1434 CU; 484-1327 CH) averaged for two years (Table 1). CR data for the F 2 population showed normal distributions in the two years (Fig. 4), with CR skewed in both yearsto high CR. Data for each individual is given in supplementary material (Table S1 and S2).
In most genotypes, CR were similar for both years evaluated, showing high correlations between them (Table 1). The date of DB was, in general, earlier in 2016 than in 2017, with the earliest on 16th Dec and 14th Dec, whereas the latest was on 4th Feb and 2nd Feb (2016 and 2017, respectively).

Heat requirements for flowering
The range of HR for flowering among the parents was 3,575 to 6,170 GDH, whereas in the segregating population it was 2,681 to 7,105 GDH (averaged for the two years) (Fig. 4). 'MB 1.37 0 , had the lowest HR of the parents, and was similar to 'Texas', whereas 'Earlygold' had the highest. Within the segregating population, most genotypes had similar HR for both consecutive years, showing a high correlation (r 2 = 0.78) between years ( Table 1). Date of flowering (FT) was, in general, earlier in 2016 than in 2017, with the exception of only three genotypes. The parental lines ('MB1.37 0 , 'Texas' and 'Earlygold') showed very similar FT in both years. Correlation between years for FT was lower than for other traits (r 2 = 0.39) ( Table 1). The earliest FT was on 21st Feb and 3rd March, and the latest 15th and 12th March (2016 and 2017, respectively). The number of days between DB to FT in the segregating population (DJD), which is the period for heat accumulation, ranged from 27 to 77 days (averaged for the two years), and there was a high correlation between both seasons (r 2 = 0.86). HR showed a bimodal distribution in 2017, whereas it was a single peak in 2016 (Fig. 4). There were single peaks for FT in both years, although it was slightly skewed to late bloom in 2017.

Correlations between traits
Highly significant correlations (R C 0.96, p B 0.01) were found between the three models used to estimate chilling accumulation for both seasons (Supplementary Table S1). Similarly, correlations between DB and these three models were very high (R C 0.99, pB 0.01). Correlations between seasons were very high for most of the traits (r 2 = 0.78-0.89), with lower values for FT and GDHF (r 2 = 0.39 and r 2 = 0.47, respectively) (Table 1). Flowering time (FT) showed a high correlation with GDHF. Correlation between FT and CR and HR traits was low for both years of study. Heat requirements (GDH) were highly and negatively correlated with CR and DB (Table 1 and Fig. 5), which means that the lower the CR, the higher the HR. As expected, HR was highly correlated with the number of days between DB and FT (DJD), and therefore the number of days for heat accumulation was highly and negatively correlated with CR and DB.

QTL analysis
Data from the eight traits under study (DB, CH, CP, CU, FT, GDH, GDHF, DJD) were used for QTL analysis. One consistent QTL per trait was identified (Table 2). QTLs were located in G1 for FT and GDHF, and in G6 for DB, the CR models (CH, CU and CP), GDH and DJD (Fig. S2).
For dormancy breaking (DB), the LOD for the QTL from G6 was 3.4 in 2016 and 3.4 in 2017 and explained 18.4 and 18.1% of the phenotypic variance, respectively. The homozygote for the peach allele increased the date of DB by 15 days compared to the almond homozygote in 2016 and by almost two days in 2017. In 2017, while both homozygous classes had similar values, there was an increase of 13 days for the heterozygous individuals. Very similar results were obtained for DJD although the LOD for the QTL in G6 was less than 2.0 (1.8) in 2017 and therefore was not considered.
For CR data similar results were obtained with the three models, although the most significant QTLs were obtained using the dynamic model, while the least was for the Utah model. For the dynamic model, a QTL was identified at the proximal end of G6 with a LOD score of 3.4 in 2016 and 3.2 in 2017, explaining, respectively, 18.4% and 17.4% of the phenotypic variance. The confidence interval for G6 QTLs spans the genomic region Pp06: 0-4.001.078 where 702 genes have been annotated. In 2016, the individuals with the homozygous peach allele needed 10 CP more to reach DB than those with the homozygous almond allele. In 2017, this difference was only four CP, but 10 with heterozygous individuals.
In the same region of G6 we also identified a QTL for GDH and DJD (only in 2016) explaining between 16.8% and 19.9% of the phenotypic variance respectively. For FT we identified a QTL at the beginning of G1 in 2017 (LOD 3.9; R 2 = 21%), but with a LOD score of 2.2 in 2016 (R 2 = 11.5%). The confidence interval spans the genomic region Pp01:0-10.521.046 bp, where 1584 have been annotated. The peach allele in homozygosis increased the FT by five days in 2016 and three in 2017, compared to the homozygous almond allele. In the same region, we also detected a QTL for GDHF in both years, with LOD values of 3.9 and 3.1 respectively.

Phenotypic data of chilling and heat requirements
It is usually assumed that almond flowers before peach. It is interesting to note that in our case, 'Texas', a late flowering almond cultivar and 'Earlygold', an early flowering peach, flower at the same time even though they have different CR and HR. The CR of the progeny are in the range of those reported for peach and almond from other Mediterranean areas (Benmoussa et al. 2017;Campoy et al. 2012;Ruiz et al. 2007). However, it must be noted that TxE is an interspecific population, and therefore results are not fully comparable to single species populations studies. The range observed in the TxE population exemplifies the difficulties for growing certain peach and almond cultivars in warm regions where annual chill accumulation is decreasing due to global warming. A large variation in the HR within the studied progeny was also found.
In Lleida (Spain) the CR calculated using the three models were very well correlated and with DB, and therefore could be used for the calculation of CR in this climatic area, as for peach (Fan et al. 2010) and apricot (Campoy et al. 2012) in colder climates. This correlation could be due to the lack of long periods of warm and fluctuating temperatures, so that chilling accumulation based on different models all steadily increased in a similar way through the two seasons. Substantial differences among the many models used have been observed in moderate mild climates (Erez et al. 1990(Erez et al. , 2000 since some of them, such as the Utah model, were developed in a cold area and are not appropriate for warmer areas. The variability of chill accumulation between years was lower with the Dynamic model than when the calculations were done with the Utah and the hoursbelow 7°C models. This may be explained by the homogenizing effect of the Dynamic model, which takes into account the synergistic effect between moderate and low temperatures for breaking dormancy (Fishman et al. 1987a, b). Other authors have previously reported similar results in apricot (Campoy et al. 2012;Ruiz et al. 2007), suggesting that the Dynamic model is optimal for the climatic conditions of Lleida and other areas of the Ebro Valley in Northern Spain.

Correlations among endodormancy and ecodormancy traits
The high negative correlations found between DJD and GDH vs. CR (CH, CP and CU) indicate that genotypes with lower CR required a longer period of heat accumulation to bloom. This has been also found in peach (Li et al. 2016) and apricot (Campoy et al. 2012). Li et al. (2016) reported a decrease of 16 days per 200 accumulated CHs, up to a threshold of approx. 950 CHs. Overall, these results suggest that in cultivars with low CRs, GDH accumulation just after CR fulfillment is less effectivethan in cultivars with higher CRs, and therefore they need more time to accomplish their HR.
A similar result was found for the correlation between DB and CR (CH, CU, and CP) against HR (GDH). Similar high correlations have been found in peach (Fan et al. 2010;Li et al. 2016;Pawasut et al. 2004;Scorza and Okie 1990) and apricot (Ruiz et al. 2007). We found a high level of variability among the progeny regarding HR, which disagrees with Linsley-Noakes and Allan (1994) who reported no differences in HRs between three nectarine cultivars with different CR. These results suggest the existence of different heat requirements among genotypes and a major genetic contribution in the control of this trait, or that this trait is not being measured accurately because the physiological base is not yet well understood. The first hypothesis is in line with the model for Douglas fir (Harrington et al. 2010), which proposes a variable threshold for the efficiency of chill and heat temperatures. However, other authors have reported contradictory results (Couvillon and Erez 1985;Guerriero Kotowski et al. 1980), which may be due to the different climate of the cited studies. There is no consensus in the literature about whether there is a clear relationship between CR and HR.
No significant correlation was found between HR (GDH) and FT, in agreement with other authors in peach (Fan et al. 2010) and apricot (Campoy et al. 2012). This result indicates that GDH is not as important as CR for determining flowering time. However, a high correlation was found between FT and GDHF which could indicate the importance of warm temperature before dormancy breaking on the flower bud formation and development, as reported previously for plum (Woznicki et al. 2019). The fact that 'Earlygold' is a low CR cultivar might also explain this result. Also Li et al. (2016) found that both the days to full bloom date and HR were negatively correlated with CH, which may indicate that less accumulated CHs could lengthen the days to full bloom date and increase the heat requirement. Together, the results indicate that CR is a major factor determining flowering time, although not the only one. Indeed, it is unclear in the literature whether heat accumulation for floral or vegetative bud break starts before or after the release of endodormancy. Recent reports have shown a positive correlation between August-September temperature and the amount and time of flowering in the following spring on plum (Døving 2009;Woznicki et al. 2019) and sweet cherry (Døving 2011). We also found a low correlation between FT and CR. This is contrary to what has been shown by other authors in peach and almond populations (Castède et al. 2014;Sánchez-Pérez et al. 2012;Fan et al. 2010). However, there is no previous literature on CR for interspecific populations. We believe that this distortion to the expected results is due to the existence of two different species in the progeny. Indeed, in the Fig. 3 can be observed that the almond 'Texas' and the peach 'Earlygold' bloom at the same time even though E has significant lower CR than T. We have identified one genomic region controlling endodormancy and ecodormancy traits (CR and GDH) located in G6 and another one controlling FT in G1, that have not been previously identified in other Prunus populations. The interspecific nature of the TxE population might explain some of the differences observed with previous data obtained in single species mapping populations. Data for GDHF, a trait highly correlated to FT, also detected a QTL in G1 colocating with that of FT, indicating either that GDHF is not related to HR but a different way of measuring FT, or that warm temperatures during the winter dormancy period, and not only after the fulfillment of the CR, are important for determining FT. The latter hypothesis would support the high negative correlation, observed in this and other studies, between CR and HR (Fan et al. 2010;Pawasut et al. 2004;Ruiz et al. 2007;Scorza and Okie 1990), since the longer the period for chill accumulation, the higher amount of GDHF likely to be accumulated. Consistent QTLs across years for CR, HR and FT have been described in various Prunus progenies. In peach, Fan et al. (2010) identified QTLs for CR, HR and FT at the latter end of G1, where the evergrowing gene (Evg) maps (Bielenberg et al. 2008), and for CR and FT in G4 and G7, where QTLs for FT have been found in various Prunus crops (Dirlewanger et al. 2012).In similar positions of G1 and G4, QTLs for CR and/or FT have been detected by Bielenberg et al. (2015) in peach, by Sánchez- Pérez et al. (2012) in almond and by Quilot et al. (2004) in an advanced backcross between P. persica cultivars and P. Davidiana. For sweet cherry, a consistent QTL in G4 for CR and FT has also been detected by Castède et al. (2014). In all cases, the QTLs for CR coincided in map position with those of FT and produced effects of similar magnitude and gene action, suggesting CR as a major cause for the FT phenotype. Where HR was studied (Fan et al. 2010;Castède et al. 2014;Sánchez-Pérez et al. 2012), the QTLs were detected at the same positions as those of CR and FT, or were not consistent over the two years and had effects generally opposite to those of the QTLs detected for CR. This suggests that HR is a minor or irrelevant factor in the determination of FT, that its measurement as GDH is inefficient, or both. Our results support these observations as we found that HR and CR detected the same QTL in both years studied. On the other hand, we did not find common QTLs for CR and FT, indicating that FT was mainly determined by factors other than those we measured. A possible explanation for this is that the TxE offspring, from the cross between two cultivars from different species with low CR, had a lower level of variation for CR than other mapping populations studied. In the cold-winter conditions of Lleida, chilling requirements could have been rapidly met, resulting in a narrower distribution of variability that the parameters used measured with low efficiency. In agreement with our results, a QTL for FT was previously found in G1 in the TxE progeny in 2012 and 2013 (Donoso et al. 2016). Here, we also identified a peak with LOD of 2.9 in 2016 for FT in G6 but it was not considered as an stable QTL as it had a LOD \ 2 in 2017 (results not shown). This does not discount that the CR may be involved in FT variability, with apparently minor effects, although the population size used could have been insufficient to detect them with a significant threshold.

Conclusions
All the models for the estimation of CR (Utah, Dynamic and Hours-below 78) worked well for the area of study, characterized by short but cold winters, with warm falls and springs. However, the Dynamic model seems the most appropriate as it reduced the year-to-year variation observed in the population. The results indicate that, although CR appears to have a more important role than HR in determining flowering time, neither factor had a major effect on this trait under the conditions of this research. For HR, the warm temperatures during endodormacy (not only after endodormancy release) may have also influenced flowering time. In summary, our data supports FT as a quantitatively inherited character with a strong genotype x environment component that is affected by both chilling and heat requirements. The observed variation in the CRs within the population studied highlights the importance and feasibility of breeding for low CRs in a new scenario of low chill accumulation due to global warming.