Two Piwis with Ago-like functions silence somatic genes at the chromatin level

ABSTRACT

Most sRNA biogenesis mechanisms involve either RNAse III cleavage or ping-pong amplification by different Piwi proteins harbouring slicer activity. Here, we follow the question why the mechanism of transgene-induced silencing in the ciliate Paramecium needs both Dicer activity and two Ptiwi proteins. This pathway involves primary siRNAs produced from non-translatable transgenes and secondary siRNAs from targeted endogenous loci. Our data does not indicate any signatures from ping-pong amplification but Dicer cleavage of long dsRNA. Ptiwi13 and 14 prefer different sub-cellular localizations and different preferences for primary and secondary siRNAs but do not load them mutually exclusive. Both Piwis enrich for antisense RNAs and show a general preference for uridine-rich sRNAs along the entire sRNA length. In addition, Ptiwi14-loaded siRNAs show a 5´-U signature. Our data indicates both Ptiwis and 2´-O-methylation contributing to strand selection of Dicer cleaved siRNAs. This unexpected function of the two distinct vegetative Piwis extends the increasing knowledge of the diversity of Piwi functions in diverse silencing pathways. We describe an unusual mode of action of Piwi proteins extending not only the great variety of Piwi-associated RNAi pathways but moreover raising the question whether this could have been the primordial one.

KEYWORDS:

• Argonaute
• piwi
• sRNA loading
• RNA interference
• siRNA
• dicer
• transgene-induced silencing
• secondary siRNAs

Introduction

RNA silencing is a term describing a broad variety of mechanisms that use short RNA molecules to regulate gene expression. These can either target already transcribed mRNAs post-transcriptionally (PTGS) or they can interfere in transcription via co-transcriptional targeting of nascent transcripts (CTGS), thus recruiting chromatin modifying complexes [1,2]. An important component of RNAi (RNAinterference) mechanisms is Argonaute proteins (Ago), which load small RNAs (sRNAs), thus creating functional complexes. Agos themselves can be phylogenetically dissected into two clades: Agos and Piwis (P-element-induced wimpy testes), the latter being discovered in Drosophila germline stem cells. Agos form the RISC (RNA induced silencing complex), with miRNAs and siRNAs both being ubiquitously expressed, whereas Piwis and piRNAs were believed to be expressed in germline cells only [3].

siRNAs are a distinct class of regulatory RNAs produced by the dsRNA (double stranded RNA)-specific ribonuclease Dicer: these have been shown across kingdoms to act either in PTGS and CTGS of protein-coding genes and structural elements such as centromeres through the life cycle [1]. In many systems, secondary siRNAs have been shown to be produced involving activity of RNA-dependent RNA polymerases (RDR).

In contrast to Dicer-cleaved siRNAs, piRNAs were mainly described to silence transposable elements during gametogenesis. However, increasing data on different piRNA mechanisms reveal an unexpected diversity of those. This diversity does not only concern piRNA targets but also their temporal/spatial occurrence and, most importantly, the piRNA biogenesis mechanisms. In mouse and Drosophila, similar mechanisms were described: single-stranded precursor 5^´-U RNAs are loaded into Piwi proteins, and this is followed by subsequent 3´-trimming of the piRNA end. This initiation is then followed by Dicer-independent amplification of piRNAs by the ping-pong mechanism involving the reciprocal cleavage of complementary ssRNA, thus generating an internal single nucleotide A preference (reviewed in [4]).

3´-Nucleotides of mature piRNAs usually carry a 2´-O-methylation, which is added after 3´-processing. Some mechanistic diversity becomes apparent as Drosophila and mouse piRNAs become 3´-trimmed after Piwi loading, but this was not reported in Caenorhabditis elegans where mature piRNAs are generated from precursors by 5´- and 3´-processing, and then the mature piRNA is loaded into Piwi [5]. In all systems, 3´-methylation of piRNAs occurs after Piwi loading.

Piwi-mediated piRNA biogenesis differs to siRNAs and miRNAs maturation by the action of Agos and Piwis [6]. Agos load duplexes of Dicer cuts and select for guide and passenger strand before screening for targets. In contrast, Piwis do not show strand selection from duplexes and show a tendency for 5´-U containing ssRNAs, the latter has been shown to be due to both biased piRNA biogenesis and selectivity of Piwis [7].

As mentioned, piRNA pathways differ extremely between species and show a lack of conservation of involved genes [8]. Studies in Drosophila demonstrated that piRNA pathway genes evolve rapidly, indicating an arms race between transposons and their cellular defence [9,10]. Also downstream mechanisms differ between species as C. elegans shows an absence of ping-pong amplification using RDR-dependent siRNAs for amplification of the initial piRNAs [11]. Moreover, a screen of non-model species revealed the absence of the piRNA system in nematode lineages other than the most prominent C. elegans: these organisms apparently use RDR-dependent siRNAs to account for transposon control [12].

Consequently, piRNA biogenesis pathways are highly diverse, and increasing evidence indicates that piRNAs are not restricted to the germline but are present in low abundance also in somatic tissues, e.g. piRNA-like sRNAs have been identified in various somatic tissues by their ping-pong signature, and increasing reports also show piRNAs regulating expression of endogenous genes in somatic cells too [13,14]. However, purely descriptive studies reporting somatic piRNA expression need to be handled with care because of miss-annotations in piRNA databases containing piRNA-sized fragments of longer RNAs showing high RNA levels in somatic tissues [15].

As a result, an ongoing discussion asks for the evolution of these multiple functions of piRNAs, and one possibility is the co-option of transposon-derived piRNAs to regulate genomic functions [16]. A piRNA analysis of several arthropod species revealed somatic piRNAs targeting transposons and mRNAs across all species, and the authors consequently scrutinize that the ancestral role of piRNAs was to protect the germline from transposons [17].

