Skip to main content

Gene co-expression network analysis reveals key pathways and hub genes in Chinese cabbage (Brassica rapa L.) during vernalization

Abstract

Background

Vernalization is a type of low temperature stress used to promote rapid bolting and flowering in plants. Although rapid bolting and flowering promote the reproduction of Chinese cabbages (Brassica rapa L. ssp. pekinensis), this process causes their commercial value to decline. Clarifying the mechanisms of vernalization is essential for its further application. We performed RNA sequencing of gradient-vernalization in order to explore the reasons for the different bolting process of two Chinese cabbage accessions during vernalization.

Results

There was considerable variation in gene expression between different-bolting Chinese cabbage accessions during vernalization. Comparative transcriptome analysis and weighted gene co-expression network analysis (WGCNA) were performed for different-bolting Chinese cabbage during different vernalization periods. The biological function analysis and hub gene annotation of highly relevant modules revealed that shoot system morphogenesis and polysaccharide and sugar metabolism caused early-bolting ‘XBJ’ to bolt and flower faster; chitin, ABA and ethylene-activated signaling pathways were enriched in late-bolting ‘JWW’; and leaf senescence and carbohydrate metabolism enrichment were found in the two Chinese cabbage-related modules, indicating that these pathways may be related to bolting and flowering. The high connectivity of hub genes regulated vernalization, including MTHFR2, CPRD49, AAP8, endoglucanase 10, BXLs, GATLs, and WRKYs. Additionally, five genes related to flower development, BBX32 (binds to the FT promoter), SUS1 (increases FT expression), TSF (the closest homologue of FT), PAO and NAC029 (plays a role in leaf senescence), were expressed in the two Chinese cabbage accessions.

Conclusion

The present work provides a comprehensive overview of vernalization-related gene networks in two different-bolting Chinese cabbages during vernalization. In addition, the candidate pathways and hub genes related to vernalization identified here will serve as a reference for breeders in the regulation of Chinese cabbage production.

Background

Chinese cabbage (Brassica rapa L. ssp. pekinensis), also known as heading cabbage or wrapping cabbage, is a leafy Brassica vegetable of the cruciferous family that originated in China with a long history of cultivation. Chinese cabbage has the characteristics of a rich variety of types, wide distribution, high yield, durability during storage and transportation, and a long supply period, and it is both highly nutritious and deeply loved by consumers. Chinese cabbage is one of the most economically important Brassica vegetable crops cultivated in Asian countries [1]. In Europe, especially Western Europe, the area of land under cultivation for Chinese cabbage has increased [2]. This indicates that the demand for Chinese cabbage throughout the year is slowly increasing. However, Chinese cabbage is susceptible to low temperatures (vernalization) and long daylight hours during the spring cultivation process, which causes it to bolt and flower quickly, thereby losing its commercial value. In contrast, in the breeding process, low temperature (vernalization) can be used to rapidly breed excellent varieties.

The transition from vegetative to reproductive growth is an important developmental step in the plant life cycle [3], and the timing of this switch is crucial for successful reproduction [4]. Vernalization, the effect of low temperature that induces and promotes flowering, is the main factor that promotes the transition from vegetative to reproductive growth in some biennial plants and annual winter plants. If plants that require low-temperature treatment do not undergo proper vernalization, flowering will be delayed by a few weeks or flower primordia will not form and will gradually decline. Different plants have different vernalization requirements depending on the developmental stage, vernalization temperature, and length of vernalization [5]. Previously, Yui and Yoshikawa [6] observed the phenomenon of low temperature promoting Chinese cabbage bolting and flowering. In the vernalization pathway, FLOWERING LOCUS C (FLC) is a key gene that controls flowering time. Many upstream genes ultimately determine bolting and flowering time by regulating the expression of FLC. FLC encodes a MADS-box transcription factor, which is a flowering inhibitor. The difference between early and late flowering depends largely on FLC allele variation [7]. FRIGIDA (FRI) is required for high FLC expression levels in Chinese cabbage and is a positive regulator of FLC [8]. Vernalization inhibits the expression of FLC and promotes flowering, and the dominant FRI allele strengthens the inhibition of FLC [9, 10]. The vernalization of Chinese cabbage also involves the expression of VIN3, VRN2, and VRN1 [11]. Among them, VRN1 and VRN2 inhibit the expression of FLC and maintain the state of vernalization. Moreover, VRN1 and VRN2 do not recover after vernalization and maintain a continuous low expression state. VIN3 participates in inhibiting the expression of FLC in early vernalization under low temperature conditions. In Chinese cabbage, Li Z et al. cloned the homologous gene BrpFLC of FLC of Arabidopsis and proved that different degrees of vernalization can reduce the transcription level of BrpFLC in different bolting-resistant cabbage varieties [12]. So far, four FLC homologous genes (BrFLC1, BrFLC2, BrFLC3, and BrFLC5) have been found and verified in Chinese cabbage [13, 14]. Recently, BrFLC5 has been proven to be a weakly regulated gene for flowering regulation in Chinese cabbage [15]. After years of research, genes including FLC, VIN3, and the VRN family are currently the most thoroughly studied genes related to vernalization in Chinese cabbage.

The transcriptome is used to study gene transcription in plant cells and the regulation of transcription overall. The application of RNA sequencing technology (RNA-Seq) has been widely used in various biological fields to explore various aspects of the life sciences. RNA-Seq has been widely used to study the related genes of many plants, including the characteristics of Arabidopsis [16], rice [17] and cucumber [18]. In a study on the vernalization of Brassica-type vegetables, Sun et al. [19] conducted a transcriptome analysis on pak choi (Brassica rapa subsp. chinensis) samples at different developmental stages after vernalized and control treatments to investigate differentially expressed genes (DEGs), and they found that Bra014527, Bra024097, and Bra035940 exhibited obvious changes after vernalization. The homologous genes of these three genes also participated in the vernalization response of Arabidopsis. Therefore, it was speculated that these genes also responded to vernalization in pak choi. Qi et al. [20] used an RNA-Seq technology to obtain information including the DEGs, functional annotations, and variable shear, of Chinese cabbage samples before and after vernalization. Four candidate genes related to flowering were screened. As an important flowering crop, it is necessary to explore the underlying molecular mechanisms of flowering induction in Chinese cabbage.

Currently, vernalization is widely applied in vegetable production, especially in leafy vegetables. Spring Chinese cabbage lose their commercial value after premature bolting as a result of low-temperature effects. The length of breeding time is also shortened due to rapid bolting and flowering caused by vernalization. Therefore, the effects of vernalization on Chinese cabbage are worth dissecting and exploring. In this study, the gradient vernalization of two different bolting Chinese cabbage accessions were used to analyze the transcriptome pattern of Chinese cabbage during vernalization. Using a weighted gene co-expression network analysis (WGCNA), specific gene co-expression networks formed in Chinese cabbage during vernalization were identified in order to find the reasons for the different bolting.

Results

RNA sequencing and gene co-expression network construction

Pearson’s correlation coefficients were used to test for biologically repeated correlations between samples. The generated cluster dendrogram was used to observe the overall correlation of the transcriptomes of the 2 Chinese cabbage accessions at different time periods (Fig. 1a). The three biological replicates from each time period and the transcriptome data both exhibited good correlation. The similarity test between the three biological replicates required the use of a principal component analysis (PCA). Using the first principal component (PC1) and second principal component (PC2), a dimensionality reduction analysis was used to analyze the similarity between each replicate (Fig. 1b). A total of 14 groups exhibited good similarity. Approximately 59.37% of the expressed genes were within the 0–5 FPKM range and 37.36% were within the 5–100 FPKM range (Fig. 1c).

Fig. 1
figure1

