Identification and QTL Analysis of Flavonoids and Carotenoids in Tetraploid Roses Based on an Ultra-High-Density Genetic Map

Roses are highly valuable within the flower industry. The metabolites of anthocyanins, flavonols, and carotenoids in rose petals are not only responsible for the various visible petal colors but also important bioactive compounds that are important for human health. In this study, we performed a QTL analysis on pigment contents to locate major loci that determine the flower color traits. An F1 population of tetraploid roses segregating for flower color was used to construct an ultra-high-density genetic linkage map using whole-genome resequencing technology to detect genome-wide SNPs. Previously developed SSR and SNP markers were also utilized to increase the marker density. Thus, a total of 9,259 markers were mapped onto seven linkage groups (LGs). The final length of the integrated map was 1285.11 cM, with an average distance of 0.14 cM between adjacent markers. The contents of anthocyanins, flavonols and carotenoids of the population were assayed to enable QTL analysis. Across the 33 components, 46 QTLs were detected, explaining 11.85–47.72% of the phenotypic variation. The mapped QTLs were physically clustered and primarily distributed on four linkage groups, namely LG2, LG4, LG6, and LG7. These results improve the basis for flower color marker-assisted breeding of tetraploid roses and guide the development of rose products.


INTRODUCTION
Roses are one of the most important ornamental plants worldwide, with a high value in various areas such as cut flowers, gardens and the food and medical industries. The genus Rosa contains approximately 200 species, but only 7-10 of them have contributed to the breeding of modern roses (Scariot et al., 2006). After a long history of hybridization, most modern roses have valuable characters that have formed from a complex genetic background at the tetraploid level. Currently, the number of modern rose cultivars has reached approximately 40,000 in number. However, the inheritance patterns of many rose traits remain unclear owing to their high heterozygosity, various ploidy levels and extreme intrageneric variation (Liorzou et al., 2016).
Flower color is one of the most significant traits of ornamental plants, and it also plays an important role in improving the aesthetic and economic values of roses. The appearance of flower color is largely affected by the various pigments synthesized by plants. The types of these pigments primarily include flavonoids, carotenoids and betalains, with flavonoids and carotenoids contributing the most to formation of flower color. Flavonols and anthocyanins are particularly important classes of flavonoids. Flavonols are highly prevalent in flower petals of various colors and play key roles in attracting insect pollinators. Flavonols primarily turn plants to yellow, while anthocyanins can provide orange, red, pink and blue colors to plants (Tanaka et al., 2008). Carotenoids primarily determine the yellow and orange colors of flowers and are also essential to human nutrition and health (Mayne, 1996).
Because of the complicated genetic background in roses, understanding the process of pigment synthesis in roses is particularly challenging. However, molecular markers and genetic mapping can provide an effective way to analyze the mechanism by which multiple genes affect variation in pigment regulation at the population level. Substantial progress has been made in the genetic mapping of tetraploid roses.
The first tetraploid rose genetic map was constructed by Rajapakse et al. (2001) using amplified fragment length polymorphism (AFLP) and simple-sequence repeat (SSR) markers. Over the last 20 years, major tetraploid populations include the double pseudo test cross population GGFC, which was used to study flower morphological traits based on AFLP and SSR markers (Guterman et al., 2002;Gar et al., 2011). The K5 population, developed from tetraploid cut roses, has been used to study the inheritance of multiple traits, including prickle density (Koning-Boucoiran et al., 2012), petal number (Koning-Boucoiran et al., 2012;Bourke et al., 2018a), resistance to powdery mildew (Koning-Boucoiran et al., 2009, 2012 and anthocyanin metabolism (Gitonga et al., 2016). The YS population, which was derived from the ancient rose 'Yunzheng Xiawei' and the modern rose 'Sun City' (Yu et al., 2015), has been shown to segregate for petal number, flower size and flower color. Lastly, a tetraploid climbing rose population (Zurn et al., 2018) has been used to study resistance to black spot.
High-throughput sequencing technologies are powerful tools that can greatly increase the density and quality of genetic linkage maps. Whole-genome sequencing (WGS) technology can be used to compare each individual genome to a reference genome and thus provide considerable nucleotide sequence information in a short time. Thus, WGS tools enable researchers to analyze the genetic mechanism of traits and finely map important genes at the genomic level (Smulders et al., 2019). These approaches have the advantages of being fast, efficient and high yielding, and after years of development, the associated costs have been dramatically reduced. In recent years, several rose high-density genetic maps based on genotyping sequencing and single nucleotide polymorphism (SNP) array technology have been constructed (Vukosavljev et al., 2016;Bourke et al., 2017;Yan et al., 2018). However, compared with these techniques, WGS can develop much more data than SNP arrays and thus, obtain more comprehensive genomic information. As a powerful tool, WGS is widely used in scientific research. Its enormous sequencing datasets provide a solid basis for analyzing individual differences, population evolution and genetic mechanisms.
Over the last 20 years, the number of marker types and overall marker density have increased substantially in rose genetic maps, with the total genome coverage also increasing (Vijay et al., 2021). However, the breeding and production of modern roses are primarily conducted with tetraploid cultivars. Differences in the ploidy level remain a major issue in rose QTL mapping, and it has resulted in low efficiency of applications based on diploid maps. Thus, with the rapid development of sequencing technology and analysis methods, QTL mapping at the tetraploid level poses significant advantages.
In this study, an ultra-high-density tetraploid rose genetic map was constructed using whole genome resequencing technology. Based on this map, we performed QTL analysis with pigment content data of flavonoids and carotenoids. Our study combined genomics and metabolomics to finely map flower-color-related loci to a genetic linkage map of tetraploid roses and analyze the genetic factors that influence pigment content in rose petals.

Mapping Population
The YS population used in this study was previously described by Yu et al. (2015). The ancient Chinese garden rose cultivar 'Yunzheng Xiawei' (YX) was used as the female parent, and the modern rose cultivar 'Sun City' (SC) was used as the male parent. Both parents are tetraploid. A total of 187 individuals were selected from among the F 1 progeny. Both parents and the resulting progeny were planted at the Xiao Tangshan Experimental Base of the National Engineering Research Center for Floriculture, Beijing, China. The parents and the YS population were grown under the same environment with regular irrigation and fertilization.

Pigment Content Extraction and Analysis
The petals of parents and offspring were collected at full bloom, immediately frozen in liquid nitrogen and stored at −80 • C. The flavonoids and carotenoids were extracted by HPLC and LC-MS as described by Wan et al. (2018).

Extraction and Analysis of Flavonoids
A total of 0.15-0.25 g samples were weighed, ground and then extracted in 0.9 mL of solvent (methanol/water/formic acid/trifluoroacetic acid, 70:27:2:1, v/v/v/v) in an ultrasonic cleaner (KQ-250DA, Jiangsu, China). The HPLC analyses were conducted using a Waters 2695 HPLC system connected with a 996-photodiode array detector (Waters, Milford, MA, United States). The supernatant was injected into an XBridge BEH C18 reversed-phase column (150 mm × 4.6 mm, 2.5 µm, Waters). The mobile phase was composed of 0.5% aqueous formic acid (A) and acetonitrile (B). The gradients were programmed as follows: 0 min, 5% B; 5 min, 10% B; 30 min, 19% B; 50 min, 40% B; 50.01-60 min, 5% B. The column temperature was 25 • C, the detection wavelengths of anthocyanin and flavonol were 520 and 350 nm, respectively. All samples were extracted in triplicate. LC-MS analysis was performed on the same HPLC system as described above using the parameters described by Wan et al. (2018).
Two-milligram flavonoid standards (cyanidin 3,5-diglucoside, pelargonidin 3,5-diglucoside, quercetin 3-O-glucoside, kaempferol 3-O-rutinoside, kaempferol 3-O-glucoside, and kaempferol) were diluted to six concentration gradients. The peak areas corresponding to each concentration of the flavonoid standards were recorded, and standard curves were drawn to obtain linear regression equations. Then, the peak area of the identified flavonoid pigment was substituted into the regression equation, and the content of the corresponding substance in the solution was calculated.

Extraction and Analysis of Carotenoids
After the samples were quickly ground, 0.2 g of each sample was weighed and promptly added to 4 mL of methanol, which was then shaken for 20 min. Afterward, 2 mL of NaCl solution (10%, w/v) was added, followed by 30 s of manual shaking. The organic phase in the upper layer was collected after the samples were allowed to stand. The mixed organic phase was then dried and saponified. A volume of 1 mL of NaCl solution (10%, w/v) was added to the saponified sample, followed by 30 s of manual shaking. Samples were allowed to stand for 5 min and then incubated with 1 mL of n-hexane/ether (3:1, v/v) for 1 h, after which the supernatant was finally evaporated to dryness and stored at −80 • C. All of the extraction processes were conducted in darkness or yellow light, and the temperature did not exceed 23 • C. The dried carotenoids were dissolved in 1 mL of organic solvent (methanol/methyl tert-butyl ether = 1:1, v/v), which was then filtered and used for HPLC analysis. Three mobile phases were used: A, methanol; B, methyl tert-butyl ether; C, ultrapure water. The gradients were set as follows: 0-2 min, 95% A + 5% C, 10 min, 95% A + 3% B + 2% C; 21 min, 95% A + 5% B; 27 min, 90% A + 10% B; 37 min, 70% A + 30% B; 40 min, 50% A + 50% B; 40.01-50 min, 95% A + 5% C. The injection volume was 20 µL; the column temperature was 25 • C, and the detection wavelength of carotenoid was 450 nm. The mass spectrometry information was obtained by HPLC-microTOF-Q. All the samples were extracted in triplicate. An LC-MS analysis was performed on the same HPLC platform as noted above, and the parameters were the same as those described by Wan et al. (2018). Carotenoid [all-E]-antheraxanthin, [all-E)]-zeaxanthin and [all-E]-lutein) were weighed and diluted with 10 mL of organic solvent (methanol: methyl tert-butyl ether = 1:1) to prepare standard solutions with six concentration gradients. The peak area corresponding to each concentration of the five types of carotenoid standards was recorded, and standard curves for the five types of standards were drawn. Corresponding linear regression equations were then obtained. The peak area of the identified carotenoids was substituted into the regression equation to calculate the content of each substance in the test solution.

