Comparative Transcriptome Analysis of Early- and Late-Bolting Traits in Chinese Cabbage (Brassica rapa)

Chinese cabbage is one of the most important and widely consumed vegetables in China. The developmental transition from the vegetative to reproductive phase is a crucial process in the life cycle of flowering plants. In spring-sown Chinese cabbage, late bolting is desirable over early bolting. In this study, we analyzed double haploid (DH) lines of late bolting (“Y410-1” and “SY2004”) heading Chinese cabbage (Brassica rapa var. pekinensis) and early-bolting Chinese cabbage (“CX14-1”) (B. rapa ssp. chinensis var. parachinensis) by comparative transcriptome profiling using the Illumina RNA-seq platform. We assembled 721.49 million clean high-quality paired-end reads into 47,363 transcripts and 47,363 genes, including 3,144 novel unigenes. There were 12,932, 4,732, and 4,732 differentially expressed genes (DEGs) in pairwise comparisons of Y410-1 vs. CX14-1, SY2004 vs. CX14-1, and Y410-1 vs. SY2004, respectively. The RNA-seq results were confirmed by reverse transcription quantitative real-time PCR (RT-qPCR). A Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs revealed significant enrichment for plant hormone and signal transduction as well as starch and sucrose metabolism pathways. Among DEGs related to plant hormone and signal transduction, six unigenes encoding the indole-3-acetic acid-induced protein ARG7 (BraA02g009130), auxin-responsive protein SAUR41 (BraA09g058230), serine/threonine-protein kinase BSK11 (BraA07g032960), auxin-induced protein 15A (BraA10g019860), and abscisic acid receptor PYR1 (BraA08g012630 and BraA01g009450), were upregulated in both late bolting Chinese cabbage lines (Y410-1 and SY2004) and were identified as putative candidates for the trait. These results improve our understanding of the molecular mechanisms underlying flowering in Chinese cabbage and provide a foundation for studies of this key trait in related species.

