Identification and analysis of functional associations among natural eukaryotic genome editing components

Paramecium strains, cultivation, and autogamy

All experiments were carried out with Paramecium tetraurelia stock 51, which may be obtained from the ATCC as ATCC30567. Cells were grown in a wheat grass powder (WGP; Pines International, Lawrence, KS) infusion medium bacterized the previous day with Klebsiella pneumoniae supplemented with 0.8 mg/l β-sitosterol. Cultivation and autogamy were carried out at 27°C as previously described²⁸.

Gene silencing

Gene silencing (knockdown) was carried out by RNA interference by feeding of Paramecium cells with E. coli expressing doubled-stranded RNA from a portion of the target gene cloned into the L4440 vector, as previously described^3,8,10,31. Sequence regions used for each of the genes silenced in DCL2/3/5-KD were as previously described^3,13. Silencing of PTCAF1 and ND7 were as previously described^10,31.

Three knockdowns of NOWA1/2 were performed; where necessary we distinguish between these knockdowns (i.e. NOWA1/2_a-KD, NOWA1/2_b-KD, NOWA1/2_c-KD). Macronuclear DNA from postautogamous cells from the first and third knockdowns was used to produce Illumina paired-end sequence libraries. Illumina sRNA libraries were produced from NOWA1/2_b-KD at the same time as those produced for sRNAs from DCL2/3-KD, DCL5-KD and control cells - the latter knockdowns were described in "Experiment 1" of 3.

The following controls were used in this manuscript:

Control 1 = feeding with an empty vector (EV) construct. This was the control for DCL2/3/5-KD and PDSG2-KD;

Control 2 = ND7-KD. This was the control for PGM-KD¹,

Control 3 = ND7-KD. This was the control for an unpublished experiment (Singh et al., in preparation);

Control 4 = ND7-KD. This was the control for EZL1-KD¹¹.

ND7 is a trichocyst discharge gene which is not expected to affect IES retention³².

Survival tests and cytological staging

Tests of survival of sexual progeny cells were performed by isolating and refeeding 30 individual cells. NOWA1/2 knockdown is generally lethal, i.e., typically all cells die within at most the time taken for 3 divisions of control cells following autogamy; for DCL2/3/5-KD 90% of cells were dead 3 days after autogamy; for PDSG2-KD after 3 days 50% of cells were dead and 47% appeared to be “sick” (i.e. dividing slower than usual); for PTCAF1-KD after 3 days 80% of cells were dead and 17% were sick.

Morphological characteristics of MAC development were used to assign the developmental stages of the control and knockdown samples: "Early", "Middle" and "Late", correspond to 50% of cells with fragmented MACs; 100% of cells with fragmented MACs; and 100% cells with fragmented MACs and a new MAC; Figure S3A,C). Note that because Paramecium cells are not completely synchronous during induction of autogamy³³ and because the underlying assessments are somewhat subjective, due caution should be exercised when comparing developmentally staged analyses between gene knockdowns.

Northern blotting

To assess the effectiveness of DCL2/3/5-KD, northern blotting was performed as in 3. The DCL2 probe was from gene position 857-1810 (ParameciumDB: GSPATG00003051001); the DCL3 probe was from gene position 518-1354 (ParameciumDB: GSPATG00027456001); the DCL5 probe was from gene position 2074-2461 (ParameciumDB: GSPATG00012932001). For PTCAF1-KD and control (ND7-KD) cells, northern blotting was performed as in^10,31.

Nowa1 antibody production and use in western blotting

A portion of Nowa1's C-terminal domain (amino acids 744-1024) fused to an N-terminal Strep-Tag was expressed in E. coli BL21-CodonPlus(DE3)-RIL cells (formerly Stratagene, currently Agilent, USA). The protein was column purified, followed by removal of the tag by HRV 3C protease cleavage. Immunization of a rabbit with this protein was performed by Eurogentec, Belgium, following their 28-day polyclonal antibody protocol, then affinity purified using the matrix bound C-terminal domain portion of Nowa1, before performing western blotting using a 1:500 antibody dilution (Figure S3B). As a loading control, western blotting was performed using 1:5000 dilution of monoclonal anti-alpha-tubulin (Sigma T8203 lot 111M4828) and 2% BSA for blocking.

Small RNA electrophoresis

For small RNA electrophoretic analysis, total RNA samples (5.0 μg) were 5'-end-labelled with [γ-³²P] ATP (5000 Ci/mmol, Amersham) by the exchange reaction of T4 polynucleotide kinase (Fermentas), then denatured and run on a 15% polyacrylamide, 7 M urea gel.

Electrophoretic mobility shift assay (EMSA)

