Comparative Transcriptome Analysis Between Cold-Tolerant and Sensitive Brassica Oleracea L. in Response to Low Temperature Stress

Background: Brassica oleracea L. occupies an important position in the annual production of vegetables. But during winter Brassica oleracea L. often suffers from low temperatures and even sub-zero temperatures. Through transcriptome analysis and identication, the pathways involved in cold tolerance of Brassica oleracea L. were analyzed and candidate genes related to cold tolerance of Brassica oleracea L. were identied. Results: Under low temperature stress, a large number of signicantly different genes were found in Zhonggan1229 (ZG, low temperature tolerance) and Yingchun (YN, low temperature sensitive). There were 3902 signicantly up-regulated genes and 5309 signicantly down-regulated genes in ZG, and 4253 signicantly up-regulated genes and 5938 signicantly down-regulated genes in YN. Among them, 1844 different genes are the specic different genes in ZG and 6089 genes are the common different genes to response the low temperature stress. By annotating the specic different genes in ZG, 26 of the top 30 enriched GO terms belonged to biological processes, 4 terms belonged to molecular functions. By annotating the common different genes, 23 GO terms belonged to biological processes, 1 GO term belonged to molecular functions, and 6 GO terms belonged to cellular components. Circadian rhythms of plants and Plant hormone signal transduction were not only signicantly enriched in the two analyzed genes, but also the effects of low temperature stress were most signicant. Among the unique different genes in ZG, 154 genes were annotated into transcription factor families, and 79 genes were up-regulated and 75 genes were down-regulated, the encoding of MYB-related proteins was the largest group. Among the different genes shared by the two varieties, 516 genes were annotated into corresponding transcription factor families, 211 genes were up-regulated and 296 genes were down-regulated, however, there were 4 genes that were up-regulated in ZG but down-regulated in YN, and 5 genes that were down-regulated in ZG but up-regulated in YN, the largest group was the protein encoding ERF. Conclusions: The results identied important genes, pathways, and transcription factors that respond to low temperature stress, provided cold tolerance gene (EC 2 ) was determined. The electrical conductivity (EC 0 ) of the double steam water was taken as the control.


