Genome-wide identification of quantitative trait loci for important plant and flower traits in petunia using a high-density linkage map and an interspecific recombinant inbred population derived from Petunia integrifolia and P. axillaris

Petunia is a very important flower in the global floriculture industry and has played a critical role as a model in plant genetic studies. Owing to limited genetic variability in commercial germplasm, development of novel petunia phenotypes and new varieties has become increasingly difficult. To enrich petunia germplasm and facilitate genetic improvement, it is important to explore genetic variation in progenitor species that may contain highly valuable genes/alleles. In this study, an interspecific recombinant inbred population (168 recombinant inbreds) derived from Petunia integrifolia × P. axillaris were phenotyped for days to anthesis (DTA), flower count (Flower_C), flower diameter (Flower_D), flower length (Flower_L), plant height (Plant_H), plant spread (Plant_S), and plant size (Plant_Z) in 2014 and 2015. Transgressive segregation was observed for all traits in both years. The broad-sense heritability on a 2-year basis varied from 0.38 (Flower_C) to 0.82 (Flower_L). Ten QTL were consistently identified in both years and by two mapping strategies [multiple QTL mapping (MQM) in MapQTL and inclusive composite interval mapping (ICIM) in IciMapping]. Major QTL explained up to 30.2, 35.5, and 47.1% of the total phenotypic variation for Plant_S, Flower_L, and Flower_D, respectively. These findings should be of significant values for introgression of desirable genes from wild petunias into commercial varieties and future genetic improvement of this important flower.


Introduction
Cultivated flowers serve a very important role in human life and health, the global economy, and the beautification and protection of the environment. Flower production has become one of the most dynamic and sophisticated sectors of the global horticulture industry. It was estimated that the worldwide production value of cultivated flowers or floricultural production reached 60 billion dollars in 2003 1 . Since then, further growth has occurred in almost every continent. It is projected that the global floriculture market will continue to grow at a compound annual growth rate of 5.4% over the period from 2016 to 2020 2 .
To sustain the global flower production industry, continuous introduction of new cultivars with improved or novel characteristics is essential. Towards this, plant breeders constantly seek to identify novel genes/alleles and combine them into new or improved cultivars. In many crops, including widely produced and used flowers and other ornamental plants 3 , the lack of genetic diversity and lack of novel genes/alleles in the commercial germplasm pool have been limiting plant breeders' progress in genetic improvement and new cultivar development 4,5 .
Identification and utilization of desirable genes/alleles from wild or progenitor species have been suggested as an effective approach to overcoming this limitation. Enormous efforts have been made in some major agronomic and horticultural crops to characterize wild and ancestor germplasm and identify favourable genes/alleles from the germplasm through phenotyping, genetic mapping, and introgression 6 . On the other hand, wild and ancestor species often perform poorly in horticultural aspects compared to elite germplasm 7 . Wild accessions may carry undesirable genes for the improvement of commercial cultivars 8,9 . When a wild species is crossed with an elite cultivar, the inferior alleles can be simultaneously dragged into cultivars, reducing the plant performance of new cultivars 10 . Numerous rounds of backcrossing are required to reduce or eliminate the inferior donor alleles from elite cultivars, which is a laborious and timeconsuming process. Genetic mapping, identification of quantitative trait loci (QTLs), and marker-assisted selection have been used to facilitate the introgression of desirable alleles from wild species to elite cultivars 11 . Over the past two to three decades, several molecular marker systems have been used in such efforts. Genotyping by sequencing (GBS) is a recently developed strategy for large-scale marker discovery 12 . It has been made possible by rapid advances in next-generation sequencing technology. With this strategy, it is possible to sequence hundreds of barcoded samples in a single sequencing lane simultaneously and to reveal single-nucleotide polymorphism (SNP) sites throughout the whole genome. The high output of SNP discovery by GBS has greatly facilitated the construction of high-density, high-resolution genetic linkage maps. GBS has been widely used to construct high-coverage linkage maps and conduct QTL analyses in multiple important agronomic crops 13 .
Garden petunia (Petunia hybrida) is a very important flower in the global floriculture production. It is cultivated all over the world and is one of the most important Solanaceae utilized for ornamental purpose 14 . Garden petunia is often among the most popular flowers planted in outdoor gardens in many countries 15 . In the United States, it ranks first in wholesale value among annual bedding plant crops 16 . Cultivated petunia originated from the cross between Petunia axillaris and Petunia integrifolia 17 . As petunias have been commercially bred with limited germplasm sources for the past 150+ years, the genetic diversity among current commercial cultivars has been low, resulting in high similarities among commercial cultivars and loss of some useful traits 8,18,19 . Several studies have indicated that wild P. axillaris and P. integrifolia carry traits that may be beneficial to commercial petunia, such as faster development rates 8 , superior freezing tolerance 19 , longer flower longevities 20 , or arthropods resistance 21 . Consequently, interest in introgressing traits from progenitor species to elite petunia cultivars has been strong 21 .
Several genetic linkage maps were developed in petunia, using restriction fragment length polymorphism markers, amplified fragment length polymorphism markers, and simple sequence repeat (SSR) markers 9,17,[22][23][24] . These genetic maps have been used to identify QTL for pollination syndrome traits (length of pistil, stigma, and corolla tube; flower scent; and corolla diameter). A SNPbased linkage map was recently reported in petunia and employed to identify QTL controlling petunia plant development rates (as well as the number of branches and flower buds and days to anthesis (DTA)) under varying temperatures 25 . All reported petunia QTL studies were conducted in the greenhouses using container-grown plants. No or few QTLs have been reported for important aesthetic traits in petunia, including plant size (Plant_Z) and flower count (Flower_C).
In this study, we (1) characterized and phenotyped seven important plant and flower aesthetic traits (DTA, Flower_C, flower diameter (Flower_D), flower length (Flower_L), plant height (Plant_H), plant spread (Plant_S), and Plant_Z) in an open field using a recombinant inbred population derived from a cross between P. integrifolia × P. axillaris for 2 consecutive years in 2014 and 2015, (2) estimated the heritability for these traits, and (3) identified and located QTL controlling these traits using a highdensity SNP bin map developed by the GBS technology.

