This study was conducted on two multi-parental diploid rose populations (TX2WOB and TX2WSE). The TX2WOB population was developed from crosses involving breeding lines (J14-3, J3-6, J4-6, and M4-4) derived from Rosa wichuraiana ‘Basye’s Thornless’ crossed with ‘Old Blush’, ‘Red Fairy’, ‘Sweet Chariot’, ‘Vineyard Song’, and ‘Little Chief’ ( Supplementary Figure 1 ) (Dong et al., 2017; Yan et al., 2019). The TX2WSE population was developed from crosses among R. wichuraiana ‘Basye’s Thornless’-derived breeding lines (TAMU7-20, TAMU7-30, J14-3, and M4-4) and ‘Papa Hemeray’, ‘Srdce Europy’, ‘Ole’, R. setigera-ARE, and R. palustris f. plena EB-ARE ( Supplementary Figure 2 ) (Young et al., 2022). In addition, these studied diploid populations are segregating for the miniature phenotype (small leaves and flowers).

Five F₁ TX2WOB rose populations (387 plants) were planted in 2010 at the Horticulture Farm of Texas A&M University in College Station, TX, USA (30.63, −96.37) in one replication and phenotyped in 2015. A subset from TX2WOB (N = 300) was planted in 2018 at the Horticulture Teaching Research and Extension Center (HortTREC) in Somerville, TX (30.524591, −96.422479) with two replications and evaluated in 2021 ( Supplementary Table 1A ). Six F₁ TX2WSE populations (N = 353) were planted in 2018 at HortTREC with two replications and phenotyped in 2021( Supplementary Table 1B ). More details on populations and field conditions are described by Rawandoozi et al. (2022) for TX2WOB and Young et al. (2022) for TX2WSE.

Temperature data were obtained from both field locations in 2015 and 2021 (Weather Underground, 2018). The data showed that the temperature varied among seasons and years. The average maximum temperature in the summer season (June, July, and August) in 2015 was hotter than in 2021 (34.1°C vs. 33.3°C). Similarly, 2015 had cooler spring and fall seasons than 2021 ( Supplementary Table 2 ; Supplementary Figure 3 ).

In 2015, the flower size parameters, including flower diameter (Diam), flower dry weight (DWT), and the number of petals (NP), were measured during three seasons [spring (April), summer (August), and fall (November)]. Phenotypic data were collected from at least three fully open flowers randomly chosen in each plant. Flower diameter (cm) was measured in the field, whereas the other two traits were measured in the lab. Flower dry weight (mg) was taken after the whole flower without the pedicel was dried for at least 3 days at 80°C. The number of petals included full-size petals and petaloids (irregularly shaped petals). In 2021, the phenotypic data for flower size were taken from four flowers per plant in the summer (June and July). Diam and NP were measured as described above. Additionally, in this year (2021), both flower fresh weight (mg) (FWT) and the number of petaloids (PD) were measured. The least-square mean (lsmean) of the phenotypic data was used in the statistical analyses and was estimated using “emmeans” v. 1.7.5 package of R (v. 4.1.2; R Foundation for Statistical Computing, Vienna, Austria).

FlexQTL has the functionality to include co-factors and was utilized in the analysis of NP. Since the major QTL for NP was located in the same genomic region harboring the DOUBLE FLOWER locus on LG3, the latter was considered as a covariate in the analysis. Thus, individuals with flowers with fewer than eight petals were considered to have single flowers, whereas those with eight or more petals were considered to have “double” flowers and were given the values of 1 or 2, respectively.

A Shapiro–Wilk test was performed to test the normality of raw and transformed data. Heritability was only estimated for traits measured in 2015 as the data was taken in three seasons using mixed models with a restricted maximum likelihood (REML) estimation method in JMP Pro version 13.2 (SAS Institute Inc., Cary, NC. USA), treating all effects as random (Littell et al., 1996). The following model was used:

where μ is the mean; σ² _FP and σ² _MP are the female (FP) and male (MP) parent variances, respectively; σ² _{Progeny(FP, MP)} is the progeny variance; σ² _Env is the environmental variance (seasons); σ² _{FP × Env}, σ² _{MP × Env}, and σ² _{Progeny × Env} are variances due to the interaction of female and male parents and progenies with the season; and σ² _error is the error variance.

The sum of the parental variances (σ² _FP and σ² _MP) was treated as an additive variance ( σ A 2 ), progeny variance [σ² _{Progeny(FP,MP)}] was as considered nonadditive variance ( σ d 2 ), and the sum of the parental and progeny variances was regarded as the genotypic variance ( σ g 2 ). The interaction of genotype [σ² _FP, σ² _MP, and σ² _{Progeny(FP,MP)}] by environment (season) was treated as the genetic–environmental variance ( σ g × e 2 ). The residual variance, confounded with progeny × environmental variance, was regarded as the error variance (σ² _error).

Broad sense heritability for each set of populations across environments was calculated as:

where E indicates the number of environments (seasons) (Holland et al., 2003; Liang et al., 2017b; Wu et al., 2019; Rawandoozi et al., 2021).

The genotype by environment variance to the genetic variance ratio was estimated as:

A genotype and genotype-by-environment (GGE) biplot analysis was employed to understand the variation due to genotype using the R package “GGEbiplots” v. 0.1.3. Pearson correlation coefficient among phenotypic traits and seasons/years was calculated.

