Nitrogen Signaling Genes and SOC1 Determine the Flowering Time in a Reciprocal Negative Feedback Loop in Chinese Cabbage (Brassica rapa L.) Based on CRISPR/Cas9-Mediated Mutagenesis of Multiple BrSOC1 Homologs

Precise flowering timing is critical for the plant life cycle. Here, we examined the molecular mechanisms and regulatory network associated with flowering in Chinese cabbage (Brassica rapa L.) by comparative transcriptome profiling of two Chinese cabbage inbred lines, “4004” (early bolting) and “50” (late bolting). RNA-Seq and quantitative reverse transcription PCR (qPCR) analyses showed that two positive nitric oxide (NO) signaling regulator genes, nitrite reductase (BrNIR) and nitrate reductase (BrNIA), were up-regulated in line “50” with or without vernalization. In agreement with the transcription analysis, the shoots in line “50” had substantially higher nitrogen levels than those in “4004”. Upon vernalization, the flowering repressor gene Circadian 1 (BrCIR1) was significantly up-regulated in line “50”, whereas the flowering enhancer genes named SUPPRESSOR OF OVEREXPRESSION OF CONSTANCE 1 homologs (BrSOC1s) were substantially up-regulated in line “4004”. CRISPR/Cas9-mediated mutagenesis in Chinese cabbage demonstrated that the BrSOC1-1/1-2/1-3 genes were involved in late flowering, and their expression was mutually exclusive with that of the nitrogen signaling genes. Thus, we identified two flowering mechanisms in Chinese cabbage: a reciprocal negative feedback loop between nitrogen signaling genes (BrNIA1 and BrNIR1) and BrSOC1s to control flowering time and positive feedback control of the expression of BrSOC1s.


Introduction
Flowering is a central event in the plant life cycle. Flowering time plasticity has evolved in relation to seasonal and developmental changes to maximize reproductive success in various environments [1]. It has emerged as one of the key traits affecting crop yield in commercial agronomic and horticultural crops because biotic and abiotic stresses can be avoided by modifying the flowering time [2,3]. Flowering is generally influenced by various hormones and nutrients. Nitric oxide (NO) is a central growth regulator known to inhibit flowering. In Arabidopsis, NO represses the expression of flowering activator genes CONSTANS (CO) and GIGANTEA (GI), which function in the photoperiod pathway, and elite variety "20" (2n = 20). CRISPR/Cas9-mediated multi-gene mutagenesis of BrSOC1s gave new insights into the flowering mechanism and demonstrated that genome editing could be applied to improve crop yield in Chinese cabbage.

Results
2.1. Differences in Vernalization-Induced Bolting Time between Lines "4004" and "50" Our previous work showed that bolting in the inbred line "4004" is induced by vernalization [12]. In this study, we compared the bolting time of two Chinese cabbage inbred lines, "4004" (early bolting) and "50" (late bolting). Seven-day-old seedlings of both lines were vernalized for 0 or 35 days, and plants were grown under long-day conditions for 60 days in soil before assessing the bolting phenotype. Under vernalization conditions, the inbred line "4004" bolted earlier than line "50" (Figure 1a). Consistent with the flowering phenotypes of the two lines, the average time to bolting was 17 days in "4004", which was approximately 13 days shorter than in line "50", exhibiting an average time to bolting of 30 days. Additionally, there were fewer rosette leaves in line "4004" than in line "50" (average of the number of leaves at bolting: 9.5 in line "4004" vs. 14.5 in line "50") ( Figure 1b). However, in the absence of vernalization, both inbred lines showed no bolting for 12 weeks (Supplementary Figure S1). Thus, our results indicated that vernalization has a significant effect on bolting in Chinese cabbage, and plants of inbred lines "4004" and "50" grown from vernalized seedlings displayed distinct phenotypes. mechanism of flowering in Chinese cabbage by performing RNA-Seq analysis of two different bolting time inbred lines ("4004" and "50" both are 2n = 20) before and after exposure to vernalization conditions. Further, we used the CRISPR system for targeted mutagenesis of BrSOC1 homologs (BrSOC1s) to introduce the late-flowering trait into an elite variety "20" (2n = 20). CRISPR/Cas9-mediated multi-gene mutagenesis of BrSOC1s gave new insights into the flowering mechanism and demonstrated that genome editing could be applied to improve crop yield in Chinese cabbage.