DNA Extraction
DNA was extracted from young leaves of the parents and the YS population using a Plant Genomic DNA Extraction Kit from TIANGEN (Beijing, China) according to the manufacturer's instructions. Total DNA was extracted from each of three biological replicates. The concentration of DNA was measured using an ultraviolet spectrophotometer. After extraction, all DNA samples were stored at −20 • C.

SNP Detection and Genotyping
Whole-genome resequencing was used to detect genome-wide variants. After detection of the samples, the genomic DNA was segmented, and the qualified gDNA was selected for library construction. Following fragmentation, a 2100 Bioanalyzer platform (Agilent Technologies, Santa Clara, CA, United States) was used for quality control, and the library was prepared after end repair and 3 A-tailing. After PCR amplification, the insert size was detected again using the Agilent 2100 Bioanalyzer. Finally, the high-quality samples were sequenced on an Illumina NovaSeq PE150 platform (Illumina, San Diego, CA, United States) using standard Illumina protocols.
The raw data were filtered as follows: (1) Cutadapt v2.3 (Martin, 2011) was used to cut the reads with N base calls at both ends or low-quality bases (quality score <20); (2) the reads with more than 5% N base calls were removed; (3) reads with more than 50% bases with a quality score <20 were removed. High-quality clean data were obtained after the quality check and filtering. Clean reads were then compared to the reference genome sequence of Rosa chinensis 'Old Blush' (Hibrand Saint-Oyant et al., 2018) using BWA software (Li and Durbin, 2010). SNPs were detected using Sentieon Haplotyper (Weber et al., 2016) based on the alignment results. The SNP markers were then genotyped under autotetraploid mode. Each locus was converted into a dosage value. The dosage values corresponding to genotypes of AAAA, AAAB, AABB, ABBB, and BBBB were 0, 1, 2, 3, and 4, respectively, with the genotypes of both parents and the offspring converted according to this rule.