Two consensus maps were developed for each diploid rose population. The TX2WOB consensus map (415 individuals) was constructed from five rose populations ( Supplementary Table 3 ). The TX2WSE was created from three rose populations (314 individuals) ( Supplementary Table 4 ).

Genomic DNA was extracted from new rose leaves using Doyle’s CTAB protocol (Doyle and Doyle, 1991). In this study, genotyping by sequencing (GBS) was accomplished using the digital genotyping procedure according to the method described by Morishige et al. (2013). Single-end sequencing was performed using an Illumina HiSeq 2500 platform. The CLC Genomics Workbench v9.0 (Qiagen, Boston, MA) was used to align the reads to the Rosa chinensis v1.0 genome (Hibrand Saint-Oyant et al., 2018). Markers were named based on their position in the rose genome and grouped into bins based on their proximity to a given restriction enzyme cut site (NgoMIV) in the reference genome (Yan et al., 2018).

For TX2WOB, before the five individual linkage maps were developed, low-quality SNP markers were eliminated. Tassel version 5 was used to remove markers if they were mapped to chromosome 0, non-biallelic, and markers that had >10% missing data. After that, a Microsoft Excel-based tool and custom R scripts were employed to remove markers with inheritance errors. Then, individual maps were developed for each population using the R package “polymapR” v. 1.1.1, which was set to remove duplicated and distorted markers (p ≥ 0.001). The datasets for each population were simplified by choosing one marker of each marker class per restriction-enzyme bin, which is defined as the region around a NgoMIV cut site, giving preference to markers that are common between populations, having little missing data, and fitting expected segregation ratios. Next, the consensus map was developed using the R package “LPmerge” v. 1.7. The R packages “LinkageMapView” v. 2.1.2 and MapChart software v. 2.32 were utilized to visualize the consensus map. Additional marker curation was performed in FlexQTL software v. 0.1.0.42 to identify and fix/remove problematic markers (singletons, double recombinations, and inheritance errors).

Regarding TX2WSE, the similar procedures described above were used to generate a linkage map, with the exception of filtering markers in PLINK v. 1.9 to eliminate Mendelian-inconsistent mistakes per population. More information about the linkage map development can be found in Rawandoozi et al. (2022) and Young et al. (2022).

QTL detection via pedigree-based analysis (PBA) was implemented through the FlexQTL software (Bink et al., 2012; Bink et al., 2014). The dataset for the TX2WOB population includes phenotypic data collected from three seasons in 2015 (spring, summer, and fall), the average of three seasons (Mean 2015), and one season in 2021 (summer), along with 1,115 SNP markers. The dataset for the TX2WSE population consists of phenotypic data from one season (summer 2021) and 866 SNP markers.

First, traits were analyzed with a mixed (additive and dominance) model. Then, due to the absence of a dominant effect, QTL analysis was performed in the only additive model at least three times with different parameter settings to ensure reproducibility (Verma et al., 2019). Markov Chain Monte Carlo (MCMC) simulations length ranged from 100,000 to 1,600,000 iterations to store 1,000 samples with thinning between 100 and 1,600.

Convergence was evaluated visually by trace and intensity plots. Twice the natural logarithm of Bayes Factors [2ln(BF)] obtained from FlexQTL software was used to determine the strength of QTLs (Kass and Raftery, 1995). The 2ln(BF) value >2, 5, and 10 indicates positive, strong, and decisive evidence, respectively (Kass and Raftery, 1995; Bink et al., 2014). This study considered major loci if those QTLs were mapped for at least two data sets with 2lnBF ≥ 5, overlapping intervals, and explaining at least 10% of the phenotypic variation. In this study, the physical positions of all mapped QTLs for a given trait(s) were compared across environments and populations. In addition, within each population, if the QTL intervals clustered in the same genomic regions for a trait, these QTLs were considered the same QTL.

The additive ( σ A ( t r t ) 2 ) , phenotypic ( σ P 2 ) , and residual ( σ e 2 ) variances were obtained for each trait from FlexQTL output and used to estimate the narrow-sense heritability (h² ), and the total phenotypic variance explained (PVE) by a QTL was calculated as follows:

where: σ A ( t r t ) 2 is the variance of the trait

where: σ A ( q t l ) 2 is the variance of QTL

QTLs were named using the Genome Database for Rosaceae QTL naming conventions (Jung et al., 2019). For example, qNP.TX2WOB-LG3.1, where q stands for QTL, “NP” is the phenotypic trait name (number of petals), “TX2WOB” or “TX2WSE” is the population name that used to develop the consensus map, “LG3” is the linkage group number, and numbers “1” or “2” are used in case there was more than one QTL detected within the same LG. A “CF” suffix was added to the QTL name to distinguish mapping based on co-factor analysis.

SNPs within the region of major QTL(s) (proximity to QTL peaks) for each trait were chosen for haplotype analysis. FlexQTL and the R package “PediHaplotyper” v. 1.0 were used to construct haplotypes (Voorrips et al., 2016). Combinations of diplotypes were used to infer haplotype effects. Diplotype effect differences were evaluated using Steel–Dwass non-parametric multiple comparison test JMP Pro version 13.2. QTL allele groups (Q or q) were assigned to haplotypes based on their effects. In the case of multi-allelic series, Q and q alleles were differentiated by an index number.

Haplotypes were traced through pedigrees records. Haplotypes that were traced to a common ancestor were considered identical by descent (IBD), whereas haplotypes that could not be traced to a known common ancestor were defined as identical by state (IBS).

