To characterize piRNAs expressed at various stages, we analyzed the small RNA sequencing data for PSG, SG, LZ, PD, and RS (Ichiyanagi et al., 2014; Inoue et al., 2017). The PSG data were obtained from samples at embryonic day 16.5, and the SG, LZ, PD, and RS data were for FACS-purified germ cells of the same individual at 21 days postnatal. After removing reads of miRNAs and cellular RNAs, such as rRNAs and tRNAs, the retained reads were mapped onto the mouse reference genome. The length distribution of the mapped piRNAs differed between the stages (Figure 1B). In PSG, the length peak was 27 nt, possibly reflecting a mixture of PIWIL2-and PIWIL4-bound piRNAs. In SG, the peak was shifted 1 nt shorter (26 nt), consistent with the predominance of PIWIL2-bound piRNAs. In LZ, the peak remained 26 nt, however, piRNAs of ≧29 nt became more abundant. In PD, the 26-nt peak became smaller, whereas the 30-nt peak became abundant. Eventually, in RS, only a single peak of 30 nt was observed, consistent with the predominance of PIWIL1-bound piRNAs. Therefore, shorter piRNAs were initially generated in male germ cells, and the population gradually shifted toward longer piRNAs during development. Notably, although long piRNAs are features of so-called “pachytene” piRNAs, our analysis revealed that such long piRNAs start to be generated at the LZ stages. To avoid ambiguity, we designated these longer piRNAs observed in spermatocytes and spermatids as “late piRNAs” (because they became expressed late) and shorter piRNAs in SG and spermatocytes as “early piRNAs” (because they were expressed early). As stated below, prospermatogonial piRNAs are distinct from these postnatal piRNAs; we designate these piRNAs as “PSG piRNAs”.

