Transcriptome and Degradome Profiling Reveals a Role of miR530 in the Circadian Regulation of Gene Expression in Kalanchoë marnieriana

Mature plants of K. marnieriana were grown in the greenhouse of the Research Institute of Subtropical Forestry (119°57′ N, 30°04′ E; Fuyang, Zhejiang, China). Tissue was collected from propagated plants of K. marnieriana grown in a growth chamber under photoperiod of 12 h light/12 h dark at 24 °C and 40% humidity. Additionally, K. marnieriana plants were subjected to the drought treatments: 10 days without watering (D10), 25 days without watering (D25) and well-watered as the control condition. For the Agrobacterium-infiltration assay, the Nicotiana benthamiana plants were grown at 22 °C under 16 h light/8 h dark condition; one-month-old plants were used for the infiltration as described by [31].

For RNA preparation, the plant tissues were collected and frozen in liquid nitrogen and stored at −80 °C before subsequent uses. Total RNA was extracted using an RNAprep Pure Plant Kit (DP441, Tiangen, Beijing, China) from the mesophyll, epidermal peel and leaves of K. marnieriana. The concentration and integrity of the total RNA were checked before library construction. A Nanodrop 2000 spectrophotometer (Thermo Fisher, CA, USA) was used to calculate the RNA concentration, and samples with more than 200 ng/uL and optical density (OD) 260/280 above 2.0 were used for small RNA sequencing. Near equal amounts of RNAs from different tissue types were mixed to construct libraries for transcriptome assembly and degradome analysis.

The small RNA sequencing libraries were prepared by a TruSeq Small RNA Library Preparation Kits (Illumina, San Diego, CA, USA) according to user’s manual from Illumina as described by [32]. Briefly, following the isolation of suitable cDNA fragments by gel electrophoresis analysis, the obtained libraries were enriched with small RNAs and loaded onto an Illumina HiSeqTM 2000 platform for sequencing by Hangzhou LC-Bio Co., Ltd. (Hangzhou, China). Total reads were trimmed of low-quality and adapter sequences using the Illumina Pipeline. The mixed RNA sample was sequenced and assembled into nonredundant unigenes using Trinity2.1.1.0 (Cambridge, MA, USA) [33]. Unigenes were tentatively identified based on the best hits against known sequences in the database. The transcriptome data of K. marnieriana were available in National Center for Biotechnology Information (TSA accession number GIXV00000000↗). All sequencing reads were deposited into the NCBI SRA database under Bioproject PRJNA684436. The identification of miRNA was initially identified according previous pipelines [34]. Genome-wide prediction of miRNAs was performed using the Shortstack3.4 (University Park, PA, USA) [35] and combined results were obtained by removing redundant miRNAs with precursor sequences. To quantify the abundance of miRNA, a transcripts per million value was defined as ‘counts of read mapped to miRNA * 1,000,000’/‘reads mapped to reference genome’, followed by differential expression analysis as described by [36]. For the degradome analysis, the extracted sequencing reads over 20 nucleotides were used to identify potentially cleaved targets by the CleaveLand 3.0 pipeline (University Park, PA, USA), and p-value of less than 0.05 was used to identify highly confident targets [37]. The Kalanchoë genome references were used to identify the annotation of target genes.

To validate the miRNA target genes, the predicted target site and miRNA precursor sequences were constructed into the Dual-Luciferase Sensors system (Addgene, 55206, Watertown, MA, USA) [31]. All constructs were verified by sequencing. The miR530 precursor was obtained by amplifying the K. marnieriana genomic DNA using the sequence specific primers (Supplementary Table S1). Then, the miRNA overexpression construct was coinfiltrated with the target sites vector in tobacco accordingly [31]. For each combination, including control and test groups, at least five biological replicates were used for efficiency analysis [38].

