A MILI-independent piRNA biogenesis pathway empowers partial germline reprogramming.

In mice, the pathway involving PIWI and PIWI-interacting RNA (PIWI-piRNA) is essential to re-establish transposon silencing during male-germline reprogramming. The cytoplasmic PIWI protein MILI mediates piRNA-guided transposon RNA cleavage as well as piRNA amplification. MIWI2's binding to piRNA and its nuclear localization are proposed to be dependent upon MILI function. Here, we demonstrate the existence of a piRNA biogenesis pathway that sustains partial MIWI2 function and reprogramming activity in the absence of MILI.

The erasure of DNA methylation during male mammalian germline reprogramming and the concomitant de-repression of transposons initiate the process of re-silencing through the production of primary piRNAs1. In mice, these piRNAs are bound to the cytoplasmic PIWI protein MILI, which is a piRNA-directed RNA endonuclease that mediates piRNA amplification of secondary piRNAs to ensure transposon repression2 (Supplementary Fig. 1). These secondary piRNAs through sequence complementarity guide the PIWI protein MIWI2 to direct de novo DNA methylation of LINE1 and IAP elements1–5 (Supplementary Fig. 1). The loss of MILI or MIWI2 results in defective de novo DNA methylation and depression of LINE1s and IAPs causing meiotic arrest3,5,6. In addition to meiotic arrest, Miwi2 mice progressively lose germ cells and are fully aspermatogenic by 9 months2,3. We found that Mili mice present a milder germ cell loss phenotype with ~50% tubules being spermatogenic at 1 year of age (Fig. 1a,b). This discrepancy in phenotype indicates that MIWI2 may have a MILI-independent role during fetal germ cell reprogramming. We therefore FACS-purified undifferentiated spermatogonia and performed whole genome bisulfite sequencing. While the loss of MILI and MIWI2 predominantly affects LINE1 and IAP methylation4,5,7 (Supplementary Fig. 2), we found in accordance with our hypothesis that Miwi2-deficiency has a more severe impact on genome re-methylation than Mili-disruption (Fig. 1c and Supplementary Fig. 2). This contrasts with a previous report8; the discrepancy likely arises from the fact that in the previous study Mili and Miwi2 spermatogonia were differentially isolated with varying degrees of somatic cell contamination8. These somatic cells have normal genomic methylation levels and thus a differential degree of contamination confounds a quantitative comparative analysis. Here we observed that the methylation of 1704 loci were dependent upon MIWI2 function, whereas 258 loci showed selective dependence upon MILI (Fig. 1c). To determine if the methylation of these loci could be attributed to the piRNA pathway, we analyzed whether piRNAs bound to MILI or MIWI2 fetal gonad RNPs map to these loci. Indeed, piRNAs from both RNPs were complimentary to these loci (Fig. 1d). We next analyzed uniquely mapping piRNAs, which would exclude transposon-associated piRNAs and found that unique piRNAs from the respective RNPs map to both MIWI2- and MILI-dependent loci (Fig. 1e). In summary, Miwi2-deficiency has a greater impact on spermatogenesis and fetal reprogramming in comparison to the loss of Mili.

Figure 1

Loss of MIWI2 results in more severe physiological and molecular phenotype than MILI-deficiency.

(a) Percentage of aspermatogenic tubules in Mili and Mili 1-year-old mice. Error bars represent a standard deviation of the mean (n=5 animals). (b) Representative images of hematoxylin and eosin stained testis cross-sections of Mili and Mili 1-year-old mice. (c) Comparisons of the percentage of CpG methylation in wild type (WT) and Mili, WT and Miwi2 as well as Mili- and Miwi2 undifferentiated spermatogonia are shown. Blue dots represent significantly differentially methylated regions (DMRs). (d & e) piRNAs (d) and uniquely mapping piRNAs (e) from MIWI2 and MILI ribonucleoprotein complexes (RNPs) mapped to loci whose methylation is dependent upon MIWI2 and MILI, respectively. Positive (red) and negative (blue) values indicate sense and antisense piRNAs, respectively. An averaged value from biological duplicates is shown (n=2 animals).

