key: cord-0305917-4levwk57
authors: Ishino, Kyoko; Hasuwa, Hidetoshi; Yoshimura, Jun; Iwasaki, Yuka W.; Nishihara, Hidenori; Seki, Naomi M.; Hirano, Takamasa; Tsuchiya, Marie; Ishizaki, Hinako; Masuda, Harumi; Kuramoto, Tae; Saito, Kuniaki; Sakakibara, Yasubumi; Toyoda, Atsushi; Itoh, Takehiko; Siomi, Mikiko C.; Morishita, Shinichi; Siomi, Haruhiko
title: Hamster PIWI proteins bind to piRNAs with stage-specific size variations during oocyte maturation
date: 2020-12-02
journal: bioRxiv
DOI: 10.1101/2020.12.01.407411
sha: d0fc0cb60aff7463f5cc24221c28752e58549c55
doc_id: 305917
cord_uid: 4levwk57

In animal gonads, transposable elements (TEs) are actively repressed to preserve genome integrity through the Piwi-interacting RNA (piRNA) pathway. In mice, piRNAs are most abundantly expressed in male germ cells, and form effector complexes with three distinct PIWI proteins. The depletion of individual Piwi genes causes male-specific sterility owing to severe defects in spermatogenesis with no discernible phenotype in female mice. Unlike mice, most other mammals have four PIWI genes, some of which are expressed in the ovary. Here, purification of PIWI complexes from oocytes of the golden hamster revealed that the size of the piRNAs loaded onto PIWIL1 changed during oocyte maturation. In contrast, PIWIL3, an ovary-specific PIWI in most mammals, associates with short piRNAs only in metaphase II oocytes, which coincides with intense phosphorylation of the protein. An improved high-quality genome assembly and annotation revealed that PIWIL1- and PIWIL3-associated piRNAs appear to share the 5′- ends of common piRNA precursors and are mostly derived from unannotated sequences with a diminished contribution from TE-derived sequences, most of which correspond to endogenous retroviruses (ERVs). Although binding sites for the transcription factor A-Myb are identified in the transcription start site regions of the testis piRNA clusters, the piRNA clusters in the ovary show no well-defined binding motifs in their upstream regions. These results show that hamster piRNA clusters are transcribed by different transcriptional factors in the ovary and testis, resulting in the generation of sex-specific piRNAs. Our findings show the complex and dynamic nature of biogenesis of piRNAs in hamster oocytes, and together with the new genome sequence generated, serve as the foundation for developing useful models to study the piRNA pathway in mammalian oocytes. Highlights - The size of PIWIL1-associated piRNAs changes during oocyte maturation - Phosphorylation of PIWIL3 in MII oocytes coincides with its association with small 19-nt piRNAs - Improved high-quality genome assembly and annotation identifies young endogenous retroviruses as major targets of piRNAs in hamster oocytes - PIWIL1- and PIWIL3-associated piRNAs share the 5′-ends of the common piRNA precursors in oocytes

Unlike mice, most other mammals have four PIWI genes, some of which are expressed in the ovary. Here, purification of PIWI complexes from oocytes of the golden hamster revealed that the size of the piRNAs loaded onto PIWIL1 changed during oocyte maturation. In contrast, PIWIL3, an ovary-specific PIWI in most mammals, associates with short piRNAs only in metaphase II oocytes, which coincides with intense phosphorylation of the protein. An improved high-quality genome assembly and annotation revealed that PIWIL1-and PIWIL3-associated piRNAs appear to share the 5´ends of common piRNA precursors and are mostly derived from unannotated sequences with a diminished contribution from TE-derived sequences, most of which correspond to endogenous retroviruses (ERVs). Although binding sites for the transcription factor A-Myb are identified in the transcription start site regions of the testis piRNA clusters, the piRNA clusters in the ovary show no well-defined binding motifs in their upstream regions. These results show that hamster piRNA clusters are transcribed by different transcriptional factors in the ovary and testis, resulting in the generation of sex-specific piRNAs. Our findings show the complex and dynamic nature of biogenesis of piRNAs in hamster oocytes, and together with the new genome sequence generated, serve as the foundation for developing useful models to study the piRNA pathway in mammalian oocytes.

Transposition of mobile DNA elements can cause severe damage by disrupting either the structural or regulatory regions on the host genome (Chuong et al., 2017; Han and Boeke, 2005; Hancks and Kazazian, 2016) . To avoid such detrimental effects, many animals have a conserved adaptive immune system known as the piRNA pathway in gonads (Iwasaki et al., 2015; Ozata et al., 2019; Pillai and Chuma, 2012) . piRNAs form effector complexes with PIWI proteins, a germline-specific class of Argonaute proteins, to guide recognition through complementary base-pairing, leading to silencing their target transposable elements (TEs) mainly in two ways: post-transcriptional silencing by PIWI-mediated cleavage of target transcripts in the cytoplasm and transcriptional silencing by PIWImediated chromatin modifications on target loci. Mutations in genes involved in the pathway can lead to sterility.

Although the mechanisms for generating piRNAs appear to largely differ among animals, the defined characteristics of piRNAs include a predominant length of 26-32 nucleotides (nt), a strong bias for uracil (U) at the 5´-ends (1U-bias), 2-Omethylation of the 3´-ends, and clustering of reads to distinct genomic locations (Iwasaki et al., 2015; Ozata et al., 2019) . The characterization of the piRNA populations in Drosophila and mouse has led to two models for the biogenesis mechanisms: the pingpong cycle (Brennecke et al., 2007; Gunawardane et al., 2007) and phased Mohn et al., 2015) , which are intimately connected. Long single-stranded precursors, often more than 10 kb in size, are derived from discrete genomic loci that are now referred to as piRNA clusters or piRNA genes (Aravin et al., 2006; Brennecke et al., 2007; Girard et al., 2006; Lau et al., 2006; Vagin et al., 2006) . Two major clusters exist to serve as genomic piRNA source loci: intergenic and genic clusters. Intergenic piRNA clusters are often located in the heterochromatic regions and comprise various types of TEs that tend to be young and potentially active, suggesting that intergenic piRNA clusters provide the host with acquired and heritable memory systems to repress TEs (Brennecke et al., 2007; Khurana et al., 2011) . Genic piRNAs mainly arise from 3′ untranslated regions (UTRs) of the protein-coding genes (Robine et al., 2009; Saito et al., 2009 ). The function of the genic piRNAs is not well understood, but some genic piRNAs 4 show significant complementarity to protein-coding genes (Saito et al., 2009; Gonzalez et al., 2015) . piRNA precursors are cleaved by either endonuclease Zucchini/mitoPLD or by the Slicer activity of PIWIs with pre-existing piRNAs to produce 5´ monophosphorylated piRNA intermediates that are loaded onto PIWI proteins (Brennecke et al., 2007; Gainetdinov et al., 2018; Gunawardane et al., 2007; Han et al., 2015; Ipsaro et al., 2012; Mohn et al., 2015; Nishimasu et al., 2012) . PIWI-piRNA complexes then initiate piRNA production that is formed by an amplification loop termed the ping-pong cycle in which reciprocal cleavage of TE and cluster transcripts by PIWI proteins generates new piRNA 5´-ends and amplifies piRNA populations while destroying TE mRNAs in the cytoplasm.

The ping-pong cycle produces two classes of piRNAs overlapped by precisely 10 nt at their 5´-ends: one class shows a strong preference for U at their 5´-ends (1U) and the second class shows a preference for adenine at nucleotide 10 (10A). The ping-pong cycle is then accompanied by the phased production of piRNAs downstream of the cleavage site, which further creates a sequence diversity of piRNAs. The 3´-ends of piRNAs are determined either by direct cleavage of Zucchini/mitoPLD (mouse Zucchini homolog) or PIWIs or by trimming piRNA intermediates by resecting enzymes (Nibbler in Drosophila, Trimmer in silkworm, and PNLDC1 in mouse) (Ding et al., 2017; Hayashi et al., 2016; Izumi et al., 2016; Nishida et al., 2018; Nishimura et al., 2018) . piRNA biogenesis is then concluded with 2-O-methylation of the 3´-ends by HENMT1 methylase, which has been hypothesized to stabilize mature piRNAs (Gainetdinov et al., 2018; Horwich et al., 2007;  Kirino and Mourelatos, 2007; Saito et al., 2007) . The extent of 3´-end cleaving/trimming and consequently piRNA length is determined by the footprint of the PIWI protein, possibly explaining the different size profiles of piRNAs associated with distinct PIWI proteins. Structural studies have shown that the 5´-and 3´-ends of the guide small RNAs, including piRNAs, are recognized by the MID-PIWI and PAZ domains of Argonaute/PIWI proteins, respectively (Matsumoto et al., 2016; Wang et al., 2008; Yamaguchi et al., 2020) .