Brassica rapa subspecies show high morphological diversity (Zhao et al., 2005). In particular, the heading Chinese cabbage (B. rapa L. subsp. pekinensis) is one of the most important Brassica vegetables in Asian countries, especially in China, Korea, and Japan (Bong et al., 2012;Cheng et al., 2014). It is a major vegetable crop produced in China, where it is considered as a key source of mineral nutrition (Wu et al., 2008). Moreover, Chinese cabbage is the main ingredient of the most popular traditional Korean side dish kimchi (Bong et al., 2012;Lee et al., 2014). The flowering Chinese cabbage (B. rapa ssp. chinensis var. parachinensis) produces elongated, tender, and thick stalks with rapid bolting as edible organs (Cheng et al., 2014;Huang et al., 2017). Pak choi (B. rapa subsp. chinensis) produces smooth dark green leaves with a prominent white midrib instead of forming a leafy head (Zhao et al., 2005).
In angiosperms, flowering is the most important developmental transition in the plant life cycle . This transition is controlled by endogenous and environmental signals Zhang et al., 2015). More than 180 genes identified by functional analyses are associated with flowering time in Arabidopsis (Fornara et al., 2010). Several of these genes form a complex regulatory network involving six key pathways, including photoperiod, vernalization, ambient temperature, age, autonomy, and gibberellin pathways (Fornara et al., 2010;Zhang et al., 2015;Huang et al., 2017). However, photoperiod and vernalization related to day length and low temperatures, respectively, have also been identified as the major pathways for the regulation of flowering time (reviewed in Song et al., 2013). Several homologs of Arabidopsis genes related to flowering in B. rapa and other plant species have been identified (Andersen et al., 2004;Greenup et al., 2009;Cheng et al., 2011;Duan et al., 2015;Xu et al., 2015). Mutations in the CONSTANS (CO), GIGANTEA (GI), and FLOWERING LOCUS T (FT), related to the photoperiod pathway, result in delayed flowering, but short days do not affect flowering time, unlike wild-type Arabidopsis (Suárez-López et al., 2001;Simpson and Dean, 2002;Moon et al., 2003). The CO gene encodes a zinc finger transcription factor that triggers the transcriptional upregulation of downstream floral integrator genes, including FT and SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1), in leaves under long-day conditions (Samach et al., 2000;Srikanth and Schmid, 2011;Zhang et al., 2015). Furthermore, FT and SOC1 act as activators of floral meristem identity genes, such as LEAFY (LFY) and APETALA 1 (AP1) (Amasino, 2005). In Arabidopsis, FRIGIDA (FRI; encoding two coiled−coil motif-containing protein) and FLOWERING LOCUS C (FLC; encoding a MADS−box transcription factor) are involved in the annual winter habit (late flowering) (Michaels and Amasino, 1999;Johanson et al., 2000;Moon et al., 2003;Amasino, 2005). Mutations in these genes result in early flowering in Arabidopsis (Michaels and Amasino, 1999;Johanson et al., 2000;Choi et al., 2011). FRI positively regulates the transcription of the flowering repressor, FLC (Choi et al., 2011). In contrast, vernalization results in early flowering via suppressing the FLC (Michaels and Amasino, 1999;Sheldon et al., 1999). Three genes, VIN3 (VERNALIZATION INSENSITIVE 3), VRN2 (VERNALIZATION 2), and VRN1 (VERNALIZATION 1), are required for the cold-mediated repression of FLC (Kobayashi et al., 1999;Gendall et al., 2001;Sung and Amasino, 2004). Su et al. (2018) reported that BrVIN3.1 and BrFLC1 are the important genetic determinants of bolting time variation in B. rapa. Recent reports have revealed that there are various orthologs of FLC and FT for flowering time variation in B. rapa and B. oleracea (Schiessl et al., 2017). In B. rapa, there are two copies of FT on chromosomes A02 and A07 and four copies of FLC on chromosomes A02, A03, and A10 (reviewed in Schiessl et al., 2017). Moreover, Shu et al. (2018) reported a major quantitative trait locus (QTL) (Ef2.1) for early flowering in "broccoli × cabbage" and identified BolGRF6 as a putative candidate for early flowering in broccoli. BrSDG8 code for encoding a histone methyltransferase is associated with bolting in B. rapa ssp. pekinensis . Huang et al. (2020) found that histone methyltransferase CURLY LEAF (CLF; Bra032169) controls the expression of flowering-related genes, and mutation in the Bra032169 caused early bolting in Chinese cabbage.
Gibberellic acid (GA), a plant hormone, plays an important role in the regulation of flowering time (Bernier, 1988). In Arabidopsis, the application of exogenous GA promotes flowering under short-day conditions (Langridge, 1957;Chandler and Dean, 1994). Mutations in genes related to GA biosynthesis and signaling can alter flowering time in Arabidopsis (Wilson et al., 1992;Sun and Kamiya, 1994;Jacobsen et al., 1996). In Brassica napus, flower initiation, flowering time, and shoot elongation are regulated by endogenous GA (Dahanayake and Galwey, 1999). Exogenous GA can stimulate flower development, while an inhibitor of GA delays or inhibits flowering in Brassica (Rood et al., 1989;Zanewich et al., 1990).
In flowering Chinese cabbage, aerial vegetative parts and floral buds are consumed. The early-bolting flowering Chinese cabbage double haploid (DH) line "CX14-1" is characterized by rapid flower stalk development and early flowering. In contrast, late-bolting heading Chinese cabbage exhibits a long period of vegetative growth before flower bud initiation. In B. rapa, longday conditions together with vernalization promote reproductive growth over vegetative growth; consequently, plants bolt before reaching the harvesting stage, which is a serious problem in this crop . Therefore, late-bolting traits are desirable than early flowering for the cultivation of spring-sown Chinese cabbage.

Plant Materials
The late-bolting DH heading Chinese cabbage (B. rapa L. ssp. pekinensis) lines "Y410-1" and "SY2004" and a vernalizationindependent early-bolting flowering Chinese cabbage (B. rapa ssp. chinensis var. parachinensis) DH line "CX14-1" were used (Figure 1). The plants were grown in the experimental field at Yuanyang (113 • 97' E and 35 • 5' N), Henan Academy of Agricultural Sciences, China. The flower primordia of these lines were used for RNA sequencing. The samples were collected, immediately frozen in liquid nitrogen, and stored at -80 • C.
Total RNA Isolation, Library Construction, RNA Sequencing, and Assembly The frozen flower bud samples were ground into a powder in liquid nitrogen. Thereafter, total RNA was isolated from 100 mg of powder using the RNeasy Mini Kit (Qiagen, Valencia, CA, United States) according to the manufacturer's guidelines. The concentration, integrity, and purity of the RNA samples were determined using a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States) and Agilent 2100 Bio Analyzer (Agilent Technologies, Palo Alto, CA, United States). RNA samples with more than seven RNA integrity numbers (RINs) were used for RNA-seq library preparation with biological replications. Nine RNA-seq libraries were constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Beverly, MA, United States) following the manufacturer's recommendations. The clustering of the indexcoded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v4-cBot-HS (Illumina, San Diego, CA, United States) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on the Illumina high-throughput sequencingby-synthesis technology platform to generate paired-end reads. The adapter sequences and low-quality reads were removed. The clean reads were then mapped to the B. rapa (Chiifu-401) v. 3.0 reference genome [Brassica database (BRAD)] 1 using HISAT2 2 .