Linkage Map Construction
In addition to the SNP markers obtained by whole genome resequencing, we also added simple sequence repeats (SSRs) (Yu et al., 2015), and the SNPs that were previously developed by specific locus amplified fragment sequencing (SLAF-seq) (Yu et al., 2021) to the current map. Linkage map construction was performed using the polyploid mapping R package polymapR (Bourke et al., 2018b).
The dosage-scored data were imported into polymapR. The markers were filtered according the following exclusion criteria: (1) non-segregating markers; (2) markers with more than 5% missing data; and (3) samples with more than 1% missing data.
The parental linkage maps were first constructed, and then the consensus map was integrated afterward. The 'cluster_SN_markers' function was first used to cluster the (1,0) markers. According to the clustering result, we chose to add (2,0) markers and use the 'bridgeHomologues' function to obtain four clusters that corresponded to the four homologs of each chromosome. The remaining types of markers were assigned to each homolog and were ordered using the MDSMap package (Preedy and Hackett, 2016). After the parental maps were constructed, all the homolog information from two parents was combined, and the markers were ordered again to construct the integrated map. Finally, the integrated genetic map was used to generate phased linkage maps to conduct QTL analysis.

QTL Analysis
The QTL analysis was conducted using composite interval mapping (CIM) with the R package polyqtlR (Bourke et al., 2018c). Phased maps and dosage marker information were imported into polyqtlR, together with the pigment content data. Identity-by-descent (IBD) probabilities for the population were estimated using polyOrigin (Zheng et al., 2020). Based on the IBD results, an initial QTL scan was performed with each set of phenotypic data. Significance thresholds of LOD values were determined using a genome-wide permutation test with 1,000 permutations (α = 0.05). The initial QTL results were outputted with peak positions and LOD values. The already-identified QTL positions were then included as genetic co-factors to reduce the amount of unwanted background variance. Further rounds of QTL scans were performed to determine other QTL positions and to check if multiple peaks within a single linkage group were indeed independent. The QTL results were confirmed and outputted when no new peaks could be identified. The visualization of QTL results was performed by polyqtlR in the R statistical computing environment (R Development Core Team, 2016). Markers located at the LOD peaks or closest to both ends of the QTL intervals were considered to be closely linked to the QTLs. The naming of QTLs was based on previously published guidelines (McCouch et al., 1997).
Eight characters related to anthocyanin content were recorded to perform further analyses ( Table 1). There were extreme values of pigment content segregation in the F 1 population, and the skewness of each of the eight characters was greater than 1, indicating a positive skew distribution. Further analysis of the frequency distribution showed that all the traits exhibited continuous variation in the mapping population (Figure 1), indicating these traits were quantitative traits controlled by polygenes and therefore suitable for QTL analysis.