Synthetic Nowa1 codon adjusted for E. coli translation (GENART synthesis) was cloned into the pMAL-c2x vector and expressed in Rosetta2(DE3)p-LysS (formerly Novagen, currently Merck) cells. Protein translation was performed at 37°C; cells were supplied with 1 mM IPTG and 0.2% glucose. Protein was purified according to the pMAL™ Protein Fusion and Purification System (NEB e8200s). EMSA was performed with the LightShift® Chemiluminescent RNA EMSA Kit (Thermo scientific, 201558) using 6.25 μM of RNA or DNA, 4 μg of MBP, and 12 μg of MBP-Nowa1.

DNA and RNA extraction, library preparation and sequencing

Macronuclear DNA was isolated from post-autogamous cells as previously described¹. For PCR analyses of IES retention, DNA was extracted from approximately 800 cells with a Nucleospin Tissue kit (Macherey-Nagel).

Total RNA was extracted from 400 ml of culture per the TRI Reagent BD protocol (Sigma) and was resuspended in RNase free water. A MirVana miRNA isolation kit (Ambion) was used for sRNA enrichment.

MAC genomic DNA and small RNA libraries were produced and sequenced per standard Illumina protocols. Paired-end Illumina DNA-seq libraries were constructed and sequenced on either the Illumina GAIIX, HiSeq 2000 or HiSeq 3000 sequencers.

Sequence data corresponding to the MAC DNA from NOWA1/2_a-KD can be obtained from the GenBank Sequence Read Archive (SRA) under the accession SRR646462. MAC DNA sequence data for NOWA1/2_c-KD, DCL2/3/5-KD, PTCAF1-KD, PDSG2-KD, ND7-KD and control 1 (empty vector) is available from the European Nucleotide Archive, respectively under the accession numbers ERS1033672, ERS1033670, ERS1033674, ERS1033673, ERS1033676, ERS1033671. sRNA sequence data for NOWA1/2_b-KD and empty vector control can be obtained from the European Nucleotide Archive accessions ERS1455004, ERS1455005, ERS1455006, ERS1455007, ERS1455008, ERS1455009, and for DCL2/3/5-KD, from ERS1455358, ERS1455359, ERS1455360.

PCR analysis of IES retention

PCR conditions and primers used to test IES retention following NOWA1/2_a-KD were the same as those in 3.

Reference genomes and IES retention scores

The following reference genomes for Paramecium tetraurelia stock 51, obtained from ParameciumDB³⁴, were used in the IES analyses and read mapping:

MAC: http://paramecium.cgm.cnrs-gif.fr/download/fasta/assemblies/ptetraurelia_mac_51.fa

MAC+IES: http://paramecium.cgm.cnrs-gif.fr/download/fasta/assemblies/ptetraurelia_mac_51_with_ies.fa (md5 checksums: dbcc54fb2987c8f60f8e765db7ed274c and 3e5b3fa65ebfaa484566a1ffddf20239, respectively).

IES retention scores (IRSs) were determined by ParTIES³⁵. For each IES, ParTIES counts mapped reads with unexcised IESs (IES+) and excised IESs (IES-) to determine an IES retention score, IRS = IES⁺ ÷ (IES⁺ + IES^-). A table of all the IES retention scores used in this manuscript is included as Dataset 1 (see Data and software availability). IES retention scores from previous knockdown experiments, i.e. PGM-KD¹, TFIIS4-KD¹², EZL1-KD (EZL1-2 silencing fragment)¹¹, DCL2/3-KD³ and DCL5-KD³, were downloaded from ParameciumDB³⁴. Unless otherwise indicated, NOWA1/2-KD IRSs are for NOWA1/2_a-KD.

IES and OES retention score regressions

To generate the IES and OES score regression matrices, shown in Figures 3A–C, Figure 4E and Figures S1A–B, we provide the program “after_ParTIES.py” (see Data and software availability), which can be run on Dataset 1 and Dataset 2 (IES and OES retention scores, respectively). This was tested using Python3.6 and relevant modules for scientific computing from the Anaconda software distribution from Continuum Analytics (version: Anaconda3-4.3.1).

Three regressions are graphically shown for IRS pairs: an ordinary least squares (OLS) regression, a LOWESS regression (non-parametric) and an orthogonal distance regression (ODR) with a linear function initialized by the slope and intercept from the OLS. OLS and ODR used standard SciPy functions. LOWESS regressions used the Python statsmodels library. ODR was performed since errors in both the "dependent" and "independent" regression variables are expected (both IRSs), not just the dependent variable as assumed by OLS. The combined effect of background IES retention and systematic IES retention errors on the regressions, and the effect of subtracting the mean of multiple control IES retention scores is shown in Figure S1A. For each pair of knockdown IES retention scores, Spearman’s rank correlation coefficient (r_s) was calculated instead of Pearson’s correlation coefficient to mitigate the effect of outliers. Reproducibility of IES retention scores between experiments is demonstrated by the strongly correlated scores of different NOWA1/2 silencing experiments (Figure S1B).

