Animal use and care

Oocytes were collected from monkeys, mice, rats, golden hamsters, rabbits, guinea pigs, Chinese hamsters, dogs, zebrafish, pigs, and goats as described below. Mouse, rat, and zebrafish oocytes were obtained in compliance with the guidelines of the Center for Excellence in Molecular Cell Science at the Chinese Academy of Sciences (CAS). Monkey oocytes were collected following the guidelines of the Institute of Neuroscience, CAS. Chinese hamster oocytes were obtained in compliance with the guidelines of Shanxi Medical University.

Monkey oocyte collection

Monkey oocytes were collected as previously described [16]. Superovulation was triggered in adult female Macaca fascicularis monkeys by intramuscular injection of 15 IU of GONAL-F recombinant human follitropin alfa (rhFSH) (Merck Serono, Germany) twice daily for 8 days and 1000 IU of human chorionic gonadotropin (hCG) (Sigma, USA) on the ninth day. Oocytes were collected by laparoscopic follicular aspiration 36 h after hCG injection. Cumulus cells were removed by adding hyaluronidase (2 mg/ml) and gently pipetting. Cumulus-free high-quality oocytes in metaphase II (MII) were visually inspected using a microscope, then cultured in HECM-9 medium at 37 °C with 5% CO2 for ~ 2 h. The MII oocytes were transferred into a droplet of acidic Tyrode’s solution (Sigma, T1788) to remove the zona pellucida, then washed twice with DPBS (Sigma, D5652). Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C until use.

Mouse, rat, and golden hamster oocyte collection

Five-week-old female mice (C57BL/6), 6-week-old Sprague Dawley (SD) rats, and 5-week-old LVG strain golden hamsters were intraperitoneally injected with pregnant mare serum gonadotropin (PMSG) (Ningbo, China) (7.5 IU, 20 IU, and 20 IU, respectively). At 48 h after PMSG injection, mice were intraperitoneally injected with 7.5 IU hCG (Ningbo, China). At 56 h after PMSG injection, rats and golden hamsters were injected with 20 IU hCG. At 20 h after hCG injection, cumulus-oocyte complexes were collected from oviducts and placed in M2 medium (Sigma, M7167). Cumulus cells were removed by adding hyaluronidase (2 mg/ml; Sigma, H3506) and gently pipetting, and MII oocytes were selected and washed twice with DPBS. Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C prior to further use.

Rabbit and guinea pig oocyte collection

Female rabbits weighing 2.5 kg (New Zealand) and 4-week-old female guinea pigs were provided by the Shanghai Songlian Lab Animal Farm. Rabbits and guinea pigs were administered 150 IU PMSG via intramuscular injection and 20 IU PMSG via intraperitoneal injection, respectively. At 56 h after PMSG injection, guinea pigs were intraperitoneally injected with 20 IU hCG; at 72 h after PMSG injection, rabbits were injected with 150 IU hCG via the ear vein. At 20 h later, animals were sacrificed, the ovaries were cut into pieces, and cumulus-oocyte complexes were released from the ovaries into M2 medium. Cumulus cells were removed by adding hyaluronidase (2 mg/ml) and gently pipetting. GV oocytes were transferred into acidic Tyrode’s solution to remove the zona pellucida, then washed twice with DPBS. Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C until use.

Chinese hamster oocyte collection

Four-week-old female Chinese hamsters were sacrificed, then the ovaries were cut into pieces and cumulus-oocyte complexes were released from the ovaries into M2 medium. Cumulus cells were removed by adding hyaluronidase (2 mg/ml) and gently pipetting. GV oocytes were transferred into fresh M2 medium. Oocytes with normal form, bright color, and no cumulus cells were selected and each was transferred to a separate PCR tube containing < 0.5 μl DPBS, then stored at − 80 °C until use.

Dog oocyte collection