Phenotypic value
Phenotypic data including mean value, mid-parents value, and data range for DTA, Flower_C, Flower_D, Flower_L, Plant_H, Plant_S, and Plant_Z of the parents and their recombinant inbred lines (RILs) and the broadsense heritability (H 2 ) estimate for each trait based on combined 2-year data are presented in Table 1. Plants of P. axillaris

Heterozygosity retention in RILs and at various marker loci
We calculated the heterozygosity level in each RIL based on the number of heterozygous marker loci out of the total number of marker loci analysed (518). The heterozygosity level in the 168 RILs ranged from 0 to 21.80% (Fig. 4), averaged to 3.03%, which is much higher than the expected 0.78% (0.5 7 ) for an F 7 RIL population. Overall,~80% of the RILs in the population had a heterozygosity level below 4.86%, and~90% of the RILs had heterozygosity level below 7.39%. We also calculated the heterozygosity level for each of the 518 marker loci by dividing the number of heterozygous RILs by the total number of RILs analysed (168). The percentage of heterozygosity per locus ranged from 0% to 23.67%. The average heterozygosity level for the 482 SNP marker loci was 2.70%, while the average heterozygosity for the 36 SSR marker loci was at least one-fold higher, reaching 6.64%. There were 12.74% of the marker loci (66) that had 5-10% heterozygosity, and 1.16% of the marker loci (6) had heterozygosity above 10%. These six loci are all of the SSR type. To show a genome-wide landscape of the heterozygosity residues, the genetic position of each marker locus with its heterozygosity level was plotted on the genetic linkage map (Fig. 5). The retained heterozygosity was not evenly distributed within and among linkage groups (LGs). LG 2, LG 3-1, LG 6, and LG 7 seem to have more regions that retained higher heterozygosity than other LGs. Heterozygosity seemed to be higher in some telomeric regions including the region from 54 to 60 cM on LG 2, the region from 0 to 8 cM on LG 6, and the region from 0 to 12 cM on LG 7.

QTL detection by multiple QTL mapping (MQM) in MapQTL
A total of 17 significant QTL in five LGs were identified for the seven petunia traits (Table 3; Fig. 6).
Two QTLs, C-Flower_D-2.1 and qFlower_L2.1, seemed to have significant QEI effects (logarithm of the odds (LOD) > 3.0). The QEI at C-Flower_D-2.1 explained 7.38% of the phenotypic variation and the QEI at qFlower_L2.1 accounted for 6.92% of the phenotypic variation. QEI effects increased the Flower_D in 2014 by 0.34 cm and the Flower_L in 2014 by 0.28 cm.
Significant QTL × QTL interactions were observed for two traits, DTA and Plant_Z ( Table 5). The interaction between two putative loci, one at 5 cM on LG 2 and the