Transcriptional relationship between samples. a Heatmap of correlation value (R square) of 42 libraries. b Principal component analysis based on all of the expressed genes, showing 14 distinct groups of samples. c Number of transcripts in the 2 Chinese cabbage accessions, based on the FPKM of different samples

After analyzing the transcriptome data of each treatment period of 2 Chinese cabbage accessions, low abundance and low variability genes were filtered out. A total of 5748 genes of ‘JWW’ and 5527 genes of ‘XBJ’ were screened out. After being log2-transformed, they were imported into the WGCNA software package for analysis. WGCNA analysis performed transcriptome data analysis in each period. Each tree branch formed a module and each leaf in the branch represented a gene, as shown in the hierarchical clustering tree (Fig. 2). Then, the tree from the dendrogram was cut into modules (clusters). Based on their correlation with vernalization and control time, sets of genes (modules) were identified. As shown in the tree dendrogram, WGCNA analysis resulted in 9 modules that were distinguishable by different colors for ‘JWW’; the number of target genes for each module ranged from 56 to 3685 (Table S1). WGCNA analysis resulted in 12 modules that were distinguishable by different colors for ‘XBJ’; the number of target genes for each module ranged from 36 to 3745 (Table S2). Each module corresponded to each period and had its correlation. Whether the correlation was positive or negative and the size of the correlation showed the degree of correlation with the target gene screened out by the transcriptome data of this period (Figs. 3 and 4a).

Fig. 2
figure2

WGCNA of gene expression in ‘JWW’ (a) and ‘XBJ’ (b) during vernalization. Hierarchical cluster trees show the co-expression modules identified by WGCNA

Fig. 3
figure3

Co-expression modules for ‘JWW’. a Relationships between modules (left) and traits (bottom). Red and blue represent positive and negative correlations, respectively, with coefficient values and p-values. b Pairwise correlation coefficients between modules. Rows and columns are the module names, numbers represent coefficient values and p-values

Fig. 4
figure4

Co-expression modules for ‘XBJ’. a Relationships between modules (left) and Traits (bottom). Red and blue represent positive and negative correlations, respectively, with coefficient values and p-values. b Pairwise correlation coefficients between modules. Rows and columns are the module names, numbers represent coefficient values and p-values

Different modules related to ‘JWW’ and ‘XBJ’ in different periods

Module-trait relationships (MTRs) were different for each vernalization and control time period. These modules contained positively and negatively related genes, and their expression levels changed at different periods. Modules with MTR > 0.7 were selected as representatives of the 2 Chinese cabbage accessions for further analysis. Five modules were selected for both ‘JWW’ and ‘XBJ’. The results revealed the following high correlations: MEbrown (r = 0.93, p = 2e− 09) in J1 days after treatment (0 DAT); MEgreenyellow (r = 0.7, p = 4e− 04) in J2 (25 DAT); MEdarkgrey (r = 0.98, p = 2e− 15) in J4 (35 DAT); MEgrey60 (r = 0.84, p = 2e− 06) in J5 (45 DAT); MEblue (r = 0.98, p = 5e− 15) in JCK (35 DAT 25 °C) (Fig. 3a); MEturquoise (r = 0.98, p = 2e− 14) in X1 (0 DAT); MEdarkgreen (r = 0.73, p = 2e− 04) in X3 (15 DAT); MEpurple (r = 0.87, p = 4e− 07) in X4 (25 DAT); MEblack (r = 0.84, p = 2e− 06) in X6 (50 DAT); and MEcyan (r = 0.99, p = 8e− 18) in XCK (25 DAT 25 °C) (Fig. 4a).

The correlations between different modules of the 2 Chinese cabbage accessions were further investigated. Based on the eigengenes of each module, some module pairs were found to be significantly positively correlated. In ‘JWW’, MEdarkturquoise was positively correlated with MEgreenyellow (r = 0.82, p = 0.001) and MEblue and MEcyan were positively correlated (r = 0.82, p = 0.002) (Fig. 3b). In ‘XBJ’, MElightyellow was positively correlated with MEdarkgreen (r = 0.83, p = 8e− 04), MEgreenyellow was positively correlated with MElightgreen (r = 0.82, p = 0.001) and MEpurple (r = 0.81, p = 0.002) and MElightgreen was positively correlated with MEcyan (r = 0.81, p = 0.002)), MElightgreen was positively correlated with MEcyan (r = 0.81, p = 0.002) (Fig. 4b). Expression gene displays were performed for each Chinese cabbage processing stage and corresponded with each module (Fig. 5). Results revealed that the enrichment and differential expression displays from the co-expression network modules exhibited similar characteristics.

Fig. 5
figure5

Gene expression levels in ‘JWW’ (a) and ‘XBJ’ (b) with their corresponding log2FPKM module values. The color gradient from blue to red indicates high to low gene expression

Biological function analysis of important co-expression network modules

GO annotations and biological function analysis were performed using 10 modules that were highly related (Figs. 6 and 7). Brassica genes were first used as queries. When the Brassica database was insufficient, Arabidopsis orthologue genes were used as queries. GO terms were derived from these annotations (Table S3; Table S4).

Fig. 6
figure6

Significant GO terms and ontological relationships (annotated from ClueGO) in ‘JWW’. The sizes of the circles represent the degree of the positive relationship between the significant GO terms. Redundant terms were grouped and presented in the same color. Each leading term, which has the highest significance, is indicated by colored font

Fig. 7
figure7

Significant GO terms and ontological relationships (annotated from ClueGO) in ‘XBJ’. The sizes of the circles represent the degree of the positive relationship between the significant GO terms. Redundant terms were grouped and presented in the same color. Each leading term, which has the highest significance, is indicated by colored font

The biological functional terms enriched in ‘JWW’ MEbrown and ‘XBJ’ MEturquoise exhibited high correlation at 0 DAT and were the largest modules (p ≤ 0.01). In the Brassica database, ‘JWW’ MEbrown and ‘XBJ’ MEturquoise were enriched together with photosynthesis, response to cytokinin, chlorophyll biosynthetic process, and response to karrikin. The differences were ribosome biogenesis, translation, and response to unfolded protein, which were enriched in ‘JWW’ MEbrown, and light harvesting in photosystem I, protein-chromophore linkage, and reductive pentose-phosphate cycle, which were enriched in ‘XBJ’ MEturquoise. In the Arabidopsis Database, photosynthesis was the most enriched functional term in ‘JWW’ MEbrown and ‘XBJ’ MEturquoise. Additionally, cellular biosynthetic process, plastid organization, and anion transport were enriched in ‘JWW’ MEbrown, while cellular response to hormone stimulus, cellular response to endogenous stimulus, and cellular response to organic substance were enriched in ‘XBJ’ MEturquoise. These results indicated that the two Chinese cabbages had a certain degree of commonality to a large extent when they were not vernalized, and that when vernalized their different biological functions and gene expression might be observable.

‘JWW’ MEgreenyellow and ‘XBJ’ MEpurple were highly correlated at 25 DAT. The most enriched biological functional term in ‘JWW’ MEgreenyellow was cell wall organization in both the Brassica and Arabidopsis databases. In ‘XBJ’ MEpurple, the most enriched biological functional term in the Brassica database was xyloglucan metabolic process, while it was cell wall organization in the Arabidopsis database. In ‘JWW’ MEgreenyellow, several important biological functional terms were enriched, including cell wall biogenesis, carbohydrate metabolic process, and phenylpropanoid metabolic process. At 25 DAT, rapid flowering in ‘XBJ’ was promoted and was highly related to MEpurple. Biological functional terms related to polysaccharide metabolism processes were enriched, including polysaccharide metabolic process, cellular polysaccharide metabolic process, cell wall polysaccharide metabolic process, glucan metabolic process, cellular glucan metabolic process, and xyloglucan metabolic process. Additionally, shoot system morphogenesis was also enriched in this module. Thus, it was speculated that polysaccharide metabolism processes were enriched at 25 DAT in ‘XBJ’ to ensure that it transitioned from vegetative to reproductive growth, which was manifested by changes in shoot system morphogenesis.

