Transcriptome and Small RNA Combined Sequencing Analysis of Cold Tolerance in Non-heading Chinese Cabbage

Background Non-heading Chinese cabbage (Brassica rapa ssp. chinensis) is an important leaf vegetable grown worldwide. However, there has currently been not enough transcriptome and small RNA combined sequencing analysis of cold tolerance, which hinders further functional genomics research. Results In this study, 63.43 Gb of clean data was obtained from the transcriptome analysis. The clean data of each sample reached 6.99 Gb, and the basic percentage of Q30 was 93.68% and above. The clean reads of each sample were sequence aligned with the designated reference genome (Brassica rapa, IVFCAASv1), and the efficiency of the alignment varied from 81.54 to 87.24%. According to the comparison results, 1,860 new genes were discovered in Pak-choi, of which 1,613 were functionally annotated. Among them, 13 common differentially expressed genes were detected in all materials, including seven upregulated and six downregulated. At the same time, we used quantitative real-time PCR to confirm the changes of these gene expression levels. In addition, we sequenced miRNA of the same material. Our findings revealed a total of 34,182,333 small RNA reads, 88,604,604 kinds of small RNAs, among which the most common size was 24 nt. In all materials, the number of common differential miRNAs is eight. According to the corresponding relationship between miRNA and its target genes, we carried out Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analysis on the set of target genes on each group of differentially expressed miRNAs. Through the analysis, it is found that the distributions of candidate target genes in different materials are different. We not only used transcriptome sequencing and small RNA sequencing but also used experiments to prove the expression levels of differentially expressed genes that were obtained by sequencing. Sequencing combined with experiments proved the mechanism of some differential gene expression levels after low-temperature treatment. Conclusion In all, this study provides a resource for genetic and genomic research under abiotic stress in Pak-choi.