Certain piRNAs are generated via the ping-pong-cycle reactions for repeat-derived RNAs, whereas the majority are generated from precursor RNAs transcribed from piRNA clusters. To study the expression dynamics of cluster-derived piRNAs, we compiled piRNA clusters based on previous reports (Kuramochi-Miyagawa et al., 2008; Li et al., 2013) and piRNAdb (https://www.pirnadb.org), yielding a list of 512 piRNA clusters. The mapped reads in each developmental stage were counted and normalized by the genomic length of a cluster, yielding reads per kilobase per million mapped reads (RPKM values). K-means clustering revealed three major groups (Figure 2A; Supplementary Table S1): 1) clusters showing high expression in PSG but low expression in postnatal stages (PSG piRNA clusters), 2) clusters showing increased expression in SG and LZ but decreased expression in PD and RS (early piRNA clusters), and 3) clusters whose expression gradually increased as development proceeded (late piRNA clusters).

Next, to gain mechanistic insights, we analyzed the lengths of piRNAs derived from each cluster at each developmental stage (Figures 2B–D). Regardless of the clusters from which they were derived, piRNA length was primarily determined by the developmental stage. For example, in PSG and SG, the piRNA lengths had a peak at approximately 27 nt, whereas the peak was at 30 nt in RS. It has been shown that piRNAs are processed by an endoribonuclease after loading onto PIWI proteins, and the structure of PIWI proteins determines the length of mature piRNAs that they bind (Kawaoka et al., 2011; Vourekas et al., 2012; Yamaguchi et al., 2020). The distribution of piRNA length is consistent with the idea that the relative number of PIWI-like proteins (PIWIL1 vs. PIWIL2) existing in the cell dictates the piRNA length. This implies that piRNAs of any cluster can bind to PIWIL1, PIWIL2, and PIWIL4 if these proteins are present in the cell.

However, we noticed that the PSG and early clusters generated shorter piRNAs than that of the late cluster in postnatal germ cells. In SG and LZ, the PSG and early clusters showed a single peak at approximately 27 nt (Figures 2B,C, green and yellow), whereas the late clusters showed a bimodal distribution with two peaks at 27 and 30 nt (Figure 2D, green and yellow). In PD, the PSG and early clusters showed two peaks at 27 and 30 nt (Figures 2B,C, orange), whereas the late clusters showed a predominant peak at 30 nt (Figure 2D, orange). In RS, the length of most piRNAs was approximately 30 nt (Figures 2B–D, red), however, a substantial fraction of piRNAs from the PSG and early clusters was approximately 27 nt in length (Figures 2B,C, red). These results suggest that piRNAs of each cluster have a preference for PIWI-like proteins to bind. The sequence features of precursor RNAs may be involved in the difference in preference between the PSG/early and late clusters.

It has been shown that piRNAs in testes at P10 are predominantly 25–28 nt in length (Aravin et al., 2007). Since such testes contain spermatogonia and preleptotene, leptotene and zygotene spermatocytes (but not germ cells of pachytene and later stages), these piRNAs are called pre-pachytene piRNAs. This seemed inconsistent with our data for P21 LZ spermatocytes, expressing longer piRNAs, which were derived from late piRNA clusters. It is conceivable that LZ spermatocytes at different ages produce piRNAs of different lengths. Therefore, we prepared LZ spermatocytes from P11 and P16 testes by cell sorting (Supplementary Figure S1) as for P21 testes (Inoue et al., 2017), and piRNAs in the LZ samples were compared (Figure 3). In P16 LZ spermatocytes, the expression levels of late piRNA clusters were almost the same, and the piRNA length distribution was also very similar. In P11 LZ spermatocytes, however, expression levels of late piRNA clusters were much lower than those in P16 and P21 LZ spermatocytes.

P11 testes do not contain pachytene spermatocytes, whereas P16 and P21 do. However, the presence of longer piRNAs in P16 and P21 LZ preparations cannot be explained by contamination, if any, of PD spermatocytes because of the high purity (>85%) of the cell fractions (Gaysinskaya et al., 2014; Inoue et al., 2017). LZ spermatocytes at P11 are developmentally distinct, as they are produced in the first wave of spermatogenesis, where spermatocytes are differentiated from non-stem-cell spermatogonia (Yoshida et al., 2006). On the other hand, at P16 and later, spermatocytes are differentiated from spermatogonial stem cells. These suggest that the first-wave LZ spermatocytes very weakly express the late piRNA cluster, whereas LZ spermatocytes in subsequent waves start to express these piRNAs. The LZ spermatocytes express Piwil1 (Hermann et al., 2018), which is consistent with the presence of PIWIL1-type longer piRNAs.

As expected, a substantial portion of piRNAs were mapped onto the retrotransposon sequences in the genome. Thus, we analyzed the time course of the abundance of retrotransposon-derived piRNAs (Figure 4A; Supplementary Table S2). In general, the expression of retrotransposon-derived piRNAs was highest in PSG, and gradually declined as development proceeded, which is consistent with the fact that piRNAs in PSG mainly work for post-transcriptional regulation and de novo DNA methylation of retrotransposons (Inoue et al., 2017).

Next, we analyzed the relationship between the evolutionary ages of retrotransposons and piRNA abundance. The average divergence from the respective consensus sequences was used to estimate the age of retrotransposons in the genome. In PSG, younger retrotransposons (with lower divergence) had more piRNAs (Figure 4B, R = –0.44, p < 10^–15). PSG piRNA clusters did not show an enrichment of young retrotransposon copies (Figure 4G). Thus, the high abundance of piRNAs from young retrotransposons is likely due to the relatively high expression of their mRNAs in prospermatogonia, which can fuel the ping-pong cycle amplification of piRNAs. This is also consistent with the role of piRNAs in the host defense system against the mobility of active retrotransposons in PSG. However, in SG and the later developmental stages, the relationship between age and piRNA abundance became weaker (Figures 4C–F).

Here, we compared the features of piRNAs expressed at various stages of male germ cell development in mice. “Pachytene piRNAs” observed in testes at ages 10 days and older are mainly long (29–31 nt) piRNAs, whereas “pre-pachytene piRNAs” observed in testes at ages 9 days and younger are mainly short (25–28 nt) piRNAs (Aravin et al., 2007). However, the present study revealed that, although longer piRNAs are highly expressed in RS, substantial amounts are also expressed at earlier stages, such as SG and LZ, particularly for piRNAs derived from the late piRNA clusters (Figure 2D). Likewise, characteristic “pre-pachytene piRNAs” (those from the PSG and early clusters) remain at the pachytene and later (PD and RS; Figures 2B,C). Since the length distribution and expression abundance gradually change in developmental stages, from SG via LZ and PD to RS (Figures 1, 2), it is inappropriate to categorize piRNAs (and clusters) into pachytene and pre-pachytene piRNAs based on their expression levels in whole testes. Here, we propose to divide them into three categories based on their expression patterns in male germ cells: PSG, early, and late piRNA clusters (see Supplementary Table S1 for details).

Leptotene and zygotene spermatocytes were prepared from mouse testes at postnatal days 11 and 16 (3 individuals respectively) according to the published protocol (Gaysinskaya et al., 2014; Inoue et al., 2017), with an exception that SH800 (SONY) was used for cell sorting (Supplementary Figure S1). RNAs were prepared from the sorted cells by isogene (Nippongene) and DirectoZol RNA prep kit (Zymo Research), quantified by Qubit (Thermo Fisher), and analyzed by TapeStation (Agilent technology). Small RNA-seq libraries were prepared from 31 to 23 ng of total RNAs (of P11 and P16, respectively) using the NEBNext Small RNA library Prep set (NEB). The libraries were sequenced on Illumina MiSeq in a 60-bp single-end run (yielding 14 million reads per sample, respectively). The sequencing data has been deposited to Gene Expression Omnibus (GEO) under the accession number, GSE200853.

The published dataset of small RNA sequencing (50-bp single-end sequencing) was retrieved from the Gene Expression Omnibus under accession numbers GSE70891 (for PSG at embryonic day 16.5) and GSE54515 (for SG, LZ, PD, and RS at 21 days postnatal). After clipping the adapter sequence, the reads were mapped to miRNAs, rRNAs, tRNAs, snoRNAs, snRNAs, and scRNAs by Bowtie2 (Langmead and Salzberg, 2012) without allowing mismatches using the following options: L 10 --norc--score-min L,0,0 --no-unal. The mapped reads were discarded, and unmapped reads were then mapped to the mouse reference genome (mm10) using Bowtie2 without allowing mismatches using the options “-L 10 --score-min L,0,0”. For the analysis of cluster-derived piRNAs, we extracted uniquely mapped reads and analyzed their lengths (10–40 bp).

The piRNA clusters used in the present study were obtained from published papers (Kuramochi-Miyagawa et al., 2008; Li et al., 2013) and piRNAdb (https://www.pirnadb.org/), and genomic coordinates were converted to mm10 using LiftOver obtained from the UCSC genome browser (Navarro Gonzalez et al., 2021). All clusters were merged using Bedtools with option -d 100, yielding 512 clusters.

To calculate the piRNA expression from clusters, uniquely mapped reads were counted for each cluster and normalized by the length of the cluster and the total number of mapped reads in the library (RPKM). The clusters were grouped based on the log expression of piRNAs at each stage using the k-means clustering method.

To calculate the expression of the retrotransposon-derived piRNAs, reads were mapped onto the mouse genome (mm10), and the RepeatMasker annotation (Navarro Gonzalez et al., 2021) of these regions was used to identify repeat-derived reads. Both sense and antisense piRNAs were counted. To estimate the ages of retrotransposons, the miliDiv values in the RepeatMasker track (per mille divergence from the consensus sequence) were averaged for all genomic copies belonging to respective retrotransposon subfamilies (or sub-subfamilies) that were annotated in the repName column.