‘JWW’ MEdarkgrey, which was highly correlated at 35 DAT, promoted rapid flowering and had many functional terms that were enriched in both databases, including response to water deprivation, response to chitin, abscisic acid (ABA)-activated signaling pathway, and response to UV-B. Additionally, response to stimulus, ethylene-activated signaling pathway, and aromatic amino acid family catabolic process, along with other terms, were positively regulated and enriched. These terms were enriched at 35 DAT during the critical vernalization period and may be the key biological functions that explain the transformation of late-bolting Chinese cabbage flowering.

MEdarkgreen, which was highly correlated with ‘XBJ’ at 15 DAT, was enriched in the functional terms nitric oxide biosynthetic process, glycolytic process, pyridine-containing compound metabolic process, sulfur amino acid metabolic process, and nitrogen cycle metabolic process, among other functional terms. The most enriched functional terms in ‘JWW’ MEgrey60 at 45 DAT included response to cold, circadian rhythm, response to temperature stimulus, and anthocyanin-containing compound metabolic process.

At 50 DAT, which was the largest vernalization period, ‘XBJ’ MEblack was enriched in functional terms related to hormones and amino acids, including response to ethylene, negative regulation of ethylene-activated signaling pathway, response to hormone, hormone-mediated signaling pathway, cellular response to hormone stimulus, amino acid export, and amino acid transmembrane transport. Additionally, reproductive growth and terms related to senescence were also enriched in this module, including positive regulation of leaf senescence, stress-induced premature senescence, and plant organ senescence.

‘JWW’ MEblue at 35 DAT at 25 °C, which was correlated with ‘JWW’ at 35 DAT in the control treatment, was enriched in the regulation of protein serine/threonine phosphatase activity, response to organic substance, hormone-mediated signaling pathway, and regulation of cellular process, among other functional terms. Notably, leaf senescence was negatively regulated and enriched in this module. Additionally, leaf senescence was positively regulated in ‘XBJ’ MEblack at 50 DAT, indicating that the leaf senescence of Chinese cabbage after vernalization may also signal bolting and flowering promotion. At 25 DAT, faster flowering was promoted in ‘XBJ’ MEcyan compared to 25 DAT at 25 °C, and ‘XBJ’ MEcyan was enriched in functional terms related to biosynthesis, including inositol biosynthetic process, aromatic compound biosynthetic process, small-molecule biosynthetic process, and wax biosynthetic process.

Hub gene selection for the ‘JWW’ and ‘XBJ’ co-expression networks

Hub genes were screened among these highly related modules across each time period. The top 20 genes that were representative of the modules were selected as they exhibited the largest “hubness” thereby providing the most detailed biological information (Figs. 8 and 9; Table S5; Table S6).

Fig. 8
figure8

Hub genes and expression profiles of ‘JWW’. a Co-expression gene networks with the greatest “hubness” in every module. Nodes are represented by dots surrounded by module colors. b log2FPKM expression profiles of the hub genes in J1, J2, J3, J4, J5, J6, and JCK. The locations of each gene correspond with A. The darker the green color, the higher the expression level

Fig. 9
figure9

Hub genes and expression profiles of ‘XBJ’. a Co-expression gene networks with the greatest “hubness” in every module. Nodes are represented by dots surrounded by module colors. b log2FPKM expression profiles of the hub genes in X1, X2, X3, X4, X5, X6, and XCK. The locations of each gene correspond with A. The darker the green color, the higher the expression level

MEgreenyellow, MEdarkgrey, and MEgrey60 were highly related modules in ‘JWW’ across vernalization periods. For MEgreenyellow, methylenetetrahydrofolate reductase 2 (MTHFR2), GDSL esterase/lipase CPRD49 (CPRD49), and amino acid permease 8 (AAP8) were enriched in amino acid transport and metabolism pathways. Carbohydrate transport and metabolism pathways were enriched in 3 genes among the 20 hub genes, including endoglucanase 10, beta-D-xylosidase 4 (BXL4), and beta-D-xylosidase 5 (BXL5). Moreover, the expression levels of MTHFR2, AAP8, BXL4, and BXL5 during vernalization were considerably higher than in the control treatment. Among the 20 genes expressed in MEdarkgrey: 3 glycosyl transferase family genes, GATL10s, and GATL17 were the hub genes of MEdarkgrey and their expression levels were the highest in ‘JWW’ J4 (35 DAT), and may be important family genes that promote faster flowering in ‘JWW’. Four AP2 domain genes, dehydration-responsive element-binding protein 1C (DREB1C), ethylene-responsive transcription factor 11 (RRF11), dehydration-responsive element-binding protein 1D (DREB1D), and ethylene-responsive transcription RAP2–13 were also hub genes found in this module. Notably, in MEdarkgrey, the B-box zinc finger protein 32 (BBX32) gene was enriched in biological functional terms related to flower development regulation. The top 20 hub genes in MEgrey60 included two 2-component response regulator-like APRR9 genes that were enriched in circadian rhythm-plant pathways.

Modules highly related to the ‘XBJ’ vernalization periods included MEdarkgreen, MEpurple, and MEblack. Of the top 20 hub genes in MEdarkgreen, 4 hub genes were involved in carbohydrate transport and metabolism, namely, BraA07g041160.3C, glucose-6-phosphate 1-dehydrogenase 3 (At1g24280), bifunctional enolase 2/transcriptional activator (ENO2), and 2,3-bisphosphoglycerate-independent phosphoglycerate mutase 1 (PGM1). In the MEgreenyellow of ‘JWW’, genes related to carbohydrate transport and metabolism pathways were also enriched, and the expression levels were notably higher than that of the control treatment, indicating that carbohydrate transport and metabolism may play an important role in the vernalization of Chinese cabbage. In MEpurple, 2 of the 20 hub genes encoded proteins and were enriched in the starch and sucrose metabolism pathway: trehalose-phosphate phosphatase A (TPPA) and sucrose synthase 1 (SUS1). Importantly, SUS1 participates in flower development. Three WRKY family genes existed as hub genes in MEblack, including WRKY transcription factor 18 (WRKY18), transcription factor 25 (WRKY25), and transcription factor 57 (WRKY57), of which, WRKY25 participated in the plant-pathogen interaction pathway. Among the top 20 hub genes in MEblack, 3 important genes were related to plant flowering, of which phophorbide a oxygenase (PAO) participated in flower development, protein twin sister of FT (TSF) regulated flower development and participated in photoperiodism and flowering, and NAC transcription factor 29 (NAC029) regulated flower development. The expression levels of these 3 genes across the ‘XBJ’ vernalization periods were significantly higher than those in the control treatment.

Validation of representative flower development-related hub genes expression

Five genes, BBX32, SUS1, PAO, TSF, and NAC029, correlated with flower development and were selected for verification by qRT-PCR. The RNA-Seq and qRT-PCR results were consistent (Fig. 10), indicating the reliability of high-throughput transcriptome sequencing. Compared with 0 DAT, the expression of SUS1 and NAC029 in ‘JWW’ at 25 DAT and 50 DAT were higher than in ‘XBJ’. The expression of TSF and BBX32 in ‘XBJ’ at 25 DAT and 50 DAT were higher than in ‘JWW’.

Fig. 10
figure10

qRT-PCR verification of the hub genes related to flower development. The bar and line graphs represent the qRT-PCR and RNA-Seq data, respectively. Data are presented as mean ± standard error (SE)

Discussion

Formation of specific co-expression networks using two bolting Chinese cabbage accessions and the WGCNA method