In the context of changing dogmas about Piwis, ciliates provide an excellent model. They use several vegetative and developmental sRNA pathways [18], and they do not harbour any Agos. Paramecium tetraurelia, for instance, contains 15 distinct Piwi proteins called Ptiwis [19]. These unicellular eukaryotes undergo sexual recombination of meiotic nuclei in order to develop somatic macronuclei (MAC) in the same cell with germline micronuclei (MIC). For further evaluation of Paramecium Piwi functions in the evolutionary context, one should be aware that the ciliate nuclear dimorphism cannot be seen as an ancestral state of multicellular species as germ-soma nuclear differentiation evolved at least twice in unicellular species [20]. As such, vegetative cells comprise functions of germline and somatic cells, rather than being exclusively somatic or germline.

Here, we characterize two Ptiwi proteins (Ptiwi13 and Ptiwi14) expressed during vegetative growth of Paramecium tetraurelia. Both are involved in transgene-induced silencing. In this mechanism, transformation of non-expressible transgenes silence endogenous gene loci and this has been shown to involve dynamic chromatin remodelling [19,21,22]. This is at first glance similar to other mechanisms in which small RNA-mediated interaction of two different genetic loci have been reported. Next to co-suppression in plants where endogenous and exogenous homologous genes are post-transcriptionally silenced by siRNAs [23], also the paramutation represents an epigenetic interaction between two different genetic loci [24]. Interestingly, paramutations in Drosophila and C. elegans can involve both piRNA and siRNA elements, respectively [25,26].

The mechanism of transgene-induced silencing in Paramecium differs from those phenomena as only truncated, non-expressible transgenes can trigger silencing [27]. Transcription of translatable and intact mRNA from transgenes appears to repress silencing and, in addition, also the deletion of these genes in F1 progeny [28]. Such transgenerational manifestations can also be observed in C. elegans and Drosophila; however, inheritance there concerns only gene silencing, not gene deletion [26,29].

In Paramecium, the precise characterization of sRNAs of transgene-induced silencing and especially of their bio-accumulation is missing. Two distinct Ptiwi proteins are involved, which is in conflict with the fact that Dicer is involved in the mechanism [22,30]. The sRNA specificity of the individual Ptiwi proteins remains unknown, and as such, their role and the origin and function of their associated sRNAs also remains elusive. The aim of this study is to dissect transgene-induced sRNAs by their loading into Ptiwi proteins to clarify about the the role and the mechanism of the two distinct Ptiwis in the mechanism in Paramecium.

Materials and methods

Cell culture, RNAi, microinjection

Paramecium tetraurelia cells (stock 51 and d4-2) were cultured as described before using Klebsiella planticola for regular food in wheat grass powder (WGP) [31]. All cultures for this study were grown at 31^°C. RNAi by feeding of dsRNA-producing bacteria was carried out as described before [32,33] using the double T7 vector L4440 in the RNAse III-deficient E. coli HT115DE3. Microinjection of the pTI-/- and FLAG fusion transgenes was carried out as described before [34].

Phylogenetic analysis

The evolutionary history was inferred using the Neighbour-Joining method with 1000 bootstrap replicates [35,36]. The optimal tree is shown in Fig. 1A. Evolutionary distances were computed using the Poisson correction method [37] and are in the units of the number of amino acid substitutions per site. Ambiguous positions were removed by the pairwise deletion option. There were a total of 1.703 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [38]. Proteins were aligned with Muscle using default parameters. Ptiwi1-15 sequences were described in [19], and we also used a curated amino acid sequence for the putative pseudogene Ptiwi04 using its paralog Ptiwi05 as a template.

Figure 1. Relationship of Ptiwi proteins and their localization in vegetative cells. A: Phylogenetic tree of Paramecium tetraurelia Ptiwi proteins (blue) in relationship to Tetrahymena Piwis (Twi, red), metazoan Piwi proteins (green) and Agos (yellow). Support values are given at nodes, see Methods for details. The amino acid sequence of the putative pseudogene Ptiwi04 was corrected using its paralog Ptiwi05. B: Localization of Ptiwi proteins in vegetative Paramecium cells injected with Ptiwi13-FLAG (top) or Ptiwi14-FLAG (bottom). Cells were analysed by indirect immunofluorescent staining using custom antibodies directed against Ptiwi13 or Ptiwi14 labelled with secondary Alexa594-conjugated antibody (red). Additionally, the cells were stained using anti-FLAG antibody. Representative overlays of Z-stacks of magnified views are presented. Other panels show DAPI (in blue), brightfield and overlay of DAPI and Alexa594 signal. White asterisk indicates the position of the macronucleus. Scale bar is 10 $\mu$m and exposure is 2 s

RNA isolation and treatment

Total RNA was isolated with TriReagent (Sigma). Integrity was checked by denaturating gel electrophoresis after DNase I (Invitrogen) digestion and subsequent purification with acid phenol. For dissection of 3ʹ-modifications by periodate oxidation, 20 $\mu$g RNA were dissolved in 17.5 $\mu$l 4.375 mM borax, 50 mM boric acid, pH 8.6, and 2.5 $\mu$l 200 mM sodium periodate were added. After 10 min incubation in the dark, 2 $\mu$l glycerol were added with another 10 min incubation. After drying in the speedvac, the pellet was dissolved in 50 $\mu$l 33.75 mM Borax; 50 mM boric acid; pH 9.5 and incubated for 90 min at 45^°C . The RNA was subsequently purified with Sephadex G-25 columns (GE).