BACKGROUND
Non-heading Chinese cabbage (Brassica rapa ssp. chinensis), also known as Pak-choi, is one of the most common cultivated vegetables in China. Its growth and productivity are often adversely affected by cold and even freezing stresses from the environment (Abe et al., 2003;Yu et al., 2010;Li et al., 2016). The response mechanism of Arabidopsis' cold acclimation and freezing resistance is not exactly the same (Rahman et al., 2020). The cold domestication of plants, which usually involves many biochemical and physiological changes, is complicated and difficult to understand (Salinas, 2002;Kaplan et al., 2004;Chinnusamy et al., 2007;Krasensky and Jonak, 2012). In Arabidopsis, the latest study found that the cold stress response was mediated by the GNOM ARF-GEF pathway (Ashraf and Rahman, 2019). After low-temperature treatment, the metabolism and transcriptome of plants are greatly affected, and the expression levels of certain genes are also regulated. Some related metabolic enzymes are inhibited, and the degree of plant metabolism are affected to some extent (Chinnusamy et al., 2007). Nowadays, genetic, biochemical, and various sequencing methods have been applied to detect how some genes cope with environmental stresses (Chinnusamy et al., 2007;Krasensky and Jonak, 2012), and many plants have been studied today, such as rice Wang et al., 2016), cotton (Zhao et al., 2012), tomato (Starck et al., 2000), potato (Shinozaki and Yamaguchi-Shinozaki, 1996), muskmelons (Wang et al., 2004), and sugarcane (Thakur et al., 2010;Anjum et al., 2011;Zhu et al., 2013). However, the response mechanism of non-heading Chinese cabbage to cold stress remains unclear, and further joint research is needed.
Nowadays, transcriptome analysis has gradually become a useful and general tool for discovering genes in multiple stress pathways, including determining the expression patterns of related genes (Martin and Wang, 2011;Hamilton and Robin Buell, 2012;Ward et al., 2012;Li et al., 2020). The related genes reported in plant secondary metabolism were discovered based on experimental methods of functional genome sequencing (Dixon, 2001;Goossens et al., 2003). In addition, RNA sequencing technology is often used to obtain complete transcriptome information from different plants, such as tea tree, chlorophytum borivilianum, and atractylodes lancea, and provide better insights into transcription or posttranscriptional force, including regulation of essential genes, during secondary metabolite biosynthetic pathways (Kalra et al., 2013;Li et al., 2015;Devi et al., 2016). Nowadays, transcriptome sequencing has been successfully used to detect expression levels of related genes in many organisms, such as rice , yeast (Nagalakshmi et al., 2008;Wilhelm et al., 2008), sweet potato , and taxus (Ge et al., 2011). With Illumina sequencing technology, millions of sequences were read at a time, and individual assembled genes were mapped into a reference transcriptome map for molecular annotation .
Studies have found that microRNA (miRNA), trans short interfering RNA (ta-siRNA), and heterochromatic short interfering RNA (hc-siRNA), all play important roles in different organisms (Axtell, 2013). Among them, miRNA, a type of endogenous small RNA, is composed of about 22 nucleotides (nt) and usually plays a negative role in regulating gene expression (Voinnet, 2009). Many studies have shown that miRNAs are often involved in plant development, hormone signaling, and abiotic stress responses (Jones-Rhoades et al., 2006;Chambers and Shuai, 2009). Generally speaking, small interfering RNAs are processed from perfect double-stranded RNA, while miRNAs are derived from single-stranded RNA transcripts, forming an imperfect double-stranded stem-loop precursor structure (Llave et al., 2002;Khraiwesh et al., 2010;Hao et al., 2012;Szittya and Burgyán, 2013). On the whole, miRNA plays a vital role in various biological and metabolic processes of plant growth and development, such as biotic (or abiotic) stresses, which can also negatively regulate the expression of target genes, by inhibiting (cutting) target mRNA or other ways (Bartel, 2004;Jones-Rhoades and Bartel, 2004;Jones-Rhoades et al., 2006;Mallory and Vaucheret, 2006;Sunkar et al., 2006;Voinnet, 2009;Wu et al., 2010). Therefore, it is important to identify miRNAs and their target genes (or miRNAs), which is essential for a better understanding of miRNA-mediated regulation of cold stress genes.
In this study, we compared the tolerance of two common and typical non-heading Chinese cabbage varieties, Suzhouqing (BcL.1) and Sijiucaixin (BcL.2), to cold stress. We found that BcL.1 is more tolerant to cold stress compared with BcL.2. We used RNA-Seq for comprehensive characterization and explored the effects of low temperature. We identified the most important genes in the low-temperature response and discussed their regulatory networks under cold stress. Furthermore, we identified conserved and novel miRNAs and their potential target genes in non-heading Chinese cabbage, and discussed the possible connections between them. Quantitative real-time polymerase chain reaction (qRT-PCR) was also used to assess the expression levels of common differentially expressed genes (DEGs) and identify those candidate genes involved in cold tolerance. This work might serve as a reference of the functional analysis of cold tolerance in non-heading Chinese cabbage.

Quantity Statistics and Venn Diagram of Differentially Expressed Genes
To study the effects of temperature on plant growth, plants (BcL.1 and BcL.2) were grown for 6 h in environments of 25 and 4 • C. Except for the temperature, the other conditions remain unvaried. 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, plant leaves are more likely to shrink or even wither than under the 25 • C treatment (Figure 1). This result corresponds to previous reports that cold stress usually downgrades the seedling vigor (Kang and Saltveit, 2002;Wang et al., 2016) and causes leaf atrophy, slows crop growth, and ultimately reduces the yield (Cruz and Milach, 2004;Oliver et al., 2007;Ruelland et al., 2009).
Afterward, to study the upregulation and downregulation of common genes shared by each group of treatments, we established a Venn diagram of differentially expressed genes. Between the 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 upregulated, and 308 common genes were downregulated (Figures 2A,B). Meanwhile, between the 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 upregulated, and 117 common genes were downregulated (Figures 2C,D). In all materials, a total of seven common genes were upregulated, and six common genes were downregulated (Figures 2E,F and Supplementary  Table 1). 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 upregulated and 2,569 downregulated genes ( Figure 3C). The annotated unigenes were then assigned to Gene Ontology (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 (Gao et al., 2015). 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 and Supplementary Table 2). To further identify the active biochemical pathways, we mapped it to the reference canonical pathways 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 (Kanehisa and Goto, 2000). A total of 24,199 unigenes were annotated based on a BLASTX search of the KEGG database (Supplementary Table 3 Table 14). Unigenes (3,700) could be classified into 23 COG categories. Out of the 3,700 unigenes, general function prediction only (679, 18.35%) was assigned to 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%) (Supplementary Figure 2A).