Phenotypic and molecular event-based RNA-Seq transcriptome analysis and WGCNA are powerful research methods [21,22,23]. WGCNA is a progressive analysis method in which variable genes are divided into co-expression modules through an unsigned network based on the gene expression patterns identified by RNA-Seq. Each module is then correlated with various traits and the gene “hubness” of each module builds the relationship between the positions of a single gene [24]. The eigengenes and hub genes of each module facilitate the establishment of the relationship between co-expressed gene clusters and concentrated traits in order to obtain clear expression patterns and screen candidate genes.

In this experiment, 2 Chinese cabbage accessions contained 21 RNA-Seq sample data points, respectively. Given that the expression of a large group of genes was affected by vernalization, WGCNA was used to construct a gene co-expression network to identify differences between modules (Fig. 2). The goal was to uncover the response mechanism of Chinese cabbage across different vernalization time periods and identify key genes. To our knowledge, there are no reports on the gene interaction networks of Chinese cabbage vernalization across different time periods. Therefore, based on the WGCNA gene co-expression network in this study, the responses of the 2 Chinese cabbage accessions to different vernalization stages were systematically analyzed at the transcriptome level.

Enrichment of different modules based on the transcriptomic differences of two Chinese cabbage accessions

Several modules that were highly related were selected for further analysis and discussion (Figs. 3 and 4). Based on the functions predicted by the modules of genes with known biological functions, the characteristics of the 2 Chinese cabbage accessions under vernalization and control treatments were analyzed and determined in order to find the reasons for their different bolting processes.

Photosynthesis, chlorophyll biosynthetic, cytokinin, and karrikin-responsive biological functional terms were enriched in the two most highly correlated modules of the 2 Chinese cabbage accessions at 0 DAT: MEbrown and MEturquoise. The contribution of photosynthesis to vegetative growth depended, to a large extent, on leaf area, chlorophyll content per leaf area, and chloroplast lifespan [25, 26]. Cytokinins are a class of hormones that regulate both the division cycle and meristem homeostasis [27, 28]. Karrikin is a plant growth regulator that promotes germination and seedling photomorphogenesis [29]. At 0 DAT, the 2 Chinese cabbage accessions were in the vegetative growth stage, thus, they exhibited obvious and consistent performance in terms of photosynthesis, chlorophyll biosynthetic, cytokinin, and karrikin. Using the vernalization of two Chinese cabbage accessions at 0 DAT as a starting period, the different performances of the two could be better analyzed and the reasons for their different bolting performances explored.

In flowering plants, includingg the model plant Arabidopsis thaliana, the shoot apical meristem (SAM) is the key determinant of overall morphogenesis [30]. We found that among highly correlated modules at 25 DAT, the MEpurple of ‘XBJ’ was more enriched in shoot system morphogenesis functional term than MEgreenyellow of ‘JWW’. Futhermore, at 25 DAT vernalization caused ‘XBJ’ to rapidly bolt and flower. This indicated that 25 DAT was a well-chosen treatment period for ‘XBJ’. Another finding was that MEpurple was enriched in many polysaccharide and sugar metabolism terms. The formation of SAM is controlled by the growth of plant cells, and the growth of plant cells is mainly controlled by the cell wall. The cell wall is a rigid structure composed primarily of polysaccharides that surround the cells and connect them together in a biological continuum [31]. This result may be explained by the formation of SAM in ‘XBJ’, which requires polysaccharides to regulate the elasticity of the cell wall to ensure cell elongation and growth. ‘XBJ’ generated energy under vernalization to ensure its successful transition from vegetative to reproductive growth. Judging from the 25 DAT enriched terms, compared with late-bolting ‘JWW’, early-bolting ‘XBJ’ was more susceptible to changes in shoot system morphogenesis and many polysaccharide and sugar metabolism functional terms under vernalization, which then promoted faster bolting and flowering.

MEdarkgrey, which was highly related to rapid flowering at 35 DAT in ‘JWW’, was enriched in the biological functional terms response to chitin, ABA-activated signaling pathway, and ethylene-activated signaling pathway. A previous study demonstrated that treatment with chitosan, a chitin derivative, in potted freesia plants caused early flowering and more flowers [32]. Another study demonstrated that endogenous ABA promoted bolting and flowering in plants after the promotion of FT and other related genes [33]. The ABA-activated signaling pathway could have been enriched during the vernalization of ‘JWW’, indicating that endogenous ABA played a certain role in the promotion of flowering in ‘JWW’. The effect of ethylene on floral transition is a complex biological process, as ethylene regulates this process by cooperating with other hormones or signal transduction pathways [34]. This finding corresponds with the ethylene-activated signaling pathway, which was enriched at 35 DAT in ‘JWW’. These biological functional terms were enriched in MEdarkgrey in ‘JWW’ and may be the key determinants of late-bolting Chinese cabbage floral transition.

In previous studies, age-dependent leaf senescence was found to be affected by developmental processes such as flowering, and the age-dependent leaf senescence phenotypes of circadian clock mutants showed significant correlation with flowering time [35, 36]. We found an interesting phenomenon in that the MEblack module, which was highly correlated with ‘XBJ’ at 50 DAT, was enriched with terms related to senescence, especially the positive regulation of leaf senescence, while in ‘JWW’ MEblue at 35 DAT at 25 °C, the negative regulation of leaf senescence term was enriched. Fifty DAT was the longest time for vernalization. At this time, early-bolting ‘XBJ’ had transitioned from vegetative growth to reproductive growth. While 35 DAT at 25 °C ‘JWW’ exhibited vegetative growth under normal growth conditions, leaves grew normally and the senescence phenomenon was in a state of resistance. With this finding, we speculated that leaf senescence was also a signal that promoted bolting and flowering of Chinese cabbage. This also laid the groundwork for our future research on the effect of vernalization on leaf senescence and flowering.

In summary, we tried to find the biological function of two different bolting Chinese cabbage accessions during the vernalization process: the early-bolting ‘XBJ’ could bolt and flower faster at 25 DAT, which was promoted by the shoot system morphogenesis and polysaccharide and sugar metabolism, while late-bolting ‘JWW’ enriched chitin, ABA, and ethylene-activated signaling pathways at 35 DAT, indicating that these regulatory pathways may promote bolting resistance in Chinese cabbage. An interesting finding was that the regulation of leaf senescence was found in the 2 Chinese cabbage-related modules, indicating that leaf senescence may be related to bolting and flowering.

Analysis of hub genes enriched in two Chinese cabbage accessions during vernalization

WGCNA was used to construct the gene co-expression networks of 2 Chinese cabbage accessions and analyze the modules that were highly related to their vernalization periods. The top 20 hub genes with the highest correlation relationships among these modules were identified to further analyze key candidate vernalization genes for Chinese cabbage with different bolting performances.

Amino acids are important constituents of proteins that play important roles in many pathways of the plant body, acting as biological stimulants under abiotic and biotic stress [37, 38]. Carbohydrates play a vital role in plant growth, reproduction, and flowering [39]. In this study, three hub genes, MTHFR2, CPRD49, and AAP8, were enriched in amino acid transport and metabolism pathways. Three other hub genes, endoglucanase 10, BXL4, and BXL5, were enriched in carbohydrate transport and metabolism pathways in ‘JWW’ MEgreenyellow. MTHFR is the least well-known enzyme in the folate-mediated 1-carbon metabolism of plants. MTHFR reactions in plants metabolize the methyl group 5,10-methylenetetrahydrofolate into serine, sugar, and starch [40]. AAP8 has been studied in seeds and siliques as an amino acid transporter and was specifically expressed in mature siliques [41]. Previous studies demonstrated that beta-D-xylosidase was widely expressed in plant flowers, siliques, and the SAM [42, 43]. The high expression levels of MTHFR2, AAP8, BXL4, and BXL5 during the vernalization of ‘JWW’ indicated that they were affected by vernalization and may have had an auxiliary promotion effect on the floral transformation of ‘JWW’. Three glycosyl transferase family genes, including two GATL10 genes and one GATL17 gene, were the main hub genes of ‘JWW’ MEdarkgrey. This finding was consistent with the results obtained from the GO biological function analysis of MEdarkgrey, indicating that sugar metabolism at 35 DAT in ‘JWW’ played an important role and promoted flower transformation. WRKY proteins, an important transcription factor superfamily involved in plant development and stress responses, have been studied in monocotyledonous and dicotyledonous plants [44]. In this study, three WRKY genes, WRKY18, WRKY25, and WRKY57, were identified in the hub genes of ‘XBJ’ MEblack. WRKY genes promote defense-related gene expression and disease resistance [45,46,47]. Vernalization is a form of low-temperature stress for Chinese cabbage. In this study, under the vernalization treatment, WRKY-related genes were enriched in ‘XBJ’ at 50 DAT.