The physiological and molecular discrepancy between Miwi2- and Mili-deficiency could indicate that a non-canonical piRNA biogenesis and partial piRNA-dependent MIWI2 function exists in the absence of MILI. To test this hypothesis, we generated total small RNA libraries from wild type and Mili E16.5 fetal testis. Albeit at low levels, a clear presence of piRNAs is found in the Mili libraries (Fig. 2a and Supplementary Fig 3a). Importantly, these are bona fide piRNAs as all the piRNAs found in absence of MILI are present in wild type fetal gonads (Fig. 2b). LINE1 and IAP piRNAs are generated in greatly reduced quantities and as expected the hallmarks of piRNA amplification are absent in the Mili libraries (Supplementary Fig. 3b-e). In the Mili mutant, piRNAs cognate to the 1704 MIWI2-dependent loci were still present but at lower levels (Fig. 2c and Supplementary Fig. 3f), with a 2- and 3-fold reduction in multi or unique-mapping piRNAs observed, respectively (Supplementary Fig. 3g). In contrast, in the Mili libraries there was a drastic 18-fold reduction in the uniquely mapping piRNAs to the 256 loci whose methylation is dependent upon MILI function (Supplementary Fig. 3g). Since piRNAs exist in Mili fetal gonocytes albeit in reduced quantities, this should result in the partial nuclear localization of MIWI2. We next analyzed the localization of MIWI2 in E16.5 Mili fetal gonocytes by confocal immunofluorescence using a rabbit polyclonal antibody (Fig. 2d and Supplementary Fig. 4a). This analysis revealed that a portion of MIWI2 indeed retains nuclear localization (Fig. 2d). Quantification of the signal in wild type and Miwi2 fetal testis revealed that the anti-MIWI2 antibody gave a non-negligible level of non-specific background staining (Supplementary Fig. 4b). To independently confirm this finding, we generated an N-terminal HA epitope tagged MIWI2 (HA-MIWI2) mouse allele (Miwi2) (Supplementary Fig. 5), permitting the use of high quality monoclonal HA antibodies. Confocal immunofluorescence of HA-MIWI2 in Mili fetal gonocytes confirmed the partial MIWI2 nuclear localization with dramatically reduced levels of non-specific background (Fig. 2e and Supplementary Fig. 4c,d).

Figure 2

piRNAs are produced and MIWI2 is partially localized to the nucleus in the absence of MILI.

(a) Length distribution of piRNAs in wild type (WT) and Mili E16.5 fetal testis. (b) piRNA expression analysis in WT and Mili E16.5 fetal testis. (c) Unique piRNAs from WT and Mili E16.5 fetal testis mapped to loci whose methylation is dependent upon MIWI2 and MILI, respectively. Positive (red) and negative (blue) values indicated sense and antisense piRNAs, respectively. An averaged value from biological duplicates is shown in panels a-c (n=2 animals). (d & e) Immunofluorescence with anti-MIWI2 antibody (d) and anti-HA antibody (e) on fetal testis of the indicated genotypes with a single fetal gonocyte shown. Scan-line profile plots (bottom panel) represent relative fluorescence intensity of MIWI2 (d) and HA (e) (green) staining in the nucleus (blue) along the analyzed trajectory (represented as arrow in d and e).

The selective loss of MILI’s endonuclease activity attenuates piRNA amplification, but leaves primary piRNA biogenesis intact, resulting in a MIWI2 RNP that can partially function to silence IAP2, whereas the complete loss of MILI function was thought to abrogate all piRNA biogenesis and MIWI2 reprogramming activity1. Here, we show that a non-canonical primary piRNA biogenesis exists in the absence of MILI which results in partial MIWI2 nuclear localization and reprogramming activity. This residual reprogramming activity is supported by the fact that more severe methylation and spermatogenesis defects are observed in Miwi2 over Mili mice. The piRNA biogenesis in Mili mice may resemble the phased primary piRNA pathway identified in Drosophila, Bombyx mori and mouse, where the binding of PIWI to cleaved transcripts from the amplification cycle coupled with Zucchini-mediated endonucleolytic cleavage results in additional piRNAs biogenesis9–12. While phased piRNA biogenesis in Drosophila can occur in the absence of piRNA amplification13, we could not detect phasing in the Mili piRNA population (data not shown), likely given the extensive 3’ trimming that occurs during mouse piRNA biogenesis14. In summary, here we show that mouse fetal primary piRNA biogenesis can exist in the absence of MILI and this residual piRNA pathway contributes to germline reprogramming.

