Xing Huang1,2†
Yongsheng Liang3†Baoqing Zhang2Xiupeng Song2
Yangrui Li1,2Zhengqiang Qin2Dewei Li2
Rongfa Chen2Zhongfeng Zhou2Yuchi Deng2Jiguang Wei1*Jianming Wu2*- 1College of Agriculture, Guangxi University, Nanning, China
- 2Key Laboratory of Sugarcane Biotechnology and Genetic Improvement (Guangxi), Guangxi Key Laboratory of Sugarcane Genetic Improvement, Sugarcane Research Institute, Guangxi Academy of Agricultural Sciences, Ministry of Agriculture, Nanning, China
- 3Nanning Institute of Agricultural Sciences, Nanning, China
Cold stress causes major losses to sugarcane production, yet the precise molecular mechanisms that cause losses due to cold stress are not well-understood. To survey miRNAs and genes involved in cold tolerance, RNA-seq, miRNA-seq, and integration analyses were performed on Saccharum spontaneum. Results showed that a total of 118,015 genes and 6,034 of these differentially expressed genes (DEGs) were screened. Protein–protein interaction (PPI) analyses revealed that ABA signaling via protein phosphatase 2Cs was the most important signal transduction pathway and late embryogenesis abundant protein was the hub protein associated with adaptation to cold stress. Furthermore, a total of 856 miRNAs were identified in this study and 109 of them were differentially expressed in sugarcane responding to cold stress. Most importantly, the miRNA–gene regulatory networks suggested the complex post-transcriptional regulation in sugarcane under cold stress, including 10 miRNAs−42 genes, 16 miRNAs−70 genes, and three miRNAs−18 genes in CT vs. LT0.5, CT vs. LT1, and CT0.5 vs. LT1, respectively. Specifically, key regulators from 16 genes encoding laccase were targeted by novel-Chr4C_47059 and Novel-Chr4A_40498, while five LRR-RLK genes were targeted by Novel-Chr6B_65233 and Novel-Chr5D_60023, 19 PPR repeat proteins by Novel-Chr5C_57213 and Novel-Chr5D_58065. Our findings suggested that these miRNAs and cell wall-related genes played vital regulatory roles in the responses of sugarcane to cold stress. Overall, the results of this study provide insights into the transcriptional and post-transcriptional regulatory network underlying the responses of sugarcane to cold stress.
Introduction
Sugarcane planting, yield, and quality are affected by abiotic factors including low-temperature stress. Previous research has shown that the response of sugarcane to cold stress is regulated by sophisticated mechanisms at the morphological, anatomical, physiological, and biochemical levels. Once subjected to cold stress, sugarcane leaves dry up, growth points and lateral buds become frostbitten, and the tilling rate and number of effective perennial stems are affected, especially at the seedling stage. Cold stress markedly inhibits sugarcane seedling growth and survival rates fall (Rasheed et al., 2010; Pang et al., 2015). The thickness of sugarcane leaves, the area of vesicular cells, the thickness of thick-walled vascular bundle tissue, and electrolyte extravasation rates are all known to be positively correlated with cold resistance (Zhu et al., 2019). Indeed, at the seedling stage, the degree of peroxidation of the leaf lipid membrane under cold stress is enhanced and occurs in concert with membrane permeability and relative conductivity. The content of soluble sugar, proline, and 3,4-methylenedioxyamphetamine (MDA) also increases (Ding et al., 2006; Wang et al., 2016a), while chlorophyll content decreases with the decrease in temperature and increased duration of cold temperature treatment. Cold stress also inhibits the activity of ribulose-1,5-diphosphate carboxylase (Rubisco) and phosphoenolpyruvate carboxylase (PEP), while SOD (superoxide dismutase) and POD (peroxidase) increase. Traditional genetic and molecular analyses have identified that COLD1-CRLK1/2-MEKK-MPK-ICE-CBF-CORs and CRPK1-14-3-3-CBF are important cold signal sensing pathways in rice and Arabidopsis (Guo et al., 2018; Shi et al., 2018). Recent RNA-seq research has revealed the presence of regulated mechanisms that are more complex in other plants (Li et al., 2016). Almost 50% of paper mulberry transcriptional factors are involved in cold stress responses (Peng et al., 2015). In addition, 1,841 differentially expressed genes (DEGs) are involved in cold stress response in maize (Li et al., 2017). A total of 47 genes are upregulated in wheat, whereas just six are downregulated in the ICE1-CBF-COR pathway as a component of cold stress response. Ten genes in the ABA-dependent pathway are upregulated, while four are downregulated across these four genotypes, confirming that ABA accumulates when plants are exposed to low temperatures (Li et al., 2018).
The regulation of gene expression via microRNA (miRNA) is a negative mechanism that occurs in response to cold stress in plants at post-transcriptional levels. In plants, miRNAs are divided into two groups, conserved, and species-specific (Esposito et al., 2020). In recent years, numerous cold-responsive miRNAs have been identified in different plant species, and some of them are involved in the CBF-dependent pathway, triggering the elimination of ROS, modifying cellular antioxidant capacity, and in the auxin signaling pathway (Megha et al., 2018). miR319 is downregulated by cold stress in rice while the overexpression of miR319 or target RNA interference leads to the upregulation of cold-responsive genes and improved cold tolerance (Thiebaut et al., 2012; Yang et al., 2013; Wang et al., 2014). miR160 and miR167 have been identified as temperature-responsive miRNAs in cotton. The miR160 family also exhibits lower expression levels at 12°C compared with the control at 25°C (Wang et al., 2016b). One hundred and seventy-two miRNAs belonging to 39 families and 126 unique examples were identified in Nelumbo nucifera (Zou et al., 2020). 172 miRNA–target pairs have been identified as being involved in plant temperature responses by integrating miRNA with mRNA expression profiles. Three conserved and 25 predicted miRNAs have also been found to be involved in response to cold stress in Brachypodium (Zhang et al., 2009).
Prior to the sequencing of the complete sugarcane genome, expressed sequence tags (ESTs) and RNA-seq were both used to determine the mechanisms involved in cold stress response. In one example, 34 cold-inducible ESTs were identified in 2003 from 1,536 ESTs in sugarcane (Saccharum sp. cv SP80-3280) exposed to cold stress (Nogueira et al., 2003). Transcriptome profiles from the low temperature-tolerant S. spontaneum clone IND 00-1037 are also analyzed, as well as 2,583 genes that are upregulated and 3,302 that are downregulated in a low temperature stress treatment (Dharshini et al., 2016). A total of 170 cold-responsive TFs in 30 families are also differentially regulated (Selvarajan et al., 2018). Transcriptome profiles of S. spontaneum under cold stress indicate that certain genes involved in transmembrane transporter activities are absent from cold-susceptible sugarcane (Jong-Won et al., 2015). Additional work on small RNA transcriptomes from sugarcane has also been possible via deep sequencing (Sternes and Moyle, 2015). However, analysis of miRNA–genes under cold stress has not been undertaken since publication of the complete sugarcane genome (Zhang et al., 2018).
The wild relative of sugarcane, Saccharum spontaneum L., is known for its adaptability to environmental stressors, particularly cold stress conditions. Thus, to gain a comprehensive understanding of the molecular mechanisms involved in the response of sugarcane seedlings to cold stress, we studied the expression profiles and relevant genes in S. spontaneum GX87-16 seedlings treated with cold stress (4°C). In addition, we integrated a DEG and differentially expressed miRNA analysis and constructed a regulation network of miRNA–mRNA. The findings from this study provide insight into the transcriptional and post-transcriptional regulation of cold tolerance and low temperature stress in sugarcane.
Methods
Plant Materials and Cold Stress Treatment
Stems of S. spontaneum GX87-16 were planted in pots and grown in an artificial incubator at 28°C under a 16-h/8-h (light/ dark) photoperiod. After 2 weeks of growth, 2 h after entering the light culture, seedlings from six pots were transferred to an artificial incubator, which had been set to 4°C for cold stress treatment. Mature leaves were collected from treatment plants at 0.5 h and 1 h after cold treatment and denoted as LT0.5 and LT1. Meanwhile, another three potted plants were sampled at the time when the six pots were transferred to a 4°C artificial incubator; these three repeats were considered as controls (CT), namely, the 0 h of cold stress. Collected leaves were frozen in liquid nitrogen and stored at −80°C prior to RNA extraction.
Construction and Sequencing of mRNA-Seq and Small RNA Libraries
Total RNA was extracted using the Trizol reagent (Invitrogen), and RNA integrity was assessed using an RNA Nano 6000 Assay Kit in an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). RNA concentrations were measured using a Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, CA, United States). Equal amounts of RNA from all samples were used for the construction of mRNA-Seq and small RNA libraries.
A series of mRNA-Seq libraries were prepared using an Illumina TruSeq RNA Sample PrepKit following the manufacturer's protocols. Sequencing was performed via paired-end reads (2 × 151 base pairs, bp) using an Illumina X Ten sequencing platform. Bands of RNA between 18 and 40 nt in length were isolated for small RNA sequencing. Libraries were prepared following the Small RNA Sample Preparation Protocol (Illumina) and sequenced on an Illumina Hiseq-2500 sequencing platform using SE50.
Processing of mRNA-Seq Data
To preprocess miRNA-Seq data, adaptor sequences were removed alongside low-quality (<Q20) bases at 5′ and 3′ ends using the software Trimmomatic (v0.30) (Bolger et al., 2014). Reads longer than 70 bp were then used for further experiments. These reads were mapped onto the sugarcane genome (Zhang et al., 2018) using the software Bowtie 2 (2.1.0) (Langmead and Salzberg, 2012) with default parameters set after the preprocessing of mRNA-Seq data. The software TopHat2 was then used to perform a sequence comparison between clean reads and the reference genome to obtain positional information for the reference genome or gene and sequence information unique to the sample. Gene expression levels were presented as fragments per kilo base per million reads (FPKM); thus, genes with FPKM values >1 were retained for further analysis.
A false discovery rate (FDR) <0.05 was used to detect DEGs. After processing with the software DEseq, we applied a more than 2-fold change as the criterion to classify DEGs between two pairs (www.bioconductor.org). Resultant P-values were then adjusted using Benjamini and Hochberg's approaches for controlling the FDR.
Functional annotations were determined using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Default parameters in the software AmiGO were used to obtain GO terms for each gene as well as analyze functional enrichment via hypergeometric tests with FDR correction. This procedure enabled us to obtain an adjusted p < 0.01 between test gene groups and the whole annotation dataset. DEGs in the KEGG pathway were analyzed using the software Cytoscape (Ideker, 2011) using the ClueGO plugin (Wouter et al., 2014).
Processing miRNA-Seq Data
In this step, clean reads were screened from raw data by removing adaptors, poly A sequences, and low-quality bases at both ends of small RNA reads. Reads were then trimmed and cleaned by removing sequences smaller than 18 nt or longer than 40 nt. Values for the sequence duplication level, GC content, Q20, and Q30 in clean data were then calculated. All downstream analyses were based on high quality clean data.
To identify known and novel miRNAs, high-quality clean reads were mapped to the S. spontaneum genome using the software Bowtie 2 (2.1.0) (Langmead and Salzberg, 2012) and searched against NCBI, Rfam, and Repbase databases to remove known classes of RNAs (i.e., mRNA, rRNA, tRNA, snRNA, and snoRNA) and repeats. The remaining reads were then used to detect both known and novel (potentially species-specific) miRNAs among the 5,940 mature miRNAs in plants using mirBase Release 19 (Ana and Sam, 2011).
We performed a differential expression analysis under two conditions using the DESeq R package (1.10.1). This approach provides a statistical analysis for determining differential expression in digital miRNA expression data using a model based on the negative binomial distribution. Thus, miRNAs with an adjusted p < 0.01 were assigned as differentially expressed.
Predicting Potential miRNA Targets: A Joint Analysis of miRNA and mRNA
The psRNATarget tool (http://plantgrn.noble.org/psRNATarget/) and software TargetFinder were used to predict S. spontaneum mRNA targets of putative miRNAs. The software packages BLAST and KOBAS2.0 were then used to compare and carry out a functional annotation of discovered target genes vs. the NR, Swiss-Prot, GO, COG, KOG, Pfam, and KEGG databases. Finally, miRNA and mRNA were related to one another using their interaction in the STRING protein database (http://string-db.org/).
All clean sequencing data used in this study are available in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under accession numbers PRJNA636260 for transcriptomics and PRJNA635765 for miRNA-seq data.
Validation of mRNA and miRNA Relative Expression by QPCR
The RNA used for RNA-seq was used for the QPCR validation of the mRNA expression level. The first-strand cDNA for each sample was prepared from 500 ng of total RNA using SuperScript II reverse transcriptase (Takara, Dalian, China) following the manufacturer's instructions, and template cDNAs were diluted 10-fold for QPCR. Primers were designed using the Primer 5 program and are listed in Supplementary Table 1. Samples and standards were run in triplicate on each plate using the SYBR Premix Ex Taq™ II kit (Takara, Dalian, China) on a StepOne™ real-time PCR system (Applied Biosystems, Foster City, USA) following the manufacturer's recommendations. QPCR was performed in a 20-μL reaction volume containing 6.8 μL of ddH2O, 10 μL of SYBR Premix Ex Taq II, 0.4 μL of ROX Reference Dye II, 0.4 μL of forward primer (10 μmol/L), 0.4 μL of reverse primer (10 μmol/L), and 2 μL of template cDNA. The PCR programs were run as follows: 30 s of pre-denaturation at 95°C, 40 cycles of 5 s at 95°C, and 30 s at 60°C, followed by steps for dissociation curve generation. Dissociation curves for each amplicon were carefully examined to confirm the lack of multiple amplicons at different melting temperatures. Relative transcript levels for each sample were obtained using the comparative cycle threshold method using the cycle threshold value of transcript Sspon.005D0014601 encoding a glyceraldehyde-3-phosphate dehydrogenase was used as a control. The 2−ΔΔCt method was used to calculate the data.
The RNA used for miRNA-Seq was used for the QPCR validation of miRNA expression. No more than 500 ng of each miRNA sample was used for miRNA first-strand preparation according to Stem-loop miRNA cDNA Synthesis kit (EnzyValley, Lot: T435A). Stem-loop primers were listed in Supplementary Table 1. QPCR was performed (Stem-loop miRNA qPCR SYBR kit, Lot: T436A) in a 20-μL reaction volume containing 3.5 μL of ddH2O, 10 μL of 2× Stem-loop miRNA qPCR SYBR ProMix, 0.5 μL of ROX Reference Dye II, 0.5 μL of forward primer (10 μmol/L), 0.5 μL of Uni-SL reverse primer (10 μmol/L), and 5 μL of template cDNA. The PCR programs were run as follows: 10 min of pre-denaturation at 95°C, 40 cycles of 5 s at 95°C, and 30 s at 60°C, followed by steps for dissociation curve generation. The 2−ΔΔCt method was used to calculate the data. The U6 snRNA was taken as the control.
Results
Gene Identification and Their Functional Annotation Analysis in This Study
A total of nine RNA-seq libraries from low temperature (LT)-treated and control (CT) seedlings were constructed for high-throughput sequencing. After the removal of low-quality reads and adapter sequences, 346,617,812 clean paired-end reads were obtained from raw data, with between 33,119,4345 and 44,605,663 obtained from each sample (Supplementary Table 2). Nearly 80% of these paired-end reads could be mapped to the reference genome, and 118,015 potential genes were identified with FPKM values between 0 and 89,290.4 (Sspon.007A0027210) (Supplementary Table 3). In total, 34,052 (28.85%), 56,311 (47.71%), 78,517 (66.53%), 80,524 (68.23%), 87,333 (74.00%), 101,433 (85.95%), and 108,262 (91.73%) genes were annotated in the COG, KOG, Swissprot, Pfam, GO, TrEMBL, and nr databases, respectively. Approximately 10.07% (11,887) of genes were annotated in the KEGG database (Supplementary Tables 3, 4, Supplementary Figure 1), with most genes significantly matched with Sorghum bicolor (70.32%).
Identifying Cold-Responsive Genes in Sugarcane
A total of 6,034 DEGs (10.19% of total genes) were detected via DEGseq analysis (p-value <0.005 and |log2 (fold change)| > 1) under cold stress conditions. Compared with CT, the LT0.5 and LT1 treatments had 2,216 and 4,581 DEGs, respectively, including 1,280 and 936; 2,349 and 2,232 upregulated and downregulated DEGs for each compared pair. When LT1 was compared with LT0.5, 454 genes, including 321 upregulated and 133 downregulated DEGs, were identified (Figure 1A, Supplementary Tables 5–7). The resultant Venn diagram suggests that 83 core DEGs were shared across all compared pairs, and the largest number of unique DEGs was 3,218 in the CT vs. LT1 comparison (Figure 1B). The heat map for all samples indicates that biological repeats reached their expected effect and that four dominant expression profile clades were found via this dynamic hierarchical clustering approach (Figure 2).
Figure 1. DEG statistics for comparisons. (A) Number of downregulated and upregulated DEGs in CT vs. LT0.5, CT vs. LT1, and LT0.5 vs. LT1, respectively. (B) Venn diagram of DEGs distributed in every sample.
Figure 2. Heat map of Saccharum spontaneum genes responding to cold stress at the transcription level. To be considered differentially expressed, a transcript must have fragments per kilobase of transcript per million mapped reads (FPKM) ≥ 0.3 in at least one sample, a 2-fold (or greater) change between samples, and P ≤ 0.05. Red indicates high expression, white indicates intermediate expression, and green indicates low expression.
Functional annotations of all DEGs for each compared pair were presented in Supplementary Tables 5–8 and Supplementary Figure 2. Comparisons between CT and LT0.5 shows that 1,744 DEGs can be assigned to GO terms; of these, the term “regulation of defense response” was the most highly enriched among both upregulated and downregulated DEGs (Figure 3A). In the 1-h cold stress treatment, 3,645 DEGs could be GO annotated and the chloroplast was most affected.
Figure 3. Annotation and functional classification of DEGs. (A) Histogram of most enriched Gene Ontology (GO) classifications of S. spontaneum DEGs under cold stress. The ordinate is the enriched GO term; the abscissa is the enriched q value of the term and the number of differential genes on the column. Different colors are used to distinguish biological processes, cellular components, and molecular functions. (B) Statistical analysis of DEGs mapped to KEGG pathways including carbohydrates, energy, and lipid metabolisms.
The resultant KEGG pathway classification (Supplementary Table 8 and Figure 3B) reveals that most significantly changed transcripts belong to metabolic pathways (especially those for starch, sucrose, pyruvate, and carbon fixation), signal transduction (plant hormones), environmental adaptation (plant–pathogen interaction), and the biosynthesis of secondary metabolites (mainly related to phenylpropanoid biosynthesis).
Protein–Protein Interaction Network Analysis in the mRNA Level Under Cold Stress
We constructed an interaction network based on the DEGs retrieved from known PPI databases. In terms of comparisons between CT and LT0.5, 1,186 DEGs can be linked via 1,100 known relationships, including activation, binding, and inhibition (Supplementary Figure 3A). The gene Sspon.004C0004740 can be identified as the core hub encoding an ATP synthase delta subunit, which was downregulated in the LT0.5 treatment. Sspon.006B0008980 (peptidyl-prolyl cis-trans isomerase CYP38), Sspon.008C0004432 (glycosyl hydrolases), Sspon.001D0013480 (fibrillarin), and Sspon.005A0004900 (ribosomal protein L23) are also important candidate DEGs that function in response to a half an hour cold stress.
After a 1-h cold stress, a total of 2,507 DEGs could be linked via 6,230 known relationships (Figure 4). Significantly, three DEGs (i.e., Sspon.004B0016440, Sspon.004B0016450, and Sspon.004B0015590) encoding the protein phosphatase 2C and the cyclic nucleotide-binding/kinase domain-containing protein (PP2C family) can all be proposed as potential hubs. Similarly, Sspon.003B0011973 (ribosomal protein L13) and Sspon.006B0008980 (peptidyl-prolyl cis-trans isomerase CYP38) also appear to be of importance to the PP2C protein family.
Figure 4. A PPI (protein–protein interaction) network of DEGs from CT vs. LT1. The interaction relationship in the STRING protein interaction database (http://string-db.org/) was then used to analyze the DEG coding protein interaction network. The size of a node in this interaction network is proportional to the degree of this node; thus, the more edges connected to a node, the greater its degree and the larger the node, and the darker the color. Nodes may be in a more central position in the network, and lines of different colors represent varied interactions. The thicker the line, the higher the score value.
In terms of comparisons between CT0.5 and LT1, 230 DEGs could be linked via 1,600 known relationships (Supplementary Figure 3B). Four TIFY family proteins, including Sspon.001B0025400, Sspon.001D0021420, Sspon.001B0025420, and Sspon.001C0021660, were calculated as the core regulated protein that interacted with other genes through 13 linkages. TIFY is a plant-specific family involved in the regulation of plant-specific biologic processes, stress, and responses to phytohormones.
Conserved and Novel miRNAs Identified in Sugarcane in This Study
To identify the miRNAs that respond to cold stress, we constructed nine structural RNA (sRNA) libraries from nine samples. These libraries yielded 110,555,326 clean reads by high-throughput sequencing, which was ~11.05 million clean reads per library (Supplementary Table 9). Mapped sRNAs were then annotated and classified into sRNA, known miRNAs, precursor RNAs, intergenic RNAs, and unknown sRNAs. To verify these results, we analyzed Pearson correlation coefficients of the miRNA expression between biological duplicates (Supplementary Figure 4); all the Pearson values for miRNA expression were >0.82.
We then obtained 20 conserved and 836 novel miRNAs belonging to 13 families from nine libraries (Supplementary Tables 10, 11). Expression profiles of identified miRNAs reveal great variation in TPM values, ranging from 83,065.85 to <10. Four miRNAs (Novel-Chr4B_41904, Novel-Chr4C_44656, Novel-Chr4C_44658, and Novel-tig00009572_93271) belong to the miR396 family and were highly expressed in all samples. Almost of all the identified miRNAs (835 of 856) were shared by all samples, indicating that these were stable in S. spontaneum, even under cold stress. sRNAs range from 18 nucleotides (nt) to 30 nt long, while the majority of mature miRNAs are 21 nt in plants. The distribution of the length of miRNA in S. spontaneum showed that 399 miRNAs were 21 nt long and 277 were 24 nt long (Figure 5A). All 16 known miRNAs were 21 nt long, and the first nucleotide was always U, which is different from the length of novel miRNAs, which were 24 nt long (Supplementary Figure 5).
Figure 5. Statistics of expressed miRNAs in this study. (A) The length distribution of miRNAs identified. (B) Number of downregulated and upregulated differentially expressed miRNAs in CT vs. LT0.5, CT vs. LT1, and LT0.5 vs. LT1. (C) Venn diagram of differentially expressed miRNAs distributed in samples.
Prediction of Targeted mRNAs in Sugarcane
A total of 496 genes were targeted by 423 miRNAs (Supplementary Tables 12, 13). A total of 15 known and 290 novel miRNAs corresponded to 65 and 318 targeted genes, respectively. The miRNA Novel-Chr5B_53523 has the largest number of targeted genes 27. A total of 453 targeted genes (accounting for 91.33%) were annotated in KEGG, GO, and other databases (Supplementary Tables 14–16). In the GO annotation, regulation of transcription (35), oxidation–reduction process (34), lignin catabolic process (14), auxin-activated signaling pathway (12), lignin biosynthesis (11), methylation (11), and cell differentiation processes (10) were enriched in biological processes. The nucleus (60), integral components of membrane (46), apoplast (15), cytoplasm (9), and plasma membrane (8) were enriched in the cellular component, and DNA binding (37), ATP binding (24), copper ion binding (14), oxygen oxidoreductase activity (14), oxidoreductase activity and oxidizing metal ions (14), iron ion binding (13), heme binding (12), lipid binding (11), DNA-binding transcription factor activity (10), O-methyltransferase activity (10), protein dimerization activity (10), and zinc ion binding (10) were the top enriched go terms for molecular function. A total of 251 genes were mapped to 132 pathways within the KEGG database; plant hormone signal transduction (29), stilbenoid diarylheptanoid and gingerol biosynthesis (12), and plant–pathogen interactions (11) were the top three annotated pathways (Supplementary Table 16).
Differentially Expressed miRNAs and Their Targets Under Cold Stress
A total 109 differentially expressed miRNAs exhibiting a more than 2-fold change were identified in S. spontaneum after cold stress, including four known and 105 novel miRNAs, respectively (Supplementary Tables 17–19). Compared with the control, the numbers of differentially expressed miRNAs after cold stress treatment for 0.5 h and 1 h were 53 (38 upregulated and 15 downregulated) and 94 (46 upregulated and 48 downregulated). When LT1 was compared with LT0.5, 28 downregulated miRNAs were identified (Figure 5B). The Venn diagram from this analysis suggests the presence of 13 core miRNAs shared by all compared pairs, including 42 unique miRNAs to CT vs. LT1, the largest number (Figure 5C).
TargetFinder was used to predict candidate targets of differentially expressed miRNAs (Supplementary Tables 20–22). In a comparison of CT and LT0.5, 10 miRNAs targeting 42 genes were identified; of these genes, 14 were targeted by miRNAs Novel-tig00084004_92871, referring to signal transduction, transcription, and energy production, such as MYB and LRR-RLK (leucine-rich repeat receptor-like protein kinase). Besides, there were three LRR-RLK genes (Sspon.006A0015270, Sspon.006B0014350, Sspon.003D0010720, and Sspon.005B0004571) were targeted by miRNA Novel-Chr6B_65233, indicating an important function in the regulation of cold stress signal transduction. Furthermore, a total of four laccase genes (Sspon.003D0006520, Sspon.003B0007420, Sspon.003C0020690, and Sspon.003C0010170) were targeted by Novel-Chr4A_40498. These 18 miRNAs were upregulated by cold stress.
A total of 16 miRNAs targeting 70 genes were identified from the comparison of LT1 and CT. The targeted genes include 14 PPR repeat family members (all targeted by Novel-Chr5C_57213), five LRR-RLKs, six scarecrow-like proteins, and four laccases.
Three miRNAs targeting 18 genes were identified in the comparison between CT0.5 and LT1. Among these, total of 11 genes encoding laccase were targeted by Novel-Chr4C_47059 which was downregulated by cold stress.
Validation of Differentially Expressed mRNAs and miRNAs by Real-time QPCR Analysis
The expression of 10 cold-responsive miRNAs and 25 DEGs was detected by QPCR. The comparisons between QPCR and mRNA-seq (or miRNA-seq) provided in Supplementary Figure 6 suggested that the result of mRNA-seq and miRNA-seq was reliable. This was also confirmed by calculating their Pearson correlation coefficients which were 0.847 and 0.852, respectively (Figure 6).
Figure 6. Correlation analysis of validation. (A) Correlation analysis of gene expression between qPCR and RNA-seq data. (B) Correlation analysis of miRNA expression between qPCR and miRNA-seq data. The X-axis represents RNA-seq (miRNA-seq) data, and the Y-axis represents the qPCR data.
Discussion
Cold stress is a major abiotic environmental factor that negatively affects plants by decelerating seed germination, root development, chlorophyll synthesis, and photosynthesis. Considering this, several protective mechanisms are employed by plants to manage adverse conditions. In sugarcane fields in low latitude and complex mountain and plain environments, cold damage has a serious impact on normal sucrose synthesis. To improve cold resistance in sugarcane, genetics and breeding research has emphasized high-quality cold-resistant varieties. A total of 6,034 DEGs were identified in this analysis as responsive to cold stress (4°C). This number is similar to that of previous reports, which proposed the presence of 2,583 upregulated and 3,302 downregulated genes from a low-temperature (10°C) treatment of S. spontaneum clone IND 00-1037 (Dharshini et al., 2016).
Transcriptional Regulatory Network Underlying Cold Stress at the Seedling Stage
Cold sensing occurs via membrane fluidity, calcium, and lipid signaling genes as well as via MAP kinases and phytohormone signaling (Zhao et al., 2017). Over the past several decades, various reports have confirmed that many transcription factors, including CBF, MYB, CAMAT, and NAC, are either directly or indirectly involved in cold stress tolerance; these were crossed with kinase, hormones, circadian cycles, Ca2+, and light to accomplish cold signal transduction (Peng et al., 2015; Li et al., 2018; Shi et al., 2018). A total of 4% of upregulated genes were related to signaling, while MYB was downregulated 2.4-fold (Dharshini et al., 2016). In maize seedlings under various stress conditions, 43 TF families containing 403 differentially expressed TFs belong to the ERF, MYB, bZIP, bHLH, WRKY, and NAC TF families (Li et al., 2017). In our study, we found that 56 DEGs belong to MYB, bZIP, GATA, and MADS-box TFs in the CT vs. LT0.5 comparison, while 179 other DEGs related to PP2C, CBL-IPK, CDPK, AMPK, calcium, cAMP, MAPK, phosphatidylinositols, plant hormones, and sphingolipids were the top enriched signaling pathways. These results emphasized that TFs, kinase, and plant hormones played vital roles in the cold signal transduction of S. spontaneum.
A total of 14 PP2Cs and 16 PP2Cs in LT0.5 and LT1 compared with CT, respectively, were either downregulated or upregulated under cold stress conditions. Three DEGs (Sspon.004B0016440, Sspon.004B0016450, and Sspon.004B0015590) encoding PP2C proteins are predicted to be potential hubs via string analysis (Figure 4). The first plant PP2C protein (ABI1) was identified as a negative regulator of ABA signaling (Schweighofer et al., 2004). An increasing number of PP2Cs have been defined in different species, and their functions have been subsequently studied. Thus, PP2C proteins were thought to be core members involved in the ABA signaling pathway (Zhu, 2016); however, more recent research showed that the PYL-PP2C-SnRK2 core ABA signaling module activates a MAPK cascade comprised of MAP3Ks MAP3K17/18, MAP2K MKK3, and MAPKs MPK1/2/7/14. These genes might regulate a range of ABA effector proteins via phosphorylation (De Zelicourt et al., 2016) and may also play major roles in the regulation of cell growth as well as cellular biotic and abiotic stress signaling (Lammers and Lavi, 2007; Fan et al., 2019). Proteins belonging to MAPK cascades have also been shown to be involved in the response of plants to cold stress (Markus et al., 2004; Tak et al., 2019), while PTP1 is known to downregulate cold stress (Liu et al., 2016). Thus, we propose that the ABA signaling pathway via PP2C is most important in terms of the cold stress responses of S. spontaneum at the seedling stage; however, future research should be conducted on this pathway.
Research has shown that both the regeneration rate and carboxylation efficiency of RuBP decline under chilling stress (Hardigan et al., 2016; Jin et al., 2017). DEG pathways shown to be enriched in GO and KEGG analyses in CT vs. LT1 have emphasized carbohydrate metabolism including how chloroplasts are related to the photosystem and starch and sucrose biosynthesis. These could also be verified by the functional annotation and PPI analysis (Supplementary Figure 3).
Plant secondary metabolism is closely related to abiotic stress responses (Zhang et al., 2017b). Flavonoids play a role in the response of higher plants to abiotic stressors (Jin et al., 2017). For example, in E. nutans, seven genes involved in phenylpropanoid biosynthesis are known to be specifically related to cold stress response (Fu et al., 2016). In terms of enzymes involved in phenylpropanoid biosynthesis, phenylalanine ammonia lyase (PAL) is one of the most relevant (Vogt, 2010). The RNA-seq data presented here show that the number of DEGs related to phenylpropanoid pathways accounts for the largest proportion of secondary metabolic pathways. Therefore, we inferred that lignin might play a key adaptive role in cold stress response.
Cold Stress-Responsive miRNAs: Sugarcane Regulatory Network
Yang et al. (2017) found that 4,002 targets were predicted for 259 miRNAs, which may be due to the absence of a sugarcane genome leading to inaccuracies in prediction. Using the published sugarcane genome as our starting point, we found 365 genes from 305 miRNAs, which is similar to reports on Astragalus membranaceus (Abla et al., 2019), Nelumbo nucifera (Zou et al., 2020), potato tubers (Ou et al., 2014), sweet potato (Xie et al., 2017), and other plants (Gupta et al., 2014).
Several previous studies on various plant species have confirmed the involvement of miRNAs in cold stress responses (Tang and Chu, 2017; Zeng et al., 2018; Song et al., 2019; Sun et al., 2019; Esposito et al., 2020). A total of 412 sugarcane miRNAs have been isolated from two cultivars, including 261 known and 151 novel miRNAs. These include nodes for ROC22 (relative cold sensitivity) and FN39 (relative cold tolerance); data indicate that 62 of these miRNAs exhibit significant cold stress-induced or depressed expression (Yang et al., 2017), which is less than that found in our study. We found that 109 miRNAs belonging to 11 families responded to cold stress, while miR444 had the most differentially expressed miRNAs (52) followed by miR169-1 (18) (22) and miR396 (13). Compared to previous studies, we also found that miR156, miR168, and miR408 are upregulated under cold stress in sugarcane (Yang et al., 2017). The miR444 in rice has been reported as a key factor relaying antiviral signaling from virus infections to OsRDR1 expression (Wang et al., 2016a) as well as the promotion of BR biosynthesis (Jiao et al., 2019). This gene is also involved in the nitrate signaling pathway (Yan et al., 2014). Increased expression of miR408 could also lead to improved Arabidopsis tolerance to cold and oxidative stress, as shown by a lower ROS level and the induction of genes related to antioxidative function (Ma et al., 2015).
We found 28, 44, and 13 miR444s in CT vs. LT0.5, CT vs. LT1, and CT0.5 vs. LT1, respectively. This is especially the case for Novel-Chr5C_55773, which exhibited induced expression after a 1-h cold stress and targeted three DEGs encoding the apyrase 3 family (Sspon.005D0023491, Sspon.005B0020100, and Sspon.005C0019251). It has been shown that these can maintain plasma membrane integrity and reduce electrolyte leakage under cold stress (Deng et al., 2015). Four LRR-RLKs were also targeted by miR444 (Novel-Chr6B_65223) and miR168 (Novel-Chr5D_60023) after a 1-h cold stress. We propose that at an early stage, cold stress reduces plasma membrane integrity by decreasing the expression of genes encoding apyrase 3 through the regulation of miR444, which also acts as membrane location receptors. These four LRR-RLK genes play an important function in cold stress signal transduction.
An integrated miRNA–gene/mRNA regulation network suggests that 25 miRNAs belong to nine families that play important roles in response to cold stress. Among these nine families, Novel-tig00084004_47811 (miR169), Novel-Chr4C_24049 (miR396), Novel-Chr5C_29143 (miR396), and Novel-tig00005469_46467 (miR159) were significant regulators. Specifically, miR169 and miR396 have been reported to play an indirect role in cold stress response by regulating the TFs (Wang et al., 2017; Sombir et al., 2020). The OsGRF4 transcription factor is targeted by miRNA396, whereas mutations disturb the GRF4-miR396 stress response network and results in the development of enlarged grains and cold tolerance enhancements in rice (Chen et al., 2019). Among 14 targets, two genes encoding MYB TFs are thought to be involved in transcriptional regulation. LRR-RLKs and phosphatidylinositol-4-phosphate 5-kinase have also been identified as regulators that function in early cold stress signal transduction.
Gene miR396 in Poncirus trifoliata functions in cold tolerance by regulating acc oxidase gene expression and modulating ethylene-polyamine homeostasis (Zhang et al., 2016). This gene also acts as a stress-responsive regulator by conferring tolerance to abiotic stressors and susceptibility to mold infection in tobacco (Chen et al., 2015). A total of 17 genes that encode pentatricopeptide repeat (PPR) family proteins were targeted by Novel-Chr5C_29143 (miR396); these genes play critical roles in all aspects of organelle RNA metabolism including development and stress defense (Zhang et al., 2017a). Another miR396 family member, Novel-Chr4C_24049, was upregulated after a 0.5-h cold stress and is thought to regulate 11 laccases−13 genes, especially Sspon.003B0004220, Sspon.003A0003410, and Sspon.007A0025870. These 11 laccases−13 genes are also differentially expressed under cold stress. Our results suggested that MiR444, miR169, and miR396 play key roles in the regulation of signal transduction, transcriptional, and lignin metabolism in sugarcane responding to cold stress.
Conclusions
The results of this study show that the multiple regulated mechanisms are involved in the response of sugarcane to cold stress. The signaling pathway mediated by PP2C and LRR-RLK is thought to play a key role in the response of sugarcane to cold stress and has not been previously reported. Most importantly, genes involved in the biosynthesis of polyunsaturated fatty acids and LEA protein, lignin, and pectin pathways were found to be directly involved in cold stress responses targeted by miR444, miR169, and miR396. These key candidate genes and miRNAs will be useful for future research on cold stress response mechanisms. This study provides new insight into the regulatory network at transcriptional and post-transcriptional levels, thus allowing us to better understand the response of sugarcane to cold stress.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
YLia, XS, and BZ: data curation. BZ, ZQ, and DL: formal analysis. RC, ZZ, and YD: investigation. YLia and JWu: funding acquisition. JWe and JWu: project administration. XH and JWu: writing—original draft. XH, YLi, and XS: writing—review and editing. All authors approved the final manuscript.
Funding
This research was supported by grants from the National Natural Science Foundation of China (31660356, 31901594, and 31360312), the Development and Regulation of Economically Important Traits in Tropical Crops (2018YFD1000500), the Guangxi Post-doctoral Science Foundation, Guangxi Science and Technology Base and Talents Special Project (GuikeAD17195100), the Fund for Guangxi Innovation Teams of Modern Agriculture Technology (nycytxgxcxtd-03-01), and the Fund Science and Technology Development of Guangxi Academy of Agricultural Sciences (Guinongke2021YT013, Guinongke2017JZ19, and Guinongke2016JM06).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.581993/full#supplementary-material
Supplementary Figure 1. Functional annotation of all genes isolated from S. spontaneum in this study. (A) Statistics for the functional annotation for all identified genes. (B) GO assignment for all identified genes. (C) Statistics of the NR annotation results mapped to other plants. (D) Unigenes were aligned to the COG (Clusters of Orthologous Groups) database to predict and categorize possible functions.
Supplementary Figure 2. Annotation statistics for DEGs in the second level of the GO database (A–C) and Annotation statistics for functional categories for DEGs from the COG database (D,E). (A,D) Statistics of DEGs from CT vs. LT0.5. (B,E) Statistics of DEGs from CT vs. LT1. (C,F) Statistics of DEGs from CT0.5 vs. LT1.
Supplementary Figure 3. DEG PPI network. Interaction relationship in the STRING protein interaction database (http://string-db.org/) to analyze DEGs coding in the protein interaction network. (A) PPI map for CT vs. LT0.5. (B) Map of PPI for CT 0.5 vs. LT1. The size of a node in this interaction network is proportional to the degree of this node; thus, the more edges connected to a node, the greater its degree, the larger the node, and the darker the color. Nodes may be in a more central position in the network and lines of different colors represent varied interactions. The thicker the line, the higher the score.
Supplementary Figure 4. Analysis of miRNA expression correlation among samples and differential expression clustering. (A) Clusters of expression patterns of differentially expressed miRNAs. The abscissa represents the sample name and the clustering result of the sample, while the ordinate represents differential miRNA and miRNA clustering results. Different columns represent different samples, while different rows represent different miRNAs. The color in this Figure represents miRNA expression level log2 (TPM) in the sample after normalization (scale number). (B) Correlation heat map of the expression levels of two samples. Pearson's Correlation Coefficient was used as an evaluation index for biological repeat correlations. The closer r is to 1, the stronger the correlation between the two replicate samples.
Supplementary Figure 5. Base preferences for miRNA. (A) Base bias of known miRNAs. (B) Base bias of novel miRNAs. The abscissa is the base position of miRNA while the ordinate is the percentage of bases A/U/C/G appearing in the miRNA at that position.
Supplementary Figure 6. QPCR validation for mRNA and miRNA. (A–C) The validation of mRNA in CT vs. LT0.5, CT vs. LT1, and LT0.5 vs. LT1. (D–F) The validation of miRNA in CT vs. LT0.5, CT vs. LT1, and LT0.5 vs. LT1. The Y-axis represents the value of relative expression which is conversed by Log2 (Fold change), the X-axis represents the ID of mRNA or miRNA.
Supplementary Table 1. List of primers used for QPCR validation.
Supplementary Table 2. Statistics of mRNA-seq data and the comparison between sequencing data and reference genome.
Supplementary Table 3. The FPKM value and Functional annotation of all identified genes from S. spontaneum in this study.
Supplementary Table 4. Statistics of functional annotation for all identified genes in this study.
Supplementary Table 5. List of DEGs in comparison between CT and LT0.5.
Supplementary Table 6. List of DEGs in comparison between CT and LT1.
Supplementary Table 7. List of DEGs in comparison between LT0.5 and LT1.
Supplementary Table 8. Statistics of KEGG annotation in comparisons.
Supplementary Table 9. Statistics of sRNA-seq reads.
Supplementary Table 10. The family and TPM value of every miRNA in each sample.
Supplementary Table 11. Statistics of miRNAs in each sample.
Supplementary Table 12. Statistics of miRNA with target and target genes.
Supplementary Table 13. List of target genes for every miRNA.
Supplementary Table 14. The function annotation of target genes.
Supplementary Table 15. List of annotation for target genes against GO database.
Supplementary Table 16. List of annotation for target genes against KEGG database.
Supplementary Table 17. The differentially expressed miRNAs in comparison of CT vs. LT 0.5.
Supplementary Table 18. The differentially expressed miRNAs in comparison of CT vs. LT1.
Supplementary Table 19. The differentially expressed miRNAs in comparison of LT 0.5 vs. LT 1.
Supplementary Table 20. Differentially expressed miRNAs under cold stress and their TPM value.
Supplementary Table 21. Corresponding miRNA-mRNA and the functional annotation of mRNA from CT vs. LT0.5.
Supplementary Table 22. Corresponding miRNA-mRNA and the functional annotation of mRNA from CT vs. LT1.
Abbreviations
MEKK, mitogen-activated protein/ERK kinase kinases; MPK, MAPK phosphatase; CBF, C-repeat cis-element binding factor; COR, cold-regulated; CRPK, calcium-regulated protein kinase; GSH, glutathione; ROS, reactive oxygen species; ERF, ethylene response factor; ARF, auxin response factors; ABA, abscisic acid; SOD, superoxide dismutase; MDA, malondialdehyde; CT, control; LT, low temperature; TF, transcriptional factor; ESTs, expressed sequence tags; FPKM, number of fragments per kilobase of transcript sequence per million base pairs sequenced; COG, clusters of orthologous groups of proteins; PP2C, protein phosphatase 2C; KOG, eukaryotic ortholog groups; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology; DEG, differentially expressed gene; LRR-RLK, leucine rich repeats-receptor-like protein kinase; LEA, late embryogenesis abundant protein; TPM, transcripts per million.
References
Abla, M., Sun, H., Li, Z., Wei, C., Gao, F., Zhou, Y., et al. (2019). Identification of miRNAs and their response to cold stress in Astragalus membranaceus. Biomolecules 9:182. doi: 10.3390/biom9050182
Ana, K., and Sam, G.-J. (2011). miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 39(Database issue), D152–D157. doi: 10.1093/nar/gkq1027
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Chen, L., Luan, Y., and Zhai, J. (2015). Sp-miR396a-5p acts as a stress-responsive genes regulator by conferring tolerance to abiotic stresses and susceptibility to Phytophthora nicotianae infection in transgenic tobacco. Plant Cell Rep. 34, 2013–2025. doi: 10.1007/s00299-015-1847-0
Chen, X., Jiang, L., Zheng, J., Chen, F., Wang, T., Wang, M., et al. (2019). A missense mutation in large grain size 1 increases grain size and enhances cold tolerance in rice. J. Exp. Bot. 70, 3851–3866. doi: 10.1093/jxb/erz192
De Zelicourt, A., Colcombet, J., and Hirt, H. (2016). The role of MAPK modules and ABA during abiotic stress signaling. Trends Plant Sci. 21, 677–685. doi: 10.1016/j.tplants.2016.04.004
Deng, S., Sun, J., Zhao, R., Ding, M., and Chen, S. (2015). Populus euphratica APYRASE2 enhances cold tolerance by modulating vesicular trafficking and extracellular atp in Arabidopsis plants. Plant Physiol. 169, 530–548. doi: 10.1104/pp.15.00581
Dharshini, S., Chakravarthi, M., Ashwin Narayan, J, Manoj, V. M., Naveenarani, M., Kumar, R., et al. (2016). De novo sequencing and transcriptome analysis of a low temperature tolerant Saccharum spontaneum clone IND 00-1037. J. Biotechnol. 231, 280–294. doi: 10.1016/j.jbiotec.2016.05.036
Ding, C., Yang, Q., Li, F., and Weifu, L. (2006). Effect of cold stress on the content of free-proline in leaf of Saccharum spontaneum L. and sclerostachya (II). J. Anhui Agric. Sci. 34, 30–33. doi: 10.13989/j.cnki.0517-6611.2006.05.010
Esposito, S., Aversano, R., Bradeen, J. M., Di Matteo, A., Villano, C., and Carputo, D. (2020). Deep-sequencing of Solanum commersonii small RNA libraries reveals riboregulators involved in cold stress response. Plant Biol. 22, 133–142. doi: 10.1111/plb.12955
Fan, K., Yuan, S., Chen, J., Chen, Y., Li, Z., Lin, W., et al. (2019). Molecular evolution and lineage-specific expansion of the PP2C family in Zea mays. Planta 250, 1521–1538. doi: 10.1007/s00425-019-03243-x
Fu, J., Miao, Y., Shao, L., Hu, T., and Yang, P. (2016). De novo transcriptome sequencing and gene expression profiling of Elymus nutans under cold stress. BMC Genomics 17:870. doi: 10.1186/s12864-016-3222-0
Guo, X., Liu, D., and Kang, C. (2018). Cold signaling in plants: insights into mechanisms and regulation. J. Integr. Plant Biol. 60, 745–756. doi: 10.1111/jipb.12706
Gupta, O., Meena, N., Sharma, I., and Sharma, P. (2014). Differential regulation of microRNAs in response to osmotic, salt and cold stresses in wheat. Mol. Biol. Rep. 41, 4623–4629. doi: 10.1007/s11033-014-3333-0
Hardigan, M. A., Crisovan, E., Hamilton, J. P., Kim, J., Laimbeer, P., Leisner, C. P., et al. (2016). Genome reduction uncovers a large dispensable genome and adaptive role for copy number variation in asexually propagated Solanum tuberosum. Plant Cell 28, 388–405. doi: 10.1105/tpc.15.00538
Ideker, T. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27, 431–432. doi: 10.1093/bioinformatics/btq675
Jiao, X., Wang, H., Yan, J., Kong, X., and Yan, Y. (2019). Promotion of BR biosynthesis by miR444 is required for ammonium-triggered inhibition of root growth. Plant Physiol. 182, 1454–1466. doi: 10.1104/pp.19.00190
Jin, J., Zhang, H., Zhang, J., Liu, P., Chen, X., Li, Z., et al. (2017). Integrated transcriptomics and metabolomics analysis to characterize cold stress responses in Nicotiana tabacum. BMC Genomics 18:496. doi: 10.1186/s12864-017-3871-7
Jong-Won, P., Benatti, T. R., Thiago, M., Yu, Q., Nora, S. G., Victoria, M., et al. (2015). Cold responsive gene expression profiling of sugarcane and Saccharum spontaneum with functional analysis of a cold inducible saccharum homolog of NOD26-Like intrinsic protein to salt and water stress. PLoS ONE 10:e0125810. doi: 10.1371/journal.pone.0125810
Lammers, T., and Lavi, S. (2007). Role of type 2C protein phosphatases in growth regulation and in cellular stress signaling. Crit. Rev. Biochem. Mol. Biol. 42, 437–461. doi: 10.1080/10409230701693342
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Li, H.-B., Li, N., Yang, S.-Z., Peng, H.-Z., Wang, L.-L., Wang, Y., et al. (2016). Transcriptomic analysis of Casuarina equisetifolia L. in responses to cold stress. Tree Genet. Genomes 13:7. doi: 10.1007/s11295-016-1090-z
Li, P., Cao, W., Fang, H., Xu, S., Yin, S., Zhang, Y., et al. (2017). Transcriptomic profiling of the maize (Zea mays L.) leaf response to abiotic stresses at the seedling stage. Front. Plant Sci. 8:290. doi: 10.3389/fpls.2017.00290
Li, Q., Byrns, B., Badawi, M. A., Diallo, A. B., Danyluk, J., Sarhan, F., et al. (2018). Transcriptomic insights into phenological development and cold tolerance of wheat grown in the field. Plant Physiol. 176, 2376–2394. doi: 10.1104/pp.17.01311
Liu, S., Chen, H., Li, X., and Zhang, W. (2016). A low-temperature-responsive element involved in the regulation of the Arabidopsis thaliana At1g71850/At1g71860 divergent gene pair. Plant Cell Rep. 35, 1757–1767. doi: 10.1007/s00299-016-1994-y
Ma, C., Burd, S., and Lers, A. (2015). miR408 is involved in abiotic stress responses in Arabidopsis. Plant J. 84, 169–187. doi: 10.1111/tpj.12999
Markus, T., Elisabeth, S., Thomas, E., Dóczi, R., Kazuya, I., Kazuo, S., et al. (2004). The MKK2 pathway mediates cold and salt stress signaling in Arabidopsis. Mol. Cell 15, 141–152. doi: 10.1016/j.molcel.2004.06.023
Megha, S., Basu, U., and Kav, N. N. V. (2018). Regulation of low temperature stress in plants by microRNAs. Plant Cell Environ. 41, 1–15. doi: 10.1111/pce.12956
Nogueira, F. T., De Rosa, V. E. Jr, Menossi, M., Ulian, E. C., and Arruda, P. (2003). RNA expression profiles and data mining of sugarcane response to low temperature. Plant Physiol. 132, 1811–1824. doi: 10.1104/pp.102.017483
Ou, Y., Liu, X., Xie, C., Zhang, H., Lin, Y., Li, M., et al. (2014). Genome-wide identification of microRNAs and their targets in cold-stored potato tubers by deep sequencing and degradome analysis. Plant Mol. Biol. Rep. 33, 584–597. doi: 10.1007/s11105-014-0771-8
Pang, X., Zhu, P., Zhou, Q., Huang, Q., Tan, Q., and Chun, L. (2015). Research progress on sugarcane response to cold stress. Guizhou Agric. Sci. 43, 39–43.
Peng, X., Wu, Q., Teng, L., Tang, F., Pi, Z., and Shen, S. (2015). Transcriptional regulation of the paper mulberry under cold stress as revealed by a comprehensive analysis of transcription factors. BMC Plant Biol. 15:108. doi: 10.1186/s12870-015-0489-2
Rasheed, R., Wahid, A., Ashraf, M., and Basra, S. M. A. (2010). Role of proline and glycinebetaine in improving chilling stress tolerance in sugarcane buds at sprouting. Int. J. Agric. Biol. 12, 1560–8530. doi: 10.3763/ijas.2009.0478
Schweighofer, A., Hirt, H., and Meskiene, I. (2004). Plant PP2C phosphatases: emerging functions in stress signaling. Trends Plant Sci. 9, 236–243. doi: 10.1016/j.tplants.2004.03.007
Selvarajan, D., Mohan, C., Dhandapani, V., Nerkar, G., Jayanarayanan, A. N., Vadakkancherry Mohanan, M., et al. (2018). Differential gene expression profiling through transcriptome approach of Saccharum spontaneum L. under low temperature stress reveals genes potentially involved in cold acclimation. 3 Biotech 8:195. doi: 10.1007/s13205-018-1194-2
Shi, Y., Ding, Y., and Yang, S. (2018). Molecular regulation of cbf signaling in cold acclimation. Trends Plant Sci. 23, 623–637. doi: 10.1016/j.tplants.2018.04.002
Sombir, R., Sonia, B., Sarita, J., and Saloni, M. (2020). Novel insights into expansion and functional diversification of MIR169 family in tomato. Planta 251, 1–17. doi: 10.1007/s00425-020-03346-w
Song, X., Li, Y., Cao, X., and Qi, Y. (2019). MicroRNAs and their regulatory roles in plant-environment interactions. Annu. Rev. Plant Biol. 70, 489–525. doi: 10.1146/annurev-arplant-050718-100334
Sternes, P. R., and Moyle, R. L. (2015). Deep sequencing reveals divergent expression patterns within the small RNA transcriptomes of cultured and vegetative tissues of sugarcane. Plant Mol. Biol. Rep. 33, 931–951. doi: 10.1007/s11105-014-0787-0
Sun, X., Lin, L., and Sui, N. (2019). Regulation mechanism of microRNA in plant response to abiotic stress and breeding. Mol. Biol. Rep. 46, 1447–1457. doi: 10.1007/s11033-018-4511-2
Tak, H., Negi, S., Rajpurohit, Y. S., Misra, H. S., and Ganapathi, T. R. (2019). MusaMPK5, a mitogen activated protein kinase is involved in regulation of cold tolerance in banana. Plant Physiol. Biochem. 146, 112–123. doi: 10.1016/j.plaphy.2019.11.012
Tang, J., and Chu, C. (2017). MicroRNAs in crop improvement: fine-tuners for complex traits. Nat. Plants 3:17077. doi: 10.1038/nplants.2017.77
Thiebaut, F., Rojas, C. A., Almeida, K. L., Grativol, C., Domiciano, G. C., Lamb, C. R. C., et al. (2012). Regulation of miR319 during cold stress in sugarcane. Plant Cell Environ. 35, 502–512. doi: 10.1111/j.1365-3040.2011.02430.x
Wang, H., Jiao, X., Kong, X., Humaira, S., Wu, Y., Chen, X., et al. (2016a). A signaling cascade from miR444 to RDR1 in rice antiviral RNA silencing pathway. Plant Physiol. 170, 2365–2377. doi: 10.1104/pp.15.01283
Wang, L. L., Zhao, T. L., Jin-Tao, G. E., and Liu, X. M. (2017). Application prospects of plant cold-stress-responsive miRNAs in cold resistance research of plants. Acta Agric. Shanghai 33, 129–134. doi: 10.15955/j.issn1000-3924.2017.06.25
Wang, Q., Liu, N., Yang, X., Tu, L., and Zhang, X. (2016b). Small RNA-mediated responses to low- and high-temperature stresses in cotton. Sci. Rep. 6:35558. doi: 10.1038/srep35558
Wang, S.-T., Sun, X.-L., Hoshino, Y., Yu, Y., Jia, B., Sun, Z.-W., et al. (2014). MicroRNA319 positively regulates cold tolerance by targeting OsPCF6 and OsTCP21 in rice (Oryza sativa L.). PLoS ONE 9:e91357. doi: 10.1371/journal.pone.0091357
Wouter, W., Carmen, A. A., Gertjan, K., Johannes, P. V., Marrije, R. B., Gemma, G. K., et al. (2014). Label-free LC-MSe in tissue and serum reveals protein networks underlying differences between benign and malignant serous ovarian tumors. PLoS ONE 9:e108046. doi: 10.1371/journal.pone.0108046
Xie, Z., Wang, A., Li, H., Yu, J., Jiang, J., Tang, Z., et al. (2017). High throughput deep sequencing reveals the important roles of microRNAs during sweetpotato storage at chilling temperature. Sci. Rep. 7:16578. doi: 10.1038/s41598-017-16871-8
Yan, Y., Wang, H., Hamera, S., Chen, X., and Fang, R. (2014). miR444a has multiple functions in the rice nitrate-signaling pathway. Plant J. 78, 44–55. doi: 10.1111/tpj.12446
Yang, C., Li, D., Mao, D., Liu, X., Ji, C., Li, X., et al. (2013). Overexpression of microRNA319 impacts leaf morphogenesis and leads to enhanced cold tolerance in rice (Oryza sativa L.). Plant Cell Environ. 36, 2207–2218. doi: 10.1111/pce.12130
Yang, Y., Zhang, X., Su, Y., Zou, J., Wang, Z., Xu, L., et al. (2017). miRNA alteration is an important mechanism in sugarcane response to low-temperature environment. BMC Genomics 18:833. doi: 10.1186/s12864-017-4231-3
Zeng, X., Xu, Y., Jiang, J., Zhang, F., Ma, L., Wu, D., et al. (2018). Identification of cold stress responsive microRNAs in two winter turnip rape (Brassica rapa L.) by high throughput sequencing. BMC Plant Biol. 18:52. doi: 10.1186/s12870-018-1242-4
Zhang, J., Xiao, J., Li, Y., Su, B., Xu, H., Shan, X., et al. (2017a). PDM3, a pentatricopeptide repeat-containing protein, affects chloroplast development. J. Exp. Bot. 68, 5615–5627. doi: 10.1093/jxb/erx360
Zhang, J., Xu, Y., Huan, Q., and Chong, K. (2009). Deep sequencing of brachypodium small RNAs at the global genome level identifies microRNAs involved in cold stress response. BMC Genomics 10:449. doi: 10.1186/1471-2164-10-449
Zhang, J., Zhang, X., Tang, H., Zhang, Q., et al. (2018). Allele-defined genome of the autopolyploid sugarcane Saccharum spontaneum L. Nat. Genet. 50, 1565–1573. doi: 10.1038/s41588-018-0237-2
Zhang, X., Teixeira Da Silva, J. A., Niu, M., Li, M., He, C., Zhao, J., et al. (2017b). Physiological and transcriptomic analyses reveal a response mechanism to cold stress in Santalum album L. leaves. Sci. Rep. 7:42165. doi: 10.1038/srep42165
Zhang, X., Wei, W., Ming, W., Zhang, H.-Y., and Liu, J.-H. (2016). The miR396b of poncirus trifoliata functions in cold tolerance by regulating ACC oxidase gene expression and modulating ethylene–polyamine homeostasis. Plant Cell Physiol. 57, 1865–78. doi: 10.1093/pcp/pcw108
Zhao, C., Wang, P., Si, T., Hsu, C. C., Wang, L., Zayed, O., et al. (2017). MAP kinase cascades regulate the cold response by modulating ICE1 protein stability. Dev. Cell 43, 618–629. doi: 10.1016/j.devcel.2017.09.024
Zhu, J. K. (2016). Abiotic stress signaling and responses in plants. Cell 167, 313–324. doi: 10.1016/j.cell.2016.08.029
Zhu, P., Pang, X., Tan, Q., Liang, C., Yan, L., Zhou, Q., et al. (2019). Effects of chilling on photosynthesis and photosynthetic pigment contents in leaves of different sugarcane varieties. Chin. J. Trop. Crops 40, 51–57.
Keywords: miRNAs, integrative analysis, cold stress, sugarcane, miRNA targeted genes
Citation: Huang X, Liang Y, Zhang B, Song X, Li Y, Qin Z, Li D, Chen R, Zhou Z, Deng Y, Wei J and Wu J (2021) Integration of Transcriptional and Post-transcriptional Analysis Revealed the Early Response Mechanism of Sugarcane to Cold Stress. Front. Genet. 11:581993. doi: 10.3389/fgene.2020.581993
Received: 10 July 2020; Accepted: 11 December 2020;
Published: 25 January 2021.
Edited by:
Rong Zhou, Aarhus University, DenmarkReviewed by:
Kranthi Varala, Purdue University, United StatesHimanshu Sharma, National Agri-Food Biotechnology Institute, India
Paula Barbim Donate, University of São Paulo, Brazil
Copyright © 2021 Huang, Liang, Zhang, Song, Li, Qin, Li, Chen, Zhou, Deng, Wei and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jiguang Wei, amlndWFuZ3dlaSYjeDAwMDQwO2d4dS5lZHUuY24=; Jianming Wu, d3VqaWFubWluZzIwMDQmI3gwMDA0MDsxMjYuY29t
†These authors have contributed equally to this work