Mammalian PIWI-piRNA pathways have been studied mainly in mice (Pillai and Chuma, 2012) . Similar to Drosophila, mice express three PIWI proteins, MIWI (PIWIL1), MILI (PIWIL2), and MIWI2 (PIWIL4) in the testis, with varying spatiotemporal expression patterns during spermatogenesis. The non-redundant role of

Piwi genes causes male-specific sterility owing to severe defects in sperm formation (Aravin et al., 2007; Aravin et al., 2008; Carmell et al., 2007; Deng and Lin, 2002; Ernst et al., 2017; Kuramochi-Miyagawa et al., 2008; Pillai and Chuma, 2012; Thomson and Lin, 2009) . Each PIWI protein associates with distinct piRNA populations; fetal prepachytene piRNAs and pachytene piRNAs. Prepachytene piRNAs are formed from TE-and other repeat-derived sequences. In contrast, pachytene piRNAs have a higher proportion of unannotated sequences with the diminished contribution of TE sequences, though they still function in TE silencing by guiding MILI and MIWI to cleave TE transcripts (De Fazio et al., 2011; Reuter et al., 2011) . A specific transcriptional factor, A-MYB, binds the promoter regions of pachytene piRNA clusters as well as core piRNA biogenesis factors, including MIWI, and initiates their transcription (Li et al., 2013) .

Some pachytene piRNA clusters are divergently transcribed from bidirectional A-Mybbinding promoters (Li et al., 2013) .

Although PIWI genes in fly and zebrafish are expressed in both testes and ovaries, mouse Piwi genes are expressed only weakly, if not at all, in the ovary, and depletion of these Piwi genes does not affect the female germline. Thus, these findings led to the assumption that the piRNA pathway does not play a role in mammalian oogenesis. Unlike mice, many other mammals possess four distinct PIWI genes (Piwil1-4), suggesting that piRNA-mediated silencing may differ between mice and mammals with four PIWI genes (Hirano et al., 2014; Sasaki et al., 2003) . However, except for mice, little is known about mammalian piRNA pathways, particularly their roles in ovaries.

Although recent studies have reported the presence of piRNA-like molecules in mammalian female germ cells, including humans (Kabayama et al., 2017; Roovers et al., 2015; Williams et al., 2015; Yang et al., 2019 ) , it is not yet known whether they play a role in the ovary because of the difficulty of technical and/or ethical application of genetic analysis. Although mice and rats belong to the Muridae family of rodents, both of which lack Piwil3, the golden Syrian hamster (golden hamster, Mesocricetus auratus) belongs to the Cricetidae family and has four PIWI genes. Golden hamsters have been used as an experimental rodent model for studying human diseases, particularly for cancer development and many different infectious diseases including COVID-19, because they display many features that resemble the physiology and pharmacological responses of humans (Hirose and Ogura, 2019; Sia et al., 2020) . In addition, methods for manipulating gene expression, including the CRISPR/Cas9 system, have been recently employed in golden hamsters (Fan et al., 2014; Hirose et al., 2020) . Herein, we analyzed PIWIassociated piRNAs in oocytes and early embryos of golden hamsters, in the hope of applying genetic analysis to the piRNA pathway in the ovary. Our analyses revealed that the size of PIWIL1-associated piRNAs changes during oocyte maturation and that PIWIL3 binds short piRNAs only at the metaphase II (MII) stage of the oocyte, which coincides with phosphorylation of the protein. With an improved high-quality genome assembly and annotation of golden hamster, we showed that PIWIL1-and PIWIL3associated piRNAs appear to share their 5´-ends. Their contents are similar to those observed with pachytene piRNAs in the mouse testis, but their targets in oocytes are mostly endogenous retroviruses. We further identified ovarian piRNA clusters, and motif search for the transcription start site regions of the piRNA clusters revealed no distinct binding motifs in their upstream regions, although A-Myb-binding motifs were enriched in the upstream regions of the testis piRNA clusters. Our study provides a basis for understanding the roles of the piRNA pathway in mammalian oocytes.

To examine the expression of PIWI genes in golden hamster gonads, we performed RNAsequencing (RNA-seq) analysis using hamster ovary, oocyte, and testis samples and analyzed the expression level of PIWI genes by calculating transcripts per kilobase million mapped (TPM) ( Figure 1A ). The estimated expression levels of Piwil1 and Piwil2 were relatively high in the testis (TPM = 24.24 and TPM = 14.07, respectively). Piwil1 was also expressed in the ovary (TPM = 4.63), while PiwiL2 is not detected in the whole ovary and appeared to be only weakly expressed in the oocyte. Piwil4 appeared to be expressed only in the testis, consistent with previous transcriptome analysis in human oocytes (Yang et al., 2019) . Interestingly, Piwil3 was highly expressed in the oocyte (TPM = 14.60). In sharp contrast, the expression of Piwil3 could not be detected in the testis. These results are consistent with previous analyses of bovine and human PIWIL3 (Roovers et al., 2015; Yang et al., 2019) . To analyze the expression levels of known PIWI-piRNA pathway factors other than PIWI genes, we also calculated the TPM values of the predicted homologous genes using RNA-seq data ( Figure S1A ). Most of the known PIWI-piRNA pathway factor homologs, including Mov10l1, Mael, MVH (mouse Vasa homolog), MitoPLD, Gtsf1, Henmt1, and Tudor domain families (Iwasaki et al., 2015) , were expressed in both testes and oocytes, suggesting that both testes and oocytes of hamsters have similar biogenesis pathways to produce piRNAs.

To confirm the expression of PIWIs in the ovary, we isolated their cDNAs from the ovary. Open reading frames (ORFs) of sequenced Piwil1 and Piwil2 cDNAs correspond well with the respective annotated gene products deposited in the Broad Institute database (MesAur1.0, Broad Institute data) ( Figure S1B ). However, to our surprise, during the cloning of the Piwil3 cDNA, we found that the large extension of the amino-terminal portion of the ORF contains 14 repeats of nucleotide sequences encoding the amino acid sequences QLQSPGAGPPRSGA ( Figure 1B ). To further confirm the expression of PIWIs in the ovary at the protein level, we produced specific monoclonal antibodies against PIWIL3 ( Figure S1C ). We also found that a monoclonal antibody that recognizes marmoset PIWIL1 (Hirano et al., 2014) cross-reacts with hamster PIWIL1 specifically among hamster PIWI proteins ( Figure S1C ). Thus, we focused our analysis on PIWIL1 and PIWIL3. Immunostaining with the antibodies produced showed that both PIWIL1 and PIWIL3 were expressed in the cytoplasm of growing oocytes in the ovary ( Figure 1C ). Western blotting with anti-PIWIL1 antibody showed a discrete band at 90 kDa in the ovary, metaphase II (MII) oocytes, and 2-cell embryos ( Figure 1D ). Western blotting with anti-PIWIL3 antibody revealed a discrete band at 130 kDa in the ovary and 2-cell embryos, which is consistent with the calculated molecular weight of the protein with the large amino-terminal extension ( Figure 1B and 1D) . Intriguingly, the size of the protein largely shifted by approximately 40 kDa in MII oocytes ( Figure 1D ).

