Mapping and identification of genetic loci affecting earliness of bolting and flowering in lettuce

Photoperiod and temperature conditions elicit different genetic regulation over lettuce bolting and flowering. This study identifies environment-specific QTLs and putative genes and provides information for genetic marker assay. Bolting, defined as stem elongation, marks the plant life cycle transition from vegetative to reproductive stage. Lettuce is grown for its leaf rosettes, and premature bolting may reduce crop quality resulting in economic losses. The transition to reproductive stage is a complex process that involves many genetic and environmental factors. In this study, the effects of photoperiod and ambient temperature on bolting and flowering regulation were studied by utilizing a lettuce mapping population to identify quantitative trait loci (QTL) and by gene expression analyses of genotypes with contrasting phenotypes. A recombinant inbred line (RIL) population, derived from a cross between PI 251246 (early bolting) and cv. Salinas (late bolting), was grown in four combinations of short (8 h) and long (16 h) days and low (20 °C) and high (35 °C) temperature. QTL models revealed both genetic (G) and environmental (E) effects, and GxE interactions. A major QTL for bolting and flowering time was found on chromosome 7 (qFLT7.2), and two candidate genes were identified by fine mapping, homology, and gene expression studies. In short days and high temperature conditions, qFLT7.2 had no effect on plant development, while several small-effect loci on chromosomes 2, 3, 6, 8, and 9 were associated with bolting and flowering. Of these, the QTL on chromosome 2, qBFr2.1, co-located with the Flowering Locus T (LsFT) gene. Polymorphisms between parent genotypes in the promotor region may explain identified gene expression differences and were used to design a genetic marker which may be used to identify the late bolting trait.