Two-sided p-values, calculated with SciPy’s spearmanr function, were used to test the null hypothesis of no correlation (for p ≤ α/m ≤ 0.01/45; i.e. with a Bonferroni correction for multiple hypothesis testing; uncorrected p-values are given in Supplementary Data S1–4), which was rejected in all the displayed correlations in the main figures, except where indicated by a circumflex accent. To facilitate assessment of the statistical significance of differences between r_s values, tables containing the two-tailed p-values are provided in Supplementary Data S5-7. p-values were calculated using methods from the R statistical package cocor version 1.1-3³⁶ that employ Fisher’s r-to-Z transformation (which are also applicable to r_s³⁷). Appropriate methods for calculating test statistics for correlations with independent groups or dependent groups with either overlapping or non-overlapping variables were used from cocor (noted in the supplementary data files). All the comparisons of r_s differences described in the text are statistically significant at α=0.01 with a Bonferroni correction for multiple hypothesis testing.

Assembly and selection of OESs

Genomic regions that are not part of the current reference Paramecium MAC+IES assembly were previously assembled¹, but we decided to use an alternative, metagenome assembler to attempt to produce a better assembly, and to minimize bacterial contamination. Reads from PGM-KD¹ not mapping to the current Paramecium tetraurelia strain 51 MAC genome and 51 MAC genome with IESs (bwa³⁸ version 0.7.12-r1039; default parameters) were assembled with IDBA-UD³⁹ version 1.0.9 (default parameters). To remove bacterial contaminants assembled during the assembly, only scaffolds with GC < 40% were selected (bacterial scaffolds peak around 58% GC; Figure S5A). This procedure resulted in a ∼20 Mb assembly consisting of 7266 scaffolds. As a small fraction of the P. tetraurelia MAC genome is currently missing from the reference assembly, this fraction was present in the ∼20 Mb assembly. To exclude the missing MAC genome regions, empty vector control (control 1) reads were mapped to the 7266 scaffolds with HISAT2⁴⁰ version 2.0.4 (parameters: “-X 600 --pen-noncansplice 0 --mp 24,8 --max-intronlen 1000), and then, using a custom Python script, only those regions with little or no sequence coverage (≤ 2.0×; compared to a mean MAC coverage of ∼100×) were extracted (Figure S5C shows an illustrative scaffold from which OESs were extracted). Coverage was smoothed using a bi-directional average of two exponentially weighted moving averages (EWMA) with a 200 bp span. This procedure extracted 4799 sequence regions (12.3 Mb).

To minimize the effects of potential genome assembly chaff (short assembled regions, esp. those that might arise from sequence errors), analyses were continued with only the low control sequence coverage regions extracted from scaffolds longer than 300 bp and with a PGM-KD sequence coverage of between 3× and 70× (Figure S5B). These 3882 regions, obtained from 5236 scaffolds, corresponding to 11.8 Mb of sequence data, are classified as putative OESs (Supplementary Data S8). Note that OES boundaries are a first approximation, since imprecise DNA elimination generates variable boundaries (as illustrated in Figure 1 and shown in Figure S5C), and many OESs situated at the ends of the scaffolds are likely incomplete due to the limitations of genome assembly (esp. since the assembly may end when untraversable repetitive regions are encountered).

Estimation of OES retention

For each knockdown, to quantify the effects of gene knockdown on OES retention, reads were mapped to the OESs (identified in the previous methods section) with HISAT2⁴⁰ version 2.0.4 (parameters: “-X 600 --pen-noncansplice 0 --mp 24,8 --max-intronlen 1000) before determining coverage, and normalized to the median coverage of the knockdown reads mapped with HISAT2 to the Paramecium tetraurelia strain 51 MAC genome plus IES assembly, i.e. OES retention score=OES_coverage/median(MAC-with-IES-genome_coverage). OES retention scores may not be as accurate as IES retention scores, both because their denominators use a measure of global, rather than local, sequence coverage and because only reads that map within OES boundaries are counted. In the latter case, read mapping may fail to detect short matches to an OES when only a short part of an OES/MDS-derived read matches, and so the shorter the OES, the greater the underestimation of its retention. Shorter OESs are more weakly retained in most knockdowns affecting genome reorganization (Figure S5D) in Paramecium, but at least some of this effect may be due to boundary-related underestimation of retention. In future, when a complete Paramecium MIC genome is assembled, the relationship between OES length and retention can be fully assessed. OES retention estimates are provided in Dataset 2. Control OES retention was subtracted from the knockdown OES retention for each OES used in the regression analyses.

sRNA mapping

sRNAs were mapped and histograms of their genomic distribution were produced as previously³ with the addition of OESs.

Sequence logos

Sequence logos were created with WebLogo⁴¹, version 3.3, with redundant sRNAs and base frequencies A=U=0.4 and C=G=0.1.

Examination of basic IES retention properties