This large increase of the PIWIL3 protein in size, together with the broad and fuzzy nature of the band on the gel ( Figure 1D ), prompted us to test whether this size shift could be due to modification of the protein with phosphorylation. Treatment of MII oocytes with calf intestinal alkaline phosphatase (CIP), an enzyme known to dephosphorylate proteins , caused a discrete band to migrate to the estimated molecular weight of 130 kDa, demonstrating that PIWIL3 is heavily phosphorylated in MII oocytes ( Figure   1E ). These results show that PIWIL1 and PIWIL3 are expressed in growing oocytes in the ovary as well as in early embryos and that PIWIL3 is modified with phosphorylation specifically at the MII stage of oocytes. Since the mouse genome lacks Piwil3 and thus the characterization of PIWIL3 protein has been delayed, our findings indicate that Piwil3 may have specific functions in female gonads.

To isolate PIWIL1-and PIWIL3-associated small RNAs from oocytes, we immunopurified the associated complexes from MII oocytes with specific monoclonal antibodies produced. We then isolated RNAs, 32 P-labeled them, and analyzed them using a denaturing polyacrylamide gel ( Figure 2A ). Intriguingly, PIWIL1 was associated with two populations of small RNAs: one with 29-30 nt and the other with 22-23 nt in MII oocytes. We then immunopurified PIWIL1-associated complexes from whole ovaries (including growing oocytes) and 2-cell embryos. The sizes of PIWIL1-associated piRNAs in whole ovaries and 2-cell embryos were 29-30 nt and 22-23 nt, respectively ( Figure 2B ). The resistance of PIWIL1-associated piRNAs in MII oocytes to periodate oxidation (NaIO4) and β-elimination reactions show that they are modified at the 3´ terminal nucleotide with a 2′-O-methyl marker ( Figure S2A ). We also isolated PIWIL1associated small RNAs from whole testes and found that small RNAs (29-30 nt) were loaded onto PIWIL1 in the testis ( Figure S2B ). These results show that piRNAs loaded onto PIWIL1 change their sizes during oocyte maturation, from the size equivalent to that observed in the testis to a mixture of long and short populations, and short piRNAs (22-23-nt). To our knowledge, this is the first study to show that the size of piRNAs loaded onto a distinct PIWI protein changes during germline development.

In sharp contrast, PIWIL3 was found to bind to a class of small RNAs with 19 nt ( Figure   2A ), which is consistent with the recent finding that human PIWIL3 associates with small RNAs (~20 nt) in human oocytes (Yang et al., 2019) . 19-20 nt RNAs in oocytes almost completely disappeared after β-elimination reactions, indicating that PIWIL3-associated piRNAs lack 2′-O-methylation at their 3′ terminal ( Figure S2C ). We failed to detect small RNAs associated with PIWIL3 in whole ovaries and 2-cell embryos ( Figure 2C ). With the caveat that this could be because of technical reasons for immunoprecipitation with the antibodies and/or the buffer conditions used, our findings suggest that PIWIL3 may bind piRNAs with 19 nt exclusively in MII oocytes but be freed from piRNAs as 'empty' PIWIL3 at the early stages of oocyte maturation and in early embryos. PIWIL3 is heavily phosphorylated only in MII oocytes, raising the possibility that phosphorylation of the protein may be required for the association with a class of short piRNAs. To test this, we performed immunoprecipitation with an anti-PIWIL3 antibody using MII oocyte lysate that had been treated with CIP and examined whether the CIP treatment affected the association of piRNAs with PIWIL3. Indeed, PIWIL3 treated with CIP was free from piRNAs ( Figure 2D ). This shows that phosphorylation of the protein is required for PIWIL3 probably either to get loaded with piRNAs or hold them or both.

Although the draft genome of the golden hamster has been sequenced (the MesAur1.0 genome), we soon came to realize that we needed much more accurate genome sequence data to further characterize these PIWI-associated piRNAs on the genome mainly because the MesAur1.0 genome sequence contains a large number of gaps (N) (17.58% of the genome; 420 Mb of the 2.4 Gb) and remains fragmented. Because TEs and other repeats in the genome are the main targets of piRNAs in many animals, the lack of accurate sequences of TEs and other repeats is a serious problem when mapping piRNAs on the genome. Accurate detection of TEs requires both full collection/classification of TE consensus sequences and high-quality genome assembly. Thus, we re-sequenced the golden hamster genome.

Details of the hamster genome assembly are shown in the Methods section. The final genome assembly is summarized in Table 1 . The nucleotide difference between the DNA Zoo MesAur1.0_HiC assembly and our assembly was 0.140%. We assessed the completeness of the genome assemblies using the BUSCO tool (Waterhouse et al., 2018) and found that our golden hamster genome assembly included 3,991 complete genes (97.2%) and 37 fragmented genes among 4,104 single-copy genes. Our new golden hamster genome allowed us to resolve several issues that had remained ambiguous. For example, although putative ancestral karyotypes of rodents in the Muridae and Cricetidae families have been partly reconstructed by traditional chromosome painting (Romanenko et al., 2012) , we compared our nearly complete genome with the mouse (Mus musculus) and rat (Rattus norvegicus) reference genomes (Methods) and identified conserved synteny blocks between the golden hamster, mouse, and rat genomes ( Figure 3A ). We inferred ancestral karyotypes by integrating synteny blocks shared between two or three species according to maximum parsimony, to minimize the amount of chromosomal rearrangement (Methods) and obtained a high-resolution ancestral karyotype of Muridae using the golden hamster genome as the outgroup as well as a precise ancestral Cricetidae karyotype ( Figure 3B ). Although our ancestral Cricetidae and Muridae karyotypes were consistent with most previous inference (Romanenko et al, 2012) , we resolved some problems: for example, it was unclear whether mouse chromosomes 3 and 4 possessed synteny blocks from a common ancestral karyotype, and our analysis demonstrated the existence of a proto-chromosome in the ancestral Cricetidae karyotype (brown blocks in Figure 3B ). We also confirmed previous speculation that mouse chromosomes 5 and 16 obtained blocks from a common Muridae proto-chromosome (light orange in Figure 3B ), and that chromosomes 10, 15, and 17 obtained blocks from a common Cricetidae protochromosome (maroon blocks in Figure 3B ). Finally, we identified two groups of mouse chromosomes (6, 8) and (8, 13) having large blocks from common Muridae protochromosomes.

Using our genome assembly, we compared the hamster genomic locus containing the Piwil3 gene with syntenic regions of the mouse and rat genomes. The Piwil3 gene is flanked by Wscd2 and Sgsm1 in the hamster genome ( Figure S3A ). The order of the two genes is conserved in the syntenic regions of the mouse and rat genomes.

We then extracted the genomic regions between these genes from our hamster genome and the reference genomes of mouse (mm10) and rat (rn6), compared the syntenic regions using dot plots ( Figure S3B ), and observed the absence of the Piwil3 gene. The validity of the hamster genomic region with Piwil3 was confirmed by checking each base in the region was covered by an ample number of long reads ( Figure S3C ). In addition, we performed a similarity search with blastn and ssearch36 using the protein-coding sequence (CDS) of hamster Piwil3 as a query and found no hits in the corresponding regions of the mouse and rat genome. The PIWIL3 gene is conserved in most mammals, including humans, suggesting that Piwil3 might have been deleted after speciation in mice and rats.

With the new golden hamster genome sequence generated, we also conducted a de novo repeat characterization and identified 177 consensus sequences of repetitive elements at the subfamily level, including 3 SINEs, 12 LINEs, 156 long terminal repeat (LTR) retrotransposons, and 2 DNA transposons. RepeatMasker analysis using our custom repeat library (RepeatMasker rodent library with newly identified 177 consensus sequences) revealed that SINEs, LINEs, LTR retrotransposons, and DNA transposons occupy 9.2%, 16.9%, 12.1%, and 1.3% of the hamster genome, respectively. The contents of TEs are equivalent to those found in mice and rats, but the fractions of SINE and LINE are, respectively, higher and slightly lower in the hamster than those observed in mice and rats ( Table 2 ). The contents of the MesAur1.0 genome (MesAur1.0_HiC) were similarly analyzed. In contrast to our assembly, a much lower proportion of young TEs were detected even using our custom repeat library ( Figure S4A , Table S4 ). This is mostly because of the high frequency of gaps (Ns) in the MesAur1.0_HiC assembly (Table 1) , which resulted in apparently less similarity between TEs and their consensus sequences.

