Transcriptome and small RNA combined sequencing analysis of cold tolerance response in non-heading Chinese cabbage

Background: Non-heading Chinese cabbage ( Brassica rapa ssp. chinensis ), as an important leaf vegetable grown worldwide. However, there is currently no 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, of which 1,613 were functionally annotated. Among them, 13 common differentially expressed genes were detected in all materials, including 7 up-regulated and 6 down-regulated. 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 RNA, among which the most common size was 24 nt. In all materials, the number of common differential miRNAs is 8. 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 that each group of differentially expressed miRNAs. Through the analysis, it is found that the distribution of candidate target genes in different materials is different. We not only use transcriptome sequencing and small RNA sequencing, but also use experiment to prove that the expression level of differentially expressed genes are obtained by sequencing. Sequencing combined with experiments proved the mechanism of some differential gene expression levels after low temperature treatment. COG: Cluster of Orthologus Groups; CYP1: Cyclophilin 1; DEGs: Differentially expressed genes; FC: Fold change; FDR: False discovery rate; GO: Gene ontology; hc-siRNA: heterochromatic short interfering RNA; KEGG: Kyoto Encyclopedia of Genes and Genomes; miRNA: microRNA; nt: nucleotides; qRT-PCR: quantitative real-time polymerase chain reaction; RNA-Seq: RNA sequencing; ta-siRNA: trans short interfering RNA


Background
One of the most common cultivated vegetables in China, non-heading Chinese cabbage (Brassica rapa ssp. chinensis), also known as Pak-choi, while its growth are often affected, from the surroundings by cold even freezing stresses [1,2]. Low temperature, well known, as an important factor, affecting growth and reducing productivity [3][4][5]. However, the cold domestication of plants, which usually involves many biochemical and physiological changes, which is complicated and difficult to understand [4][5][6][7]. The metabolism and transcriptome of plants, which can be greatly affected by cold stress, after low temperature treatment, the expression levels of certain genes will also be regulated, as we all know, some related metabolic enzymes will be inhibited, and the degree of plant metabolism will be affected to some extent [4]. Nowadays, genetic, biochemical, and various sequencing methods, which have been applied to respond to many genes that can cope with corresponding environmental stress [8][9][10]. To better understand the response mechanisms of cold stress, many plants have been studied today, such as rice [11,12], cotton [13], tomato [14], potato [15], muskmelons [16], and sugarcane [17][18][19]. However, the response mechanisms of non-heading Chinese cabbage to cold stress, which is not very clear, and further joint research is needed.
Nowadays, transcriptome analysis, which has gradually become a useful and general tool for discovering genes in multiple stress pathways, including determining the expression patterns of related genes [20][21][22]. The related genes reported in plant secondary metabolism, which can be discovered based on experimental methods of functional genome sequencing [23,24]. In addition, RNA sequencing technology, which 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 post-transcriptional force, including regulation of essential genes, during secondary metabolite biosynthetic pathways [25][26][27]. Nowadays, transcriptome sequencing, which has been successfully used to detect expression levels of related genes in many organisms, such as rice [28], yeast [29,30], sweet potato [31], and taxus [32].
With Illumina sequencing technology, read millions of sequences at a time, and map individual assembled genes into a reference transcriptome map for molecular annotation [33].
Studies have found that microRNA (miRNA), trans short interfering RNA (ta-siRNA), and heterochromatic short interfering RNA (hc-siRNA), all of which play important roles in different organisms [34]. Among them, miRNA is composed of about 22 nucleotides (nt), which is a type of endogenous small RNA, and usually also plays a negative role in regulating gene expression [35]. Many studies have shown that miRNAs, which are often involved in plant development, hormone signaling, and abiotic stress responses [36,37].
Generally speaking, small interfering RNAs are processed from perfect double-stranded RNA, while miRNAs are derived from single-stranded RNA transcripts, which form an imperfect double-stranded stem-loop precursor structure [38][39][40][41]. On the whole, miRNA, which 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 [35,37,[42][43][44][45][46]. 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 on global change. We identified the some 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 explored their interactions. Quantitative real-time polymerase chain reaction (qRT-PCR), which was also used to assess the expression levels of common differentially expressed genes (DEGs), and identify those candidate genes involved in cold tolarence. This work might serve as a reference that the functional analysis of cold tolarence in non-heading Chinese cabbage.

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][49][50].

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 GOslim 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 7 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.

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  (Table 1S). In addition, we performed cluster analysis on all screened differentially expressed genes ( Figure 3C; Figure 4C; Figure 5C; Figure 6C).

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. R 2 : 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 upregulated 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.
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 plantpathogen 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).