Flavonol Content Variation
The flavonols in the petals of parents and F1 progeny were determined; the results are shown in (Table 2), and the frequency distribution is shown in (Figure 2). The contents of 20 flavonols and total quercetins, total kaempferols and total flavonols showed transgressive segregation. The coefficient of variation of 23 traits related to flavonol contents ranged from 44.85% to 191.52% and showed continuous variation in the F1 population, indicating that the traits were suitable for QTL analysis.

Carotenoid Content Variation
Carotenoids in the petals of parents and F 1 progeny were also determined ( Table 3). There were 16 carotenoids detected in the TABLE 1 | Statistical analysis of anthocyanins in the petals of parents and F 1 progeny.

Traits
Compound   petals of the male parent SC, and no carotenoids were found in the female parent YX. A total of 16 carotenoid compounds were identified in the F1 progeny. The carotenoid contents in the petals of different progenies ranged from 0.00 to 467.12 µg/g FW, and the coefficient of variation ranged from 166.86 to 821.44%. In total, 21 traits were recorded for further analysis.
The frequency distribution showed that all 21 traits exhibited continuous variation in the mapping population (Figure 3), indicating that these traits were quantitative traits controlled by polygenes and therefore suitable for QTL analysis.

Sequencing Data and SNP Markers
In total, 33.48 and 43.85 Gb of raw data were obtained for the paternal and maternal lines, respectively. On average, 12.27 Gb of raw data were obtained for each offspring. The sequencing data were of sufficient quality, with an average Q30 percentage of 90.95%. The mapping rate of the male and female parents were 97.79 and 98.18%, respectively, indicating high-quality mapping results. After variant calling, a total of 16,038,514 initial SNPs were obtained from whole genome resequencing. These SNPs were first filtered, and those that were homozygous, had missing data or were of insufficient sequencing depth were removed. However, there were still too many SNPs remaining after screening (10,467,012). Thus, we used a sampling method to divide each chromosome into 10,000 regions according to the number and order of loci. Intermediate sites were selected within each region. Finally, 10,000 markers were obtained for each of the seven chromosomes. The pre-developed 440 SSR markers and 17,257 SNP markers obtained by SLAF-seq were then added, genotyped and imported into polymapR. During the mapping process, the markers were further filtered following the criteria within polymapR, and 78,438 markers were thus eliminated.  Finally, a total of 9,259 high quality markers were obtained and mapped successfully. Nine different marker segregation types were identified, namely, 1 × 0, 0 × 1, 2 × 0, 0 × 2, 1 × 1, 1 × 2, 2 × 1, 1 × 3, and 2 × 2 (Figure 4). Type 1 × 1 segregation, spanning 2,746 markers, was the most common marker type, comprising 29.67% of all the markers mapped. The single-dose markers, type 1 × 0 and type 0 × 1 segregation patterns corresponding to 1,563 and 1,728 markers, respectively, together comprised 35.54% of all the markers. Type 2 × 2 and type 1 × 3 segregation patterns