The expression levels of miRNAs were examined by qRT-PCR in conjunction with a Mir-X™ miRNA FirstStrand Synthesis kit (Cat. 638313, TaKaRa, Dalian, China) according to the user’s manual. The target gene expression analysis was performed using a PrimeScript kit (Cat. RR037Q, TaKaRa, Dalian, China). The gene-specific primers (Supplementary Table S1) were designed by Primer Express 3.0.1 (Applied Biosystems, Foster City, CA, USA). The RT-qPCR was performed on an ABI PRISM 7300 Real-Time PCR System (Foster City, CA, USA). For qRT-PCR experiments, three biological replicates were used, and each was repeated three times as technical replicates. The data were analyzed with the 2–ΔΔCT method [39].

K. marnieriana is a closely related species to K. fedtschenkoi and K. laxiflora according to its morphology and biological characteristics [40] (Supplementary Figure S1). Using the sequence of 18S rRNA, we determined that this accession was closely related to other two reported K. marnieriana plants (Supplementary Figure S1). To assess the drought responses in K. marnieriana, we performed the acidity titration measurements between dawn and dusk under different watering conditions. We found that, under the regularly watering condition, the difference value (ΔH+) of titratable acid in mature leaves of K. marnieriana averaged 72 nmol H+ per kg, and the acidity levels were significantly increased under the 10 days (D10) and 25 days (D25) of drought treatments (Supplementary Figure S1C), indicating that the drought treatments enhanced the production of leaf titratable acids.

To reveal the miRNAs and their targets that were responsive to drought, we designed an experiment of drought treatment by withholding water. Leaf epidermis and mesophyll tissues were separated and collected for sampling at the D10 and D25 conditions (Figure 1). We performed a combination of transcripts and small RNAs sequencing analyses and focused on the miRNAs and their target genes. The K. laxiflora genome and transcriptome assembly were both used as the references for miRNA gene identification. The transcriptome assembly was obtained by using the mixed RNA samples contained 26,336 unigenes with the N50 of 1434 bp (TSA GIXV00000000↗). We consolidated the prediction results to remove redundant miRNAs. In total, we obtained 94 miRNAs, of which 29 were supported by the genome assembly of K. fedtschenkoi (Figure 1; Supplementary Dataset S1). To predict the targets of miRNA, we performed the Parallel Analysis of RNA Ends (PARE, also known as degradome sequencing) to reveal the potential targeted transcripts of miRNAs in K. marnieriana (Figure 1). The potential targets of miRNA were identified based on the frequency of RNA ends. We determined that the integrative analyses of differential expression of miRNAs and their targets were sufficient to inform the CAM-related regulation of gene expression in Kalanchoë.

To uncover the miRNAs, we filtered the low-quality sequences and obtained 18–25 bp small RNA (sRNA) sequencing reads for the bioinformatics identification pipeline (Supplementary Table S2). For the miRNA identification from transcriptome assembly, the secondary structure and the reads-mapping profile were evaluated to generate the high-confident miRNAs. We found that the 29 genome-supported miRNA displayed a more uniform precursor length and higher expression levels comparing to the transcriptome-based miRNAs (Figure 2A). The mature sequences of miRNAs were annotated using plant miRNAs databases. We discovered that eight miRNAs from the 29 genome-supported miRNAs were not homologous to other plant miRNAs and were designated as Kalanchoë-specific miRNAs (Figure 2B). To further examine these Kalanchoë-specific miRNAs, we analyzed the secondary structure and sRNA coverage of precursors; all eight miRNAs displayed canonical features of miRNA genes (Figure 2B).

To identify miRNA genes responsive to drought, we performed expression analysis to reveal differentially expressed miRNAs. The expression level of mature miRNAs at each condition was quantified (Supplementary Dataset S2) and used for subsequent statistical analysis. We found that 55 miRNAs were differentially expressed in epidermis, which was more than twice the number in mesophyll (Figure 3A). Only four miRNAs were found in both tissues (Figure 3B); 28 miRNAs were discovered under the D25 condition in epidermis (Figure 3B). These results suggested that the drought stress might cause substantial changes of small RNAome, which could subsequently regulate the gene expression. We then focused on the Kalanchoë-specific miRNAs expression patterns under the drought conditions. We found that six Kalanchoë-specific miRNAs were responsive to the drought treatments (Figure 3C). Particularly, Kal-miR-C1878, significantly upregulated in both mesophyll and epidermis, and Kal-miR-C689, which showed opposite patterns between mesophyll and epidermis under the drought conditions (Figure 3C).