The custom library substantially improved the detection ability of recently active TEs, as represented by a higher proportion of young elements, for example, those with low (< 5.0) Kimura two-parameter (K2P) divergence from the consensus ( Figure   4A ). Hamster-specific subfamilies of B1 SINE were recently active, and 36,000 copies of young full-length B1s are present in the assembly. We identified three groups of the LINE-1 (L1) family, two of which were recently highly active (L1-4_MAu and L1-5_MAu), and the genome harbors at least 1,000 young full-length copies ( Figure 4B ).

The most active LTR retrotransposons belong to the ERV2/ERVK superfamily, including IAP (for ERV classification, see Kojima 2018; Vargiu et al., 2016) . We identified 20 families of ERV2 that contain an internal portion between LTRs, and 11 of them possess a clear reverse transcriptase domain. There are over 13,000 LTR sequences and 1,600 internal portions of the recently active elements in the genome. Among them, we found three recently expanded families: IAP1E_MAu, ERV2-5_MAu, and ERV2-7_MAu, which accounted for 29.9%, 14.7%, and 26.8% of the very young (<2.0 K2P divergence) LTR retrotransposons, respectively ( Figure 4C ). It is likely that not only IAP but also other ERV2/ERVK elements are currently active in the hamster genome. Together, these results show that our effort to re-sequence the golden hamster genome significantly improved annotations, especially for recently active TEs, which are potential piRNA targets.

To characterize piRNAs loaded onto PIWIL1 or PIWIL3 in oocytes, we performed small RNA sequencing using piRNA samples immunopurified with an anti-PIWIL1 or an anti-PIWIL3 antibody from oocytes. For PIWIL1-associated piRNAs, we also sequenced piRNA samples immunopurified with an anti-PIWIL1 antibody from testes. Replicates were highly correlated with each other (R 2 > 0.9) (data not shown); therefore, we merged the reads. First, we performed a length distribution analysis of the obtained reads ( Figure   5A -E). We observed peak signals at the size range of 29-30 nt in the testes and ovaries, 29 nt and 23 nt in MII oocytes, and 23 nt in 2-cell embryos of PIWIL1-associated piRNAs and 19 nt in MII oocytes of PIWIL3-associated piRNAs, confirming that the length of PIWIL1-and PIWIL3-associated piRNAs sequenced is consistent with that observed on the gels. Then we divided into two groups using the length of PIWIL1-associated piRNAs in MII oocyte; Oocyte long piRNAs (L1 OoL-piRNAs) and Oocyte short piRNAs (L1 OoS-piRNAs) for further analysis. Sequencing also revealed that the piRNA populations in all samples showed a strong 1U bias, a conserved characteristic for piRNAs ( Figure   5A -E).

We then mapped piRNAs to the new hamster genome (hamster.sequel.draft-20200302.arrow.fasta). Among the PIWI-associated piRNA reads, 50.0~67.4% of the reads were uniquely mapped to the genome and 6.2~43.3% were mapped multiple times ( Figure 6A , upper panel). Then, we annotated each piRNA read to analyze the genomic region from which the piRNA reads were generated. Approximately 79.05~89.84% of the reads mapped to unannotated intergenic regions, and only 7.07~13.19% originated from TEs ( Figure 6A , lower panel). This profile is similar to that of pachytene piRNAs associated with mouse MIWI, in which ∼70% is derived from intergenic regions and ∼24% from TEs (Reuter et al., 2011) .

The characteristics of PIWIL1-piRNAs were virtually indistinguishable between the testis and ovary. To determine whether this was due to the common piRNA sequences, we calculated reads per million mapped reads (RPM) for each piRNA sequence and compared the correlation between testis and ovary samples. The RPM of piRNAs greatly varied between the testis and ovary samples (R 2 = 2.08e-06). In addition, when we checked the top ten piRNA sequences with the most abundant read counts, none of the sequences were common between the testis and ovary (data not shown). These results show that testis and ovary PIWIL1-piRNAs possess distinct sequences ( Figure   6B ).

We found that most piRNAs in hamster female gonads were derived from LTR retrotransposons, including endogenous retroviruses (ERVs), compared to PIWIL1-piRNAs in the hamster testis, which corresponded to both LINEs and LTR retrotransposons ( Figure S5 ). This observation is consistent with the fact that the main targets of the piRNA pathway in mouse testes are L1 and IAP elements. In the mouse genome, L1 is the most active TE family, as represented by an increasing accumulation of their young copies ( Figure S4B ). There are also young LTR retrotransposons in mice, among which IAP and MERVL (ERV3/ERVL) families are highly active, while ERV2/ERVK families, except IAPs, showed much lower proportions among the young elements ( Figure S4C ). Interestingly, most LTR retrotransposons from which piRNAs in hamster female gonads were derived were evolutionally younger judged by the K2P divergence ( Figure 4C ). This suggests that hamster oocyte piRNAs were generated from the higher activity of recent transposition. In contrast, testis PIWIL1-piRNAs target both LINEs and LTR retrotransposons. Together with the diversity in PIWI-associated piRNA sequences of the oocyte and testis ( Figure 6B ), this supports the notion that PIWI-piRNA target TEs are different for male and female gonads. Figure 2 and 5 show that piRNAs loaded onto PIWIL1 in MII oocytes consist of two distinct populations, short piRNAs (L1 OoS-piRNAs) and long piRNAs (L1 OoL-piRNAs). Shorter piRNAs (18-20 nt) were loaded onto PIWIL3 in MII oocytes. These findings prompted us to test whether these piRNAs are derived from the same genomic loci and whether L1 OoS-piRNAs and/or PIWIL3-associated piRNAs are processed products of L1 OoL-piRNAs by cleaving and/or trimming their 3´-ends. To this end, we asked whether the genomic positions of the 5´-or 3´-ends of the L1 OoL-piRNAs differ from those of L1 OoS-piRNAs or PIWIL3-piRNAs and calculated the probabilities of the 5´-or 3´-ends of L1 OoL-piRNAs coinciding with the 5´-or 3´-ends of L1 OoS-piRNAs or PIWIL3-piRNAs (Gainetdinov et al., 2018) . In MII oocytes, L1 OoL-piRNAs were much more likely to share 5´-ends with L1 OoS-piRNAs and PIWIL3-piRNAs (Oo PIWIL3) than with 3′-ends (0.65 for 5′-ends versus 0.03 for 3′-ends and 0.54 for 5′-ends versus 0.02 for 3′-ends, respectively; Figure 7A ). These results suggest that approximately half of these three piRNA populations share the same 5′ cleaved piRNA precursors, thereby sharing the same piRNA clusters as piRNA sources. Thus the size differences among these piRNAs account for their 3´-end formation of these populations that may be determined by the footprint of piRNA-binding PIWI proteins, probably because of either the structures of partner PIWI proteins (in particular, the distance between MID-PIWI and PAZ domains) or conformational changes in the partner proteins (see Discussion).

The finding that piRNA populations in oocytes share the same 5′ cleaved piRNA precursors suggests little ping-pong amplification among them. We calculated the ping-pong signature of each piRNA and found that they had little or no ping-pong signatures ( Figure 7B ). This result indicates that most of the piRNAs expressed in oocytes are not involved in the ping-pong cycle. Recently, a novel class of small RNAs with 21-23 nt, termed short PIWI-interacting RNAs (spiRNAs), was identified in mouse oocytes (Kabayama et al., 2017) . They are produced by the ping-pong cycle driven by the MILI slicer activity, and thus, small piRNAs found in hamster oocytes are distinct from spiRNAs.