Quantitative Real-Time-Polymerase Chain Reaction Validation of the Candidate Differentially Expressed Genes Responsive to Cold Tolerance
To test the reliability of the transcriptome sequencing results, qRT-PCR analysis was used. In this study, 13 common candidate DEGs were selected and detected in all treatments by qRT-PCR analysis (Figure 7). The results of transcriptome sequencing 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, almost all expression levels analyzed by qRT-PCR were highly consistent with the transcriptome sequencing results. These results also confirmed the reliability of transcriptome sequencing data (Figure 7). Through qRT-PCR analysis, it was found that there was only one downregulated gene (Brassica_rapa_new gene_1,153), 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), which were collected, sequenced, and analyzed. Total reads (34,182,333) 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 was the most abundant class among both clean and unique reads in all groups (Figure 8 and Table 1).

Analysis of Known miRNA
In order to obtain the details of the miRNA matched on each sample, the abovementioned reads were mapped to the reference sequence, which are compared with the specified range of sequences in miRBase, including the secondary structure of the known miRNAs on the match, and the information on the sequence, length, and number of occurrences of the miRNA in the present invention. When the miRNA developed into a mature body from the precursor, the process was completed by dicer digestion. The specificity of the cleavage site makes the miRNA mature sequence the first base. There was a strong bias,  so the first base distribution of miRNAs of different lengths was also carried out, in addition to the base distribution statistics of the miRNAs. As shown in Table 2, the number of miRNAs in the C4 group was the highest, 338,668, and the number of miRNAs in the C1 group was the least, 129,932. However, in the C3 group, the types of miRNAs were the most, 1,774; the C1 group had the least types, 1,244. In Figure 9, we listed the secondary structure of the 10 known miRNAs on the match.

Predicted New miRNA
The signature hairpin structure of miRNA precursors can be used to predict new miRNAs. As shown in Table 3, the number of miRNAs was the highest in the C4 group, 103,663, and the number of miRNAs was the least in the C1 group, 60,344. However, in the C3 group, the types of miRNAs were the most, 2,985; the C1 group had the least types, 2,318. We listed the secondary structure of the 10 predicted new miRNA on the match (Figure 10).

Screening and Identification of Differential miRNAs
The correlation analysis of gene expression levels between samples was carried out to test the reliability of the experimental results and the rationality of sample selection. In Figure 11C, R 2 , the square of the Pearson correlation coefficient, was basically at 0.772-1, indicating that the similarity of expression patterns between samples is higher. In Figure 11, the TPM density distributions of miRNA under different experimental conditions were compared. Finally, by using volcano plots, we inferred the overall distribution of differential miRNA. Differential miRNAs were screened based on fold changes in levels and corrected significance levels (padj/q value). In the C2 and C1 groups, 20 differential miRNAs were upregulated, and 16 differential miRNAs were downregulated. In the C4 and C3 groups, 31 differential miRNAs were upregulated, and 47 differential miRNAs were downregulated. In the C3 and C1 groups, 44 differential miRNAs were upregulated, and 33 differential miRNAs were downregulated. In the C4 and C2 groups, 37 differential miRNAs were upregulated and 47 differential miRNAs were downregulated (Figure 12).