Discussion
Significant phenotypic variation was observed for all seven traits (DTA, Flower_C, Flower_D, Flower_L, Plant_H, Plant_S, and Plant_Z) in the F 7 RIL population derived from the cross between P. integrifolia and P. axillaris (Table 1 and Fig. 3). All traits exhibited certain degrees of transgressive segregation in the RIL population. These results suggest that P. axillaris and P. integrifolia possess a very different genetic background. This is probably because P. axillaris and P. integrifolia evolved from two different ecogeographic isolations 27 and have developed their own pollination syndromes, which have restricted natural gene flows between the two species 28 . Similar transgressive segregation was previously reported for Flower_D, Flower_L, and DTA in P. integrifolia × P. axillaris F 2 populations 8,9 . Traits with higher H 2 across different growing environments and/or seasons can be reliably selected in breeding. In this study, Plant_H (H 2 = 0.70), Flower_D (H 2 = 0.79), and Flower_L (H 2 = 0.82) were found to have high heritability. Similar results were reported in previous petunia studies 8,9,29 . Four traits, including Flower_C (H 2 = 0.38), DTA (H 2 = 0.40), Plant_S (H 2 = 0.58), and Plant_Z (H 2 = 0.58), have moderate H 2 , suggesting that these traits are more sensitive to environmental changes. Moderate H 2 was observed for DTA in previous studies 8,9 . Similar H 2 estimates were also reported for Plant_S, Plant_Z, and Flower_C in a P. axillaris × Petunia exserta RIL population 30 .
In Theoretically the average heterozygosity in an F 7 inbred population should be 0.78%, as the heterozygosity reduces by half for each cycle of inbreeding. In this study, we observed a much higher level of heterozygosity (3.03%) in the F 7 interspecific petunia population (Figs. 4 and 5). Higher levels of heterozygosity were also observed in several telomeric regions of the LGs. This phenomenon was observed in several previous studies 31,32 . The biological meaning for retaining these heterozygous segments in petunia remains to be understood. One hypothesis might be that these segments are important to petunia growth and development, and when they become homozygous or fixed, petunia plants may have reduced fitness.
To identify consistent, useful QTL in petunia, we phenotyped a large number of RILs over two growing seasons and replicated each RIL several times in each growing season. By using MQM in MapQTL, a total of 17 putative QTL controlling seven important petunia plant and flower traits were observed. Twelve of these QTL were also confirmed by ICIM in the IciMapping software. Ten   25 used the same population, linkage map, and mapping algorithm (MQM) as we did in this study to identify and locate QTL for the traits Flower_D and DTA. The distinct difference between Guo et al. 25 and this study was that the population was previously phenotyped in an artificial growing environment (greenhouses with precise temperature and light control) in a temperate climate while the population was phenotyped in open fields in a subtropical climate in this study. Guo et al. 25 detected the QTL FD2.1 at 30.49 cM to 31.78 cM on LG 2 for Flower_D. A QTL (C_Flower_D-2.1) was detected for petunia Flower_D in the present study by using the ICIM strategy. C_Flower_D-2.1 is located at 28.25 cM to 29.75 cM on LG 2, very close to FD2.1. We projected the confident intervals of FD2.1 and C_Flower_D-2.1 to the P. axillaris genome and found one common scaffold Peaxi162Scf0007 (Supplemental Table S1). This scaffold contained 149 genes (Supplemental Table S2); two of the genes (Peax-i162Scf0007g00237 and Peaxi162Scf01617g) are predicted to be involved in plant auxin synthesis. It is known that auxins play an important role in regulating plant flower initiation and organ growth 33,34 (Supplemental Table S3); these two genes may be considered as candidate genes for petunia Flower_D regulation. Two other genes (Peax-i162Scf0007g02848 and Peaxi162Scf0007g02859, coding for an F-box protein and the developmental-regulator-ULTRAPETALA, respectively) (Supplemental Table S3) are predicted to be involved in floral meristem initiation or expansion in petunia, tomato, or Arabidopsis 35,36 . These two genes may be considered as candidate genes for petunia Flower_D regulation as well.
One of the Flower_L QTL, qFlower_L2.1, was located at 43.74 cM (LG 2) when the 2014 phenotyping data were analysed but shifted to 46.03 cM (LG 2) when the 2015 phenotyping data were analysed. Significant QEI (LOD = 5.31) was also observed for this QTL. We suspect that this QEI might have played some role in the QTL peak position shift between years or environments. Comparing QTL from Guo et al. 25 with those from the present study led to the recognition that two pairs of QTL for DTA, qDTA1.1 and DTA1.1 25 and qDTA6.1 and DTA6.1 25 , were mapped onto the same LGs, but in different regions in the LGs (around 21 cM for qDTA1.1 and around 31 cM for DTA1.1; around 8 cM for qDTA6.1 and around 19 cM for DTA6.1). It remains to be determined whether these QTL represent different loci or have shifted positions between studies. The trait DTA had low H 2 and was prone to be influenced by environment conditions, and the two studies were completed in very different growing environments (greenhouses with very good temperature and light control in a temperate climate 25 vs. open fields in a subtropical climate). These factors might contribute to change of QTL peak position between studies.
Results from this study showed that P. integrifolia possesses beneficial alleles that can increase Flower_C, while P. axillaris carries favourable alleles for high Flower_C, increased Flower_L and Flower_D, and larger Plant_Z. These alleles are of significant value for introgression into commercial petunia cultivars. Gene introgression from wild species to commercial cultivars has been seldom reported in ornamental plants, but it has been practised frequently in tomato 37 , rice 38 , and other crop species. It is estimated that 30-50% of major QTL from wild species or progenitor species could be beneficial for commercial breeding 35 . The major and environmentally stable QTL identified in this study can be very useful for further genetic improvement of petunia. The QTL information could facilitate the development of traitassociated molecular markers and accelerate markerassisted transfer of favourable QTL from wild species to commercial cultivars while minimize dragging of chromosomal fragments containing deleterious genes (linkage drag) into commercial germplasm.
The three QTL-rich segments in LG 1, LG 2, and LG 3 (Fig. 6) indicate the presence of linkage among different alleles or pleiotropic alleles that control two or multiple traits. This phenomenon has been observed in other plant species, such as Brassica napus 39 , sorghum 40 , and sweet cherry 41 . In an F 7 P. axillaris × P. exserta population, the presence of QTL-rich chromosomal regions were also observed LG 1, LG 2, and LG 4 42 . These results indicated that QTL-rich chromosomal regions may be common in petunia and may have additional value for petunia breeding.