We found that more than half of PIWIL1-and PIWIL3-associated piRNAs are likely to share 5´-ends of precursor transcripts, suggesting that the primary source of piRNAs is also shared among them. To test the possibility that PIWIL1-and PIWIL3-associated piRNAs share piRNA clusters from which they are derived, we next focused on the identification of piRNA clusters in hamster testes, ovaries, MII oocytes, and 2-cell embryos. We identified piRNA clusters using proTRAC (version 2.4.3) under the following conditions as described (Yang et al., 2019) with some modifications: (1) more than 75% of the reads were appropriate to the length of each piRNA; (2) more than 75% of the reads exhibited the 1 U or 10 A preference; (3) more than 75% of reads were derived from the main strand; and (4) -pimin option with 21, 28, and 18 for oocyte PIWIL1-piRNAs and PIWIL3-piRNAs, respectively. We identified 101, 89, 74, 55, and 61 piRNA clusters in testis PIWIL1-piRNAs, ovary PIWIL1-piRNAs, MII oocyte PIWIL1-piRNAs, 2-cell PIWIL1-piRNAs, and oocyte PIWIL3-piRNAs, respectively. In the testis, both unidirectional clusters, in which piRNAs map to only one strand, and bidirectional clusters, in which the polarity of piRNA production switches between plus and minus strands, were identified, although unidirectional clusters were dominant ( Figure 8A and Figure S7A ). We found that the majority of the piRNA clusters identified in female gonads were unidirectional ( Figure 8B and Figure S6A ).

Interestingly, approximately 90% of the piRNA clusters corresponding to PIWIL1-associated piRNAs in testis, ovary, and MII oocytes were sex-specific, and only eight piRNA clusters were identified commonly in both male and female gonads ( Figure   8C ). These results support the idea that transcription of the primary source of piRNA (piRNA clusters) plays a key role in the production of sex-specific piRNAs ( Figure 6B and 8C). We next examined overlaps among ovaries, MII oocytes, and 2-cell embryo PIWIL1 piRNA clusters and found that they shared 31 clusters ( Figure 8D ). PIWIL1-piRNAs in 2-cell embryos shared approximately 85% of the source clusters with those in MII oocytes, while they shared approximately 60% of the source clusters with those in the ovary, which suggests that the genomic regions where piRNAs are derived are altered during oocyte maturation. PIWIL3-associated piRNAs shared 77% of the clusters with PIWIL1-associated piRNAs in MII oocytes. These results are consistent with the findings that PIWIL1-and PIWIL3-associated piRNA populations are likely to share the same 5′ cleaved piRNA precursors ( Figure 7A ). Finally, we found that piRNAs loaded onto PIWIL3 in MII oocytes shared a large number of piRNA clusters with piRNAs loaded onto PIWIL1 in 2-cell embryos, suggesting that they may have common targets ( Figure   8E ).

Given that the transcription of piRNA clusters greatly differs in male and female gonads, we analyzed the motif sites surrounding the potential transcription start sites of the testis and ovary piRNA clusters. We first extracted sequences surrounding the transcriptional start sites of unidirectional and bidirectional piRNA clusters, predicted from small RNA-seq and RNA-seq mapping data. We used these sequences and performed motif searches using MEME v.5.1.0, which can discover motifs enriched in the given dataset. The results of the MEME analysis show that the A-Myb binding site is significantly represented in the bidirectional piRNA clusters in the testis transcriptional 

In this study, we generated a new golden hamster genome, which revealed the presence of young and possibly transposition-competent TEs in the genome. This also allowed us to characterize piRNAs in golden hamster oocytes. Intriguingly, the size of PIWIL1-associated piRNAs changes during oocyte maturation. In sharp contrast, PIWIL3 binds to piRNAs only in MII oocytes, and the size of loaded piRNAs is shorter (19 nt).

The change in the size of PIWIL1-associated piRNAs during oocyte maturation may necessitate unloading of the 3´-ends of the long piRNAs from PIWIL1 to shorten the long PIWIL1-associated piRNAs to short 22-23 nt either by trimming or by direct cleavage to determine their new 3′-ends. Alternatively, PIWIL1-associated short piRNAs may not be the processed products of loaded long piRNAs, but they may be produced by a mechanism similar to that produced by long piRNAs and then loaded onto either newly produced PIWIL1 or PIWIL1, which has unloaded piRNAs. Because it is thought that the size of small guide RNAs is determined by the footprint of small RNA-binding Argonaute/PIWI proteins (Wang et al., 2008) , the change in the size of PIWIL1-associated piRNAs implies a change in the structure of the protein that accommodates short piRNAs. The question is how the change in the structure of PIWIL1 is induced to either unload long piRNAs or reload short piRNAs or to conclude the production of short piRNAs, but not long piRNAs, from intermediates, of which the 5´-ends are likely anchored within the MID-PIWI domain of PIWIL1. It is known that the release of the 3´-end of the guide small RNA from the PAZ domain of some bacterial Argonaute proteins occurs during target recognition, which is accompanied by conformational changes in the PAZ domain (Kaya et al., 2016; Sheng et al., 2014) . A recent study has also shown that disengagement of the small RNA 3´-end from the PAZ domain occurs during the mammalian Argonaute activity cycle (Baronti et al., 2020) . Thus, it is tempting to speculate that conformational changes in PIWIL1, either upon target recognition of long piRNA-loaded PIWIL1 or by hitherto unknown mechanisms, may occur to conclude the production of short piRNAs.

In other words, PIWIL1 could switch between states of structural rearrangements to produce two types of piRNAs. In this context, it is interesting that short piRNAs are loaded onto PIWIL3 when the protein is heavily phosphorylated. Indeed, we demonstrated that phosphorylation is required for PIWIL3 to associate with piRNAs. It is known that the loading of small guide RNAs onto Argonaute proteins appears to be affected by phosphorylation, although the underlying mechanisms are poorly understood (Treiber et al., 2019) . Phosphorylation could induce changes in the structure of PIWIL3 so that the protein is loaded with processed intermediates of piRNA precursors and the footprint of the protein determines the size of loaded piRNAs. Alternatively, but not mutually exclusive, PIWIL3 may need to be phosphorylated to stably hold piRNAs. Our findings suggest that hamster PIWI proteins in the oocyte can alternate between states (allostery). It will be interesting to model the structural changes in the PIWI protein using published data on structures of PIWI proteins and other Argonaute proteins.

We found that PIWIL3 binds piRNAs only in MII oocytes but appears to be free from piRNAs at other stages of oocyte maturation. There are precedents for piRNAindependent functions of PIWI proteins. Mouse PIWIL1 (MIWI) was found to bind spermatogenic mRNAs directly, without using piRNAs as guides, to form mRNP complexes that stabilize mRNAs essential for spermatogenesis (Vourekas et al., 2012) .

Recent studies have also shown that human PIWIL1 (HIWI) functions as an oncoprotein via piRNA-independent mechanisms (Li et al., 2020; Shi et al., 2020) . Although Argonaute/Piwi proteins tend to be destabilized when small RNAs are not loaded onto them (Elkayam et al., 2012; Haase et al., 2010; Kobayashi et al., 2019; Smibert et al., 2013) , these studies suggest that PIWI proteins may be stable in some conditions in the absence of piRNAs to play a role in some cellular functions.

In mouse testes, the main target TEs in the piRNA pathway are those of LINE1 and IAP family members. Two distinct piRNA populations are present in mouse testes: prepachytene piRNAs are enriched in TE-derived sequences (approximately 80% of the total) (Aravin et al., 2008) and pachytene piRNAs have a higher proportion of unannotated sequences, with diminished contribution from TE-derived sequences (approximately 25%) (Aravin et al., 2006; Girard et al., 2006) . These piRNAs guide PIWI proteins to target LINE1 and IAP elements and repress them either by cleaving their transcripts in the cytoplasm or by modifying their chromatin loci in the nucleus (Iwasaki et al., 2015; Ozata et al., 2019) . Recent studies also support a model in which TE-derived piRNAs serve as a guide to selectively target young L1 families for de novo DNA methylation (Molaro et al., 2014) or H3K9me3 modification (Pezic et al., 2014) in fetal testes. However, the Slicer activity of both MIWI and MILI is still required to maintain TE silencing in the mouse testis after birth, suggesting that continuous cleavage of TE transcripts by the Slicer activity is essential for repressing TEs in mouse testes (Reuter et al., 2011; De Fazio et al., 2011) . We found that the contents of PIWIL1-and PIWIL3associated piRNAs in hamster oocytes are similar to those observed in mouse pachytene piRNAs. However, the major target TEs in the piRNA pathway in hamster oocytes appear as endogenous retroviruses (ERVs), including ERV2 families. Indeed, our results are concordant with the fact that 41.5% of recently active LTR retrotransposons are accounted for by ERV2 families such as ERV2-7_MAu and ERV2-5_MAu. The differences in piRNA targets between testes and oocytes suggest the possibility that the activity of IAP and these ERV2 elements might be distinctively controlled and their young copies in the genome might have jumped in different types of germline cells.