Discussion
As a convenient tool for transcriptome analysis, RNA-Seq, has been increasingly targeted at species that do not have access to genomic sequences [53][54][55]. Here, we used RNA-Seq to measure gene expression levels in non-heading Chinese cabbage 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 datas from each groups, we found that between the BcL.1-25 and BcL.1-4 groups, the most differentially expressed genes (6,208 DEGs), which were enriched, and the number was almost 6 times that of the other groups, and the difference is the largest. At the same time, the ribosome pathway, which is worth mentioning, the most abundant genes (287 genes), which are also shown, far more than 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 [56,57]. This also suggests that the ribosomal pathway, which might be involved in cold stress.
RNA-Seq searched for some low temperature-related differentially expressed genes. qRT-PCR, a practical method, often used to identify gene expression levels [58]. To confirm the reliability of the RNA-Seq results, we performed qRT-PCR experiment to identify them.
According to qRT-PCR results, the expression patterns of all unigenes, which were consistent with the transcriptome sequencing data, showed that our experimental results were reliable.
small RNAs are short, non-coding RNAs, which usually 19-25 nt in length, and two protruding sizes of 21 and 24 nt, respectively [59]. In general, miRNA corresponds to the 21 nt class of small molecule RNA, research have also found that small RNAs shown a wide range of functions, including heterochromatin formation, gene silencing, and DNA methylation [60,61]. By small RNA sequencing, the result shows that 24 nt length, which is the most abundant category, among pure and unique reads in all groups (Figure 8). This result was highly consistent with previous studies on A. thaliana [ 62], Oriza sativa [ 63], Medicago truncatula [ 64], and Populus trichocarpa [ 65].
Further, the known miRNAs are analyzed, and new miRNAs are predicted. We discover that these selected known miRNAs, which have at least two stem-loop structures that are obtained by self-folding ( Figure 9A-9J), suggesting that there might be protein or chromosomal binding sites. Meanwhile, these predicted miRNAs, which have at least two stem-loop structures that are obtained by self-folding ( Figure 10A-10J), suggesting that there might also be protein or chromosomal binding sites. The functions of these miRNAs and their target genes, which were comprehensively analyzed, 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 annotated in biological processes, and the KEGG pathway, which was also used to describe the corresponding metabolic pathways. After low temperature treatment, most target genes, which were enriched in the metabolic pathways of starch and sucrose, which was basically consistent with the results of the above-mentioned transcriptome study. This indicated that a series of biosynthetic and metabolic pathways, which might be induced, after low temperature treatment [66]. Here, a combination of transcriptome and small RNA sequencing, which 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.

Conclusions
In this study, a total of 63.43 Gb clean data, which was obtained from the transcriptome analysis. Based on the comparison results, a total of 1,860 new genes were discovered, 13 common DEGs, which were detected in all treatments, including 7 up-regulated and 6 down-regulated. Some of the DEGs, which were involved in response to stress, especially in terms of abiotic stress, such as freezing, cold and salt stresses. Meanwhile, we used qRT-PCR experiment to confirm changes in the expression levels of these genes. Further, we performed miRNA sequencing analysis on the same material. We found that the results revealed a total of 34,182,333 small RNA reads, 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 is 8. The number of known mature miRNAs and the number of precursors, which were 110 and 80, respectively; the number of predicted novel miRNA matures and the number of precursors, which were 75 and 84, respectively. According to GO and KEGG enrichment analysis, single-organism 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 plantpathogen interaction pathway. Additionally, our findings highlight 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, which used in this study are non-heading Chinese cabbage varieties, hours, respectively. The leaves of plants, which were harvested, immediately frozen in liquid nitrogen, and stored at -80°C for experimental used.

Transcriptome sequencing
Following the manufacturer's instructions to extract total RNA from the samples, by using RNAiso Plus reagent (Takara Bio, Dalian, China). RNA samples, which were examined, by using a spectrophotometer and electrophoresed on a 1% agarose gel. cDNA libraries construction and transcriptome sequencing, which were performed by the Biomarker Technologies (Beijing, China). The using Illumina HiSeq X Ten sequencing platform and Pe150 mode sequencing, all clean reads, which were subsequently mapped to the Brassica rapa reference genome sequence (IVFCAASv1) (http://brassicadb.org). The clean reads of each sample, which were sequence aligned with the designated reference genome. Gene expression analysis, which was performed, based on the comparison results; differentially expressed genes, which were identified, based on their expression levels in different samples, and their functional annotation and enrichment analysis were aslo performed [67].