Accession codes

Small RNA sequencing data that support the findings of this study have been deposited in www.ebi.ac.uk/ena with the accession code PRJEB19792. WGBS data that support the findings of this study have been deposited in www.ebi.ac.uk/arrayexpress with the accession code E-MTAB-5561. MIWI2 and MILI RNP-associated piRNA libraries were used from previously published study2.

Note: Any Supplementary Information and Source Data files are available in the online version of the paper. The imaging data that support the findings of this study are available from the corresponding author upon reasonable request.

Online Methods

Mouse strains and experimentation

The Mili and Miwi2 mouse alleles were previously produced in the O’Carroll laboratory2,15 and maintained on a mixed 129/C57Bl/6 genetic background. The Miwi2-tdTomato transcriptional allele was used to isolate spermatogonia (Carrieri C and O’Carroll D, unpublished observations). The Miwi2 allele was created using CRISPR/Cas9 gene editing technology16. The sequence encoding a HA-tag was inserted after the starting ATG codon of Miwi2 to generate the Miwi2 allele that results in an N-terminal HA-MIWI2 fusion protein. An sgRNA (5’ ACCCAGGAACATGAGTGGAC) targeting Miwi2 in vicinity of the starting ATG codon was used. A DNA oligonucleotide with HA and BamHI restriction site flanked by 60bp of 3’ and 5’ Miwi2 genomic sequence was used as a template for homology directed repair. The sgRNA (12,5 ng/ul), DNA oligonucleotide (12,5 ng/ul), as well as Cas9 mRNA (25 ng/ul) were injected into the cytoplasm of two pronuclei stage C57Bl/6 mouse embryos. Manipulated embryos were cultured overnight and embryos that had proceeded to the 2 cell stage were transferred into the oviduct of pseudopregnant females. The presence of Miwi2 allele in mice was screened by digesting PCR products amplified by primers (Fw 5’ACAGCCACACCGTCTCTTTT; Rv 5’CAGGATAGCCAAAGGAAGGA) flanking Miwi2 ATG with BamHI. The Miwi2 allele was also confirmed by sequencing. For the analysis of progressive germ cell loss, Mili and Mili mice were aged to one year, sacrificed and testes collected. Each genotype group contained five animals to ensure adequate power to detect and describe differences between the genotypes. E16.5 fetal testes were used for immunofluorescence analysis and small RNA library preparation. No blind experiments or randomization of samples or animals were employed in this study. All mice used in this study were male; and bred and maintained in EMBL Mouse Biology Unit, Monterotondo, and subsequently in the Centre for Regenerative Medicine, Edinburgh. All procedures were done in accordance to the current Italian legislation (Art. 9, 27. Jan 1992, nu116) under license from the Italian health ministry or the UK Home Office regulations, respectively.

Spermatogonia isolation

Isolated postnatal day 14 testes from mice of the respective genotypes that also carried the Miwi2-tdTomato transcriptional allele in heterozygosity were dealbulginated and digested with collagenase (0,5 mg/ml, Sigma) at 32°C for 10 minutes shaking in 1ml of DMEM media (Life technologies) supplemented with Pen/Strep (Life technologies), NEAA (Life technologies), sodium pyruvate (Life technologies) and sodium lactate (Sigma). After digestion seminiferous tubules were collected by centrifugation at 1000 rpm for 5 min. Single cell suspension was prepared by further digesting seminiferous tubules in 1 ml of 0,05% trypsin (Life technologies) and DNase (Sigma), added to 0,05 mg/ml of final concentration, at 32°C for 8-10 minutes shaking. After digestion trypsin was neutralized by adding 100 μl of FCS. Cells were spun and resuspended in PBS with 3% FCS. Germ cells were stained with following antibodies: c-Kit 1:800 (eBioscience 25-1171), CD45-biotin 1:400 (eBioscience 13-0451), CD51-biotin 1:100 (Biolegend 104104), streptavidin-qDot 1:50 (eBioscience 93-4317), CD9 1:50 (eBioscience 11-009) in PBS with 3% of FCS and 0,01% sodium azide. SYTOX Blue (Life technologies) was used as a living dye. Target population was selected by excluding cell doublets and CD45pos and CD51pos cells by initially choosing to analyze cells with low side scatter and negative for live cell dye SYTOX Blue. Among CD45neg and CD51neg cells, a target population was identified as c-Kitneg and having a live fluorescence of tdTomato and being CD9pos. Cells were sorted with an 85μm nozzle using a FACSAria II SORP (BD Biosciences). This strategy results in the isolation of undifferentiated spermatogonia with a purity of greater than 98%.