Recent studies have shown that ERVs, which constitute a large portion (8~10%) of mammalian genomes (Mandal and Kazazian, 2008) This also implies that spermatogonial dysfunction, observed in PIWI-deficient mice, may not simply be due to deregulation of TEs but due to transcriptional rewriting. It will be interesting to see if the lack of PIWIL1 or PIWIL3 leads to dysfunctions in hamster oocytes.

We found that nearly 80% of piRNA clusters corresponding to PIWIL1associated piRNAs in testis, ovary, and MII oocytes of hamsters were sex-specific. This is in agreement with previous studies with total small RNA-seq of ovaries, indicating that piRNAs in human and bovine ovaries mostly share common piRNA clusters with piRNAs in testes (Roovers et al., 2015; Williams et al., 2015) . However, a recent study showed that only about 30% of human oocyte piRNA clusters overlapped with the human testis-piRNA clusters, proving that testes and oocytes differentially express piRNA clusters (Yang et al., 2019) . Our study, together with that of human oocyte piRNAs, suggest a model in which the expression of oocyte and testis piRNA clusters have different functions with distinct targets. Transcriptional factors that regulate their expression are also distinct, though further analysis, including ATAC-seq to define transcription start sites of piRNA precursors, will be needed to identify the transcriptional factors responsible for piRNA clusters in hamster oocytes. In addition, we demonstrated that nearly half of the two populations of PIWIL1-associated piRNAs in oocytes share common clusters and that nearly half of PIWIL3-associated piRNAs in MII oocytes map to clusters from which PIWIL1-associated piRNAs are derived. These results suggest the possibility that a considerable portion of oocyte piRNA cluster transcripts are processed by a common biogenesis pathway but are loaded onto distinct PIWI proteins.

In summary, we have demonstrated that piRNAs are abundantly expressed in hamster oocytes and their size changes during oocyte maturation. Given the recent development of methods to produce genome-edited hamsters (Fan et al., 2014; Hirose et al., 2020) , our findings have set the stage for understanding the role that the piRNA pathway plays in mammalian oocytes. Our newly assembled golden hamster genome will also greatly promote the use of golden hamster as a model to understand human disease.

Ovaries were collected from 4-week-old golden hamsters that were obtained from the Sankyo Labo Service Corporation, Inc. MII oocytes were collected from 8-week-old golden hamsters that were injected with 40 U of equine chorionic gonadotropin (Serotropin; ASUKA Pharmaceutical Co., Ltd., Tokyo, Japan) at estrus. Two-cells were collected from 8-week-old golden hamsters that were injected with equine chorionic gonadotropin at estrus and mated with adult male hamsters.

Total RNA for PIWIL3 5′ RACE was extracted from the ovaries of 8-week-old golden hamsters using ISOGEN (Nippon Gene) and RNeasy (Qiagen). 5′ RACE was performed using the SMARTer RACE 5/3 kit (Takara Bio) according to the manufacturer's instructions. The PCR-amplified PIWIL3 5′ RACE fragments were cloned into the pRACE vector and sequenced.

To construct the golden hamster PIWI gene expression vectors, PIWI genes were amplified by PCR using gene specific primers and 3-week-old hamster testis and ovary cDNA. The PCR products were digested with SpeI and NotI, and then cloned into pEF4-3xFlag vector.

Anti-PIWI monoclonal antibodies were produced as described previously Nishida et al. 2007) is available for all these applications.

Western blot analysis was performed as described previously (Miyoshi et al. 2010) . Onetenth of protein from ovaries, proteins from 15 oocytes and 2-cell embryos, and one-tenth of protein from the testes were subjected to SDS-PAGE. Culture supernatants of antimarmoset PIWIL1 (MARWI) (1A5-A7) hybridoma cells (Hirano et al. 2014 ) and antihamster PIWIL3 (3E12) hybridoma cells (1:500) and mouse monoclonal anti-β-tubulin (1:5000, DSHB, E7) were used.

The ovaries from 8-week-old wild type female golden hamsters were embedded without fixation for frozen sections. Freezing blocks were sliced to 8 µL and treated with an anti-

Alexa488-conjugated anti-mouse IgG (Molecular Probes) was used as the secondary antibody. Fluorescence was observed with an IX72 fluorescence microscope (Olympus).

Immunoprecipitation details have been described previously (Saito et al. 2009 ). Briefly, the ovaries were homogenized using TissueLyser II (QIAGEN) and oocytes and 2-cell embryos in which the zona pellucida were eliminated using polyvinyl acetate with acetic 

Periodate oxidation and β-elimination were performed as described previously (Hirano et al. 2014; Kirino and Mourelatos 2007; Ohara et al. 2007; Saito et al. 2007; Simon et al. 2011) . A 100 µL mixture consisting of purified RNAs and 10 mM NaIO4 (Wako, 199-02401) was incubated at 4°C for 40 min in the dark. The oxidized RNAs were then incubated with 1 M Lys-HCl at 45°C for 90 min to achieve β-elimination.