sRNA sequencing and analyses

For siRNA sequencing, 17–25 nt small RNA fractions were isolated by denaturating PAGE and subjected to standard small RNA library preparation using the NEB Next small RNA sequencing Kit (NEB, Frankfurt a.M., Germany). The procedure includes 3´-OH and 5^´-monophosphate-specific ligation steps, and we tried to lower 3^´-2^´-O-me biases by 18 hours 3^´-ligation at ${16}^{\circ }$C. After 10 PCR cycles, the libraries were gel-purified and sequenced on the HiSeq 2500 using the Rapid Mode with 28 cycles. Reads were de-multiplexed, and adapter sequences were trimmed using Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim galore/) that uses Cutadapt [39] with a stringency cut-off of 10. For analysis of reads, we used normalized counts and converted these values to transcripts per million (TPM), which we also refer to as sRNA accumulation. For the analyses specific to endogenous clusters, a TPM value greater than one was termed to be present in Ptiwi IPs. We used the RAPID pipeline to obtain the normalized counts, implementing the KnockDown Corrected Scaling (KDCS) method [40]. Hierarchical clustering of data sets was performed with complete linkage using an Euclidean distance measure, and heatmaps were created using R/Bioconductor package gplots (v3.0.1.1). Shown statistical analyses were performed on average of reads of dupli- or triplicates of experiments, including error bars based on calculated variance. Data are deposited at the European Nucleotide archive, ENA, Acc. Nr. PRJEB38766, and information of sequencing depth and mapping statistics are provided in Table S1.

sRNA signatures

Sequence logos of 23nt sRNAs were generated using WebLogo3 [41] with error bars twice the height of the correction for small sample size. Probabilities for overlapping reads from aligned sRNA reads were calculated using the small RNA signature analysis tool in Galaxy [42]. sRNAs from 17 to 25 nt were mapped to each region of interest, allowing no mismatches or multimapper in bowtie [43], and overlaps from 1 to 25 nucleotides were calculated. Plots for read length distribution and coverage were created using Geneious Prime 2020.1.2.

Antibodies, western blots, immunostaining

