Fine mapping and analysis of candidate genes for qFT7.1, a major quantitative trait locus controlling flowering time in Brassica rapa L

qFT7.1, a major QTL for flowering time in Brassica rapa was fine-mapped to chromosome A07 in a 56.4-kb interval, in which the most likely candidate gene is BraA07g018240.3C. In Brassica rapa, flowering time (FT) is an important agronomic trait that affects the yield, quality, and adaption. FT is a complicated trait that is regulated by many genes and is affected greatly by the environment. In this study, a chromosome segment substitution line (CSSL), CSSL16, was selected that showed later flowering than the recurrent parent, a rapid-cycling inbred line of B. rapa (RcBr). Using Bulked Segregant RNA sequencing, we identified a late flowering quantitative trait locus (QTL), designated as qFT7.1, on chromosome A07, based on a secondary-F2 population derived from the cross between CSSL16 and RcBr. qFT7.1 was further validated by conventional QTL mapping. This QTL explained 39.9% (logarithm of odds = 32.2) of the phenotypic variations and was fine mapped to a 56.4-kb interval using recombinant analysis. Expression analysis suggested that BraA07g018240.3C, which is homologous to ATC (encoding Arabidopsis thaliana CENTRORADIALIS homologue), a gene for delayed flowering in Arabidopsis, as the most promising candidate gene. Sequence analysis demonstrated that two synonymous mutations existed in the coding region and numerous bases replacements existed in promoter region between BraA07g018240.3C from CSSL16 and RcBr. The results will increase our knowledge related to the molecular mechanism of late flowering in B. rapa and lays a solid foundation for the breeding of late bolting B. rapa.


