Two inbred lines, 12S3 and 12S4, which exhibit significant differences in bolting, were selected as the parents. Line 12S3, with extreme resistance to bolting, was used as the female and recurrent parent, while the early bolting line 12S4 was used as the male parent to develop a BC₁ and an F₂ population (Qian et al., 2017). The two parental lines, the derived F₁ line, and the 148 BC₁ individuals (BC_1a) were planted in a field in spring 2015 for primary mapping. In addition, 200 BC₁ progenies (BC_1b) and 150 F₂ progenies were planted in the same location in spring 2020 in natural conditions for validation of QTLs and narrowing down of the QTL regions. All of these materials were developed and tested by the Spinach Research Group, Institute of Vegetables and Flowers (IVF), and Chinese Academy of Agricultural Sciences (CAAS).

Each individual plant was visually inspected daily, and the bolting date was determined as the date that a stem of a plant was seen to be at least 5 cm in length (Fu et al., 2020). The bolting time (BOT) was then determined as the period from the sowing date to the bolting date. The phenotypic data of all plant materials in these experiments were analyzed with Excel 2013 (Microsoft Corp., Redmond, United States) for calculating the mean, standard error (SE), and coefficient of variation (CV) in each line or population.

At the four true-leaf stage, fresh young leaves were collected from each plant of the F₁, BC_1a, BC_1b, F₂ populations, and parents, immediately frozen in liquid nitrogen, and stored in a –80°C freezer. Genomic DNA was extracted from each plant using the cetyltrimethyl ammonium bromide (CTAB) method (Murray and Thompson, 1980). The DNA concentration and quality were assessed using a ND-2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States) and 1.0% agarose gel electrophoresis, respectively.

Specific length amplified fragment sequencing (SLAF-seq) is an efficient method of large-scale genotyping developed on the basis of high-throughput sequencing technology and reduced representation library (RRL). In brief, an SLAF pilot-experiment was first designed to improve the efficiency of SLAF-seq, which considered the uniform distribution and avoided the duplication of SLAFs. Next, according to the pre-experiment, the SLAF library was conducted as follows The genomic DNA from each sample was completely digested by the two restriction enzymes - RsaI and HaeIII (New England Biolabs, NEB). After digestion, the DNA fragments were repaired with adenine and duplex tag-labeled sequencing adapters. Twenty polymerase chain reaction (PCR) cycles were used to enrich the concentration of fragments and the PCR products were then purified and pooled. The sample was performed by 2% agarose gel electrophresis (120 V, 60 min). After gel purifcation, DNA fragments of 364–414 bp were excised and diluted for paired-end sequencing. Finally, the selected SLAFs were sequenced on an Illumina High-seq 2500 sequencing platform (Illumina, Inc.; San Diego, CA, United States).

The analysis of SLAF-markers followed the procedures described by Sun et al. (2013). All SLAF paired-end reads were clustered on the basis of sequence similarity, which was detected by BLAST (-tileSize = 10, -stepSize = 5). Sequences with over 95% identity were grouped in one SLAF locus. SLAFs with two to four tags were deemed as polymorphic SLAFs.

The SNP molecular markers were obtained from 4080 SLAF markers from the spinach high-density genetic map constructed by Qian et al. (2017), following which a total of 300 KASP primers was designed by the LGC company (Shanghai, China), and the slow bolting parent 12S3 and early bolting parent 12S4 were tested (Liu et al., 2021b). A subset of KASP primers were selected and used to genotype the BC_1a.

For the KASP assays, each sample contained 2.5 μL 2 × KASP Master mix, 0.07 μL KASP Assay mix, and 2.5 μL genomic DNA diluted to 20–30 ng/μL. The reaction system was as follows: 94°C for 15 min, 10 cycles of 94°C for 20 s and 61°C (0.6°C drop per cycle) for 60 s and a further 26 cycles of 94°C for 20 s and 55°C for 60 s. An additional three cycles of 20 s at 94°C and 60 s at 55°C were executed if the results of the initial KASP thermal cycles did not acquire sufficiently defined genotype clusters. In addition to DNA samples, two no-template controls (NTCs) were included on each 384-well PCR plate. All plates were read below 40°C in a 7900 HT Fast Real-Time PCR System (Applied Biosystems), and the data were analyzed using SDS2.3 software (supplied by Applied Biosystems) (Semagn et al., 2014).

The SNP markers were selected with no segregation distortion, and markers with more than 25% missing data were also excluded. The valid markers were then used to construct the linkage map from the BC_1a population using JoinMap 4.0 software (Van Ooijen, 2006). All markers were firstly grouped based on a threshold of LOD = 3.0, while all other settings were left at their default values.