Previously, each newly silenced Paramecium gene involved in IES excision appeared to generate a different distribution of IES retention^3,11,12. From additional knockdowns we have subsequently performed, it can now be seen that the IES retention distributions for NOWA1/2 knockdown, TFIIS4 knockdown and DCL2/3/5 triple knockdown are very similar (Figure 2A; see Methods for the IES retention score calculation). Given the role of TFIIS4 and Dcl2/3/5 in producing and processing of RNAs, and the suggested role of Nowa1 in RNA interactions (see “Nowa1/2 is involved in RNA scanning, may mediate scnRNA- and iesRNA-targeted DNA cutting and excision” later in the results), we classify these proteins as "RNA-associated" proteins involved in Paramecium MAC genome development. The IES retention score distributions of the EZL1 knockdown and PTCAF1 knockdown also have similar forms, though the effect of the former knockdown is slightly stronger (Figure 2B). Since both Ezl1 and PtCAF1 affect histone modifications^10,11, we classify these proteins as those involved in "histone modification".

Figure 2. IES retention scores and their relation to basic IES properties.

(A–D) For the sake of clarity, histogram bars with more than 10000 IESs have been truncated. The number of IESs binned in the left-most bars are as follows: (A) NOWA1/2-KD - 17468, TFIIS4-KD - 14152, DCL2/3/5-KD 18915 (Figure S3B and S3D provide Western and Northern blots, for NOWA1/2-KD and DCL2/3/5-KD, respectively); (B) PTCAF1-KD - 12614, EZL1-KD - 8047; (C) DCL2/3-KD - 37410, DCL5-KD - 35019, DCL2/3/5-KD - 22136; (D) Control - 43294; PDSG2-KD - 20730. Note that, as judged by survival tests in comparison to those previously published 8, PDSG2-KD may have been of reduced efficiency. (E,F) IES retention scores (scales on right y-axes) vs. IES length over shorter (E) (≤ 200 bp) and longer (F) (≤ 1000 bp) IES length scales. The IES length distribution is shown in grey in the background (scale on left y-axes). Lines are the exponentially-weighted mean averages (EWMA) with spans of 5 (E) or 50 (F) bp. (G–H) Mean base frequencies of positions 1-3 after the TA repeat relative are plotted relative to IES retention score for IESs from the third (45-55 bp) length peak (only the highest frequency bases are shown — see (H) for the base color scheme; similar trends are visible for ~10 bp windows surrounding longer IES length peaks). IES length is restricted to control for IES end base frequencies variation associated with large scale (26-1000 bp) IES length variation. EWMA lines for spans of 10 data points (intervals of 0.01) are plotted for data within 2 standard deviations (dotted vertical lines) of the mean retention score (dashed vertical line). Mean IES lengths (bp; counting only one of the two TA pointers) vs. IES retention score are potted in the bottom subgraphs (navy diamonds). See Figure S2 for end base frequency graphs of EZL1-KD and TFIIS4-KD and for end base frequency graphs of shorter IESs (26-36 bp). Equivalent graphs for DCL2/3-KD, DCL5-KD and PGM-KD can be seen in 15.

As we previously reported, the effects of either DCL2/3 knockdown or DCL5 knockdown on IES retention alone are quite modest^3,15. Even when their IES retention scores are summed, it is apparent that the knockdowns of DCL2/3 and DCL5 influence considerably fewer IESs, and to a lesser degree than the PGM knockdown. For DCL2/3/5, NOWA1/2 and TFIIS4, given that the gene knockdowns are efficient (as judged by depletion of their mRNAs or proteins; Figure S3B,D and¹²), we think the similarity of their IES retention score distributions indicates that the upper bounds of IES retention for these knockdowns have been approached. Consequently, these knockdowns may indicate most of the Paramecium IESs requiring scnRNAs and iesRNAs for their excision.

Though an interesting idea, we find no evidence of a proposed missing class of sRNAs^11,12 responsible for targeting and excision of additional IESs. Instead, we observe that after DCL2/3/5-KD the levels of IES-matching sRNAs are negligible (Figure 5C). We therefore infer that most IES excision in Paramecium does not require IES-targeting sRNAs. This is consistent with past studies reporting that most tested IESs were non-maternally controlled^27,28, i.e. that scnRNAs alone are insufficient for most IES excision^3,15. On the other hand, from IES retention following DCL2/3/5-KD we infer that many IESs require some scnRNAs or iesRNAs to cleanly remove all their copies (the exact number of affected IESs depending on the chosen retention score threshold, e.g. 19043 IESs are affected by DCL2/3/5-KD above a threshold of 0.1 retention).

An important observation from the DCL2/3/5 knockdown is the considerably stronger IES retention compared to the separate DCL2/3 and DCL5 knockdowns (Figure 2C); for example, at an IES retention score > 0.1, for DCL2/3/5-KD 19043 IESs are affected vs. 4815 IESs for DCL2/3-KD, and 5384 IESs for DCL5-KD. This suggests some co-operation between scnRNAs and iesRNAs in IES excision. From the DCL2/3/5 triple silencing, it is clear that multiple gene co-silencing may be necessary to gauge the full contributions of individual genes to IES excision. This is a critical consideration in gene function investigations in Paramecium tetraurelia with its abundance of paralogs that arose from successive whole genome duplications⁴².