Plant material
The RIL mapping population consisted of 168 individuals and was developed by crossing P. integrifolia (PI 28546, from the USDA Ornamental Plant Germplasm Center, Columbus, OH) and P. axillaris (PI 28546; USDA Ornamental Plant Germplasm Center) and selfing their progeny for seven generations (F 7 ) following a single seed descent procedure. This mapping population was previously described by Guo et al. 25 . P. axillaris exhibits an apical dominance growth habit, long internodes, long floral tubes, and large floral limbs. In contrast, P. integrifolia has a creeping growth habit, short internodes, short floral tubes, and small floral limbs. Two weeks later, the germination trays with young seedlings were transferred to a greenhouse where the air temperature was maintained between 25°C and 30°C. After 12 days, six seedlings were individually transplanted to 72-cell trays (66.04 cm in length × 33.02 cm in width) filled with a commercial soilless substrate (Fafard ® 3B; Conrad Fafard, Agawam, MA, USA). The seedlings were grown in the same greenhouse until they were ready to be transplanted to ground beds. Seedlings were fertilized twice weekly using a commercial water-soluble fertilizer containing 15% (w/w) total nitrogen, 5% phosphate (P 2 O 5 ), and 15% potassium (K 2 O) (Peters ® Excel; Everris NA, Dublin, OH, USA). Two weeks later, all seedlings were moved to a shade house with 30% shade and kept there for 1 week to acclimate the seedlings to the outdoor environment. After the acclimation, four seedlings per RIL and parent were transplanted to mulched, raised ground beds in the UF/GCREC experimental farm (central Florida; N 27 o 45", S 82 o 13"). The ground beds were fumigated with Pic-Clor 60 (60% chloropicrin and 40% 1,3-dichloropropene) at 45 kg per 1000 m 2 1 month prior to transplanting. Transplanted petunia plants were irrigated with a drip irrigation system 30 min a day. Each plant received 8 g of controlled-release fertilizer Osmo-cote® (The Scotts Miracle-Gro Company, Marysville, OH, USA). During the petunia growing season (late February to mid-June), the daily average air temperature ranged from 11°C to 28°C in 2014 and from 6°C to 28°C in 2015. The total precipitation during the growing season was 42.39 cm in 2014 and 46.30 cm in 2015. The experiments in both years followed a randomized complete block design, with four replicates and one RIL plant per experimental unit.