The BC_1b and F₂ populations were used to confirm and narrow down the predicted region, and other KASP markers were developed based on the SNP variation between the two parents around the initial QTL area. The QTLs for bolting were also detected using QTL IciMapping 4.2 software (Meng et al., 2015) based on the phenotype of 148 BC_1a individuals. The Composite Interval Mapping of ADDitive QTL (ICIM-ADD) method was used for QTLs. The parameters were as follows: a step in 1 cM, probability in stepwise regression of 0.001, and LOD = 3.0. The final QTLs were named based on the method of McCouch et al. (1997): “ q” + the English abbreviation of the trait + the chromosome number + “ . “ the QTL number.

Based on the of spinach genome annotations (version Sp75) in SpinachBase,¹ the genes related to bolting and flowering within the identified interval were selected for further analysis. The full-length RNA was extracted at the 12-leaf-stage and the promoters of the candidate genes were sequenced between the two parents. The specific primers were designed by Primer3 plus² (Table 2). The candidate genes were cloned and the sequences were aligned by MUSCLE software.³ Finally, the gene structure was elucidated based on the re-sequenced result.⁴

Quantitative real-time PCR was employed to evaluate the expression of the candidate genes from the seedling to bolting stages. Leaf tissue of the 12S3 and 12S4 lines was collected at 6, 9, 12, 15, and 18 weeks until both parents began bolting in spring of 2021 (12S3 bolted 18 weeks after planting; 12S4 bolted 15 weeks after planting) (Figure 7B), and the total RNA was extracted using a Plant Total RNA Mini Kit (GeneBetter Biotech, Beijing, China⁵). The cDNA was synthesized from 500 ng total RNA with a TranScript One-Step gDNA Removal and cDNA Synthesis Kit (TransGen Biotech, Beijing, China⁶). Three independent biological and three technical replicates of each period were performed and analyzed. The synthesized cDNA was subjected to quantitative real-time (qRT)-PCR analysis using a QuantStudio™ 12 K Flex Real-Time PCR System (Applied Biosystems) with SYBR Fast qPCR Mix (TaKaRa⁷). The reaction mixture contained 70 ng template cDNA, 0.2 μM of gene-specific primer (Table 2), 0.2 μM ROX Reference Dye II, 3.4 μL ddH₂O, and 5 μL 2 × SYBR Fast qPCR Mix in a 20 μL volume. The qRT-PCR was performed at 95°C for 30 s, followed by 35 cycles of 95°C for 15 s and 60°C for 1 min. The relative expression was calculated using the 2^–ΔΔCT method. SpActin was used as the reference gene (Lee and Zeevaart, 2002).

In 2015, the bolting time of line 12S3 and 12S4 was on average 62 (60–65) and 46.5 (45–48) days, respectively, indicating differences in the bolting time of the parents. From the BC_1a line, the bolting time ranged from 48 to 66 days, with an average of 54.5 days. The bolting time of 155 BC_1b and the 123 F₂ individuals was from 47 to 65 days, with a mean value of 57.6 days in BC_1b and 55.8 days in the F₂ populations. The phenotypic traits are summarized in Table 1 and Supplementary Table 1. Moreover, these segregating populations showed continuous variation in bolting time, suggesting that the bolting trait has a quantitatively inherited character in spinach (Figure 1 and Supplementary Table 2).

In our previous study (Qian et al., 2017), a total of 4080 SLAF markers for 148 BC_1a individuals were acquired by SLAF-seq, and the linkage groups were coded with six linkage groups (P01–P06) in a total length of 1125.97 cM, which matches the spinach chromosome numbers (Supplementary Table 3). After the exclusion of missing and disqualified data, 130 BC_1a individuals were finally used to map the QTLs (Supplementary Table 3). Combining the SLAF high-density genetic map with bolting time in 130 BC_1a progenies, a major QTL (named qBT1.1), which contributed 49.07% of the phenotypic variance (PVE) (Table 1), was identified at the interval 15.82–18.97 cM on LG3 between two adjacent SLAF markers (Marker 2552708 and Marker 1611427), with an LOD score of 16.39. Based on the spinach genome Sp75 (Xu et al., 2017), this QTL was mapped on chromosome 1 in the region of 47.72–50.61 Mb (Figure 2).

A total of 181 informative SNP markers (Liu et al., 2021b) and 147 BC_1a individuals were selected for KASP-based linkage analysis (Supplementary Table 4). Based on the 181 KASP markers and screened 127 BC_1a plants, two QTLs (qBT1.2 and qBT1.1) were mapped to 103.5–105.5 cM and 163.5–166.0 cM on LG3 (Figure 3) and were located at 41.44–42.02 Mb and 46.76–49.12 Mb on chromosome 1, respectively. The LOD scores were 3.46 and 13.39, explaining 8.90 and 40.86% PVE, separately (Table 1). In 2020, the same 181 KASP markers were used in 185 BC_1b and 112 F₂ populations and they were co-located in the same area between KM706861 and KM3309304 (Figure 3). The LOD score was 10.5623 and explained 41.99% PVE in BC_1b individuals while the figures were 19.9202 and 51.20% in F₂ plants, respectively (Table 1).