The narrow sense heritability (h ²) varied among traits, environments, and populations. It was low for PD (0.16–0.39); low to moderate for Diam (0.31–0.56), FWT (0.38–0.46), and DWT (0.40–0.57); and low to moderately high for NP (0.38–0.74) ( Table 1 ). The broad sense heritability was moderately high to high (0.75–0.87) for Diam, DWT, and NP ( Supplementary Table 6 ).

Diam exhibited moderately high broad sense heritability (H ² = 0.75), moderate to strong correlations among environments (r = 0.52–0.87) ( Supplementary Table 7 ), and a moderate ratio of G×E to G variance ( σ g × e 2 / σ g 2 = 1.00 ) ( Supplementary Table 6 ). These results suggest the existence of genotype sensitivity to the environment. However, the other two traits, DWT and NP, showed high to very high broad sense heritability (H ² = 0.82 and 0.87, respectively), strong correlations among environments (r = 0.70–0.93), and minimal genotype by environment interaction ( σ g × e 2 / σ g 2 = 0.65 and 0.43, respectively). These results were further supported by the greater amount of variation explained by the first principal component (PC1) (82.61% and 92.73%) for DWT and NP as compared to Diam (80.83%) ( Supplementary Figure 7 ). Furthermore, the GGE biplots for Diam showed that the summer vector was far from the other two seasons, indicating that the high temperatures in the warm season (average monthly summer temperature ~29°C) discriminate genotypes differently than the cool seasons (average monthly spring and fall temperature ~20–22°C). In addition, cool seasons for this trait discriminated genotypes similarly, supported by the short distance between them. This finding is confirmed by the correlations among these three seasons ( Supplementary Table 7 ). In contrast, the GGE biplots showed for NP and DWT that fall is more discriminating and discriminates differently than spring and summer seasons.

In general, for both NP and DWT, the σ g × e 2 / σ g 2 ratio were<1, suggesting that we should be able to make progress in selecting individuals across different seasons as the genetic variance is greater than the G×E variance. In comparison, the selection for Diam is season dependent. Consequently, the selection for Diam during both cool and warm seasons would be preferred to develop a stable flower diameter.

Among traits, the correlation results showed that flower weight, whether it was DWT (TX2WOB, 2015) or FWT (both populations in 2021), was positively correlated with Diam (r=0.27–0.47) and NP (r= 0.53–0.85) ( Supplementary Table 8 ). This finding was expected as the flower weight reflected the combination of Diam and NP. A weak and negative correlation was found between Diam and NP, reflecting a weak tendency of the larger flowers to have fewer petals in the studied rose populations.

In addition, high positive correlations were found between FWT and NP with PD in TX2WOB (r= 0.78 and 0.79) and TX2WSE (r= 0.89 and 0.69). Lastly, there were strong positive correlations between the 2 years’ data (2015 and 2021) in TX2WOB for Diam (r=0.66), fresh and dry flower weight traits (r=0.76), and NP (r=0.93) ( Supplementary Table 8 ).

Several QTLs associated with flower size traits were identified across the two diploid rose populations and six phenotypic data sets (spring, summer, fall, and the mean of 2015 for TX2WOB and summer 2021 for both TX2WOB and TX2WSE), which were evaluated over two locations in Texas ( Table 1 ).

For TX2WOB, FlexQTL software detected nine QTLs associated with Diam on all linkage groups (LGs), except LG7, across five seasons ( Table 1 ; Figure 1 and Supplementary Figure 8 ). A major QTL (qDiam.TX2WOB-LG1) was consistently detected on the distal end of LG1 across two cool seasons (spring and fall) and mean 2015 data sets between 68.2 and 69.8 cM [60.44–60.95 Mbp on the Rosa chinensis genome v1.0 (Hibrand Saint-Oyant et al., 2018)] with mode (peak) at 69.0 cM ( Table 2 ; Figure 1 and Supplementary Figure 8 ). This QTL showed decisive evidence, high posterior intensity, and the phenotypic variation explained (PVE) by this QTL ranged from 67% to 80% ( Tables 1 , 2 ). In this study, qDiam.TX2WOB-LG1 passed our inclusion criteria.

Two minor QTLs were detected on LG2 and LG6. qDiam.TX2WOB-LG2.1 was mapped at the middle part of LG2 between 49.7 and 57.9 cM (48.36–56.70 Mbp) in four data sets with strong evidence, and PVE ranged from 11% to 15% ( Tables 1 , 2 ; Figure 1 ). qDiam.TX2WOB-LG6.2 was detected between 58.4 and 61.8 cM (61.67–63.22 Mbp) in two data sets with positive evidence and PVE of 9%–12%. Most other QTLs identified in this population were environment specific.

Three QTLs were mapped in the TX2WSE population for Diam with positive/strong evidence on LG 2, 3, and 7 using one data set (summer 2021) ( Table 1 ; Figure 2 and Supplementary Figure 8 ). qDiam.TX2WSE-LG2 was between 106.52 and 114.97 cM (71.30–72.31 Mbp) and PVE of 12%. qDiam.TX2WSE-LG3 and qDiam.TX2WSE-LG7 were at the proximal ends of LG3 and LG7, respectively, with PVE of 12%–21% ( Table 2 ).