The diurnal expression of gene expression is fundamental to CAM, and we ask whether the differentially expressed miRNAs under drought stress were involved in the circadian expression of genes. Through the Parallel Analysis of RNA Ends (also known as PARE sequencing or degradome), we predicted the targets of miRNAs by using the transcripts from the Kalanchoë genome database. In total, 381 miRNA-target regulations (p < 0.05) were identified (Supplementary Table S3). To further search for the miRNA-target regulation related to CAM pathway, we evaluated the diurnal transcriptome dataset from K. fedtschenkoi leaves [19]. We focused on genes that displayed the expression peaks at morning and dusk periods (Supplementary Figure S2A; Supplementary Dataset S3). Interestingly, we found that PLASMA MEMBRANE PROTON ATPASE 2 (Kaladp0008s0304) was identified as a target of miR369, which was revealed as a key regulator of stomatal movement in K. fedtschenkoi [19]. In the dusk-peaked gene cluster, we showed that two candidate targets of Kal-miR_C1817, including a MYB transcription factor and POF like gene, displayed high expression levels in the daytime (Supplementary Figure S2B). In the morning-peaked gene cluster, we found that two TZP-like genes were the targets of miR530 (Supplementary Figure S2B). TZP genes, also known as BLUS3-like genes, were essential regulators of circadian clock and blue light signaling pathways [41,42]. We hypothesized that the miR530-TZP regulatory module might be involved in directing the diurnal gene expression patterns in Kalanchoë.

To investigate the miR530-TZP regulatory module, we identified and validated the genomic locus of miR530 through cloning and sequencing. The precursor of miR530 contains a typical stem-loop structure with abundant accumulation of small RNAs at the mature and star regions (Figure 4A). Based on the deep sequencing analysis, we found that the mature miR530 levels in mesophyll tissues were not significantly changed under the drought treatments (Figure 4B), while more than 20-fold reduction was revealed in the epidermis under D10 and D25 conditions (Figure 4C). To validate the targets of miR530, the Dual-Luciferase Transient Expression System was employed to evaluate the efficacy of miRNA and its targets [38]. This analysis showed that both targeting sites of TZP1 and TZP2 were regulated by miR530 (Figure 4C), suggesting those transcripts were both cleaved by miR530 in Kalanchoë.

In order to reveal the roles of miR530 in the regulation circadian gene expression, we performed time-course analysis of miR530, TZP1 and TZP2 expression. Under the normal light–dark cycling condition, we showed that the miR530 expression peaked in the middle of the night and dropped to the lowest level before dawn (Figure 5A); this pattern negatively correlated with expression of TZP1/2 which peaked before the onset of light (Figure 5A). The expression pattern of miR530 under normal light–dark cycle was tested by JTK_cycle, which displayed a rhythmic expression (p-value = 0.003, period = 12, amplitude = 0.2) [43]. A dramatic reduction of TZP1/2 after lighting for 10 min, suggested that expression of TZP1/2 was also regulated by the light signaling pathway (Figure 5A). We continued to monitor the expression profiles for two consecutive days under the constant light conditions. We found that miR530 retained the similar accumulation pattern in Day 1, but the peaked expression shifted about nine hours in Day 2 (Figure 5B,C). Alternatively, the expression of TZP1/2 remained a single peak expression with slight shifts in Day 1 but displayed disrupted expression patterns in Day 2 (Figure 5B,C). Taken together, these results suggest that the expression of miR530 is potentially under the regulation of circadian rhythm, and both light signaling and circadian clock pathways might be involved in the regulation of TZP1/2 expression.

All data supporting reported results can be found from the links related to the manuscript.