Differences in Vernalization-Induced Bolting Time between Lines "4004" and "50"
Our previous work showed that bolting in the inbred line "4004" is induced by vernalization [12]. In this study, we compared the bolting time of two Chinese cabbage inbred lines, "4004" (early bolting) and "50" (late bolting). Seven-day-old seedlings of both lines were vernalized for 0 or 35 days, and plants were grown under long-day conditions for 60 days in soil before assessing the bolting phenotype. Under vernalization conditions, the inbred line "4004" bolted earlier than line "50" (Figure 1a). Consistent with the flowering phenotypes of the two lines, the average time to bolting was 17 days in "4004", which was approximately 13 days shorter than in line "50", exhibiting an average time to bolting of 30 days. Additionally, there were fewer rosette leaves in line "4004" than in line "50" (average of the number of leaves at bolting: 9.5 in line "4004" vs. 14.5 in line "50") ( Figure 1b). However, in the absence of vernalization, both inbred lines showed no bolting for 12 weeks (Supplementary Figure S1). Thus, our results indicated that vernalization has a significant effect on bolting in Chinese cabbage, and plants of inbred lines "4004" and "50" grown from vernalized seedlings displayed distinct phenotypes.  Extensive bolting variation in two Chinese cabbage inbred lines, "4004" and "50", in response to vernalization. (a) Bolting phenotypes of the early-bolting line "4004" and the late-bolting line "50" following vernalization for 0 or 35 days. Non-vernalized seedlings (0 days, 0D Ver) were grown at 23 • C for 30 days. Vernalized seedlings (35 days, 35D Ver) were grown at 5 ± 1 • C and 12 h light/12 h dark photoperiod for 35 days and then transferred to a growth room at 23 • C for 20 days. Scale bar = 5 cm. (b) Statistical analysis of the days to bolting and the number of leaves in the two inbred lines. The leaves were counted when the plants started bolting. All experiments were performed in three independent biological replicates (n = 10 per replicate). Statistically significant differences are indicated with asterisks (* p < 0.05, *** p < 0.001; Student's t-test).

Transcriptome Profiling
To investigate the molecular mechanisms determining the bolting time in Chinese cabbage, we performed comparative RNA-Seq analyses of lines "4004" and "50" vernalized for 0 (0D Ver) and 35 days (35D Ver) (Supplementary Figure S2). RNA isolates were prepared from the shoots of both inbred lines, and cDNA libraries were sequenced on the Illumina HiSeq 2000 platform in triplicate (primer sets used in this study, Supplementary Table S1). The base quality of raw reads was evaluated using a series of pre-processing steps, and 56-76% of the raw reads were retained (Supplementary Table S2), indicating the use of relatively stringent pre-processing analysis criteria. Then, clean reads were mapped using the B. rapa reference genome (http://brassicadb.org/brad/, accessed on 27 April 2021, version 1.5) (Supplementary Table S3). Approximately 70-84% of the clean reads mapped to the coding sequences of genes, indicating that the transcriptome sequences were of good quality. Although pre-processing eliminated a large number of reads, the final mapped reads were suitable for further functional transcript analysis.

Identification of Differentially Expressed Genes (DEGs) and Functional Pathway Enrichment
Analysis Comparing Lines "4004" and "50" The DEGs were identified based on the following criteria: normalized read count ≥ 500, log 2 (fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05. The comparison between the two inbred lines ("4004" vs. "50") identified a total of 1710 DEGs, most of which (1229, 70%) were specific to the non-vernalized (0D Ver) conditions, whereas a less significant proportion of 10% (179) was detected under vernalization conditions (35D Ver) ( Figure 2a); 302 DEGs (17.5%) were found under both 0D Ver and 35D Ver conditions. This suggested that the heterogeneity between the two inbred lines based on the up-and down-regulation of DEGs was mostly inherent and not much affected by vernalization.

Identification of B. rapa Nitrate Reductase Genes and Relative Expression Levels Derived by
Comparing Lines "4004" and "50" Among the DEGs that were up-regulated in line "50" relative to their expression in line "4004", there were several genes that were highly associated with "Nitrogen metabolism" pathways. Therefore, we decided to examine DEGs involved in the NO signaling pathway. To study the functional relevance of nitrogen for flowering in Chinese cabbage, we applied the Arabidopsis research results to Brassicaceae plants. Previous studies had demonstrated that NIA and NIR were associated with flowering in Arabidopsis. Specifically, the nia1 nia2 double knockout mutant of Arabidopsis displayed an early-flowering phenotype and antagonistic gene expression compared with SOC1, indicating a negative regulation mechanism for flowering [5,37]. Based on this information, we created a schematic diagram showing the negative regulation of flowering by nitrogen ( Figure 3a) and assumed that the genes participating in N metabolism, NIA and NIR, may play an important role in controlling bolting time in lines "4004" and "50". We examined whether the DEGs BrNIA and BrNIR that we identified in Chinese cabbage are homologs of Arabidopsis genes by performing a phylogenetic tree analysis based on the predicted amino acid sequences using tblastx. Arabidopsis proteins AtNIA1 (At1g77760) and AtNIR1 (At2g15620) were highly homologous with the corresponding B. rapa proteins ( Figure 3b). Specifically, BrNIA and BrNIR genes were identified as DEGs by comparing the two lines with and without vernalization (Supplementary Table S4  In the qPCR analysis, the expression level of each from "4004" on day 0 was defined as "1". Error bars represent standard error (SE) of three replicates. Statistically significant differences are indicated with asterisks (* p < 0.05, *** p < 0.001; Student's t-test). (d) Nitrogen content in lines "4004" and "50" at each vernalization time point (0 and 35 days). Three black dots indicate each biological replicate (n = 3), bar graphs mark the average of biological triplicates, error bars mean ± SE of biological triplicates, and the different letters indicate a significant difference between samples (p < 0.05, one-way ANOVA followed by Tukey post-hoc test).

Identification and Analysis of Flowering Time Genes Differentially Expressed between the Two Inbred Lines
Next, we analyzed the DEGs enriched in the flowering pathway. The DEGs analysis (read count > 500, log2(fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05) for genes upregulated in "4004" vs. "50" detected only 12 flowering time-related DEGs (0D Ver: 6 DEGs, and 35D Ver: 6 DEGs) (Supplementary Table S5). Interestingly, the three BrSOC1s genes related to flowering time expression control were identified as up-regulated flowering time DEGs in line "4004" under 35D Ver conditions, and they also had the tendency that their expression levels were inversely correlated to those of BrNIA1 and BrNIR1 (Supplementary Table S4). However, it was difficult to broadly identify flowering (b) Phylogenetic relationship between Arabidopsis and Chinese cabbage nitrogen reductase-encoding genes. The phylogenetic tree was created using the neighbor-joining method. The scale bar shows evolutionary distances based on amino acid substitutions. (c) Quantitative RT-PCR (qPCR) results of nitrogen reductase genes (Bra015656 and Bra015227) compared to that of RNA-seq results. The actin (BrACT2) was used as an internal control. In the qPCR analysis, the expression level of each from "4004" on day 0 was defined as "1". Error bars represent standard error (SE) of three replicates. Statistically significant differences are indicated with asterisks (* p < 0.05, *** p < 0.001; Student's t-test). (d) Nitrogen content in lines "4004" and "50" at each vernalization time point (0 and 35 days). Three black dots indicate each biological replicate (n = 3), bar graphs mark the average of biological triplicates, error bars mean ± SE of biological triplicates, and the different letters indicate a significant difference between samples (p < 0.05, one-way ANOVA followed by Tukey post-hoc test).
To confirm the reliability of RNA-Seq data, the expression levels of two genes, BrNIA1 (Bra015656) and BrNIR1 (Bra015227) (shown in the phylogenetic tree), were analyzed by qPCR ( Figure 3c). The results of RNA-Seq and qPCR analyses were generally consistent. Overall, Bra015227 (BrNIR1) and Bra015656 (BrNIA1) had higher expression levels in line "50" than in line "4004" regardless of the vernalization status. Moreover, vernalization (35D Ver) increased the BrNIR1 expression level by a 2-fold change in line "50" relative to the level in line "4004". In general, the expression levels of NO signaling genes were well correlated with the bolting phenotypes of lines "4004" and "50". This suggested that the effect of nitrate signaling on the flowering pathway in Chinese cabbage was similar to that in Arabidopsis.
Then, we examined the nitrate content in lines "4004" and "50" based on the differences in the NO signaling gene expression as previously described [38]. The nitrogen content was substantially higher in line "50" than in line "4004" measured at the 0 and 35 days vernalization time points with a more than 10-fold difference (Figure 3d). This result was consistent with the previous finding that nitrogen starvation rapidly suppressed the NO signaling genes NIA (Os02g0770800, Os08g468100) and NIR (Os01g0357100) in rice [39]. Our findings supported the notion that BrNIA and BrNIR were directly involved in nitrate/nitrite assimilation in Chinese cabbage. Thus, we concluded that nitrogen signaling genes were related to the flowering pathway in Chinese cabbage as in the Arabidopsis system.

Identification and Analysis of Flowering Time Genes Differentially Expressed between the Two Inbred Lines
Next, we analyzed the DEGs enriched in the flowering pathway. The DEGs analysis (read count > 500, log 2 (fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05) for genes upregulated in "4004" vs. "50" detected only 12 flowering time-related DEGs (0D Ver: 6 DEGs, and 35D Ver: 6 DEGs) (Supplementary Table S5). Interestingly, the three BrSOC1s genes related to flowering time expression control were identified as up-regulated flowering time DEGs in line "4004" under 35D Ver conditions, and they also had the tendency that their expression levels were inversely correlated to those of BrNIA1 and BrNIR1 (Supplementary  Table S4). However, it was difficult to broadly identify flowering time-related DEGs based on just 12 DEGs. We reasoned that we could increase the number of flowering time DEGs by adjusting the expression levels via normalized read counts because flowering time genes were expressed at relatively low levels. We observed the pattern of increasing the number of flowering time-related DEGs while lowering the expression level conditions. There were 34 and 15 flowering time-related DEGs without and with vernalization, respectively, using the following criteria: read count > 50, log 2 (fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05 (Figure 4a). To assess the molecular function of flowering time DEGs with known roles in the flowering pathway, we examined the ratio of flowering enhancer and repressor genes expressed in both inbred lines with the read count ≥ 50 (Figure 4b). Under 0D Ver conditions, the numbers of up-regulated flowering enhancers vs. repressors among flowering time-related DEGs were similar (16 vs. 17, respectively) in the comparison of "4004" vs. "50" using up-regulating conditions (left graph). Under 35D Ver conditions, there were more up-regulated enhancer genes than repressor genes (9 vs. 2), and the repressor pattern also coincided with their bolting characteristics (1 vs. 3) (right graph). Furthermore, we identified the flowering time DEGs according to the distribution of flowering pathways. Under 0D Ver conditions, the Circadian/Light/Photoperiod (C/L/P) pathway accounted for more than 50% of up-regulated Ft DEGs in the comparison of "4004" vs. "50" (Figure 4c, left graph). In contrast, under 35D Ver conditions, the DEGs belonging to the C/L/P pathway were selectively decreased (Figure 4c, right graph). These results indicated that vernalization diminished the number of up-regulated flowering enhancers in the comparison of "4004" vs. "50" for flowering time DEGs. This applied especially to the flowering time DEGs of the C/L/P pathway, which were more affected by vernalization than those of other flowering pathways. pathway were selectively decreased (Figure 4c, right graph). These res vernalization diminished the number of up-regulated flowering comparison of "4004" vs. "50" for flowering time DEGs. This applie flowering time DEGs of the C/L/P pathway, which were more affecte than those of other flowering pathways.

Validation of the Expression of Flowering Time Genes Involved in the Flowering Regulatory Network
To validate the RNA-Seq data of DEGs involved in the flowering pathway, we conducted a qPCR analysis of flowering time genes exhibiting a 2-fold difference in expression levels between lines "4004" and "50" before and after vernalization (Supplementary  Table S5, * marked) using the same RNA samples that were used for the RNA-Seq analysis ( Figure 5). Most C/L/P pathway genes were induced by vernalization in both inbred lines, except BrTEM1 (a flowering repressor). The expression levels of C/L/P pathway flowering time genes, BrCIR1 and BrTEM1, were highly correlated between the RNA-Seq and qPCR analyses, although the differences for the BrCRY2 and BrPRR5 expression levels between the two inbred lines were smaller with the qPCR method than with the RNA-Seq method. Overall, the expression profile of BrCIR1 in the two inbred lines, but not that of BrTEM1, was consistent with the bolting phenotypes specific for the inbred lines. Next, we examined the expression of BrVRN1, which was the only flowering time DEG in the V pathway. The expression level of BrVRN1 varied substantially between lines "4004" and "50" before vernalization, which was consistent with their respective bolting phenotypes. Moreover, the results of qPCR and RNA-Seq analyses of BrVRN1 were similar, although there was no significant difference after vernalization. Interestingly, among the flowering time DEGs belonging to the gibberellin and Integrator pathways, the BrSOC1 homologs (BrSOC1-1, BrSOC1-2, and BrSOC1-3) reached remarkably high expression levels in line "4004" upon vernalization; however, in line "50", all BrSOC1 genes were barely expressed, regardless of vernalization. The BrGID1C gene was differentially expressed between the two inbred lines without vernalization but showed no significant difference in expression under vernalization conditions, which was consistent with the RNA-Seq data. The expression of BrGID1C in both inbred lines was also consistent with their respective bolting phenotypes, especially under non-vernalization conditions.  We also tested the expression levels of other Ft genes from both inbred lines (Supplementary Figure S4). The genes of a major flowering determinant, BrFLC1 and BrFLC3, were substantially repressed in both inbred lines after vernalization. Expression of most Ft genes (BrCCA1, BrCOL1-2, BrCOL9, BrVIN3A, BrGA1, BrTOC1, and BrLHY) was induced by vernalization regardless of the inbred line. The BrSPA3 and BrGI genes were not affected by vernalization. We hypothesized that the BrSOC1s, which displayed a distinct difference between "4004" and "50" after vernalization, might have a more important role in regulating the flowering of Chinese cabbage than other Ft genes. Based on RNA-seq and qPCR analyses, we identified BrSOC1s as the most important genes in vernalization. Then, we employed the CRISPR/Cas9 system to create BrSOC1s knockouts, which we used to examine whether these homologs were critical for the flowering phenotype in response to vernalization in Chinese cabbage. For this experiment, we choose the elite line "20" with a bolting time that was typical for BrSOC1s knockouts. In addition, we attempted to edit multiple BrSOC1 genes of (three homologs: BrSOC1-1/BrSOC1-2/BrSOC1-3) at once by locating the protospacer adjacent motif (PAM) sequence and designing single-guide RNA (SgRNA) for the 1st exon of the three BrSOC1s (Figure 6a). Then, we performed the BrSOC1s editing using the Chinese cabbage transformation and regeneration method in the inbred line "20" as previously reported [40]. After co-cultivation of Chinese cabbage and Agrobacterium tumefaciens carrying the CRISPR/Cas9 vector system (No. 1 in Figure 6b), we confirmed the shoot and root regeneration in a selection medium (No. 2 and 3 in Figure 6b). Finally, we generated BrSOC1s-edited T 0 plants (No. 4 in Figure 6b). We conducted next-generation sequencing (NGS) of T 0 plant samples, which confirmed BrSOC1s-specific gene editing (Supplementary Figure S5) based on the insertion/deletion (indel) proportion. The genotype proportion of BrSOC1-1 (Bra004928) with +1 bp insertion was 93.7% and that of BrSOC1-2 (Bra000393) with −1 bp deletion and +1 bp insertion combined was 87.4%. In contrast, the BrSOC1-3 (Bra039324) genotype proportion with +1 bp insertion was 34.0% (Figure 6c). In addition, "20" CRISPR/Cas9-brsoc1-1/1-2/1-3 T 1 (B36-1) plants did not display a bolting phenotype until 60 days after the 35 days of vernalization, unlike the "20" WT ( Figure 6d). The average time to bolting was approximately 17 days in vernalized WT line "20" plants, but no bolting phenotype was observed within 100 days in the CRISPR/Cas9-mediated BrSOC1s mutagenesis T 1 plants after the vernalization. In addition, the average number of leaves at the time of bolting was approximately 15 in WT line "20", whereas the leaves could not be enumerated in the genome-edited brsoc1s because of leaf senescence, which already occurred during the long growth period (>60 days after 35 days of vernalization) (Figure 6e). Finally, we found that the high proportion of gene-edited BrSOC1-1 and BrSOC1-2 was sufficient to introduce the late-bolting trait in Chinese cabbage.

Backcross Breeding for Brsoc1/2/3 Mutagenesis Controlling the Late-Bolting Trait
We harvested 20 seeds from the T 0 plant and performed NGS analysis, along with flowering phenotype observation in T 1 plants obtained from 20 seeds. We determined that T 1 plants, which were typically edited with BrSOC1-1 and BrSOC1-2 as in B36-1-9 (brsoc1-1/1-2) and BrSOC1-1, BrSOC1-2, and BrSOC1-3 as in B36-1-18 (brsoc1-1/1-2/1-3), displayed late bolting (Figure 7a). However, T 1 plants were difficult to use in the experiment because the Cas9 gene was integrated all plants. We then analyzed T 2 generations of the two T 1 plants. In addition, Chinese cabbage lines were selected that had the Cas9 gene removed (Supplementary Figure S6). line "20", whereas the leaves could not be enumerated in the genome-e because of leaf senescence, which already occurred during the long growth days after 35 days of vernalization) (Figure 6e). Finally, we found that the hig of gene-edited BrSOC1-1 and BrSOC1-2 was sufficient to introduce the late-b Chinese cabbage.  flowering phenotype observation in T1 plants obtained from 20 seeds. We determined that T1 plants, which were typically edited with BrSOC1-1 and BrSOC1-2 as in B36-1-9 (brsoc1-1/1-2) and BrSOC1-1, BrSOC1-2, and BrSOC1-3 as in B36-1-18 (brsoc1-1/1-2/1-3), displayed late bolting (Figure 7a). However, T1 plants were difficult to use in the experiment because the Cas9 gene was integrated all plants. We then analyzed T2 generations of the two T1 plants. In addition, Chinese cabbage lines were selected that had the Cas9 gene removed (Supplementary Figure S6). Three black dots indicate each biological replicates (n = 3), the graph bars represent the averages of biological triplicates, error bars mean ± SE of biological triplicates, and different letters indicate a significant difference between samples (p < 0.05, one-way ANOVA followed by Tukey post-hoc test).
To validate that the CRISPR/Cas9-mediated mutagenesis of brsoc1s was linked to the late-bolting trait, the B36-1-9 line was backcrossed (B36-1-9 BC). The indel frequency of BrSOC1s was restored to about 50% in B36-1-9 BC T1 plants compared to before being B36-1-9 based on the NGS analysis (Supplementary Table S6). The B36-1-9 BC T1 plants were BrACT2 was used as an internal control to quantify the transcripts. The expression level of a gene in the "20" wild type (WT) on day 0 (no vernalization) was defined as "1". Values are mean ± SD. Statistically significant differences are indicated with asterisks (* p < 0.05, ** p < 0.01, *** p < 0.001; Student's t-test). (c) Nitrogen content in "20" WT and B36-1 T 2 #18 at both vernalization time points (0 and 35 days). Three black dots indicate each biological replicates (n = 3), the graph bars represent the averages of biological triplicates, error bars mean ± SE of biological triplicates, and different letters indicate a significant difference between samples (p < 0.05, one-way ANOVA followed by Tukey post-hoc test).

Loss-of-Function BrSOC1s Showed Reciprocal Regulation with Nitrogen Signaling Genes and Nitrogen Content and Preserved Positive Feedback Regulation of the Original BrSOC1s Genes
To assess the molecular characteristics of the genome-edited brsoc1s mutants, we performed a qPCR analysis of BrSOC1s and nitrogen signaling genes and determined the nitrogen content in T 2 plants based on the previous results comparing the two inbred lines, "4004" and "50". Firstly, we determined the expression levels of original BrSOC1 genes based on the observation that the positive feedback regulation is closely connected to flowering time in Arabidopsis [41]. The WT and B36-1 plants displayed different expression levels and patterns of the BrSOC1s. Under 0D Ver, when the SOC1 expression remained at a basal level, there was no significant difference in the BrSOC1s expression levels of the genome-edited plants compared to those in the WT, although BrSOC1-2 was decreased in the B36-1s. Upon vernalization, the expression levels of all three BrSOC1s were markedly induced in the WT line "20" with up to 30-fold change, whereas the B36-1 genome-edited plants did not respond to vernalization, displaying significantly lower expression levels (Figure 7b, upper graphs). The loss-of-function NIA1 in Arabidopsis had antagonistically increased expression levels of SOC1. Therefore, we examined the nitrogen signaling genes BrNIA1 and BrNIR1 in B36-1s plants. Remarkably, the expression levels of these genes differed from those in the WT. Without vernalization, the expression level of BrNIR1 was lower in the B36-1s mutants than in the WT, but there was an approximately 4-fold (B36-1#9) or 10-fold (B36-1#18) expression level increase in response to vernalization. Similarly, the expression level of another nitrogen signaling gene, BrNIA1, was up to 40-fold increased after vernalization in B36-1s mutants, compared to the expression level without vernalization (Figure 7b, lower graphs). The significant expression level variations in the nitrogen signaling genes in B36-1s plants led us to the analysis of the nitrogen content. Interestingly, the nitrate content was higher in the B36-1#18 mutant than in the WT (4-fold higher in the 0D Ver sample and 9-fold in the 35D sample). The nitrate content in WT "20" plants was reduced in response to vernalization with a 2-fold change, whereas the B36-1#18 plants displayed no significant variation (Figure 7c). In addition, transcript levels of BrCIR1 and BrCRY2 were also decreased in the B36-1s than that of WT "20" under 35D Ver conditions (Supplemental Figure S8). Although the nitrate content did not completely match the tendency of the nitrogen reductase expression level observed by comparing WT and B36-1s, the data suggested that BrSOC1s played a central role in flowering and reversely modulated the expression levels of nitrogen signaling genes in Chinese cabbage. Additionally, this result indicated a BrSOC1s-specific control mechanism by positive feedback, which appeared to be conserved between Chinese cabbage and Arabidopsis.

Discussion
Although the knowledge of the flowering mechanism is of fundamental importance for improving crop productivity, the molecular basis of flowering time control has not yet been elucidated in Chinese cabbage. In this study, we accessed the molecular variation in flowering time between early-and late-bolting inbred lines of Chinese cabbage via RNA-Seq and CRISPR/Cas9 technologies.
We previously detected 223 flowering time genes in the early-bolting Chinese cabbage inbred line "4004", of which approximately 20% were regulated by vernalization [12]. In the present study, we performed transcriptome profiling of two inbred lines, "4004" and "50", exhibiting different bolting times. However, comparative transcriptome profiling of these two lines was difficult because of two reasons. Firstly, most of the DEGs identified between the two lines were detected in the absence of vernalization, and only a few DEGs were identified after vernalization ( Figure 2). Thus, the DEGs identified between the two lines before vernalization were responsive to vernalization. For instance, the expression of BrFLCs was highly down-regulated in both lines after vernalization, despite the significant difference between the two lines prior to vernalization, and the expression of BrTEM1, BrGID1C, and BrVRN1 was also substantially decreased after vernalization in line "4004" but not in line "50" (Supplementary Table S5). Secondly, because of the low read counts in our analysis, there was a substantial number of flowering time genes that failed to meet the criteria to be categorized as DEGs in response to the two vernalization conditions or based on the differences between the two inbred lines. Finally, only 5% of the flowering time genes (12 out of 223) were differentially expressed between line "4004" and line "50".
This suggested the possibility that other regulatory pathway genes controlled the bolting time in these two inbred lines.
The identification of DEGs involved in N metabolism is an interesting finding because nitrate or N in plant tissues affects the flowering time [42,43]. Therefore, BrNIR and BrNIA are good candidates for determining the flowering mechanism in Chinese cabbage. The RNA-Seq and qPCR analyses revealed that the expression levels of BrNIR and BrNIA were up-regulated in line "50" grown with or without vernalization ( Figure 3C). NIR is the main enzyme responsible for the production of NO and affects a wide range of physiological and developmental processes in plants, including flowering timing in Arabidopsis [5,44,45]. In Arabidopsis, flowering is facilitated when plant growth occurs under low-nitrate conditions compared with growth under high-nitrate conditions. Furthermore, the expression of flowering time genes is regulated by nitrate availability; for example, low nitrate conditions repress flowering repressors and activate flowering enhancers [43,46]. Because nitrate availability affects the flowering pathway and the timing of the vegetativeto-reproductive phase transition, we surmise that increased expression of BrNIR and BrNIA in line "50" contributes to flowering repression in Chinese cabbage, consistent with the study in Arabidopsis [5]. The differentially expressed flowering time genes that affect the N metabolism may comprise a potential regulatory network linked to N signaling and the flowering pathway in Chinese cabbage, as suggested by the findings in Arabidopsis. Our data suggest that BrCIR1 and BrSOC1s are responsive to vernalization ( Figure 5). The MYB-related protein CIR1, also known as REV2, is involved in circadian regulation and acts as a negative regulator of flowering in Arabidopsis [47]. Furthermore, nitrate treatment significantly induces CIR1 expression [6]. This led us to speculate that CIR1 is a critical factor for the N-mediated regulation of flowering time. Surprisingly, BrCIR1 expression was higher in line "50" than in line "4004" after vernalization ( Figure 5). Expression of the floral integrator gene SOC1 also increases under low-nitrate conditions, leading to accelerated flowering [46]. The difference in the expression level of BrSOC1s between the two inbred lines was also correlated with the transcript levels of N metabolism genes BrNIR and BrNIA (Figures 3 and 5).
CRISPR/Cas9-mediated mutagenesis is a reliable system employed in functional studies on multiple-gene clusters of unknown functions and is becoming one of the most effective tools to create desirable phenotypes in crops relying on transformation technology. Therefore, it is the preferable tool for verifying that the regulation and function of BrSOC1s are associated with N signaling in the Chinese cabbage flowering pathway. There are three SOC1 genes in Chinese cabbage, but it is not known whether all of them are involved in flowering time control. However, their involvement in flowering is considered to be highly likely because their response to vernalization appears to be robust ( Figure 5), and the double and triple mutants have a delayed flowering phenotype ( Figure 6). Moreover, the knockout mutants of SOC1 that were backcrossed with WT line "20" appear to indicate that BrSOC1 and BrSOC2 are concurrently involved in controlling flowering time in Chinese cabbage (Supplementary Figure S7).
In the BrSOC1 triple mutant, the expression of the BrNIR1 and BrNIA1 genes was increased after vernalization (Figure 7b), suggesting the possibility that the SOC1 gene inhibited NO signaling during vernalization treatment, as previously shown in Arabidopsis [5]. It can be presumed that the NO signal and SOC1 form a regulation loop that inhibits each other (Figure 8). In addition, SOC1 expression was not increased in the brsoc1s triple mutant even after vernalization. This may occur if the brsoc1s mutant transcript is unstable, but it appears to be less likely because it was similar or slightly lower before the vernalization process. Therefore, it is possible that BrSOC1 expression is regulated by positive feedback after vernalization. These two regulation mechanisms-suppression of negative regulators (NO signal) and positive feedback of itself-can quickly induce the SOC1 gene expression and strongly increase it after vernalization. More detailed research will be needed to confirm our model. Thus, our results suggest that NO signaling genes play an essential role in the regulation of flowering time in Chinese cabbage by differentially regulating the expression of the key flowering enhancer SOC1 (Figure 8).
brsoc1s triple mutant even after vernalization. This may occur if the brsoc1s mutant transcript is unstable, but it appears to be less likely because it was similar or slightly lower before the vernalization process. Therefore, it is possible that BrSOC1 expression is regulated by positive feedback after vernalization. These two regulation mechanismssuppression of negative regulators (NO signal) and positive feedback of itself-can quickly induce the SOC1 gene expression and strongly increase it after vernalization. More detailed research will be needed to confirm our model. Thus, our results suggest that NO signaling genes play an essential role in the regulation of flowering time in Chinese cabbage by differentially regulating the expression of the key flowering enhancer SOC1 (Figure 8). In summary, comparative transcriptome profiling of two inbred lines identified genes with a potentially critical role in regulating the bolting time in Chinese cabbage. Additionally, our data suggest that genes involved in nitrate metabolism also function in the flowering pathway in Chinese cabbage. In summary, comparative transcriptome profiling of two inbred lines identified genes with a potentially critical role in regulating the bolting time in Chinese cabbage. Additionally, our data suggest that genes involved in nitrate metabolism also function in the flowering pathway in Chinese cabbage.

Plant Materials and Growth Conditions
Inbred lines of Chinese cabbage (Brassica rapa ssp. pekinensis), "4004" (early-bolting line), "50" (late-bolting line), and "20" (early-bolting line, elite variety) were used in this study. Seeds were obtained from NongHyup Seed in Korea (Gyeonggi-do, Anseong, Korea). Seeds were sown in sterilized soil and placed in a growth room maintained at 23 • C under long-day conditions (16 h light/8 h dark). After 1 week, vernalization of seeds was initiated in a cold room at 5 ± 1 • C with a 12 h light/12 h dark photoperiod for 35 days. Seeds not subjected to vernalization served as a control. After vernalization, plants were transferred to the growth room and grown for 40 days. The flowering phenotype was assessed based on the rosette leaf numbers and the days to bolting, which were recorded when the length of the floral axis was ≥0.5 cm. All experiments were performed in three independent biological replicates (n = 10 per replicate). To conduct RNA-seq analysis, shoot samples from five independent plants without vernalization (0 days) and with vernalization (35 days) were collected at the same time point in the light/dark cycle; samples were collected in three independent biological replicates (Supplementary Figure S1a).

Total RNA Isolation, Library Construction, and RNA-Seq Analysis
Total RNA was extracted from 250 mg shoot tissue of lines "4004" and "50" subjected to vernalization for 0 days (D0_4004 and D0_50) or 35 days (D35_4004 and D35_50) in three biological replicates using RNAiso Plus TRIzol (TaKaRa, Shiga, Japan). For line "4004", we integrated duplicate RNA-seq data that we previously generated; therefore, only one set of samples was newly analyzed in this line. The purified total RNA was analyzed with a NanoDrop 1000 spectrophotometer (Thermo Fisher, Waltham, MA, USA). Then, 6 µg RNA was used for cDNA library construction after passing an RNA quality evaluation. Subsequently, cDNA sequencing libraries were constructed using the TruSeq RNA Library Prep Kit (Illumina, Inc., San Diego, CA, USA) and sequenced on the Illumina HiSeq™2000 platform (Illumina, Inc., San Diego, CA, USA) to generate 101 bp paired-end reads. All raw sequences were deposited in the GenBank Gene Expression Omnibus (GEO) under the accession numbers GSE106444 (for previous "4004" RNA-seq duplicates) and GSE139375.

Transcriptome Assembly and Mapping
Raw sequence reads were analyzed using the RNA-seq parameter of the Illumina pipeline. The raw reads were filtered using several trimming steps to remove adaptor contamination, low-quality sequences, and reads containing stretches of Ns. The DynamicTrim and LengthSort programs of SolexaQA (v. 1.13) were used on reads with a Phred quality score of 31 (Q ≥ 20) and a minimum length of 25 bp [48] (Supplementary Figure S1b).

Identification and Analysis of DEGs
Gene expression data of lines "4004" and "50" obtained in this study (a total of eight samples) and of four samples of line "4004" obtained previously [12] (a total of 12 samples) were used for the identification of DEGs. To identify DEGs between two vernalization treatments and between the two inbred lines, raw reads were normalized using the DESeq library in R (v3.5.1) [50]. Genes with normalized read count ≥ 500, log 2 (fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05 were designated as differentially expressed. Furthermore, flowering time-related DEGs were analyzed with read count > 50 for the purpose of gathering the number of DEGs (with the same criteria: log 2 (fold-change) ≥ 1, and adjusted p-value (FDR) < 0.05). The expression of DEGs in lines "4004" and "50" (4004 vs. 50) was normalized relative to their expression in line "50" treated with vernalization for 0 days (D0_50) to identify up-and down-regulated genes, which were depicted as Venn diagrams drawn using Microsoft PowerPoint and as bar graphs generated using Microsoft Excel. Flowering enhancer/repressor and pathway were analyzed as described previously [12].

Real-Time Reverse Transcription PCR (qPCR) Analysis
To validate the results of RNA-Seq analysis, selected genes were analyzed by qPCR in three biological replicates. Total RNA was treated with DNaseI (Thermo Fisher Scientific, Waltham, MA, USA) to remove genomic DNA contamination, and 2 µg total RNA was used for first-strand cDNA synthesis using PrimeScript TM RT Master Mix (TaKaRa, Shiga, Japan), according to the manufacturer's instructions. The cDNA was resuspended in nuclease-free water and used for qPCR analysis with a CFX Connect™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). BrActin (Bra022356) was used as an internal control for normalizing mRNA levels. All primer sets used in this study are listed in Supplementary  Table S1.

Bioinformatics Analysis
Amino acid sequences of AtNIA1 (At1g77760) and AtNIR1 (At2g15620) proteins were used for phylogenetic tree analysis. BrNIA1s and BrNIRs amino acid sequences were constructed using tblastx from transcript id of RNA sequencing data. The phylogenetic tree of NIA1 and NIR1 between Arabidopsis and Chinese cabbage was constructed using the neighbor-joining method in Molecular Evolutionary Genetic Analysis (MEGA; version 7) for whole protein sequences.

Nitrate Assay
The nitrate assay protocol was based on the salicylic acid method [51,52]. The assays were performed using three independent B. rapa leaves in the same position as previously described in an Arabidopsis methods [38]. Weighed leave samples (up to 100 mg) in a 1.5 mL tube were frozen in liquid nitrogen and milled into powder using a grinder. Then, 1 mL of sterilized deionized water was added into the tube, followed by boiling for 25 min. The samples were centrifuged at 13,000 rpm for 10 min at 4 • C, and 0.1 mL supernatant was transferred into a 15 mL conical tube (Hyundai Micro, Seoul, Korea). Next, 0.4 mL of 5% salicylic acid-sulfate acid solution (5 g salicylic acid (Sigma-Aldrich, St. Louis, MO, USA) in 100 mL sulfate acid (Junsei, Tokyo, Japan) was added, the sample was mixed thoroughly by inverting, the reaction incubation was conducted at 23 • C for 20 min and 9.5 mL of 8% NaOH solution was added. Before measurement, standard curves were checked using KNO 3 − (Sigma-Aldrich, St. Louis, MO, USA; 10 to 120 mg/L concentration), and regression analysis was conducted based on the standard curve. The OD 410 value was determined after cooling the tube to room temperature; 0.1 mL of sterilized deionized water was used as a control instead of the supernatant.

CRISPR/Cas9-Mediated Mutagenesis of BrSOC1-1/1-2/1-3 and Genetic Transformation
Cas-Designer (http://www.rgenome.net/cas-designer/, accessed on 27 April 2021) was used to design sgRNA. The sequences of the Chinese cabbage SOC1 genes (Bra000393, Bra004928, and Bra039324) were compared to find a region with matching sequences. Then, sgRNAs were designed (5 -TCTGATCATCTTCTCTCCTAAG-3 ) in the corresponding region to target the three BrSOC1 genes with one sgRNA. The sgRNA was synthesized with a restriction enzyme sequence for cloning, and it was cloned into the pHAtC vector [53] treated with the restriction enzyme AarI (CACCTGC (4/8)ˆ). The recombined vector was transformed into Agrobacterium tumefaciens LBA4404 and then used for transformation of Chinese cabbage plants.
The Chinese cabbage elite line "20" (NongHyup Seed. Gyeonggi-do, Anseong Korea) was used for the transformation of Chinese cabbage according to the previously published method [40]. For transformation, hypocotyl was incubated for~1-2 days in darkness, and then it was cut into segments with a length of 0.5 to 1 cm. Co-culture was performed with the transformed Agrobacterium tumefaciens in the dark for 2 days. After washing, the explants were cultivated in a callus induction medium (MS salts, including 3% sucrose, 5 mg/L benzyl adenine (BA), 1 mg/L naphthaleneacetic acid (NAA) and 300 mg/L cefotaxim) in the dark for 3 days. The induced calli were transferred to a shoot induction medium (MS salt, including 3% sucrose, 10 mg/L BA, 1 mg/L cefotaxime, and 10 mg/L hygromycin). Once the shoot was obtained, it was cultured in a root-inducing medium (MS salt, including 3% sucrose, 0.1 mg/L NAA, and 0.1 mg/L gibberellin).

Backcross Breeding for Brsoc1s
All organs, except the pistil, were removed from the "20" buds that were not opened; then, the Brsoc1/2/3 (B36-1-9, T 1 ) pollen was used for crossing. The obtained seeds were sown, and genomic DNA was extracted from each plant. The indel analysis was performed using NGS to confirm that the Brsoc1/2/3 mutation was hetero. Several confirmed plants were selected, and selfing was performed to obtain seeds of the next generation.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ijms22094631/s1, Figure S1: Bolting phenotypes of two Chinese cabbage inbred lines without vernalization; Figure S2: Experimental design for RNA-seq analysis; Figure S3: Identification of flowering time (Ft) genes based on gene ontology (GO) terms; Figure S4: Results of qPCR and RNAseq analysis of flowering time (Ft)-related genes in the inbred lines "4004" and "50" grown with or without vernalization; Figure S5: Multiple sequence alignment of insertion/deletion (indel) variants of three BrSOC1 homologs; Figure S6: Cas9 and Hygromycin (Hyg.) PCR analysis in B36-1 T2 plants; Figure S7: Backcrossing (BC) of brsoc1s gene-edited Chinese cabbage; Figure S8: Quantitative PCR analysis in BrSOC1s-edited Chinese cabbage plants grown with or without vernalization; Table S1: Primers used for RT-qPCR analysis; Table S2: Cleaned short reads statistics after trimming of RNA sequencing data; Table S3: Statistics of reads mapping to B. rapa coding sequence database; Table  S4: BrNIA1 a nd BrNIR1 groups sequencing data;  Data Availability Statement: All raw sequences were deposited in the GenBank Gene Expression Omnibus (GEO) under the accession numbers GSE106444 (for previous "4004" RNA-seq duplicates) and GSE139375.