For the TX2WOB population, four QTLs associated with DWT were identified ( Table 1 ; Figure 1 and Supplementary Figure 9 ). Of these, major QTLs on LG1 and LG3 were consistently detected with either strong or decisive evidence at the same positions in at least two seasons and showed high PVE ( Tables 1 , 2 ). Hence, these two QTLs were considered for further analysis. qDWT.TX2WOB-LG1 was consistently mapped across three data sets in 2015 at the distal end of LG1 between 68.2 and 69.8 cM (60.44–60.95 Mbp) and mode at 69.0 cM with large effects (PVE, 40%–41%) ( Table 2 ; Figure 1 ).

The second major QTL on LG3, qDWT.TX2WOB-LG3, was identified at the middle part of LG3 with intervals ranging from 42.4 to 54.0 cM (30.61–39.34 Mbp) across three seasons and the mean in 2015. Peaks for this QTL co-localized at positions between 43.0 and 51.9 cM, with PVE ranging from 24% to 38% ( Table 2 ; Figure 1 ). The other two QTLs on LG2 and LG5 were mapped once and considered environment specific.

One major QTL for FWT was mapped with decisive evidence using TX2WOB and one season (summer 2021) on LG3, qFWT.TX2WOB-LG3 ( Table 1 ; Figure 1 ). The QTL interval ranged from 42.4 to 46.7 cM (30.61–33.57 Mbp) and PVE of 23% ( Table 2 ).

On the other population, TX2WSE, two QTLs were mapped for FWT on LG3 and LG5 using one data set (summer 2021) ( Tables 1 , 2 ; Figure 2 ). qFWT.TX2WSE-LG3 was a major QTL mapped with an interval ranging from 25.38 to 33.53 cM (30.15–33.83 Mbp), decisive evidence, and PVE of 46%.

For TX2WOB, one to three QTLs for NP were detected per environment on LG3, two of these passed our inclusion criteria ( Table 1 ; Figure 1 and Supplementary Figure 10 ). The first major LG3 QTL that was mapped at the middle part (qNP.TX2WOB-LG3.2) was common to all five data sets examined across 2 years, with PVE of 34%–76% ( Tables 1 , 3 ; Figure 1 ). Peaks for this QTL co-localized across all environments, having their mode at 43.0, 46.2, and 46.7 cM, and their interval between 42.4 and 46.7 cM (30.61–33.57 Mbp), except for fall 2015, which had wider intervals spanning 42.4–51.9 cM (30.61–36.23 Mbp) ( Table 3 ). The second major QTL on LG3, qNP.TX2WOB-LG3.1, was located upstream qNP.TX2WOB-LG3.2, at 25.4–31.7 cM (18.88–22.06 Mbp) and identified in the fall and the mean 2015 data sets with PVE up to 37% ( Table 3 ; Figure 1 ). The LG2 QTL was mapped at the upper part of the LG and was environment specific in this study, as it appeared only once in the dataset.

As for TX2WSE, one major QTL (qNP.TX2WSE-LG3) with decisive evidence and large effect (58% PVE) was identified at the middle part of LG3 between 25.38 and 33.53 cM (30.15–33.83 Mbp) with a peak at 31.88 cM ( Tables 1 , 3 ; Figure 2 and Supplementary Figure 10 ).

Moreover, the posterior intensity of qNP.TX2WOB-LG3.2 and qNP.TX2WSE-LG3 was generally >1 ( Table 3 ), which implies that FlexQTL assigned two QTLs within the same QTL interval. This was supported by Bayes factors values, which showed evidence of more than one QTL ( Table 1 ). However, the downstream analysis of “iqtl.out” file from FlexQTL output showed that the distance between these QTLs was very short (ranging from 1.7 to 2.9 cM) (data not shown). FlexQTL is not equipped to distinguish QTLs within this short distance.

Since the major QTLs of NP (qNP.TX2WOB-LG3.2 and qNP.TX2WSE-LG3) colocalized with the DOUBLE FLOWER locus on LG3 in all data sets in both populations, this trait was considered for co-factor analysis to assess the effect of this locus of masking the detection of additional QTLs with minor effects. The results from the co-factor analysis were similar to those of the original data set in TX2WOB. However, the major QTL on LG3 (qNP.TX2WOB-LG3.CF) had a larger interval between 39.0 and 46.7 cM (27.65–33.57 Mbp) in most data sets, 25.4 and 51.9 cM (18.88–36.23 Mbp) in fall, and 25.4 and 41.3 cM (18.88–29.38 Mbp) in summer 2015 with PVE of 12%–28% ( Tables 1 , 3 ; Supplementary Figure 11 ). In addition, the minor QTL on LG2 was mapped in the middle part of the LG. However, for TX2WSE, a new QTL was mapped upstream to qNP.TX2WSE-LG3 between 18.39 and 29.88 cM (27.80–29.08 Mbp) on LG3 (qNP.TX2WSE-LG3.CF) with PVE of 41% ( Table 3 ). The QTL mapping results using the DOUBLE FLOWER locus as a covariate helped identify an additional major QTL on LG3 for TX2WSE and another minor QTL on LG2 for TX2WOB. In addition, a new chromosomal region was discovered on LG3 for TX2WOB that was confounded with one or both major QTLs (qNP.TX2WOB-LG3.1 and qNP.TX2WOB-LG3.2).

Regarding PD, one QTL on LG3 with decisive evidence and large effect was mapped for each population using one data set (summer 2021) and was considered for the downstream analysis ( Tables 1 , 2 ; Figures 1 , 2 and Supplementary Figure 12 ). qPD.TX2WOB-LG3 was detected between 46.2 and 51.9 cM (33.14–36.23 Mbp) with PVE of 17% in TX2WOB ( Table 2 ), whereas qPD.TX2WSE-LG3 was between 29.24 and 31.88 cM (31.09–34.36 Mbp) and PVE of 42% ( Table 2 ).