Introduction
The economically important crop, Brassica rapa, has long been cultivated worldwide, mainly as a vegetable foodstuff, such as Chinese cabbage and Pak-choi, and to a lesser extent to produce fodder and oilseed, such as turnip rape and yellow sarson (Carpio et al. 2011). Among the agronomic traits in B. rapa, flowering time (FT) is important because it affects the yield of seeds and the harvested crop's commercial quality (Wu et al. 2012), thus determining their growing season and cultivation area (Xiao et al. 2019). The regulation of FT is complex, involving multiple genes (quantitative trait loci (QTLs)) and is markedly affected by environmental conditions, making it a challenge to identify linked markers or related genes for marker assisted selection (MAS)-based breeding (Liu et al. 2016).
Many QTLs related to FT in B. rapa have been identified in the last 30 years, mainly based on a wide range of biparental mapping populations or natural populations (Teutonico and Osborn 1995;Osborn et al. 1997;Ajisaka et al. 2001;Nishioka et al. 2005;Yang et al. 2007;Lou et al. 2007Lou et al. , 2011Li et al. 2009aLi et al. , 2015Yuan et al. 2009;Zhao et al. 2010;Kakizaki et al. 2011;Xiao et al. 2013Xiao et al. , 2019Zhang et al. 2015;Liu et al. 2016;Xi et al. 2018;Wang et al. 2018b;Fu et al. 2020;Kaur et al. 2021). Most of these populations are traditional primary populations, including F 2 , BC 1 , doubled haploid (DH), and recombinant inbred lines (RILs). However, these populations are only useful to detect QTLs with relatively large effects, because the segregants often have complicated backgrounds because of the presence of large parent-derived chromosomal fragments, thus QTLs Communicated by Isobel AP Parkin.
1 3 with minor effects might be missed. By contrast, a wider range of QTLs can be identified using advanced mapping populations, including near-isogenic lines (NILs) and chromosome segment substitution lines (CSSLs) (Nadean and Frankel. 2000). Moreover, secondary F 2 or F 3 populations can be produced by from backcrossing a target NIL or CSSL with the recurrent parent, making them more suitable for fine mapping and positional cloning of a target QTL (Yano. 2001). Over 75 CSSL libraries in 17 major crops have been constructed in the last three decades, which have made significant contributions to QTL characterization (Balakrishnan et al. 2019), despite the development of these population being labor and time intensive. However, currently, only a few CSSLs are available for B. rapa Wang et al. 2018b). Most of the above-mentioned QTLs identified in primary populations have not yet been fine mapped, mainly because of a lack of optimal genetic material, and very few studies resulted in the actual cloning of the gene responsible for flowering in B. rapa. To the best of our knowledge, to date, only a few genes, such as BrVIN3.1, BrFLC1, BrFLC2 (Su et al. 2018;Jeong et al. 2019), and BraELF6 , have been successfully cloned and subsequently confirmed by transformation in B. rapa and in Arabidopsis thaliana, respectively. Furthermore, transcriptional analysis (RNA-sequencing (RNA-seq)) based on CSSLs (or NILs) also has advantages for identifying candidate genes underlying QTLs (Keurentjes et al. 2007).
In Arabidopsis, six major pathways were identified that regulate flowering time, including the photoperiod pathway, the vernalization pathway, the gibberellin pathway, the autonomous pathway, the ambient temperature pathway, and the age pathway (Fornara et al. 2010). Flower locus T (FT), as the core gene of the flowering pathway, is located downstream of the critical gene of the photoperiod pathway CONSTANS (CO). CO activates the expression of FT, and the bulk of the FT protein in leaves moves via the phloem to the shoot apical meristem (SAM), where it combines with FLOWERING LOCUS D (FD) and promotes plant early flowering (Abe et al. 2005;Corbesier et al. 2007). The FT family includes five other genes: TERMINAL FLOWER 1 (TFL1), MOTHER OF FT AND TFL1 (MFT), TWIN SISTER OF FT (TSF), BROTHER OF FT AND TFL1 (BFT), and ARABIDOPSIS THALIANA CENTRORADIALIS (ATC ) . FT, TSF, and MFT facilitate the transition to reproductive development and flowering, whereas TFL1, BFT, and ATC suppress this process (Yoo et al. 2004;Wickland and Hanzawa, 2015). In addition, TFL1 competes with FT for FD to reduce the expression of the SAM key gene LEAFY (LFY), thus suppressing early flowering in Arabidopsis (Wigge et al. 2005;Zhu et al. 2020). ATC is a TFL1-like gene in Arabidopsis that is homologous with CEN-TRORADIALIS (CEN) and was first identified in Antirrhinum (Bradley et al. 1996). ATC inhibits flowering and regulates inflorescence morphology (Banfield and Brady 2000). ATC encodes a protein that is 77% similar to Antirrhinum CEN and 67% similar to TFL1 at the amino acid level. ATC , CEN, and TFL1 overexpression showed similar phenotypes in wildtype Arabidopsis (Ratcliffe et al. 1998;Mimida et al. 2001). ATC acts in a non-cell autonomous manner to inhibit flowering in Arabidopsis and is specifically expressed in the hypocotyl, through long distance movement from the SAM. LFY and AP1 are the critical factors in the SAM. ATC might inhibit the expression of LFY and AP1 to delay flowering (Hempel et al. 2000;Huang et al. 2012;Gao et al. 2017). Many ATC homologous genes have been identified in different species, such as GoCEN-Dt in cotton, Hordeum vulgare CENTRORADIALIS (HvCEN) in barley, SELF-PRUNING (SP) in tomato, and ZEA CENTRORADIALIS (ZCN) in maize (Pnueli et al. 1998;Danilevskaya et al. 2010;Comadran et al. 2012;Liu et al. 2018).
Recently, the identification and isolation of genes underpinning QTLs associated with agronomic traits in crops have been accelerated significantly because of the emergence of rapid and inexpensive next-generation sequencing (NGS)based technologies combined with plant genetics (Nguyen et al. 2019). Traditional QTL mapping integrated or combined with QTL-seq and RNA-seq is a highly efficient and accurate approach for QTL mapping and validation, which enables the identification of candidate genes associated with agronomic traits of interest, which has been widely utilized in diverse crops (Lu et al. 2014;Berenguer et al. 2015;Gelli et al. 2017;Li et al. 2017;Shu et al. 2018;Wang et al. 2018a;Liu et al. 2019;Park et al. 2019;Wen et al. 2019;Huang et al. 2021).
Previously, a set of CSSLs was developed by our group using the rapid-cycling inbred line of Brassica rapa (RcBr) and the Chinese cabbage inbred line 08A061 as the recipient and donor parents, respectively ). Among the developed CSSLs, CSSL16 showed significantly later flowering than the recurrent parent, RcBr. To map and identify the candidate gene(s) responsible for this late flowering, we developed a secondary F 2 population derived from a cross between CSSL16 and RcBr. Bulked Segregant RNA-Seq (BSR-Seq) identified a QTL, qFT7.1, which was validated using conventional QTL mapping. Ultimately, recombinant analysis narrowed down qFT7.1 to a 56.4-kb interval on chromosome A07, allowing the candidate genes to be identified. Our findings represent a benchmark to further determine the molecular mechanism of late flowering in B. rapa.