Histology

Testes were fixed in Bouin’s fixative (Sigma) overnight at 4°C and paraffin embedded. 8 µm thick sections were stained with hematoxylin and eosin using routine methods.

Immunofluorescence

Testes were freshly embedded into OCT compound (Sakura), 8 μm sections were cut, fixed in 4% paraformaldehyde (Sigma) 10 minutes at room temperature and then permeabilized for 10 minutes at room temperature in 0.1% Triton-X100 (Sigma). Subsequently sections were blocked for 30 minutes at room temperature in 10% normal donkey serum (Sigma), 1% BSA (Sigma) and 0.1M glycine (Sigma). Primary antibody incubation was done overnight at 4°C in the blocking buffer. The following antibodies were used at indicated dilutions: rabbit polyclonal anti-MIWI2 1:2002, rabbit monoclonal anti-HA 1:500 (Cell Signaling 3724). Appropriate donkey anti-rabbit AlexaFluor488 secondary antibody (Life technologies) was used in a dilution 1:1000. DAPI (5 μg/μl) (Life technologies) was used to stain DNA. Samples were mounted with ProLong Gold antifade reagent (Life technologies). Leica TCS SP5 confocal microscope was used to acquire images. Five images were acquired per each biological replicate per genotype. After acquisition images were processed with ImageJ and Adobe Photoshop CS5 computer programs. All images corresponding to the same experiment were acquired and processed applying the same settings. To determine MIWI2 and HA immunofluorescence intensity in the nuclei of gonocytes, fetal testes cross-sections were stained with MIWI2 and HA antibodies as described above. Fetal testes from two embryos were used per each genotype analyzed as biological duplicates, except a biological triplicate was used for Miwi2 fetal testis. All the data was collected in 2-3 independent immunofluorescence experiments. Acquired stacks of images were merged using ImageJ computer program. ROI (region of interest) manager tool in ImageJ was used to quantify relative fluorescence intensity in the nuclei, where ROIs were manually defined, avoiding cytoplasmic signal.

Bisulfite conversion

For bisulfite conversion DNA was isolated from sorted cells as described above. Whole genome bisulfite sequencing libraries were generated using a post-bisulfite adaptor tagging (PBAT) method as previously described17 using 10 cycles of PCR amplification. Libraries were sequenced using Illumina HiSeq 2000. Three biological replicates (cells independently isolated from three animals) were generated per genotype to ensure adequate power to detect statistically significant differences between the respective genotypes.

Bisulfite sequencing analysis

Raw sequence reads were trimmed to remove both poor quality calls and adapters using Trim Galore (v0.4.1, www.bioinformatics.babraham.ac.uk/projects/trim_galore/, Cutadapt version 1.8.1, parameters: --paired). Trimmed reads were aligned to the mouse genome (GRCm38) in paired-end mode to be able to use overlapping parts of the reads only once. Alignments were carried out with Bismark v0.14.418 with the following set of parameters: paired-end mode: --pbat. Reads were then deduplicated with deduplicate_bismark selecting a random alignment for position that was covered more than once. CpG methylation calls were extracted from the deduplicated mapping output ignoring the first 6 bp of each read to reduce the methylation bias typically observed in PBAT libraries using the Bismark methylation extractor (v0.14.4) with the following parameters: a) paired-end mode: --ignore 6 --ignore_r2 6; b) single-end mode: --ignore 6. CpG methylation calls were analysed using R and SeqMonk software (www.bioinformatics.babraham.ac.uk/projects/seqmonk/). 50 adjacent CpG running window probes were generated and percentage of methylation determined for probes containing at least 5 reads and 3 CpG on the pooled replicate data. The methylation level was expressed as the mean of individual CpG sites. Graphing and statistics were performed using Seqmonk and RStudio.

