Comparative Analysis of miRNA Expression Profiles between Heat-Tolerant and Heat-Sensitive Genotypes of Flowering Chinese Cabbage Under Heat Stress Using High-Throughput Sequencing

Heat stress disturbs cellular homeostasis, thus usually impairs yield of flowering Chinese cabbage (Brassica campestris L. ssp. chinensis var. utilis Tsen et Lee). MicroRNAs (miRNAs) play a significant role in plant responses to different stresses by modulating gene expression at the post-transcriptional level. However, the roles that miRNAs and their target genes may play in heat tolerance of flowering Chinese cabbage remain poorly characterized. The current study sequenced six small RNA libraries generated from leaf tissues of flowering Chinese cabbage collected at 0, 6, and 12 h after 38 °C heat treatment, and identified 49 putative novel miRNAs and 43 known miRNAs that differentially expressed between heat-tolerant and heat-sensitive flowering Chinese cabbage. Among them, 14 novel and nine known miRNAs differentially expressed only in the heat-tolerant genotype under heat-stress, therefore, their target genes including disease resistance protein TAO1-like, RPS6, reticuline oxidase-like protein, etc. might play important roles in enhancing heat-tolerance. Gene Ontology (GO) analysis revealed that targets of these differentially expressed miRNAs may play key roles in responses to temperature stimulus, cell part, cellular process, cell, membrane, biological regulation, binding, and catalytic activities. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis identified their important functions in signal transduction, environmental adaptation, global and overview maps, as well as in stress adaptation and in MAPK signaling pathways such as cell death. These findings provide insight into the functions of the miRNAs in heat stress tolerance of flowering Chinese cabbage.


Introduction
Flowering Chinese cabbage (Brassica campestris L. ssp. chinensis var. utilis Tsen et Lee) belongs to the Brassicaceae family and is mainly grown and consumed in southern China [1]. This vegetable crop is valuable for human diet due to its high soluble fiber, favorable taste, richness in vitamin catalytic and binding activities. The comparative study identified differentially expressed miRNAs and their involvements in controlling heat tolerance in flowering Chinese cabbage.

Plant Materials, Growth Conditions, and Total RNA Isolation
Two genotypes of flowering Chinese cabbage (Brassica campestris L. ssp. chinensis var. utilis Tsen et Lee),  and Liuye 50 (HS), were grown in a growth chamber at Guangzhou University at 28/22 • C for 14/10 h (day/night). Plants at the five-leaf stage were transferred into another growth chamber at 38 /29 • C (14/10 h) for heat treatments. Samples were collected from the fully expanded upper leaves of HT and HS plants after 0 (control), 6 and 12 h of heat treatments. Collected tissues were flash-frozen immediately in liquid nitrogen, and then stored at 80 • C until RNA isolation [22,23]. Trizol RNA extraction kit (Invitrogen, Waltham, MA, USA) was used for total RNA isolation following the manufacturer's protocol.

Construction and Sequencing sRNA Libraries
RNA was isolated from three biological replicates at each time point and then three RNA samples per treatment were combined into one tube in an equal amount for library construction. All the cDNA libraries were constructed using TruSeq Small RNA Preparation Kit following the manufacturer's protocol (Illumina, San Diego, CA, USA). In brief, RNA 3 -and RNA 5 -adapters were ligated to total RNA, cDNA constructs were created using reverse transcription after PCR, and then samples of small cDNA fragments of different lengths (18-30 nt) were run on 6% denaturing polyacrylamide gel by electrophoresis [18,24,25]. The final cDNA libraries were sequenced using Illumina HiSeq at the Beijing Genomics Institute (BGI, Shenzhen, China).