Plant materials and trait measurement
Our laboratory previously constructed a set of CSSLs populations from a cross between RcBr as the recipient parent, 1 3 which is an extremely early-flowering and vernalizationindependent inbred line, and 08A061 as the donor parent, which is a very late-flowering and vernalization-dependent Chinese cabbage inbred line ). In the CSSL population, one line, CSSL16, with a late flowering phenotype, was chosen for backcrossing with the recurrent parent, RcBr. The secondary F 2 CSSL16/RcBr population was subsequently constructed using self-pollination.
The phenotypic analyses were carried out at the Experiment Station of Shenyang Agriculture University, Shenyang, China (41.8°N, 123.4°E). Four flowering indices were assessed to calculate the phenotypic scores of the individuals: The bolting index (BI), days to reach a 5/10-cm-high elongated floral stalk (DE5/DE10), and FT (Liu et al. 2016). CSSL16 and RcBr were evaluated under four environmental conditions (E1, E2, E3 and E4; Table S1). A secondary F 2 population consisting of 500 individuals was sown in the greenhouse in September 2019, which were utilized for BSR-Seq. Whereas 300 individuals were used for conventional QTL analysis in March 2020. The progeny of the recombinant individuals screened from the secondary F 2 population (2200 individuals), along with parental controls, was grown in March 2021 for fine mapping of the identified QTL. All plants were sown directly into 10-cm pots without providing any extra vernalization, as described in our previous study ).

RNA isolation and extreme pool construction
For BSR-Seq analysis, we constructed two extreme pools, an L-pool (late-flowering pool) and an E-pool (early flowering pool), by mixing the same amounts of RNA from 25 late-flowering or 25 early-flowering plants, according to the phenotypic scores of 500 F 2 individuals recorded in the fall of 2019. About 50 days after sowing, the apex leaves of each pool were sampled and subjected to RNA isolation. The TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract the total RNA. An Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), a NanoDrop spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA), and 1% agarose gel were used to assess the quantity and quality of the extracted RNA. RNA (1 μg) with an RNA integrity number (RIN) > 7 was then processed for next generation sequencing library construction (NEBNext ® Ultra™ RNA Library Prep Kit for Illumina®; NEB, Ipswich, MA, USA).
According to the basic principle of the ED value, the occurrence frequency of the four bases A, T, C, and G at the SNV site was statistically different in the population, and the corresponding base frequency of the two trait groups was calculated by distance. To eliminate the background noise, the ED value of each different SNV site was raised to the power of five, termed ED 5 (Su et al. 2016). All ED 5 values were sorted, and the different SNV sites corresponding to the top 1% of ED 5 values were screened, and then, specific chromosome segments were located according to the distribution of the different SNV sites.