Linkage Map Construction
The tetraploid rose genetic linkage map was constructed using polymapR with a total of 9,259 markers, including 175 SSRs and 9,048 SNPs. The final linkage map contained seven linkage groups (LGs), which is consistent with the number of haploid chromosomes observed in rose. The detailed map information is shown in Table 4 and Figure 5. The total length of the integrated map was 1285.11 cM, with an average distance of 0.14 cM between markers, indicating the high density of the map. The length of each LG ranged from 169.40 to 205.20 cM. LG2 and LG3 had the most markers, 1,669 and 1,653, respectively, with a minimum average distance of 0.11 cM.
LG7 was the smallest linkage group, which contained a total of 949 markers separated by an average distance of 0.19 cM. Between two adjacent markers, the average percentage of gaps (≤5% total map length) was 99.9%, indicating the high coverage and even marker distribution of the linkage map.

QTL Analysis
A total of 67 anthocyanins, flavonols and carotenoids were analyzed by a whole-genome QTL analysis, and the detailed information is presented in Table 5. Significant QTLs were detected for 33 compounds, including four anthocyanins, 18 flavonoids and 11 carotenoids. Across the traits, 46 QTLs were detected, explaining 11.85-47.72% of phenotypic variation. The QTLs were primarily distributed on four linkage groups: LG2, LG4, LG6 and LG7 (Figure 6). Among them, the number of QTLs detected on LG6 (20) was the highest, with all of them contributing to variation in the flavonoid content; the numbers of QTLs on LG4 and LG7 were 15 and 8, respectively, and most of them were associated with variation in the carotenoid content. Three QTLs were identified on LG2. Four QTL were detected for the anthocyanin content, all of which were distributed on LG6. The phenotypic variation rate ranged from 31.31 to 45.97%, which shows that all of them were major QTLs (i.e., contributing to more than 10% of the variation). The confidence intervals of the four loci highly overlapped. Among them, qA9-6-1 had the highest contribution rate, accounting for 45.97% of variation. The molecular marker closely linked to qA9-6-1 was Chr6_36524322. Additionally, 23 QTLs were detected for flavonol contents in the three linkage groups. Among them, 16 QTLs were found on LG6. On LG2 and LG4, there were three and four QTLs, respectively. The phenotypic variation rate ranged from 11.85 to 47.72%. Among the 23 QTLs, 12 and eight QTLs explained more than 20 and 30% of phenotypic variation, respectively. Among the detected loci, qF15-6-1 had the rate of highest contribution, and its closely linked molecular marker was Chr6_16827958. A total of 19 QTLs were identified to be associated with the content of carotenoid, which were distributed in two linkage groups. The number of QTLs on LG4 and LG 7 were 11 and eight, respectively, and they explained from 14.45 to 26.38% of the phenotypic variation. The 19 loci had highly overlapping confidence intervals. Most had the same QTL peaks and linked markers, indicating that many carotenoids were controlled by the same QTLs. Among the loci detected, qC24-4-1 had the highest contribution rate to phenotypic variation, and its closely linked molecular marker was Chr4_56663092.