MII oocytes were collected from 12- to 24-month-old beagles on the estimated ovulation day through continuous observation. The collected oocytes were transferred into M199 medium containing 25 mM HEPES. Cumulus cells were removed by adding hyaluronidase (2 mg/ml) with gentle pipetting. The oocytes were then washed twice in DPBS and each oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C until use.

Zebrafish oocyte collection

Female AB strain zebrafish were anesthetized with 0.16 mg/ml Tricaine for 2–3 min, then transferred to phosphate-buffered saline (PBS) in a 35-mm plastic petri dish. Eggs were removed by gently pressing on the belly of each fish. Each egg was transferred to a separate 1.5-ml tube and stored at − 80 °C.

Pig oocyte collection

Pig ovaries were obtained from a local slaughterhouse and transported to the laboratory in 0.9% NaCl (containing 100 IU/ml penicillin) within 2 h by Nanjing Zhushun Biotechnological Co., LTD. Antral follicles on the ovary surfaces were released by aspirating with an 8-gauge needle attached to a 10-ml syringe and washing three times in TCM 199 medium (Gibco, 11,150,059). Good quality cumulus-oocyte complexes with competent oocytes containing three cumulus cell layers were selected and cultured in TCM 199 medium supplemented with 0.1% polyvinyl alcohol (PVA), 3.05 mM D-glucose, 0.91 mM sodium pyruvate, 0.5 IU/mL follicle-stimulating hormone (FSH), 0.5 IU/mL luteinizing hormone (LH), 10 ng/mL epidermal growth factor (EGF), 75 μg/mL penicillin, 50 μg/mL streptomycin, 0.57 mM cysteine, and 10% follicular fluid. The cumulus-oocyte complexes were cultured in an incubator at 39 °C with 5% CO2 for 44 h. Cumulus-oocyte complexes were then transferred to M199 medium and cumulus cells were removed by adding hyaluronidase (2 mg/ml) and gently pipetting. MII oocytes with the first polar body were selected and transferred into acidic Tyrode’s solution to remove the zona pellucida, then washed twice with DPBS. Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at -80 °C until use.

Goat oocyte collection

Goat ovaries were obtained from a local slaughterhouse and transported to the laboratory in 0.9% NaCl (containing 100 IU/ml penicillin) within 2 h by Nanjing Zhushun Biotechnological Co., LTD. Antral follicles on the surface of ovaries were released by aspirating with an 8-gauge needle attached to a 10-ml syringe and washed three times in TCM 199 medium (Gibco). Good quality cumulus-oocyte complexes with competent oocytes containing three cumulus cell layers were selected and cultured in TCM 199 medium supplemented with 10% fetal bovine serum (FBS, Gibco), 1 μg/ml 17 β-estradiol, 24.2 mg/L sodium pyruvate, 0.5 IU/ml FSH, 0.5 IU/ml LH, 10 ng/ml EGF, 100 mM cysteamine, and 200 mM cystine. The cumulus-oocyte complexes were cultured in an incubator at 38.5 °C with 5% CO2 for 24 h. Cumulus-oocyte complexes were then transferred into M199 medium and cumulus cells were removed by adding hyaluronidase (2 mg/ml) and gently pipetting. MII oocytes with the first polar body were selected and transferred into acidic Tyrode’s solution to remove the zona pellucida, then washed twice with DPBS. Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C until use.

Human oocyte collection

Pre-matured oocytes in the germinal vesicle (GV) stage were donated by patients undergoing a cycle of intracytoplasmic sperm injection (ICSI). Oocytes were cultured in Continuous Single Culture medium (CSC, Irvine Scientific, CA, USA) with 10% Serum Substitute Supplement (SSS, Irvine Scientific) at 37 °C with 5% CO2. The presence of the germinal vesicle was confirmed in each oocyte under a macroscope 24 h later. The GV oocytes were then transferred into a droplet of acidic Tyrode’s solution (Sigma, T1788) to remove the zona pellucida, then washed twice with Dulbecco’s phosphate-buffered saline (DPBS) (Sigma, D5652). Each naked oocyte was transferred to a separate PCR tube containing < 0.5 μl DPBS and stored at − 80 °C prior to use.