As we previously proposed, IES excision/recognition efficiency and the need for particular genome reorganization proteins vary as a function of IES end bases and lengths¹⁵. A trend of increasing IES retention with IES length was shown for the knockdowns of DCL2/3, TFIIS4 and EZL1^3,11,12,15. It can now be seen that, with the exceptions of PGM-KD and DCL5-KD, most knockdowns leading to IES retention show this trend (Figure 2E,F). A new observation with respect to IES length is that IES retention for the knockdowns of DCL2/3/5, TFIIS4 and NOWA1/2 also has an obvious 10-11 base periodicity for ~45-130 bp IESs (Figure 2E), mirroring that of the IES length peaks (which are proposed to reflect the constraints of DNA twist on transposase excision¹). This signal is not apparent in knockdowns with weaker effects (DCL2/3-KD, DCL5-KD and PDSG2-KD), and only weakly visible in the knockdowns with stronger effects (PTCAF1 and EZL1). We previously also observed that the terminal IES base frequencies relative to IES length show a periodic, 10-11 bp signal (Figure 2A of 15). As in knockdowns of DCL2/3, DCL5 and PGM¹⁵, for all the new knockdowns we examined, the end base frequencies of longer IESs (e.g. 45-55 bp) vary in a similar manner with respect to their IES retention scores (Figure 2G–H; Figure S2), consistent with the idea that some IESs have a greater requirement for the products of these genes than others¹⁵.

Relationships among genes involved in genome editing inferred from IES retention

While similar IES retention distributions suggest different gene knockdowns might affect retention of particular IESs in a similar manner, examination of pairwise correlations is necessary to clearly establish this. We therefore determined all the possible pairwise correlations between IES retention scores of different gene knockdowns (Figure 3). From this figure, it can be seen that the DCL2/3/5, NOWA1/2, TFIIS4 (group 1), and histone modification genes (group 2) not only have the most similar IES retention score distributions within each of these two groups, but also the strongest IES retention score correlations. Though not as strong as those within the DCL2/3/5, NOWA1/2, TFIIS4, and histone modification groups, positive correlations are also evident between these groups (Figure 3; e.g. Spearman’s rank correlation coefficient, r_s=0.74 between DCL2/3/5-KD and PTCAF1-KD). The correlations between PGM knockdown and other gene knockdowns are generally weak — not much greater than the correlations between the control (ctrl3) and other gene knockdowns. Such weak correlations are an expected consequence of the ability of PiggyMac to excise most IES copies without scnRNAs/iesRNAs and without chromatin state changes brought about by Ezl1 and PtCAF1.

Since IES retention increases with IES length for many of the knockdowns influencing IES excision (see Figure 2E,F), we examined the retention score correlations between the abundant IESs corresponding to the first length peak (e.g. ≤ 35 bp; ∼15600 IESs; the IES length peaks can be seen in the background of Figure 2E) compared to longer IESs (e.g. ≥ 500 bp; ~380 IESs; patterns of correlations between these length intervals are intermediate to those in Figure 2E,F). For short IESs there are good correlations between all the possible pairs of IES retention scores from DCL2/3/5-KD, NOWA1/2-KD, TFIIS4-KD, PTCAF1-KD and EZL1-KD. This could be because the protein products of all the genes involved in these knockdowns are cooperating in the excision of short IESs. For the subset of short IESs, the right-skewed IES retention score distributions of these knockdowns also all resemble each other (diagonal of Figure 3B). We previously noted that DCL5-KD leads to greater retention of short IESs than DCL2/3-KD, i.e. iesRNAs are marginally more important for short IES excision than scnRNAs and short IESs are generally weakly epigenetically controlled^3,15. This is reflected by the stronger correlations in IES retention score for short IESs between DCL5-KD and the knockdowns of NOWA1/2, TFIIS4, PTCAF1 and EZL1, relative to the correlations between DCL2/3-KD and the same genes.

When focusing on longer IESs (e.g., those > 500 bp) the IES retention score correlations between the knockdowns of histone-associated genes, and also those between DCL2/3/5, NOWA1/2, TFIIS4 and histone-associated genes weaken substantially. Compared to their short counterparts the retention score distributions of longer IESs also look very different, with those of DCL2/3, DCL2/3/5, NOWA1/2, TFIIS4 closer to uniform, and those of the histone modification group closer to normal. The forms of the histone modification group IES retention score distributions now also become more like that of PGM-KD. These patterns are consistent with the chromatin state changes brought about by EZL1-KD and PTCAF1-KD exerting more influence on long IESs than short ones, with a much lesser requirement for iesRNAs, and generally some requirement for scnRNAs (i.e. long IESs are generally more strongly epigenetically controlled). At longer scales, the effect of PTCAF1-KD and EZL1-KD on IES retention (e.g. at ≥ 500 bp, r_s=0.58) does not correlate as well as at the shorter scale (e.g. at ≤ 35 bp, r_s=0.88), likely reflecting that the products of these genes do not have exactly the same consequences on chromatin state at this scale. The decrease in IES retention score correlations between either DCL2/3-KD or DCL2/3/5-KD and histone modification groups largely reflects that iesRNAs have little influence on long IES excision.