Additionally, five genes related to flower development, BBX32, SUS1, PAO, TSF, and NAC029, were expressed in ‘JWW’ MEdarkgrey and ‘XBJ’ MEpurple and MEblack. B-box (BBX) zinc finger proteins play critical roles in both vegetative and reproductive plant development [48]. A previous study proved that BBX32 in Arabidopsis is a clock gene that interacts with COL3 and enables COL3 to bind to the FT promoter, thereby promoting the transcriptional regulation of flowering time [49]. A separate study demonstrated that BrBBX32 binds to BrAGL24 in Chinese cabbage through the B-box domain, which regulates flowering time [50]. In this study, BBX32 was enriched in ‘JWW’ at 35 DAT, indicating that vernalization induced BBX32 expression and promoted the flowering transition of Chinese cabbage. A previous study showed that sucrose levels increased in the leaves and SAM of Arabidopsis exposed to strong radiation, thereby promoting bolting and flowering by increasing FT expression and inducing SUS1 expression [51]. Sugar levels regulate plant flowering [52, 53]. In this study, vernalization induced SUS1 levels, indicating that SUS1 can be used as a candidate gene for Chinese cabbage vernalization. TSF is the closest homologue of FT and transgenic plants that overexpress TSF exhibit premature flowering [54]. The high expression levels of TSF in MEblack, which was highly correlated with ‘XBJ’ at 50 DAT, indicated that ‘XBJ’ began reproductive growth at this time; TSF was also continuously expressed. PAO is a chloroplast envelope-bound Rieske-type iron-sulfur oxygenase. The degradation of chlorophyll in Arabidopsis is related to PAO activities [55, 56]. In this study, the expression of PAO reached its maximum level at 50 DAT in both Chinese cabbage accessions, indicating that Chinese cabbage was gradually aged during vernalization. NAC family transcription factors play a role in leaf senescence [57]. One previous study found that multiple NACs played regulatory roles in flowering [58]. In this study, under the vernalization treatment, the expression of NAC029 in both Chinese cabbage accessions increased and was significantly upregulated in ‘XBJ’. This finding demonstrated that the early-bolting ‘XBJ’ accession could more quickly adapt to vernalization and after the formation of the SAM, its leaves gradually aged. When we performed biological functions on the highly correlated modules of two Chinese cabbage treatment periods, we found that the positive and negative regulation of leaf senescence was enriched in the vernalization and non-vernalization periods, which was consistent with the high expression of PAO and NAC029 found here. This collectively proved that the vernalization process and the aging mechanism have a connection, although whether it promotes the flowering transition remains to be determined.

Conclusion

Vernalization as an important trait that is directly linked to production potential. It is necessary to elucidate the regulatory mechanisms involved in differently-bolting Chinese cabbage varieties. In this study, a WGCNA was conducted using RNA data from two Chinese cabbage accessions that bolt differently in order to reveal the key pathways and hub genes that cause bolting and flowering during vernalization. The results revealed that shoot system morphogenesis and polysaccharide and sugar metabolism induce early-bolting ‘XBJ’ to bolt and flower faster; chitin, ABA and ethylene-activated signaling pathways were enriched in late-bolting ‘JWW’; and leaf senescence and carbohydrate metabolism pathways were found to be enriched in the two Chinese cabbage-related modules, indicating that these may be related to bolting and flowering. Additionally, five genes related to flower development, BBX32 (binds to the FT promoter), SUS1 (increases FT expression), TSF (the closest homologue of FT), PAO, and NAC029 (plays a role in leaf senescence), revealed different vernalization mechanisms. The findings of this study provide a comprehensive overview of vernalization-related gene networks in Chinese cabbage and uncovered candidate hub genes in the vernalization process that can be utilized in future breeding research.

Methods

Chinese cabbage accessions

Two Chinese cabbage accessions, the late-bolting Jin Wawa (JWW) and early-bolting Xiao Baojian (XBJ), were provided by the Chinese Academy of Agricultural Sciences located in Beijing, China. The two materials were highly inbred lines. After plants were grown in a nursery greenhouse under normal conditions for 32 d, the vernalization experiment began and lasted for 50 d (Fig. 11a). The treatment conditions were 4 °C for the vernalization treatment and 25 °C for the control treatment.

Fig. 11
figure11

a Treatment process of two Chinese cabbage accessions: (a), (b), and (c) are ‘JWW’; (d), (c), and (f) are ‘XBJ’; (a) and (d) are before vernalization (0 DAT); (b) and (c) are during vernalization (25 DAT); (c) and (f) are after flowering. b ‘JWW’ flowering time promoted by vernalization. c ‘XBJ’ flowering time promoted by vernalization

Sample selection and RNA sequencing

Based on the timing of flowering caused by vernalization (Fig. 11b; c), the following samples were collected for RNA-Seq: from ‘JWW’, selected samples included J1 (0 days after treatment (DAT)), J2 (25 DAT), J3 (30 DAT), J4 (35 DAT), J5 (45 DAT), J6 (50 DAT), and JCK (35 DAT 25 °C), and from ‘XBJ’, selected samples included X1 (0 DAT), X2 (10 DAT), X3 (15 DAT), X4 (25 DAT), X5 (40 DAT), X6 (50 DAT), and XCK (25 DAT 25 °C). Three biological replicates were collected for each sample. The vernalization treatment, sample collection method, RNA-seq period selection, and the detailed methods used for data processing are described in our previous study [59].

Gene co-expression network construction and visualization

The RNA-Seq data were analyzed to construct gene co-expression networks using the R package, WGCNA [60]. Based on the following criteria, Fragments Per Kilobase of transcript per Million mapped reads (FPKM) ≥ 1, and a variation of FPKM: cv ≥ 0.5 and cv ≤ sd (genes number)/mean (genes number) (‘sd’ represents the standard deviation of the sample, and ‘mean’ represents the calculated average of the sample),genes of the 2 Chinese cabbage accessions were screened for co-expression network construction. From ‘JWW’, 5748 co-constructed genes were screened out. The following parameters were used to identify each gene module: weighted network, unsigned; hierarchical clustering tree, dynamic hybrid tree cut algorithm; power, 5; and minimum module size, 30; minimum height for merging modules, 0.29995. From ‘XBJ’, 5527 co-constructed genes were screened out. The following parameters were used to identify each gene module: weighted network, unsigned; hierarchical clustering tree, dynamic hybrid tree cut algorithm; power, 5; minimum module size, 30; and minimum height for merging modules, 0.3131.

To describe the most common gene expression models in each module, module eigengenes were used. Module eigengenes are the first major components of the expression matrix and are used to summarize the module overview and feature data. Pearson’s correlation coefficients were used to calculate the correlation between the module characteristic genes and the degree of vernalization of the two Chinese cabbage accessions. A heat map was drawn according to the correlation coefficients. The depth of color represents the correlation between the module and the degree of vernalization.