DNA extraction and marker development
The CTAB method, with minor modifications (Murray and Thompson 1980), was used to extract total DNA from the two parental lines and individuals of the secondary F 2 populations. The PCR reaction volume and amplification were same as those described previously ). The primers for the InDel markers were designed using Primer Premier 5.0 software (Premier Biosoft, San Francisco, CA, USA) based on sequence variations of the target region identified by BSR-seq and whole-genome re-sequencing between RcBr and CSSL16. In addition, some primers were used that amplified simple sequence repeat (SSR) markers, which were used previously in our laboratory. The primer information is shown in Table S2.

QTL analysis and fine mapping
The BSR-seq-identified QTL for FT was confirmed using classical QTL mapping, assisted by polymorphic markers. The secondary F 2 populations utilized for conventional QTL analysis consisted of 300 individuals sown in March 2020. QTL mapping was performed using composite interval mapping (CIM) implemented in Windows QTL Cartographer 2.5 (Silva et al. 2012). The determined logarithm of odds (LOD) value for putative QTL declaration was determined after 1000 permutation tests at a significance level of P < 0.05 and a threshold of 3.0.

3
The progeny of recombinant individuals screened from a larger secondary F 2 population (2200 individuals) sown in a greenhouse in March 2021 was used for fine mapping of the identified QTL. The means of the homozygous recombinant phenotype of the progeny were analyzed using Student's t-test in SPSS v17.0 (IBM Corp., Armonk, NY, USA) and compared with that of the control (RcBr) at a significance level of P < 0.01.

Whole-genome resequencing
Total DNA was extracted from young leaves. Its quality was determined using 0.8% agarose gel electrophoresis, and it was quantified using an ultraviolet spectrophotometer. The Illumina NovaSeq platform was used for 2 × 150 bp pairedend sequencing. The raw data were cleaned using Adap-terRemoval (version 2) (Schubert et al. 2016), and high quality reads were compared with the Brassica rapa V3.0 reference genome. Single nucleotide polymorphisms (SNPs) and InDels were detected using the GATK software (Van der Auwera et al. 2013) and analyzed using the ANNOVAR software (McKenna et al. 2010;Wang et al. 2010). SNP and InDel loci were used for population-specific locus analysis. The ED value was used to identify difference regions between the two parents and was displayed as an image.