Figure 3. IES retention score correlations and regressions.

(A) For the purposes of regression, for each IES, negative control IES retention scores have been subtracted and zeroed when the subtraction yields a negative value, i.e. IRS_regression=max(IRS_experiment - mean(IRS_controls), 0) (controls=controls 1-4). This calculation was performed to reduce the effects of natural, background IES retention and systematic errors in IES retention score estimation as much as possible (see Figure S1, which examines these effects on regression). Hexagonal binning of IES retention scores was used to generate the plots in the lower triangular matrix. Red lines are for OLS regression, orange lines for LOWESS, and grey lines for ODR. IES retention score distributions given along the diagonal are shown to a maximum of 4000 IESs except for PGM-KD (maximum 8500 IESs). For each plot in the lower triangular matrix (M_i,j), Spearman’s rank correlation coefficient (r_s) is given in the corresponding position diagonally opposite in the upper triangular matrix (M_j,i). Note r_s=0.91 for biological replicates of NOWA1/2-KD (Figure S1B). (B) as in (A) but for short IESs (≤ 35 bp). IES retention score distributions given along the diagonal are shown to a maximum of 1500 IESs except for PGM-KD (maximum 3200 IESs). (C) as in (A) but for long IESs (≥ 500 bp). IES retention score distributions given along the diagonal are shown to a maximum of 50 IESs except for PGM-KD (maximum 75 IESs).

Relationships among genes involved in genome editing inferred from OES retention resemble those of long IESs

To date, information about the effects of gene silencing upon Paramecium’s other eliminated sequences (OESs), the remaining deleted DNA currently not annotated as IESs, has typically come from Southern blotting or dot blotting, e.g. 8,10,11,13,31, though more recent studies have also begun to investigate the overall effects of gene knockdown on these regions^11,12. Most gene knockdowns that lead to IES retention also lead to retention of transposons within OESs^8,10–13,31. For each of the knockdowns used to investigate IES retention, we examined OES retention by mapping the knockdown reads to genomic regions which do not match MDSs and IESs from the reference P. tetraurelia assemblies that were extracted from a PGM-KD DNA assembly (see Methods for more details). Since PGM-KD leads to the most IES retention, like previous reports^11,12, we assume that much of the remaining DNA eliminated from the MIC genome can be recovered from this knockdown.

Examination of the distributions of OES retention (Figure 4A–D), generated similarly to those for IES retention (Figure 2A–D), shows that while some features of OES retention resemble those of IES retention, e.g. the forms of the NOWA1/2-KD, TFIIS4-KD and DCL2/3/5-KD OES retention distributions are roughly similar, there are also important differences. OES retention due to DCL5-KD is very weak, hardly differing from that of a negative control, and DCL2/3/5-KD OES retention is weaker than that of DCL2/3-KD. Compared to the very strong IES retention following PGM-KD versus all other knockdowns (Figure 2D), PGM-KD OES retention is not as pronounced (Figure 4D). Like IESs, in part this may be because the strength of OES retention in other knockdowns is greater, but it is also possible that DNA cutting enzymes other than PiggyMac may be involved in OES elimination.

Figure 4. OES retention and OES retention correlations.

(A–D) For the sake of clarity, histogram bars with more than 2000 OESs have been truncated. The number of OESs binned in the left-most bars are as follows: (A) NOWA1/2-KD - 1238, TFIIS4-KD - 1230, DCL2/3/5-KD 2812; (C) DCL2/3-KD - 473, DCL5-KD - 3788, DCL2/3/5-KD - 2812; (D) Control - 3882; PDSG2-KD – 3062; PGM-KD – 709. (E) OES retention score correlations and regressions. For the purposes of regression, for each IES the mean of two negative control OES retention scores (control 1 and control 2) have been subtracted, and zeroed when the subtraction yields a negative value, i.e. ORS_regression=max(ORS_experiment – mean(ORS_controls), 0). This calculation was performed to reduce the effects of natural, background OES retention and systematic errors in OES retention score estimation as much as possible (but note that when OES retention scores are very low, some residual background retention may lead to some correlation, e.g. between DCL5-KD and control 3). Hexagonal binning of OES retention scores was used to generate the plots in the lower triangular matrix. Red lines are for OLS regression; orange lines for LOWESS and grey lines for ODR. For each plot in the lower triangular matrix (M_i,j), Spearman’s rank correlation coefficient (r_s) is given in the corresponding position diagonally opposite in the upper triangular matrix (M_j,i).