Functional Annotation and Novel Gene Discovery
The functions of assembled genes and novel gene discovery were performed by BLAST searches against the NR [National Center for Biotechnology Information (NCBI) non-redundant protein sequences 3 ], Swiss-Prot 4 , gene ontology (GO) 5 , Clusters of Orthologous Groups of proteins (COG) 6 , Pfam 7 , and Kyoto Encyclopedia of Genes and Genomes (KEGG) 8 databases. After prediction of the amino acid sequences of the new genes, HMMER (Eddy, 1998) was used for comparisons with the Pfam database to obtain annotation information.

Quantification of Gene Expression and Identification of Differentially Expressed Genes
Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). The FPKM value for each gene was quantified according to the length of the gene and read count mapped to this gene. Differentially expressed genes (DEGs) between early-bolting and late-bolting Chinese cabbage lines were identified using the R package DEGseq (Wang et al., 2010). The resulting p-values were adjusted using the Benjamini and Hochberg method to control the false discovery rate. Genes with an adjusted p-value < 0.01 found by DEseq were identified as DEGs.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis
Gene ontology enrichment analysis of the DEGs was implemented by the GOseq R packages based on the Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs. KEGG (Kanehisa et al., 2008) is a database resource for understanding highlevel functions and utilities of the biological system from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other highthroughput experimental technologies (see text foot note 8). We used KOBAS (Mao et al., 2005) software to test the statistical enrichment of differential expression genes in KEGG pathways. A GO enrichment analysis of DEGs was performed using the GOseq R package based on the Wallenius non-central hypergeometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs.

cDNA Synthesis and qPCR Validation
A total of 1 µg of RNA was used for cDNA synthesis using SuperScript III following the manufacturer's protocol (Invitrogen, Gaithersburg, MD, United States). Then, the FPKM values for nine randomly selected genes were validated by reverse transcription quantitative real-time PCR (RT-qPCR) using the LightCycler 480II (Roche, Mannheim, Germany). A total of 45 ng of cDNA was used as a template for RT-qPCR with gene-specific primers (Supplementary Table 1) using 2 × SyGreen Mix (qPCRBIO Lo-ROX) (PCR Biosystems, London, United Kingdom). Thermocycling conditions were 95 • C for 5 min, 45 cycles of 95 • C for 10 s, 60 • C for 10 s, and 72 • C for 15 s. At the end of the PCR cycles, the Ct values were analyzed using LightCycler 480II software (Roche). The efficiency of each gene-specific primer was determined using pooled cDNA samples. The expression of each gene was normalized using the

Transcriptome Sequencing of Early-and Late-Bolting Chinese Cabbage Floral Buds
We performed transcriptome analyses of two late-bolting heading Chinese cabbage (B. rapa L. ssp. pekinensis) DH lines (Y410-1 and SY2004) and one early-bolting flowering Chinese cabbage DH line (CX14-1) (B. rapa ssp. chinensis var. parachinensis) with three biological replications using Illumina sequencing technology. The raw reads generated by RNA-seq were deposited in the NCBI "Sequence Reads Archive" (SRA) under the accession number PRJNA605481. After removal of low-quality sequences, adapters, and ambiguous reads, a total of 721.49 million clean paired-end reads were obtained ( Table 1). The clean reads for each sample totaled 8.81 Gb, and the Q30 base percentage was 92.71% or greater ( Table 1). The clean reads were aligned with the B. rapa (Chiifu-401) v. 3.0 reference genome (BRAD, see text foot note 1), and the efficiency of the alignment ranged from 86.20 to 89.54%. Furthermore, variable splicing prediction, a gene structure analysis, and new gene discovery were performed based on the comparison. The clean reads were then assembled into 47,363 transcripts and 47,363 genes, of which 3,144 were predicted as novel genes (Supplementary Table 2).