Analysis of hub genes in the gene co-expression network

Hub genes are good representatives of each co-expression module and have important biological significance in the system analysis. Hub genes are genes with the most connection points in each module, and their height is represented by the kME value. The kME value is based on the Pearson correlation coefficient between the expression level and module eigengenes. The kME and eigengene connectivity of each gene are calculated by signedKME, which includes the edge and node characteristics. The genetic network map, which was drawn according to the kME values, was created using Cytoscape software [61].

A gene ontology (GO) enrichment analysis was conducted on the genes using GOseq R software [62] and ClueGO [63]. The Brassica Database v3.0 IDs [64] were used as search queries for the GOseq R software annotations; GO terms with FDR values < 0.01 were selected for output. The TAIR10 IDs were used as search queries for the ClueGO annotations. ClueGO is a cytoscape plugin for visualizing large gene clusters in a functionally grouped network that can analyze both single clusters and compare them based on their specificity and the same aspects of multiple cluster functions. The ClueGO network was set to ‘medium’ and its connectivity was based on a kappa score of 0.4. GO terms with p ≤ 0.01 were considered to be significant. Other parameters were based on the original ClueGO values. Gene functions were annotated based on the Swiss-Prot, KOG/COG, KO, Pfam, and Nr NCBI databases.

Quantitative Real-Time PCR (qRT-PCR) and the evaluation of candidate hub gene expression

Five hub genes were selected to evaluate their expression levels by qRT-PCR analysis. Gene-specific primers were designed using Primer v5.0. Actin was used as an internal control for gene expression (Table S7). The Bio-Rad CFX96 RT-PCR Detection system (Bio-Rad, Hercules, CA, USA) and SYBR Green II PCR Master mix (Takara, Nojihigashi, Kusatsu, Japan) were used for the qRT-PCR reactions. The gene expression data were analyzed using the 2-ΔΔCt method [65]. SPSS v19.0 (SPSS, Chicago, IL, USA) was used to conduct a one-way analysis of variance (ANOVA) with Duncan’s multiple range post-hoc test and a significance threshold of p < 0.05. Results were visualized using Sigmaplot v10.0 (Systat Software Inc., San Jose, CA, USA).

Availability of data and materials