Introduction
The phase transition from vegetative to the reproductive stage is a critical step in the life cycle of flowering plants that affects the successes of flowering, fruit set, and seed production. In lettuce (Lactuca sativa L.), the transition from the vegetative to reproductive stage occurs in the meristem concealed within the rosette. During the transition, the vegetative shoot apical meristem cap is elongated, and then, the microscopic floral primordia are formed (Chen et al. 2018a, b). Subsequently, rapid elongation of the stem internodes occurs, a process called bolting, and the inflorescence expands followed by flowers opening over the course of several days. Since lettuce is a crop grown for its leaf rosettes, early bolting is an undesirable trait. Determining the genetic and environmental factors affecting bolting would increase the ability to develop cultivars with greater resistance to premature bolting. Communicated by Herman J. van Eck.
The vegetative-flowering transition is regulated by a complex network of genetic and environmental factors. In the well-studied Arabidopsis (Arabidopsis thaliana (L.) Heynh) model, more than 100 genes have been implicated in the control of flowering time (Komeda 2004;Putterill et al. 2004;Amasino 2010;Bouché et al. 2016). Internal regulatory pathways converge with pathways integrating environmental signals to control the timing of meristem differentiation. Endogenous factors influence floral induction and include the gibberellin (GA) pathway, circadian clock, and plant age (Fornara et al. 2010;Song et al. 2015;Bao et al. 2020). Environmental cues influencing flowering time include day length, light quality, ambient temperature, and vernalization (Michaels and Amasino 2000;Samach and Coupland 2000;Blázquez et al. 2003;Song et al. 2013;Romera-Branchat et al. 2014;Capovilla et al. 2015;Cho et al. 2017).
Lettuce is a facultative long-day plant, and cultivars and accessions exhibit wide variability in the rate of transition to flowering (Ryder 1988;Lafta et al. 2017). Very early cultivars (Ryder 1996) start to bolt and flower seven weeks after planting in long days in a greenhouse while very late cultivars may flower as late as four months after planting. Besides genetic factors, environmental factors substantially affect the rate of transition to the reproductive phase (Ito et al. 1963;Jenni et al. 2013). Earlier bolting (less than 90 days) can be induced in certain genotypes by increased day length (Waycott 1995) or by high ambient temperatures (Rappaport et al. 1956), thus leading to a strong genotype x environment (GxE) interaction (Lafta et al. 2017;Wien 2020).
For lettuce production, a prolonged vegetative phase is beneficial for accumulation of leaf biomass and to maintain quality. Premature bolting results in elongated cores and/or cracked heads that are not suitable for the market. In addition, the initiation of flower signaling causes biochemical changes to form accumulation of latex in the leaves which makes them undesirably bitter (Simonne et al. 2002). Premature bolting can lead to partial or complete yield loss, particularly in early fall planting periods in hot production regions. Therefore, the genetic basis for regulation of bolting and flowering has been the focus of numerous studies over the years (Thompson 1943;Rappaport et al. 1956;Ryder 1983;Ryder and Milligan 2005;Han et al. 2021). Ryder (1983Ryder ( , 1988 identified two genes controlling early flowering (Ef) Ef-1 and Ef-2, with earliness being incompletely dominant. The author later described four more Ef genes in additional crosses (Ryder and Milligan 2005). More recently, multiple quantitative trait loci (QTLs) for bolting and flowering time, as well as for stem length, which correlates with earliness of bolting, have been identified on lettuce chromosomes 2 and 7 (Jeuken and Lindhout 2004;Hartman et al. 2012Hartman et al. , 2013Han et al. 2021) in crosses between cultivated lettuce and its wild relatives. In a recombinant inbred line (RIL) population developed from crossing late and early bolting lettuce cultivars and evaluated in four summer experiments, three QTLs were identified on chromosomes 2, two QTLs on chromosome 5, and one each on chromosomes 7 and 8 (Jenni et al. 2013). Linkage mapping performed on another RIL population developed from crossing two cultivars with similar time to bolting revealed QTLs on chromosomes 1, 4, and 5 (Mamo et al. 2019). Perhaps due to the polygenic nature of the trait, transgressive segregation was observed in several bi-parental populations (Silva et al. 1999;Jenni et al. 2013;Mamo et al. 2019). A genome-wide association study (GWAS) identified several single nucleotide polymorphism (SNP) markers linked to bolting (on chromosomes 1, two on 7, and one that was not mapped) and flowering (on chromosomes 5, 7, 8 and one that was not mapped; Kwon et al. 2013). Although bolting and flowering are closely related processes, mapping of each phenotype separately resulted in only partial overlap of the associated SNPs. Another GWAS study performed on ~500 accessions identified a major QTL for developmental rate on chromosome 7, while minor QTLs with substantial QTL x environment interaction were detected on most of the other chromosomes (Sthapit Kandel et al. 2020).
In lettuce, the homolog of the Arabidopsis FLOWERING LOCUS T (FT) gene had only one copy and LsFT expression levels increased with higher ambient temperature and induced earlier differentiation to flowering (Fukuda et al. 2011). The response of LsFT to heat treatment was also characterized in FT RNAi plants (Chen et al. 2018a). In these lines, there was also a reduction in expression of putative downstream genes (based on homology to Arabidopsis AP2, AP3, and LFY). Another study showed association between the bolting time and the expression levels of additional Arabidopsis homologs in lettuce: LsFVEL, LsFLDL, LsLDL, LsLFYL, LsAP1L (Fukuda et al. 2017). The lettuce homolog of AtSOC1 (LsSOC1), an important floral integrating gene, was found to be expressed in lettuce shoot apical meristems during the reproductive phase. Silencing of this gene greatly delayed flowering (Chen et al. 2018a, b). LsSOC1 was also shown to greatly increase upon heat treatment which promoted flowering.  found many genes that were differentially expressed between two lettuce cultivars with different bolting time. The authors implicated a group of MADS-box and GA genes as flowering regulators (Ning et al. 2019). GA and related genes were shown to be related to bolting (Fukuda et al. 2009;Umetsu et al. 2011) specifically in response to high temperatures (Liu et al. 2018). The similarities in expression and function of these central genes suggest the main structure of the signaling pathway to be conserved between lettuce and Arabidopsis. There are many examples where homologs of Arabidopsis genes function differently in crop species (Blümel et al. 2015) and 1 3 paralogs that broadened natural variation in flowering time (Blackman et al. 2010). In other cases, several genes have been found in crops that have no Arabidopsis homologs. Therefore, focusing solely on Arabidopsis homologues could lead to identifying new players or overlook the central genes that are involved in regulation of bolting and/or flowering in lettuce.
While studying functions of gene homologs and molecular genetic manipulation sheds light on the molecular mechanism of flowering in lettuce, their use for traditional breeding is limited. For conventional breeding that uses crossing and selection, information about beneficial alleles that exist in plant accessions, including marker-assisted selection (MAS) (GuoYou et al. 2010;Simko 2013), is most needed.
To identify novel genes for lettuce flowering time, we used a recombinant inbred lines (RIL) population segregating for flowering time. The RIL population was created from a cross between PI 251246 and cv. Salinas (Yoong et al. 2016). Parental line PI 251246 is a primitive type lettuce accession cultivated for seed oil that skips the rosette stage and flowers within six weeks after germination (Grube and Ryder 2004). Parental line cv. Salinas, an iceberg type, is relatively late to bolting, 95 days under long-day conditions and ~ 132 days under short-day conditions (Waycott 1995).
The wide range of flowering phenotypes and environmental sensitivities observed in the RIL population, and careful evaluation in controlled environments, revealed multiple environment-specific QTLs and interactions among them. Fine mapping in conjunction with gene sequencing and expression studies allowed identification of candidate genes and the development of a MAS assay that can be used by breeders to confer bolting resistance.

RIL population, plant growing conditions, and phenotyping of phenological stages
A F 7 RIL population, consisting of 161 families, was evaluated under controlled environmental conditions. The RIL was derived from single seed descent following a cross between PI 251246 and cv. Salinas. For the first 14 days, plants were grown in controlled environment chambers under a photoperiod regime consisting of 8 h light and 16 h dark at 20 °C for 14 days to allow seedling establishment. Light was provided by a combination of fluorescent and incandescent electric bulbs. After 14 days, plants were transferred to growth chambers under four light and temperature conditions: long day and low temperature (LD-LT; 16 h light /8 h dark at 20 °C), long days and high temperature (LD-HT; 16 h light /8 h dark at 35 °C), short days and low temperature (SD-LT; 8 h light / 16 h dark at 20 °C), and short days and high temperature (SD-HT; 8 h light/ 16 h dark at 35 °C). In each of these four environmental conditions, the phenological stage of three plants per RIL and both parents was evaluated weekly using a rating scale of 1 = rosette; 2 = bolting-visible internode elongation; 3 = visual buds; 4 = expanded inflorescence; 5 = flowering-opening of first flower; 6 = more than half of buds flowered; 7 = open involucres. Also evaluated was the number of days to bolting (D BLT ; days after germination to reaching stage 2) and days to flowering (D FLT ; days after germination to reaching stage 5). The population-wide response to environmental conditions was statistically evaluated with two-tailed Student's t-test using the population mean and the variance that were calculated from the values of all RILs grown at each environment. The experiment-wise p-value lower than 0.05 was considered statistically significant; thus, the test-wise error rate was adjusted for multiple comparisons using Bonferroni correction. It has been shown that D BLT and D FLT share both common and unique pathways (Kwon et al. 2013). To potentially identify regulation of flowering distinct from bolting, the residuals from the linear regression of D FLT to D BLT (D BFr ) were calculated for each RIL and used alongside other traits as an independent phenotype for QTL analysis.

Genotyping and linkage map construction
DNA was extracted from each RIL family (DNeasy Plant Mini Kit,Cat No. 69104,Qiagen Inc.,Valencia,CA), digested with the EcoT22I enzyme, and subjected to genotyping by sequencing (GBS; Elshire et al. 2011) by the Cornell Institute of Biotechnology (http:// www. biote ch. corne ll. edu/ brc/ genom ics-facil ity). SNP calling and filtering were performed using the TASSEL SNP discovery pipeline (Glaubitz et al. 2014). A linkage map based on informative markers was generated using an online version of MSTmap (http:// mstmap. org/). Fully linked markers, with no recombination between them, were collapsed to a single combined marker. The linkage map for QTL analysis was re-generated to verify the correct placing and centimorgan (cM) distance determination of the combined markers and the original ones.

QTL mapping and model selection
QTL mapping and model selections were performed in R using the R/qtl package (Broman et al. 2003). Simple interval QTL mapping (SIM) was performed using the Haley-Knott regression method (Feenstra et al. 2006). The significance LOD threshold for each environmental condition was determined by 1,000 permutations. QTLs were deemed significant when exceeding α = 0.05 LOD threshold and suggestive when exceeding α = 0.1 LOD threshold. QTL naming followed the nomenclature convention using q + trait 1 3 abbreviation + chromosome number + QTL number on the particular chromosome. Trait abbreviations followed those suggested by Han et al. (2021).
All significant and suggestive QTLs were considered in the model fitting for D BLT and D FLT in each environmental condition separately. Initially, additive and interactive effects of all QTLs were considered in the model. Subsequently, a backward selection was applied to eliminate QTLs out of the model when they had no significant additive effect or interaction. The Bayesian information criterion (BIC; Broman and Speed 2002) with a δ value of 2 to minimize false positives was calculated for each evaluated model. Finally, the best model was selected based on the lowest BIC, i.e., the model explaining the highest variance with the least number of parameters.

Fine mapping
To enable fine mapping of the qBLT7.2-qFLT7.2 QTL region, an F 2 population was developed from a cross between two RILs (RIL 68 and RIL 69) with high similarity of marker alleles across the whole genome, but harboring contrasting marker alleles in the region of interest.
Genomic polymorphisms between RIL 68 and RIL 69 in the region of interest were identified by SNP mining of RNA-seq data generously provided by Fei-Yian Yoong (Yoong et al. 2016). To identify additional polymorphisms, Sanger sequencing was performed (Eton Bioscience Inc., San Diego, California) on PCR amplification products that targeted genes and inter-gene segments within the region of interest. Polymorphisms identified along the qBLT7.2-qFLT7.2 region were used to design high-resolution DNA melting (HRM) markers (Simko 2016; Supplementary  Table 1). Testing for recombination and genotyping of all recombinants with the HRM markers was performed with the LightScanner instrument (BioFire Diagnostics LLC, Salt Lake City, Utah). PCR conditions and florescent detection were according to manufacturer's instructions. Specifically, PCR annealing temperature was 66 °C, and elongation temperature was held at 72 °C for 30 s and LightScanner temperature scan ranged between 72 and 96 °C.
To identify recombinants, 2,094 F 2 plants, all derived from a single F 1 plant from the RIL 68 × RIL 69 cross, were evaluated. Of these, 773 plants were found to have recombination between two co-dominant markers flanking the QTL region. Recombinant plants were genotyped with 18 additional markers, designed throughout the qBLT7.2 / qFLT7.2 region. The phenological development of recombinants was evaluated weekly on plants grown under long days and mild temperature (18-25 °C) conditions in which the parental RILs (RIL 68 and RIL 69) significantly differed in their flowering time. Differences between parent genotypes and F 2 plants used in the fine mapping study were statistically evaluated with two-tailed Student's t-test using the mean and the variance calculated from the values of at least five plants. The experiment-wise p-value lower than 0.05 was considered statistically significant; thus, the test-wise error rate was adjusted for multiple comparisons using Bonferroni correction.

qBFr2.1 allele effect validation
In order to validate the effect of the P allele (originating from PI 251246 parent) on flowering time, F 3 plants from the RIL 68 × RIL 69 cross were evaluated. Plants selected by HRM markers were homozygous for P alleles in qBFr2.1 and for S alleles (originating from cv. Salinas parent) in other detected D FLT QTLs. For further validation, a series of back crosses of the F 3 plants to cv. Salinas was conducted, selecting for the P allele with the HRM marker LZ001. F 3 and backcross plants were evaluated for D FLT in the greenhouse during spring under ~ LD-LT conditions.

Gene expression assays
To examine the range of gene expression responses of flowering related genes, the RIL parental lines (PI 251246 and cv. Salinas, along with the RILs used for fine mapping, RIL 68 and RIL 69, were used. The plants were grown under the four photoperiod and temperature conditions (LD-LT, LD-HT, SD-LT, SD-HT) and evaluated weekly using the phenological rating scale. In previous experiments a ~ 14-day difference between appearance of a first bud and a first flower was consistently observed. Furthermore, the appearance of the bud demonstrated that the transition from a vegetative to a reproductive phase occurred. Therefore, in this experiment D BUD (days after germination to reaching stage 3) was used to quantify developmental progress. The response to environmental conditions of the genotypes used in the gene expression study was statistically evaluated with two-tailed Student's t-test using the mean and the variance calculated from the values of at least nine plants. The experiment-wise p-value lower than 0.05 was considered statistically significant; thus, the test-wise error rate was adjusted for multiple comparisons using Bonferroni correction.
At selected time points, leaf and stem tip samples were harvested for RNA extraction. Leaf samples were taken from the distal part of the second-to-newest fully expanded leaf in the rosette. Stem tip samples included the apical meristem with minimal wrapping of young leaves. To improve uniformity of results, all samples were collected between hour two and three of the growth chamber light period. For each sample, three independent pools of tissue were collected from three individual plants and flash-frozen in liquid nitrogen. Total RNA was extracted from each sample using the RNeasy Plant Mini Kit (Cat No. 74904, Qiagen Inc., Valencia, CA)) according to the manufacturer's instructions.
A total of 31 candidate genes within the qBLT7.2 / qFLT7.2 QTL region were identified using the lettuce reference genome (ver. 4; Reyes-Chin-Wo et al. 2017). Genes whose predicted function was deemed unlikely to be involved in reproductive transition regulation were excluded from further consideration, resulting in 12 candidate genes retained. To better understand the effects of the QTL alleles on the reproductive transition regulatory network, 15 genes known to be involved in flowering time regulation in Arabidopsis were tested for expression (Supplementary Table 2). Of these, eleven have been identified in lettuce Huo et al. 2016;Ning et al. 2019) including MADS-box genes. The MADS-box genes assessed in this work are arbitrarily named MAD1-7 because of their high sequence homology. Additionally, three normalization genes (PP23, PP21, and TIP4) were chosen from the literature based on their low constitutive expression in various tissues and in different experimental conditions. PP23, PP21, and TIP4 encode for Protein Phosphatase 2a subunit 3, Protein phosphatase 2A-1, and TIP41-like protein, respectively (Czechowski et al. 2005;Borowski et al. 2014;Sgamma et al. 2016).
The GenomeLab GeXP Analysis System (SCIEX, USA) was used to assay mRNA levels in tissues with multiplex primer design (Supplementary Table 2). Gene expression primers for multiplexed reactions were designed using GenomeLab eXpress Profiler software. The multiplexed expression panel consisted of 27 candidate genes plus the three normalization genes. Complementary primer binding sites for both genotypes were verified by sequencing. The cDNA for GeXP was synthesized from 50 ng of total RNA using the GenomeLab GeXP Start Kit. PCR and multiplex detection were performed according to manufacturer's instructions (Hayashi et al. 2007). Gene expression was quantified as area under the peak. All peaks under a background noise threshold of 500 counts were removed, and remaining peaks were normalized to the intensity of the internal control Kan r gene. Since TIP4 was not stably expressed in the experiment, in order to select optimal genes for normalization, the full panel of 30 genes was subject to an expression stability test following the GeNorm algorithm (Vandesompele et al. 2002). The geometric mean of expressions levels of the selected three most stable genes, PP21, PP23, and 7b_8, was used for normalization of expression levels of all target genes.
Lactuca sativa FLOWERING LOCUS T gene expression was independently evaluated by RT-real-time PCR, using the same RNA samples as were used for the GeXP analysis. Reverse transcription was performed using 170 ng of total RNA with the High-Capacity RNA-to-cDNA™ Kit (Cat no: 4387406, Thermo Fisher Scientific Inc., Foster City, California). The normalization gene PP23 was selected for real-time PCR based on having lowest variation in the GeXP panel. Target and normalization genes were detected in the same well by a duplex TaqMan reaction using reporter dyes FAM and VIC, respectively, and custom-designed primers and probes (Supplementary  Table 2), with an annealing temperature of 60 °C (Custom TaqMan® Gene Expression Assays, Thermo Fisher Scientific Inc., Foster City, California). Color compensation to account for differences in using different dyes was applied to the raw data. Dose curve experiments revealed amplification efficiency for both genes and genotypes ranging from 1.98 to 2.02; thus, an efficiency of 2 was used for calculation of target to reference concentration ratio (conc. ratio = 2 −(Cp target−Cp ref) ). Each sample was tested in triplicate, and their average was used for further statistical analysis.
Gene expression differences between the genotypes and controlled environment condition were compared using Tukey's Honest Significant Difference (HSD) test. For comparison of environmental condition effect on gene expression, values of the four environmental conditions at a specific time point were compared by Tukey's HSD test for each genotype. For comparisons of genotype effect on gene expression, values of the four genotypes at a specific time point were compared by Tukey's HSD for each condition. An indication of a significant increase or decrease was made by Student's t-test, comparing gene expression values of the specific genotype and environmental condition between the time points described in the text. Significance was determined at α = 0.05 corrected for the number of comparisons in each time point.

Candidate gene homology, sequencing, and polymorphism annotation
Candidate genes were identified by positioning GBS markers to the lettuce reference genome (Reyes-Chin-Wo et al. 2017) and using blastp to search the Arabidopsis Information Resource (TAIR) protein database with the predicted lettuce protein sequence. Reciprocally, the protein of the best Arabidopsis match was used to search the lettuce TSA database in NCBI using tblastn (McGinnis and Madden 2004).
Candidate genes and their promotor regions were sequenced using PCR amplification of overlapping gene segments. The PCR products were purified and subjected to Sanger sequencing (Eton Bioscience Inc., San Diego, California). To better understand the biological significance of differences discovered in promotor regions, an in silico annotation based on known Arabidopsis promotor motifs was performed using AgRIS (Davuluri et al. 2003) and PLACE (Higo et al. 1999).

QTLs associated with day length and temperature
Under all environmental conditions, D BLT were lowest for PI 251246 and highest for cv. Salinas with no transgressive segregation observed among RILs. In contrast, transgressive segregation for D FLT was observed with several RILs families displaying both earlier and later D FLT compared to PI 251246 and cv. Salinas, respectively (Fig. 1).
Both D BLT and D FLT of RILs grown under HT were significantly shorter (p < 0.0001, t-test) than under LT (Table 1). Day length did not significantly affect D BLT , but D FLT was significantly shorter (p < 0.0001, t-test) in LD compared to SD. The range of days between bolting and flowering also differed among the four environmental conditions, and both HT and LD significantly (p < 0.0001, t-test) reduced the observed range. These results support the idea that D BLT and D FLT share both common and unique pathways (Kwon et al. 2013;Chen et al. 2018b) and hence support the use of residuals of the linear correlation between them, D BFr , as an additional phenotypic trait.
QTL analysis using the RIL population revealed a total of two bolting QTLs across the four growing conditions, accounting for 13.6-42.9% of the total phenotypic variance explained (PVE) (Fig. 2, Table 2). One significant flowering QTL was identified, which accounted for 24.0-51.9% of the PVE in various conditions. In addition, three QTLs were identified for D BFr , which accounted for 13.0-37.4% of the PVE.
A major QTL for D BLT , qBLT7.2, was found under LD-LT conditions while the same locus showed only a minor effect at LD-HT and SD-LT and no observable effect at SD-HT conditions. This locus also had a large effect on D FLT , qFLT7.2, in LD conditions under both high and low temperature (PVE = 52.0% and 47.3% under LD-LT and LD-HT, respectively), but had a lower effect under SD-LT conditions (PVE = 24.1%). The QTL effect on D BFr was relatively small under LD-LT and SD-LT conditions (13% and 9% of PVE, respectively), but high under LD-HT (PVE = 37.5%). These results show that under LD conditions alleles at this locus affect bolting and flowering differently depending on the ambient temperature.
Three other QTLs associated with D BLT , D FLT and D BFr were identified, one on chromosome 8 (qBLT8.1), one each on chromosome 2 (qBFr2.1), and chromosome 6 (qBFr6.1). The allele contributing to the delayed flowering phenotype at qBFr2.1 originated from the early flowering parent PI 251246.

QTL models
To identify QTLs that explained the highest amount of phenotypic variation for D BLT and D FLT , Bayesian information criterion were used (Broman and Speed 2002). The PVE explained by the models ranged from 13.5% to 83.8%, of a single QTL model for bolting under SD-HT and a five QTL model for flowering under LD-LT, respectively (Table 3). An identical five QTL model was the best fit for both flowering under LD-HT and LD-LT; all other conditions tested had unique models (Table 3, Supplementary Table 3).
The models show that QTLs which were detected only at a suggestive (0.1 ≥ α ≥ 0.05) level using SIM can have a significant effect on flowering through interaction with other QTLs (i.e., epistasis). Furthermore, the extent of their effect and interaction is affected by the light and temperature conditions. The most obvious condition-dependent interaction was seen between qFLT2.1 and qFLT7.2 (Fig. 3), where at LD-HT and SD-LT alleles in these positions have constant effects, while at SD-HT only qFLT2.1 alleles affect D FLT . At LD-LT, the qFLT2.1 alleles have an effect on flowering only when the homozygous combination of alleles originating from cv. Salinas is present at the qFLT7.2 locus, thus demonstrating epistatic effect. In this population, under all environmental conditions, the alleles maintained directionality of their effect on bolting and flowering.

Fine mapping of qBLT7.2 / qFLT7.2 region
Under field and controlled environment studies, qBLT7.2 / qFLT7.2 was consistently observed to have the largest PVE of all QTLs detected. Therefore, this locus was selected for   Table 1), five closely linked markers, LZ431, LZ481, LZ261, LZ299, and LZ425 displayed complete association with the flowering phenotype (Fig. 4). The phenological stage of F 2 plants heterozygous in these markers was significantly higher (p < 1e −7 , t-test) than that of plants with homozygous S allele in all closely linked markers. The QTL region of interest was re-defined based on these markers and included a 1.15-Mbp-long region (chromosome 7: 207,702,411-209,071,520 bp) which contained 31 predicted genes. Genes whose function, predicted by Arabidopsis homology annotations (Reyes-Chin-Wo et al. 2017), was unlikely to be involved in development regulation (such as glycolysis or chloroplast targeted genes) were not pursued for further analysis. The predicted annotation of twelve genes further studied in the fined mapped region was verified by Arabidopsis protein homology (Supplementary Table 2), and three genes were found homologous to genes known to be involved in flowering time regulation. CONSTANS LIKE-9 (COL9) was found to be the best matching homolog for the Lsat_1_v5_gn_7_97240.1 lettuce gene. In Arabidopsis, COL9 delays flowering time and is a FT, SOC1, and CO repressor (Cheng and Wang 2005). The full coding region and part of the UTRs of COL9 in cv. Salinas and PI 251246 were sequenced using both genomic and cDNA. A single polymorphism between sequences of the two genotypes was found in the fourth intron, but none in the coding region. The COL9 gene-specific maker used for fine mapping, LZ431, mapped to the center of the QTL region and matched exactly with the flowering time phenotype (Fig. 4d). PHYTOCHROME C (PhyC) was found to be the best match for Lsat_1_v5_gn_7_96941.1. PhyC encodes a phytochrome protein that is known to be involved in flowering time regulation in Arabidopsis (Pearce et al. 2016).
Sequencing of PhyC revealed six SNP polymorphisms in the coding sequence between cv. Salinas and PI 251246, all in the first exon. A gene-specific marker, LZ425, based on a SNP in the first exon of this gene, matched exactly with the flowering time phenotype (Fig. 4d). Of the six described SNPs, five were silent mutations, but one in position 118 (cv. Salinas, reference genome) causes a change of isoleucine to phenylalanine in PI 251246. This mutation is located within the histidine kinase A domain (SM00388), as annotated by SMART database (Letunic and Bork 2018), and serves as dimerization and phosphoacceptor domain of histidine kinases (Vierstra and Davis 2000). The isoleucine residue, corresponding to isoleucine in position 901 in Arabidopsis, is conserved within the PhyC genes in angiosperms (supplementary information). It is also conserved in histidine kinase A domains of more diverse organisms, including prokaryotes (it is either isoleucine, leucine, or alanine in over 60% of matched sequences in the SMART database). While the functional effect of this amino acid change is not known, it could alter the gene function in PI 251246.
AT4G14540.1, which is a nuclear factor Y, subunit B3, transcription regulator involved in photoperiod regulation (Kumimoto et al. 2008), was the best matching homolog for Lsat_1_v5_gn_7_96781.1. However, a gene-specific marker was not fully linked to the flowering phenotype so this gene was not considered to be a likely candidate.

Expression of candidate genes from qBLT7.2 / qFLT7.2 region
To identify genes whose gene expression pattern corresponds to QTL effects in various conditions, PI 251246, cv. Salinas, RIL 68, and RIL 69 were sampled at multiple time points during development under the same environmental conditions used for QTL mapping. Of the four genotypes, PI 251246 was the earliest to bolt in all environmental conditions, bolting 21 days after germination. Soon after bolting, buds were visible on these plants (35 D BUD in SD-LT and 28 D BUD in the other three conditions). The three other genotypes bolted and flowered later than PI 251246 (Fig. 5). In these genotypes, D BLT was significantly lower in HT conditions compared to LT (p < 0.0001, t-test), while day length did not have a consistent effect on D BLT . What seems to be an exception in cv. Salinas is probably caused by low light levels due to aged bulbs at the beginning of the growth period in the SD-LT growth chamber. Only the cv. Salinas plants in the chamber were etiolated, and therefore, accurate Instead, the presented data are corresponding to the first observation of stem internode elongation that is closely associated with bolting. While PI 251246 D BUD was early and relatively unchanged in the four conditions, RIL 69 was also early, but mean D BUD ranged 28-44 across the four environmental conditions. RIL 68 and cv. Salinas had significantly (p < 0.001, t-test) higher D BUD than PI 25146 and RIL 69, and both displayed a wide range of mean D BUD in the four environmental conditions (71-109 and 67-124 for RIL 68 and cv. Salinas, respectively). For all genotypes, both HT and LD led to lower D BUD compared to LT (p < 0.0001, t-test) and SD (p < 0.05, t-test), respectively. These data indicate that tested temperature conditions had a larger effect than day length on development rate.
Gene expression studies focused on identifying genes which displayed expression patterns consistent with D BLT and D BUD phenotypes in response to specific environmental conditions, or genes whose expression pattern differed between PI 251246 and RIL 69, with P allele (originating from PI 251246), and cv. Salinas and RIL 68, with S allele (originating from cv. Salinas). PI 251246 and RIL 69 reached budding stage in only a few weeks and therefore had only the earlier sampling points available for gene expression assays.
Of the twelve targeted genes located in the qFLT7.2 locus (Supplementary Table 2), most did not display any consistent significant patterns between environmental conditions or genotypes ( Supplementary Fig. 1). Interestingly, many genes in the study displayed an increase in gene expression in tips of RIL 68 between 42 and 49 days after germination under HT conditions region (7b_1, 7b_4, 7b_5, 7b_6, 7b_7, 7b_8, 7b_10, and 7b_11, putatively annotated as U-box protein, COL9, unknown protein1, P24, unknown protein2, unknown kinase, PhyC, and FAR1, respectively), which is not seen under LT conditions. Additionally, several genes (7b_1, 7b_3, 7b_4, 7b_5, 7b_7, and 7b_11, putatively annotated as U-box protein, HNRNPA1_3, COL9, unknown protein1, unknown protein2, and FAR1, respectively) also displayed a decrease in expression in RIL 69 tips under SD-LT during the 42 days of sampling. 7b_1 (U-box protein) had significantly lower expression under LD-HT early in the experiment (day 21 or 28) in cv. Salinas, PI 251246, and RIL 68 and significantly higher expression under SD-LT later in the experiment.
Notably, expression of COL9 (Lsat_1_v5_gn_7_97240.1) in cv. Salinas and RIL 68 leaves was significantly higher in HT conditions compared to LT conditions, in which D BUD was higher. On the other hand, in PI 251246 and RIL 69 expression of this gene was more effected by day length, as plants grown under LD displayed significantly higher expression than SD grown plants. In addition, at 14 days after germination leaf expression was higher in PI 251246 and RIL 69 compared to cv. Salinas and RIL 68 (Fig. 6). Additionally, PhyC (Lsat_1_v5_gn_7_96941.1) gene expression was significantly higher in cv. Salinas and RIL 68 compared to PI 251246 and RIL 69 in tips under all conditions (Fig. 7). In leaves, expression of the PhyC gene significantly increased in PI 251246 and RIL 69 at 21 days after germination in SD-LT, in which D BUD is higher, compared to the other environmental conditions.
Based on gene expression, homology, and sequencing, the strongest candidate genes within the fine mapped qBLT7.2 / qFLT7.2 region include PhyC and COL9. The lettuce PhyC homolog found in the QTL region has a mutation causing a change of a conserved isoleucine to phenylalanine in a domain that may influence the protein function. Moreover, the correspondence between

QTL qDBFr2.1
Under SD-HT conditions, the QTL associated with the residuals from the linear regression of D FLT to D BLT , designated as qDBFr2.1, mapped to chromosome 2 and accounted for 14.4% of the phenotypic variation observed. In this QTL, the late flowering allele was contributed by the PI 251246 parental line (Table 2). When considering multi-QTL models, qDBtFr2.1 was included in two of the four environmental conditions for D BLT and all four conditions for D FLT (Table 3). Within this genomic region, only a single gene, Lsat_1_v5_gn_2_17881, was found to have an Arabidopsis homolog known to be related to flowering, the lettuce Flowering Locus T (LsFT) gene. FT is a widely conserved gene in dicots and monocots and is a key integrator of floral induction signaling pathways (reviewed in: Ballerini and Kramer 2011). To validate the location of this gene in the QTL, a gene-specific HRM marker (LZ001) mapped to the center of the qDBFr2.1 QTL.

LsFT gene expression
At 14 days after germination, when all plants were grown under SD-LT, gene expression levels of LsFT were assessed on leaves and tips of cv. Salinas, PI 251246, RIL 68, and RIL 69. PI 251246 had the highest LsFT expression in leaves, with lower expression in RIL 69, and very little to no expression was observed in RIL 68 and cv. Salinas (Fig. 8). In tips, LsFT expression was highest in cv. Salinas, but only lowly expressed in PI 251246, RIL 69 and RIL 68 (Fig. 8).
In leaves of PI 251246 sampled 21 days after germination, expression of LsFT increased about fourfold higher than observed at 14 days after germination under LD-HT (Fig. 9). LsFT expression of leaves in RIL 69 increased by four to five fold at 21 and 28 days after germination under LD-HT, LD-LT, respectively (Fig. 9). In cv. Salinas, LsFT expression increased about tenfold above that observed at 14-days under LD-HT (Fig. 9). LsFT expression remained essentially unchanged in RIL 68 up through 80 days after germination (Fig. 9).
For tips, LsFT increased sharply in PI 251246 by 21 days under all four environmental conditions (Fig. 9). In RIL 69, LsFT increased under LD-LT and SD-LT (but no data for LD-HT, SD-HT) at 28 and 42 days after germination, respectively (Fig. 9). In RIL 68 and cv. Salinas, LsFT expression fluctuated under all environmental conditions, but remained essentially unchanged across all sampling dates (Fig. 9).

LsFT sequencing
Sequencing of genomic and cDNA revealed no polymorphisms between cv. Salinas and PI 251246 in the coding sequence of the gene. There were, however, numerous SNPs and indels in introns and the promotor region. The region up to 600 bp upstream of the transcription start site had relatively few polymorphisms. However, the region between 660 and 1,186 bp upstream of the start codon had many SNPs and indels in PI 251246 compared to cv. Salinas. Based on in silico motif annotation (Davuluri et al. 2003) and (Higo et al. 1999), two GATA motifs (agatag or agataa) and a SORLIP1 motif (agccac) were found to in cv. Salinas but were absent in PI 251246. In Arabidopsis, these motifs are related to light response (Terzaghi and Cashmore 1995;Adrian et al. 2010). None of the predicted motifs were found only in PI 251246 while absent in cv. Salinas (Fig. 10).

QTL qBFr2.1 effect validation
To validate the effect of the delayed flowering phenotype contributed by the P allele from early bolting parent PI 251246, we phenotypically evaluated F 3 plants from the RIL 68 × RIL 69 cross which were homozygous for P alleles in the qBFr2.1 region. When the plants were grown in the greenhouse under ~ LD-LT conditions, they flowered during the same week or one week later than cv. Salinas. In a series of backcrosses of F 3 plants to cv. Salinas, selection for the P allele was done using the LZ001 HRM marker. BC 2 plants flowering ranged from 92 to 128 days with an average of 111.3, while cv. Salinas plants flowered after 100 days under LD-LT (Table 4). Furthermore, in BC 1 S 1 families, which are expected to segregate, D FLT ranged 85-153. In each BC 1 S 1 family, at least two of the ten plants evaluated flowered more than two weeks later than cv. Salinas plants. These results demonstrate the late flowering effect of the LsFT (qBFr2.1) P allele, derived from the early flowering PI 251246 parent, in the background of cv. Salinas. Moreover, since the effect of this QTL is much higher in warmer temperatures (Table 2), the benefit of the allele is expected to be even greater under such conditions.
Since marker LZ001 is located within the first intron of the LsFT gene, is tightly linked to qBFr2.1, and distinguishes between S and P alleles, its utility for markerassisted selection was tested on a panel of diverse lettuce cultivars. Cultivars Amazona, Anais, Ancora, Reine des Glaces, Costa Rica, Darkland, Eruption, Flashy Troutback, Green Towers, Hearts Delight, Iceberg, La Brillante, Lolla Rossa, PI 491224, Pacific, Salad Crisp, Sniper, Tiber, Tinto, and Valmaine all had the same marker allele as cv. Salinas. A Lactuca serriola accession (11G99) had the same P allele HRM profile as PI 251246. Only cv. Sentry had a HRM pattern distinct from those for both S and P alleles, indicating that this cultivar harbors an allele different from all other tested genotypes.

Discussion
This work identified multiple QTLs affecting bolting and flowering time, which is not surprising considering the complex regulatory network controlling flowering time in the plant kingdom (reviewed in : Komeda 2004;Putterill et al. 2004;Andrés and Coupland 2012). The phenotypic divergence and genetic distance between the RIL population parents allowed detection of multiple QTLs. The current finding that most QTLs were involved in epistatic interactions strengthens the previous results documenting epistatic effect of QTLs on flowering time (Ryder and Milligan 2005;Schwartz et al. 2009).
The critical effect of day length and temperature on plant development (Koornneef et al. 1998;Song et al. 2013) and lettuce flowering (Ito et al. 1963;Jenni et al. 2013;Lafta et al. 2017) is well established. This study aimed to dissect the effect of these environmental factors and their combinations over genetic regulation of development. The QTL detected under the SD-HT condition were not detected under the other three environmental conditions tested, suggesting a different bolting and regulatory network (Table 2). While under both LD conditions and also in SD-LT one major  (Kumar et al. 2012;Fernández et al. 2016). Similarly, in fava bean (Vicia faba) under SD-HT many small QTLs were revealed compared to LD conditions (Catt et al. 2017). Despite the extremely divergent bolting and flowering times of cv. Salinas and PI 251246, transgressive segregation was observed for D FLT . In two detected QTLs and one suggestive QTL, qBFr2.1, sqBFr3.1, and sqFLT9.1, late flowering alleles were contributed by the very early flowering parent. Presence of these alleles, and their interaction with other loci, could explain the transgressive segregation. If these alleles are not yet present in domesticated lettuce cultivars, they can contribute to development of lettuce breeding lines and cultivars with late flowering. For bolting, on the other hand, no transgressive segregation was observed, and either the qBLT7.2 allele from PI 251246 determined early bolting or, under SD-HT, the PI 251246 allele in qBLT8.1. QTLs for flowering, stem length, and plant development were previously found on this chromosome (Hartman et al. 2012;Jenni et al. 2013;Kwon et al. 2013;Sthapit Kandel et al. 2020).
To date, all published research regarding QTLs of bolting and/or flowering time-related traits, whether using bi-parental population or diversity panel for GWAS analysis, identified QTLs in chromosome 7 (Jeuken and Lindhout 2004;Hartman et al. 2012;Jenni et al. 2013;Kwon et al. 2013;Lee et al. 2021;Sthapit Kandel et al. 2020;Han et al. 2021;You, Still, Rosental, Hayes and Simko-unpublished data]. Although exact comparisons of QTL intervals is not always possible, the consistent detection supports the importance of chromosome 7 for lettuce flowering control. The current findings of qFLT7.2, and its high effect QTL explaining a large percentage of the variance in bolting and flowering time, confirm the contribution of this region, specifically in LD conditions. Bolting-flowering QTLs on chromosome 2 have been mapped in other RIL populations from inter-and intraspecific crosses and qBFr2.1 consistently had a significant additive effect in high temperature (HT) conditions (Hartman et al. 2012;Jenni et al. 2013). Furthermore, an L. serriola accession used in previous bolting studies (Hartman et al. 2012) had the same P allele for marker LZ001, linked to the qBFr2.1, as PI 251246.
The two QTLs focused on in this study, qBFr2.1 and qFLT7.2, had epistatic relations, and the effect of the qFLT7.2 QTL was specific to the environmental condition. In SD-HT, when qFLT7.2 had no additive effect, no epistatic effect was observed in the models and interaction plot (Fig. 3). However, in LD condition, the QTL P allele determined early flowering while S allowed qBFr2.1 QTL alleles to affect D FLT . Since D FLT in the F 2 plants heterozygous in qFLT7.2 was early, similar to the P genotype, S can be defined as the recessive allele and could perhaps indicate loss of function. This suggests that in environmental conditions and genotypes where qFLT7.2 is functional and has an effect, it overrides any effect qBFr2.1 may have. Consequently, only in conditions or alleles where qFLT7.2 is not functional, qBFr2.1 (and perhaps other loci, as seen in SD-HT) conveys their effect. In other words, qBFr2.1 is not 'induced' only under SD-HT, but rather is manifested in these conditions due to lack of qFLT7.2 effect. This observation is supported by the gene expression pattern of LsFT, the putative gene underlying qBFr2.1 QTL. During most of development LsFT expression seems to be governed more by qFLT7.2 allele than by the qBFr2.1 allele. Specifically, RIL 68 expression in tips and leaves is similar to cv. Salinas and not like PI 251246 and RIL 69, with whom it shares P allele in qBFr2.1. The latter genotypes have higher LsFT expression which increases in the first weeks of development and have an early flowering phenotype. This suggests that qFLT7.2 is an LsFT expression inducer and supports the detailed work describing the effect of LsFT expression to induce flowering in lettuce (Fukuda et al. 2011(Fukuda et al. , 2017. In accordance with FT's known role as floral integrator (Kim et al. 2013), other genetic and environmental factors influence FT expression levels. In this study, under HT conditions, FT had higher expression, as found previously in lettuce (Fukuda et al. 2011;Han et al. 2016) and Arabidopsis (Schwartz et al. 2009;Kim et al. 2012). In Chrysanthemum (Chrysanthemum morifolium) leaves, on the other hand, heat caused reduction in FT-like gene expression and delayed flowering (Nakano et al. 2013). However, the increase under HT conditions in leaves is not translated to an increase in tips in genotypes with qFLT7.2 S allele. The lack of epistasis in the SD-HT QTL model stands in contrast with the apparent influence of qFLT7.2 on LsFT gene expression in this condition, as well. 1 3 The only time point when cv. Salinas LsFT gene expression matched its late flowering allele in qBFr2.1 locus was in the tips of 14-day old plants. At this point, all plants were grown in SD-LT conditions, considered to be least inductive of flowering. The exclusive high expression in cv. Salinas's tips, not seen in leaves, may initiate an induction of flowering in the absence of qFLT7.2. The inconsistency in LsFT expression between leaves and tips, also observed in later stages in cv. Salinas and RIL 68, is especially interesting given previous reports that FT moves from leaves to the meristem within the tip where its protein is activated (Corbesier et al. 2007;Putterill and Varkonyi-Gasic 2016).
FT has a long promoter with many cis regulatory elements which affect its expression (Schwartz et al. 2009;Adrian et al. 2010;Liu et al. 2014), such as those identified to be different between cv. Salinas and PI 251246. The lettuce FT gene is located in a relatively gene-poor region, with only one gene detected within a 3 Kbp upstream transcription start site with no others for 116 Kbp, as seen in other plant species (Adrian et al. 2010). The modified elements in PI 251246, GATA motif, and SORLIP1 motif, may influence the expression pattern differences seen between tips of cv. Salinas and the P allele genotypes at 14 days.
As opposed to the most obvious candidate gene underlying the qBFr2.1 QTL, fine mapping was required to shorten the candidate list for qFLT7.2. After narrowing the interval to a 1.15 Mbp region, two genes remained plausible candidates. Still, the possibility that both genes have a role in the control of flowering in lettuce, either together or each at different conditions, cannot be excluded. Lsat_1_ v5_gn_7_96941.1 was identified as an Arabidopsis PhyC homolog. PhyC belongs to a family of phytochrome factors (Sharrock and Quail 1989), which have important roles in multiple aspects of plant growth and development, including germination, chloroplast development, photomorphogenesis, shade avoidance, and photoperiod-dependent flowering (Kendrick and Kronenberg 1994;Franklin et al. 2003;Nishida et al. 2013). In Arabidopsis, PhyC is not required for long-day induction of flowering time (Franklin et al. 2003). It is, however, involved in day length perception and delay of flowering in SD, together with PhyA and PhyB (Monte et al. 2003). In wheat and barley, on the other hand, PhyC was found to have a key role in LD-induced flowering (Nishida et al. 2013;Pearce et al. 2016).
Phytochrome family proteins bind a billin chromophore and share a basic structure (Pham et al. 2018), with distinct and partially overlapping roles in development control and signal transduction. The N terminus of phytochromes has a kinase domain, and kinase activity was shown to be important for plant photo-response (Shin et al. 2016). Gene sequencing of both parents revealed a mutation in PI 251246 located in a conserved amino acid within the histidine kinase A domain. While the effect of the change on protein function is not known, it could possibly effect photoperiod-dependent flowering time regulation.
The closest Arabidopsis homolog of Lsat_1_v5_ gn_7_97240.1 is COL9. The CO-like gene family threegroup structure is conserved in various plant groups (Wang et al. 2019). Although most genes in the family promote flowering in Arabidopsis, COL9 was found to suppress FT, SOC1, and CO and delay flowering in Arabidopsis (Cheng and Wang 2005). Functional homologs of COL9 affect flowering time in Arabidopsis (Cheng and Wang 2005), rice (Wu et al. 2018), maize (Maldonado et al. 2019) and were also identified in Ambrosia artemisiifolia (Mátyás et al. 2019). In this study, COL9 expression was higher in early flowering genotypes with the P allele, which could suggest a flowering promoting role. On the other hand, expression in S allele genotypes was higher in LT conditions, in which flowering is later. This contrasts with a previous study in which LsCOL9 expression was found to significantly increase under flowering inducing heat conditions (Liu et al. 2018). These discrepancies may be due to unidentified polymorphisms in the genotype's promotors, leading to differing sensitivity to the environmental conditions. So the effect of COL9 gene expression on the flowering phenotype is not yet well defined. Taken together, either LsCOL9 or LsPhyC or both could be qBLT/ qFLT7.2 causal gene involved in lettuce flowering.
In summary, this work contributes to efforts of breeding cultivars less prone to premature bolting through improved resolution and understanding of the major qBLT/qFLT7.2 QTL. This QTL was detected in multiple experimental populations and has a determining effect on bolting and flowering rates in most lettuce production environments.
Another main contribution is the identification of the qBFr2.1 allele from PI 251246 not currently present in the rosette-forming lettuce genepool, as demonstrated with marker LZ001 on a diverse group of cultivars. Since the primitive oil seed cultivar (PI 251246) allele has a significant effect in delaying flowering in SD-HT conditions and in delaying bolting and flowering in certain field conditions (You, Still, Rosental, Hayes and Simko-unpublished data), introgression of the allele into commercial germplasm may lead to development of breeding lines with delayed bolting. The allele may be particularly useful when breeding cultivars specifically adapted for lettuce production in fall plantings in the southwestern USA, an environment especially susceptible to premature bolting and flowering. The LZ001 marker can be used to facilitate breeding by MAS.