Mapping and characterization QTLs for phenological traits in seven pedigree-connected peach families

Environmental adaptation and expanding harvest seasons are primary goals of most peach [Prunus persica (L.) Batsch] breeding programs. Breeding perennial crops is a challenging task due to their long breeding cycles and large tree size. Pedigree-based analysis using pedigreed families followed by haplotype construction creates a platform for QTL and marker identification, validation, and the use of marker-assisted selection in breeding programs. Phenotypic data of seven F1 low to medium chill full-sib families were collected over 2 years at two locations and genotyped using the 9 K SNP Illumina array. Three QTLs were discovered for bloom date (BD) and mapped on linkage group 1 (LG1) (172–182 cM), LG4 (48–54 cM), and LG7 (62–70 cM), explaining 17–54%, 11–55%, and 11–18% of the phenotypic variance, respectively. The QTL for ripening date (RD) and fruit development period (FDP) on LG4 was co-localized at the central part of LG4 (40–46 cM) and explained between 40 and 75% of the phenotypic variance. Haplotype analyses revealed SNP haplotypes and predictive SNP marker(s) associated with desired QTL alleles and the presence of multiple functional alleles with different effects for a single locus for RD and FDP. A multiple pedigree-linked families approach validated major QTLs for the three key phenological traits which were reported in previous studies across diverse materials, geographical distributions, and QTL mapping methods. Haplotype characterization of these genomic regions differentiates this study from the previous QTL studies. Our results will provide the peach breeder with the haplotypes for three BD QTLs and one RD/FDP QTL to create predictive DNA-based molecular marker tests to select parents and/or seedlings that have desired QTL alleles and cull unwanted genotypes in early seedling stages.