Western blot

Pooled MII oocyte samples (100 mouse oocytes, 75 rat oocytes, and 44 golden hamster oocytes) were each lysed in 8 μl of RIPA buffer (Beyotime, P0013C) with 1/100 proteinase inhibitor (Sigma, P8340) added. Samples were incubated on ice for 10 min, then 2 μl of 4 × protein loading buffer was added. Samples were incubated at 95 °C for 5 min, then resolved with 8% SDS-PAGE. Samples were then semi-dry transferred from the gel onto a polyvinylidene fluoride membrane and stained with anti-Dicer (Abcam, ab259327, 1:1000) or anti-tubulin (Proteintech, 66,031–1-Ig, 1:2000) antibodies overnight at 4 °C. Finally, membranes were incubated with appropriate horseradish peroxidase-conjugated secondary antibodies at room temperature for 1 h and developed using Eblot.

NaIO4 oxidation of oocyte RNA

Total RNA was extracted from five MII or five GV oocytes from each species with TRIzol Reagent (Invitrogen, 15,596,026) after the addition of spike-in oligos. The RNA pellets were dissolved in 4 μl water, and 2 μl of each resulting sample was saved as a control. The remaining 2 μl of each RNA sample was treated with an oxidation reaction mixture, containing 1.5 μl 100 mM NaIO4 (Sigma, S1878), 2 μl of 5 × borate buffer (150 mM borax and 150 mM boric acid at pH 8.6), and 4.5 μl water. The reactions were incubated in the dark for 30 min at room temperature. RNA was then precipitated by adding 600 μl 100% ethanol with 25 μg linear acrylamide and 30 μl sodium acetate (3 M, pH 5.2) and incubating at − 30 °C for 30 min. Samples were centrifuged at 16,200 × g for 15 min at 4 °C. Pellets were washed twice with 75% ethanol and dissolved in 2 μl water for library construction.

Generation of polyclonal antibodies for PIWIL3 and PIWIL1

Polyclonal antibodies used for immunoprecipitation were produced in rabbits by GL Biochem (Shanghai, China). Peptides with the following sequences were chosen as antigens to produce antibodies against PIWIL3 in guinea pig, rabbit, dog, pig, and goat, respectively: CQDLVVNTREKLRHVRHSKTG, C-EAQAVRALEAPQLHAREAE, CRHLEPQRHLEPQRHLEPQR, CAPEPAEPQPSEVARASEVT, and C-RTDPPLSFADLVRRGTAAQ. A peptide with the sequence C-HDLGVNTRQNLDHVKESK was used as the antigen to produce the anti-PIWIL1 antibody. Cysteine residues were included at the N terminal of the synthetic peptides for later coupling. The rabbits used for antibody production underwent seven rounds of immunization with the appropriate antigen prior to blood collection. Antibodies were then purified with peptide affinity chromatography.

Immunoprecipitation of antibodies in transfected HEK293T cells

HEK293T cells were transiently transfected with plasmids encoding Flag-tagged PIWIL1, PIWIL2, or PIWIL3 using Lipofectamine 2000 (Invitrogen, 11,668,019) following standard procedures. Cells were harvested 36 h after transfection and incubated in lysis buffer (50 mM Tris–HCl [pH 7.4], 1% Polyoxyethylene (40) nonyl phenyl ether [NP-40], 150 mM NaCl, 0.5 mM Dithiothreitol [DTT], and 1 × proteinase inhibitor cocktail [Sigma]) for 10 min on ice. Cell lysates were centrifuged at 12,000 × g for 10 min at 4 °C. After removing 40 μl of supernatant to save as input, 400 μl of the supernatant of each sample was mixed with 3 μg of antibody specifically recognizing PIWIL3 coupled with Protein A agarose, then incubated overnight at 4 °C with gentle rotation. After washing the agarose with wash buffer (50 mM Tris–HCl [pH 7.4], 0.005% NP-40, and 150 mM NaCl) four times, immunoprecipitation (IP) pellets were boiled at 95 °C for 5 min with 1 × protein loading buffer. Anti-Flag-HRP (Sigma, A8592) was used to detect the target proteins via western blot.