Examination of the correlations in OES retention in Figure 4E reveals a number of important effects: (i) between PGM-KD and other knockdowns, OES retention correlation is greater than general IES retention score correlations (Figure 3A); (ii) NOWA1/2-KD correlates most strongly with TFIIS4-KD and DCL2/3-KD; PTCAF1-KD and EZL1-KD correlate the best with each other, and moderately with DCL2/3-KD, DCL2/3/5-KD, TFIIS4-KD and NOWA1/2-KD; (iii) though a clear role for the Pdsg2 protein has yet to be defined, OES retention following PDSG2-KD correlates most strongly with PGM-KD; (iv) as OES length increases (e.g. > 4 kb) the PTCAF1-KD and EZL-KD OES retention scores become very narrowly distributed in contrast to those of DCL2/3-KD, NOWA1/2-KD and PGM-KD, suggesting that chromatin state changes brought about by PtCaf1 and Ezl1 may have a very general effect on OES retention. From these observations and those of the effects of previously reported analyses of the overall sequence complexity of retained OES DNA following gene silencing^11,12, it follows that OES elimination mostly does not need iesRNAs, whereas scnRNAs, together with TFIIS4 and Nowa1/2 are necessary for OES elimination. The stronger correlations in OES retention scores of other knockdowns relative to those of PGM-KD suggests greater cooperation between PiggyMac, scnRNAs and other protein factors in OES elimination than in general IES excision. Overall, the patterns of OES retention correlation for the different knockdowns resemble those of retention correlations for long IESs, suggesting most OESs may be treated like long IESs.

Nowa1/2 is involved in RNA scanning, and may mediate scnRNA- and iesRNA- targeted DNA cutting and excision

As the Nowa1/2 proteins have previously been implicated in removal of maternally controlled IESs³¹, but their precise role in Paramecium genome development has not yet been established, we wished to obtain additional insight into the function of these proteins, and their relation to other molecules involved in this process. Nowa1/2 proteins contain two distinct domains: an unstructured N-terminal domain consisting of alternating blocks of ‘FRG’ repeats and ‘GGWG’ repeats, and a structured C-terminal domain which is sufficient to route Nowa1 from the maternal MAC to the zygotic nuclei³¹. Nowa1/2's ‘FRG’ repeats resemble the ‘RGG’ boxes present in many hnRNP and other RNA-binding proteins, raising the possibility that they could facilitate interaction between Piwi-bound sRNAs and nascent transcripts on the old and the new MAC. One way Nowa1/2 may mediate these interactions is via its "GGWG" repeats, which resemble the "WG/GW" repeats found in Argonaute binding proteins^43,44 (Argonaute proteins are members of the Piwi protein superfamily⁴⁵). "WG/GW" repeat-containing proteins in Tetrahymena (Wag1 and CnjB) have been shown to interact with its Twi1 Piwi protein⁴⁶. In principle, Nowa1/2 could interact with Paramecium’s sRNA-bound Piwi proteins in both the old and the developing new MAC, and during the internuclear transport of the sRNAs. Like PtCAF1-KD, NOWA1/2-KD also substantially decreases H3K9me3 and H3K27me3 in the old and developing new macronuclei¹⁰.

To gain more insight into the role of Nowa1/2 in Paramecium genome reorganization we analyzed the effects of NOWA1/2-KD upon IES retention and the sRNA population during macronuclear development (Figure 5). We previously showed that in control and DCL5-KD cells, MAC genome-matching scnRNAs decrease during development, consistent with RNA scanning, and that there is a progressive increase in iesRNAs³. It can now also be seen that OESs give rise mainly to 25 nt scnRNAs and, compared to IESs, relatively few iesRNAs (Figure 5A; see also Figure S3G,H).

Figure 5. Analysis of the effects and properties of NOWA1/2-KD.

(A) Length distribution of small RNAs from control cells and NOWA1/2-KD cells (see Methods for details of histogram construction and normalization). Vector=vector backbone only. Figure 2SA shows the cytological characterization of the cells these sRNAs were isolated from. (B) Radioactively end-labelled total RNA from NOWA1/2-KD cells and control cells (cells fed with E. coli producing dsRNAs corresponding to the plasmid L4440 with no sequence target in the Paramecium genome; the gel portion showing the control sRNAs was previously shown in Figure 2A of 3 as "Control 1; Experiment 1"). (C) sRNA length distribution for DCL2/3/5-KD; the peak at 33 nt corresponds to an unknown sRNA family (e.g. 5'-GGAUCUAUCGUAUAGUGGUUAGUACCUGAGGCU-3') that matches a handful of locations in the MAC genome, which we have also be observed in other controls and knockdowns (e.g. from 8). See Figure 2SE for cytological characterization. (D) To examine the relationship between maternal control and IES retention scores (values determined in Table S1), maternal control scores were obtained from the maximal IES retention observed in Figure 6 of 28. Linear regression lines are shown along with Pearson's r and the two-tailed p-values for hypothesis testing with a null hypothesis that the regression slope is zero. (E) Retention of IESs tested by PCR for control and NOWA1/2-KD cells. NOWA1/2-KD IES retention scores are given to the right. Additional bands indicate the presence of heterodimers of ssDNAs with excised and unexcised IESs. (F) Electrophoretic mobility shift assay with MBP (maltose-binding protein), MBP-Nowa1 (alternating lanes) and oligonucleotides, visualized on a native polyacrylamide gel.