For Diam, only two QTL genotypes (QQ and Qq) were predicted at the peaks of qDiam.TX2WOB-LG1 and qDiam.TX2WOB-LG2.2, where q and Q were associated with low and high Diam, respectively ( Supplementary Figures 13A, B ). Progenies were categorized into QQ and Qq groups with a mean Diam of 3.4 and 2.8 cm for qDiam.TX2WOB-LG1, and 3.5 and 3.4 cm for qDiam.TX2WOB-LG2.2, respectively. In contrast, the three genotype classes for qDiam.TX2WSE-LG2 had Diam means of 3.6, 3.3, and 3.1 cm for QQ, Qq, and qq classes, respectively ( Supplementary Figure 13C ). Generally, progenies with the favorable allele (Q) were associated with increasing Diam and were more prevalent in TX2WOB dataset compared to TX2WSE population.

Regarding the flower dry weight for TX2WOB, two QTL genotypes were predicted for both major QTLs. QTL genotypes had an average DWT of 6.4 and 3.2 mg for qDTW.TX2WOB-LG1 and 8.5 and 7.2 mg for qDTW.TX2WOB-LG3, for progenies having the QQ and Qq genotypes, respectively ( Supplementary Figures 14A, B ). The two QTL genotypes classes at qFTW.TX2WSE-LG3 have estimated mean FWT of 280.7 and 141.3 mg for Qq and qq genotype classes, respectively ( Supplementary Figure 14C ). No individuals were in the qq genotype class in TX2WOB ( Supplementary Figures 14A, B ). Therefore, more favorable alleles (Q) are associated with increasing flower weight in TX2WOB as compared to TX2WSE.

Two QTL genotype groups were identified at the LG3 QTLs for NP across the two populations. For TX2WOB, QTL genotypes had an average NP of 26.7 and 11.4 for qNP.TX2WOB-LG3.1 and 35.5 and 15.4 for qNP.TX2WOB-LG3.2 for progenies having the Qq and qq genotypes, respectively ( Supplementary Figures 15A, B ). In contrast, progenies were categorized into Qq and qq groups with a mean NP of 25.8 and 7.4 for qNP.TX2WSE-LG3 ( Supplementary Figure 15C ). In summary, unfavorable alleles (q) associated with decreasing petals were common in TX2WSE (approximately 50% of progenies). The opposite was true in TX2WOB.

As for PD, three QTL genotypes classes were estimated at qPD.TX2WOB-LG3 with a mean of PD 2.9, 2.5, and 2.0 for QQ, Qq, and qq classes, respectively ( Supplementary Figure 16A ). In contrast, individuals were classified into Qq and qq groups with an average of 3.0 and 0.6 for qPD.TX2WSE-LG3, in which half of the individuals had had q-allele associated with fewer PD ( Supplementary Figure 16B ).

Thus, lacking or unbalanced representation of QTL genotype classes in most studied traits due to the small sample size of some populations hindered our ability to make robust conclusions about QTL gene action.

This analysis considered two major QTLs that were mapped on LG3 in TX2WOB, which showed high PVE and posterior intensity. The four SNPs (25.4–31.7cM) spanning ~3.2 Mbp in the qNP.TX2WOB-LG3.1 region were chosen for haplotyping ( Figure 7A ). Three SNP haplotypes were discovered in which G1 and G2 were associated with increasing NP and assigned to the Q-allele, while G3 was linked to decreasing the trait (q-allele) ( Figure 7A ). G2G1 was the most common diplotype (216 individuals) ( Figure 7B ). The analysis of the haplotype effects showed that the effects of G1 and G2 could not be differentiated based on G2G1 and G1G1. Meanwhile, G1 showed a greater effect than G3 (G1G1 to G1G3). The pedigree information indicated that the three haplotypes were IBS. G1 was inherited from ‘Old Blush’, PP-M4-2 ‘Little Chief’, or ‘Violette’ ( Figure 7A ); G2 was from PP-J14-3, ‘Ducher’, or R. wichuraiana; and G3 came from PP-M4-4, ‘Old Blush’, or R. wichuraiana.

For qNP.TX2WOB-LG3.2, eight SNPs (42.4–46.7cM) spanning ~3 Mbp in the QTL region were selected for the analysis ( Figure 7C ). Three SNP haplotypes were detected. H1 and H2 increased NP (Q-allele), while H3 was linked to decreasing the trait (q-allele). H1 was the most common haplotype ( Figures 7C, D ). The diplotype analysis revealed the presence of three QTL-alleles with different effects on NP and ordered as Q ₁ (H1) > Q ₂ (H2) > q (H3). H1 was inherited from various sources (PP-J14-3, ‘Old Blush’, ‘Little Chief’, ‘Red Fairy’, or ‘Violette’) ( Figure 7C ), H2 came from ‘Old Blush’ or ‘Little Chief’, and H3 originated from ‘Ducher’, R. wichuraiana, WOB26, or PP-M4-2.

