Integrated RNA-Seq and sRNA-Seq Analysis Identifies Chilling and Freezing Responsive Key Molecular Players and Pathways in Tea Plant (Camellia sinensis)

Abstract

Tea [Camellia sinensis (L) O. Kuntze, Theaceae] is one of the most popular non-alcoholic beverages worldwide. Cold stress is one of the most severe abiotic stresses that limit tea plants' growth, survival and geographical distribution. However, the genetic regulatory network and signaling pathways involved in cold stress responses in tea plants remain unearthed. Using RNA-Seq, DGE and sRNA-Seq technologies, we performed an integrative analysis of miRNA and mRNA expression profiling and their regulatory network of tea plants under chilling (4℃) and freezing (-5℃) stress. Differentially expressed (DE) miRNA and mRNA profiles were obtained based on fold change analysis, miRNAs and target mRNAs were found to show both coherent and incoherent relationships in the regulatory network. Furthermore, we compared several key pathways (e.g., 'Photosynthesis'), GO terms (e.g., 'response to karrikin') and transcriptional factors (TFs, e.g., DREB1b/CBF1) which were identified as involved in the early chilling and/or freezing response of tea plants. Intriguingly, we found that karrikins, a new group of plant growth regulators, and β-primeverosidase (BPR), a key enzyme functionally relevant with the formation of tea aroma might play an important role in both early chilling and freezing response of tea plants. Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) analysis further confirmed the results from RNA-Seq and sRNA-Seq analysis. This is the first study to simultaneously profile the expression patterns of both miRNAs and mRNAs on a genome-wide scale to elucidate the molecular mechanisms of early responses of tea plants to cold stress. In addition to gaining a deeper insight into the cold resistant characteristics of tea plants, we provide a good case study to analyse mRNA/miRNA expression and profiling of non-model plant species using next-generation sequencing technology.

Fig 1. Overview of DE mRNAs, DE miRNAs and length distribution of the small RNA among libraries.

The histograms showing the number of DE mRNAs (A) and miRNAs (D) among libraries. The Venn diagrams showing the overlaps among four comparisons of DE mRNAs (B) and miRNAs (E). (C) refers to the length distribution of the small RNA in different libraries.

Fig 2. Hierarchical clustering of DE mRNAs and DE miRNAs among libraries.

Hierarchical clustering of the DE mRNAs (A) and DE miRNAs (B), using the DGE and sRNA-Seq data derived from five samples (CT4, CT8, FT4, FT8 and CK samples) based on log₁₀ (RPKM+1) and log₁₀ (TPM+1) values, respectively.

Fig 3. Distribution of differentially expressed transcription factors.

The histograms showing the number of up- or down-regulated transcription factors in CT4 (A), CT8 (B), FT4 (C) and FT8 (D) samples, respectively.

Fig 4. Biological process network of GO term enrichment for differentially expressed (DE) mRNAs.

Over-represented GO terms for DE mRNAs in CT4 (A), CT8 (B), FT4 (C) and FT8 (D) samples, respectively. The node size represents the number of genes associated to a given GO term and node filled color reflects the adjusted P-value. End nodes were indicated by blue label and green border.

Fig 5. KEGG pathway analysis of differentially expressed (DE) mRNAs.

A total of 10, 17, 10 and 6 significantly enriched KEGG pathways were identified in CT4, CT8, FT4 and FT8 samples, respectively. Line color of edge represents the P-value of pathway.

Fig 6. Numbers of miRNA member in each family in C. sinensis.

Fig 7. miRNA-mRNA correlation network.

Circuits in red, blue, yellow and green ellipse are indicating CT4, CT8, FT4 and FT8 samples, respectively. The triangles and circles in the network are indicating miRNAs and target mRNAs, respectively. Down-regulated mRNAs and miRNAs are shown as green and up-regulated mRNAs and miRNAs are shown as red.

Fig 8. qRT-PCR validation for DE miRNAs and mRNAs.

Twelve DE mRNA and six DE miRNA were selected for the quantitative RT-PCR analysis, including AHK4 (A), DREB1B (B), DREB1C (C), NAC100 (D), ARF6 (E), BPR (F), GI (G), GST3 (H), PP2A13 (I), HEXO1 (J), HT1 (K), ATH1 (L), csn-miR156.17(M), csn-miR398.5(N), csn-miR164.4(O), csn-miR167.6(P), csn-miR164.1(Q) and csn-miR171.2(R). The GADPH gene and 5.8S rRNA were chosen as the endogenous control for mRNA and miRNA qRT-PCR, respectively.

Fig 9. Differential expression of carbohydrate metabolism and plant hormone signal transduction related genes.

The histograms showing the number of up- or down-regulated carbohydrate metabolism (A) and plant hormone signal transduction (B) related genes in CT4, CT8, FT4 and FT8 samples, respectively.