Since it is known that smaller IESs can be nested within larger IESs⁴⁷, we searched for IESs embedded within OESs. Visual inspection of DNA reads mapping to the assembled OESs from different gene knockdowns confirms the presence of shorter IESs in OESs, often with matching iesRNAs. We previously noted that, though both scnRNAs and iesRNAs are enriched at IES ends, scnRNAs overlap IES ends, whereas iesRNAs map very specifically within the limits of IESs³. Given the patterns of OES retention shown in Figure 4E, we therefore suggest that iesRNAs primarily facilitate precise IES excision, whereas scnRNAs facilitate both precise IES excision and imprecise OES elimination.

Compared to control cells, scnRNAs accumulate and iesRNAs are reduced after NOWA1/2 silencing (Figure 5A, B). scnRNA accumulation is consistent with inhibition of RNA scanning when NOWA1/2 is silenced, i.e. MDS scnRNAs are not removed. We do not think that NOWA1/2-KD affects RNA scanning by interfering with transport of scnRNA-bearing Piwis (Ptiwi01/09⁹), as we observed no differences in localization of GFP-tagged Ptiwi09 with and without silencing NOWA1/2 (Figure S4). On the other hand, the pronounced mobility shift in an electrophoretic mobility shift assay upon addition of synthetic MBP-Nowa1 protein (Figure 5F) indicates Nowa1 may interact with RNA, RNA:DNA hybrids, and, to a lesser degree, DNA. Hence NOWA1/2-KD-induced interference in RNA scanning in the old MAC could be due to disruption of interactions normally mediated by Nowa1/2 between nucleic acids and sRNA-carrying Ptiwi proteins.

From the electrophoretic visualization of sRNAs (Figure 5B) it can be seen that iesRNAs are reduced following NOWA1/2-KD (since Figure 5A shows only proportions, it is not possible to deduce whether absolute levels of sRNAs change from this figure). A decrease in iesRNAs was also observed following DCL2/3 co-silencing, and, to a lesser degree for the independent DCL2 and DCL3 silencings³. For the lower quantity of iesRNAs in NOWA1/2-KD cells, there are two possible explanations, both of which may be valid: (i) Nowa1/2 is directly involved in iesRNA production or stability; (ii) Nowa1/2 depletion inhibits scnRNA-dependent IES excision due to the inhibition of RNA scanning, which in turn inhibits iesRNA production. Teasing apart these possibilities is a possible future experimental avenue to gain additional insight into the function of the Nowa1/2 proteins.

NOWA1/2 silencing was previously reported to lead to retention of maternally controlled IESs, and indeed maternally controlled IES retention levels (following microinjection of complementary IESs into the old MAC) correlate well with the IES retention scores (Figure 5D). However, following NOWA1/2 silencing we also observe low levels of retention of some non-maternally controlled IESs by PCR analyses and in DNA-seq data (Figure 5D, E; Table S1). As the non-maternally controlled IESs that are weakly retained following NOWA1/2-KD are also weakly retained following DCL2/3-KD, we infer that they are weakly epigenetically controlled (Table S1). Therefore, due to the sensitivity of IES retention detection in DNA-seq data, some "non-maternally controlled" IESs may actually be weakly epigenetically controlled (it is also possible that with more sensitive methods these IESs would be classified as weakly maternally controlled). From Figure 3, we infer that Nowa1/2 is not only involved in scnRNA-dependent IES excision (epigenetically controlled IES excision), but also iesRNA-dependent IES excision, since IES retention scores of NOWA1/2-KD are positively correlated with those of both DCL2/3-KD and DCL5-KD, and strongly correlated with those of DCL2/3/5-KD.

From these observations of the effects of NOWA1/2-KD, we propose that Nowa1/2's WG/GW (Argonaute-binding) repeats could mediate IES targeting by both Ptiwi01/09-bound scnRNAs (both old and new MAC) and iesRNAs bound to other Ptiwi proteins⁴⁸ in the new MAC. The apparent failure to remove MAC genome-matching scnRNAs after NOWA1/2 silencing (Figure 5A) is consistent with the involvement of Nowa1/2 in RNA scanning in the old macronucleus, as expected if these proteins mediate pairing of scnRNAs with nascent maternal non-coding transcripts. In principle, Nowa1/2 proteins could interact with such transcripts via their FRG (RNA-binding RG boxes) repeats. From the correlations in Figure 4, we infer that Nowa1/2 proteins may operate predominantly in concert with scnRNAs and not iesRNAs during OES elimination.