Tools for Genetic Studies in Polyploids
Polyploid species follow a very complicated pattern of meiosis that differs from the pattern in diploid organisms. Therefore, it is very difficult to conduct genetic studies of polyploids. As more research in polyploid species is conducted, the necessary tools are also being developed, thus, rendering the genetic analysis of polyploid populations progressively more reliable.
The number of tools for polyploid linkage mapping is continuously increasing, from the well-known TetraploidMap (Hackett et al., 2007), to recently developed R packages, such as PERGOLA (Grandke et al., 2017), polymapR (Bourke et al., 2018b), and MApoly (Mollinari and Garcia, 2019), among others. More tools have been developed to focus on a greater number of population types and multiple ploidy levels. For example, PolymapR, which was used in this study, can perform linkage analysis for polysomic triploids, tetraploids and hexaploids, as well as segmental allotetraploid populations (Bourke et al., 2018c). The polyqtlR R package is a reliable tool for QTL analysis using identity-by-descent (IBD) probabilities in F 1 populations of outcrossing in heterozygous polyploid species. With the continual analysis of polyploid genetic model systems Sun et al., 2020), the software tools for polyploid genetic research will also provide a better understanding of polyploid inheritance.

High Density Linkage Map Construction
Linkage map construction for roses has primarily been focused at the diploid and tetraploid levels. The first rose genetic linkage map was constructed by Debener and Mattiesch (1999) using a diploid Rosa multiflora population. Subsequently, more diploid rose maps using different markers have been published (Debener et al., 2001;Crespel et al., 2002;Dugo et al., 2005;Yan et al., 2005Yan et al., , 2018Linde et al., 2006;Hibrand Saint-Oyant et al., 2008;Remay et al., 2009;Spiller et al., 2011;Li et al., 2019;Lopez Arias et al., 2020). Studies of tetraploid populations are relatively less common than those of diploid rose populations. Therefore, our study aimed to increase the quality of tetraploid rose genetic map by improving the map length, density and coverage, thus, providing a solid foundation for the fine mapping of quantitative traits, map-based gene cloning and functional gene mining of tetraploid roses in the future.
Most of the published linkage maps of roses are based on traditional molecular markers, such as AFLPs, restriction fragment length polymorphisms (RFLPs) and SSRs. However, markers developed by sequencing technology have gradually replaced traditional markers. With the advantages of higher coverage and substantially higher marker number, the length and density of genetic maps can be significantly improved. Yan et al. (2018) generated the first diploid rose genetic map based on genotyping by sequencing technology data. They mapped 20 SSRs and 3,507 SNPs onto a total map length of 892 cM, with an average gap of 0.25 cM. However, they used Fragaria vesca as the reference genome, and that may have reduced the accuracy of SNP calling. Li et al. (2019) constructed a diploid rose genetic map using RAD-seq technology. The map was improved in terms of its length, which was 1027.4 cM, but the marker number and map density decreased, with an average distance of 0.96 cM between adjacent markers. Compared with reduced-representation genome sequencing (RRGS), wholegenome resequencing has a wider genome coverage and can thus, obtain more comprehensive genomic information. WGS technology has been widely used in crops, such as peanut (Agarwal et al., 2018), watermelon (Guo et al., 2019), cucumber , and others. However, this approach has not been widely used in rose research.
In previous research, the YS population was used to construct a genetic map using SLAF-seq technology (Yu et al., 2021) ( Table 6). A total of 6,842 markers were contained in a 1,158.90 cM map, with an average distance of 0.18 cM between markers. However, the maximum gap was 29.16 cM, indicating that the distribution of markers was uneven. Based on the previous map, we made further improvements using WGS technology. In this new ultra-high-density genetic map, a total of 9,259 molecular markers were well-distributed over seven LGs. The map has achieved substantial improvements in both length and density, with a total length of 1,285.11 cM and an average marker spacing of 0.14 cM. This improved map can serve as a strong tool for rose genetic research, as well as the marker-assisted breeding of tetraploid roses.