Background
Low temperature stress is an abiotic stress frequently occurring in agricultural production, which affects the growth and development of plants. After a long time of evolution, various plants have formed a set of molecular regulatory system to deal with cold stress. By regulating a series of physiological and biochemical reactions, they can adapt to low temperature stress [1,2]. Chinese scientists revealed that bio lms, calcium signaling, endocytosis and antioxidant capacity determine plant cold resistance at physiological and cellular levels [3][4][5][6]. At present, studies have been carried out on the low temperature stress of many plants, including Arabidopsis, rice, rape, tea, wheat, tobacco, sorghum and other important crops. Therefore, the research to the plant under low temperature damage effect and its mechanism, to explore the mechanism of plant resistance to freezing and its prevention measures, not only for cultivation has a theoretical signi cance, but also for development to foster more varieties of cold resistance of plant, increase the cold region in China plant biological yield and economic yield has extensive application value.
Brassica oleracea L. is an important vegetable crop widely planted by China and other countries, prefer mild and moist, adequate light. They are hardy and have the ability to adapt to high temperatures, and the suitable temperature for them to growth is 15-20°C. When the eshy stem expands, it is easy to brosis in high temperature above 30°C. In addition to this, the selection of soil is not very strict, but suitable for cultivation in humus-rich clay loam or sandy loam. Brassica oleracea L. seedlings must be vernalized at 0-10°C and then licked, blossomed and bore fruit under long sunshine and moderate temperature. Since the "13th Five-Year Plan", according to statistics, the planting area of Brassica oleracea L. in China is about 900,000 hm 2 [7], and the cultivated area is still expanding in recent years, and due to its characteristics of different maturity, suitable for different seasons, wide adaptability and strong stress resistance, Brassica oleracea L. occupies an important position in the annual production of vegetables [8]. But low temperature was one of the most important limiting factors for the wintering cultivation of Brassica oleracea L. in open elds, during winter Brassica oleracea L. often suffers from low temperatures and even sub-zero temperatures. The low temperature conditions make its yield and quality decline. Therefore, it is of great signi cance to reveal the cold resistance mechanism of Brassica oleracea L. for the study of winter cultivation and cold resistance breeding.
With the completion of the Human Genome Project, the rst-generation sequencing method has been unable to meet the requirements of large-scale genome sequencing such as deep sequencing and repeat sequencing, which has promoted the birth of the second-generation sequencing technology, its most remarkable characteristic is high throughput, low cost. Currently, Roche's 454 technology, Illumina's Solexa technology and ABI's SOLiD technology are the mainstream sequencing technologies [9]. RNA-seq is mainly based on the second-generation sequencing platform, with a whole new and complete set of library building, sequencing and analysis systems [10]. Since it was proposed in 2006, RNA sequencing technology has developed rapidly and been widely used. In 2015, it replaced microarray as the main transcription technology [11][12][13]. At present, it has been applied to analyze the cold stress response mechanism of many important cash crops, such as rice, rapeseed, wheat and corn. It is of great signi cance for the in-depth analysis of the molecular mechanism of cold-tolerance and the breeding of cold-tolerance varieties. In addition, it can also provide technical reference and basic information for the study of the stress mechanism of other species.
In the present study, we used the Illumina sequencing platform to sequence the transcriptomes of one cold-tolerant (ZG) and one cold-sensitive (YN) Brassica oleracea L. genotypes. And we compared the transcriptome differences between the two different varieties of Brassica oleracea L. genotypes. The differentially expressed genes (DEGs) were screened and identi ed for the two different varieties of Brassica oleracea L.genotypes. Then pathway enrichment analysis of the DEGs and gene ontology (GO) enrichment analysis were used to indicate the biological processes involved in the response of Brassica oleracea L. to cold stress. The goals of this study were to identify and characterize potential candidate genes involved in cold stress tolerance in Brassica oleracea L.. This should contribute to our understanding of the cold response mechanisms in Brassica oleracea L., and therefore assist efforts to improving cold tolerance in Brassica oleracea L. in future breeding programs.

Plant materials and growth conditions
High cold-tolerance Brassica oleracea L. genotype (Zhonggan1229) and low cold-tolerance Brassica oleracea L. genotype (Yingchun) were provided by Vegetable Breeding-Laboratory of Chinese Academy of Agricultural Sciences. Two varieties were selected to be treated at room temperature of 22°C(control) and under cold stress of 4°C (treatment) for 24h. The two Brassica oleracea L. varieties include the coldtolerant Zhonggan1229 (ZG), the cold-sensitive Yingchun (YN). Among them, the normal temperature treatment and low temperature treatment of YN are expressed by YC and YT, the normal temperature treatment and low temperature treatment of ZG are expressed by ZC and ZT.
The two materials were grown in a 1:1:1 mixture of peat, perlite and vermiculite, Culture and growth were placed in an arti cial light incubator under controlled conditions: 100 soil 20 mol photos, light 100 soil, 16/8 h day/night time, temperature ratio 25°C-20°C day/night, humidity of 60 ± 20%. After 45 d of culture, when the well-grown rape seedlings grew 6-8 true leaves, the Brassica oleracea L. seedlings with the same growth status were selected and treated at a low temperature of 4°C.
Part of the physiological indexes of the two cultivars Relative electrical conductivity To measure relative electrical conductivity (REC), fresh 0.1 g leaves were placed in a coarse test tube, 10 mL double distilled water was added, and soaked in an oscillator at 150 r/min for 12 h. The conductivity (EC 1 ) was measured by Thermo-A212 conductometer after standing. The test tube was put into the boiling water bath for 20 min to kill the plant tissue. After being taken out, it was cooled quickly with tap water and the electrical conductivity (EC 2 ) was determined. The electrical conductivity (EC 0 ) of the double steam water was taken as the control.