The 5′-end of the RNAs were labeled with [gamma-32 P] ATP (Perkin Elmer) and T4 polynucleotide kinase (ATP: 5-dephosphopolynucleotide 5'-phosphotransferase, EC 2.7.1.78). The labeled RNAs were cleaned using Micro Bio-Spin™ column 30 (Bio-Rad) and ethanol precipitation. The precipitated RNAs were subjected to SDS-PAGE with 15% urea.

Raw PacBio reads (Table S1 ) were assembled into contigs using the FALCON genome assembler, which is widely used for processing long reads (Chin et al., 2016) . To correct assembly errors in the FALCON contigs, we applied the Racon module (Vaser et al., 2017 ) three times. To obtain a chromosome-scale genome assembly, we aligned all contigs to the 22 golden hamster chromosomes (MesAur1.0_HiC, DNA Zoo) ( Figure S7 ).

We attempted to fill 3,797 gaps in the reference chromosomes using the FALCON contigs and error-corrected reads. We generated error-corrected reads using the FALCON assembler, which aligned PacBio raw reads to each other, obtained the consensus sequence of aligned reads using multiple alignments, and then output the consensus sequences as error-corrected reads. To fill unsettled gaps, we aligned error-corrected reads to the gaps using the minimap2 program (Li, 2018) and manually inspected the results to determine whether the gaps were spanned by more than one error-corrected read (Table S2 ). Finally, we polished the draft assembly using the PacBio raw reads and Arrow program.

We used the FALCON genome assembler version 2018.08.08-21.41-py2.7-ucs4-beta (Chin et al. 2016 ) with default parameter settings to assemble PacBio reads. To correct errors in the assembly, we applied the RACON module (version 1.4.10; Vaser et al. 2017) three times with default parameter settings.

We aligned all contigs in the assemblies to the golden hamster reference assembly (DNA Zoo release MesAur1.0_HiC.fasta.gz) using MUMmer 4.0.0beta2 software (Marçais et al. 2018) with the nucmer program, and the following parameters: mum min cluster = 100, max gap = 300. We also used the minimap2 version 2.13 (Li, 2018) with default parameter settings.

We polished the draft genome using the Arrow program (version 2.3.3; https://github.com/PacificBiosciences/GenomicConsensus) with default parameter settings. We used the QUAST 5.0.2 tool (Gurevich et al. 2013) to calculate mismatch and indel ratios for our golden hamster genome with respect to the DNA Zoo Hi-C genome, both before and after genome polishing (Table S3) .

We lifted gene structures and other genome annotations from the golden hamster reference assembly (MesAur1.0 release 100) to our golden hamster assembly. Such crossassembly mapping typically requires an annotation file in a standard format (e.g., GFF3; https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md), chain alignment (Kent et al. 2003) , and a program capable of mapping annotations based on chain alignment (e.g., liftOver) (Speir et al. 2016 

First, we identified synteny blocks that were shared between two or three species, that is, hamster, mouse, and/or rat (details provided in Figure S8 ). Next, we inferred ancestral karyotypes of Muridae and Cricetidae ( Figure 3B ) by integrating synteny blocks according to maximum parsimony, minimizing the total number of chromosomal rearrangements such as chromosome fusions, chromosome fissions, and translocations.

We excluded inversions from chromosomal rearrangements, which were markedly more prevalent than the other rearrangements. We compared our nearly complete golden hamster genome with the mouse (Mus musculus) sand rat (Rattus norvegicus) reference genomes. Figure 4A shows conserved synteny blocks between the golden hamster, mouse, and rat genomes. We inferred ancestral karyotypes by integrating synteny blocks shared between two or three species according to maximum parsimony to minimize the amount of chromosomal rearrangement. We obtained a high-resolution ancestral karyotype of Muridae, including mice and rats, using the golden hamster genome as the outgroup as well as a species ancestral Cricetidae karyotype ( Figure 3B ).

We first used the RepeatModeler ver. 2.0 (Flynn et al, 2020) coupled with LTR_retriever ver. 2.8 (Ou and Jiang, 2018) for de novo identification of repetitive elements in the hamster genome. Among the initial repeat candidates obtained, 282 elements, covering ~37% of the genome in total, were used for subsequent in-depth characterization. In the accurate identification process, we conducted a BLASTN search and collected 80-100 sequences along with the 6-kbp flanking sequences of each element. The sequences were aligned with MAFFT ver 7.427 (Katoh and Standley, 2013) followed by manual curation, and a refined consensus sequence was constructed to be used for the next round of blastn search. This process was repeated for three rounds at maximum until the consensus sequence reached the 5' and 3' termini. We finally constructed 177 consensus sequences at the subfamily level and classified them based on the sequence structure and a comparison with known elements using CENSOR (Jurka et al. 2005) and RepeatMasker ver. 4.1.0. We finally constructed a custom repeat library by adding 177 novel consensus sequences to the original rodent repeat library. RepeatMasker analysis was conducted for genome annotation using the cross-match engine with the sensitive mode (-s). Young (i.e., recently active) complete TEs were selected based on the <5.0 K2P divergence from the consensus sequence and the full-length insertions, although ignoring the lack of sequence homology in up to 50 bp of the 5-terminal of LINEs and 3 bp at the edge of other TEs.

PIWIL1 and PIWIL3-associated piRNAs were prepared as described previously (Hirano et al. 2014 RNAs. PIWIL3-piRNAs with 18-20 nt were selected for further analysis.

The RNA-seq library in oocytes was prepared using the SMART-Seq® Stranded Kit 

Sequence logos were generated using the motifStack R package (http://www.bioconductor.org/packages/release/bioc/html/motifStack.html). Small RNA sequences were aligned at the 5´-end, and nucleotide bias was calculated per position.

Annotation of genome mapped reads was determined as described previously (Iwasaki et al., 2017) with some modifications. We examined the overlap between read-mapped genomic regions and feature track data. Each feature data was obtained from and UCSC LiftOff pipeline for protein-coding genes. Reads were assigned to a feature when the length of its overlap was longer than 90% of the small RNA. The priority of the feature assignment was defined to avoid any conflict of the assignment. For this study, the priority was in the following order: transposon, repeat, miRNA, rRNA, tRNA, snRNA, snoRNA, and protein-coding genes (exons and introns). Any unassigned regions were regarded as unannotated regions.

Prediction of PIWI-piRNA targets derived from TE regions was determined as described previously (Hirano et al. 2014; Iwasaki et al. 2017 ) with some modifications. First, we eliminated piRNA reads that mapped to tRNAs or rRNAs. The extracted reads were aligned to consensus sequences of transposons (Rodentia custom library), allowing two mismatches. The alignment was performed using Bowtie because of the large number of obtained sequences.

Prediction of piRNA clusters was performed using proTRAC version 2.4.3 (Rosenkranz et al. 2012 ) under the following conditions as described previously (Yang et al., 2019) :

(1) more than 75% of the reads that were appropriate to the length of each piRNA; (2) more than 75% of the reads exhibited the 1 U or 10 A preference; (3) more than 75% of reads were derived from the main strand; and (4) -pimin option with 21, 28, and 18 for oocyte PIWIL1-piRNAs and PIWIL3-piRNAs, respectively.

We analyzed the motif sites surrounding the transcription start sites of the testis and ovary piRNA clusters. We first extracted sequences surrounding transcriptional start sites of unidirectional and bidirectional piRNA clusters predicted from small RNA-seq and RNAseq mapping data. For the detection of bidirectional piRNA clusters, we first determined the number of reads per base in the cluster based on the results of RNA-seq with TopHat version 2.1.1 (Kim et al., 2013) . We next detected the direction of each base site and a region in which the same direction was contiguous by more than 200 bp was identified.

If (+) or (-) occupies more than 75% of the cluster, the cluster is designated as the direction.

If the number of reads was less than five, the direction was the same as the previous base (Trapnell et al., 2010) . The direction of the sequence strands was the same as in the transcripts. We then used these sequences and performed motif searches using MEME version.5.1.0 (Bailey et al., 2009) . Tomtom version 5.1.1 (Bailey et al., 2009 ) was used for motif comparison.

To visualize the read density obtained from smRNA-seq snd RNA-seq, we created a BigWig file by using HOMER version 4.11 (Heinz et al. 2010 ) and displayed the Integrative Genomics Viewer (IGV) version 2.4.1. (Robinson et al., 2011) . The normalized expression level of each sample was calculated using reads per million reads (RPM).

The accession number for the deep-sequencing datasets reported in this paper is PRJDB10770 in DDBJ.

We thank all members of the Siomi laboratory, especially Drs. Hirotsugu Ishizu and Akihiko Sakashita, for the discussions and comments on this study. We are also grateful to Yukiteru Ono for assisting the bioinformatics. S. Morishita thanks Drs. Erez Brown blocks show the existence of a proto-chromosome in the ancestral Cricetidae karyotype. We confirmed the previous speculation that mouse chromosomes 5 and 16 obtained blocks from a common Muridae proto-chromosome (light-orange), and that chromosomes 10, 15, and 17 obtained blocks from a common Cricetidate protochromosome (maroon blocks). We also identified two groups of mouse chromosomes (6, 8) and (8, 13) having large blocks from common Muridae proto-chromosomes. Figure S1 . Related to Figure 1 .

(A) Expression levels of homolog genes in hamster testis, ovary, and MII oocyte, which is known to correspond to the PIWI-piRNA pathway in mouse and Drosophila.

Expression levels are normalized by transcripts per kilobase million mapped (TPM).

(B) A phylogenetic tree of PIWI genes in hamster, human, mouse, and Drosophila. Their evolutionary distances were calculated using the maximum likelihood method.

(C) Western blotting was performed on 293T cells, which were transfected with Piwil1

and Piwil3 genes, respectively. The results show that anti-MARWI (marmoset PIWIL1) and hamster PIWIL3 antibodies specifically recognize their targets. and Tmem1 was conserved between hamsters and humans.

(B) From our hamster genome and the reference genomes of mouse (mm10) and rat (rn6), we retrieved and compared the three genomic regions between Wscd2 and Sgsm1. The

Piwil3 encoding region is present in the focal region of the hamster genome but is absent in the mouse and rat genomes. The Piwil3 region is deleted in the mouse genome and is replaced with another sequencing in the rat genome.

(C) We aligned to PacBio long reads to the assembled contig with Piwil3 (the red-colored region). In the lower portion, the respective blue and red colored lines show the alignments of reads in the plus and minus strands. We observed that each base position was covered by an ample number (ten or more) of long reads, and the read coverage was nearly even, thereby confirming the accuracy of the assembled contig. x-axis and the mouse (red) and rat (green) genomes on the y-axis. Statistics of raw and "polished" long genome sequencing reads from the golden hamster sample using the PacBio Sequel system. Mismatch and indel ratios of our golden hamster genome assembly before and after polishing using PacBio long reads. 

Ishino et al.

Ishino et al. the simultaneous degradation of endo-siRNAs and piRNAs in early embryos of mice 28,32 , these 20-nt os-piRNAs declined more rapidly than the 30-nt piRNAs in human morula embryos (Fig.6a) and correlated with the disappearance of HIWI3 mRNA during human embryo development ( Supplementary Fig. 4a Fig. 7 Comparison of the Argonaute family proteins associated with small RNAs in human and mous RNA (endo-siRNA), and PIWI-interacting RNAs (piRNA) are the three types of small RNAs that bi preferredfirst nucleotide for all the types of piRNAs, and adenine (A) is the preferred tenth nucleo HIWI3-associated 20-nt oocyte short piRNAs (os-piRNAs) are shorter and do not bear the 2 ′-O-methy expression of long piRNAs ( ≥ 25 nt) was relatively lower than the expression of os-piRNAs (~20 nt) in (~22 nt) in mouse oocytes 

A novel class of small RNAs bind to MILI protein in mouse testes

A piRNA pathway primed by individual transposons is linked to de novo DNA methylation in mice

Developmentally regulated piRNA clusters implicate MILI in transposon control

MEME SUITE: tools for motif discovery and searching

Base-pair conformational switch modulates miR-34a targeting of Sirt1 mRNA

Discrete small RNA-generating loci as master regulators of transposon activity in Drosophila

LTR retrotransposons transcribed in oocytes drive species-specific and heritable changes in DNA methylation

MIWI2 is essential for spermatogenesis and repression of transposons in the mouse male germline

Phased diploid genome assembly with single-molecule real-time sequencing

Regulatory activities of transposable elements: from conflicts to benefits

Piwi is a key rgulator of both somatic and germline stem cells in the Drosophila testis

A slicer-mediated mechanism for repeat-associated siRNA 5' end formation in Drosophila

QUAST: quality assessment tool for genome assemblies

Probing the initiation and effector phases of the somatic piRNA pathway in Drosophila

Noncoding RNA. piRNA-guided transposon cleavage initiates Zucchini-dependent, phased piRNA production

LINE-1 retrotransposons: modulators of quantity and quality of mammalian gene expression?

Roles for retrotransposon insertions in human disease

Genetic and mechanistic diversity of piRNA 3'-end formation

Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities

Small RNA profiling and characterization of piRNA clusters in the adult testes of the common marmoset, a model primate

Acrosin is essential for sperm penetration through the zona pellucida in hamsters

The golden (Syrian) hamster as a model for the study of reproductive biology: Past, present, and future

The Drosophila RNA methyltransferase, DmHen1, modifies germline piRNAs and single-stranded siRNAs in RISC

The structural biochemistry of Zucchini implicates it as a nuclease in piRNA biogenesis

A Drosophila fragile X protein interacts with components of RNAi and ribosomal proteins

Deep sequencing and high-throughput analysis of PIWI-associated small RNAs

PIWI-Interacting RNA: Its Biogenesis and Functions

Identification and Functional Analysis of the Pre-piRNA 3' Trimmer in Silkworms

Repbase Update, a database of eukaryotic repetitive elements

Roles of MIWI, MILI and PLD6 in small RNA regulation in mouse growing oocytes

MAFFT multiple sequence alignment software version 7: improvements in performance and usability

A bacterial Argonaute with noncanonical guide RNA specificity

Evolution's cauldron: duplication, deletion, and rearrangement in the mouse and human genomes

Adaptation to P element transposon invasion in Drosophila melanogaster

SnapShot: Vertebrate transposons

MUMmer4: A fast and versatile genome alignment system

<200-1885-3-PB.pdf>

Crystal Structure of Silkworm PIWI-Clade Argonaute Siwi Bound to piRNA

Molecular mechanisms that funnel RNA precursors into endogenous small-interfering RNA and microRNA biogenesis pathways in Drosophila

Noncoding RNA. piRNA-guided slicing specifies transcripts for Zucchini-dependent, phased piRNA biogenesis

Two waves of de novo methylation during mouse germ cell development

Gene silencing mechanisms mediated by Aubergine piRNA complexes in Drosophila male gonad

Hierarchical roles of mitochondrial Papi and Zucchini in Bombyx germline piRNA biogenesis

Structure and function of Zucchini endoribonuclease in piRNA biogenesis

PNLDC1, mouse pre-piRNA Trimmer, is required for meiotic and post-meiotic male germ cell development

The 3' termini of mouse Piwi-interacting RNAs are 2'-O-methylated

LTR_retriever: A Highly Accurate and Sensitive Program for Identification of Long Terminal Repeat Retrotransposons

PIWIinteracting RNAs: small RNAs with big functions

StringTie enables improved reconstruction of a transcriptome from RNA-seq reads

piRNA pathway targets active LINE1 elements to establish the repressive H3K9me3 mark in germ cells

piRNAs and their involvement in male germline development in mice

Murasaki: a fast, parallelizable algorithm to find anchors from multiple genomes

Miwi catalysis is required for piRNA amplification-independent LINE1 transposon silencing

Integrative genomics viewer

Piwi proteins and piRNAs in mammalian oocytes and early embryos

proTRAC--a software for probabilistic piRNA cluster detection, visualization and analysis

A regulatory circuit for piwi by the large Maf gene traffic jam in Drosophila

Pimet, the Drosophila homolog of HEN1, mediates 2'-O-methylation of Piwi-interacting RNAs at their 3' ends

Identification of eight members of the Argonaute family in the human genome

SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation

Structure-based cleavage mechanism of Thermus thermophilus Argonaute DNA guide strand-mediated DNA target cleavage

PIWIL1 promotes gastric cancer via a piRNA-independent mechanism

Pathogenesis and transmission of SARS-CoV-2 in golden hamsters

Recognition of 2'-O-methylated 3'-end of piRNA by the PAZ domain of a Piwi protein

Casein kinase II phosphorylates the fragile X mental retardation protein and modulates its biological properties

Homeostatic control of Argonaute stability by microRNA availability

The UCSC Genome Browser database: 2016 update

The biogenesis and function of PIWI proteins and piRNAs: progress and prospect

Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation

Regulation of microRNA biogenesis and its crosstalk with other cellular pathways

A distinct small RNA pathway silences selfish genetic elements in the germline

Classification and characterization of human endogenous retroviruses; mosaic forms are common

Fast and accurate de novo genome assembly from long uncorrected reads

Defective germline reprogramming rewires the spermatogonial transcriptome

Mili and Miwi target RNA repertoire reveals piRNA biogenesis and function of Miwi in spermiogenesis

Structure of the guide-strand-containing argonaute silencing complex

BUSCO Applications from Quality Assessments to Gene Prediction and Phylogenomics

Discovery and Characterization of piRNAs in the Human Fetal Ovary

Crystal structure of Drosophila Piwi

Single-cell CAS-seq reveals a class of short PIWI-interacting RNAs in human oocytes

The authors declare that they have no conflict of interest. cannot distinguish between whether end 2′-O-methylation are strictly ey are also present in other cells in the gle-cell small RNA-seq method protudy of small RNAs in rarer cell types, s. To improve the sensitivity of small ecting small RNAs expressed at a low loped CAS-seq, which employs Cas9eterodimers and IVT linear ampli fis during ampli fication. Compared to a A-seq method 59 , CAS-seq provides a mapping ratio and can reliably detect small RNAs expressed at low levels, siRNAs. CAS-seq strictly relies on the uence for spCas9/sgRNA recognition er. Fortunately, the most frequently NA-seq, such as TruSeq (Illumina, d with the TGG, which is the same S-seq. Other Cas family proteins equences may be used to circumvent ns for custom-designed adapters. In NA/DNA hybrid substrates by Cas9 spectrum of spCas9, which could also other applications. oocytes were approximately 20 nt in than the known piRNAs in the testes 9 Table S3