Candidate gene sequence analysis
The annotation information of genes in the candidate region was obtained from the Brassica rapa database (BRAD; http:// brass icadb. cn/#/ Downl oad/) and The Arabidopsis Information Resource (TAIR; http:// www. arabi dopsis. org/). The specific primers (Table S2) to amplify the full-length coding sequences and promoter sequences were designed using Primer Premier 5.0. A Gel Extraction Kit (CWBIO, Beijing, China) was used to purify the PCR products, which were ligated into a pGEM-T Easy Vector (Promega, Madison, WI, USA), followed by sequencing at GENEWIZ. The sequences were aligned using the DNAMAN software (Lynnon Biosoft, San Ramon, CA, USA), and the structure of the candidate gene was displayed using online software (http:// gsds. cbi. pku. edu. cn/).

Quantitative real-time reverse transcription PCR (qRT-PCR) analysis of candidate gene expression
The expression level of the candidate gene was detected using qRT-PCR. Total RNA of RcBr and CSSL16 from roots, leaves, cotyledons, hypocotyls, stems, flowers, and the SAM were isolated using an RNA extraction kit (Aidlab Biotechnologies Co., Ltd., Beijing, China). The RNA was then reverse transcribed to cDNA. The cDNA was used as the template for the quantitative real-time PCR (qPCR) step of the qRT-PCR protocol (reaction volume: 25 µL, comprising 2 × Ultra SYBR Mixture, 2 µL of diluted cDNA, 1 µL of 0.2 μM forward and reverse primers, and 21 µL of RNasefree water). The reaction conditions were: initial denaturation at 95 °C for 10 min, followed by 40 cycles of 95 ℃ for 15 s and 60 °C for 1 min. This was followed by meltingcurve analysis: 95 °C for 15 s, 60 °C for 1 min, 95 °C for 15 s, and 60 °C for 15 s. The 2 −△△Ct method was used to analysis the relative expression level. Cycle threshold (C t ) values were shown as the means of three independent biological replicates. Each sample was analyzed as three independent technical replicates. QuantStudio™6 Flex Manager software (Livak and Schmittgen 2001) was used to analyze the data. Primer Premier 5.0 was used to design gene-specific primers (shown in Table S2), with the Actin gene being used as the internal control .

Genotypic and phenotypic characterization of RcBr and CSSL16
To detect the segment that had introgressed from 08A061 on ten chromosomes of CSSL16, the two parental lines, RcBr and CSSL16, were genotyped using whole-genome resequencing. A total of 67,507,156 and 64,540,134 high-quality reads were detected in RcBr and CSSL16, and the clean data were compared with the reference genome. This identified 18,894 SNPs and 4,469 InDels on the ten chromosomes between the two parental lines. According to the ED value calculation, the variation was mainly distributed in chromosomes A02 (physical location 1,775,235-2,512,196 bp) and A07 (physical location 15,350,379-16,648,887 bp) (Fig. 1). The total substituted segment derived from the donor parent, 08A061 was approximately 2.04 Mb, and the background recovery rate was about 99.42% (351.10/353.14). To help the interpretation of the nature of the introgression in the recurrent parent, flowering-related genes of two substituted segments are listed in Table S4.
RcBr and CSSL16 showed a significant difference in FT under multiple environments (Fig. 2a) (Fig. 2b). The two parental lines also showed significant differences in DE5, DE10, and BI under all four growth conditions (Table 1). In conclusion, RcBr and CSSL16 showed significant differences in all flowering-related traits.

Identification and validation of qFT7.1
Through BSR-seq analysis, we were able to map 45,425,180 and 41,697,006 clean reads to the B. rapa reference genome from the E-pool and L-pool, respectively. A total of 218,944 SNPs in the E-pool and 209,924 SNPs in the L-pool were identified. All ED 5 values were sorted, the top 1% of the ED 5 values was used as the threshold, and the different corresponding SNV sites were screened. The distribution of different SNV sites confirmed the candidate region. This candidate region of the QTL for FT was located on chromosome A07, starting at 15,486,952 and ending at 16,546,846; thus, the candidate interval covers about 1.06 Mb (Fig. 3). The candidate QTL underlying FT in this region was designated as qFT7.1.
To validate qFT7.1, identified by BSR-seq analysis, we carried out conventional QTL analysis with 300 F 2 individuals in March 2020. The frequency distribution of the F 2 population presented a normal distribution and showed Scale bars = 5 cm; **Significant at P < 0.01 the phenomenon of transgressive segregations (Fig. S2). A total of 13 polymorphic markers (Table S2) were screened from the difference interval (the donor segment of 08A061), which were detected using whole genome re-sequencing (Chromosomes A02 and A07), and all polymorphic markers were used for classical QTL mapping. One QTL with a LOD value of 32.2, explaining 39.9% of the phenotypic variation, was found to control FT and was located between marker InDel714 and InDel716, corresponding to a physical position of 15,539,588 bp to 16,499,043 bp on chromosome A07 (Fig. 4b). However, we could not detect any QTL on chromosome A02. Thus, the conventional QTL analysis confirmed the QTL qFT7.1, which was identified via BSRseq analysis.

Fine mapping of qFT7.1
The candidate QTL, qFT7.1, was preliminary mapped to a 1.06 Mb candidate region on chromosome A07. A larger F 2 population consisting of 2200 individuals were used to refine the position of qFT7.1. Recombinant plants were screened with markers InDel723 and InDel716, and a total of 19 recombinant individuals were obtained. The homozygous recombinant progeny was divided into ten groups according to their genotypes. The mean value of the homozygous progeny (DE5, DE10, and FT) was compared with RcBr at the P < 0.01 level.
Recombinants L1 and L2 both showed no difference with RcBr, while L3 and L4 had the opposite genotype to L1 and L2. Recombinants L3 and L4 both showed significant differences with RcBr, but no difference with CSSL16; thus, these groups placed the QTL in a region upstream of markers InDel705 and InDel708, respectively. In the same way, the results for recombinants L5 and L6 demonstrated that the QTL was located upstream of markers InDel706 and InDel707, and the results for recombinant L7 placed the QTL in a region upstream of marker InDel702. Furthermore, recombinant L8 identified that the QTL was located upstream of marker InDel715, whereas recombinant L9 was significantly different from CSSL16 and similar to RcBr, which delimited the QTL to a region downstream of marker InDel714. Recombinant group L10 showed a significant difference with RcBr, thus group L10 placed the QTL in a region downstream of marker InDel723. By positioning of the recombinant groups, qFT7.1 was finally narrowed down to a 56.4-kb interval between marker InDel714 and InDel715, representing the physical position of 15,539,588 bp to 15,595,959 bp on chromosome A07 (Fig. 4c). 19.23 ± 1.12 24.32 ± 0.65** 23.24 ± 1.73 29.96 ± 2.04** 38.78 ± 2.78 51.55 ± 3.37** 7.00 ± 0.00 5.00 ± 0.00** Fig. 3 The distribution of the ED 5 value of differential SNPs on Brassica rapa chromosomes according to BSR-Seq analysis. BSR-Seq-based distribution of SNPs on chromosomes. The x-axis shows the 10 B. rapa chromosomes, and the y-axis shows the ED 5 values of the filtered SNPs; the horizontal line is the threshold of the top 1% Fig. 4 Fine mapping of qFT7.1. a Graphic representation of the genotype of chromosome A07 encompassing qFT7.1. Sequence variations of chromosome A07 between the two parental lines detected by whole genome re-sequencing. The orange region represents the confidence interval of qFT7.1 identified by BSR-seq. b Traditional QTL mapping was preformed to validate the QTL qFT7.1. Physical map of the qFT7.1 region on chromosome A07. Traditional QTL analysis was used to validate qFT7.1, which was preliminary located between marker InDel714 and InDel716. The position of markers is shown on the x-axis, and the LOD value is shown on the y-axis. The LOD value of qFT7.1 was 32.2, which explained 39.9% of the phenotypic variation. c Fine mapping of qFT7.1. The phenotype and genotype of the ten homozygous recombinant groups and the two parental lines (RcBr and CSSL16) used for fine mapping. The marker genotypes of RcBr are shown as black bars and those of CSSL16 are shown as white bars; DE5, DE10, and FT data appear as means ± SD. Significant differences for the traits of the recombinant compared with those of the parents are indicated using superscript letters (a, b). Student's t test was used to distinguish significant difference at P < 0.01. d Structure of the BraA07g018240.3C coding region. Whole genome re-sequencing and cloning detected sequence variation of BraA07g018240.3C between the two parental lines; the black regions represent exons, and the straight lines represent intron; two base variation mutations were identified in the first and second exons

Candidate gene annotation
According to the B. rapa reference genome database annotation, nine genes were annotated to the 56.4-kb region (BraA07g018220.3C-BraA07g018300.3C). The detailed information for these genes is shown in Supplementary  Table 3. All the genes were compared with Arabidopsis homolog genes and analyzed for their function. We found that BraA07g018240.3C is homologous with Arabidopsis gene At2g27550, a key gene regulating FT. This gene was an Arabidopsis CENTRORADIALIS homolog (ATC ) gene, which belongs to FLOWERING LOCUS T (FT) family and encodes a protein similar to TERMINAL FLOWER1 (TFL1), the overexpression of the gene encoding which leads to a similar phenotype as TFL1 overexpression. The encoded protein from the identified gene might inhibit the expression of critical flowering genes LFY and AP1, acting in a non-cell autonomous manner to delay flowering (Huang et al. 2012;Zhu et al. 2020).

Expression analysis by qRT-PCR
The candidate gene expression level was detected using qRT-PCR, which indicated that BraA07g018240.3C was specifically expressed in the root and hypocotyl, and not in other tissues. The expression level showed significant differences in the hypocotyl, but not in the root, between RcBr and CSSL16 (Fig. 5). In addition, we could not detect any expression of the other paralogs of BraA07g018240.3C in different tissues using qRT-PCR, including in the SAM and hypocotyl, which was verified using RNA-seq (Table S5). The expression of the other eight genes was also detected in the root, stem, leaf, flower, and hypocotyl, with BraA07g018270.3C and BraA07g018300.3C showing significantly different expression levels in flowers (Fig. 6). The signals of each flowering pathway were collected in the SAM and were used together to determine the FT. AP1 and LFY are both main inflorescence meristem genes and play a central role in the flowering regulatory network (Wellmer and Riechmann 2010). To verify the most likely candidate gene, the expression of LFY and AP1 homologous genes in the SAM were detected in RcBr and CSSL16. LFY homologous genes included BraA02g043220.3C and BraA06g025360.3C, and the AP1 homologous gene included BraA02g018970. 3C,BraA07g030470.3C,and BraA07g034100.3C in B. rapa. Five specific primers,, were used to analyze the expression levels of the AP1 and LFY genes (Table S2). The results showed that the expression levels of all the LFY and AP1 homologous genes were significantly different between RcBr and CSSL16, with all the genes being downregulated in CSSL16 (Fig. 7). Higher expression of LFY and AP1 resulted in Arabidopsis premature flowering; therefore, the protein encoded by ATC (BraA07g018240.3C) possibly downregulates the expression of LFY and AP1, which would lead to delayed flowering, similar to the function of ATC in A. thaliana. In conclusion, we predicted that BraA07g018240.3C was the most likely candidate gene. To further confirm the candidate gene, we analyzed the sequence variations of BraA07g018240.3C between the two parental lines.

Sequence variation analysis of BraA07g018240.3C
A map-based cloning method was used to clone the sequence of BraA07g018240.3C from the two parents, and the sequencing results were analyzed using DNA-MAN software. To identify variations in the candidate gene sequence, a specific primer, 24-C, was used to detect CDS sequence variation (Table S2). The full length gene for BraA07g018240.3C is 1600 bp, starting at 15,558,430 and ending at 15,560,029, including three introns and four exons. The CDS sequence had an A to T mutation at the 12 th base in first exon and a base T to C change at the 32nd base in the second exon; however, both mutations were synonymous (Fig. 4d). Two specific primers, QG-1 and QG-22 (Table S2), were used to amplify the promoter. The result showed many changes in the promoter regions of BraA07g018240.3C between the two parental lines (Fig.  S1).

Discussion
In this study, we employed BSR-seq based on secondary F 2 populations derived from CSSL16 and RcBr to identify the QTL qFT7.1 (Fig. 3), which is responsible for late flowering in B. rapa, and was further validated using classical QTL mapping (Fig. 4c). QTL-seq (BSR-seq) combined with classical QTL mapping has proven to be a powerful tool to identify major QTLs controlling traits of interest in a variety of crops (Lu et al. 2014;Berenguer et al. 2015;Gelli et al. 2017;Shu et al. 2018;Park et al. 2019;Wen et al. 2019;Huang et al. 2021). Thus, most of the populations utilized for QTL-seq are preliminary populations, such as F 2 , DHs, and RILs. QTL-seq is utilized mainly to detect major QTLs with large effects, whereas QTLs with minor effects might not be detected by QTL-Seq, for which traditional QTL mapping is more suitable. The background is basically the same as the recurrent parent, and the positioning accuracy was relatively high.
Ultimately, qFT7.1 was fine mapped to a 56.4-kb interval, between the two InDel markers, InDel714 and InDel715, (Fig. 4c) and a physical position of 15,539,588 to 15,595,959 on chromosome A07. In our previous studies, we did not detect any QTLs in the candidate region basing on F 2 , RIL, and CSSLs derived from the identical parents, RcBr and 08A061 Liu et al. 2016). The CSSLs were constructed using 166 InDel and SSR markers that were distributed relatively evenly on the ten chromosomes; however, a low marker density is likely to lead to small introgression segments being missed. We did not detect any introgression segments on chromosome A07 for CSSL16 based on a limited number of markers (data not shown); therefore, we re-sequenced the two parental lines to identify the possible segment derived from 08A061 (Fig. 1). The nature of the bolting or flowering trait in B. rapa is complex and highly influenced by environmental factors, which perhaps is another reasons why no flowering-related QTL could be detected in the F 2, RIL, and CSSLs from the A07 candidate region. Until now, no other flowering-related QTLs have been detected in the candidate interval of qFT7.1 on chromosome A07 in B. rapa, allowing us to identify the Arabidopsis ATC homologous gene for the first time, which is of a great significance to breed late flowering varieties of B. rapa.
BraA07g018240.3C was homologous with the Arabidopsis gene ATC (At2g27550), which belongs to the FT family and acts systemically to inhibit floral initiation. We found that BraA07g018240.3C was expressed specifically in the hypocotyl and root, which was consistent with the results of Huang et al. (2012): in Arabidopsis, the ATC gene is mainly expressed in vascular tissues, but not in the apex. In the present study, we could not detect any expression (the number of reads was zero) of BraA07g018240.3C based on BSR-seq, which proved that BraA07g018240.3C is specially expressed in the root and hypocotyl of B. rapa further. According to our results, the expression of AP1 and LFY genes in the SAM was significantly higher in RcBr than in CSSL16, which was consistent with the results reported by Liu et al (2009) and Kaneko-Suzuki et al. (2018), allowing us to speculate that ATC might inhibit AP1 and LFY expression positively via long distance transport in the SAM, which then delays flowering in B. rapa. Some studies indicated that ATC homologous proteins have similar functions and their upregulated expression delays flowering. Our study also found that the BraA07g018240.3C expression level in CSSL16 was higher than that in RcBr, which was consistent with the prediction that upregulated expression of this ATC -like gene would delay the flowering; however, this function needs to be verified in transgenic plants. In addition, we found no significant difference in the expression of BraA07g018240.3C in the roots between the two parents, possibly because of the high expression of the ATC gene in the hypocotyl (Fig. 5). Previous studies have shown that in Arabidopsis, the ATC gene is specifically expressed in vascular tissue, while in rice and maize it is mainly expressed in the leaf and stem (Lazakis et al. 2011;Huang et al. 2012;Kaneko-Suzuki et al. 2018). The expression of the ATC gene was significantly downregulated under short day conditions, but showed no significant expression difference under long day conditions (Huang et al. 2012). Taken together, the results showed that the lack of a significant difference in expression of BraA07g018240.3C in the roots between the two parents might be caused by photoperiod changes under natural conditions. We cloned and sequenced the candidate gene BraA07g018240.3C promoter and CDS sequence. This identified two synonymous mutations in the CDS and some mutated bases in the promoter sequence. Changes in the promoter region might influence gene expression. For example, tomato SELF PRUNING 5G is a critical gene for FT, and Zhang et al. (2018b) found that changes to the promoter region resulted in delayed flowering. The candidate gene for a major QTL controlling tomato weight, fw2.2, also had changes in its promoter sequence, which influenced fruit weight (Nesbitt and Tanksley. 2002). In B. rapa, Su et al (2021) identified sequence variations in the promoter of BrHISN2, which conferred cold-dependent expression on BrHISN2, resulting in leaf yellowing. Thus, the significantly different expression levels of BraA07g018240.3C between the two parental lines might be caused by changes in the promoter region.
In conclusion, we fine mapped a novel QTL for FT, qFT7.1, to a physical interval of 56.4 kb on chromosome A07. The CSSL16 allele at the qFT7.1 locus regulates the FT negatively at the bolting stage of B. rapa. According to 1 3 the sequence and expression level, the most likely candidate gene, BraA07g018240.3C, encodes a TFL1-like protein.
These findings could aid our understanding of the mechanisms underlying flower formation and provide a genetic resource for B. rapa crop improvement studies.