Proline content
In order to measure the proline content (Pro), 0.2 g of leaves were placed in the test tube, 5 mL 3% sulfonyl salicylic acid was added, boiling water bath for 10 min, 2 mL extraction solution, 2 mL acetic acid and 2 mL acidic indoline were absorbed, mixed and then boiling water bath for 30 min, cooled, added 4 mL toluene, shaken for 30 s, stood for a while, absorbed the upper proline toluene solution and measured the absorption value at 520 nm.

Soluble sugar content
In order to measure the soluble sugar content, 0.2 g leaves were put into a centrifuge tube, and 15 mL distilled water was added. After being taken out and cooled, the leaves were ltered into a 50 mL volumeter bottle, and the volume of distilled water was constant. 0.25 mL extract, 0.25 mL distilled water and 2.5 mL anthranone reagent were put in the test tube, mixed and then bathed in boiling water for 10min. After cooling, the light absorption value at 620 nm was determined.
RNA extraction,construction and Illumina sequence of cDNA library Sample after the extraction of total RNA in eukaryotes, with Oligo (dT) mRNA magnetic beads enrichment, to join in mRNA from fragmentation buffer makes its fragments short clips, mRNA after fragments as templates, again in six bases random primers (random hexamers) synthetic cDNA rst chain, and add buffer, dNTPs, RNase H and DNA polymerase I synthesize cDNA second chain, After puri cation by QiaQuick PCR kit and elution by EB buffer, terminal repair, alkali addition A, sequencing connector, and then agarose gel electrophoresis to recover the large and small fragments of the target, PCR ampli cation was carried out, and the whole library was completed. Illumina HiSeq TM 2500 was used for sequencing the constructed library.