Background
Peaches and nectarines [Prunus persica (L.) Batsch] are deciduous fruit trees belonging to the Rosaceae family. These are native to China and grown throughout the world in a wide range of environments. The gross production value of peaches and nectarines in 2016 was $825 million in the United States and $17,107 million globally [1].
Breeding of woody perennial crops is not an easy task since their long juvenility periods and large plant size makes maintaining large populations in the field expensive [2]. The use of marker-assisted breeding (MAB) provides a tool to do an early selection of seedlings, identify superior parents, improve the selection of elite alleles for essential traits, and stack desirable alleles [3,4]. This strategy is pertinent for perennial fruit tree to reduce breeding operational costs [3].
Ripening date in peach trees is a crucial element for extending the production season and developing cultivars that ripen throughout the harvest season. Also, the ripening process is involved in the regulation of several metabolic pathways such as blush, sugar/acid balance, and the flesh softening in peach fruits [31]. Narrow sense heritability (h 2 ) for ripening date ranges from high to very high (0.79-0.94) [15,32,33]. The major QTL for controlling RD was mapped on LG 4 at~44 cM in the Prunus T × E reference map, and a putative candidate gene was located at~10.5 Mbp on the peach genome sequence v.1 [30,31,34,35]. This QTL explained~50 to 98% of the phenotypic variability. The RosBREED project has verified this locus is significant in the U.S. breeding programs [18]. Likewise, a QTL for RD on chromosome 4 was detected in apricot, sweet cherry [31], and almond [36].
Currently, DNA-based tests for a few breedingrelevant traits have been developed and used in the peach marker-assisted selection application, including maturity date (G4mat) [39], quality traits, and fruit bacterial spot resistance. Thus, work is needed to develop DNA tests for BD and FDP traits and to validate SNP-based DNA test (G4mat) for ripening date to enable their use in the TX and other breeding programs [3,[40][41][42].
The objectives of this study are to identify new and/or validate the major QTLs previously reported for bloom date, ripening date, and fruit development period through pedigree-based analysis approach (PBA) using Texas peach/nectarine germplasm. Also, to estimate QTL genotypes for important breeding parents and to identify predictive SNP marker(s) associated with desired QTL alleles. Results from this research will facilitate the design of DNA tests linked to these QTL(s) or genes for routine use for marker-assisted breeding.

Phenotypic data analysis
The mean BD value ranged from 42.3 ± 3.9 (CA11) to 50.2 ± 9.45 (TX13), and a maximum range of 51, with the number of observations between 82 in CA11 and 143 in overall mean (Additional file 1: Table S1). In our study, BD distribution varied across environments and overall mean (Additional file 2: Fig. S1). The CA environments were skewed towards the lower values, whereas the TX exhibited multimodal profiles with two or more peaks in both environments. This was expected as some of the higher chill genotypes had delayed bloom in the lower chill Texas site compared to California (5 40 vs.~1090 chilling hours) [43]. Normal distribution was seen in the overall mean of BD. RD exhibited an average between 129.2 ± 16.7 (TX12) and 157.4 ± 17.7 (CA11), with greater (87.0) and lower (59.5) RD ranges in the overall mean and TX13 data sets, respectively. CA and the overall mean data sets were slightly skewed towards the higher values, while the TX data sets were skewed towards the lower values. On average, fruit ripened approximately 17 days later at Fowler, CA than at College Station, TX. FDP mean values ranged from 81.2 ± 16.9 (TX12) and 115.3 ± 16. 9 (CA11) with FDP range from~67 (CA12) to 91 (TX13) days. The minimum number of observations (59) was recorded for CA11 compared to 138 observations for the CA12 and overall mean data sets. Similar to RD, FDP for CA data sets were slightly skewed towards higher values compared to the TX environments, which skewed towards lower values while the overall mean showed normal distribution. Fruit had development periods that were 23 days longer, on average, at Fowler, CA, than at College Station, TX. This was an effect of cooler temperatures during early fruit development in March and April for CA11 and CA12 (~15 and 9°C) relative to TX12 and TX13 (21 and 18°C).
Among these traits, a strong correlation was found between RD and FDP (r = 0.91) (Additional file 1: Table  S2), and a moderately weak correlation was observed between FDP and BD (− 0.45). The negative correlations between BD and FDP suggest that the earlier blooming genotypes experience a delay in the rate of fruit development due to cooler temperatures. A weak correlation was found between RD and BD traits (− 0.14).

Genotype by environment interactions
The genotype × environment interaction (G × E) is the differential sensitivity of genotypes to different environments. If such interaction exists, the selection would be complicated and result in genetic gains reduction in a breeding program. Understanding the G × E interactions is key to increasing the efficiency of marker-assisted selection for complex traits [44].
In this study, RD and FDP showed very high broadsense heritability (H 2 = 0.95 and 0.96, respectively), strong correlations among environments (r = 0.91), and minimal G × E variance ( σ 2 gÂe =σ 2 g ratio = 0.20) (Additional file 1: Table S3 and S4) whereas BD trait, showed highbroad-sense heritability (H 2 = 0.88), strong correlation among environments (r = 0.83) and a moderate genotype by environment interaction ( σ 2 gÂe =σ 2 g ¼ 0.70). All traits had comparable PC2 values and ranged from 5.5 to 6.8 (Additional file 1: Table S5), implying that the environments equally discriminate the populations for these traits. Finally, the minimal G × E effect of RD and FDP is supported by the relatively similar length of the environmental vectors in the GGE biplots, especially within the same location, indicating a high correlation among them and equal discriminatory ability of the four environments (Additional file 2: Fig. S2). Also, the distance between the environmental vectors was closer between CA11 and CA12, and between TX12 and TX13 for RD and FDP, respectively, illustrate that genotypes responded similarly in these two environments. This is confirmed by the highest positive correlations between CA11 and CA12 (r = 0.87, RD and 0.84, FDP) and between TX12 and TX13 (r = 0.79, RD and 0.89, FDP) for RD and FDP) (see Additional file 1: Table S4).
For BD, the sharper angle and less distance were observed between CA12 and TX12, TX12 and TX13, and CA12 with TX13 (Additional file 2: Fig. S2), indicating a stronger correlation between these environments (r = 0.73, 0.75, and 0.65) (Additional file 1: Table S4). The best discrimination of BD among genotypes was observed in the CA11 environment indicated by the longer vectors for these environments (Additional file 2: Fig.  S2). Also, the environment CA11 was far from the other three environments and showed less correlation coefficient. However, the low number of observations of this environment (82) may have affected the correlation and G × E results.
Three QTLs were mapped for BD on three linkage groups (LG1, 4, and 7) across the four environments (site × year combinations) and their overall mean. The QTL on LG1 was at the distal end and showed strong to decisive evidence in all data sets (Table 1 and Additional file 2: Fig. S3). The QTL on LG4 was mapped in three environments (except CA11) and the overall analysis, showing positive and decisive evidence. At the same time, the QTL on LG7 was seen in only two environments and the overall analysis with decisive evidence. FlexQTL software found one to two candidate QTLs for RD and FDP depending on the environment; however, only the QTL on the middle part of G4 passed our inclusion criteria. (Table 1 and Additional file 2: Fig. S4 and S5).
For BD, the proportion of phenotypic variation explained (PVE) ranged from 17 to 54%, 11 to 55%, 11 to 18% for LG1, LG4, and LG7, respectively ( Table 2). The highest posterior QTL intensity (0.96) showed in LG1 for BD-mean, and the lowest intensity (0.21) was found in LG4 for BD-TX12. The highest additive effect (~10 days) was in LG4 for BD-TX13, and the lowest (~2 days) showed in LG1, 4, and 7 for BD-CA12. The QTL on LG1 was co-localized across all data sets with an interval between 172 and 182 cM (peaks, 174, 176, and 178 cM), and the physical position of this chromosomal region was 43,058,300 -45,586,061 bp on the peach genome sequence v2.0, ( Table 2, Fig. 1a, and Additional file 1: Table S6). Likewise, peaks of QTL on LG4 of three data sets, except CA12, clustered at mode 50 cM, with an interval between 48 and 54 cM and physical chromosomal position between 11,956,738 -13,633,831 bp. Regarding LG7, the peaks co-localized at either 64 or 66 cM with an interval from 62 to 70 cM and physical chromosomal position between 15,513,277 -17,226,623 bp on the peach genome sequence v2.0 ( Table 2 and Fig. 1b, Additional file 1: Table S6).
The proportion of phenotypic variation explained by RD QTL on LG4 ranged between 46 and 75% ( Table 2). The highest posterior QTL intensity (1.80) and the highest additive effect (~19 days) were found in CA12. In most environments, the observed high intensity (greater than one) implies that FlexQTL assigned two QTLs within the same QTL interval with an average distance between them of 1.0 cM across all sampled models. This distance is very short to be genetically meaningful for population sizes. This QTL had mode at either 44 or 45 cM, overlapping intervals from 40 to 46 cM across all data sets, and the physical chromosomal position between 10,396,616 to 11,298,736 on the peach genome sequence v2.0 ( Table 2, Fig. 1c, and Additional file 1: Table S6). The proportion of phenotypic variation explained by FDP QTL on LG4 ranged between 40 and 71% ( Table 2). The highest posterior QTL intensity (1.60) was for CA12 and the lowest (0.79) for TX12. The highest additive effect (~20 days) was found in TX13. Likewise, this QTL had a mode at either 44 or 45 cM, overlapping intervals from 40 to 46 cM across all data Markov chain Monte Carlo (MCMC) run length, phenotypic mean (μ), phenotypic variance (σ 2 P ), residual variance(σ 2 e ), additive variance(σ 2 A ), narrow-sense heritability (h 2 ), the linkage groups (LG) that QTLs were mapped on 2ln(BF). Bayes Factor, a measure quantifies the support from the data for the number of QTLs in the model (QTL evidence), after pair-wise model comparison (1/0, 2/1, and 3/2) such as 'one-QTL model' vs. 'zero-QTL sets, except TX12, and has a physical chromosomal position between~10,396,616 to 11,298,736 bp of the peach genome sequence v2.0 (Table 2 and Additional file 1: Table S6). Like RD, the high intensity that is noticed in most data sets indicates two tightly linked QTLs within the QTL interval, and the gap between them averaged to 1.4 cM across all sampled models. So, the distance is also too short to be genetically dissected in these studied population sizes.
QTL associated haplotypes, number of QTL-alleles, their effect, predictive markers, and sources On LG1, 11 SNPs in the predicted qBDG1 region (172.23-182.34 cM) (Additional file 1: Table S7), chosen for haplotyping, revealed eight SNP haplotypes across the seven parents in which H8 was a common haplotype ( Table 3). The estimation of the diplotype effect identified families of two parents (Y434-40 and 'Victor') were segregating for this QTL. The results also discovered multiple Q-alleles of various effects associated with H1 to H7, and only one q-allele was linked to low phenotypic values associated with H8.
The examination of the haplotype /diplotype effects (Fig. 2a) revealed that the effect of H7 and H1could not differentiated when comparing H5H7<>H5H1 and H8H1<>H8H7. Likewise, the effects of H5 and H8 could not be differentiated when comparing H5H1 to H8H1 and H5H7 to H8H7. Also, H7 had a larger effect than H8 and H3 in the comparison H8H7<>H8H8 and H8H7<>H8H3, respectively. The effect size of H1 was greater than H2 and H3 when comparing H8H1 to H8H2 and H8H3. In general, H8 had a smaller effect than H1, H2, H3, H6, and H7, when comparing H8H8 to H8H1, H8H2, H8H3, H8H6, and H8H7. Hence, H1 and H7 had similar and the largest effects, and both coined as Q1, then followed by H3, H6, H2, and H8, which were represented as Q2, Q3, Q4, and q, respectively. However, the under-representation of QTL genotypes hindered the estimation of H4 and H5 effects. All of these haplotypes could be differentiated from H8 by various pairs of adjacent SNP markers by contrasting either ABor BA-alleles for 1) snp_1_46757382 and ss_ 135737 to BB of H8, or 2) ss_128625 and ss_128603 to AA of H8, and 3) ss_129512 and ss_128603 to also AA of H8 (Table 3 and Additional file 1: Table S7). Breeding parents 'Galaxy', Y426-371, Y435-246, Y434-40, and TX2B136 were considered as founders in this study and the sources of these SNPs were unknown because their ancestors were not available for genotyping. On the other hand, the Q-allele (H5) of 'Victor' was inherited from F_Goldprince, and the q-allele (H8) of both 'Victor' and TXW1490-1 was inherited from Fla3-2 through 'TropicBeauty'.
On LG4, there were 13 SNP markers in the BD QTL region (47.83 to 54.54 cM) (Additional file 1: Table S7) selected for haplotyping. That revealed five SNP haplotypes in the seven parents. H1 and H3 were the most common haplotypes ( Table 3). Families of four parents (Y435-246, Y426-371, 'Galaxy', and 'Victor') were heterozygous for this QTL. H2 and H3 were associated with the Q-allele while H1, H4, and H5 with the q-allele.
The examination of the haplotype/diplotype effects in Fig. 2b revealed that H3    effects/magnitudes of some haplotypes on BD, e.g., H5H3 (qQ) had a larger effect than H3H2 (QQ). That could be explained by several reasons such as the presence of interaction with other loci, H5 (q) having a smaller effect on decreasing BD among the other haplotypes (H1 and H4) associated with decreasing BD, or H2 (Q) having less magnitude on increasing BD. The low number of diplotype observations or high variance within a diplotype class might also have caused these issues. More than one predictive SNP marker associated with H2 and H3 (Q-allele) were identified (Table 3). A-allele at ss_415301 (50.09 cM) along with three more SNP markers distinguished H3, whereas the A-allele at ss_ 414387 (48.43 cM) and the other two SNP markers were unique for H2. In contrast, H1, H4, and H5 (q-allele) were distinguished by two adjacent BB-alleles at ss_ 414387 and ss_415301. The H3 Qallele was found in TX2B136, 'Flordaprince', F_TXW1490_1, Y426-371, and F-Goldprince while the H2 Qallele came from Y435-246 and 'Galaxy'. The qalleles were found in Y435-246, Y426-371, 'Galaxy', Y434-40, Fla3-2, and 'TropicBeauty'.
15 SNP markers in the predictive QTL region for both RD and FDP traits (42.33 to 45.19 cM) (Additional file 1: table S7), on the middle part of LG4, were picked for haplotype analyses. FlexQTL implies this genomic region had more than one QTL within the same interval Results discovered four SNP haplotypes associated with RD and FDP across the seven parents of which H1 QTL alleles for each parent cultivar are presented with ♀ and ♂ for maternal and paternal parent sources, respectively. Parents that are heterozygous for the QTL are in bold. Allele(s) for predictive SNP marker(s) associated with Q or q-alleles for increasing or decreasing a given trait, respectively, are shown inunderscored bold. Q/q of different effect magnitude are indicated by subscript numbers. The identity of the SNP markers and their physical and genetic location is given in Additional file 1: Table S7 was common ( Table 3). Families of five parents (Y426-371, Y434-40, 'Galaxy', 'Victor', and TXW1490-1) were segregating in this region. The diplotype analysis revealed the presence of four statistically distinct phenotypic classes (Fig. 3 a and b). H3 had a larger effect than H1 and H4 when comparing H4H3<>H4H1 and H1H3 <> H4H1, respectively. Likewise, H2 showed a smaller effect than H1 on both RD and FDP when comparing H1H1<>H1H2 and from H4H1<>H4H2 just in FDP not RD as their effects could not be differentiated (Fig. 3 a and b). Thus, the effect size of haplotypes can be ordered as H3 > H4 > H1 > H2 that is differentiated by Q1, Q2, q1, and q2, respectively. The major finding in this study was the presence of multiple QTL alleles of different effects for a single locus. That may explain why the Bayes Factor values and high intensities of most data sets of this study suggested the presence of two QTLs.