Transposon analysis in bisulfite sequencing data

Repeat locations for a pre-defined set of repeat classes of interest were extracted from the pre-masked repeatmasker libraries (mm10 - Dec 2011 - RepeatMasker open-4.0.5, http://www.repeatmasker.org/species/musMus.html). Repeat instances overlapping annotated genes in the Ensembl gene set were removed to avoid mixing signals from genic expression with specific expression of repetitive sequences. Methylation levels at the repeat instances were quantitated by summing up all methylation calls and non-methylation calls for all instances of each class of repeat and calculating the percentage of methylated Cs over all Cs.

Small RNA libraries

Isolated E16.5 fetal testes were homogenized in 1 ml of Qiazol (QIAGEN). Extracted total fetal testes RNA was used to generate small RNA libraries using NEBNext Multiplex Small RNA Library Prep Set for Illumina (Set 1) (NEB) following manufacturer’s instructions with 20 cycles of PCR amplification. Small RNA libraries were generated from two biological replicates (RNA isolated from fetal testes from two animals) per genotype to ensure adequate confidence in observations.

Small RNA sequencing analysis

All samples were initially aligned against Rfam19 in order to filter out tRNA sequences (tRNA hits with an alignment identity score > 90% were excluded from the rest of the analysis). Analysis of the filtered samples was then performed using SequenceImp20. Input reads were first trimmed from the 3’ adapter with reaper (using default configuration for read geometry without barcode) and de-duplicated with tally, which are both part of the Kraken suite of tools20. The length distribution of all cleaned reads between 18 and 32 nucleotides was recorded in order to check for depletion of piRNA sequences between the wild type and Mili conditions. Cleaned reads were later filtered by length (retaining only 24-32 nucleotide long sequences) and aligned against the Mouse genome (Ensembl release 66) allowing up to 2 mismatches and reporting up to 20 hits per sequence, when analyzing for all mappers. In the case of uniquely mapped sequences, bowtie call from the SequenceImp pipeline was tweaked using the parameter -m = 1 (parameter -k was set to 1 for both the unique and all mappers cases). BAM output files from the alignment step were intersected with BED files containing the coordinates of 1704 and 258 loci, whose methylation is dependent on MIWI2 and MILI, respectively. PiRNA counts within each locus were calculated as the average number of fragments aligning against the locus, divided by the size of the locus region in 1 kb units. Expression densities were limited within the interval (-100, 100) in order to filter out outliers from 3 overexpressed loci and thus increase densities resolution for all loci. With regards to the piRNA differential expression analysis, a custom database of all 26-31 nucleotide long unique sequences found across all wild type and Mili replicates was initially built. Each sample replicate was then aligned against this database and expression levels of all matching sequences were quantified between the two conditions. As for the quantification of LINE1 and IAP repeats, the analysis was performed using the ‘features’ step of SequenceImp for repeat elements, allowing up to 3 mismatches and correcting the read counts to the number of genome mapping reads. The ping-pong signatures and 1U-10A content of the LINE1 and IAP elements were also calculated as part of this step. In all cases, normalization was performed based on the total number of reads of transcripts that remained unchanged between the two conditions. Specifically, the types of transcripts that were used for normalisation (based on the official Ensembl genebuild annotation) were: miRNAs, rRNAs, snRNAs, snoRNAs, processed_transcripts, aa_tRNAs, Mt_tRNAs, other_tRNAs, miscRNAs and RNA_repeats. Significance was assessed using the BootstRatio algorithm21.

Code availability

Computer code used to analyze small RNA sequencing and whole genome bisulfite sequencing data is available upon request from A.J.E./D.V. and R.B./W.R., respectively.