Peptides corresponding to the amino acids 684–698 and 449–463 of Paramecium Ptiwi13 (C-DDAPPQARKNNKSPY) and Ptiwi14 (C-QNWMQRLTAEIGDK), respectively, were used for immunization of rabbits. Purification of antibodies from serum was performed by coupling the respective peptides to SulfoLink coupling resin (Thermo Scientific) and following the manual instructions. Purified antibodies were tested by dotblot assays (Fig S2). Western blots were carried out as described previously [44] using indicated antibodies diluted 1:250 in 5% milk/TBST. Indirect immunofluorescence staining was carried out as previously described [45]. Cells were permeabilized in 2.5% Triton X-100 and 1% formaldehyde for 30 min followed by fixation in 4% formaldehyde and 1.2% Triton X-100 for 10 min. After blocking in 3%BSA/TBST, the cells were incubated in primary antibody diluted 1:200 in 3%BSA/TBST under mild agitation overnight at 4^°C. After washing and incubation with 1:2500 Alexa Fluor 568 F(ab’)2 fragment of goat anti-Rabbit IgG (H + L) (Thermo Scientific # A-21069), the cells were stained with DAPI and mounted in VECTASHIELD (VectorLaboratories). Images were acquired using Zeiss Axio Observer with ApoTome. For expression of tagged Ptiwis, the respective orf was cloned into Paramecium FLAG-Vectors pPXV containing three FLAG sequences either at the N or C terminus (kind gift of M. Valentine and J. Van Houten, Vermont, USA) as described in [46]. Injected clones were screened by single-cell PCR, and positives were grown for cell fixation and protein isolation. Macronuclei were isolated as described [47], and protein was isolated by adding preheated Laemmli sample buffer with subsequent boiling for 5 min.

Ptiwi-immunoprecipitation

Transgenic Ptiwi lines harbouring the pTI-/- transgene and a single Ptiwi-FLAG fusion construct (described above) were used for Ptiwi IPs using monoclonal anti-FLAG M2 (Sigma). Our procedure follows the protocol by [48] for developmental Ptiwis with the following modifications. 500 k cells of a single transgenic line were grown and harvested by snap freezing in 2 ml lysis buffer. 1 ml of the lysate was broken in a Dounce homogenizer and 1 ml was sonified until also MACs were destroyed. After addition of 50 $\mu$l Anti-FLAG M2 Magnetic Beads (Sigma) and incubated over night by gentle agitation. After washing beads with wash buffer and re-suspended in 100$\mu$l. 10$\mu$l were used for western controls by addition of 2.5$\mu$l Laemmli sample buffer and subsequent boiling for 2 min. 90$\mu$l were extracted with TriReagent LS (Sigma) according to the manufacturer’s recommendation.

Results

Ptiwi13 and 14 prefer different subcellular localizations

Fig. 1A shows the evolutionary relationship between the 15 Paramecium tetraurelia Ptiwi proteins. The phylogenetic tree reveals ciliate Piwi proteins clustering with metazoan Piwis, which are clearly separated from Agos. The two Piwis involved in transgene-induced silencing, Ptiwis 13 and 14, show some degree of similarity but no close relationship. They do not appear to be the result of one most recent genome duplication of which the Paramecium tetraurelia genome has undergone at least three. Ptiwi 14 has an ohnolog of the most recent WGD, Ptiwi08, which shows developmental expression [49,50]. In addition, a recent study identified orthologs of Ptiwi 13 and 14 in most of the species of the Paramecium aurelia complex and in addition in Paramecium caudatum and Paramecium bursaria [51]. An analysis of the catalytic domain (Fig S1) shows that the catalytic DEDH tetrad [52] is present in Ptiwi13 and 14, suggesting that both are capable of slicer activity.

To clarify the subcellular localization, we raised antibodies against specific peptides corresponding to both Ptiwis for immunolocalization (Methods and Fig S2) and additionally used FLAG-tagged Ptiwi transgenes, which we used for later Ptiwi IPs as well. Figure 1B indicates clear cytosolic Ptiwi13 signals in stainings with specific and FLAG antibodies. This cytosolic signal appears a bit structured, likely due to fixation-induced binding of soluble proteins to ER membranes. Ptiwi13 custom antibodies reveal additional MAC signals in ca. 20% of cells as shown in Fig S2, and also the FLAG antibodies do not show an absence of Ptiwi13 signal in the MAC. We conclude that Ptiwi13 has a predominant localization in the cytosol but can also appear in the MAC. Ptiwi14 staining with custom- and FLAG-Abs shows mainly MAC signals and only faint signals in the cytoplasm. We conclude that both Ptiwis have different sub-cellular localization preferences in MAC and cytosol, but for both Ptiwis, also less intense signals in the respective other compartment are apparent. This is supported by the comparison of total and MAC protein Western blots shown in Fig S2 and by an in silico analysis of the amino acid sequences by ngLoc method, which is a Bayesian classification method to predict localization of proteins [53]. According to the multi-localization confidence score (MLCS above 20), the algorithm predicts that both Ptiwis shuttle between nucleus and cytosol, where the evidence for Ptiwi13 is higher (Fig S3). Slight differences in the ratio of Mac and cytosolic signals between FLAG and custom Abs may be due to increased levels of antigenic sites, resulting from the over-expression of FLAG-Ptiwi constructs. As these data exclude that the over-expression or the FLAG tag causes false positive localization, we proceeded with Ptiwi-IPs of these fusion proteins.

Ptiwi13 and 14 have different loading preferences for endogenous and exogenous sRNAs

For Ptiwi IPs, we injected FLAG-tagged Ptiwi13 and 14 transgenes, respectively, into a transgenic RNAi-strain harbouring the pTI-/- transgene (Fig. 2A). The latter contains a GFP marker and, additionally, a truncated version of the endogenous ND169 gene causing silencing of the endogenous locus. Western blots of aliquots of the lysates/pulldowns verified the successful immunoprecipitaion and the absence of soluble proteins present in the supernatant (Fig. 2B). As a first insight, the trimmed read length distribution of Ptiwi IP reads shown in Fig. 2C reveal a clear 23nt peak, which is the predominant siRNA read length in Paramecium. We then mapped reads to different classes of RNA templates and quantified them relative to the respective abundance in Ptiwi overexpression lines to limit the effects of an individual Ptiwi overexpression to stabilization of individual RNA species. Please note here that Ptiwi over-expression may cause unspecific binding of abundant RNA species.

Figure 2. Analysis of sRNAs in Ptiwi immunoprecipitations. A: Experiment overview. A single cell was injected with the pTI-/- transgene. After establishment of a stable line, the cells were injected with FLAG-Ptiwi13/14 constructs, respectively (green). B: Control Western blots for the IPs using anti-FLAG Abs for Ptiwi detection and anti-GFP (Sup.-Supernatant, IP-Immunoprecipitation). Two different setups of the IPs used sonication (S) and douncing (D) for cell lysis, the latter remains MAC structure but permeabilized. C: Total read length distribution of Ptiwi IPed reads after adapter trimming. Average of reads from three IP replicates is shown. D: Relative enrichment of RNA reads in Ptiwi IPs mapping to different categories of genomic templates. Average of reads from three IP replicates was calculated, and the enrichment in reference to individual Ptiwi overexpressing lines is shown. * p-value < 0.005

Fig. 2D indicates that Ptiwi13 enriches for sRNAs of exogenous precursors such as food bacteria and mitochondria. This is in agreement with a previous report that Ptiwi13 is also involved in exogenously triggered RNAi when paramecia are fed with dsRNA producing bacteria [19]: it has later been shown that Paramecium also converts exogenous ssRNA of the food bacteria such as rRNA and mRNA into siRNAs [19,54]. In addition, Ptiwi13 also enriches for fragments of rRNA and snRNA.

In contrast, Ptiwi14 IPs show accumulation of small RNAs produced from all protein coding genes (all MAC genes) and a subset of previously characterized siRNA producing genes (SRCs, small RNA clusters) of the Paramecium genome [55]. Small RNAs from ncRNAs are clearly underrepresented in Ptiwi14 IPs as well as fragments of snoRNA and tRNAs. The latter are also not found in Ptiwi13 IPs.

In summary, our data indicates both Ptiwis not to be redundant but with different localization and loading preferences.

2^° siRNAs are enriched but not exclusively found in Ptiwi14

Both Ptiwis have been earlier shown to be necessary for efficient transgene-induced silencing [19]. Figure 3A shows the genomic structure of the endogenous ND169 gene involved in trichocyst discharge. This gene becomes silenced on the chromatin level when cells are injected with a truncated form of this gene shown below: the pTI-/- transgene shows two deletions: one on the 5ʹ-coding region (ND-1) and the 3ʹ-coding region including the 3^´-UTR (ND-2) [22].

Figure 3. Ptiwi13 and 14 load transgene associated sRNAs. A: Detailed scheme of the endogenous ND169 locus (top) and the truncated transgene (bottom). Introns are numbered and brackets symbolize specific junctions. Shaded regions are not part of the transgene (ND-1 and ND-2). B: Coverage tracks of siRNAs mapped to the endogenous ND169 locus. siRNAs were separated by their direction (sense/antisense). Regions accounted for 1^° siRNAs (NDgene, yellow) and 2^° siRNAs (ND-1 and ND-2, grey) are highlighted. Coverage track in log scale is shown for one replicate each, while numbers on the right indicate untransformed sense and antisense coverage. C: Read length distribution of 1^° siRNAs and 2^° siRNAs from Ptiwi IPs and total RNA from pTI-/- injected cells. Data are shown as proportion of reads mapping to the NDgene locus and D: the ND-1 locus. E: Relative enrichment of RNA reads in Ptiwi IPs mapping to different regions of the transgene. Enrichment was calculated for the average of reads of three replicates in reference to individual Ptiwi overexpressing lines. * p-value < 0.005

Mapping sRNA reads to the endogenous ND169 is shown in (Fig. 3B): siRNAs mapping to the regions called ND-1 and ND-2 result from the endogenous ND169 gene only, as these regions are not present in the transgene. siRNAs mapping to ND-1 and ND-2 therefore represent 2^° siRNAs.

They are therefore 2^° siRNAs produced from the targeted gene after triggering by 1^° siRNAs produced from the transgene [22]. Those regions existing in the transgene and the endogenous gene (called NDgene in the following) consist of both 1^° and 2^° siRNA; however, as the abundance of 2^° siRNAs is more than 10-fold lower compared to 1^°, the NDgene regions show predominantly 1^° RNAs [22]. Figure 3B also shows the reads from Ptiwi IPs. These maps already rebut one of our first hypothesis on the question why two different Ptiwis may be involved in this mechanism: both Ptiwi IPs show reads mapping to the NDgene and to the ND-1/2 region. As such, they do not mutually exclusively load 1^° siRNAs and 2^° siRNAs. Analysing quality and quantity of these sRNAs, Fig. 3C and Fig. 3D show that both Ptiwi bound sRNAs are of predominant 23nt length. Ptiwi14 significantly enriches more 2^° siRNAs (Fig. 3E). The finding that silencing of the ND169 gene by the transgene was shown to occur on the chromatin level [22] may make sense in this context as 2^° siRNAs may then be produced from nascent transcripts in the nucleus and loaded by nuclear Ptiwi14.

Ptiwi13 and 14 specifically load 23nt antisense siRNAs

To characterize the nature of Ptiwi loaded sRNAs, we had first a look at the ratio of 23nt reads to other read lengths. Comparing total RNA from pTI-/- transgenic cells to IPs, Fig. 4A shows that both Ptiwis specifically load 23nt sRNAs; however, the ratio of 23nt to other lengths varies between different RNA species. For several RNAs, e.g. food bacteria, rRNA etc., one can see that Ptiwis specifically select 23nt sRNAs among many other RNAs. It seems likely that fragments of these RNAs are produced by different mechanisms creating several lengths of sRNA of which Ptiwis enrich for 23nt sRNAs.

Figure 4. Asymmetric modification and Ptiwi selectivity contributes to accumulation of antisense siRNAs. A: Number of 23nt reads mapping to each indicated genomic feature and transgene regions accounting for 1^° and 2^° siRNAs were related to the total number of reads of other sizes. Calculation is shown for the mean of Ptiwi IPs and RNA from pTI-/- injected cells as control (total RNA). B: Antisense ratio of reads from Ptiwi IPs calculated by merging three replicates each. C: Difference in the antisense ratio of reads to the antisense ratio of respective control was calculated. Data is shown for reads from knockdown of Ptiwis in duplicates and IPs (triplicates) mapping to the indicated transgene regions. D: Antisense ratio of small RNAs from pTI-/- transgene samples (total RNA, untreated) and small RNAs treated with sodium periodate (+NaIO₄). Average of reads from two replicates is shown

This appears different for transgene-associated 1^° and 2^° siRNAs, which show almost identical ratio of 23nt siRNAs in total RNA and IPs, suggesting that distinct biogenesis mechanisms contributes to more precise sRNA cleavage. Going more into detail with these transgene associated siRNAs, both Ptiwis load predominantly antisense RNAs as shown in Fig. 4B. We compared this to Ptiwi knockdowns, in which the individual Ptiwis are silenced by introduction of dsRNA by feeding bacteria. In the transgene (pTI-/-) background, silencing of Ptiwi13 or 14 causes rescue of the ND169 silencing phenotype, meaning that cells can eject trichocysts again [22]. Comparing the antisense ratio of Ptiwi knockdown and IPs to each other, our data indicate a certain decrease of the antisense ratio in knockdowns and an increase in IPs (Fig. 4C). Its worth to note here that Ptiwi13 is recursive: being involved in the dsRNA feeding pathway may be less efficient compared to Ptiwi14 silencing but has been shown to be efficient for reporter RNAi rescue in several instances [19,22]. These changes in the antisense ratio are only moderate, and it is either possible that both Ptiwis are redundant or that also other factors contribute to strand selection. In many systems, 2ʹ-O-methylation was shown to occur in context of Piwi-associated sRNAs [56]. Using periodate oxidation of RNA and subsequent library preparation, we can show that both 1^° and 2^° are resistant to periodate thus likely to be methylated at the 3ʹ-end (Fig S4). Moreover, our data indicate that predominantly antisense sRNAs are modified in this manner (Fig. 4D), suggesting that this modification contributes to strand selection and stabilization.

Ptiwi14 loaded siRNAs have a 5ʹ-uridine preference

As the current data implicates that both Ptiwis select strands from Dicer cuts rather than amplify sRNAs in a ping-pong manner, we followed this idea by analysing the sequence logos of transgene-associated sRNAs. Figure 5A shows logos of 1^° and 2^° sRNAs obtaining 5ʹ-uridine preference. To decide whether these RNAs should result from a Dicer cut, one should see an A preference at position 21 for the non 5ʹ-U reads. These are shown in Fig. 5A, but one cannot identify such a Dicer signature nor a ping-pong signature (an A at position 10 of the non 5ʹ-U reads). Ptiwi IPs (Fig. 5B) reveal that the 5ʹ-U preference of total RNA is mainly due to Ptiwi14, whose RNAs show a much stronger 5ʹ-U preference compared to Ptiwi13. Unfortunately, lack of Dicer or ping-pong logos do not allow for further conclusions about the biogenesis mechanisms.

Figure 5. Sequence logos of 1^° and 2^° siRNAs. A: Sequence logos of 23nt antisense reads from pTI-/- injected cell lines mapping to transgene regions. Logos for either all sequences or the ones without 5ʹ-U are shown. B: Sequence logos of 23nt antisense reads from Ptiwi13 IP (left) and Ptiwi14 IP (right). C: Overlap predictions of small RNAs from 17 to 25nt of untreated, total RNA and the same small RNAs treated with sodium periodate (+NaIO₄). Z scores for overlapping pairs are included

We therefore additionally analysed reads for their overlapping signature: transgene 1^° siRNAs show a peak at 21nt overlaps, which fit to 23nt Dicer cuts (Fig. 5C). The 21nt overlaps are prominent but not dominant in Fig. 5C, which is likely due to the strand selection by Ptiwis, which causes degradation of the passenger strand. Fig. 5C in addition shows that we cannot identify any 21nt read overlaps in periodate-treated samples in agreement with the hypothesis of strand-specific methylation.

Interpreting this as an argument for Dicer cleavage, this is contrary to the missing Dicer signature in sequence logos, which would have been an an A-preference at position 21 in non-5´-U reads. We have to consider that the observed 5ʹ-U preference is not that strong compared to 5ʹ-Us in Paramecium scnRNAs for instance [30,57], and thus the complementary 21-As on the passenger strand might not be detectable. In addition, Fig S5 shows that the 5ʹ-U preference is still pronounced in periodate-treated RNAs, thus indicating that nucleotide preference and methylation co-occur on the same molecules.

sRNA uridine content contributes to strand selection

When analysing the sequence logos, not only the 5ʹ-U preference was observed but some logos suggested that stabilized strands are rich in uridines (Fig S5, Fig S6). We followed this by analysis of the antisense ratio along the transgene and endogenous ND169. Figure 6A shows again that most areas show dominant antisense preference for 1^° and 2^° siRNAs: an exception is the promoter proximal region (called ND-5), which shows almost 50/50 strand distribution. We consequently calculated uridine and adenosine content of these regions (Fig. 6B), revealing that the ND-5 region is different from the other regions as it shows a much higher uridine content on the sense strand. We therefore asked whether this could be seen in sRNAs, too. Also in sRNAs, the ND-5 region shows a different behaviour compared to other regions (Fig. 6C), and we consequently calculated the U-content of 23nt sRNAs of (i) in silico diced RNA, (ii) total transgene siRNAs and (iii) siRNAs of Ptiwi IPs. As demonstrated in Fig. 6D, the exceptional sense bias of the promoter proximal ND-5 region correlates to the enrichment of U-rich sRNAs, mainly by Ptiwi14. The analysis further reveals first that for all regions, the U-content of the more abundant antisense RNAs is higher than the sense siRNAs, and second for almost all regions for 1^° and 2^°, the U-content of Ptiwi IPed sRNAs is above in silico diced sRNAs. Thus, strand selection by Ptiwi13 and Ptiwi14 seems to include a general selection for uridine-rich sRNAs.

Figure 6. U-content analysis of 23nt sRNAs in Ptiwi IPs. A: Sense/antisense reads from pTI-/- injected cells mapping to indicated regions in 50nt windows. B: Percentage of adenine and uridine of the sense RNA transcript of each region. C: Heatmap of reads mapping to the indicated regions separated by their direction. D: X axis shows U-content of reads found in Ptiwi IPs, while density represents number of reads. In silico data is generated by counting U-content of 23mers of the DNA sequence for each region

Transgene-induced silencing mimics endogenous siRNA accumulation

We finally had a closer look at endogenous siRNAs. We have recently described siRNA-producing loci in the Paramecium genome, showing read length preference of 23nt [55]. Figure 7A shows also for these endogenous clusters a predominant overlap of 21nt, indicating Dicer to be involved at least in the majority of them. Surprisingly, most of the endogenous clusters can be identified also in IPs of Ptiwi13 and 14 (Fig. 7B, C and Fig S7). As shown in Fig S8, those SRC-derived, Ptiwi-loaded siRNAs with high antisense ratios also show a higher U-content in antisense reads. This again leads to the hypothesis of loading preferences of Ptiwis for uridine-rich siRNAs. Supplement Fig. S9 shows that both Ptiwis load SRC small RNAs independent of the expression level of genes: as Paramecium SRCs correlate with both, silent and high expressed genes, the function of these sRNAs is hardly understood, but the Ptiwi-IP data here suggest that Ptiwis do not dissect between sRNAs from silent or high expressed genes. In these analyses, we will miss any trans-acting mechanisms as mapping with only one mismatch allowed will likely result in cis acting correlations only. It is likely that also trans actions of SRC-produced siRNAs could occur, which is difficult to analyse as we do not know about the target recognition of Ptiwi-bound sRNAs. Correlating gene expression level with sRNA antisense ratio of SRCs located in protein coding genes, we also cannot identify a correlation between mRNA and antisense ratio of small RNAs in SRCs (Suppl. Fig S10). Few silent genes, however, show a strict bias of siRNAs, which could be examples of classical sRNA silenced genes but as this analysis was only possible for few genes (Suppl. Fig. S10). As a result of the comparison of transgene-associated silencing and endogenous small RNA producing loci, both appear to have common genetic requirements with endogenous siRNA accumulation pathways and appears therefore as a suitable model to study endogenous siRNA accumulation. Transgene-induced silencing appears therefore as a suitable model to study endogenous siRNA accumulation.

Figure 7. Endogenous sRNAs in Ptiwi IPs. A: Calculation of overlaps of endogenous 17–25nt sRNAs isolated from wildtype RNA including Z-score for overlapping pairs. B: sRNAs mapping to SRCs (small RNA clusters) in the Paramecium genome are analysed by their presence in Ptiwi IPs. Venn diagrams of the two IPs of sonified (S) and dounced (D) lysates with numbers of SRCs detected and the proportion of each fraction of the total found SRCs. C: Three examples of endogenous loci shown by coverage tracks of unnormalized data (surface antigen 51A, cluster22 and SRC 204)

Discussion

The dissection between Agos and Ptiwis was originally not only based on sequence similarity but on their spatial and temporal activity in germline and somatic cells. The most important distinction between both was due to their action in strand selection: Agos load RNAse III generated duplexes, whereas Piwis load longer ssRNA and generate their own sRNAs. Although many recent reports were published describing non-canonical functions of Piwis in somatic cells targeting non-transposable elements, the latter aspect of strand selection activity remains an important difference between the two groups.

Ciliates belong to the SAR-supergroup, which are equidistant to animals, plants and fungi [58] and especially Paramecium is different to many other species, not only metazoans, as the genome does not contain Agos but 15 Ptiwis [19]. Thus, Paramecium offers exciting possibilities for evolutionary comparison of RNAi mechanisms and the mode of action of individual components.

Transgene-associated siRNAs are Dicer products

We started this work here based on the surprising finding that two different Ptiwis are involved in a process where a non-expressible transgene silences a genetic locus at the chromatin level [19,22]. There are a couple of possible hypotheses why two distinct Ptiwis are necessary for this mechanism. First, a logical idea would be the action of ping-pong amplification. Although this would have been supported because Ptiwi13 and Ptiwi14 likely have slicer activity, it is in conflict with the fact that Dicer1 is necessary to produce at least the 1^° siRNAs (Fig S11). A second hypothesis would have been that both Ptiwis distinguish between 1^° and 2^° siRNAs. Our data clearly shows that both hypotheses are not true because we do not see any ping-ping signature but a 21nt overlap of reads, thus strongly suggesting that Dicer cuts at least the majority of siRNAs. Ptiwi IPs disprove also the second hypothesis because both Ptiwis load both 1^° and 2^° siRNAs, however, in different quantities. But what is the function of both Ptiwis then? Could they have the very same function being redundant? This seems not very likely because our data shows some discrete differences between Ptiwi13 and Ptiwi14: (i) Ptiwi13 loads more sRNAs from exogenous templates, (ii) both have different sub-cellular localization preferences, (iii) Ptiwi14 loads more 2^° siRNAs and finally (iv) Ptiwi14 shows a much stronger preference for 5ʹ-U RNAs. It seems therefore more likely that both Ptiwis have indeed distinct and specialized functions in this mechanism. Concerning the subcellular localization, our data resulting from IF, Western and in silico analysis support the idea that both Ptiwis can be found in both cytosol and MAC but with different preferences. Our data indicate to be more Ptiwi13 in the cytsol and more Ptiwi14 in the MAC. It seems likely that both are involved in a shuttling process between cytosol and nucleus, reminiscent of the nematode Ago NRDE-3, which localizes in the cytosol and redistributes to the nucleus when bound to 2^° siRNAs from the feeding pathway [59]. Also Arabidopsis Ago4 assembles with siRNAs in the cytosol and is then transported in the nucleus [60]. Among these examples for nuclear import of sRNA loaded Ago/Piwis, this would make sense if we think about the trigger for transgene-induced silencing. It has been shown that explicitly non-expressible transgenes induce RNAi [28,33], tempting that a quality control mechanism is involved in this process to dissect which transgene can produce translatable mRNA. Such a processes, e.g. nonsense-mediated RNA decay, works in the cytosol as translation is an efficient way to dissect between defect and translatable mRNAs. As we have previously shown that transgene-induced silencing works on the chromatin level, this cytosolic signal needs to be transported into the nucleus.

Ptiwis select strands from dicer products

Our data indicate that both Ptiwis select for strand-specific sRNAs from Dicer cut duplexes, which represents a non-canonical function of Piwi proteins. This finding is fostered by several aspects. Dicer knockdown reduces all sRNAs [22,30], and in addition we have shown here that bulk transgene siRNAs show 21nt overlaps of 23nt RNAs. Our study also brings strand-asymmetry in association with individual properties of sRNAs, apparently contributing to strand selection and stabilization by Ptiwis: 5ʹ-U preferences, U-content and 3ʹ-methylation. 5ʹ-U preferences have been frequently described, e.g. for the Piwi lacking Arabidopsis Ago1 [61]. It is also quite reminiscent to the strong 5ʹ-U preference of Drosophila Piwi (and weaker in Aubergine), which act together with Ago3 in the ping-pong amplification of piRNAs [62]. However, we cannot identify any sRNAs with an A-preference at position 10, which would result from such a mechanism. 5ʹ-Nucleotide preferences were also reported for the developmental Ptiwis in Paramecium showing a strong 5ʹ-UNG signature [30,48,57]. Recent evidence from in vitro dicing experiments show that this signature is due to cleavage preference of the involved Dicer-like enzymes rather than due to preferential Ptiwi loading [63]. The authors speculate about a co-evolution of Dicer-like enzymes to produce sRNAs targeting germline-specific DNA with strong sequence bias at its ends. This seems likely for the particular need to target conserved sequence-ends, but for the control of endogenous gene expression, accumulation of such conserved sequence features in siRNAs would not make sense. Interestingly, the scnRNA mechanism holds another interesting aspect to discuss as it involves also a special kind of 2^° sRNAs (iesRNAs) transcribed from already excised and circularized IESs. However, 1^° and 2^° developmental sRNAs have been shown to be loaded in distinct Ptiwis. This is contrary to our vegetative mechanism; this makes sense because scnRNAs and iesRNAs have different biogenesis mechanisms and different properties [57,63]. Vice versa, we may conclude from this that 1^° and 2^° siRNAs in transgene-induced silencing have identical biogenesis mechanisms, meaning that the same mechanism act on the transgene and the endogenous gene. Similar to transgene-associated Ptiwis, the knockdown of Ptiwi10 and Ptiwi09 during meiosis resulted in accumulation of duplexes and the authors concluded that these Ptiwis could be responsible for strand selection in reminiscence of Agos [48]. This is further supported by the involvement of Dicer/Dicerl-like proteins in the biogenesis of these two sRNA classes [57], thus indicating that the two developmental Piwi proteins act also more like Agos, which is similar to our finding here. Ciliates may in general use Piwis in an Ago manner, which is surprising not only for the vegetative Ptiwis described here but even more for the massive elimination of transposons and transposon-derived sequences during development, which is the hall mark of Dicer-independent piRNA in other organisms.

Nucleotides as a biochemical reason and the resulting thermodynamic behaviour of a sRNA duplex are only individual aspects for strand selection. It has been demonstrated that many protein factors in addition to Agos contribute to strand selection, allowing for dynamic adaptation of the miRNA system in response to challenges to adapt gene expression [64]. Two aspects need to be taken into account: phosphorylation by de novo RDR initiation and availability of RISC targets. Concerning the latter, it has first been shown in plants that the alteration of the target RNA binding quality alters miRNA abundance [65]. It seems likely that such a parameter also contributes to the strand selection in our example, as for instance the ND-5 region still does not show a clear sense bias, which could be explained if indeed antisense strands are preferred due to available targets by a sense (m)RNA. As mentioned above, also phosphorylation of RDR transcripts may play a role in strand selection: Tetrahymena Dicer2 was shown to be physically coupled with the RDRC: in vitro, Dicer2 cleaves discrete siRNAs from the 5´-triphosphorylated ends of dsRNA only [66]. As also the transgene-induced silencing here employs RDR activity, cyclic and phased de novo RDR activity on the sense transcript could produce duplexes with triphosphorylated antisense ends. If Ptiwis would select for those, the RDR products would be preferentially loaded.

Conclusion

Our study is another evidence for the extreme diversification of small RNA amplification not only in ciliates and the data raises the question whether ciliates use their Piwis totally different to other species or whether this function of Piwis could be the primordial one and later split into Agos and Piwi mode of actions. Vice versa to the ciliate Piwis in strand selection, yeast Ago has been demonstrated to load single-stranded RNAs, which become trimmed into pri-RNAs which are Dicer independent [67]. This means that S. pombe Ago can also process longer ssRNA into functional small RNAs. Such a Piwi-like function may fit to the position of S. pombe Ago between mammalian Agos and Piwis in our phylogenetic reconstruction in Fig. 1A. This non-canonical Ago function in yeast together with our data in ciliate Piwi function, let us assume that Piwis and Agos in unicellular eukaryotes are more diverse in their activity than expected. One may hypothesize that yeast and ciliates owning only Agos or Piwis, respectively, use those in a more flexible way compared to species harbouring both Agos and Piwis. This flexibility along with the absence of miRNAs in ciliates [55] may also be compared with the highly conserved miRNA loading Agos in vertebrates, which show a much higher degree of conservation compared to siRNA loading Agos of e.g. nematodes and insects, which are still in an arms race with viral adaption [68]. Although this comparison of conserved miRNA loading Agos in mammals and evolutionary flexible antitransposon and antiviral Agos/Piwis in plants, nematodes and single-celled organisms make sense at first glance, recent studies also demonstrate antiviral RNAi in interferon deficient mammalian cells [69], so depending on the extent of this, also the mammalian RNAi mechanisms still needs to adopt to new pathogens.

As until now, no ping-pong amplification has been demonstrated in any ciliate, further research has to clarify whether this is absent in ciliates at all. From the evolutionary point of view, Ago-like usage of Piwis seems surprising in unicellular eukaryotes as Agos have been demonstrated on bacteria already [70], although their function is less understood. Further studies need to clarify whether Agos may have been depleted in ciliates or, on the other hand, if strand selection of dsRNA by Piwis may have been the primordial function.

Acknowledgments

This work was supported by grants from the German research Council (DFG) to MS (SI1379/3-1), MHS SCHU (3140/1-1) and MJu MJ (CRC894). We are grateful to Dominique Furrer and Mariusz Nowacki for sharing RIP experience, Kaz Mochizuki for advice in Ptiwi IPs, and Megan Valentine and Judy VanHouten for sharing the FLAG fusion vector. Finally thanks to Christoph Kellner for scripts on U-counts and Anke Behnke for help in phylogenetic analyses.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by grants from the German research Council (DFG) to MS (SI1379/3-1), MHS SCHU (3140/1-1) and MJu MJ (CRC894).