The SLAF-based and KASP-based QTLs in the two strategies indicated that qBT1.1 was a stable locus that could be used for fine-mapping and cloning. The SNP variations were explored in the sequences at this region of 12S3 and 12S4. The raw reads were first-filtered by fastp 0.12.0 (Chen et al., 2018) and the alignment data were obtained on the spinach genome Sp75 (Xu et al., 2017) by BWA 0.7.17-r1188 (Li and Durbin, 2009). The vcf files were finally generated by Samtools/Bcftools 0.1.19 – 44428 cd (Li et al., 2009). To further refine the mapping region, the KASP markers were developed from 40 to 51 Mb on chromosome 1 by the file.

After eliminating the invalid segregation data, 20 efficient KASP markers were designed for fine mapping (Supplementary Table 5). In the expanded BC₁ population (BC_1b), 185 individuals were obtained in 2020, and then 16 recombinant plants were ultimately acquired by the KASP genotyping. A major QTL (qBT1.1) for bolting time was verified between KMBL53 (31.0 cM) and KM3309304 (31.5 cM) (Figure 4). Furthermore, in the F₂ population in 2020, the 112 plants were used to map the QTLs for bolting time, and nine recombinant individuals were obtained. A major QTL in the F₂ populations was also fine-mapped in the interval of KMBL53 (32.4 cM) and KM3309304 (33.1 cM). In conclusion, a QTL named qBT1.1 was detected with a 0.56-Mb region between KMBL53 (47.56 Mb) and KM3309304 (48.12 Mb) on chromosome 1 (Figure 4).

A total of 68 genes were located in a 560 kb region based on the spinach genome (version Sp75)⁸ (Supplementary Figure 1 and Supplementary Table 6). Among these genes, two genes were unannotated, 23 genes encoded various enzymes, and four genes had transmembrane structure. Two genes Spo04942 and Spo04967 were found to be homologs of bolting and flowering genes in Arabidopsis (Table 2, Figure 5, and Supplementary Table 6), and thus could potentially be the candidate genes controlling the spinach bolting trait.

Spo04942 is MADS-box transcription factor that is homologous to Arabidopsis FLC and the sugar beet FLC homolog FLC-LIKE 1, thus was renamed as SpFLC. In the coding area, there were two synonymous SNP variations and one non-synonymous SNP variation that led to the change from tyrosine (12S3) to asparagine (12S4) at position 98 in the domain area (Figure 6). In the 2-kb upstream non-coding sequences, no difference was found between the two parents. The gene Spo04967 may encode a zinc finger protein similar to CONSTANS-LIKE 14, which belonged to the COL family, and thus was named as SpCOL14. With Sanger sequencing, no variation was found in the coding region of SpCOL14 between 12S3 and 12S4. However, from the 2-kb upstream non-coding region, variations of about 900 bp were found between 12S3 and 12S4 (–730 bp to –1653 bp) (Figure 6) that may affect the expression of SpCOL14.

Except for SpFLC and SpCOL14, no other genes were found to be related to bolting and flowering genes, such as FT, SOC1, and FLOWERING LOCUS D (FD) that were reported in other plants (Hideki et al., 2021). Interestingly, in this 0.56 Mb region, we found some transcription factors that regulate the various stages of plant growth and development, including Spo04911, which is an NAC domain-containing protein that plays a role in regulating plant growth and stress resistance, and Spo04943, which is a MADS-box factor with an AGAMOUS-LIKE 9 homolog to floral organ development. The factors may play some minor roles in affecting bolting regulation, but the gene regulatory network will be identified in the future with technological progress.

We assessed the expression patterns using qRT-PCR analysis between 12S3 and 12S4 at different growth stages to further assess the two candidate genes (Table 2). 12S4 bolted after 15 weeks, while the bolting time of 12S3 was after 18 weeks (Figure 7B). Although the expression of these two candidate genes at the five stages were significantly different, a similar change trend of expression levels of SpFLC was found both in 12S3 and 12S4, with a high level observed in both at week 9. However, the expression level of SpCOL14 between the 12S3 and 12S4 plants was associated with differences in the phenotype; in 12S4 (Figure 7A), the expression of SpCOL14 showed a high level at week 9, while 12S3 showed high expression at week 12 (Figure 7C). These results suggested that SpCOL14 could potentially be the key candidate gene controlling bolting in spinach.