Functional Annotation and Classification
The putative functions of assembled genes were annotated by searches against public databases, including GO, COG, KEGG Orthology (KOG), KEGG, eggNOG, PFAM, NR, and SWISS-PROT. Among them,47,214,40,497,35,083,32,944,and 24,390 genes were annotated to the NR, eggNOG, PFAM, SWISS-PROT, and KOG databases, respectively ( Table 2). In addition, 2,564 of 3,144 novel genes were functionally annotated (Supplementary Table 2). Furthermore, we performed GO, COG, and KEGG pathway analyses to illustrate the biological functions of the Chinese cabbage floral bud transcriptomes. A GO term enrichment analysis was performed to identify terms in three general categories, biological process (BP), molecular function (MF), and cellular component (CC)  Table 6-8 and Figure 6A). These results indicated that early and late bolting might be associated with DEGs in these functional subgroups. According to a COG functional annotation analysis, 14,411 genes were classified into 25 COG categories ( Figure 6B). The predominant COG categories represented in the Chinese cabbage floral bud transcriptomes were G (carbohydrate transport and metabolism), R (general function prediction only), and T (signal transduction mechanisms) ( Figure 7B). Several genes in these categories were differentially expressed, which might contribute to flowering time differences between the early-and late-bolting Chinese cabbage DH lines.

Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis of the Differentially Expressed Genes
Different gene products interact with each other to exert biological functions, and pathway annotation analyses are therefore useful to predict the functions of gene. KEGG is a resource for the systematic analysis of gene function and genomic information databases, providing information about genes and expression patterns from a whole network perspective (Kanehisa et al., 2004). We therefore analyzed the metabolic pathways of the DEGs between early-and late-bolting Chinese cabbage lines using the KEGG database. Total 2,430 and 2,417 DEGs were assigned to 124, 128, and 127 pathways for Y410-1 vs. SY2004, Y410-1 vs. CX14-1, and SY2004 vs. CX14-1, respectively (Supplementary Table 9). Moreover, these pathways involved 15,408 genes, which were different from the DEGs FIGURE 4 | MA plot of differential expression for two samples. Each point represents the differentially expressed genes (DEGs). The x axis indicates log 2 (FPKM), and the y axis indicates log 2 (FC). The red and green dots represent significantly up-and downregulated genes, respectively; while black dots represent non-differentiated genes. FPKM, fragments per kilobase of transcript per million fragments mapped; FC, fold change.
assigned to pathways, indicating that some genes contribute to more than one KEGG pathway ( Table 2 and Supplementary  Table 9). For example, the novel gene predicted in this study "Brassica_rapa_newGene_2848" is involved in homologous recombination, mismatch repair, nucleotide excision repair, and DNA replication. However, six pathways, including plantpathogen interaction, starch and sucrose metabolism, carbon metabolism, biosynthesis of amino acids, phenylpropanoid biosynthesis, and plant hormone signal transduction, contained over 100 DEGs between early-and late-bolting lines (SY410-1 vs. CX14-1 and SY2004 vs. CX14-1) ( Figure 6C)

Expression Patterns of Differentially Expressed Genes Related to Plant Hormone and Signal Transduction
The expression levels of most of the DEGs related to ABA receptors were downregulated, except PYL8 and PYL9, while most of the DEGs encoding ABA-insensitive five-like protein were upregulated in late-bolting Y410-1 compared with earlybolting flowering Chinese cabbage DH lines (Figure 8 and Supplementary Table 9). Similarly, the expression levels of most of the DEGs related to ABA receptors were downregulated, except for PYL8 and PYL9. In contrast, most of the DEGs related to ABA-insensitive 5-like proteins were upregulated in late-bolting compared with early-bolting flowering Chinese cabbage DH lines (Figure 8 and Supplementary Table 10). The expression levels of DEGs encoding gibberellin receptors GID1A and GID1B were upregulated in both late-bolting (Y410-1 and SY2004) heading Chinese cabbage compared with the early-bolting flowering Chinese cabbage DH line (Figure 8 and Supplementary Table 10).
The expression levels of JA signaling genes, like jasmonateamino synthetase (JAR1, a GH3 family of protein), were downregulated in both late-bolting heading Chinese cabbage DH lines (Y410-1 and SY2004) compared with the earlybolting flowering Chinese cabbage DH line (Figure 8 and Supplementary Table 10). Moreover, the expression levels of DEGs encoding TIFY proteins were upregulated in the earlybolting flowering Chinese cabbage DH line (CX14-1) compared with both late-bolting heading Chinese cabbage lines (Figure 8  and Supplementary Table 10).
The unigenes related to auxin signal transduction, including auxin-responsive proteins, auxin transporter-like proteins, auxin responsive factor (ARF), and auxin-induced proteins, were differentially expressed between late-bolting heading Chinese cabbage and early-bolting flowering Chinese cabbage DH lines. The expression levels of 34 of these DEGs were upregulated in late-bolting heading Chinese cabbage Y410-1 compared with early-bolting flowering Chinese cabbage DH lines. Likewise, 27 and 33 DEGs were down-and upregulated, respectively, in another late-bolting DH line, SY2004 (Figure 8 and Supplementary Table 10). However, among these DEGs, 27 upregulated and 15 downregulated were common to both latebolting lines (Supplementary Table 10).

