Genotypic and phenotypic characterization of RcBr and CSSL16
To detect the segment which had introgressed from 08A061 on ten chromosomes of CSSL16, the two parental lines, RcBr and CSSL16, were genotyped using whole‑genome re-sequencing. A total of 67,507,156 and 64,540,134 high quality reads were detected in RcBr and CSSL16, and the clean data was compared to reference genome. This identified 2,999,421 SNPs and 627,670 InDels on the ten chromosomes between the two parental lines. According to the ED value calculation, the variation was mainly distributed in chromosome 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 B. rapa whole genome size is about 353.14 Mb (B. rapa V3.0), and the background recovery rate was about 99.42% (351.10/353.14).
RcBr and CSSL16 showed a significant difference in FT under multiple environments (Fig.2a). Under E1 growth conditions, the FT of RcBr (75.47 ± 3.78) was earlier than CSSL16 (98.19 ± 2.07) by about 23 days, under E2 conditions, the FT of RcBr (52.26 ± 1.27) was earlier than CSSL16 (63.43 ± 1.28) by about 11 days, under E3 conditions, the FT of RcBr (44.21 ± 3.55) was earlier than CSSL16 (61.69 ± 3.88) by about 17 days, and under E4 conditions, the FT of RcBr (38.78 ± 2.78) was earlier than CSSL16 (51.55 ± 3.37) by about 13 days (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 difference 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 ED5 values were sorted, the top 1% of ED5 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 the QTL, qFT7.1, identified by BSR-seq analysis, we carried out conventional QTL analysis with 300 F2 individuals in March 2020. A total of 13 polymorphic markers (Table S2) were screened from the difference interval (donor segment of 08A061), which was 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 BSR-seq 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 F2 population consisting of 2200 individuals was 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 were divided into ten groups according to their genotypes. The mean value of the homozygous progeny (DE5, DE10 and FT) were compared with RcBr at 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 to a region upstream of markers InDel705 and InDel708. 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, and physical position of 15,539,588 bp to 15,595,959 bp on chromosome A07 (Fig. 4c).
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 the genes is shown in Supplementary Table 3. All the genes were compared with Arabidopsishomolog 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 homologue (ATC) gene which belongs to FLOWERING LOCUS T (FT)family and encodes a protein similar to TERMINAL FLOWER1 (TFL1), the overexpression of which leads to a similar phenotype as TFL1. The encoded protein from the identified gene might inhibit the expression of critical flowering genes LEAFY (LFY) and APETALA1 (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 results 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). The expression of the others eight genes was also detected in the root, stem, leaf, and flower, with BraA07g018270.3C and BraA07g018300.3C showing significantly different expression levels in flower (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, RT-22, RT-36, RT-97, RT-47, and RT-41 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 result in Arabidopsis premature flowering, therefore, 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
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 12th base in first exon and a base T to C change at the 32nd base in the second exon, however, mutations were both 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.3Cbetween two parental lines (Fig. S1).