As for TX2WSE, eight SNPs were also chosen to conduct haplotype analysis in the qNP.TX2WSE-LG3 region (25.38–33.53 cM) spanning ~3.7 Mbp ( Figure 8A ). Seven SNP haplotypes were identified, in which J1, J2, and J3 were associated with the increasing NP and were given Q-alleles. J4, J5, J6, and J7 were related to decreasing the trait (q-allele) ( Figure 8A ). The diplotype analysis revealed that J1 and J2 (Q-alleles) had a similar effect when comparing J2J7 to J1J7 ( Figure 8B ). No difference was observed between J5 and J6 (J3J5 to J3J6) or between J7 and J6 (J1J7 to J1J6 and J5J7 to J5J6). J1 had a greater effect than J3 (J1J6 to J3J6), and the latter also had a greater effect than J4 and J5 (J3J5 to J4J5 and J3J6 to J5J6). Thus, three QTL alleles of different effects were identified at this locus. The haplotype effects order was [J1 =J2]>J3>[J5=J6=J7] corresponding to Q ₁, Q ₂, and q, respectively. However, the under-representation of some diplotype compounds hindered our ability to conclude the magnitude of the J4 (q-allele) effect compared to other haplotypes associated with q-alleles ( Figure 8B ). The highest and lowest diplotypes were observed in J1J6 (36.6) and J4J5 (6.8).

The sources of J1 were MEV, R36, or ‘Violette’, J3 was PP-J14-3, while J2 was a recombinant haplotype and originated from recombination events between the parents of ‘SC’ (‘LC and ‘Violette’) ( Figure 8A ). J4 came from R. wichuraiana, J5 from ‘Old Blush’ and PP-M4-4, J6 from HIA or ‘Sweet Chariot’, and J7 from ‘Sweet Chariot’.

In addition, the interplay between qNP.TX2WOB-LG3.1 and qNP.TX2WOB-LG3.2 was assessed ( Supplementary Figure 18 ). The analysis revealed one dose of either Q- or q-alleles at the qNP.TX2WOB-LG3.2 increased/decreased NP more than those at qNP.TX2WOB-LG3.1 ( Supplementary Figure 18 ). This indicates that qNP.TX2WOB-LG3.2 has a larger effect than qNP.TX2WOB-LG3.1. This result was supported by diplotype effects (~41 vs. ~15) ( Figures 7B, D ) and QTL genotype (~18 vs. 15) ( Supplementary Figures 15A, B ) analyses for qNP.TX2WOB-LG3.2 and qNP.TX2WOB-LG3.1, respectively. The effect of qq-genotypes at both loci could not be determined for the same reason mentioned above. Accordingly, future QTL mapping studies using broader and more diverse germplasm are crucial to enhance the representation of three QTL genotype classes/diplotype combinations.

Haplotype analysis was conducted on the LG3 PD QTL mapped in each population. In the TX2WOB population, five SNP markers between 46.2 and 51.9 cM spanning ~3.0 Mbp in the qPD.TX2WOB-LG3 were selected for haplotype analysis ( Figure 9A ). Five distinct SNP haplotypes were identified, of which K1 and K2 were associated with increasing PD (Q-allele), while K3, K4, and K5 were linked to decreasing PD (q-allele) ( Figure 9A ). The assessment of the haplotype/diplotype effects showed the effect of K1 seemed to be similar to that of K2 (K5K2 to K5K1) ( Figure 9B ). Likewise, K3 had the same effect as K4 (K3K1 to K4K1). K2 had a greater effect on increasing PD than K3 by comparing K2K1 to K3K1. In addition, K4 increased PD more than K5 (K4K1 to K5K1). Thus, the haplotypes were ordered from higher to lower effects, [K1=K2]>[K3=K4]>K5 and were assigned the Q, q ₁, and q ₂ QTL alleles, respectively.

K2K1 showed the highest PD (~3.7), K4K5 had the lowest (~2.0), and K1 was the most prevalent haplotype ( Figure 9B ). K1 was traced back to various sources (PP-J14-3, ‘Old Blush’, ‘Red Fairy’, or ‘Little Chief’) ( Figure 9A ). K2 came from ‘Red Fairy’ and ‘Violette’, K3 came from R. wichuraiana, K4 from ‘Old Blush’, while K5 was inherited from ‘Ducher’, ‘R-Wich ‘Old Blush’, or PP-M4-2.

In the TX2WSE population, six SNP markers between 29.24 and 31.88 cM spanning ~3.3 Mbp were used in the haplotype analysis of qPD.TX2WSE-LG3. Five distinct SNP haplotypes were identified ( Figure 9C ). L1 was the only haplotype associated with increasing PD (Q-allele), and L2, L3, L4, and L5 were associated with decreasing the trait (q-allele). The results also showed that L1 had a larger effect than L3 and L4 based on L1L2 to L3L2, L1L4 to L3L4, and L1L2 to L4L2 diplotypes ( Figure 9D ). L2 and L4 (q-alleles) had a similar effect when comparing L1L2 to L1L4 and L3L4 to L3L2. Likewise, L5 increased PD more than L2 (L1L5 to L1L2). Thus, there were two QTL alleles of different effects associated with decreasing PD, q ₁ (L5) > q ₂ (L2=L4). However, the effect of L3 (q-allele) was not determined due to the low representation of some diplotype combinations ( Figure 9D ).

The difference between the highest (L1L5) and lowest (L3L2) was ~3 petaloids. The pedigree map showed that L1 was inherited from PP-J14-3, MEV, R36, and ‘Violette’ ( Figure 9C ). The source of L2 was either ‘Sweet Chariot’ or HIA, while L3 and L5 came from R. wichuraiana and ‘Sweet Chariot’, respectively. Lastly, ‘Old Blush’ and PP-M4-4 were the sources for L4.