The raw data have been submitted under BioProject number PRJNA615255 to the Sequence Read Archive (SRA) database at NCBI (https://www.ncbi.nlm.nih.gov/sra/PRJNA615255).

Abbreviations

JWW:

Jin Wawa

XBJ:

Xiao Baojian

DAT:

Days after treatment

RNA-seq:

RNA sequencing

FPKM:

Fragments Per Kilobase of transcript per Million mapped reads

WGCNA:

Weighted gene co-expression network analysis

GO:

Gene ontology

ABA:

Abscisic acid

DEGs:

Differentially expressed genes

PCA:

Principal components analysis

PC1:

First principal component

PC2:

Second principal component

MTRs:

Module-trait relationships

SAM:

Shoot apical meristem

qRT-PCR:

Quantitative real-time PCR

References

  1. 1.

    Laczi E, Apahidean AS, Apahidean AI, Varga J. Agrobiological Particularities of Chinese Cabbage: Brassica rapa L. ssp. pekinensis (Lour) Hanelt and Brassica rapa L. ssp. chinensis (Lour) Hanelt. Bull Univ Agri Sci Vet Med Cluj-Napoca. 2010;67(1):508.

    Google Scholar 

  2. 2.

    Eniko L, Alexandru SA, Jolan V, Alexandra IA. Research concerning the bolting of Chinese cabbage (Brassica campestris var. pekinensis (Lour.) Olson) in early crops in polyethylene tunnels. Acta Univ Sapientiae Agri Environ. 2011;3:152–60.

    Google Scholar 

  3. 3.

    Sung S, Amasino RM. Remembering Winter: toward a molecular understanding of Vernalization. Annu Rev Plant Biol. 2005;56(1):491–508. https://doi.org/10.1146/annurev.arplant.56.032604.144307.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Araki T. Transition from vegetative to reproductive phase. Curr Opin Plant Biol. 2001;4(1):63–8. https://doi.org/10.1016/S1369-5266(00)00137-0.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Kim DH, Sung S. Genetic and Epigenetic Mechanisms Underlying Vernalization. Arabidopsis Book. 2014;12(12). https://doi.org/10.1199/tab.0171.

  6. 6.

    Yui S, Yoshikawa H. Bolting resistant breeding of Chinese cabbage. 1. Flower induction of late bolting variety without chilling treatment. Euphytica. 1991;52(3):171–6. https://doi.org/10.1007/BF00029393.

    Article  Google Scholar 

  7. 7.

    Sheldon CC, Rouse DT, Finnegan EJ, Peacock WJ, Dennis ES. The molecular basis of vernalization: the central role of FLOWERING LOCUS C (FLC). Proc Natl Acad Sci U S A. 2000;97(7):3753–8. https://doi.org/10.1073/pnas.97.7.3753.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Xiao X, Wu C, Xu Z, Yang Y, Fan S, Wang H. Molecular cloning, characterization and expression analysis of bolting-associated genes in flowering Chinese cabbage. Genes Genomics. 2015;37(4):357–63. https://doi.org/10.1007/s13258-014-0264-z.

    CAS  Article  Google Scholar 

  9. 9.

    Lin S, Wang J, Poon S, Su C, Wang S, Chiou T. Differential regulation of FLOWERING LOCUS C expression by Vernalization in cabbage and Arabidopsis. Plant Physiol. 2005;137(3):1037–48. https://doi.org/10.1104/pp.104.058974.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Takada S, Akter A, Itabashi E, Nishida N, Shea DJ, Miyaji N, Mehraj H, Osabe K, Shimizu M, Takasakiyasuda T. The role of FRIGIDA and FLOWERING LOCUS C genes in flowering time of Brassica rapa leafy vegetables. Sci Rep. 2019;9(1):13843. https://doi.org/10.1038/s41598-019-50122-2.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Zhang S, Li X, Sun C, He Y. Epigenetics of plant vernalization regulated by non-coding RNAs. Hereditas. 2012;34(7):829–34. https://doi.org/10.3724/sp.j.1005.2012.00829.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Li Z, Zhao L, Cui C, Kai G, Zhang L, Sun X, Tang K. Molecular cloning and characterization of an anti-bolting related gene (BrpFLC) from Brassica rapa ssp. pekinensis. Plant Sci. 2005;168(2):407–13. https://doi.org/10.1016/j.plantsci.2004.08.012.

    CAS  Article  Google Scholar 

  13. 13.

    Schranz M, Quijada P, Sung S, Lukens L, Amasino R, Osborn T. Characterization and effects of the replicated flowering time gene FLC in Brassica rapa. Genetics. 2002;162(3):1457–68.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Xiao D, Zhao J, Hou X, Basnet R, Carpio D, Zhang N, Bucher J, Lin K, Cheng F, Wang X, Bonnema G. The Brassica rapa FLC homologue FLC2 is a key regulator of flowering time, identified through transcriptional co-expression networks. J Exp Bot. 2013;64(14):4503–16. https://doi.org/10.1093/jxb/ert264.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Xi X, Wei K, Gao B, Liu J, Liang J, Cheng F, Wang X, Wu J. BrFLC5: a weak regulator of flowering time in Brassica rapa. Theo Appl Genet. 2018;131(10):2107–16.

    CAS  Article  Google Scholar 

  16. 16.

    Bhargava A, Clabaugh I, To JP, Maxwell BB, Chiang YH, Schaller GE, Loraine AE, Kieber JJ. Identification of Cytokinin-responsive genes using microarray meta-analysis and RNA-Seq in Arabidopsis. Plant Physiol. 2013;162(1):272–94. https://doi.org/10.1104/pp.113.217026.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Shen C, Li D, He R, Fang Z, Xia Y, Gao J, Shen H, Cao M. Comparative transcriptome analysis of RNA-seq data for cold-tolerant and cold-sensitive rice genotypes under cold stress. J Plant Biol. 2014;57(6):337–48. https://doi.org/10.1007/s12374-014-0183-1.

    CAS  Article  Google Scholar 

  18. 18.

    Li Y, Tian S, Yang X, Wang X, Guo Y, Ni H. Transcriptomic analysis reveals distinct resistant response by physcion and chrysophanol against cucumber powdery mildew. PeerJ. 2016;4:e1991. https://doi.org/10.7717/peerj.1991.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Sun M, Qi X, Hou L, Xu X, Zhu Z, Li M. Gene Expression Analysis of Pak Choi in Response to Vernalization. PLoS One. 2015;10(10).

  20. 20.

    Qi X, Townsley BT, Aguilarmartinez JA, Yin L, Gao X, Hou L, Gao M, Li M. Cloning and characterization of the LFY homologue from Chinese cabbage (Brassica rapa subsp. pekinensis). Hortic Environ Biotechnol. 2015;56(6):821–9. https://doi.org/10.1007/s13580-015-0066-5.

    CAS  Article  Google Scholar 

  21. 21.

    Shao Z, Zhang P, Lu C, Li S, Chen Z, Wang X, Duan D. Transcriptome sequencing of Saccharina japonica sporophytes during whole developmental periods reveals regulatory networks underlying alginate and mannitol biosynthesis. BMC Genomics. 2019;20(1):975. https://doi.org/10.1186/s12864-019-6366-x.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Sari E, Cabral AL, Polley B, Tan Y, Hsueh E, Konkin D, Knox RE, Ruan Y, Fobert PR. Weighted gene co-expression network analysis unveils gene networks associated with the Fusarium head blight resistance in tetraploid wheat. BMC Genomics. 2019;20(1):925. https://doi.org/10.1186/s12864-019-6161-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Li L, Long Y, Li H, Wu X. Comparative Transcriptome analysis reveals key pathways and hub genes in rapeseed during the early stage of Plasmodiophora brassicae infection. Front Genet. 2020;10:1275. https://doi.org/10.3389/fgene.2019.01275.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Dai R, Xia Y, Liu C, Chen C. csuWGCNA: a combination of signed and unsigned WGCNA to capture negative correlations. bioRxiv. 2018:288225.

  25. 25.

    Offermann S, Peterhansel C. Can we learn from heterosis and epigenetics to improve photosynthesis. Curr Opin Plant Biol. 2014;19:105–10. https://doi.org/10.1016/j.pbi.2014.05.010.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Apel K, Kloppstech K. The effect of light on the biosynthesis of the light-harvesting chlorophyll a/b protein : evidence for the requirement of chlorophyll a for the stabilization of the apoprotein. Planta. 1980;150(5):426–30. https://doi.org/10.1007/BF00390180.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Ferreira FJ, Kieber JJ. Signaling: Cytokinin Signaling. Curr Opin Plant Biol. 2005;8(5):518–25. https://doi.org/10.1016/j.pbi.2005.07.013.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Zhao Y. The role of local biosynthesis of auxin and cytokinin in plant development. Curr Opin Plant Biol. 2008;11(1):16–22. https://doi.org/10.1016/j.pbi.2007.10.008.

    CAS  Article  PubMed  Google Scholar 

  29. 29.

    Baldrianova J, Cerný M, Novak J, Jedelský P, Diviskova E, Brzobohatý B. Arabidopsis proteome responses to the smoke-derived growth regulator karrikin. J Proteome. 2015;120:7–20. https://doi.org/10.1016/j.jprot.2015.02.011.

    CAS  Article  Google Scholar 

  30. 30.

    Eric M, Elliot M, Victoria G, Tobias M. Morphogenesis in plants: modeling the shoot apical meristem, and possible applications. Workshop Evolvable Hardw. 1999.

  31. 31.

    Sassi M, Traas J: New insights in shoot apical meristem morphogenesis: isotropy comes into play. Plant Signal Behav 2015, 10(11):0, e1000150, DOI: https://doi.org/10.1080/15592324.2014.1000150.

  32. 32.

    Salachna P, Agnieszka Z. Effect of chitosan on plant growth, flowering and corms yield of potted freesia. J Ecol Eng. 2014;15(3):97–102.

    Google Scholar 

  33. 33.

    Riboni M, Galbiati M, Tonelli C, Conti L. GIGANTEA enables drought escape response via Abscisic acid-dependent activation of the Florigens and suppressor of overexpression of constans1. Plant Physiol. 2013;162(3):1706–19. https://doi.org/10.1104/pp.113.217729.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Iqbal N, Khan NA, Ferrante A, Trivellini A, Francini A, Khan MIR. Ethylene role in plant growth, development and senescence: interaction with other Phytohormones. Front Plant Sci. 2017;8:475.

    PubMed  PubMed Central  Google Scholar 

  35. 35.

    Kim H, Kim H, Vu Q, Jung S, McClung C, Hong S, Nam H. Circadian control of ORE1 by PRR9 positively regulates leaf senescence in Arabidopsis. PNAS. 2018;115(33):8448–53. https://doi.org/10.1073/pnas.1722407115.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Song Y, Yang C, Gao S, Zhang W, Li L, Kuai B. Age-triggered and dark-induced leaf senescence require the bHLH transcription factors PIF3, 4, and 5. Mol Plant. 2014;7(12):1776–87. https://doi.org/10.1093/mp/ssu109.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  37. 37.

    Heuer B. Influence of exogenous application of proline and glycinebetaine on growth of salt-stressed tomato plants. Plant Sci. 2003;165(4):693–9. https://doi.org/10.1016/S0168-9452(03)00222-X.

    CAS  Article  Google Scholar 

  38. 38.

    Maini P, Bertucci BM. Possibility to reduce the effects of the viruses with a biostimulant based on amino acids and peptides. Agro Food Industry Hi Tech. 1999;10(5):26–8.

    CAS  Google Scholar 

  39. 39.

    Ellithy M, Reymond M, Stich B, Koornneef M, Vreugdenhil D. Relation among plant growth, carbohydrates and flowering time in the Arabidopsis Landsberg erecta × Kondara recombinant inbred line population. Plant Cell Environ. 2010;33(8):1369–82.

    CAS  Google Scholar 

  40. 40.

    Roje S, Wang H, Mcneil SD, Raymond RK, Appling DR, Shacharhill Y, Bohnert HJ, Hanson AD. Isolation, characterization, and functional expression of cDNAs encoding NADH-dependent Methylenetetrahydrofolate Reductase from higher plants. J Biol Chem. 1999;274(51):36089–96. https://doi.org/10.1074/jbc.274.51.36089.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Okumoto S, Schmidt R, Tegeder M, Fischer WN, Rentsch D, Frommer WB, Koch W. High affinity amino acid transporters specifically expressed in xylem parenchyma and developing seeds of Arabidopsis. J Biol Chem. 2002;277(47):45338–46. https://doi.org/10.1074/jbc.M207730200.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Goujon T, Minic Z, Amrani AE, Lerouxel O, Aletti E, Lapierre C, Joseleau J, Jouanin L. AtBXL1, a novel higher plant (Arabidopsis thaliana) putative beta-xylosidase gene, is involved in secondary cell wall metabolism and plant development. Plant J. 2003;33(4):677–90. https://doi.org/10.1046/j.1365-313X.2003.01654.x.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Minic Z, Rihouey C, Do CT, Lerouge P, Jouanin L. Purification and characterization of enzymes exhibiting β-d-Xylosidase activities in stem tissues of Arabidopsis. Plant Physiol. 2004;135(2):867–78. https://doi.org/10.1104/pp.104.041269.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Chen F, Hu Y, Vannozzi A, Wu K, Cai H, Qin Y, Mullis A, Lin Z, Zhang L. The WRKY transcription factor family in model plants and crops. Crit Rev Plant Sci. 2018;36:311–35.

    Article  Google Scholar 

  45. 45.

    Chen C, Chen Z. Potentiation of developmentally regulated plant defense response by AtWRKY18, a pathogen-induced Arabidopsis transcription factor. Plant Physiol. 2002;129(2):706–16. https://doi.org/10.1104/pp.001057.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Li S, Fu Q, Chen L, Huang W, Yu D. Arabidopsis thaliana WRKY25, WRKY26, and WRKY33 coordinate induction of plant thermotolerance. Planta. 2011;233(6):1237–52. https://doi.org/10.1007/s00425-011-1375-2.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Jiang Y, Liang G, Yu D. Activated expression of WRKY57 confers drought tolerance in Arabidopsis. Mol Plant. 2012;5(6):1375–88. https://doi.org/10.1093/mp/sss080.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Huang J, Zhao X, Weng X, Wang L, Xie W. The Rice B-Box Zinc Finger Gene Family: Genomic Identification, Characterization, Expression Profiling and Diurnal Analysis. PLoS One. 2012;7(10).

  49. 49.

    Tripathi P, Carvallo M, Hamilton EE, Preuss S, Kay SA. Arabidopsis B-BOX 32 interacts with CONSTANS-LIKE3 to regulate flowering. Proc Natl Acad Sci U S A. 2017;114(1):172–7. https://doi.org/10.1073/pnas.1616459114.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Guanpeng MA, Zhao D, Wang T, Zhou L, Guilian LI. BBX32 interacts with AGL24 involved in flowering time control in Chinese cabbage (Brassica rapa L. ssp. pekinensis). Notulae Botan Horti Agrobotanici Cluj-napoca. 2018;47(1):34–45.

    Article  Google Scholar 

  51. 51.

    King RW, Hisamatsu T, Goldschmidt EE, Blundell C. The nature of floral signals in Arabidopsis. I. Photosynthesis and a far-red photoresponse independently regulate flowering by increasing expression of FLOWERING LOCUS T (FT). J Exp Bot. 2008;59(14):3811–20. https://doi.org/10.1093/jxb/ern231.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Lebon G, Wojnarowiez G, Holzapfel B, Fontaine F, Vaillantgaveau N, Clement C. Sugars and flowering in the grapevine (Vitis vinifera L.). J Exp Bot. 2008;59(10):2565–78. https://doi.org/10.1093/jxb/ern135.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Seo PJ, Ryu JY, Kang SK, Park C. Modulation of sugar metabolism by an INDETERMINATE DOMAIN transcription factor contributes to photoperiodic flowering in Arabidopsis. Plant J. 2011;65(3):418–29. https://doi.org/10.1111/j.1365-313X.2010.04432.x.

    CAS  Article  PubMed  Google Scholar 

  54. 54.

    Yamaguchi A, Kobayashi Y, Goto K, Abe M, Araki T. TWIN SISTER OF FT (TSF) acts as a floral pathway integrator redundantly with FT. Plant Cell Physiol. 2005;46(8):1175–89. https://doi.org/10.1093/pcp/pci151.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Pružinska A, Tanner G, Aubry S, Anders I, Moser S, Muller T, Ongania K, Krautler B, Youn J, Liljegren SJ. Chlorophyll breakdown in senescent Arabidopsis leaves. Characterization of chlorophyll Catabolites and of chlorophyll catabolic enzymes involved in the Degreening reaction. Plant Physiol. 2005;139(1):52–63. https://doi.org/10.1104/pp.105.065870.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Chung DW, Pružinska A, Hortensteiner S, Ort DR. The role of Pheophorbide a Oxygenase expression and activity in the canola green seed problem. Plant Physiol. 2006;142(1):88–97. https://doi.org/10.1104/pp.106.084483.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Guo Y, Gan S. AtNAP, a NAC family transcription factor, has an important role in leaf senescence. Plant J. 2006;46(4):601–12. https://doi.org/10.1111/j.1365-313X.2006.02723.x.

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    Sablowski R, Meyerowitz EM. A homolog of NO APICAL MERISTEM is an immediate target of the floral homeotic genes APETALA3/PISTILLATA. Cell. 1998;92(1):93–103. https://doi.org/10.1016/S0092-8674(00)80902-2.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Dai Y, Zhang S, Sun X, Li G, Yuan L, Li F, Zhang H, Zhang S, Chen G, Wang C, Sun R. Comparative Transcriptome analysis of gene expression and regulatory characteristics associated with different Vernalization periods in Brassica rapa. Genes. 2020;11(4):392. https://doi.org/10.3390/genes11040392.

    CAS  Article  PubMed Central  Google Scholar 

  60. 60.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559. https://doi.org/10.1186/1471-2105-9-559.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Lopes CT, Max F, Farzana K, Donaldson SL, Quaid M, Bader GD. Cytoscape web: an interactive web-based network browser. Bioinformatics. 2010;18:18.

    Google Scholar 

  62. 62.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):1–12.

    Article  Google Scholar 

  63. 63.

    Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman W, Pages F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3. https://doi.org/10.1093/bioinformatics/btp101.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Zhang L, Cai X, Wu J, Liu M, Grob S, Cheng F, Liang J, Cai C, Liu Z, Liu B. Improved Brassica rapa reference genome by single-molecule sequencing and chromosome conformation capture technologies. Horticult Res. 2018;5(1):1–11.

    Google Scholar 

  65. 65.

    Dai Y, Yuan L, Zhang S, Wang J, Xie S, Zhao M, Chen G, Sun R, Wang C. Comprehensive Evaluation for Cold Tolerance in Wucai (Brassica campestris L.) by the Performance Index on an Absorption Basis (PIabs). Agronomy. 2019;9(2):61.

    CAS  Article  Google Scholar 