RNA immunoprecipitation of PIWI proteins in oocytes

Protein A Magnetic Beads (NEB, S1425S) were washed twice with RIP lysis buffer (50 mM Tris–HCl [pH 7.4], 150 mM NaCl, 1 mM EDTA, 0.5% NP-40, 0.5 mM DTT, and 0.1 U/μl RNase Inhibitor) containing 0.4 U/μl proteinase inhibitor cocktail (Sigma). For immunoprecipitation of PIWI protein for each sample, the washed beads derived from 2.5 μl of the original bead slurry were mixed with 1.5 μg antibody to a total volume of 10 μl and incubated on a rotator at room temperature for 0.5 h. Coupled beads were washed with RIP buffer before use. Oocyte samples were pooled for extraction for some species, namely guinea pig (six GV oocytes), goat (10 MII oocytes), pig (10 MII oocytes), and rat (nine MII oocytes). The pooled samples were lysed by adding 30 μl, 30 μl, 40 μl, or 22.5 μl cold RIP lysis buffer, respectively, in a DNA low-bind tube (Eppendorf, Germany) and incubating on ice for 10 min. Samples were then centrifuged at 16,200 × g for 5 min at 4 °C. Each 7.5 μl sample of oocyte lysate was transferred to a PCR tube that contained antibody-coupled beads for RIP assays. RNA was extracted with TRIzol from the remaining oocyte lysate in each tube after distribution for IP reactions. For IP reactions, after rotating for 5 h at 4 °C, beads were washed with RIP lysis buffer by rotating at room temperature for 5 min; this was repeated three times and followed by another wash with PBS. The beads were resuspended in 2 μl water and incubated at 75 °C for 5 min. After cooling on ice, the tubes were placed on a magnetic stand (Thermo Scientific, 12321D) for 10 s, and the supernatants were recovered for RNA extraction and sequencing library construction.

Small RNA library construction

Small RNA libraries were constructed from single oocytes as previously described [33]. Briefly, each single oocyte was lysed and incubated at 72 °C for 3 min, then cooled on ice. After 3′ adapter ligation, samples were incubated with 5 U of lambda exonuclease and 25 U of 5′ de-adenylates. Small RNAs were reverse-transcribed after 5′ adapter ligation. Two rounds of amplification were conducted to obtain the libraries, which were recovered on a 6% polyacrylamide gel.

mRNA library construction

mRNA library construction was conducted as previously described [33]. Briefly, a single oocyte was lysed and incubated at 72 °C for 3 min, then cooled on ice. Next, 1 μl of a 1/500,000–1/50,000 dilution of the ERCC RNA Spike-In Mix (Invitrogen, 4,456,740) was added to each sample. After reverse transcription and PCR pre-amplification, cDNA was purified and 3 ng used for a tagmentation reaction with Tn5 transposase. Amplified libraries were recovered and sequenced as described above.

Phylogenetic species tree

A phylogenetic tree was built for the 12 representative vertebrate species with TimeTree [54] and visualized with ggtree [55].

Sources of RNA sequences and genome assemblies