Cluster Analysis of Differential miRNAs
The clustering pattern of differential miRNA expression under different experimental conditions was determined by using differential miRNA cluster analysis. For each comparison combination, a set of differential miRNAs was obtained and was 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 was the highest in the C2 group, while C4 had the lowest. In the C4 group, while some miRNAs had comparatively lower expression levels, others, such as bra-miR9557-3p and bra-miR9557-5p, were expressed at very high levels (Figure 13).

Venn Diagram of Differential miRNAs
Next, we show more intuitively the common and unique differences of each comparison combination. When the number of miRNAs was greater than or equal to two and less than or equal to five, the number of differential miRNAs, which was obtained by comparison in each group, can be counted and plotted as a Venn diagram (Figure 14). In Figure 14A, there were 17 common differential miRNAs; while in Figure 14B, there were 38 common differential miRNAs. In all combinations, eight common differential miRNAs are shown in Figure 14C ( Table 4 and Supplementary Table 18).

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. Through the 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 and Supplementary Table 19). 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 (Supplementary Figure 4A and Supplementary Table 21). 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) were  Figure 5A and Supplementary Table 23). Finally, 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 (Supplementary Figure 6A and Supplementary Table 25).

the most representative GO terms (Supplementary
Through the KEGG pathway enrichment analysis, between the C2 and C1 groups, the starch and sucrose metabolism pathway was the largest, with 42 genes, followed by the plantpathogen interaction pathway, with 36 genes (Figure 15B and Supplementary Table 20). Between the C4 and C3 groups, the starch and sucrose metabolism pathway, was the largest, with 74 genes, followed by the plant-pathogen interaction pathway, with 67 genes (Supplementary Figure 4B and Supplementary  Table 22). Between the C3 and C1 groups, the plant-pathogen interaction pathway, with 68 genes, and the RNA transport pathway, also containing 68 genes, are both the pathways with the most genes (Supplementary Figure 5B and Supplementary  Table 24). Between the C4 and C2 groups, the plant-pathogen interaction pathway was the largest, with 86 genes, followed by the starch and sucrose metabolism pathway, with 78 genes (Supplementary Figure 6B and Supplementary Table 26).
FIGURE 13 | Cluster analysis of differentially expressed miRNA. The above figure is the overall hierarchical clustering chart. Clustering is performed by the log 10 (TPM + 1) value. Red indicates high expression miRNA, and blue indicates low expression miRNA.

DISCUSSION
As a convenient tool for transcriptome analysis, RNA-Seq, has been increasingly targeted in species that do not have access to genomic sequences Liu et al., 2013;Soetaert et al., 2013). Here, we used RNA-Seq to measure the gene expression levels in non-heading Chinese cabbage of different varieties and under low-temperature treatment. To obtain all the DEGs from RNA-Seq data, the expression of all genes was analyzed depending on the RPKM. By comparing the data from each groups, we found that between the BcL.1-25 and BcL.1-4 groups, the most differentially expressed genes (6,208 DEGs) were enriched, and their number was almost six times that of the other groups, and the difference was the greatest. At the same time, the ribosome pathway is worth mentioning, and the most abundant genes (287 genes) were also shown, far more than those of the other groups. The results showed that cold stress could affect the genes involved in the expression of these pathways. Previous reports have also detected these rich pathways, partially reflecting the credibility of our results (Kreps et al., 2002;Pang et al., 2013). This also suggests that the ribosomal pathway might be involved in cold stress.
RNA-Seq was searched for some low temperature-related differentially expressed genes. qRT-PCR was used to identify gene expression levels (Taylor et al., 2010). To confirm the reliability of the RNA-Seq results, we performed qRT-PCR experiment to identify them. According to the qRT-PCR results, the expression patterns of all unigenes were consistent with the transcriptome sequencing data, showing that our experimental results were reliable.
Small RNAs are short, non-coding RNAs, usually 19-25 nt in length, and two protruding sizes of 21 and 24 nt, respectively (Kim et al., 2012). In general, an miRNA corresponds to a 21 nt class of small molecule RNA. Research have also found that small RNAs showed a wide range of functions, including heterochromatin formation, gene silencing, and DNA methylation (Lippman and Martienssen, 2004;Penterman et al., 2007). By small RNA sequencing, the result showed that a 24-nt length is the most abundant category, among pure and unique reads in all groups (Figure 8). Our result was highly consistent with previous studies on A. thaliana (Rajagopalan et al., 2006), Oriza sativa , Medicago truncatula (Szittya et al., 2008), and Populus trichocarpa (Puzey et al., 2012).
Furthermore, the known miRNAs are analyzed, and new miRNAs are predicted. We discovered that these selected known miRNAs have at least two stem-loop structures that are obtained by self-folding (Figures 9A-J), suggesting that there might be protein or chromosomal binding sites. Meanwhile, these predicted miRNAs have at least two stem-loop structures that are obtained by self-folding (Figures 10A-J), suggesting that there might also be protein or chromosomal binding sites. The functions of these miRNAs and their target genes were comprehensively analyzed and might provide new insights into miRNA-mediated epigenetic control of Pak-choi under low temperature.
To better understand the functional roles of these predicted miRNA targets, the target genes, which were functionally FIGURE 14 | Venn diagram of differential miRNA. The large circles represent each comparison combinations, and the sum of the number in each large circle represents the total number of differential miRNAs from the comparison combinations, and the overlapping portions of the circles represent the number of common differential miRNAs between the combinations. (A) Venn diagram of differential miRNA between C4 vs C3 and C2 vs C1 groups. (B) Venn diagram of differential miRNA between C4 vs C2 and C3 vs C1 groups.  annotated in biological processes, and the KEGG pathway was also used to describe the corresponding metabolic pathways. After low-temperature treatment, most target genes were enriched in the metabolic pathways of starch and sucrose, and was basically consistent with the results of the abovementioned transcriptome study. This indicated that a series of biosynthetic and metabolic pathways might be induced after low-temperature treatment (Hu et al., 2016). Here, a combination of transcriptome and small RNA sequencing was used to analyze the cold tolerance of non-heading Chinese cabbage. This study might promote further molecular regulation mechanisms of cold tolerance in Pak-choi.