Download references

Acknowledgements

We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

Funding

This work was funded by the National Key Research and Development Program of China (grant number 2017YFD0101802) and the China Agricultural Research System (CARS-23-A-02). This work was performed at the Key Laboratory of Biology and Genetic Improvement of Horticultural Crops, Ministry of Agriculture, Beijing, China.

Author information

Affiliations

Authors

Contributions

Conceived and designed the experiments: SZ2 (Shujiang Zhang), RS, and CW. Performed the experiment: YD and XS. Wrote the manuscript: YD and XS. Analyzed data and submitted data: FL, SZ1 (Shifan Zhang) and HZ. Modified the manuscript: GL, LY, and GC. All of the authors have read and approved the final manuscript.

Corresponding author

Correspondence to Shujiang Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Log2FPKM values of the 5748 variable genes in ‘JWW’ (J for ‘JWW’). Table S2. Log2 FPKM values of the 5527 variable genes in ‘XBJ’ (X for ‘XBJ’). Table S3. Significant GO terms for each gene set of ‘JWW’. Table S4. Significant GO terms for each gene set of ‘XBJ’. Table S5. Annotation information of the hub genes of each module in ‘JWW’. Table S6. Annotation information of the hub genes of each module in ‘XBJ’. Table S7. qRT-PCR specific primers for hub gene related to flower development.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Dai, Y., Sun, X., Wang, C. et al. Gene co-expression network analysis reveals key pathways and hub genes in Chinese cabbage (Brassica rapa L.) during vernalization. BMC Genomics 22, 236 (2021). https://doi.org/10.1186/s12864-021-07510-8

Download citation

Keywords

  • Chinese cabbage
  • Gradient-vernalization
  • RNA sequencing
  • Weighted gene co-expression network analysis
  • Hub genes
\