Validation of RNA-seq Data by RT-qPCR
We further tested the reliability of FPKM expression patterns (determined by RNA-seq) by RT-qPCR. Eight DEGs were randomly selected, and their relative expression levels were quantified using the same RNA samples extracted from earlyand late-bolting Chinese cabbage lines for RNA sequencing. The results confirmed that the expression patterns of the analyzed unigenes were consistent with the FPKM expression pattern obtained by RNA-seq (Figure 10).

DISCUSSION
Late bolting is an important economic trait in spring-sown lines, and late-bolting varieties have been developed and characterized for off-season production to meet the annual demand for Chinese cabbage (Yang et al., 2007;Su et al., 2018). Nevertheless, little is known about the genes and pathways associated with flowering time in flowering Chinese cabbage (B. rapa ssp. chinensis var. parachinensis) and late-bolting heading Chinese cabbage (B. rapa var. pekinensis). RNA sequencing technology has been effectively utilized for transcriptome analyses in a wide range of texa (Wang et al., 2009;Lu et al., 2010). This high-throughput sequencing technology can be used to determine global gene expression differences between populations or species with phenotypic differences and responses to environmental stress (Mortazavi et al., 2008;Wang et al., 2009;Miao and Luo, 2013). In the present study, we performed a comparative transcriptome analysis of early-bolting Chinese cabbage DH lines by RNA-seq. We obtained 721.49 million clean paired-end reads assembled into 47,363 unigenes, including 3,144 genes predicted as novel ( Table 1 and Supplementary Table 2). We detected 12,932 (Y410-1 vs. CX14-1) and 12,711 (SY2004 vs. CX14-1) DEGs between two sets of early-and late-bolting Chinese cabbage (Figure 2A and Supplementary Tables 4, 5).
A functional analysis of assembled unigenes revealed that the predominant COG categories were carbohydrate transport and metabolism, general function prediction only, and signal transduction mechanisms ( Figure 6B). Several unigenes in these categories were differentially expressed and might explain the difference in flowering times in early-bolting Chinese cabbage DH lines. Moreover, we identified 3,144 novel unigenes, of which 2,564 were functionally annotated (Supplementary Table 2).
A KEGG pathway enrichment analysis of DEGs between early-and late-bolting DH lines indicated that plant hormone and signal transduction (Ko04075) and starch and sucrose metabolism (Ko00500) were the most enriched pathways in both comparisons (Supplementary Table 9). These results suggest that unigenes related to plant hormones, signal transduction, and sucrose metabolism are involved in the regulation of flowering time in these Chinese cabbage DH lines.
Previous studies of mutants related to GA biosynthesis or signal transduction have revealed that GA can alter the flowering time (Wilson et al., 1992;Peng and Harberd, 1993;Sun and Kamiya, 1994;Jacobsen et al., 1996;Peng et al., 1997;Andres et al., 2014). Moreover, the DELLA domain protein RGA (repressor of ga1-3), GAI (GA insensitive), and RGA-like1 (RGL1 and RGL2) act as negative regulators of the GA signaling pathway (Yu et al., 2004;Silverstone et al., 2007). Our results revealed that two DEGs encoding DELLA proteins RGA1 (BraA06g040430) and RGA2 (BraA09g023210) are upregulated in late-bolting heading Chinese cabbage "Y410-1" (Figure 8 and Supplementary Table 10) and downregulated at the flowering stage compared with the vegetative stages in earlybolting flowering Chinese cabbage. These results suggest that the downregulation of negative regulators of the GA signaling pathway might increase GA levels, and the lack of function of these genes may trigger early bolting. Similar results have been reported by Huang et al. (2017), who reported that RGA1 and RGA2 are downregulated in flowering Chinese cabbage.
Several previous studies have indicated that the phytohormone jasmonic acid (JA) also regulates flowering time in Arabidopsis (Zhai et al., 2015). The F-box protein COI1 (coronatine insensitive 1) degrades JAZ (contains TIFY and Jas domains) repressors (Hoo and Howe, 2009). Transcript levels of NaJAZd and NaJAZh are upregulated in the early floral stages of NaJAZi-silenced plants due to the high JA content in the flowers (Li et al., 2017). Our results also showed that DEGs encoding TIFY proteins are more highly expressed in the early-bolting flowering Chinese cabbage DH line (CX14-1) than in both late-bolting heading Chinese cabbage lines (Figure 8 and Supplementary Table 10). These results suggest that the elevated expression of these genes in early flowering Chinese cabbage might be related to a high JA content.
Another steroidal phytohormone, brassinosteroid (BR), promotes flower induction in plants (Matsoukas, 2014). The BR signaling genes BZR1 (brassinazole-resistant1) and BES1 suppress the expression of important genes related to BR biosynthesis (CPD, constitutive photomorphogenesis, and dwarfism; DWF4, and DWARF4) by binding to their promoter regions in Arabidopsis (Wei et al., 2017). Moreover, BR-deficient/BRinsensitive mutants show delayed flowering time (Domagalska et al., 2007;Zhu et al., 2013). We also found that DEGs related to BZR1, BZR2, and BRI1 (brassinosteroid insensitive 1) were more highly expressed in both late-bolting heading Chinese cabbage DH lines than in early-bolting flowering Chinese cabbage (Figure 8 and Supplementary Table 10). These results suggest that the upregulation of BR signaling genes might affect BR biosynthesis and cause late bolting in late-bolting heading Chinese cabbage DH lines.
Starch and sucrose play important roles in flowering (Turnbull, 2011;Cho et al., 2018). Previous reports revealed that trehalose-6-phosphate acts as a signal molecule for flowering initiation in different plant species, including Arabidopsis thaliana (Sheen, 2014), grape (Caspari et al., 1998), and citrus (Shalom et al., 2014). Besides, Micallef et al. (1995) found that an increased level of endogenous sucrose promotes flowering in tomato. Wahl et al. (2013) reported that TREHALOSE-6-PHOSPHATE SYNTHASE 1 (TPS1) is required for the regulation of flowering time in A. thaliana. In the present study, alpha,alphatrehalose-phosphate synthase 5 (BraA03g047730) was upregulated in both late-bolting heading Chinese cabbage lines compared with early-bolting flowering Chinese cabbage lines (Figure 9). Razi et al. (2008) demonstrated that BoFLC alleles segregate independently from flowering time alleles in Brassica oleracea. However, Yuan et al. (2009) reported that variation in BrFLC1 is linked to flowering time in B. rapa. Moreover, Xie et al. (2015) performed a QTL analysis and showed that Br2 is the key determinant of BrFLC2 and a candidate flowering time locus in B. rapa. A transposon insertion in the coding sequence of BrFT2 located on a QTL on chromosome A07 (region Br5) causes late flowering . Therefore, we further analyzed the flowering time-related genes, especially FLC and FT, based on FPKM expression values. Unigenes, such as FLK (BraA03g031700) and FLD (BraA03g034300), showed higher expression levels in both late-bolting lines than in earlybolting lines, while FLT (BraA02g016700) showed the opposite pattern (Supplementary Table 11). With regard to FT genes, the expression levels of three unigenes encoding proteins related to flowering time control, such as FCA (BraA01g020520), FY (BraA02g004930), and FPA (BraA09g036880), were higher in both late-bolting lines than in early-bolting lines ( Supplementary  Table 11). Nonetheless, none of these unigenes were differentially expressed between early-and late-bolting lines.

DATA AVAILABILITY STATEMENT
Data is available at NCBI SRA accession: PRJNA605481.

AUTHOR CONTRIBUTIONS
XZ and YY conceived and designed the experiments. YZ, SY, and ZW performed the experiments. HS, LL, LN, and MH-U-R prepared the figures and tables. XW and MR drafted the work or revised it critically for important content. All authors contributed to the article and approved the submitted version.