Exogenous RNAi mechanisms contribute to transcriptome adaptation by phased siRNA clusters in Paramecium

Abstract

Extensive research has characterized distinct exogenous RNAi pathways interfering in gene expression during vegetative growth of the unicellular model ciliate Paramecium. However, role of RNAi in endogenous transcriptome regulation, and environmental adaptation is unknown. Here, we describe the first genome-wide profiling of endogenous sRNAs in context of different transcriptomic states (serotypes). We developed a pipeline to identify, and characterize 2602 siRNA producing clusters (SRCs). Our data show no evidence that SRCs produce miRNAs, and in contrast to other species, no preference for strand specificity of siRNAs. Interestingly, most SRCs overlap coding genes and a separate group show siRNA phasing along the entire open reading frame, suggesting that the mRNA transcript serves as a source for siRNAs. Integrative analysis of siRNA abundance and gene expression levels revealed surprisingly that mRNA and siRNA show negative as well as positive associations. Two RNA-dependent RNA Polymerase mutants, RDR1 and RDR2, show a drastic loss of siRNAs especially in phased SRCs accompanied with increased mRNA levels. Importantly, most SRCs depend on both RDRs, reminiscent to primary siRNAs in the RNAi against exogenous RNA, indicating mechanistic overlaps between exogenous and endogenous RNAi contributing to flexible transcriptome adaptation.

Figure 1.

(A) Overview of the small RNA cluster (SRC) generation workflow. The first row visualises the different serotypes according to a transcriptome analysis done in (27). (B) Length distribution of SRCs. (C) Number of serotype specific SRCs (y-axis) detected in the WT serotype samples (replicates were merged), stratified according to the predominant small RNA length (dicer call), where N means that no predominant length could be found.

Figure 2.

(A) Serotype specific SRCs expressed in the WT serotype (replicates were merged) samples were overlapped with annotated regions. Each annotated element is counted only once (distinct counting) and the number of elements of the different types is shown on the y-axis for all four serotypes. (B) sRNA reads (log10) in SRCs overlapping different genomic annotations (rows) and restricted to small RNA length (x-axis) for 51A serotype. (C) Length distribution of sense and antisense sRNAs mapping to protein-coding genes in 51A serotype. (D) Set intersection plots for SRCs overlapping with genes across serotypes. Genes consistently overlapping in the four WT serotypes, are called Genes associated with SRCs (GSRCs). (E) Boxplot of mRNA expression (y-axis, log2 TPM) in the four WT serotypes of GSRCs and other expressed genes. (F) Boxplot of total sRNA reads (y-axis; log2) in the EEJ and introns of GSRCs. (G) A histogram of the number of GSRCs (y-axis) plotted against the overlap percentage (in basepairs) of a GSRC by one or more SRCs (x-axis) is shown here. Based on the overlap percentage, GSRCs were split into three classes as shown (vertical lines) for further analysis. *P < 0.05 with a two-tailed Wilcoxon test.

Figure 3.

(A) Total number of phased clusters observed in all WT serotypes. (B) Example IGV screenshot of identified phased cluster, C909, annotated as a gene (ID: PTET.51.1.G0170152). (C) Set intersection plots for phased SRCs identified in all WT serotypes. (D) Number of phased clusters in each class of GSRC (as explained in Figure 2).

Figure 4.

(A) Box plots of the total small RNA read counts (y-axis; log2) from the WT, and mutant samples (Mutants: 51A-Rdr1, and 51A-Rdr2). The P-values indicated are from two-tailed wilcoxon test. (B) Set intersection plots for Differentially Expressed (D.E.) SRCs identified by DESeq2 in the mutant samples. We name the intersection groups as shown here (D - Downregulated, U - Upregulated, 1, 2 represents the mutant 51A-Rdr1 and 51A-Rdr2, respectively). For these intersection groups, we show the overlap percentage of GSRC by SRCs (C), the percentage of phased SRCs (D) and the sense/antisense ratio (y-axis; log2) in WT (E) here. In (E) the P-values (one-tailed wilcoxon test) are indicated for only groups with statistically significant differences.

Figure 5.

(A) Boxplots of total sRNA fold change (y-axis; log2 mutant/WT) in each mutant is shown. Boxplots are grouped based on whether a SRC is predicted as an unphased, or phased loci. (B) Same as (A), but shows the mRNA fold change (y-axis; log2 mutant/WT). The P-values indicated (A and B) are from two-tailed wilcoxon test. (C) A heatmap of all the phased SRCs, total small RNA read counts (log2) of all the WT replicates is shown.

Figure 6.

(A) A histogram of the number of GSRCs (y-axis) plotted against the Pearson correlation of mRNA expression and total sRNA (x-axis) of all WT replicates. Scatter plots of an example GSRC with high positive (B) and negative (C) correlation with mRNA expression (see r values in the plot). (x-axis; log2 TPM) against total sRNA accumulation (y-axis; log2 TPM) are shown.