Collecting phenotype data
In each year, the large RIL population as well as their parents was phenotyped for seven plant and floral traits of the most significance to the use of petunia as a bedding and garden plant, including DTA, Flower_C, Flower_D, Flower_L, Plant_H, Plant_S, and Plant_Z. DTA were calculated as follows: DTA (days) = date of first flower anthesis − date of seed sowing. When approximate 50% of progeny came into flowering, all flowers on each plant were counted weekly for 7 weeks to obtain Flower_C data. Three fully opened flowers per RIL and parent were randomly selected to collect data for Flower_D and Flower_L. Flower_D was measured from one side of the flower's petal to the opposite side; Flower_L was measured from the base of the flower's calyx to the top of the flower's corolla. Plant_H and Plant_S were recorded near the end of the growing season (early June to mid-June).
Plant maximal spreads were measured along the longest axis between two opposite margins of the plant. While plant minimal spread was taken along a straight line between two plant margins that were perpendicular to the maximal Plant_S. Two directional Plant_S (plant maximal and minimal spread) were collected and then averaged to represent Plant_S. Plant_H was measured from the bed surface to the highest point of the plant. The value of Plant_Z was determined by Plant_H, plant maximal spread, and plant minimal spread and calculated using the formula: Plant_Z (m 3 ) = [π × (plant maximal spread ÷ 2) × (plant minimal spread ÷ 2) × Plant_H].

Statistical analysis
The statistical software JMP Pro 10.0.2 (SAS institute Inc., Cary, NC, USA) was used to calculate progeny distribution for each trait studied and Pearson's correlation coefficients and to estimate the broad-sense heritability for each trait. All board-sense heritability (H 2 ) estimates were calculated using the following statistical model: y ijk = µ + G i + E j + G i × E j + B h(j) + ε ijk , where y ijk represents the measured phenotypic value of the studied trait for individual plant ijk , µ the population mean value for the specific trait, G i the genetic effect, E j the environment effect, G i × E j the effect of interactions between genotype and environment, B h(j) the block effect, and ε ijk the random error. All components (G i , E j , G i × E j , B h(j) , and ε ijk ) in this model were treated as random effects.

Calculation of heterozygosity in RILs and at marker loci in LGs
The heterozygosity level of each RIL and at each marker locus was calculated using the marker genotyping data described by Guo et al. 25 . The genotyping data consist of data from 482 SNP and 36 SSR markers. Molecular marker genotypes were categorized into either being heterozygous or homozygous. The level of heterozygosity (%) in each RIL was calculated by dividing the total number of heterozygous marker loci in each RIL by the total number of marker loci analysed. The level of heterozygosity at each marker locus was obtained by dividing the total number of heterozygous RILs by the total number of RILs analysed. The resulting data were plotted in the software Matplotlib 43 to show a genome-wide landscape of heterozygosity retention with and among LGs.

QTL identification and analysis
The genetic linkage map described by Guo et al. 25 was used for QTL identification and localization in this study. The genetic map contained 518 bins (482 SNPs and 36 SSRs) spanning a total genetic distance of 220.2 cM across petunia's seven chromosomes. Molecular markers in this genetic map could be located to 620 scaffolds, 0.74% of the total number of scaffolds in the assembled P. axillaris genome (https://solgenomics. net/organism/Petunia_axillaris/genome) 42 . Nevertheless, these 620 scaffolds contain 747,650 kb of nucleotides, which is approximately 53.4% of the P. axillaris genome (1.4 Gb). The software MapQTL 6.0 44 was employed for QTL analysis. Putative QTL regions were first determined by interval mapping and the resulting highest scored markers were then highlighted and labelled. These stamped markers were subsequently treated as cofactors and run in MQM. The LOD thresholds for putative QTL were determined by permutation tests (1000 times per run) with the significant threshold at 95th percentile of LOD scores. Only QTL with a LOD score more than the LOD threshold value were declared and retained in the analysis and reported here. To verify the QTL detected in MapQTL and to estimate the QEI and QTL × QTL (epistasis) effects, the ICIM model software IciMapping 4.1 45 was used; the LOD threshold value 3.0 was used to declare significant QTL and QEI, and the LOD cutoff of 5.0 was used to declare the presence of epistasis.