In this study, narrow sense heritability was low to moderate for Diam, consistent with the prior studies (Liang et al., 2017a; Liang et al., 2017b). In addition, the moderate narrow sense heritability estimated for DWT and FWT agreed with those previously reported for DWT (Liang et al., 2017a; Liang et al., 2017b). Regarding NP, moderate to moderately high narrow sense heritability was found in this and previous studies (Gitonga et al., 2014; Liang et al., 2017a; Liang et al., 2017b). Lastly, low to moderate narrow sense heritability was observed for PD.

Moreover, all flower size traits, Diam, DWT, and NP showed a moderately high to high broad sense heritability (0.75, 0.82, and 0.87) ( Supplementary Table 6 ) as previously reported (Gitonga et al., 2014; Liang et al., 2017a; Liang et al., 2017b). This study highlighted the important additive and non-additive genetic components for the studied traits.

The G×E interaction was varied from low ( σ g × e 2 / σ g 2 ratio, 0.43–0.65) for NP and DWT to moderate ( σ g × e 2 / σ g 2 ratio, 1.00) for Diam. These results indicate that DWT and NP have a strong genetic component with less environmental influence than Diam (Gitonga et al., 2014; Liang et al., 2017a), and hence, the selection for these two traits can be made irrespective of seasons.

The research conclusively showed that all examined traits, over the two evaluated years, experienced a significant reduction when subjected to heat stress (average summer monthly temperature, approximately 28–28.6°C). However, these traits increased under cooler conditions (average spring and fall monthly temperature, approximately 20–22°C) ( Supplementary Table 2 ; Supplementary Figure 3 ). This aligns with the findings of previous studies (Gitonga et al., 2014; Greyvenstein et al., 2014; Liang et al., 2017a; Liang et al., 2017b).

DWT was moderately positively correlated with Diam and petal number over all data sets in both populations ( Supplementary Table 8 ), as previously reported (Liang et al., 2017b). Therefore, this suggests that the increase in DWT appears to be due to an increase in Diam and/or NP.

In contrast, the correlation between the Diam and the petal number was low in this diploid (Liang et al., 2017b) and in tetraploid rose populations (Yu et al., 2021), indicating that these traits are inherited independently as reflected in the major QTLs for NP and Diam being on different LGs. Furthermore, a high positive correlation between NP and PD in both populations was anticipated as both traits result from transformations from stamens (Dubois et al., 2010).

FlexQTL software detected one major QTL on LG3 associated with PD that was consistently mapped across both populations. This QTL was between 33.14 and 36.23 Mbp, and 31.09 and 33.83 Mbp for qPD.TX2WOB-LG3 and qPD.TX2WSE-LG3, respectively, with PVE 17%–41%. To our knowledge, this is the first QTL report for PD.

Therefore, future studies using data over several environments (seasons, years, and locations) from more diverse germplasm are essential to test robust LG3 QTL across environments and identify whole genetic pathways that regulate PD in roses.

Haplotype analyses for the major QTLs and those common among populations revealed several SNP haplotypes and predictive SNP marker(s) associated with increasing/decreasing QTL alleles. In addition, multiple functional alleles with different effects were found for most traits of this study.

According to pedigree information, for Diam, the ancestors ‘Old Blush’ and R. wichuraiana from J14-3 were common sources of Q-alleles across populations, while ‘Violette’ was the common source of Q-allele over LG1 and LG2 loci of TX2WOB. In addition, ‘Violette’ and PP-J14-3 were the common sources for Q-alleles for NP, PD, FWT, and DWT over the two populations. Moreover, ‘Little Chief’, a very small diameter flowered rose, was the common source of q-allele for Diam across populations, while R. wichuraiana, a single (five petals) flowered species rose, through J14-3 was a common source of q-allele for most loci associated with NP, PD, FWT, and DWT across populations.

Alternatively, the genetic information of estimated diplotype effects can be utilized for future selection, such as A7A4 (LG1), B1B5 (LG2) for TX2WOB, and C2C1 (LG2) for TX2WSE for larger Diam; regarding flower weight, D7D4 (LG1) and E1E1 (LG3) for TX2WOB, while F2F7 and F1F7 (LG3) for TX2WSE; and NP, G2G1, G1G1, and H1H1of both major QTLs on LG3 for TX2WOB, whereas J1J6, J2J7, and J1J7 for TX2WSE. These results will aid breeders in parental selection to develop cultivars with high ornamental quality, consequently increasing the market value of roses. Eventually, major QTLs for these traits could be used to develop high-throughput DNA tests for routine use in a DNA-informed breeding program.

In this study, QTLs of flower size traits were clustered at the same genomic region on LG1 and LG3 in the rose genome ( Figure 10 ). The first QTL cluster was at the lower part of LG1 (~60.44–60.95 Mbp) for Diam (qDiam.TX2WOB-LG1) and DWT (qDWT.TX2WOB-LG1) ( Figure 10A ). The haplotype results showed that both loci are in the coupling phase with each other, supported by a positive correlation between the traits.