Transcriptome assembly and functional annotation
The raw data of the transcriptome sequencing, which 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 Trinity program [68].
Firstly, the certain short reads with overlap regions, which were assembled into longer contiguous sequences for each library. Then, the distance of different contigs, which 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 TGI Clustering tool [69]. The new genes discovered, which were performed with NR [70], Swiss-Prot [71], GO [72], COG [73], KOG [74], Pfam [75], and KEGG [76] databases, using BLAST [77] software. Using KOBAS2.0 [78] to obtain the KEGG orthology result of the new gene for sequence alignment. After predicting the amino acid sequence of the new gene, using HMMER [79] software to align with the Pfam database, in order to obtain the annotation information of the new gene.

Differentially expressed genes analysis
Following fragments per kilobase of exon per million fragments mapped reads (FPKM) method, the expression level of unigene, which was calculated [80]. The ratio of the FPKM values (using 0.001 instead of 0 if the FPKM was 0), which were taken as the fold-changes in the expression of each gene to identify DEGs between each groups. The false discovery rate (FDR) control method, which was used to identify the threshold of the p-value in multiple tests, in order to compute the significance of the difference in transcript abundance [81]. In this result, only fold change with |log 2 (case_FPKM/control_FPKM)| ≥ 1, and an FDR ≤ 0.001, which 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 [82]. Meanwhile, the DEGs, which were annotated with GO and KEGG databases.

Validation of DEGs with qRT-PCR
qRT-PCR, which 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, which was synthesized, using a PrimeScript™ II First Strand cDNA synthesis kit (Takara Bio, Dalian, China), according to the manufacturer's protocol. The primers were designed, using the software Premier 5.0, and listed in Table 1. The quantified expression levels of the tested genes were normalized against the housekeeping genes Cyclophilin 1 (CYP1) 19 [83]. The qRT-PCR assays, which were performed with three biological, and technical replicates. Each reaction, which was performed in 20 μL reaction mixtures containing a diluted cDNA sample as template, SYBR Premix Ex Taq (2×) (TaKaRa, Kyoto, Japan) and gene-specific 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, which was performed according to previous report [84]. The comparative Ct value method, which was adopted to analyze the relative gene expression, according to a previous analysis, RNA expression levels relative to actin gene, which were calculated as 2 -ΔΔCT [ 84,85].

small RNA sequencing
High-throughput sequencing (such as Illumina HiSeqTM2500 / MiSeq and other sequencing platforms) sequenced raw image data files, which are converted into sequenced reads by base calling analysis, we call them raw data or raw reads. As a result, it is stored in FASTQ (abbreviated as fq) file format, which contains the sequence information of read sequences and its corresponding sequencing quality information. Referring to the standard definition of miRNA [86,87], the candidate target gene of miRNA, which was compared as the query sequence with the Brassica rapa database (http://brassicadb.org/brad/). The control samples and inoculation samples, which were mixed for small RNA libraries construction, respectively. According to the reported procedures, the construction of small RNA libraries, which were completed [88]. Three micrograms of total RNA per sample, which was used as input material for the small RNA library [89]. small RNA library construction and small RNA deep sequencing, which were proceeded following the detailed protocol provided by the genome sequencing company (Novogene, China).

Bioinformatic analysis of sequence data
Raw data, which were first processed through custom Perl and Python scripts to obtain clean data. The clean data, which were mapped to the reference sequence in miRBase21.0 20 by Bowtie [90], without mismatch to look for known miRNAs. Then, the other reads, which were integrated to predict novel miRNAs using the available miREvo [91] and miRDeep2 [92] software. The miRNA counts as well as base bias, which were identified by using custom scripts. Then, the miFam.dat (http://www.mirbase.org/ftp.shtml) was used to look for families of known miRNAs. The novel miRNA precursor was submitted to Rfam (http://rfam.sanger.ac.uk/search/) 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), which was used to estimated miRNA expression levels [93]. DEG seq R package, which was used to analyze the differential expression of two samples with the criterion of Q < 0.01 and |log 2 (fold-change)| > 1 [94].

Construction of degradome libraries
Target genes of candidate miRNAs, which were verified by degradome sequencing by using total RNA same to the RNA used for small RNA sequencing library construction, following the published parallel analysis of RNA Ends (PARE) protocol [95]. 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, which was performed to predict target genes of miRNA [96].
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 snnotation suites. Subsequently, Revigo tool (http://revigo.irb.hr), which was implemented for enrichment analysis of the target genes. The KOBAS software, which was used to test the statistical enrichment of the target gene candidates in the

Consent for publication
Not applicable.

Availability of data and materials
The datasets supporting the conclusions of this manuscript are included within the article and its additional files.

Competing interest
The authors declare no any competing interest.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.