CONCLUSION
In this study, a total of 63.43 Gb clean data were obtained from the transcriptome analysis. Based on the comparison results, a total of 1,860 new genes were discovered in Pak-choi, and 13 common DEGs were detected in all treatments, including seven upregulated and six downregulated. Some of the DEGs were involved in response to abiotic stresses, such as freezing, cold, and salt stresses. Among them, the cold stress response is more obvious, so we used qRT-PCR experiment to confirm changes in the expression levels of these genes. Furthermore, we performed miRNA sequencing analysis on the same material. We found that the results revealed a total of 34,182,333 small RNA reads and a total of 88,604,604 kinds of small RNA, among which the most common size was 24 nt. In all materials, the number of common differential miRNAs was eight. The number of known mature miRNAs and the number of precursors were 110 and 80, respectively. The number of predicted novel miRNA matures, and the numbers of precursors were 75 and 84, respectively. According to the GO and KEGG enrichment analysis, singleorganism cellular process in the biological process, membrane in the cellular component, and protein binding in the molecular function were almost all the most representative level GO terms in all data sets, while the starch and sucrose metabolism pathway was almost all the largest, followed by the plant-pathogen interaction pathway. Our findings highlighted the significance of cold signaling in Pak-choi and might provide a foundation for subsequent research under abiotic stress in the future.