Mapping Flower Color Related Traits
The improvement of flower color is of substantial significance to the ornamental value and economic value of roses. A substantial amount of research has been conducted on pigment contents, compositions and variations in roses. Schmitzer et al. (2009) studied the pigment composition of modern rose petals and found that the main anthocyanins were cyanidin 3,5dioxyglucoside and geranium 3,5-dioxyglucoside, and the secondary anthocyanins were cyanidin 3-O-glucoside and geranium 3-O-glucoside. These four anthocyanins were also detected and mapped in F1 population in this study, and all of them were located on LG4. The flavonols detected in this study, kaempferol 3-O-rutinoside, kaempferol 3-O-glucoside and kaempferol, are widely found in the petals of R. chinensis, Rosa gallica, Rosa Rugosa, and Rosa damascena (Kumar et al., 2009;Sarangowa et al., 2014). This proved that the results of this study are precise and accurate. Compared with anthocyanins and flavonols, studies on carotenoids in Rosa petals have been studied less frequently.
In this study, we detected a total of 46 QTLs related to pigment content traits based on the newly developed linkage map. The phenotypic variation associated with these QTLs was greater than 10% for all the traits, indicating that all of them were major QTLs. The LOD thresholds obtained from the QTL analysis were high, ranging from 6.20 to 13.26. This may result in the omission of some QTLs with smaller effects but can effectively reduce the detection of false-positive QTLs, which is necessary for further fine mapping and mapbased cloning.
Many studies have mapped flower color traits of diploid and tetraploid rose populations, but most have focused on chromatism. Few studies have analyzed anthocyanins and mapped major QTLs across several linkage groups. Henz et al. (2015) analyzed total anthocyanins in diploid rose progeny under different environments, detecting two major QTLs on ICM2 and ICM6 in all the environments. A marker derived from a bHLH (basic helix-loop-helix) gene was mapped to the center of ICM6, which corresponds to our current results (Figure 6). Gitonga et al. (2016) also detected QTLs on ICM6 in the tetraploid K5 population. In our study, all of the variation in content of anthocyanin compounds was mapped to LG6, with peak positions around 65 cM. This is consistent with previous studies and showed that bHLH is significantly associated with anthocyanin biosynthesis. Moreover, most flavonols also mapped to LG6, and the confidence intervals overlapped with the range of anthocyanins. It is assumed that these traits are controlled by TABLE 6 | Comparison between maps constructed by the specific locus amplified fragment sequencing (SLAF-seq) and whole-genome sequencing (WGS) methods across linkage groups (LG).

Method
Number of markers Total distance (cM) Average distance (cM) Max gap (cM) LG1 LG2 LG3 LG4 LG5 LG6 LG7 Total the same genes in the flavonoid metabolism pathway and have similar genetic bases. Few studies have mapped loci associated with carotenoids in roses. Schulz et al. (2016) presented a genome-wide association study (GWAS) for tetraploid roses and mapped anthocyanin and carotenoid contents. However, because of the lack of a published rose genome, they were limited to mapping the traits to the Fragaria and Prunus genomes. In early studies, Yu et al. (2021) used a previous map developed by SLAF-seq technology to locate QTLs associated with carotenoid contents and found that there were QTL clusters on both LG4 and LG7 (Table 6), which is also consistent with the current results. However, in previous studies, QTLs were analyzed by interval mapping. Thus, the results were not as accurate or comprehensive as the results of composite interval mapping in this study.
QTLs are easier to map to the same region if the traits are highly correlated with each other (Guan et al., 2015). The QTLs mapped in this study were clustered, and the QTLs for different carotenoids were located in roughly the same regions of LG4 and LG7. Santos and Simon (2002) also found that multiple QTLs that controlled the content of β-carotene and ζ-carotene were clustered together, which may be owing to the fact that related phenotypic traits are regulated by one or a few pleiotropic gene. Thus, they can effectively improve the utilization of genes and reduce the loss of favorable allele combinations through recombination.
In this study, the tetraploid rose genetic map with the highest density was constructed by combining genomic and metabolomic techniques, and QTL mapping of metabolites related to flower color was conducted. Integrating genomics, metabolomics and other omics data to analyze the synthesis pathways of secondary metabolites can assist modern bioengineering technology to improve the content of related pigments in ornamental plants and promote the comprehensive utilization of edible and medicinal plant materials. The results of this study provide the basis for the fine mapping of flower color genes in roses and thus, suggest new objectives for the production of natural metabolites and the development of rose products.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BC performed the efinal data analysis and wrote the manuscript. HW collected and analyzed the phenotypic data. YH, CY, and LL helped to improve the method and edited the manuscript. CY, HP, and QZ participated in supervision and helped to revise the manuscript. All authors contributed to the article and approved the submitted version.