The genomic region located between 60.44 and 60.95 Mbp on chromosome 1 of the Rosa chinensis genome v1.0 (Hibrand Saint-Oyant et al., 2018), according to data derived from the GDR database, encompasses 87 candidate genes ( Supplementary Table 9 ). Among those, five genes appear to directly affect plant growth and/or floral development. A number of the cytochrome P450 was reported in this genomic region. A study by Liu et al. (2015) reported CYP715A1 (At5g52400), the lone representative of P450 family in Arabidopsis thaliana, underlining its fundamental role in regulating petal expansion, emission of volatiles, and flower maturation. Furthermore, actin-depolymerizing factor 7 (ADF) was found in this region, which is known to be a key regulator of F-actin organization, flowering, and cell and organ expansion in Arabidopsis. This F-actin was reported to be required for various cellular processes, including cell division and expansion (Kost et al., 1999). Similarly, a study (Li et al., 2009) on cotton (Gossypium hirsutum) reveals that ADF genes (GhADF6 and GhADF8) contribute to petal development and found that the low levels of the ADF genes were detected at the early stages of petal development, while a surge in the gene expression was observed at the later stages of petal development. F-box protein, one of the super protein families, has also been identified in this region. Zhang et al. (2019) pointed out the diverse functionalities of the F-Box genes, including roles in development processes such as plant hormonal signal transduction, floral development, circadian rhythms, responses to biotic and abiotic stresses, and others. The other candidate gene in this QTL interval is FAR1-related sequence 5, which is recognized for its vital role in plant growth and developmental processes, including shoot meristem and floral development in Arabidopsis (Li et al., 2016; Ma and Li, 2018). Lastly, two SEUSS-like 2 genes were found in this genomic region on chromosome 1, which was reported to have a role in the negative regulation of the floral homeotic gene AGAMOUS (AG) in Arabidopsis, which is responsible for specifying stamen and carpel identity in plants. Franks et al. (2002) suggested that SEUSS mutant causes the partial transformation of floral organs and a mild decrease in the number of floral organs. Furthermore, double mutation of SEUSS and AG led to narrow sepals and petals and reduced plant height.

The second cluster found in TX2WOB was at the middle region of LG3 (~30.61–39.34 Mbp) harbored QTLs for the FWT (qFWT.TX2WOB-LG3), DWT (qDWT.TX2WOB-LG3), NP (qNP.TX2WOB-LG3.2), and PD (qPD.TX2WOB-LG3) ( Figure 10B ). Similarly, in TX2WSE, a single cluster of major QTLs was identified for FWT, NP, and PD at the same region in LG3 (~30.15–34.36 Mbp) ( Figure 10C ). In addition, the haplotype results revealed that DWT/FWT and NP shared identical haplotypes associated with increasing/decreasing phenotypic values and their sources, which also partially coincided with PD. Thus, the increase in flower weight seems to be due to an increase in flower diameter and the number of petals or petaloids. This finding was consistent with positive correlations that were found among these traits.

This QTL region on LG3 overlaps with the DOUBLE FLOWER locus (33.24–33.53 Mbp). The candidate gene in this region, APETELA2/TOF (33.23–33.24 Mbp), is proposed to control the switch from simple to double flower and the number of petals within double flowers (Hibrand Saint-Oyant et al., 2018). Moreover, this specific region was previously reported to be associated with the QTL for number of petals (Debener and Mattiesch, 1999; Linde et al., 2006; Koning-Boucoiran et al., 2012; Roman et al., 2015; Bourke et al., 2018; Hibrand Saint-Oyant et al., 2018; Schulz et al., 2021; Yu et al., 2021).

Interestingly, this middle region on LG3 is considered a hot spot in the rose genome as major or candidate genes associated with key morphological traits have been identified such as the continuous flowering locus (RoKSN) (Iwata et al., 2012), and prickle density (TTG2, TESTA TRANSPARENT GLABRA2) (Hibrand Saint-Oyant et al., 2018; Zhou et al., 2020), and QTLs associated with resistance to black spot and cercospora diseases in roses (Rawandoozi et al., 2022; Rawandoozi et al., 2023) have also been reported in this genomic region.

Three additional clusters for minor QTLs were found in TX2WOB. The region between 48.36 and 64.33 Mbp on LG2 was associated with Diam, DWT, and NP. Yu et al. (2021) also mapped the LG2 QTL for Diam (qfdia-2-2), which co-localized with the minor LG2 QTL for Diam (48.36–56.76 Mbp) in this study and NP from another study (Bourke et al., 2018). Schulz et al. (2021) reported two candidate genes that are involved in the regulation of plant growth on LG2 between 65.6 and 68.7 Mbp. The first candidate gene LONGIFOLIA is known to be responsible for regulating cell elongation in Arabidopsis (Lee et al., 2006). The second, RhNAC100, was discovered in rose petals, and it was found to be an ethylene-inducible NAC transcription factor. A study showed that the suppression or silencing of RhNAC100 led to a substantial increase in petal size (Pei et al., 2013).

In addition, there was a cluster on LG3 between minor QTL for Diam across both populations (15.44–27.65 Mbp) and the second major QTL for NP (qNP.TX2WOB-LG3.1). The distal end of LG5 (85.60–85.70 Mbp) had QTL for DWT and Diam.

Overall, the co-localization of LG1 or LG3 between two or more traits may result from either a single locus with the pleiotropic effect or QTL clusters corresponding to a tight linkage between distinct loci. This has coincided with the positive correlations among the most studied traits.

Further studies of the genetic basis for these flower size traits are needed using larger and more diverse germplasm. In addition, phenotyping in multiple seasons/years is required for Diam, as it showed moderate G×E and petaloids. This would help to test QTL stability across environments and genetic backgrounds before being utilized in a breeding program. Additional work on these populations would be needed to identify the gene(s) responsible for Diam, flower weight, and petaloids.