Identification of Conserved and Novel miRNAs
After removing poly-A tags, no-insert tags, adapter sequences, oversized insertion tags, 5'-primer contaminants, and small tags (sequences beyond 15-30 nucleotides or without 3' primers), remaining sequence reads were further analyzed using the Bowtie2 web program to determine the length distribution of the sRNAs by mapping the clean reads to other sRNA databases and to the reference genome [26].
The unique sRNAs were aligned to known ncRNAs in the Rfam database (http://www.sanger.ac. uk/science/tools/rfam) to remove snoRNA, rRNA, snRNA, tRNA, and scRNA using NCBI BLASTN. Perfectly matched reads were excluded from further analysis; and remaining sequences were compared to Brassica database (http://brassicadb.org/brad/) to determine mismatched and matched sequences. These reads with no more than three mismatched nucleotides were considered as conserved miRNA candidates, whereas those with more than three unmatched sequences were considered as putative novel miRNA. miRBase software [17] was used for further prediction of the novel miRNAs.

Analysis of Differentially Expressed miRNAs
The levels of miRNA expression were compared between HT and HS genotypes after the heat treatments to determine differentially expressed miRNAs in flowering Chinese cabbage. To predict heat tolerance associated miRNAs in flowering Chinese cabbage, the fold-change of miRNA was determined as the ratio of miRNA expressions between HT and HS lines. The false discovery rate (FDR) was adjusted by analyzing significant p-value thresholds in different tests [18]. The normalized miRNA expression level was used to determine the fold changes (log 2 ratio) of miRNA expression in each sample. To avoid calculation error, miRNA expression level was normalized and converted to transcripts per million (TPM) from 0 to 0.01 in all libraries. The minimum criteria for comparative analysis of low expression of miRNAs was adjusted as if miRNA had normalized expression of < 1 in all libraries. The normalization equation is as follow: Normalized expression = actual miRNA read counts/total counts of clean reads × 10 6 . The fold-change values and p-values were calculated using normalized data, and the fold-change values were used to generate a scatter plot: Fold-change = log 2 (N 2 /N 1 ). The p-value was determined as follows: where x and y denote total clean sRNA reads, while N 1 and N 2 are the normalized miRNA expression levels in the control and the treatment, respectively.

Prediction of miRNA Secondary Structure
The Zuker folding algorithm implemented in Mfold (http://mfold.rna.albany.edu/?q=mfold) was used to determine the secondary structures of miRNAs using default parameters [27]. Minimal folding free energy (MFE) and minimal free energy index (MFEI) parameters were used to differentiate the miRNAs from other sRNA sequences. Moreover, sRNA sequences sustaining Meyers guidelines were assumed as prospective miRNAs with following considerations: mature miRNAs have no more than one bulge and the bulge size is not higher than two, mismatch sequences should be less than three, and high MFEI values and high negative MFE values must be in predicted secondary structures and properly fold into stem-loop hairpin structures [28].

Target Prediction of Differentially Expressed miRNAs
To identify and analyze differentially expressed miRNAs, the software packages TargetFinder [29] and TAPIR [30] were used as described earlier [31]. To achieve reliable results with a confidence interval, only common binding sites that were predicted by both tools were chosen for further analysis.

GO and KEGG Prediction of miRNA-Related Regulatory Pathways
Gene Ontology (GO; http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (http://www.genome.jp/kegg/pathway.html) databases were analyzed to predict the key regulatory pathways using a corrected p-value (≤ 0.05) with a threshold derived from a hypergeometric test [32]. GO and KEGG analyses were classified sequence reads into different groups based on their function, including biological processes, cellular components, and molecular functions.

Validation of miRNAs Using RT-qPCR
To validate the sequencing data, four novel miRNAs and four conserved miRNAs-whose expressions were either down-or up-regulated-were selected randomly for RT-qPCR as previously described [33]. RT-qPCR used the same RNA samples used for sequencing that were collected from 0, 6, and 12 h of heat-treated heat-tolerant and heat-sensitive genotypes. Reverse transcription was carried out using RNA-tailing and primer-extension reverse transcription (RT)-PCR. For each quantitative real-time (qRT)-PCR, 2 µl template cDNA was mixed with 10 µl 2× SYBR Green PCR master mix (Takara Bio Inc., Kusatsu, Japan) and 5 pmol each of the forward and reverse primers in a final volume of 20 µl. The amplification program started at 95 • C for 5 min, then for 40 cycles of 95 • C for 15 s, 60 • C for 20 s, and 70 • C for 20 s, followed by a thermal denaturation step to generate the dissociation curves for verification of amplification specificity. Specific primers used in this study were listed in Table S1 with U6 as the internal control. Relative miRNA expression levels were quantified using the previously described 2 -∆∆CT method [34].

Statistical Analysis
Statistical analysis was conducted using the statistical product and service solution (SPSS) software version 22.0 (IBM Corp., Chicago, USA). The significant differences between treatments were determined using Student's t-test or Tukey's test for multiple comparisons after one-way analysis of variance at significance level of p < 0.05 or p < 0.01. All the results were presented as means ± SEM.

Analysis and Classification of sRNAs Sequence
Sequencing of six sRNA libraries constructed using samples from 0, 6, and 12 h of heat stress yielded 28 To examine which type of sRNAs is associated with heat tolerance of flowering Chinese cabbage, Rfam database was searched to categorize these clean sequence reads into different classes (Table 1). rRNA reads were 0.36%-1.26% in HT genotypes and 0.69% in HS genotypes. tRNA read was higher (0.26%-5.38%) in HT genotypes than in HS genotypes (0.09%-0.16%). In addition, snRNA was 0.06%-0.12% for HT genotypes and 0.03%-0.07% for HS genotypes; snoRNA was 0.03% for both HT and HS genotypes. After removal of tRNA, rRNA, snRNA, and snoRNA, 22,222,402, 13,150,565, and 17,741,890 unique reads from the HT genotype and 25,021,953, 13,857,771, and 24,056,406 unique reads from the HS genotype treated with 0, 6, and 12 h under high temperature were mapped, respectively.
The sequence read length distribution patterns of sRNAs were similar between the HT and HS libraries in general, but a higher proportion of sRNAs were observed in the HS genotype (79.1%) than in the HT genotype (50.7%). The lengths of sRNAs were mainly from 21 to 24 nt ( Figure 1) with the predominant length of 21 nt, followed by 24, 23, and 22 nt in both the HT and HS genotypes (Table S3).

Known miRNAs from Flowering Chinese Cabbage
To identify the known miRNAs from flowering Chinese cabbage, sRNA sequences were searched against the known miRNAs in miRBase (release 17.0) to find these reads with no more than three mismatched nucleotides as known miRNAs. The lengths of the known miRNAs usually ranged from 21-24 nt. Sixty-two small RNAs were identified to have identical sequences to B. campestris in miRBase and were considered as known bra-miRNAs (Table S4). Although many miRNAs had a relatively low expression, they showed significantly differential expression between HT and HS genotypes. The top ten abundantly expressed miRNAs were bra-miR398-3p, bra-miR168a-5p, bra-miR396-5p, bra-miR168b-5p, bra-miR171e, bra-miR160a-5p, bra-miR159a, bra-miR162-3p, bra-miR171a, and bra-miR156a-5p.
genotype treated with 0, 6, and 12 h under high temperature were mapped, respectively.
The sequence read length distribution patterns of sRNAs were similar between the HT and HS libraries in general, but a higher proportion of sRNAs were observed in the HS genotype (79.1%) than in the HT genotype (50.7%). The lengths of sRNAs were mainly from 21 to 24 nt ( Figure 1) with the predominant length of 21 nt, followed by 24, 23, and 22 nt in both the HT and HS genotypes (Table  S3).

Known miRNAs from Flowering Chinese Cabbage
To identify the known miRNAs from flowering Chinese cabbage, sRNA sequences were searched against the known miRNAs in miRBase (release 17.0) to find these reads with no more than three mismatched nucleotides as known miRNAs. The lengths of the known miRNAs usually ranged from 21-24 nt. Sixty-two small RNAs were identified to have identical sequences to B. campestris in miRBase and were considered as known bra-miRNAs (Table S4). Although many miRNAs had a relatively low expression, they showed significantly differential expression between HT and HS genotypes. The top ten abundantly expressed miRNAs were bra-miR398-3p, bra-miR168a-5p, bra-miR396-5p, bra-miR168b-5p, bra-miR171e, bra-miR160a-5p, bra-miR159a, bra-miR162-3p, bra-miR171a, and bra-miR156a-5p.

Identification of Novel miRNAs
Novel miRNA candidates were identified on the basis of MFE value, the miRNA/miRNA* duplex and secondary structure of precursor sequences. The stem-loop hairpin secondary structures were predicted from precursor sequences, suggesting that most of the identified miRNAs rely on the 5′ arm of the hairpin structure. To determine if novel miRNAs were involved in heat tolerance in flowering Chinese cabbage, all mappable sRNA sequences were searched against the Brassica database and the miRBase to eliminate previously known miRNAs. Any sRNAs that could be exactly mapped to the reference genome but not as conserved miRNAs were assumed to be novel miRNA

Identification of Novel miRNAs
Novel miRNA candidates were identified on the basis of MFE value, the miRNA/miRNA* duplex and secondary structure of precursor sequences. The stem-loop hairpin secondary structures were predicted from precursor sequences, suggesting that most of the identified miRNAs rely on the 5 arm of the hairpin structure. To determine if novel miRNAs were involved in heat tolerance in flowering Chinese cabbage, all mappable sRNA sequences were searched against the Brassica database and the miRBase to eliminate previously known miRNAs. Any sRNAs that could be exactly mapped to the reference genome but not as conserved miRNAs were assumed to be novel miRNA candidates. To increase the accuracy of novel miRNA prediction, the miRNA/miRNA* criterion was evaluated. A total of 49 novel miRNA candidates were identified from the HT and HS libraries of flowering Chinese cabbage (Table S5). The mean MFE value predicted for pre-miRNAs was -42.13 kcal/mol, ranging from -23.6 to -217.2 kcal/mol. The identified novel miRNA length varied from 21 to 24 nt. Novel-mir09, novel-mir112, novel-mir125, novel-mir149, novel-mir187, novel-mir202, and novel-mir248 had striking secondary structures with lower MFE values and were considered as key putative miRNAs ( Figure 2). Of all the novel miRNAs, novel-mir202, novel-mir225, novel-mir255, novel-mir248, novel-mir187, novel-mir170, and novel-mir99 had the highest expression levels. candidates. To increase the accuracy of novel miRNA prediction, the miRNA/miRNA* criterion was evaluated. A total of 49 novel miRNA candidates were identified from the HT and HS libraries of flowering Chinese cabbage (Table S5). The mean MFE value predicted for pre-miRNAs was -42.13 kcal/mol, ranging from -23.6 to -217.2 kcal/mol. The identified novel miRNA length varied from 21 to 24 nt. Novel-mir09, novel-mir112, novel-mir125, novel-mir149, novel-mir187, novel-mir202, and novel-mir248 had striking secondary structures with lower MFE values and were considered as key putative miRNAs ( Figure 2). Of all the novel miRNAs, novel-mir202, novel-mir225, novel-mir255, novel-mir248, novel-mir187, novel-mir170, and novel-mir99 had the highest expression levels.

Differential expression profiling of known and novel miRNAs
To investigate the heat induced miRNAs in flowering Chinese cabbage, differentially expressed miRNAs were identified between heat treated (6 and 12 h heat treatments) and non-heat treated (0 h control) HT and HS genotypes (Figure 3). A total of 43 known and 49 novel miRNAs were differentially expressed in at least one of two time points (6 and 12 h) of heat treatments.

Differential Expression Profiling of Known and Novel miRNAs
To investigate the heat induced miRNAs in flowering Chinese cabbage, differentially expressed miRNAs were identified between heat treated (6 and 12 h heat treatments) and non-heat treated (0 h control) HT and HS genotypes (Figure 3). A total of 43 known and 49 novel miRNAs were differentially expressed in at least one of two time points (6 and 12 h) of heat treatments.
Among the heat induced known miRNAs, 21 (4 up-regulated and 17 down-regulated) were differentially expressed in both HT and HS genotypes, whereas nine (1 up-regulated and 8 down-regulated) were differentially expressed only in the HT genotype (Table 2) and 13 (2 up-regulated and 11 down-regulated) were differentially expressed only in the HS genotype (Table S6). Among the novel miRNAs, 6 (1 up-regulated and 5 down-regulated) were differentially expressed in both HT and HS genotypes, 14 (3 up-regulated and 11 down-regulated) were differentially expressed in HT genotype only, and 29 (7 up-regulated and 22 down-regulated) were differentially expressed in HS genotype only (Table S7). The five differentially expressed known miRNAs and novel miRNAs that were common in both HT and HS genotypes are listed in Table 3. Among the heat induced known miRNAs, 21 (4 up-regulated and 17 down-regulated) were differentially expressed in both HT and HS genotypes, whereas nine (1 up-regulated and 8 downregulated) were differentially expressed only in the HT genotype (Table 2) and 13 (2 up-regulated and 11 down-regulated) were differentially expressed only in the HS genotype (Table S6). Among the novel miRNAs, 6 (1 up-regulated and 5 down-regulated) were differentially expressed in both HT and HS genotypes, 14 (3 up-regulated and 11 down-regulated) were differentially expressed in HT genotype only, and 29 (7 up-regulated and 22 down-regulated) were differentially expressed in HS genotype only (Table S7). The five differentially expressed known miRNAs and novel miRNAs that were common in both HT and HS genotypes are listed in Table 3.

Functional Annotation of miRNA Target Genes
GO analysis on the putative miRNA target genes further divided them into three function categories: biological process, molecular function, and cellular component (Figure 4; Table S9). In the biological process category, the most enriched GO terms include biological regulation, cellular process, metabolic process, response to temperature stimulus, regulation of biological process, and single organism process. Interestingly, these target genes might play a significant role in diversifying biological processes such as signaling, localization, response to stimulus, and developmental process. In cellular component category, the most enriched GO terms include cell, cell part, membrane, macromolecular complex, organelle, and organelle part. The majority of molecular functions include nucleic acid binding transcription factor activity, catalytic activity, and transporter activity and binding.  KEGG analysis predicted the potential pathway networks that were involved in response to heat stress in flowering Chinese cabbage including transport and catabolism, signal transduction, folding, sorting and degradation, environmental adaptation, replication and repair, transcription, and translation ( Figure 5; Table S10). KEGG analysis predicted the potential pathway networks that were involved in response to heat stress in flowering Chinese cabbage including transport and catabolism, signal transduction, folding, sorting and degradation, environmental adaptation, replication and repair, transcription, and translation ( Figure 5; Table S10).
The KEGG pathway analysis also predicted various miRNA-regulated pathways based on the miRNA targeted genes identified and demonstrated enriched target genes including genes for stress response, stress tolerance, and stress adaptation [32,35]. MEKK1 and SnRK2 were up-regulated and their expression might play a key role in responses to cold/salt stresses, pathogen infection, and drought tolerance in HT genotype. SUMM2 and WRKY33 were expressed in both HT and HS genotypes and might be involved in regulation of cell death defense response and camalexin synthesis, respectively ( Figure S1). The KEGG pathway analysis also predicted various miRNA-regulated pathways based on the miRNA targeted genes identified and demonstrated enriched target genes including genes for stress response, stress tolerance, and stress adaptation [32,35]. MEKK1 and SnRK2 were up-regulated and their expression might play a key role in responses to cold/salt stresses, pathogen infection, and drought tolerance in HT genotype. SUMM2 and WRKY33 were expressed in both HT and HS genotypes and might be involved in regulation of cell death defense response and camalexin synthesis, respectively ( Figure S1).

Validation of Differentially Expressed miRNAs by RT-qPCR
To confirm the results from the RNA sequencing, four novel and four conserved miRNAs were randomly selected to represent both down-and up-regulated miRNAs for RT-qPCR. A strong correlation (r 2 = 0.843) was detected between RT-qPCR and RNA-seq data (Figure 6), indicating a good agreement in expression levels (down-regulation or up-regulation) between RT-qPCR and RNA sequencing data for the selected miRNAs, confirming that the differentially expressed miRNAs of B. campestris ssp. Chinensis predicted by RNA sequencing were real.

Validation of Differentially Expressed miRNAs by RT-qPCR
To confirm the results from the RNA sequencing, four novel and four conserved miRNAs were randomly selected to represent both down-and up-regulated miRNAs for RT-qPCR. A strong correlation (r 2 = 0.843) was detected between RT-qPCR and RNA-seq data (Figure 6), indicating a good agreement in expression levels (down-regulation or up-regulation) between RT-qPCR and RNA sequencing data for the selected miRNAs, confirming that the differentially expressed miRNAs of B. campestris ssp. Chinensis predicted by RNA sequencing were real.

miRNA in Flowering Chinese Cabbage
Vegetable plants have to cope with diverse environmental stresses that may trigger numerous gene regulatory mechanisms, such as orchestration of gene expression regulation at a post-transcriptional level, and re-establishing and restoring cellular homeostasis, reducing cell-cycle regulation, and adaptive growth response [36,37]. Since the plant miRNA may involve in many functions in response to stress environments, it is an integral topic in functional genomic research. To date, numerous reports have revealed miRNA as key regulatory molecules that participate in plant metabolisms, stress responses, tissue development and many other functions [38,39]. The functional involvement of plant miRNAs in response to abiotic stresses was originally suggested by surveys of NCBI expressed sequence tags and profiling miRNA expression after challenging plants with certain stress stimuli to predict miRNA targets [40]. Flowering Chinese cabbage has much fewer miRNAs registered in miRBase compared to other vegetable crops and model plant Arabidopsis [11,41].
Using high-throughput sequencing, we identified 43 previously reported and 49 putative novel miRNAs in the HT and HS genotypes of flowering Chinese cabbage after heat treatments. sRNAs showed wide variation in sequence length, with 21 to 24 nt as most abundant, which agree with several previous reports [11,18]. Previously, we identified 41 conserved and 18 novel miRNAs in a flowering Chinese cabbage, Youlv 501, after heat treatment using high-throughput sequencing [21]. Likewise, using computational and deep sequencing methods, 24 novel and 161 known miRNAs from 51 families were identified in HT and HS genotypes of Brassica oleracea L. var italic [11]. High-throughput sequencing has been successfully used to identify novel and known miRNAs in response to stress in various species including soybean [42], maize [43], Populus euphratica [44], and rice [45]. In B. rapa, 21 novel miRNAs belonging to 19 miRNA families that were identified using NGS, bra-miR1885b.3 and bra-miR5718 were reported to be involved in response to heat [18]. Likewise, 221 conserved and 125 novel miRNAs were reported to play key functions in plant development and growth, metabolism, and stress responses in B. rapa [46]. Using a comparative genomics approach, 126 novel miRNAs were identified in B. juncea and were predicted to participate in regulation of different biological processes in response to drought, salinity, and high temperature stresses [20]. These findings reveal that more novel miRNAs with important regulatory functions in flowering Chinese cabbage can be identified using improved reference genome sequences.

miRNA-Targeted Gene-Networks Involved in Response of Flowering Chinese Cabbage to Heat Stress
Upon heat stress, several genes that are responsive to abiotic stress in plants can be directly targeted by miRNAs. Target mRNA cleavage depends upon its complementary miRNA. It has been shown that conserved miRNA miR398a targeted BracCSD1 gene, and there is an inverse relationship between miRNA and its target in Brassica rapa [18], indicating that the gene is heat-sensitive and the miRNA exhibits a heat-inhibitive function in response to heat stress. Furthermore, it has been revealed that several plants used more phosphate and nitrogen under heat stress than normal growth conditions, and this phenomena may be due to miRNA repression and roles of the target genes in adaptation and regulation of Pi starvation [58]. Moreover, miR156h and miR156g targeted SPL -family genes that were involved in floral transition and vegetative phase of B. rapa [19], and miR156a regulated SPL2, SPL3, SPL9, and SPL10 targets, but SPL2 showed a significant down-regulation under heat stress.
In conclusion, in this extensive study, we analyzed heat-responsive miRNAs from HT and HS genotypes of flowering Chinese cabbage after heat stress treatments and identified 49 novel and 43 known miRNAs that possess important regulatory roles in regulating heat stress responses by targeting numerous genes. The relationships between the target genes and the miRNAs in response to heat stress revealed from this study could help us design new breeding tools using biotechnological approaches for genetic improvement of heat tolerance of flowering Chinese cabbage.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/11/3/264/s1, Figure S1: KEGG pathway map for MAPK signaling for heat-tolerant (HT) and heat-sensitive (HS) genotypes of flowering Chinese cabbage; Table S1: List of primers used in this study; Table S2: Analysis of small RNA sequences from the libraries of HT and HS flowering Chinese cabbage; Table S3: The length distribution and percentage of detected reads from heat-responsive total small RNAs of the flowering Chinese cabbage HT and HS genotype by sequencing; Table S4: Expression profiles of heat-responsive conserved miRNAs in HT and HS genotypes of flowering Chinese cabbage at 6 h and 12 h; Table S5: Novel miRNAs identified from HT and HS flowering Chinese cabbage genotypes under heat stress; Table S6: Differentially expressed heat-responsive known miRNAs in flowering Chinese cabbage HT and HS genotypes at 6 h and 12 h; Table S7. Differentially expressed heat-responsive novel miRNAs in flowering Chinese cabbage HT and HS genotypes at 6 h and 12 h; Table S8: Potential targets of differentially expressed miRNAs between HT and HS Flowering Chinese cabbage under heat stress; Table S9; Gene Ontology (GO) annotation of target genes of heat stress-responsive microRNA in the HT and HS libraries of flowering Chinese cabbage; Table S10; KEGG classification of target genes of heat stress-responsive microRNAs in the HT and HS libraries of flowering Chinese cabbage.