Discussion
In this study, the bloom date was moderate to highly heritable (0.44-0.82) as has been previously reported [15,[24][25][26][27] in a range of germplasm, indicating that expression of bloom date is not heavily influenced by environmental effects which were supported by G × E results. Narrow sense heritability was moderate to high for RD (0.59 to 0.83) as was found in previous studies [15,18,26,[46][47][48]. FDP also has an important additive genetic component as indicated by a high to very high (0.65 to 0.82) estimated narrow-sense heritability reported in this and previous studies [15,[24][25][26][27].
The QTL at the middle region of LG4 for BD mapped between ss_413934 and ss_419614, in the interval between 12 and 13.6 Mb, and PVE ranged from 11 to 55%. This QTL overlaps with the BD QTL on LG4 (qFD4.2) at nearest markers ss_417840 and ss_440116 (13.1 to 16.0 Mb) reported by Hernández Mora, et al. [15].
Lastly, the QTL at the distal end of LG7 was flanked by ss778568 and snp_7_17628094, spanned from 15.5 to 17.2 Mb and explained~11 to 18% of BD phenotypic variation. This finding agreed with Romeu, et al. [30] who found a QTL for BD on LG7 at the nearest marker ss_779224 (15.7 Mb), which was close to our QTL peaks (ss_780816 (16.3 Mb) and ss_779362 (15.7 Mb)). Moreover, this region overlapped with the QTL (15.4 to 19.4 Mb; PVE~60%) reported by Fan, et al. [29].
The only one of the three QTLs was detected in CA11 is probably due to that this environment had a low number of phenotypic data (82 records). The G × E for BD in the studied populations may result from the response of the high-chill seedlings to the lack of chill hours that delayed the blooming period.
In summary, this study provides more evidence that three mapped QTLs for BD on LG1, 4, and 7 are major loci for controlling BD and were supported by other studies using low-and medium-chill germplasm and biparental family mapping. It was also supported by the polygenic nature of BD inheritance. Additional QTLs for BD were also reported on LG2 [15,49], LG3 [17,30], LG6 [15,30,50], and LG8 [15,17]. Thus, further studies using more diverse germplasm will be important to continue to characterize additional QTLs and candidate genes to identify the genetic pathway regulating the BD in peach.
The examination of haplotype/diplotype effects uncovered the high prevalence of a few haplotypes, e.g., H8 (qallele), H3 (Q-allele), and H7 (q-allele) on LG1, 4, and 7, respectively, reflecting the relatively narrow genetic base of peach germplasm. Also, the results revealed the presence of multiple Q-alleles of different effects for the QTL on LG1 (Q1, Q2, Q3, and Q4) along with only one q-allele. In general, the small family sizes and consequently the low/lack representation of various compound diplotypes (e.g., 6 to 9 observations in some diplotypes of LG1) hindered the ability to make conclusions on the haplotype effects (H4 and H5) or the interplay among the three mapped QTLs for BD.
One QTL associated with RD and FDP was mapped at the middle part of LG4 (10.4-11.3 Mb) with PVE 46-75% and 40-71%, respectively. This specific genomic region was reported as associated with RD trait previously by Nuñez-Lillo, et al. [35] (~10.9 Mb), Romeu, et al. [30] (~10.7 Mb), Frett [18]  . This held true using early-, mid-, and late-maturing populations. The co-localization between QTLs for RD and FDP was supported by the strong correlation (r = 0.87) (data not shown) between these traits in this study as well as previous work [6,15].
Also, all data sets, except TX12, showed decisive evidence (BF ≥ 10) with high intensity for the presence of a second QTL on LG4. This could be explained by that TX12 had higher temperatures during the critical fruit development months (March and April) [51] compared to other sites. The higher temperatures accelerated RD and shortened FDP in this environment, which minimized the phenotypic variation as mentioned earlier (Additional file 1: Table S1).
Furthermore, the haplotype analysis of this chromosomal region revealed multiple predictive loci (ss_ 410398, ss_410794, and ss_412662) for decreasing and increasing for either RD or FDP. Likewise, examining the relative effects of haplotypes and estimated QTL genotypes revealed a series of QTL alleles of different effect at this locus that we coined Q1, Q2, q1, and q2. The use of multi-parent populations for finding multiple functional alleles of different effect was also reported for two acidity QTLs/genes in apple by Verma, et al. [52] and for the blush QTL in peach using the current germplasm by Rawandoozi, et al. [16]. In our germplasm, the RD QTL on LG4 co-localized with a QTL for soluble solids concentration (SSC) and blush reported by Rawandoozi, et al. [16]. These co-localizations had also been reported by other studies [15,34,53]. A pleiotropic effect of the RD has been reported on several quality traits [15,34,35,39]. Co-factor analysis could be useful in future studies to account for one trait when analyzing another, e.g., accounting for RD for analyzing SSC or blush traits.
Overall, additional QTL mapping through pedigreebased analysis across a wider range of breeding germplasm is needed to identify and characterize additional QTLs to understand the whole genetic pathway controlling RD and FDP traits. Moreover, larger family sizes would ensure better representation of QTL genotype classes for estimating QTL effects and allow improved downstream analysis in case of multiple QTL alleles of different effects at a single locus and/or gene by gene interaction.
At the genomic region of the detected QTLs for these traits, candidate genes have been reported. For BD, the QTL interval (43,058,300 -45,586,061 bp) of LG 1, the most promising candidate genes for the major QTL affecting blooming time and chilling requirement in LG1 were the Dormancy-associated MADS-box (DAM) genes within the evergrowing (evg) locus in peach, apricot, and almond [29,54,55].
Prupe.1G531600 (DAM5) and Prupe.1G531700 (DAM6) genes were identified as potential candidate genes of lateral bud endodormancy release in peach [29,56,57]. Prupe.1G531500 gene is described as MADS-box protein short vegetative phase (SVP) and it plays a role in controlling meristem development during the vegetative phase and flower development as well as in floral meristem determination [58]. Prupe.1G549600 and Prupe.1G548000 genes are described as agamous-like MADS-box proteins AGL11 and AGL12, respectively. AGL11 is a vital gene to control ovule identity and associated placental tissues in Arabidopsis [59]. While a MADS-box gene AGL12 regulates root development and flowering transition in Arabidopsis [60]. Prupe.1G554100 (AGL80) is also a member of the MADSbox family of genes. In Arabidopsis, AGL80 was found to be involved in female gametophyte development [61].
Likewise, many candidate genes have been reported within the interval (11,956,633,831 bp) of LG4. Prupe.4G208000 is described as a Forkhead-associated (FHA) domain-containing protein (DDL) that plays an important role in plant growth and development [62]. Prupe.4G197000 gene was proposed to link to auxin synthesis and response which is known to be involved in fruit set and ripening [63]. Prupe.4G202200, Fertilization Independent Endosperm (FIE) polycomb group protein, in Arabidopsis thaliana FIE regulates endosperm and embryo development and suppresses flowering during embryo and seedling development [64]. Prupe.4G207300 (uclacyanin) is associated with pollen grain development in rice [65]. Prupe.4G205500 (early nodulin-like protein 1) gene is reported to be engaged in determining the reproductive potential in Arabidopsis [66]. In the QTL region (15,513,226,623 bp) of LG7, Prupe.7G130900, CURLY LEAF (CLF) gene, is associated with the repression of FLOWERING LOCUS T (FT) gene and other flowering-time genes during the vegetative growth of the plant [67]. Prupe.7G153400 gene is described as a ATP-dependent DNA helicase (DDM1), the importance of this gene was previously reported for DNA methylation in genes and transposable elements [68]. Prupe.7G133100 (Zeaxanthin epoxidase) gene has been identified to play an important role in resistance to stresses, seed development, and dormancy in Arabidopsis [69].
Within the RD/FDP locus on LG4 (10,582,092 to 11,298, 736), a list of candidate genes has been previously reported in this region. NAC072 (Prupe.4G816800) is the candidate gene for controlling the ripening date in peach [39]. Also, there are three other genes proposed to be involved in the determination of RD/FDP in peach. Prupe.4G79900 gene is needed for normal embryo development in Arabidopsis and maize [70,71]. Prupe.4G179800 gene is described as Early nodulin-like protein 1 and PtNIP1in Arabidopsis and loblolly pine, respectively [72]. It is expressed in immature zygotic and somatic embryos of developing seeds. Prupe.4G179200 gene with functional annotation Purine permease 10 in Arabidopsis and OsPUP7 in rice [73], and showed a flowering delay in rice. Finally, Prupe.4G185800 [74] and Prupe.4G187100 [75] genes that were reported to be associated with the regulation of the anthocyanin biosynthetic pathway in peach. Hence, these results confirming the pleiotropic effect of the RD on several quality traits, including blush that was previously reported [15,16,34,35,39].

Conclusions
The pedigree-based analysis was successfully used as a statistical method for discovering and validating QTLs. Four QTLs associated with three important phenological traits were validated using low-medium-chill peach/nectarine germplasm. Two minor QTLs were also identified. This approach increases the genetic background explored, improves statistical power, and allows the simultaneous detection and validation of QTLs.
QTLs for BD on LG1, 4, and 7 were verified, and the SNP haplotypes associated with increasing or decreasing BD were identified. A single QTL with multiple QTL alleles of different effects was detected on the central part of LG4 for both RD and FDP. Our findings would help breeding programs make crossing decisions to pick the combination of parents that have SNP haplotypes associated with lowering BD to produce progeny with better adaptation to subtropical environments like Texas or increasing BD to ensure better adaptation to temperate environments, whereas the results of RD and FDP will facilitate better targeting for specific ripening periods. Ultimately, the SNP haplotypes associated with these QTLs could be converted into easy-to-use high throughput markers (e.g., Simple Sequence Repeat (SSR), Kompetitive Allele-Specific PCR (KASP), and Sequence Characterized Amplified Region (SCAR) markers) to routinely use in MAB. In general, this approach would save time and resources, particularly for fruit breeders since perennial wood species have long juvenility periods, and large populations are expensive to maintain in the field.

Plant materials
Briefly, we included in this study 143 seedlings from seven related F 1 families derived from seven parents descending from 12 founders. The parents are all cultivated germplasm that has been developed by the Stone On the other hand, 'Galaxy', Y435-246, Y424-40, and Y426-371 were developed by the USDA Stone Fruit Breeding program from 'Armking' and germplasm from Rutgers University (New Brunswick, NJ), the University of Florida (Gainesville, FL), and Georgia USDA Stone Fruit Breeding program (Byron, GA). Seedlings and parental genotypes were grown in College Station, TX (30°37′41.60″N, 96°22′27.38″W), and Fowler, CA (36°38′21.37″N, 119°42′20.51″W). Full details on plant materials and plot establishment and design can be found in Rawandoozi, et al. [16].

Phenotypic evaluations
Phenotypic data were taken at both locations across 2 years (2011-2012 in CA, and 2012-2013 in TX) on individual trees for three phenological traits, bloom date (BD), ripening date (RD), and fruit development period (FDP). The date of first (10% blossoms open) and full bloom (60 to 80% of the blossoms open) were visually assessed in the field and recorded for each tree. Ripening date was determined when 20% of fruits are pickable by visually inspecting the presence of a few soft fruits in the field for maturity two times per week. Both full bloom and ripening dates were converted to Julian days (0-365). FDA is difference in days between BD and RD.

Heritability and G × E
Variance components for the studied traits were estimated using a linear mixed model with the residual maximum likelihood (REML). Results from REML were used to estimate the broad-sense heritability across the environments, as explained by Rawandoozi, et al. [16]. The R package GGEBiplots version 0.1.1 was used to estimate the variations due to genotypes and G × E. Pearson's correlation coefficients were also estimated among phenotypic traits within and across the environments using R software version 4.0.3.

Genotyping and linkage map
Plant samples were genotyped using the IPSC 9 K SNP Array for Peach [11], and SNP data were curated following the workflow described by Vanderzande, et al. [76]. After filtration, a total of 1487 informative SNPs were distributed over eight chromosomes using a conversion factor in which every 1 Mb corresponded to 4 cM [76].

QTL mapping
FlexQTL software (version 0.1.0.42) with an additive genetic model conducted by Markov Chain Monte Carlo (MCMC) simulation was used for QTL mapping. The analysis was run at least three times on each data set. Different prior and maximum QTL numbers were used in each run to reach effective chain size (ECS) ≥ 100 for the mean, variance of the error, number of QTLs, and QTL variance, as recommended to draw reliable and accurate conclusions [13,77]. MCMC length ranged from 100,000 to 3600,000 iterations to store one thousand samples with a thinning between 100 and 3600. Convergence was evaluated visually via trace and intensity plots [13]. Twice the natural logarithm of Bayes Factors [2ln(BF)] obtained from FlexQTL software used as evidence for presence and number of QTLs [78]. The 2ln(BF) value greater than 2, 5, or 10 indicate positive, strong, and decisive evidence, respectively. In this study, loci were considered if QTL had 2lnBF ≥ 5 or that 2 ≤ 2lnBF < 5 for at least two data sets, the QTLs with overlapping intervals of at least 2 cM on the same linkage group, and explained at least 10% of the phenotypic variation.
The additive (σ 2 AðtrtÞ Þ, phenotypic ðσ 2 P Þ, and residual ð σ 2 e Þ variances were obtained from FlexQTL output to estimate the narrow-sense heritability (h 2 ), and the proportion of phenotypic variance explained (PVE) as follows: : σ 2 AðqtlÞ isthevarianceofQTL The QTL nomenclature in this study described by Rawandoozi, et al. [16] is a modification of that of Fan et al. [29].

Haplotypes analysis
SNPs within the significant QTL interval were considered for haplotype analysis using the FlexQTL software and PediHaplotyper package of R [19]. Haplotype effects were determined from combinations of diplotypes by comparing the effects of the H1|H2 and H1|H3 diplotypes. The nonparametric multiple comparison Steele-Dwass test (P < 0.05) was used to assess the significance of differences using JMP Pro Version 13.2 (SAS Institute Inc., Cary, NC, 2016) as described by Rawandoozi, et al. [16].