Quantity statistics and veen diagram of differentially expressed genes
To study the effects of temperature on plant growth, plants (BcL.1 and BcL.2) were grown for up to 6 hours in environments of 25 °C and 4 °C. Except for materials and temperature, the other conditions remain the same. Then, we observed that low temperatures have an important effect on plant phenotype. Whether it is BcL.1 or BcL.2, under 4 °C treatment, leaves of plant are more likely to shrink or even wither than 25 °C treatment (Figure 1). This result corresponds to previous reports that cold stress usually lowers the seedling vigor [12, 47], leaf atrophy, slows crop growth and ultimately reduces yield [48-50].
Afterward, to study the up-regulation and down-regulation of common genes shared by each group of treatments, we established a venn diagram of differentially expressed genes. Between G0 (BcL.2-25 vs BcL.2-4) and G2 (BcL.1-25 vs BcL.1-4) groups, a total of 313 common genes were up-regulated and 308 common genes were down-regulated ( Figure 2A and 2B); while between G1 (BcL.1-25 vs BcL.2-25) and G3 (BcL.1-4 vs BcL.2-4) groups, a total of 344 common genes were found to be up-regulated, 117 common genes were down-regulated (Figure 2C and 2D); in all materials, a total of 7 common genes were up-regulated and 6 common genes were down-regulated (Figure 2E and 2F; Table S1). We speculated that these DEGs might help increase the potential application value of non-heading Chinese cabbage under cold stress.
Functional annotation and classification
Between the BcL.1-25 and BcL.1-4 groups, 6,208 DEGs (p < 0.05) were detected, including 3,639 up-regulated and 2,569 down-regulated genes (Figure 3C). The annotated unigenes were then assigned to Gene Othology (GO) terms for functional classification. Three main categories (biological process, molecular function, and cellular component) of GO classification were analyzed separately to investigate their functional distribution. To simplify the functional distribution of plants, we assigned the annotated sequences to GO-slim terms to obtain a “thin” version of classification [51]. Cellular process (GO:0009987, 3,589 genes) and metabolic process (GO:0008152, 3,424 genes) in the biological process, cell part (GO:0004464, 5,170 genes) and cells (GO:0005623, 5,169 genes) in the cellular component and binding activity (GO:0005488, 2,808 genes) and catalytic activity (GO:0003824, 2,279 genes) in the molecular function were the most representative level 2 GO terms in all three data sets, respectively (Figure 3A; Table S2). To further identify the active biochemical pathways, we mapped it to the reference canonical pathways, which in the Kyoto Encyclopedia of Genes and Genomes (KEGG). KEGG is thought to provide a basic platform for systematic analysis of gene function in terms of the network of gene products [52]. A total of 24,199 unigenes were annotated based on a BLASTX search of the KEGG database (Table S3): 263 biosynthesis pathways were predicted and classified into five categories. Of which, the ribosome pathway was the largest, containing 287 genes (287 out of 1,586, 18.10%) (Figure 3B; Figure S1A; Table S10). The annotated unigenes were categorized into different functional groups based on the Cluster of Orthologus Groups (COG) database (Table S14). 3,700 unigenes could be classified into 23 COG categories. Out of 3,700 unigenes, general function prediction only (679, 18.35%) was assigned into the COG category of general function prediction, which represented the largest functional group of the 23 COG categories, followed by translation, ribosomal structure and biogenesis (433, 11.70%), transcription (329, 8.90%), replication, recombination and repair (293, 7.92%) and signal transduction mechanisms (290, 7.84%) (Figure S2A).
Between the BcL.2-25 and BcL.2-4 groups, 964 DEGs (p < 0.05) were detected, including 437 up-regulated and 527 down-regulated genes (Figure 4C). Through GO enrichment stratification analysis, cellular process (GO:0009987, 552 genes) and metabolic process (GO:0008152, 534 genes) in the biological process, cell part (GO:0004464, 741 genes) and cells (GO:0005623, 741 genes) in the cellular component and binding activity (GO:0005488, 429 genes) and catalytic activity (GO:0003824, 433 genes) in the molecular function, which were the most representative level 2 GO terms, respectively (Figure 4A; Table S4). Through KEGG pathway enrichment analysis (Table S5), the protein processing in endoplasmic reticulum pathway was the largest, containing 26 genes (26 out of 251, 10.36%) (Figure 4B; Figure S1B; Table S11). Through COG classification of differentially expressed gene. Out of 554 unigenes (Table S15), general function prediction only (109, 19.68%) was assigned into the COG category of general function prediction, which represented the largest functional group of the 21 COG categories, followed by posttranslational modification, protein turnover, chaperones (59, 10.65%), amino acid transport and metabolism (49, 8.84%), carbohydrate transport and metabolism (45, 8.12%) and signal transduction mechanisms (35, 6.32%) (Figure S2B).
Between the BcL.1-25 and BcL.2-25 groups, 1,448 DEGs (p < 0.05) were detected, including 1,052 up-regulated and 396 down-regulated genes (Figure 5C). Through GO enrichment stratification analysis, cellular process (GO:0009987, 729 genes) and metabolic process (GO:0008152, 695 genes) in the biological process, cell part (GO:0004464, 988 genes) and cells (GO:0005623, 989 genes) in the cellular component and binding activity (GO:0005488, 556 genes) and catalytic activity (GO:0003824, 465 genes) in the molecular function were the most representative level 2 GO terms, respectively (Figure 5A; Table S6). Through KEGG pathway enrichment analysis (Table S7 ), the ribosome pathway was the largest, containing 96 genes (96 out of 366, 26.23%) (Figure 5B; Figure S1C; Table S12). Through COG classification of differentially expressed gene. Out of 681 unigenes (Table S16), general function prediction only (124, 18.21%) was assigned into the COG category of general function prediction, which represented the largest functional group of the 23 COG categories, followed by translation, ribosomal structure and biogenesis (109, 16.01%), transcription (51, 7.49%), replication, recombination and repair (50, 7.34%) and amino acid transport and metabolism (44, 6.46%) (Figure S2C).
Between the BcL.1-4 and BcL.2-4 groups, 1,479 DEGs (p < 0.05) were detected, including 854 up-regulated and 625 down-regulated genes (Figure 6C). Through GO enrichment stratification analysis, cellular process (GO:0009987, 750 genes) and metabolic process (GO:0008152, 697 genes) in the biological process, cell part (GO:0004464, 1,056 genes) and cells (GO:0005623, 1,057 genes) in the cellular component and binding (GO:0005488, 610 genes) and catalytic activity (GO:0003824, 536 genes) in the molecular function, which were the most representative level 2 GO terms, respectively (Figure 6A; Table S8). Through KEGG pathway enrichment analysis (Table S9) , the DNA replication pathway was the largest, containing 13 genes (13 out of 329, 3.95%) (Figure 6B; Figure S1D; Table S13). Through COG classification of differentially expressed gene. Out of 681 unigenes (Table S17), general function prediction only (144, 17.98%) was assigned into the COG category of general function prediction, which represented the largest functional group of the 23 COG categories, followed by carbohydrate transport and metabolism (74, 9.24%), posttranslational modification, protein turnover, chaperones (67, 8.36%), replication, recombination and repair (65, 8.11%) and amino acid transport and metabolism (58, 7.24%) (Figure S2D).
Clustering and functional enrichment of DEGs in all treatments
Among them, 13 DEGs were detected in all treatments, including 7 up-regulated and 6 down-regulated. Some of the DEGs, which were involved in response to stress (GO:0006950) and stimulus (GO:0050896), as well as response to abiotic stress (GO:0009628), such as freezing (GO:0050826), cold (GO:0009409) and salt (GO:0009651) stress. Some of the DEGs response to growth hormone (GO:0060416) and water deprivation (GO:0009414) (Table 1S). In addition, we performed cluster analysis on all screened differentially expressed genes (Figure 3C; Figure 4C; Figure 5C; Figure 6C). Nearly all differentially expressed genes are up-regulated between BcL.1-25 (N-25) and BcL.1-4 (N-4) groups; between BcL.2-25 (B-25) and BcL.2-4 (B-4) groups, most of the genes were up-regulated, while the BcL.2-25 (B-25) group showed more up-regulated genes; between BcL.1-25 (N-25) and BcL.2-25 (B-25) groups, the performance of up-regulated genes and down-regulated genes, which is very similar to that between BcL.1-4 (N-4) and BcL.2-4 (B-4) groups (Figure S3).
qRT-PCR validation of the candidate DEGs responsive to cold tolarence
To test the reliability of the transcriptome sequencing results, qRT-PCR analysis, which were used to perform. In this study, 13 common candidate DEGs, which were selected and detected in all treatments by qRT-PCR analysis (Figure 7). The results of transcriptome sequencing, which were compared with the results of qRT-PCR experiments. Our results showed that even if the fold changes in the expression levels of certain genes detected by transcriptome sequencing and qRT-PCR analysis did not match, but almost all expression levels analyzed by qRT-PCR, which were highly consistent with the transcriptome sequencing results. These results, which also confirmed the reliability of transcriptome sequencing data (Figure 7). Through qRT-PCR analysis, it was found that there was only one down-regulated gene (Brassica_rapa_newGene_1153), and its expression level was different from the RNA-Seq data (Figure 7).
Overview of small RNA sequencing data
In this study, these samples included C1 (BcL.1-25), C2 (BcL.1-4), C3 (BcL.2-25) and C4 (BcL.2-4) were collected, sequenced and analyzed. 34,182,333 total reads were generated, and 8,836,042 unique reads were isolated. After removing low-quality reads, the length distribution of the small RNAs (18-35 nt) revealed that a length of 24 nt, which was the most abundant class among both clean and unique reads in all groups (Figure 8; Table 2).
Analysis of known miRNA
To obtain the details of the miRNA matched on each sample, the above-mentioned reads, which mapped to the reference sequence are compared with the specified range of sequences in miRBase, including the secondary structure of the known miRNAs on the match. Information on the sequence, length, and number of occurrences of the miRNA in the present invention. When the miRNA develops into a mature body from the precursor, the process is completed by dicer digestion. The specificity of the cleavage site makes the miRNA mature sequence sequence the first base. There is a strong bias, so the first base distribution of miRNAs of different lengths is also carried out, in addition to the base distribution statistics of the miRNAs. As shown in Table 3, the number of miRNAs in the C4 group, which was the highest, 338,668; the number of miRNAs in the C1 group, which was the least, 129,932. However, the C3 group had the largest number of miRNAs in the 1,774 species; the C1 group had the least miRNA species, which was 1,244 species. We listed the secondary structure of the 10 known miRNAs on the match (Figure 9).
Predicted new miRNA
The signature hairpin structure of miRNA precursors, which can be used to predict new miRNAs. As shown in Table 4, the number of miRNAs was the highest in the C4 group, 103,663; the number of miRNAs was the least in the C1 group, 60,344. However, the C3 group had the largest number of miRNAs in the 2,985 species. The C1 group had the least miRNA species, which was 2,318 species. We listed the secondary structure of the 10 predicted new miRNA on the match (Figure 10).
Screening and identification of differential miRNAs
To test the reliability of the experimental results and the rationality of sample selection, the correlation analysis of gene expression levels between samples, which was carried out. R2: The square of the pearson correlation coefficient is basically at 0.772-1 (Figure 11C), indicating that the similarity of expression patterns between samples is higher. By comparing the TPM density distribution of miRNA under different experimental conditions, so that the TPM distribution between different experimental conditions, which can be checked as a whole (Figure 11A; Figure 11B). Finally, by using volcano plots, to infer the overall distribution of differential miRNAs. Differential miRNAs, which were screened based on fold changes in levels and corrected significance levels (padj / qvalue). In the C2 and C1 groups, 20 differentially miRNAs, which were up-regulated and 16 differentially miRNAs were down-regulated; in the C4 and C3 groups, 31 differentially miRNA, whichs were up-regulated and 47 differentially miRNAs were down-regulated; in the C3 and C1 groups, 44 differentially miRNAs, which were up-regulated and 33 differentially miRNAs were down-regulated; in the C4 and C2 groups, 37 differentially miRNAs, which were up-regulated and 47 differentially miRNAs were down-regulated (Figure 12).
Cluster analysis of differential miRNAs
The clustering pattern of differential miRNA expression under different experimental conditions, which was determined by using differential miRNA cluster analysis. For each comparison combination, a set of differential miRNAs, which is obtained and will be used for hierarchical clustering analysis. The number of miRNAs with high expression levels of C1, C2, and C3 was higher than that of the C4 group. In addition, the number of highly expressed miRNAs, which was highest in the C2 group, while C4 was the smallest. Some miRNAs in the C4 group have lower expression levels in other groups, but the highest expression in the C4 group, such as bra-miR9557-3p, bra-miR9557-5p, and so on (Figure 13).
Venn diagram of differential miRNAs
In order to more intuitively show the common and unique differences of each comparison combination. When the number of miRNAs, which is greater than or equal to 2 and less than or equal to 5, the number of differential miRNAs, which obtained by comparison in each group can be counted and plotted as a venn diagram (Figure 14). In Figure 14A, there are 17 common differential miRNAs; while in Figure 14B, there are 38 common differential miRNAs. In all combinations, 8 common differential miRNAs are shown in Figure 14C (Table 5; Table S18).
Enrichment analysis of differential miRNA candidate target genes
After obtaining the differentially expressed miRNAs between the groups, according to the correspondence between the miRNA and its target genes, we performed GO and KEGG enrichment analysis on the set of target genes of each group of differentially expressed miRNAs.
GO enrichment stratification analysis, between the C2 and C1 groups, single-organism cellular process (GO:0044763, 896 genes), membrane (GO:0016020, 747 genes), and protein binding (GO:0005515, 994 genes) were the most representative GO terms in biological process, cellular component, and molecular function, respectively (Figure 15A; Table S19); between the C4 and C3 groups, single-organism cellular process (GO:0044763, 1,748 genes), membrane (GO:0016020, 1,431 genes), and protein binding (GO:0005515, 1,973 genes) were the most representative GO terms (Figure S4A; Table S21); between the C3 and C1 groups, single-organism cellular process (GO:0044763, 1,841 genes), membrane (GO:0016020, 1,511 genes), and protein binding (GO:0005515, 2,041 genes), which were the most representative GO terms (Figure S5A; Table S23); between the C4 and C2 groups, single-organism cellular process (GO:0044763, 1,880 genes), membrane (GO:0016020, 1,515 genes), and ion binding (GO:0005515, 2,198 genes) were the most representative GO terms (Figure S6A; Table S25).
Through KEGG pathway enrichment analysis, between the C2 and C1 groups, the starch and sucrose metabolism pathway, which was the largest, 42 genes, followed by the plant-pathogen interaction pathway, 36 genes (Figure 15B; Table S20); between the C4 and C3 groups, the starch and sucrose metabolism pathway, which was the largest, 74 genes, followed by the plant-pathogen interaction pathway, 67 genes (Figure S4B; Table S22); between the C3 and C1 groups, the plant-pathogen interaction pathway, 68 genes, RNA transport pathway, also containing 68 genes, they are both the pathway with the most genes (Figure S5B; Table S24); between the C4 and C2 groups, the plant-pathogen interaction pathway, which was the largest, 86 genes, followed by the starch and sucrose metabolism pathway, 78 genes (Figure S6B; Table S26).