Plant Growth Environment and Treatment Conditions
The materials used in this study were non-heading Chinese cabbage varieties, Suzhouqing (BcL.1) and Sijiucaixin (BcL.

Transcriptome Sequencing
Following the manufacturer's instructions, we used RNAiso Plus reagent (TaKaRa, Dalian, China) to extract total RNA from the samples. RNA samples were examined by using a spectrophotometer and electrophoresed on a 1% agarose gel. cDNA library construction and transcriptome sequencing were performed by Biomarker Technologies (Beijing, China). Transcriptome sequencing was done by Beijing Biomarker Biotechnology 1 Co., Ltd. By Illumina HiSeq X 10 sequencing platform and Pe150 mode sequencing, all clean reads were subsequently mapped to the Brassica rapa reference genome sequence (IVFCAASv1) 2 . The clean reads of each sample sequence were aligned with the designated reference genome. Gene expression analysis was performed based on the comparison results; differentially expressed genes were identified based on their expression levels in different samples, and their functional annotation and enrichment analysis were also performed (Wang et al., 2017).

Transcriptome Assembly and Functional Annotation
The raw data of the transcriptome sequencing were purified by trimming adapters, removing reads containing poly-N, and rejecting the low-quality data (quality value ≤ 10 or unknown nucleotides larger than 5%) to get the clean reads. Meanwhile, the proportion of nucleotides with quality values greater than 30 (Q30) and GC content of the clean data were calculated. Then, all of the clean reads were assembled, using the Trinity program (Grabherr et al., 2011). First, the certain short reads with overlap regions were assembled into longer contiguous sequences for each library. Then, the distance of different contigs was recognized, mapping the clean reads, based on the paired-end information, to obtain the sequence of the transcripts. Finally, the unigenes were obtained, performing the sequence of potential transcript to the TGI Clustering tool (Pertea et al., 2003). The new genes discovered were performed with NR (Deng et al., 2006), Swiss-Prot (Apweiler et al., 2004), GO (Ashburner et al., 2000), COG (Tatusov et al., 2000), KOG (Koonin et al., 2004), Pfam (Finn et al., 2013), and KEGG (Kanehisa et al., 2004) databases, using BLAST (Altschul et al., 1997) software. KOBAS2.0  was used to obtain the KEGG orthology result of the new gene for sequence alignment. After predicting the amino acid sequence of the new gene, we used the HMMER (Eddy, 1998) software to align with the Pfam database and obtained the annotation information of the new gene.

Differentially Expressed Gene Analysis
Following fragments per kilobase of exon per million fragments mapped reads (FPKM) method, the expression level of unigene was calculated (Mortazavi et al., 2008). The ratio of the FPKM values (using 0.001 instead of 0 if the FPKM was 0) was taken as the fold changes in the expression of each gene to identify DEGs between each groups. The false discovery rate (FDR) control method was used to identify the threshold of the p-value in multiple tests and to compute the significance of the difference in transcript abundance (Reiner et al., 2003). In this result, only fold change with | log 2 (case_FPKM/control_FPKM)| ≥ 1, and an FDR ≤ 0.001 were taken as the threshold for significantly differential expression. The log 2 -transformed FPKM value for DEGs was applied to generate heat map by MeV 4.7 (Howe et al., 2011). Meanwhile, the DEGs were annotated with the GO and KEGG databases.

Validation of Differentially Expressed Genes With Quantitative Real-Time Polymerase Chain Reaction
qRT-PCR was used to confirm the expression of common differentially expressed genes in Pak-choi. Total RNA was extracted from each sample. The first-strand cDNA was synthesized, using a PrimeScript TM II First Strand cDNA synthesis kit (TaKaRa Bio, Dalian, China), according to the manufacturer's protocol. The primers were designed, using the software Beacon Designer 7.9, and listed in Table 5. The quantified expression levels of the tested genes were normalized against the housekeeping genes Cyclophilin 1 (CYP1) . The qRT-PCR assays were performed with three biological and technical replicates. Each reaction was performed in 20-µl reaction mixtures containing a diluted cDNA sample as template, SYBR Premix Ex Taq (2×) (TaKaRa, Kyoto, Japan) and genespecific primers. Conditions for quantitative analysis were as follows: 95 • C for 3 min, then 40 cycles (95 • C 30 s, 60 • C 30 s), and 72 • C for 30 s. qRT-PCR was performed according to a previous report (Song et al., 2016). The comparative Ct value method was adopted to analyze the relative gene expression according to a previous analysis and RNA expression levels relative to the actin gene were calculated as 2 − CT (Pfaffl, 2001;Song et al., 2016).

Small RNA Sequencing
By high-throughput sequencing (such as Illumina HiSeqTM2500/MiSeq and other sequencing platforms), sequenced raw image data files were converted into sequenced reads by base calling analysis. The raw data, containing the sequence information of read sequences and the corresponding sequencing quality information, were stored in FASTQ (abbreviated as fq) file format. Referring to the standard definition of miRNA (Allen et al., 2005;Schwab et al., 2005), the candidate target gene of miRNA was compared as the query sequence with the Brassica rapa database 3 . The control samples and inoculation samples were mixed for small RNA library construction, respectively. According to the reported procedures, the construction of small RNA libraries were completed (Sunkar et al., 2008). Three micrograms of total RNA per sample was used as input material for the small RNA library (Huang et al., 2018). Small RNA library construction and small RNA deep sequencing proceeded following the detailed protocol provided by the genome sequencing company (Novogene, China).

Bioinformatic Analysis of Sequence Data
The raw data were first processed through custom Perl and Python scripts. The clean data were mapped to the reference sequence in miRBase21.0 by Bowtie (Langmead et al., 2009), without mismatch to look for known miRNAs. Then, the other reads were integrated to predict novel miRNAs using the available miREvo (Wen et al., 2010) and miRDeep2 (Friedländer et al., 2011) software. The miRNA counts as well as base bias were identified by using custom scripts. Then, the miFam.dat 4 was used to look for families of known miRNAs. The novel miRNA precursor was submitted to Rfam 5 to look for Rfam families.

Venn Diagrams of Known miRNAs and Novel miRNAs
Normalization formula (normalized expression = mapped read count/total reads * 1,000,000) was used to estimate miRNA expression levels (Zhou et al., 2010). DEG seq R package was used to analyze the differential expression of two samples with the criterion of Q < 0.01 and | log 2 (fold change)| > 1 (Storey, 2003).

Construction of Degradome Libraries
Target genes of candidate miRNAs were verified by degradome sequencing by using total RNA, the same as the RNA used for small RNA sequencing library construction, following the published parallel analysis of RNA End (PARE) protocol (German et al., 2009). The data analysis was processed, following the procedure instructions (Novogene, China).

Target Gene Prediction and Annotation for Known and Novel miRNAs
The psRobot_tarin psRobot was performed to predict target genes of miRNA (Wu et al., 2012). To further explore the detailed molecular mechanism of miRNAs in Pak-choi response to cold stress, the target transcripts of differentially expressed miRNAs were analyzed by GO and KEGG functional annotation suites. Subsequently, the Revigo tool 6 was implemented for enrichment analysis of the target genes. The KOBAS software was used to test the statistical enrichment of the target gene candidates in the KEGG pathways (Mao et al., 2005).

DATA AVAILABILITY STATEMENT
The full data sets have been submitted to NCBI Sequence Read Archive (SRA) under Accession SUB8670775, Bioproject: PRJNA682276.

AUTHOR CONTRIBUTIONS
JW completed the relevant experiments and wrote the manuscript. QZ, XY, and XH modified and approved the manuscript. XY and XH interpreted the results and coordinated the research. All authors have read and finalized the draft, agreeing that the manuscript will be submitted to this prestigious journal.