Identi cation and analysis of DEGs
The mapped reads method for gene expression is used in the calculation of FPKM (Fragments Per Kilobase of transcript), and the calculation formula is as follows: FPKM method can eliminate the in uence of gene length and sequencing quantity difference on the calculated gene expression, and the calculated gene expression quantity can be directly used to compare the gene expression difference between different samples.
To identify DEGs across samples or groups, the edgeR package (http://www.r-project.org/) was used. We identi ed genes with a fold change≥2 and a false discovery rate (FDR)<0.05 in a comparison as signi cant DEGs. DEGs were then subjected to enrichment analysis of GO functions and KEGG pathways.

Functional annotation of differentially expressed genes
Gene Ontology (GO) is an international standardized gene functional classi cation system which offers adynamic-updated controlled vocabulary and a strictly de ned concept to comprehensively describeproperties of genes and their products in any organism. GO has three ontologies: molecular function, cellular component and biological process. The basic unit of GO is GO-term. Each GO-term belongs to a type of ontology, pvalue≤0.05 was taken as the threshold, and the GO term satisfying this condition was de ned as the GO term that was signi cantly enriched in the differentially expressed transcripts. The main biological functions of differentially expressed transcripts can be determined by GO functional signi cance enrichment analysis.
KEGG is the main public database on pathway [62]. Through Pathway signi cant enrichment, the main biochemical metabolic pathways and signal transduction pathways involved in differentially expressed transcripts can be identi ed. Pathways with qvalue≤0.05 were de ned as pathways that were signi cantly enriched in the differentially expressed transcripts.
Veri cation of DEG expression with quantitative real-time PCR In order to verify the results of RNA-seq data, 13 DEGs were randomly selected for quantitative real-time PCR (qRT-PCR), and the results of these genes were consistent with the RNA-seq results (Fig 10). Thirteen potential functional genes were selected for qRT-PCR, among them, the rst four genes were uniquely differentially expressed in ZG in response to low temperature stress, and the last nine were shared differentially expressed genes in response to low temperature stress in the two cultivars. The actin gene of Brassica oleracea L. was used as the internal control in this study. The RNA extracted from the sample leaves was transcribed by cDNA. The cDNA was diluted vefold, and each 2 μl of the diluted cDNA was added into the 18μl PCR mixture. qRT-PCR was performed with 1 min pre-heating at 95°C, followed by 40 cycles of 95°C for 20 s, 55°C for 20 s, and 72°C for 30 s. The relative expression levels were calculated via the 2 -ΔΔCt method. Each sample was analyzed in four biological replicates. The primer sequences used for PCR are listed (Addition le 3: Table S3).

Result
Changes of physiological indexes of the two cultivars during cold stress As can be seen from the picture given( Fig. 1), after 24 hours of low temperature treatment at 4°C, the leaves of ZG plants were slightly withered and yellow, and the stems were also slightly drooping.
Compared with ZG, the leaf wilting and yellowing of YN treated at 4°C were more obvious, and the drooping of stem was more serious. It can be preliminarily inferred that ZG has better cold tolerance than YN. By the bar chart as you can see (Fig. 2), through 4°C low temperature processing, after the two varieties of MDA content and soluble sugar content has the varying degree rise, the MDA content in ZG has risen by 206%, soluble sugar content increased by 166%, the MDA content of YN increased by 104%, soluble sugar content increased by 63%. However, under low temperature stress, the changes of REC and free proline content were different. After the low temperature treatment at 4°C, the REC of ZG plant leaves decreased by 7%, and the content of free proline decreased by 38%. The REC of YN leaves increased by 18% and the content of free proline increased by 33%. That's roughly the same as previous studies [14][15][16].
Transcriptome sequencing results and quality assessment RNA was extracted from the leaves of the plants after two temperature treatments of YN and ZG, and two replicates were performed for each RNA samples. The RNA samples were used to construct cDNA library.
Illumina HiSeq TM 2500 was used for sequencing the constructed library. A total of 390 million clean reads were generated and an overview of the sequencing results is shown in Table 1 and (Addition le 1: Table  S1). The base number with mass value Q≤20 accounted for 98.28%-98.69% of the whole reads, and the base number with mass value Q≤30 accounted for 94.47%-95.66% of the whole reads. In addition, 74.62% clean reads successfully mapped the reference genome.
As can be seen from the Fig. 3, the proportion of Low_quality (the proportion of low quality reads in the total number of reads) in the four experimental groups was around 1%, while the proportion of High_quality_clean (the proportion of the original sequence data obtained after the removal of impurities in the total number of reads) was above 98%. This shows that the data obtained by sequencing is true and reliable.

Distribution of differentially expressed genes (DEGs)
Under low temperature stress, a large number of signi cantly DEGs were found in ZG and YN. There were 3902 signi cantly up-regulated genes and 5309 signi cantly down-regulated genes in ZG under low temperature stress. However, there were 4253 signi cantly up-regulated genes and 5938 signi cantly down-regulated genes in YN under low temperature stress (Fig. 4, Addition le 2: Table S2). Among them, 1844 DEGs are the speci c differential genes in ZG in response to low temperature stress, which is speculated to be the reason for its strong cold resistance. And there were 6,089 DEGs common DEGs in response to low temperature stress in the two varieties, and these genes were presumed to be the genes that responded to low temperature stress commonly found in Brassica oleracea L. varieties (Fig. 5).
In this study, 1844 speci c DEGs in response to low temperature stress in ZG and 6098 common DEGs in response to low temperature stress in ZG and YN were analyzed and annotated. Among the 1844 genes, 791 genes were up-regulated and 1053 genes were down-regulated. Among the 6089 genes, 2611 genes were up-regulated in both cultivars, and 3418 genes were down-regulated in both cultivars. In addition, 38 genes were up-regulated in YN but down-regulated in ZG, and 22 genes were down-regulated in YN but upregulated in ZG.

Annotation of differentially expressed genes
The GO enrichment analysis was performed on the speci c DEGs in ZG in response to low temperature stress. As can be seen from the Fig. 6, 26 of the top 30 enriched GO terms belonged to biological processes, and only 4 terms belonged to molecular functions. The most signi cant enrichment in biological processes, is GO: 0009755 (hormone-mediated signal pathway), there are 61 genes involved in this process, and the most signi cant enrichment in molecular functions is GO: 0016740 (transferase activity), there are 145 genes involved in this process. Presumably the genes that can resist cold stress may be mostly biological processes terms.
The GO enrichment analysis was also performed on the common DEGs in response to low temperature stress in the two varieties. Among the GO enrichment results of these genes, 23 GO terms belonged to biological processes, 1 GO term belonged to molecular functions, and only 6 GO terms belonged to cellular components. The most signi cant enrichment in biological processes is: GO: 0009628 (response to abiotic stress), there are 540 genes involved in this process; The most signi cant enrichment in cellular components is: GO: 0009521 (optical system), there are 33 genes involved in this process; The most signi cant enrichment in molecular functions is: GO: 0016758 (transferase activity, transfer hexosyl), there are 72 genes involved in this process. The above results suggested that most of the universal genes in Brassica oleracea L. resistance to low temperature stress belong to biological processes terms.
After KEGG pathway-analysis of the speci c differential genes in ZG in response to low temperature stress and the common differential genes in response to low temperature stress in the two varieties (Fig.  7). Circadian rhythms of plants and Plant hormone signal transduction was signi cantly enriched in the two analyzed genes, and the effects of low temperature stress were most signi cant. In particular, circadian rhythm pathway is an important pathway found in this experiment (Fig. 8). The expression of genes involved in these two pathways also changed signi cantly in response to low temperature stress( Table 2, Table 3). For example, ARF7 that found in ZG encodes auxin response factor 7-like in the ARF transcription factor family, which is up-regulated under low temperature stress, In addition, ABF3 encoded ABSCISIC ACID-INSENSITIVE 5-like protein 6 in the bZIP transcription factor family, in both ZG and YN, which showed an up-regulated trend in both cultivates under low temperature stress.
Analysis of transcription factor gene family classi cation of differentially expressed genes under low temperature stress In order to analyze the working mechanism of Brassica oleracea L. genes in response to low temperature stress, the selected two DEGs were classi ed and analyzed as transcription factors (Fig. 9). Among the unique differentially expressed genes in ZG in response to low temperature stress, a total of 154 genes were annotated into transcription factor families. Among the 154 genes, 79 genes were up-regulated and 75 genes were down-regulated, the encoding of MYB-related proteins was the largest group, with 15 genes involved, followed by proteins encoding BHLH, ERF, bZIP, C2H2, and DOF, etc. Among the 15 genes, 8 genes were up-regulated and 7 genes were down-regulated.
However, after analyzing the common DEGs that in two varieties. A total of 516 genes were annotated into corresponding transcription factor families. Among the 516 genes, 211 genes were up-regulated and 296 genes were down-regulated. There were also 4 genes that were up-regulated in ZG but downregulated in YN, and 5 genes that were down-regulated in ZG but up-regulated in YN, the largest group was the protein encoding ERF, with 59 genes, followed by proteins encoding bHLH, bZIP, MYB and NAC, etc. However, after analyzing the differentially expressed genes shared by the two varieties in response to low temperature stress, the largest group was the protein encoding ERF, with 59 genes, followed by proteins encoding BHLH, BZIP, MYB, NAC, etc. Among the 59 genes, 32 genes were up-regulated to low temperature stress in both cultivars, 26 genes were down-regulated to low temperature stress in both cultivars, and 1 gene was up-regulated in YN but down-regulated in ZG in response to low temperature stress.

Discussion
Low temperature is an important factor restricting the cultivation of Brassica oleracea L.. At present, there are few studies on the cold tolerance mechanism of Brassica oleracea L. [17][18][19][20], and the differences of cold tolerance among different varieties have not been su ciently analyzed and demonstrated. In this study, two varieties with different cold tolerance were used: low-temperature sensitive type (YN) and lowtemperature tolerant type (ZG). Considering that it is common for Brassica oleracea L. to experience a temperature of about 4°C during winter overwintering cultivation, only one temperature treatment was carried out. In order to verify the results of RNA-seq data, 13 DEGs were randomly selected for qRT-PCR, and the results of these genes were consistent with the RNA-seq results.

Analysis of differentially expressed genes under low temperature stress
In this study, two varieties were treated at 4°C for 24h, and then several representative physiological indexes were determined [21][22][23]. These physiological indices indicated that ZG was more cold tolerant than YN. Through RNA-seq results, there were 9211 DEGs in response to low temperature in ZG, and 10191 DEGs in response to low temperature in YN. In the two varieties contain a large number of genes make response to low temperature. Interestingly, the number of DEGs in response to low temperature stress in YN was slightly higher than that in ZG, Therefore, it can be inferred that it is very important to analyze the cold-resistance mechanism of uniquely DEGs in ZG. All the DEGs of the two varieties were analyzed, and 6098 DEGs were found in both varieties. These genes were speculated to be involved in the common mechanism of response to low temperature stress in Brassica oleracea L.. In addition, 1844 DEGs were only speci cally expressed in ZG, and these genes were speculated to be speci c genes that could improve cold tolerance of Brassica oleracea L..
The GO enrichment of genes from biological processes, cellular components and molecular functions was carried out, and it was found that most of the DEGs, both common to the two varieties and speci c to ZG, were concentrated in terms of biological processes, with little involvement in terms of the other two. In addition, among the top 30 enriched GO terms, 5 terms belonging to biological processes are involved in both genes. These 5 terms are GO: 0009719 (response to endogenous stimulus), GO: 0007165 (signal transduction), GO: 0023052 (signaling), GO: 0044700 (single organism signaling) and GO: 0001101 (response to acid chemical). These ve terms are speculated to be involved in the basic mechanism of drought resistance of Brassica oleracea L. It is worth noting that 4 GO terms belonging to molecular functions only appear in ZG speci c differentially expressed genes, they are respectively: GO: 0016740 (transferase activity), GO: 0019787 (ubiquitin-like protein transferase activity), GO: 0015399 (primary active transmembrane transporter activity) and GO: 0015405 (p-p-bond-bromo-driven transmembrane transporter activity), it can be speculated that these four terms are important functional terms to enhance the cold tolerance of Brassica oleracea L.. Both the same terms with different genes involved and the different terms with the same genes involved can re ect the response degree of varieties to chilling injury [24]. On the other hand, after KEGG pathway enrichment of the two genes, it was found that among the top 20 enrichment pathways, there are two common enrichment pathways: ko04712 (Circadian rhythm-plant) and ko04075 (Plant hormone signal transduction), and all the genes that were commented into the transcription factor family were involved in these two pathways. It is speculated that these two pathways play an important role in the cold resistance mechanism of Brassica oleracea L.. Therefore, It is important to collate the focused pathways and analyze the differentially expressed genes involved in the pathways to explore the mechanism of cold resistance in Brassica oleracea L.. And the terms jointly enriched by these two varieties, the special terms of ZG with strong cold tolerance and the genes involved in these functional terms may be important basis for the improvement of Brassica oleracea L..

Plant hormone signal transduction
Under low temperature stress, plant growth hormone is involved in plant growth and development, signal transduction and other aspects of stress, such as ABA, GA and IAA [25,26]. In this study, a large number of DEGs related to plant hormones appeared under low temperature stress.
In previous studies, the EIN3 is validated and its mutant EIN3-1 in ethylene signal transduction pathway, found three in Arabidopsis thaliana genome and EIN3 higher genetic sequence similarity, respectively named EIL1, EIL2 and EIL3. EIN3 and EIL1, EIL2 can be complementary EIN3-1 phenotype, that they are involved in the ethylene signal transduction [27,28]. However, EIL1 has been proved to be the second transcriptional activator in regulating sulfur de ciency response in Arabidopsis thaliana, second only to SLIM1/EIL3, and EIL3 has also been shown to play an important role in resisting abiotic stress and regulating the expression of ethylene synthesis genes in mulberry, and is up-regulated in response to abiotic stress [29,30]. In this study, EIF2 did not appear, EIF1 and EIF3 did appear in both varieties, but EIF1 showed a down-regulated trend while EIF3 showed an up-regulated trend. Therefore, it can be inferred that EIL3 also plays a positive regulatory role in the response of Brassica oleracea L. varieties to low temperature stress.
In this study, ABF3 and ABF1 appeared in both varieties, and both genes were up-regulated under low temperature stress while involved in regulating hormone signal transduction. The role of ABF1 and ABF3 in stress signal regulation network has been studied in Arabidopsis thaliana, ABF1 was mainly involved in response to low temperature and ABA stress. ABF3 was mainly induced by ABA, high salinity, low temperature, heat and oxidative stress [31][32][33][34].
In addition, ERF2 overexpression has been proved to improve the cold tolerance of tomato and rice, and the expression of ERF15 and ERF16 is regulated by the JA signaling pathway in tomato plants [35,36]. In this experiment, ERF2 and ERF15 appeared in the two varieties. ERF15 showed an up-regulation trend in response to low temperature stress in the two varieties, while ERF2 showed the opposite trend.
In addition to the above contents, some other genes only appeared in ZG under low temperature stress and showed an up-regulated trend: ARF7, TGA1, TGA3, and TGA7. At present, there are few research results on ARF7, and it is basically con rmed that ARF7 is involved in the growth and development of ower buds and pericarp of plants [37,38]. However, in this study, ARF7 was annotated as an auxin response factor. For TGA1, TGA3 and TGA7, it has only been con rmed that they are widely involved in the disease resistance response of plants, and there is no direct connection with plant hormone signal transduction [39][40][41][42].

Circadian rhythm-plant
Photoperiod changes and day-night alternation in nature play an important and complex role in plant growth and development. Studies have shown that circadian rhythm is involved in the response of plants to low temperature stress [43,44]. The effect of biological clocks on plants has also been demonstrated [45].
In down-regulated genes present in two cultivars, only PIF3 was involved in both plant circadian subtitles and plant hormone signal transduction. Recent studies have shown that PIF3 is involved in regulating plant hormone synthesis and response to low temperature stress [46,47]. It is necessary to further study the role of this gene in the cold tolerance mechanism of Brassica oleracea L..
In Arabidopsis thaliana, loss of the central clock genes, CIRCADIAN CLOCK ASSOCIATED 1 (CCA1) and LATE ELON-GATED HYPOCOTYL (LHY), results in both reduced feuptake and photosynthetic efficiency, whereas CCA1 overexpression confers the opposite effects [48,49]. In this study, LHY gene and CCA1 gene appeared in the differential genes in response to low temperature stress, which participated in the process of plant circadian rhythm and belong to the MYB related family. The difference is that LHY plays a role in both cultivars, while CCA1 plays a role in ZG, and both of them showed an up-regulated trend after cold injury. Its suggesting that the chilling tolerance of LHY gene on Brassica oleracea L. is more basic, while CCA1 is more targeted to improve the chilling tolerance of Brassica oleracea L.. Also worth mentioning is COL2 gene, which shows an up-regulated trend in ZG but a down-regulated trend in YN under low temperature stress. In addition, CO and COL1 genes were found in both cultivars and belonged to the same gene family as the above COL2 gene. Studies have con rmed that CO gene is involved in the photoperiodic regulation pathway of plants, and plays an important role in it [50,51]. And the function of COL1 and COL2 genes, which are homologous to the CO gene in Arabidopsis, was also preliminarily identi ed. COL1 gene was proved to shorten the circadian rhythm, but the function of the COL2 gene is not yet understood [52]. CO and COL1 genes are also involved in the circadian rhythm of the two varieties of plants, and both of them show an up-regulated trend under low temperature stress. There have been some studies on CO-like gene family [53,54], and it is speculated that this gene or this gene family may be one of the reasons for the strong cold tolerance of ZG. How COL2 gene and CO-like gene family participate in the cold tolerance mechanism of Brassica oleracea L. can be studied in the future. Now it has been proved that low temperature induced anthocyanin accumulation under light requires HY5 and HYH, and it has also been found that the reduction of endogenous gibberellin is conducive to low temperature induced anthocyanin accumulation. The gibberellin-degrading enzyme encoding gene GA2ox1 was also up-regulated at low temperature and was dependent on HY5/HYH [55][56][57]. In this study, HY5 gene appeared in the circadian rhythm pathways of both varieties. HYH gene and HY5 gene showed an up-regulation trend in response to low temperature stress in the two varieties. However, CDF1, which is only speci c in ZG, was studied in tea tree, non-heading Chinese cabbage, tomato and other crops [58-60], and it was found that CDF1 may interact with FKF1 in non-heading Chinese cabbage, and then participate in the photoperiod owering pathway. However, CDF1 did not play an obvious role in regulating photoperiod owering pathway in tea and tomato. In this study, CDF1 gene only appears in ZG and participates in the circadian rhythm pathway, which also presents an up-regulated trend in response to low temperature stress. So, this gene may be an important gene of ZG with strong cold tolerance.

Conclusion
In this study, Illumina HiSeqTM 02500 was used to carry out ZG and YN cold resistance of different varieties of the transcriptome analysis, found a lot of DEGs, including 1844 DEGs for cold resistance is strong in ZG speci c genes in response to low temperature stress, there were 6089 DEGs of two varieties of genes in response to low temperature stress. GO enrichment and KEGG pathway analysis were performed for these two genes, these differentially expressed genes involved in response to low temperature stress is involved in many biological processes, cell component and molecular function, including participation in circadian rhythm-plant and all kinds of plant hormone synthesis, also associated with a variety of common transcription factors, the gene is associated with the cold resistant mechanism of Brassica oleracea L., need further study and analysis, the results of this study was to Brassica oleracea L. follow-up for cold resistance breeding research provides cold resistance gene resource.

Consent for publication
Not applicable.

Availability of data and materials
Extra data has been appended as supplementary Tables. The accession number for sequence data generated in this study is PRJNA726022 available at https://dataview.ncbi.nlm.nih.gov/object/PRJNA726022?reviewer=vkvnco25irsb6g7pqo3hc657ra.

Competing interests
The authors declare that they have no competing interests.    Figure 1 Morphological changes of ZG(a) and YN(b) seedlings treated at 4°C were compared. Plants on the left side of each gure were normal temperature control(22°C), while plants on the right side were low temperature treatment at 4°C.  Pie chart of the lter reads in ZG under normal temperature treatment. The sample of ZG treated at 22°C(ZC) was taken as an example. The red part represents the proportion of reads containing joint sequences to the total number of reads, the yellow part represents the proportion of reads containing unknown bases to the total number of reads, the green part represents the proportion of low-quality reads to the total number of reads, and the blue part represents the proportion of data obtained after removing impurities from the original sequence data to the total number of reads.

Figure 4
Cluster diagram of differential expression patterns between ZG and YN groups before and after treatments (a). Based on gene expression, hierarchical clustering was performed on the relationship between samples and genes, and the clustering results were presented using heat maps. The gene expression amount of each sample was calculated with base 2 as the logarithm value, and the hierarchical clustering analysis was conducted for different samples and genes. Each column represents a sample, and each row represents a gene. The expression of genes in different samples is shown in different colors. The redder the color means the higher the expression level are, and the bluer the color means the lower the expression level are. Edger software was used to analyze the differences in gene expression between ZG and YN groups before and after treatment (b). FDR and Log2FC were used to screen the differential genes under the screening conditions of FDR<0.05 and |log2FC|>1.
Page 24/30 Figure 5 Venn diagram of genetic statistics of differences between ZG and YN groups before and after treatment.
Four samples were compared in pairs, and different sets represented the number of differentially expressed genes in different parts.

Figure 6
The top 30 Go terms of the speci c DEGs of ZG responds to low temperature stress (a) and the common DEGs of the two cultivars in response to low temperature stress (b). The plot is arranged according to pvalue from small to large. The three colored columns represent the GO terms of biological processes, molecular functions and cellular components, and the length of the columns represents the number of genes enriched in the terms.

Figure 7
The top 20 pathway enrichment of the speci c DEGs of ZG responds to low temperature stress (a) and the common DEGs of the two cultivars in response to low temperature stress (b). The bubble diagram is made by arranging the enrichment pathway qvalue from small to large. The closer the color of the bubble is to red, the smaller the qvalue of the pathway is, and the higher the enrichment degree of genes in the pathway is. The size of the bubble in the pathway re ects the number of genes involved in the pathway, and the larger the bubble is, the more genes involved in the pathway.  Supplementary Files