miRNA sequence annotations were downloaded from miRbase (Version 22) (http://www.mirbase.org/). mRNA, lncRNA, rRNA, snoRNA, and tRNA sequences were downloaded from Ensembl Genes (www.ensembl.org, version 100) and the Genomic tRNA Database (http://lowelab.ucsc.edu/GtRNAdb). The reference genome sequences, repeat annotations (RepeatMask files), and genomic annotations (Gene Transfer Format) were downloaded from the US National Center for Biotechnology Information (NCBI) RefSeq database. The accession or version numbers were as follows: GRCh38.p13 (human), Macaca_fascicularis_5.0 (monkey), OryCun2.0 (rabbit), Cavpor3.0 (guinea pig), GRCm38.p6 (mouse), Rnor_6.0 (rat), MesAur1.0 (golden hamster), CriGri-PICRH-1.0 (Chinese hamster), CanFam3.1 (dog), ARS1 (goat), Sscrofa11.1 (pig), and GRCz11 (zebrafish). The consensus TE sequences were downloaded from Repbase [56]. RepeatModeler2 was also used to predict and classify consensus TEs using the default parameters [57].

Oocyte small RNA-Seq

High-throughput RNA sequencing was performed on a NovaSeq 6000 (PE150). Because small RNAs typically range from 18–40 nt in length, we only used the R1 FASTQ files for subsequent analyses. Cutadapt was used to clip adaptor sequences and to filter out low-quality reads [58]. Reads that failed to match the adaptor sequences or that were shorter than 17 nt in length were discarded. Redundant sequences were collapsed as high quality reads (Gordon A, Hannon GJ: Fastx-toolkit. FastQ/A short-reads pre-processing tools, unpublished) and mapped to the reference sequences using bowtie [59]. The resulting reads were classified as known miRNAs, tRNA-derived small non-coding RNAs (tsRNAs), rRNA-derived small non-coding RNAs (rsRNAs), small snoRNAs, lncRNAs, and mRNAs, successively. Reads that could not be mapped to these known small RNAs were used in successive piRNA and endo-siRNA prediction. Sequences that were not annotated as any of the RNA categories described above were classified as “others.”

Oocyte mRNA-Seq

Trimmomatic [60] was used to clip adaptor sequences and to filter out low-quality reads. The high-quality reads were mapped to reference genomes using STAR [61]. Gene expression levels were calculated as transcripts per kilobase of exon model per million mapped reads (TPM) with StringTie [62].

piRNA sequence logos

The sequence logos for nucleotide positions 1–10 in piRNAs were calculated with weblogo [63].

Ping-pong signatures

The score for a distance of x nt was calculated as follows:

where Mi is the number of reads with piRNA 5′ ends located at position i and N i+x is the number of reads with piRNA 5′ ends located at position i + x (opposite strand). When x = 0, the 5′ ends of two piRNAs had no overlap. When x = 10, the 5′ ends of two piRNAs had a 10-nt overlap (ping-pong). For analyses including multi-mappers, reads were apportioned based on the number of times they could be aligned to the genome. In ping-pong analyses, overlaps at nt positions 1–9 and 11–20 were used as the background to calculate Z-scores.

Prediction of piRNAs and piRNA clusters

Small RNAs that could not be mapped to known annotations were used for piRNA cluster prediction using previously described methods [64]. To remove potential noise, we used stringent parameters in proTRAC (v2.4.2) to identify piRNA clusters in these 12 species (-1Tor10A 0.75 -1Tand10A 0.5). Due to differences in the length distribution of small RNA sequences in these species, we set the piRNA length limit to 17–23 nt (-pimin 17 -pimax 23) in rabbits and guinea pigs. In mice, rats, Chinese hamsters, and zebrafish, we set the piRNA length limit to 24–32 nt (-pimin 24 -pimax 32). In humans, monkeys, golden hamsters, dogs, goats, and pigs, we set the piRNA length limit to 17–32 nt (-pimin 17 -pimax 32). Other parameters in proTRAC were not changed, and the default parameters were used. Unknown reads of 17–32 nt in length located in these clusters were defined as piRNA candidates.

Prediction of endo-siRNAs and endo-siRNA clusters

Endo-siRNA clusters were defined as previously described [16]. Briefly, we used proTRAC (v2.4.2) to identify endo-siRNA clusters in mouse oocytes (-pimin 17 -pimax 24 -1Tor10A 0.25 -1Tand10A 0.1 -pdens 0.01 -clstrand 0.5). The resulting coordinates of endo-siRNA clusters are listed in Additional file 2: Table S18. Compared to piRNAs, endo-siRNAs have no preference for 10A, and the preference for 1U is less significant. Therefore, reads that could not be mapped to known small RNAs (miRNAs, tsRNAs, rsRNAs, and sn/snoRNAs) and were 17–24 nt in length were used to predict candidate endo-siRNA clusters, which were classified as those that met two criteria: (1) more than 75% of the reads were 21–23 nt in length and (2) at least four sequences in the cluster had a 2-nt overhang and the percentage of sequences having a 2-nt overhang in the cluster was > 10%. Reads 17–24 nt in length that could be mapped to these endo-siRNA clusters were predicted to be endo-siRNAs.

miRNA families

We compared miRNA expression levels in the 12 representative vertebrate species using miRNA families. miRNA family annotations were downloaded from TargetScan (Release 7.2) [65]. When small RNAs were mapped to the miRNA families, there were 0 mismatches allowed. Expression levels of miRNA families were qualified by size factor [66].

Synteny analysis

Synteny analysis was performed using liftOver in the UCSC toolkit (http://genome.ucsc.edu/cgi-bin/hgLiftOver). Human (hg38) piRNA clusters were aligned to monkey (macFas5), rabbit (oryCun2), guinea pig (cavPor3), mouse (mm10), rat (rn6), golden hamster (mesAur1), Chinese hamster (criGriChoV2), dog (canFam3), pig (susScr11), goat (oviAri4), and zebrafish (danRer11) genomes using the command “liftOverminMatch = 0.1”. piRNA-producing clusters in the human genome that had been successfully identified in another genome were considered syntenic. Loci producing piRNAs from a syntenic location in another species were considered to be evolutionarily conserved.

PhyloP scores calculation

ClustalW2 was used for multiple sequences alignment [67]. We used the likelihood ratio test (LRT) method and CONACC mode from phyloP to compute conservation scores for each site in the alignment [68]. We then divided the site equally into 100 bins according to their coordinates and calculated average phyloP scores for each bin. The values represent -log p-values under a null hypothesis of neutral evolution. The sites predicted to be conserved are assigned positive scores while sites predicted to undergo accelerated evolution are assigned negative scores.

Heatmap and principal component analyses

A heatmap was plotted using the “pheatmap” package in R [69]. PCA was performed using the prcomp function in R and visualized using the “ggplot2” package [70].

Small-RNA spike-ins and normalization

Three unmethylated and three methylated small-RNA spike-ins were synthesized by Integrated DNA Technologies (IDT) and pooled. The sequence of each spike-in can be found in our previous study [16]. To quantify the expression levels of small RNAs in the RNA IP assays and NaIO4 oxidation data, we normalized the read counts of small RNAs to those of the spike-ins.

Mapping piRNA and RNA-Seq data to TEs

To measure the expression levels of TE-related piRNAs, piRNAs were mapped to consensus TEs. For analyses including multi-mappers, reads were classified based on the number of times they could be aligned to the TEs. We quantified the TE-related piRNA expression levels in reads of exon model per million mapped fragments (RPM). Oocyte RNA-Seq data were mapped to the consensus TEs using bowtie2 [71]. TE mRNA expression levels were quantified in fragments per kilobase of exon model per million mapped fragments (FPKM).

TE genomic features

BED files containing repetitive elements predicted by RepeatMasker were downloaded from NCBI RefSeq. The genomic features of average TE divergence and TE copy numbers were calculated based on the BED files. Pearson’s correlation coefficients for piRNA, TE mRNA abundance, and genomic features were calculated using the “corrplot” package in R [72].