Integrated small and long RNA sequencing reveals piRNA mediated transposon repression during human oogenesis

Domain: biology
Published: 2026-03-11
Editorial outcome: minor_revision

Abstract
The piwi-interacting RNAs (piRNAs) are essential for controlling transposable elements (TEs) activity in germ cells, yet their expression dynamics and functions during human oogenesis remain poorly understood. Here, we simultaneously profile small and long RNAs in individual human oocytes across four developmental stages. piRNAs, especially PIWIL3-associated short piRNAs, are the predominant small non-coding RNAs during human oogenesis. A marked increase in short-piRNAs coincide with a global reduction of TE expression, particularly LINE-1 and endogenous retroviruses (ERVs). In contrast, PIWIL1- and PSIWIL2-associated long piRNAs correlate with the downregulation of certain specific ERV subfamilies. Genomic analyses reveal that highly productive piRNA clusters have evolved asymmetric antisense insertion bias toward LINE-1 and ERVs, enabling TE families-specific regulation. Together, our study provides a valuable single-cell dataset of small and long RNA co-expression landscapes in developing human oocytes and reveals coordinated yet distinct roles of different PIWI/piRNA classes in TE repression during human oogenesis, with short-piRNAs acting as the primary and broad-spectrum suppressors, and long-piRNAs providing coordinated ERV-specific repression.

Introduction
Transposable elements (TEs) are mobile genetic elements that are capable of invading and propagating within host genomes. While TE activity has played a major role in shaping genome evolution, their mobilization also poses a threat to genomic integrity by inducing DNA double-strand breaks and facilitating illegitimate recombination1, 2. Based on their mobilization mechanisms, TEs are categorized into two major classes3, 4. Class I elements, or retrotransposons, mobilize via an RNA intermediate and include both long terminal repeat (LTR) and non-LTR retrotransposons. LTR retrotransposons primarily consist of endogenous retroviruses (ERVs) and ERV-derived solo LTR elements. Non-LTR retrotransposons comprise long interspersed nuclear elements (LINEs) such as L1, CR1, and L2; short interspersed nuclear elements (SINEs) including SINE1/7SL, SINE2/tR NA, and SINE3/5S rRNA; and a composite element SINE-VNTR-Alu (SVA). Class II elements, or DNA transposons, mobilize via a DNA intermediate. Although the vast majority of TEs, all DNA transposons included, are extinct, several retrotransposon families remain active in the human population5. These active elements can mobilize autonomously or rely on trans-acting factors encoded by other autonomous elements. Notably, subfamilies such as L1PA1, SVA, and AluY exhibit insertional polymorphisms among human individuals, reflecting ongoing transpositional activity6, 7. Among ERVs, HERV-K is the only lineage known to exhibit polymorphic insertions in modern humans8, 9. Additional support for their ongoing activity comes from somatic insertions observed in neuronal and cancer cells10-15, indicating their potential to mobilize post-zygotically and contribute to somatic mosaicism.

While TE repression is effectively maintained in somatic cells, genome-wide epigenetic reprogramming during primordial germ cell formation erases most DNA methylation marks, creating a vulnerable window for TE activation16-18. In females, this epigenetic state persists until the initiation of follicle growth19. Unchecked TE activity during this critical period can lead to genomic instability, meiotic arrest, and infertility20, 21. The piRNA pathway is a major defense mechanism against TE activity in the germline. PIWI proteins, a subfamily of ARGONAUTE, associate with piRNAs to mediate TE suppression via both post-transcriptional cleavage and epigenetic silencing22-31. In Drosophila and zebrafish, disruption of PIWI proteins leads to sterility in both males and females due to uncontrolled TE activity22, 32-36. In mice, the loss of Piwil1 (Miwi), Piwil2 (Mili), or Piwil4 (Miwi2) results in TE dysregulation and male infertility, whereas female oogenesis remains largely unaffected29, 32, 37, 38. Recent transcriptomic analyses of mammalian oocytes have revealed substantial variability in small RNA composition across species39. Interestingly, the Piwil3 gene is absent from the mouse genome and mouse oocytes primarily express endo-siRNAs rather than piRNAs, which are processed by a unique Dicer isoform (DicerO) and functionally substitute for piRNAs in TE suppression40. In contrast, human and most other mammalian oocytes predominantly express PIWIL3-associated short piRNAs (short-piRNAs), which account for most of the small RNAs in these cells39, 41, 42. Functional studies using golden hamsters as a model have demonstrated that loss of the piRNA pathway compromises oocyte competence and leads to early embryonic arrest43-45. Additionally, knockout studies of individual PIWI proteins suggest their non-redundant roles in TE suppression and gametogenesis46.

However, the temporal dynamics of piRNAs expression during human oogenesis and the functional relationship between distinct piRNA populations and TE regulation remain poorly understood. To address this gap, we profiled both small and long RNA transcriptomes in individual human oocytes across four key developmental stages. Our analyses reveal a global repression of TEs following the primary follicle stage, with short-piRNAs serving as the predominant and broad suppressors and long-piRNAs providing coordinated, ERV-specific repression. These findings highlight the distinct yet coordinated roles of piRNA populations in TE silencing and provide new insights into the regulatory landscape of human oocyte development.

Results
Profiling small and long RNAs in developing human oocytes
A total of 30 human oocytes spanning four developmental stages, primordial follicle (PrF), primary follicle (PF), germinal vesicle (GV), and metaphase II (MII), were collected. Small and long RNA transcriptomes were simultaneously profiled at the single-oocyte level using our previously established protocol39, 43, 46 (Fig. 1a, Materials and Methods). To ensure high data quality for downstream analysis, we applied stringent quality-control criteria to both the small RNA and long RNA datasets. For small RNA data, we removed samples with a useful read rate below 60% after 3’ adaptor trimming and short-read filtering, a mapped read rate below 40%, or rRNA/tRNA-derived reads exceeding 30%. For long RNA data, we excluded samples with a mapped read rate below 80%, mitochondrial-derived reads exceeding 10%, or rRNA-derived reads exceeding 10%. Following quality control filtering, 23 oocytes were retained for further analyses, comprising 5 PrF, 9 PF, 4 GV, and 5 MII stage oocytes (Supplementary Data 1). On average, 5.88 million raw reads per small RNA library and 7.66 million reads per long RNA library were obtained (Supplementary Data 2). Among these retained oocytes, an average of 87.88% of small RNA-seq reads were deemed usable, with a genome mapping rate of 71.17%. Long RNA-seq libraries showed mapping rates ranging from 87.84% to 95.38%. Across all oocytes, we detected an average of 170 miRNAs, 113,861 piRNAs, and 16,145 genes (Supplementary Data 2).

Global downregulation of TEs during human oogenesis
Unsupervised principal component analysis (PCA) of protein-coding genes expression successfully segregated oocytes into four clusters, each corresponding to a different developmental stage (Fig. 1b and Supplementary Data 3). Gene Ontology (GO) enrichment analysis of stage-specific feature genes (558 for PrF, 570 for PF, 356 for GV, and 738 for MII) revealed key biological processes associated with each stage (Supplementary Fig. 1). Genes related to RNA and protein biosynthesis were significantly overrepresented in PrF and PF stages, while pathways involved in mitotic cell cycle regulation and nuclear division were enriched in GV and MII stages, reflecting ongoing oocyte growth and maturation47.

To assess TE expression dynamics during oogenesis, we mapped long RNA-seq reads to representative human TE sequences. PCA of TE expression revealed that PrF and PF oocytes clustered together, whereas GV and MII oocytes formed separate groups (Fig. 1c and Supplementary Data 4). A progressive decline in TE-derived transcript levels was observed across developmental stages, with a marked reduction occurring after the PF stage (Fig. 1d). Specifically, 56 TE subfamilies were significantly downregulated between the PF and GV stages, while 12 subfamilies were upregulated (Fig. 1e and Supplementary Data 5). From the GV to MII stages, most TE subfamilies exhibited reductions in expression, although this reduction was less pronounced compared to the PF-to-GV transition (Fig. 1e, f).

L1 experienced the most substantial repression during oocyte development, with an average 35-fold reduction from the PrF to MII stages (Fig. 1f, g). Within the L1 clade, the primate-specific L1PA subfamily displayed the strongest downregulation, particularly L1PA2, which underwent a dramatic 129-fold downregulation (Fig. 1f). Approximately two-thirds of ERV subfamilies followed similar patterns to L1, showing progressive downregulation throughout oogenesis (Fig. 1f). For instance, HERV-L, a human specific ERV, exhibited the most substantial repression with an overall 89-fold decrease, followed by HERV-K, which demonstrated a 30-fold reduction from PrF to MII.

Despite the general trend of TE repression, approximately one-third of ERV subfamilies exhibited increasing expression from the PrF to GV stages, with eight subfamilies showing overall upregulation across oocyte development (Fig. 1f). Among these, HERVS71 exhibited the most notable increase, with a 3.4-fold upregulation. In comparison to robust silencing observed for L1 and most ERVs, Alu and SVA elements underwent relatively moderate repression, with fold changes of 4.7 and 5.4 from PrF to MII stage, respectively. Collectively, these findings indicate a hierarchical repression pattern during human oogenesis, where L1 and the majority of ERVs undergo extensive silencing, while certain ERVs, Alu and SVA elements remain relatively resilient.

piRNAs are the predominant sncRNAs during human oogenesis
Analyses of small RNA composition demonstrated that piRNAs were the most abundant sncRNAs species, accounting for 80.51% of total sncRNAs at the PrF stage and increasing to 84.19% at MII (Fig. 2a). PCA of piRNA expression profiles effectively distinguished PrF and PF oocytes, while GV and MII stages exhibited highly similar expression patterns (Fig. 2b). Among the PIWI genes, PIWIL3 was the most abundantly expressed, while PIWIL1 and PIWIL2 exhibited relatively lower expression levels, and PIWIL4 was nearly undetectable (Fig. 2c). PIWIL3 immunoprecipitation from PrF- and PF-stage oocytes, coupled with small RNA sequencing, confirmed that short-piRNAs (17–24 nt) are specifically associated with PIWIL3 (Fig. 2d). Oxidation assays further showed that these short-piRNAs lack 2′-O-methylation at the 3′ ends, as their signals were largely eliminated following oxidation (Fig. 2e). In contrast, long-piRNAs (25-35 nt), likely associated with PIWIL1 and/or PIWIL2, were resistant to oxidation (Fig. 2e), indicating 3’ terminal 2’-O-methylation. These biochemical characteristics of short- and long-piRNAs are consistent with observations in human GV and MII oocytes39, 41. The proportion of short-piRNAs increased from 44.88% at PrF to 73.21% at GV, remaining at 71.28% at MII (Fig. 2a), suggesting a predominant role of short-piRNAs in human oocyte compared to long-piRNAs. Genomic coordinate analysis showed that short- and long-piRNAs were far more likely to share 5’ ends than 3’ ends across all stages (PrF: 0.51 vs. 0.02; PF: 0.68 vs. 0.01; GV: 0.66 vs. 0.01; MII: 0.63 vs. 0.01) (Supplementary Fig. 2), suggesting that both populations are derived from the same 5’-cleaved precursors, consistent with findings in golden hamster43.

To investigate whether phased piRNA biogenesis occurs during human oocyte development, we quantified the 3′-to-5′ and 5′-to-5′ distances between piRNA ends mapped to the same genomic strand and evaluated phasing enrichment using Z-scores48. To validate this analytical strategy, we reanalyzed previously published Tdrkh⁺/⁻ and Tdrkh⁻/⁻ mouse testis datasets, in which robust piRNA phasing has been reported49. As expected, we observed a strong 1-nt peak in the 3′-to-5′ distance distribution, broad enrichment at ~30–40 nt in the 5′-to-5′ distance distribution, and a pronounced 5′-terminal uridine (1U) bias in downstream phased piRNAs (Supplementary Fig. 3a-d). Applying the same analytical framework to piRNAs derived from developing human oocytes revealed significant enrichment at a 1-nt separation in the 3′-to-5′ distance distribution (Z-scores 3.6–6.8) (Supplementary Fig. 3e), indicative of phased piRNA processing. Notably, short-piRNAs displayed stronger phasing signals than long-piRNAs. Downstream phased piRNAs showed a pronounced 5′-terminal U bias, particularly among short-piRNAs (Supplementary Fig. 3e and Supplementary Fig. 3g-h), consistent with canonical phased piRNA signatures in mammals. However, no prominent peaks were observed in the 5′-to-5′ distance distribution (Supplementary Fig. 3f). Unexpectedly, enrichment at a 0-nt separation exceeded that at a 1-nt separation in the 3′-to-5′ distance distribution (Supplementary Fig. 3e), indicating frequent 1-nt overlaps between adjacent piRNAs. These 1-nt-overlapping downstream piRNAs retained a strong 5′-terminal U bias, whereas their upstream counterparts exhibited marked enrichment for 3’-terminal uridine (Supplementary Fig. 3i, j). This distinctive overlap pattern may result from frequent 3′ non-templated uridylation of piRNAs, a possibility that warrants further investigation through perturbation of enzymes involved in piRNA 3′-end tailing.

Unlike in mouse oocytes, endo-siRNAs were undetectable at any sequenced developmental stage of human oocytes, aligning with observations in golden hamsters and other mammals39, 43-45 (Fig. 2a). miRNAs, which constituted 0.34% of total sncRNAs in the PrF stage oocytes, exhibited a modest increase after the GV stage (Fig. 2a). miRNA expression profiles grouped the oocytes into two distinct clusters: PrF/PF and GV/MII (Fig. 3a and Supplementary Data 6). Feature analyses identified 11 miRNAs highly enriched in PrF and PF stage oocytes, while 58 miRNAs were upregulated at GV/MII stage oocytes (Fig. 3c). Reactome pathway enrichment analysis revealed that these miRNAs are involved in pathways related to RUNX2, ALK, and TP53, as well as growth factor receptor signaling, second messenger systems, and MAPK signaling (Fig. 3d). Among the ARGONAUTE (AGO) genes, AGO2 was predominantly expressed, with expression progressively increasing during human oogenesis, while AGO1 and AGO3 exhibited consistently low expression (Fig. 3b). AGO4 was barely detectable in any analyzed stage, in line with its established distinctive role in spermatogenesis50, 51. Although miRNAs were expressed at lower levels in oocytes than in somatic cells, their stage-specific expression, the progressive upregulation of AGO2, and enrichment of their targets in critical signaling pathways suggest that miRNAs may contribute to oocyte development and warrant further investigation. Overall, these findings underscore the dominant role of piRNAs, particularly short-piRNAs, in human oogenesis.

Conserved TE preference of short- and long-piRNAs
To explore the potential contribution of piRNAs to global TE downregulation, we mapped short- and long-piRNAs separately to representative human TE sequences on both sense and antisense strands. Remarkably, approximately 80% of long-piRNAs were associated with ERVs across all developmental stages, compared to around 50% of short-piRNAs (Fig. 4a, b). This preference was largely driven by antisense-aligned long-piRNAs. At the PrF stage, 80.64% of all TE-aligned long-piRNAs were antisense ERV-related, compared to only 5.86% for sense ERV-related, resulting in an antisense-to-sense ratio greater than 13. Although this ratio decreased as development progressed, the bias toward antisense persisted, with a 4.7-fold preference at the MII stage. In contrast, short-piRNAs exhibited a broader association across various TE superfamilies, accompanied by a more balanced antisense-to-sense ratio, generally ranging between 2 to 3 throughout oogenesis (Fig. 4a, b).

Notably, the preference for ERV association among long-piRNAs was conserved across species. In cynomolgus monkey (Macaca fascicularis) and golden hamster MII oocytes, 72.49% and 69.69% of long-piRNAs, respectively, were ERV-related, compared to 53.21% and 43.87% of short-piRNAs (Fig. 4c, d and Supplementary Fig. 4a, b). In addition, the antisense-to-sense strand ratio for ERV-related long-piRNAs exceeded 10 in golden hamster oocytes, whereas there was almost no bias for ERV-related short-piRNAs (Supplementary Fig. 4a, b). Collectively, these results suggest that TE-related antisense long-piRNAs exhibit a conserved, preferential targeting of ERVs across species, indicating their specialized role in ERV repression. In contrast, TE-related antisense short-piRNAs exhibit broader TE targeting, likely reflecting their more generalized silencing function across multiple TE families.

The widespread role of short-piRNAs in TE downregulation
Analyses of TE-aligned short-piRNA expression revealed three distinct groups corresponding to the PrF, PF, and GV/MII stages (Fig. 5a and Supplementary Data 4). A marked upregulation of short-piRNA was observed from PrF to GV stages during oogenesis (Fig. 5b, c). From PrF to PF, 162 TE subfamilies exhibited significant increases in short-piRNAs expression, nearly twice as many as those showing downregulation. From PF to GV, 130 TE subfamilies showed upregulated short-piRNAs, while only seven were downregulated (Fig. 5c and Supplementary Data 7). Notably, the upregulation of short-piRNAs preceded the onset of substantial TE reduction, which was primarily initiated after the PF stage (Fig. 1d, e), indicating a temporal delay in the repressive effects of short-piRNAs on TE activity at the mRNA levels.

Among TE superfamilies, L1-related short-piRNAs exhibited the most pronounced upregulation throughout oogenesis, aligning with L1 being the most substantially downregulated TE superfamily (Fig. 5d, Fig. 1g, and Supplementary Fig. 5a). Following L1, ERV-related short-piRNAs displayed the second most prominent upregulation, although considerable heterogeneity was detected across different ERV subfamilies, mirroring the varied repression patterns observed in these elements. In addition, Alu- and SVA-related short-piRNAs showed relatively modest increases, mainly during the PF-to-GV transition, consistent with their comparatively moderate downregulation (Fig. 5d and Supplementary Fig. 5a). A robust inverse correlation between TE downregulation levels and the upregulation levels of corresponding short-piRNAs further supports the contribution of short-piRNAs to global TE repression (Fig. 5e).

The Ping-Pong amplification mechanism is essential for the post-transcriptional TE repression via PIWI/piRNA-mediated cleavage of target RNAs, generating complementary piRNA pairs with a characteristic 10-nt overlap at their 5’ ends52, 53. This process also exhibits an uracil (U) preference at the first nucleotide (1U) of antisense piRNAs, and an adenine (A) bias at the tenth nucleotide (10A) in sense. Indeed, high 1U and 10A preferences were observed in retrotransposons-associated piRNAs (sense strand 1U: 89.14%, 10A: 62.85%; antisense strand 1U: 87.86%, 10A: 44.31%, Supplementary Fig. 5b, c). Among TE superfamilies, Alu- and SVA-associated piRNAs exhibited the strongest Ping-Pong signatures, followed by moderate activity for L1 and weak activity for ERV (Supplementary Fig. 5d). We next examined Ping-Pong activity at the level of individual TE subfamilies. Specifically, short-piRNAs related to most of the L1, Alu, and SVA subfamilies showed strong Ping-Pong signatures (Fig. 5f), with the intensity of these signatures correlating with their progressive downregulation throughout oocyte development. However, this signature was less prevalent among ERV-related short-piRNAs, with notable exceptions such as THER1A and THER1B. Notably, approximately two-thirds of L1 subfamilies exhibited strong Ping-Pong activity (Group A) and consistently showed higher 1U and 10A proportions, greater antisense and sense piRNA abundance, and elevated transcript levels compared with those with low signals (Group B) (Supplementary Fig. 6a, b). Interestingly, genomic analysis of the top 20 piRNA clusters revealed that Group B L1 subfamilies carried more antisense insertions than Group A, while sense insertions were comparable (Supplementary Fig. 6c). Given the lower abundance of both antisense and sense piRNAs for Group B, these results suggest that primary processing and PIWI loading preferentially enhance piRNA biogenesis from Group A L1 subfamilies. These findings indicate that weak Ping-Pong activity in certain L1 subfamilies likely reflects their lower transcript abundance and inefficient primary piRNA processing. Collectively, these observations suggest that the PIWIL3-mediated Ping-Pong pathway may contribute to the extensive regulatory influence of short-piRNAs on TE repression during human oogenesis.

Although the majority of piRNAs mapped to various TE families, a notable minority (5-10%) originated from coding or other exonic non-TE genomic regions (Supplementary Fig. 5b, c). Unlike TE-related piRNAs, these piRNAs lacked the characteristic 10A bias (Supplementary Fig. 5c), suggesting they may participate in gene regulation through Ping-Pong–independent mechanisms.

Long piRNAs play a collaborative role in ERV repression
TE-aligned long-piRNA expression segregated oocytes into three distinct groups: PrF, PF, and GV/MII stages (Fig. 6a and Supplementary Data 4). However, long-piRNAs were less effective than short-piRNAs in distinguishing PrF from PF stage oocytes. Unlike the widespread upregulation of short-piRNAs during oogenesis, TE-aligned long-piRNAs showed a declining trend in expression (Fig. 6b, c). Differential expression analysis revealed a substantial decrease in TE-aligned long-piRNAs between the PF and GV stages, with 58 TE subfamilies exhibiting significant downregulation of their corresponding long-piRNAs, compared to 13 subfamilies that were upregulated (Fig. 6c and Supplementary Data 8). Despite the overall decline, 18 ERV subfamilies and two derivative LTR elements exhibited significant upregulation of their associated long-piRNAs, aligning with the preferential ERV targeting by long-piRNAs (Fig. 6d and Fig. 4a, b).

The Ping-Pong signature was either undetectable or substantially reduced for most TE-aligned long-piRNAs comparing with short-piRNAs (Supplementary Fig. 7a), including those ERVs or solo LTRs with upregulated long-piRNAs (Supplementary Fig. 7b). Analysis of TE-related piRNAs using the adult human ovary dataset from a previous study54 further confirmed that long-piRNAs exhibited weaker Ping-Pong signatures than short-piRNAs for SVA and Alu elements, whereas L1 and ERV showed very low or undetectable Ping-Pong signatures (Supplementary Fig. 8).

To examine the cooperative contributions of short- and long-piRNAs to TE repression during human oogenesis, we assessed the overlap among TEs with significantly altered mRNA levels (upregulated or downregulated), TEs with upregulated short-piRNAs, and those with upregulated long-piRNAs (Fig. 6e). Of the 96 downregulated TE subfamilies, 61.46% (59 out of 96) were associated with either upregulated short- or long-piRNAs. Among these, 93.22% (55 out of 59) were associated with upregulated short-piRNAs, while 20.34% (12 out of 59), most of which belonged to ERV subfamilies, were linked to upregulated long-piRNAs. L1HS and L1PA2 represent TE subfamilies that are primarily regulated by short-piRNAs with downregulation levels correlated inversely with a gradual increase in corresponding antisense and sense short-piRNAs and accompanied by prominent Ping-Pong signatures (Fig. 6f, Supplementary Fig. 9, and Supplementary Fig.10). For HERV-K, only sense short-piRNAs were markedly increased, and this increase was significantly correlated with reduced HERV-K transcript levels (Supplementary Fig. 9d and Supplementary Fig. 10d). Among the TE subfamilies associated with long-piRNAs, 66.67% (8 out of 12) also showed upregulation of short-piRNAs, suggesting a coordinated role for both piRNA classes in repressing these elements. HERV-L, one of the most strongly downregulated ERV subfamilies, showed an inverse correlation with both antisense and sense short- and long-piRNAs during oocyte development (Fig. 6g and Supplementary Fig. 10c), accompanied by a weak but detectable Ping-Pong signature. All these observations suggest a dominant role for short-piRNAs in broad TE repression, with long-piRNAs contributing to ERV-specific repression.

In addition, we identified eight TE subfamilies, including HERVS71, LTR25, MER11D, LTR1C1, LTR17, LTR6B, LTR14A, and LTR3B, that remained transcriptionally active during human oogenesis (Fig. 1f and Fig. 6e). Compared with TE subfamilies that were downregulated, these upregulated TE subfamilies showed higher transcript abundance but produced markedly fewer corresponding antisense and sense short- and long-piRNAs, accompanied by weaker Ping-Pong signatures and reduced 1U and 10A biases (Supplementary Fig. 11a). Genomic distribution analyses within the top 20 piRNA clusters further revealed substantially fewer antisense insertions for these upregulated TE subfamilies relative to downregulated TE subfamilies, consistent with their reduced antisense piRNA levels (Supplementary Fig. 11b). Notably, seven of these elements, with the exception of LTR1C1, are specific to the human genome. The limited or absent repression of these evolutionarily younger elements may reflect a lag in the adaptation of the piRNA/PIWI system to those recently emerged transposon threats.

Highly productive piRNA clusters show antisense bias
Primary piRNA biogenesis is tightly associated with specific genomic loci across various species25, 53, 55. To determine the genomic origins of piRNAs expressed during human oogenesis, we employed a proximity-based clustering strategy adapted from the piClust algorithm56. Using a threshold of at least four uniquely mapped piRNAs within a 2.5-kilobase (kb) window and merging overlapping regions, we identified 14,603 piRNA-generating clusters. Remarkably, the top 20 most highly expressed clusters accounted for 91.72 % of all uniquely mapped piRNAs (Fig. 7a, b, and Supplementary Data 9) and generated both short- and long-piRNAs, with short-piRNA production generally increasing as oogenesis progressed (Fig. 7c). These clusters displayed dynamic expression patterns during oocyte development, with their contributions to TE-related piRNA population varying by stage (Fig. 7a, b, and Supplementary Data 9). For example, the cluster chr4:78106722–78396420|+ ranked second for L1 antisense short-piRNA generation at the PrF and PF stages, but became the predominant source at the GV/MII stages. Similarly, the cluster chr12:114437233–114625671|− produced the second-largest number of ERV antisense short piRNAs at the PrF stage and emerged as the primary contributor after the PF stage. Notably, only two pairs of bidirectionally transcribed clusters were identified among the top 20 piRNA clusters (Supplementary Fig. 12), and expanding the analysis to the top 100 clusters revealed two additional pairs. This indicates that, in contrast to pachytene piRNAs in mouse testis57, only a small fraction of human oocyte piRNA clusters arises from bidirectionally transcribed loci, consistent with observations in golden hamster oocytes58.

TE enrichment analysis revealed a strong bias toward ERV fragments within the top 10 clusters, particularly for the antisense strand, when compared to genomic background regions (Fig. 7d and Supplementary Fig. 13 and Supplementary Data 10). This asymmetric insertion pattern facilitated the preferential generation of antisense-aligned piRNAs targeting ERVs (Supplementary Fig. 13b). This phenomenon is reminiscent of the flamenco/COM locus in Drosophila, which generates antisense piRNAs to silence gypsy, Idefix, and ZAM elements53, 59, 60. Interestingly, although L1 elements did not show a pronounced difference in antisense versus sense insertion frequencies, the top 10 piRNA clusters still produced markedly more antisense than sense piRNAs (Fig. 7d and Supplementary Fig. 13a). This suggests selective primary processing, in which L1 antisense piRNAs may be preferentially generated or more efficiently loaded onto PIWI proteins. In contrast, Alu elements were not enriched in these top clusters compared to genomic background regions (Fig. 7d and Supplementary Fig. 13c), consistent with the relatively low production of Alu-derived piRNAs and their modest repression during oocyte development.

The most productive piRNA cluster, spanning approximately 320 kb between WDR1 and ZFN518B on chromosome 4, contributed 34.05% of all uniquely mapped piRNAs. This unidirectionally transcribed cluster contained abundant antisense insertions of L1 and ERVs (Fig. 8a) and generated numerous short- and long-piRNAs targeting these elements (Supplementary Fig. 13a, b). Similarly, the second most productive cluster, also on chromosome 4 (~300 kb between CCNG2 and CXCL13), contributed 15.41% of uniquely mapped piRNAs and displayed a comparable enrichment for antisense insertion of L1 and ERVs (Supplementary Fig. 13a, b and Supplementary Fig. 14a). By contrast, the third most productive cluster featured minimal L1 context but was enriched for ERVs, producing predominantly antisense piRNAs targeting ERVs (Supplementary Fig. 13a, b and Supplementary Fig. 14b). This suggests that this cluster primarily functions as a specialized regulator of ERVs, distinct from the broader TE repression roles observed in the top two clusters.

SVA, a relatively young TE element with approximately 4,000 annotated copies in the human genome, was represented by antisense insertions in one of the top 20 clusters. This insertion was accompanied by the production of SVA-targeting antisense piRNAs (Fig. 7d, Fig. 8b, and Supplementary Fig. 13d), suggesting that this cluster may serve as a SVA suppressor. Conversely, L2 elements, which are largely inactive in the human genome, were absent from the top clusters (Supplementary Fig. 13e). These patterns suggest that evolutionary pressures preferentially shape the composition and activity of piRNA clusters based on the functional relevance and activity of specific TE families. Taken together, these findings indicate the unique evolutionary dynamics of human piRNA clusters and their adaptive specialization in repressing diverse and evolving TE populations during oocyte development.

Methods
Ethics statement
The study was approved by the Reproductive Study Ethics Committee of Shanghai Zhongshan Hospital (approval number: B2018-244). Written informed consent was obtained from all participants prior to oocyte donation, and no compensation was provided, following clinical protocols of the Department of Assisted Reproduction at Zhongshan Hospital, affiliated with Fudan University, Shanghai, China. All human participants were female, as the study focused on oocyte development. No sex or gender analysis was conducted because it was not applicable. The original tissue samples and oocytes were fully consumed and are not available for further use.

Human oocyte collection
Human follicles were isolated from fresh ovarian cortex obtained from patients undergoing ovarian cyst or cystic teratoma resection, following previously described procedures78. Briefly, ovarian cortical tissues were fragmented into pieces measuring 30–600 μm and digested with 1 mg/mL collagenase type IA (Sigma-Aldrich) and 8 IU/mL DNase I (Sigma-Aldrich) in 10 mL of M199 medium (Gibco) at 37°C for 30 minutes with gentle agitation. Follicles were subsequently isolated and evaluated under a stereomicroscope (Nikon SMZ 1000), and their diameters were measured using an inverted microscope (Nikon Ti2).

Ovarian stimulation and oocyte collection were performed as previously described41, 78. GV oocytes were collected following stimulation. MII stage oocytes were obtained after overnight incubation to facilitate maturation, characterized by the presence of a distinct first polar body under microscope observation.

For all oocytes at all developmental stages, granulosa cells were removed via incubation with 1 mg/mL hyaluronidase. The zona pellucida was eliminated using Acidic Tyrode's Solution (Millipore), yielding denuded oocytes, which were then washed and transferred to PCR tubes containing 4 μL of cell lysis buffer (0.2% Triton X-100, Sigma; 1 U/μL RiboLock RNase inhibitor, Thermo Fisher). The oocytes were stored at -80°C until further processing.

Single-oocyte small-RNA library construction
Single-oocyte small-RNA libraries were constructed as previously described43. Briefly, lysates from the single oocyte were incubated at 72 °C for 3 minutes to release and unfold the small RNAs. After 3’ adapter ligation, 5 U of lambda exonuclease and 25 U of 5′ deadenylates were used to remove the excess 3’ adapter. Subsequently, the 5′ adapter was ligated, followed by reverse transcription and pre-amplification. For final amplification, 1 μL of the pre-amplified product was used as a template. Libraries were size-selected by electrophoresis on a 6% polyacrylamide gel, and 130–160 bp DNA fragments were extracted to recover small RNA libraries. Sequencing was performed on the HiSeq X Ten platform (Illumina) using 2 × 150 bp reads. Oligonucleotide sequences are provided in the Supplementary Data 11.

Single-cell mRNA library construction
Single-oocyte mRNA was captured and amplified using Smart-seq2, as previously described43, 79. In brief, single oocytes or embryos were incubated at 72°C for 3 min to release the RNAs. For every single oocyte or embryo, 0.125 μL 1/100,000 dilution of the ERCC RNA Spike-In Mix (Invitrogen, Cat# 4456740) was used. Following reverse transcription and PCR pre-amplification, cDNA libraries were purified by Agencourt Ampure XP beads (0.8:1 ratio) (Beckman Coulter, Cat#A63881). 2 ng of purified cDNA was subjected to a tagmentation reaction with Tn5 transposase80. Sequencing was carried out on the Illumina HiSeq X Ten platform with 2x150 bp reads. Oligonucleotide sequences are provided in the Supplementary Data 11.

Immunoprecipitation
Five PrF or PF oocytes were lysed on ice for 10 min in lysis buffer containing 50 mM Tris-HCl (pH 7.4), 150 mM NaCl, 1 mM EDTA, 0.5% NP-40, 0.5 mM dithiothreitol, 0.1 U μL−1 RNase Inhibitor (Thermo Fisher Scientific) and 0.4 U μL–1 proteinase inhibitor cocktail (Sigma-Aldrich). Dynabeads Protein G (Thermo Fisher Scientific) were washed once with lysis buffer and incubated with 1 μg of homemade anti-PIWIL3 or rabbit IgG at room temperature for 40 min. The oocyte lysates were mixed with antibody-coupled beads and rotated gently at 4 °C for 5 h. The beads were then washed three times with lysis buffer and once with PBS, and RNA was eluted from the beads by adding 4 μL of nuclease-free water and followed by incubation at 72 °C for 5 min with gentle agitation. The supernatant containing the RNA was collected using a magnetic holder and used for small RNA library construction.

NaIO oxidation of oocyte RNAs
Total RNAs were extracted from five PrF and PF oocytes using TRIzol Reagent (Ambion) supplemented with 8.14 × 10−3 pmol of spike-in oligos. RNA pellets were dissolved in 4 μL of water. Half of the RNA was reserved as a non-oxidation control, and the remaining half was treated with an oxidation reaction mixture containing 1.5 μL of 100 mM NaIO4 (Sigma-Aldrich), 2 μL of 5× borate buffer (150 mM borax, 150 mM boric acid, pH 8.6), and 4.5 μL of water. The reaction mixtures were incubated in the dark at room temperature for 30 min, followed by RNA precipitation with 600 μL of 100% ethanol, 25 μg of linear acrylamide, and 30 μL of 3 M sodium acetate (pH 5.2). Samples were incubated at −30 °C for 30 min, centrifuged, and washed twice with 75% ethanol. The RNA pellets were dissolved in 2 μL of water and used for small RNA library preparation.

Data processing of small RNA-seq and long RNA-seq
For small RNA-seq, raw read1 sequences were processed using the fastx_toolkit (v.0.0.14) for quality filtering and adapter trimming. Reads lacking the 3’ adapter sequence or shorter than 17 bp post-trimming were discarded. Reads of 17-40 bp were mapped to the human genome (hg19) using Bowtie (v.1.2) with parameters --k=100 --v=0.

For long RNA-seq, raw pair-end reads were processed using TrimGalore (v.0.5.0) for adapter trimming and quality filtering. The remaining reads were mapped to the human genome (hg19) and representative TE sequences using STAR (v2.7.1a) with parameters --winAnchorMultimapNmax 10000 --outFilterMultimapNmax 10000. Gene abundance was quantified using RSEM (v.1.2.31) with default settings.

For multiply mapped piRNAs or long RNA-seq reads, each read was assigned to all possible mapped loci and categorized by TE superfamily in the following order: LINE/L1 → SINE/Alu → LTR/ERV → LTR/others → SVA → others. If a read or piRNA species mapped to multiple TE subfamilies within the same superfamily, it was counted in each TE family. When summarizing results at the superfamily level (e.g., LINE/L1), each piRNA or read was counted only once. This uniform assignment strategy was applied to all TE-related analyses, including TE transcript abundance, TE-associated short and long piRNAs, and strand-specific subsets, to ensure consistent quantification and accurate cross-dataset comparison.

Categorization of small non-coding RNA
Mapped small RNA-seq reads were further aligned to annotations of small non-coding RNAs in the following order: pre-miRNAs, tRNAs, rRNAs, snoRNAs, and snRNAs, using Bowtie (v.1.2) with parameters --a --v=0 --norc. Reads that precisely matched the 5’ start site of miRNAs and 3’ ends with less than 3 bp deletions or additions derived from pre-miRNAs were considered for miRNA quantification.

Unassigned sequences that were not categorized as pre-miRNA, tRNA, rRNA, snoRNA, or snRNA, with lengths ranging from 17–40 nt, were identified as piRNA clusters using a modified piClust with MinReads=4 and Eps=2500 bp. Overlapping clusters were merged using mergeBed from bedtools (v.2.30.0). Reads mapping within these clusters were defined as piRNA candidates. piRNA candidates were subsequently aligned to representative TE sequences using Bowtie with the parameters --a --v = 2. Based on alignment orientation, piRNAs were classified as TE-aligned antisense piRNAs (mapping to the opposite strand of TEs) or TE-aligned sense piRNAs (mapping to the same strand of TEs). Reads were further categorized by TE superfamily in the following order: LINE/L1, SINE/Alu, LTR/ERVs, LTR/others, SVA, and others. If a read mapped to multiple TE subfamilies within a superfamily, it was counted simultaneously for each subfamily.

Cell Clustering and Feature Gene Analysis
Cell clustering based on gene, TE, miRNA, and piRNA expression levels was performed using principal component analysis in Seurat (v.2.3.4). Feature genes for each stage were identified using FindAllMarkers from Seurat, and enriched pathways were analyzed using Gene Ontology website.

Differential Expression Analysis
Differential expression of TEs and their aligned piRNAs was analyzed using DESeq2, with thresholds of fold change ≥ 1.50 or ≤0.67 and p-value ≤0.05. TEs or aligned piRNAs with average expression levels below 100 RPM across all samples were excluded from this analysis.

Identification of Ping-Pong signature
Ping-Pong signatures were detected by calculating the Z score of 10-nucleotide overlaps of piRNAs from opposite strands. A Z score > 1.96 was considered significant57 (p < 0.05).

Sources of Annotations
Human (hg19) and cynomolgus monkey (Macaca fascicularis, MacFas5) genome sequences were obtained from the UCSC Genome Browser, while the golden hamster genome (MesAur1.0) was retrieved from Ensembl. For human and cynomolgus monkey, the representative sequences of TEs were retrieved from RepBase, and genomic TE annotations were sourced from UCSC. For the hamster, TE sequences and annotations were obtained following previously published protocols43, 46. Gene annotations for hg19 were sourced from GENCODE (version 19). Functional RNA annotations for human were obtained from miRBase (pre-miRNAs and miRNAs), Genomic tRNA Database (tRNAs), NCBI GenBank (18S and 28S rRNAs), and Ensembl (5S and 5.8S rRNAs, snoRNA and snRNAs).

Statistics and Reproducibility
No formal statistical method was used to predetermine the sample size. Sample sizes were determined based on prior studies to ensure reliable measurement of experimental parameters according to field standards39, 41, 43-45, 81. Samples were randomly allocated, and blinding was implemented during library construction. A total of 7 (PrF), 12 (PF), 4 (GV), and 7 (MII) oocytes were collected at each stage, of which 5 (PrF), 9 (PF), 4 (GV), and 5 (MII) were retained after quality filtering and included in downstream analyses. Data analyses were conducted using R version 3.6.3. Results are presented as mean ± standard deviation (s.d.), with statistical tests and p values provided in figure legends.

Data Availability
The sequencing data generated in this study have been deposited in the NCBI GEO database under accession number GSE296363 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE296363]. Publicly available datasets used in this study include small-RNA sequencing data for golden hamsters from GSE16952843 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE169528] and cynomolgus monkey MII oocytes from GSE9521841 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE95218].

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The Piwi-piRNA pathway is an evolutionarily conserved germline specific repressor of transposable elements (TEs). In their manuscript \"Integrated small and long RNA sequencing in single oocytes reveals piRNA-mediated transposon repression during human oogenesis,\" Zhang et al., extending their previous analysis of human piRNAs (Ref 41), characterized piRNAs during human oocyte maturation at single cell levels. The results suggest a model in which short-piRNAs (17-24nt) act as the primary and broad-spectrum TE suppressors and long-piRNAs (25-35nt) provide coordinated endogenous retrovirus (ERV)-specific silencing. Overall, the results are extensively presented and the main conclusions are supported by the data. Studies of this sort are good references and resources for further comparisons, However, once the points below are addressed I feel that the manuscript may be suitable for publication in Nature Communications.\nBelow we list several areas that in our opinion warrant improvement. We also suggest additional analyses that, if fruitful, may elevate the biological significance of the study.\n1. By definition, piRNAs are small RNAs loaded onto PIWI proteins. They previously examined PIWI-associated small RNAs (Ref 41). However, in this study, they sequenced total small RNAs. How do they know “short-piRNAs” are indeed piRNAs and associated with PIWIL3? How many of piRNAs identified here are identical to those identified in the previous study (Ref 41)? Do they exhibit no characteristics of endo-siRNAs with 2nt overhangs and complete complementarity?\n2. Figure 1f (line 136-144): The authors show that L1PA2 has undergone a 128.90-fold downregulation while HERV-L and HERV-K have exhibited a 89.43-fold and a 29.50-fold repression, respectively. The authors should show expression levels of corresponding piRNAs (just like Figure 4e). This will help the reader understand the relationship between expression levels of piRNAs and degree of TE repression.\n3. Figure 1f: Some group of ERVs are upregulated during oocyte development. What might distinguish this group of TEs from other TEs which are downregulated during oocyte development. Are these TEs located in piRNA clusters identified in this study?\n4. Line 166: “3’-O-methylation” should be “2’-O-methylation at the 3’ end.” The authors should describe methods used to examine the presence/absence of this modification. This reviewer cannot find methods used in Methods.\n5. Figure 4f: It is peculiar that the authors do not cite the paper of Ishino et al. (NAR 49, 2700-2720, 2021) and discuss the discrepancy between the two studies. Ishino et al. conclude that most of the piRNAs expressed in oocytes are not involved in the ping-pong cycle.\n6. Figure 4f: L1 subfamilies are largely divided into two subgroups in terms of the z-score of ping-pong signature (high vs low). Why? Is the low z-score due to the primary processing? Is there any peculiarity of piRNAs in this group?\n7. Figure 6a-d (line 311-313): The authors claim a strong bias toward L1 and ERV fragments within the top 20 clusters. However, 11 out of the top 20 clusters show L1 antisense bias, which does not seem statistically strong. Could these data imply that piRNA clusters are differentially expressed at each stage of oocyte development? Do piRNA cluster share common promoters/transcription factor binding sites in a stage-specific manner, for example?\n8. Figure 6e: How do the authors know that the piRNA cluster is unidirectionally transcribed? Are short-piRNAs produced by the ping-pong cycle derived from unidirectionally transcribed cluster transcripts?", "reviewer_2": "Integrated small and long RNA sequencing in single oocytes reveals piRNA-mediated transposon repression during human oogenesis\nThe manuscript describes the expression, by nucleotide sequencing, of small non coding RNAs, mRNAs and transposon elements in human oocytes of 4 different developmental stages: primordial follicle, primary follicle, germinal vesicle and metaphase II. It is concluded that short piwi interacting RNAs (piRNAs) act as primary suppressors of transposon activity, while long piRNAS coordinate silencing of specific ERV transposons.\nIn general, the manuscript provides a nice descriptive dataset of expressed piRNA and transposon sequences at different stage of human oocytes. In that respect it is of value to the scientific community.\nFunctional data that demonstrate the activity of certain piRNAs, or how possible transposon activity would affect the quality of the oocyte or a possible developing embryo are lacking however. This limits the conclusions that can be drawn from the data, and as such it is merely a descriptive manuscript. There may be correlation, but there is no demonstration of a causation.\nMajor comments.\nThe analysis of piRNA, PIWI protein genes and transposon expression of human oocytes at different developmental stages is solid, well described and well presented in clear figures. The sequencing, though solid and extensive, is not novel. Similar data have been published previously eg PMID: 34340992; PMID: 33622565; PMID: 25818294; PMID: 31358756.\nFunctional data are lacking in the manuscript. In that sense it is descriptive, and conclusions drawn are rather speculative. Different interpretations of the data are possible, and how presence of absence of piRNAs would affect the quality and functionality of oocytes remains unknown.\nAlthough in the introduction it is mentioned that in mice oocytes transposon activity is thought to be regulated via miRNAs and a specific Dicer isoform, crucial information here that the mouse genome does not harbour Piwil3 should be mentioned.\nIn the results section it is stated that 30 oocytes were collected. The value of this number is limited, as only 23 oocytes were used. The numbers is Fig1A ae misleading: PrF n=7: only 5 were used, PF n=12: only 9 were used, MII n=7: only 5 were used.\nIt is not clear whether the somatic cells surrounding the oocytes (granulosa and cumulus) were removed. For the GV and MII oocytes it is evident that these were removed but for the oocytes from primordial follicles this is less clear. How certain are the authors that all somatic cells were indeed removed? Indeed in figure 1a, the oocytes are presented with somatic cells present.\nFig 1c, and 2d, the colored areas are a bit misleading, just the colored dots would be more informative and realistic.\nIn the different violin plots such as Fig 1g, 2c, 2e, 4b, 4d,5f, 5g indications of which data sets differ statistically from others is essential.\nThe use of the word ‘repression’ is suggestive, or misleading, as it indicates and active mechanism. For instance in line 150 where it is stated that Alu and SVA elements underwent repression with fold changes of 4.71 and 5.36. What is the indication that these elements were repressed, rather than having reduced expression? From these data there is no indication of repression. Similar for ‘robust silencing’ in line 149. There are no experimental data indicating that the expression is actively silenced, instead of expression simply being at a lower level.\nThe manuscript is concentrated on indicating that piRNAs are important for reducing transposon activity. At the same time, expression of miRNAs, even at low levels, is suggested to ‘ contribute to oocyte development’ line 188. Although it is understandable that piRNAs are linked with suppression of TE activity, the data presented in this manuscript do not exclude that piRNA, similar to miRNAs, have another role besides repression of TE activity that is beneficial for acquiring oocyte developmental competence.\nOn the same note, expression of several TE increase during oocyte development (line 148; fig 1f). What is the proposed function of this increased expression. Could it be important for acquisition of developmental competence?\nFig 2g is difficult to read, individual lines/arrows cannot be distinguished. As such the figure is not informative.\nFigure 3c, ‘ Monkey’ is a general term, it should be further specified; eg Macaca fascicularis, or whichever species was used here.\nLine 267 The Ping-Pong signature was either undetectable or 267 substantially diminished for most TE-aligned long-piRNAs. How does this vibe with other publications such as PMID: 25818294?\nThe conclusions are too strong. For instance line 352, a role for piRAs has not been provided, but suggested. By examining piRNA and TE expression levels, it is by no means demonstrated that these affect each other. Line 352, the results do not ‘emphasize the suppressive effects of piRNAs on TE’, but instead provide a correlation between piRNAs and TEs.\nLine 364-371: It should be made clear that (control of) TE activity in germ cells is different from TE activity in somatic cells.\nMinor\nLine 34 ‘(long-piRNAs)’ can be removed.", "reviewer_3": "Please see uploaded PDF (to preserve formatting)."}, "author_response": "We sincerely thank the reviewer for the thorough and constructive evaluation of our manuscript. We are encouraged that the reviewer found our results comprehensive, our main conclusions well supported by the data, and our study potentially valuable as a reference for future research. We also greatly appreciate the reviewer’s recognition of the novelty of our analyses and of the biological model we propose for the roles of short- and long-piRNAs in TE regulation in human developing oocytes.\nWe further thank the reviewer for the insightful comments and professional suggestions, which have substantially strengthened the biological significance and overall quality of our study after incorporation. We have carefully revised the manuscript in response to all points. Below, we provide a detailed, point-by-point response, highlighting the corresponding revisions made to the text, figures, and supplementary materials where appropriate.\n1. By definition, piRNAs are small RNAs loaded onto PIWI proteins. They previously examined PIWI-associated small RNAs (Ref 41). However, in this study, they sequenced total small RNAs. How do they know “short-piRNAs” are indeed piRNAs and associated with PIWIL3? How many of piRNAs identified here are identical to those identified in the previous study? Do they exhibit no characteristics of endo-siRNAs with 2nt overhangs and complete complementarity?\nResponse: We thank the reviewer for raising these important points. In our previous work, we identified that ~19-nt piRNAs associate with PIWIL3 and lack 2’-O-methylation at their 3’ ends (1, 2). In the revised manuscript, we performed PIWIL3 immunoprecipitation (PIWIL3-IP) from human primordial and primary follicle oocytes, followed by sequencing of the associated small RNAs. This analysis confirmed that PIWIL3 specifically binds short-piRNAs with a size peak at 17–24 nt (Fig. 2d). Furthermore, we verified that most of short-piRNAs lack 2′-O-methylation at their 3′ ends, as their signals were largely eliminated after oxidation treatment (Fig. 2e). In contrast, long-piRNAs (25-35 nt) were resistant to oxidation (Fig. 2e), indicating 3’ terminal 2’-O-methylation. These biochemical characteristics of short- and long-piRNAs are consistent with observations in human GV and MII oocytes (1, 2). The new analyses have been incorporated into (Fig. 2d, e) and described in the corresponding Results section (Line 171-179) of the revised manuscript. The corresponding experimental details have been provided in the Materials and Methods section (Lines 528-551).\nFigure 2. piRNAs are the predominant small non-coding RNAs during human oogenesis.\n(d) Length distribution of PIWIL3-associated small RNAs and total small RNAs from input oocytes. The PIWIL3-bound fraction is enriched for short-piRNAs (17-24 nt), whereas input samples contain both short (17-24 nt) and long-piRNAs (25-35 nt). Approximately 2.5 oocytes were pooled per sample for the IP assay.\n(e) Length profiles of small RNAs before and after oxidation treatment in PrF and PF oocytes. Short-piRNAs (17-24 nt) were largely depleted upon oxidation, whereas long-piRNAs (25-35 nt) were largely retained, confirming the absence of 2’-O-methylation in short-piRNAs and its presence in long-piRNAs. Approximately 2.5 oocytes were pooled per sample for each oxidation assay.\nNext, we compared the short-piRNAs identified in this study with previously reported PIWIL3-IP and oxidation-treated libraries from human GV and MII oocytes (1, 2). We found that approximately 50% of short-piRNAs overlapped with PIWIL3-associated sequences. Consistently, short-piRNAs were more abundant in non-oxidized than in oxidized libraries (Rebuttal Figure 1).\nRebuttal Figure 1. Comparison of short-piRNAs in this study with previously reported PIWIL3-IP and oxidation-treated libraries from human GV and MII oocytes.\n(a) Proportion of short-piRNAs identified in GV and MII oocytes from this study that overlap with PIWIL3-associated sequences.\n(b) Comparison of short-piRNAs identified in GV and MII oocytes from this study with previously reported oxidized and non-oxidized small RNA libraries. In both (a) and (b), PIWIL3-IP, oxidation, and non-oxidation treatments were performed using pooled GV and MII oocytes (1, 2). Statistical significance was evaluated using the Wilcoxon test.\nNext, we examined the sequence characteristics of these short-piRNAs. To validate the accuracy of our analytical pipeline for detecting 2-nt 3’ overhangs, we first applied it to analyze endo-siRNAs in mouse oocytes. As expected, we observed a clear signature of 2-nt 3′ overhangs (overhang length = 2) and perfect base-pairing (overhang length = 0), characteristic of endo-siRNAs (Rebuttal Figure. 2a). Applying the same pipeline to short-piRNAs identified in this study revealed no obvious 2-nt 3′ overhangs or perfect complementarity signature (Rebuttal Figure. 2b). Instead, these short-piRNAs exhibited canonical piRNA features, including a strong 1U and 10A bias (Supplementary Fig. 4c), a 10-nt Ping-Pong signature (Supplementary Fig. 4d).\nRebuttal Figure 2. Analysis of 3′ overhang signature for endo-siRNAs in mouse MII oocytes and short-piRNAs in human oocytes.\n(a) Distribution of endo-siRNAs (20-24 nt) with varying 3’-end overhang lengths in mouse MII oocytes. Data represent pooled endo-siRNAs from 5 oocyte samples.\n(b) Distribution of short-piRNAs (17-24 nt) with varying 3’-end overhang lengths in PrF, PF, GV, and MII oocytes. Data represent pooled short-piRNAs derived from 4-9 oocytes across different developmental stages. In both (a) and (b), Z-scores for overhang lengths of 0 nt (Z0) and 2 nt (Z2) are shown and calculated with overhang lengths of 1 and 3-20 nt used as background.\nTogether, these results demonstrate that short-piRNAs characterized in our study are bona fide PIWIL3-associated piRNAs that lack 2′-O-methylation at their 3’ ends across multiple oocyte stages.\n2. Figure 1f (line 136-144): The authors show that L1PA2 has undergone a 128.90-fold downregulation while HERV-L and HERV-K have exhibited a 89.43-fold and a 29.50-fold repression, respectively. The authors should show expression levels of corresponding piRNAs (just like Figure 4e). This will help the reader understand the relationship between expression levels of piRNAs and degree of TE repression.\nResponse: We analyzed the relationship between TE expression changes and the corresponding levels of short- and long-piRNAs, with strand specificity, focusing on L1PA2, HERV-L, and HERV-K. Our analysis revealed that the downregulation of L1PA2 was significantly correlated with an increase in both sense and antisense short-piRNAs (Supplementary Fig. 9b). In contrast, L1PA2-related sense and antisense long-piRNAs decreased during human oocyte development. For HERV-L, its downregulation was accompanied by an increase in both short- and long-piRNAs, independent of strand orientation (Supplementary Fig. 9c), suggesting a coordinated role for both piRNA classes in repressing this element. For HERV-K, only sense short-piRNAs showed a marked increase from the PrF to the MII stage, and this increase was significantly correlated with reduced HERV-K transcript abundance (Supplementary Fig. 9d). In addition, we examined the expression changes of L1HS, AluYb8, and SVA_A, along with their corresponding short- and long-piRNAs (Supplementary Fig. 9a, e, and f).\nCollectively, these results support the conclusion that short- and long-piRNAs may play distinct, TE subfamily-specific regulatory roles during human oocyte maturation. These new analyses have been incorporated into the revised Supplementary Fig. 9 and are described in the corresponding Results section (Line 316-329) of the revised manuscript.\nSupplementary Figure 9. Relationship between piRNA abundance and the extent of TE downregulation.\n(a-f) Scatter plots showing the correlations between TE expression levels and the abundance of their corresponding short- and long-piRNA in developing human oocytes for L1HS (a), L1PA2 (b), HERVL (c), HERVK (d), AluYb8 (e), and SVA A (f). Pearson correlation coefficients (R) and p-values (p) are shown in each panel.\n3. Figure 1f: Some group of ERVs are upregulated during oocyte development. What might distinguish this group of TEs from other TEs which are downregulated during oocyte development. Are these TEs located in piRNA clusters identified in this study?\nResponse: We thank the reviewer for this insightful question. We examined multiple characteristics of upregulated and downregulated TE subfamilies, including related short- and long-piRNA abundance, Ping-Pong Z-scores, 1U and 10A ratios, and TE transcript levels. We found that upregulated TE subfamilies displayed significantly higher transcript levels than downregulated TE subfamilies (Supplementary Fig. 10a). However, apart from this higher-level in transcript abundance, all other piRNA-related signatures were markedly lower in upregulated TE subfamilies—most notably, the abundance of corresponding antisense short- and long-piRNAs. Consistently, upregulated TE subfamilies exhibited weaker Ping-Pong signatures, reflected by reduced 10-nt overlap ratios and 1U/10A biases.\nSupplementary Figure 10. Comparative characteristics of up- and down-regulated TE subfamilies.\n(a) Comparison of up- and down-regulated TE subfamilies based on related piRNA abundance, Ping-Pong Z-scores, 1U and 10A ratios, and TE transcript levels.\n(b) Boxplots showing the density of upregulated and downregulated TE insertions within the top 20 piRNA clusters, separated by strand orientation (Antisense vs. Sense). For both (a) and (b), statistical significance was calculated using a t-test.\nWe next analyzed the genomic distribution of upregulated and downregulated TE subfamilies within the top 20 piRNA clusters, which collectively generate over 90% of uniquely mapped piRNAs. Although upregulated TE subfamilies were present within these clusters, their antisense insertions were significantly less frequent than those of the downregulated TE subfamilies, consistent with their reduced antisense short- and long-piRNA levels (Supplementary Fig. 10b). Notably, seven out of the eight upregulated TE subfamilies are human-specific elements. Given their relatively high transcript abundance and reduced antisense piRNA expression, we propose that these upregulated elements likely represent evolutionarily younger TE subfamilies for which the piRNA/PIWI pathway has not yet fully adapted to achieve efficient repression—potentially due to insufficient antisense insertions within highly productive piRNA clusters.\nThese new analyses have been incorporated into the revised Supplementary Fig. 10 and are described in the corresponding Results section (Line 333-341) of the revised manuscript, which reads as:\n“Compared with TE subfamilies that were downregulated, these upregulated TE subfamilies showed higher transcript abundance but produced markedly fewer corresponding antisense and sense short- and long-piRNAs, accompanied by weaker Ping-Pong signatures and reduced 1U and 10A biases (Supplementary Fig. 10a). Genomic distribution analyses within the top 20 piRNA clusters further revealed substantially fewer antisense insertions for these upregulated TE subfamilies relative to downregulated TE subfamilies, consistent with their reduced antisense piRNA levels (Supplementary Fig. 10b)”.\n4. Line 166: “3’-O-methylation” should be “2’-O-methylation at the 3’ end.” The authors should describe methods used to examine the presence/absence of this modification. This reviewer cannot find methods used in Methods.\nResponse: We thank the reviewer for pointing out this error. It has been corrected to “2’-O-methylation at the 3’ end” in the revised manuscript.\nIn the original version, we did not perform oxidation treatment followed by small RNA sequencing to assess the presence or absence of 2’-O-methylation, as our previous studies had already shown that short-piRNAs in GV and MII oocytes lack this modification, whereas long-piRNAs possess it (1, 2). The initial description was therefore based on these published findings. Moreover, oocytes at the early follicular stages yield extremely limited amounts of RNA, and human oocytes are particularly precious, making oxidation-based library construction technically challenging. Nevertheless, to further strengthen the rigor of the present study, we collected oocytes at the primordial and primary follicular stages and performed oxidation sequencing. Comparison between oxidized and non-oxidized samples revealed that short-piRNAs were largely depleted upon oxidation, whereas long-piRNAs were largely retained, particularly in PF-stage oocytes (Fig. 2e). These data confirm that, in oocytes at the primordial and primary follicular stages, PIWIL3-associated short-piRNAs lack 2′-O-methylation at their 3′ ends, whereas long-piRNAs possess this modification, consistent with our previous observations in GV and MII oocytes (1, 2). The incomplete depletion of short-piRNAs is likely due to the extremely limited input RNAs (approximately 2.5 oocytes per reaction), which may reduce the efficiency of the oxidation treatment.\nFigure 2. piRNAs are the predominant small non-coding RNAs during human oogenesis.\n(e) Length profiles of small RNAs before and after oxidation treatment in primordial follicle oocytes. Short-piRNAs (17-24 nt) were largely depleted upon oxidation, whereas long-piRNAs (25-35 nt) were largely retained, confirming the absence of 2’-O-methylation in short-piRNAs and its presence in long-piRNAs. Approximately 2.5 oocytes were pooled per sample for each oxidation assay.\nThese new results have been included in the revised Figure 2e and described in the Results section (Line 173-179) in the revised manuscript. The corresponding experimental details have been provided in the Materials and Methods section (Lines 541-551).\n5. Figure 4f: It is peculiar that the authors do not cite the paper of Ishino et al. (NAR 49, 2700-2720, 2021) and discuss the discrepancy between the two studies. Ishino et al. conclude that most of the piRNAs expressed in oocytes are not involved in the ping-pong cycle.\nResponse: We thank the reviewer for pointing out this important reference. We note a key methodological difference between the two studies. Ishino et al. analyzed all piRNA sequences associated with PIWIL1 and PIWIL3 in hamster testis, ovary, MII oocytes, and 2-cell embryos, and reported no appreciable A enrichment at the 10th position or 10-nt overlaps between sense and antisense piRNAs. Based on these genome-wide analyses, they concluded that PIWIL1- and PIWIL3-bound piRNAs exhibit minimal or no Ping-Pong activity. In contrast, our study specifically focused on TE-related piRNAs. Because most genomic TE copies are truncated or inactive, we mapped piRNAs to consensus TE sequences downloaded from Repbase, which represent reconstructed ancestral active elements. This approach enables a more targeted evaluation of potentially active TEs and their associated piRNA signatures, which may not be detectable when analyzing the whole piRNA population.\nTo directly address this discrepancy, we reanalyzed the data of Ishino et al. by using both strategies: (i) all piRNA sequences, (as in Ishino et al.), and (ii) piRNAs mapped to consensus TE sequences (as used in this study and our previous work) (1, 3, 4). Consistent with Ishino et al., the nucleotide bias analysis using all piRNA sequences showed very weak or no 10A enrichment (Rebuttal Figure 3a). In contrast, the analysis based on consensus TE sequences revealed a modest but detectable 10A bias across all examined samples (Rebuttal Figure 3b). Consistently, piRNAs mapped to the sense and antisense strands of consensus TEs exhibited an obvious 10-nt overlap, indicating the presence of a low-level Ping-Pong signature (Rebuttal Figure 3c).\nTogether, these results suggest that the apparent discrepancy between our findings and those of Ishino et al. primarily arises from differences in analytical strategy.\nRebuttal Figure 3. Ping-Pong signature analysis for the dataset of Ishino et al.\n(a) Nucleotide enrichment analysis of all piRNAs associated with PIWIL1 and PIWIL3 in hamster testis, ovary, MII oocytes, and 2-Cell embryonic stage.\n(b) Nucleotide enrichment analysis for consensus TE-aligned piRNAs associated with PIWIL1 and PIWIL3 in hamster testis, ovary, MII oocytes, and 2-Cell embryonic stage.\n(c) Ping-Pong overlap distributions of consensus TE-aligned sense and antisense piRNAs associated with PIWIL1 and PIWIL3 in hamster testis, ovary, MII oocytes, and 2-Cell embryonic stage. Z-scores for 10-nt overlaps (Z10) were calculated by using positions 0-9 and 11-20 as background and are shown in each panel.\n6. Figure 4f: L1 subfamilies are largely divided into two subgroups in terms of the z-score of ping-pong signature (high vs low). Why? Is the low z-score due to the primary processing? Is there any peculiarity of piRNAs in this group?\nResponse: To investigate the differences between L1 subfamilies with high (Group A) and low (Group B) Ping-Pong signatures (Supplementary Fig. 5a), we systematically analyzed nine features for Group A and B, including related short- and long-piRNA expression levels, overall 1U and 10A ratios, Ping-Pong Z-scores of short- and long-piRNAs, and the corresponding transcript abundances. L1 subfamilies with stronger Ping-Pong signatures (Group A) exhibited consistently higher values across all examined features than those with weaker Ping-Pong signatures (Group B) Supplementary Fig. 5b). As the reviewer noted, short-piRNAs in Group A displayed more pronounced Ping-Pong signatures than those in Group B, and a similar but weaker trend was also observed for long-piRNAs. This pattern was further supported by a higher 10A ratio observed in Group A. Moreover, the transcript abundance of L1 subfamilies in Group A, as well as the expression levels of their corresponding sense and antisense piRNAs, was significantly higher than in Group B.\nAll these observations suggest that the elevated levels of both antisense and sense piRNAs, together with the greater availability of TE transcripts as Ping-Pong substrates in Group A, promote more efficient Ping-Pong amplification.\nSupplementary Figure 5. Comparative features of L1 subgroups with high and low Ping-Pong signatures\n(a) Classification of L1 subfamilies into Group A (high Ping-Pong Z-scores) and Group B (low Ping-Pong Z-scores).\n(b) Comparison of nine features between Group A and Group B L1 subfamilies, including related short- and long-piRNA expression levels, overall 1U and 10A ratios, Ping-Pong Z-scores of short and long-piRNAs, and the corresponding transcript abundances.\n(c) Distribution of sense and antisense insertions of L1 subfamilies from Group A and B within the top 20 piRNA clusters. For both (b) and (c), statistical significance was calculated using a t-test.\nFinally, we examined the genomic distribution of L1 subfamilies from Group A and B within the top 20 piRNA clusters. Interestingly, Group B showed a higher proportion of antisense insertions compared to Group A, while sense insertions displayed no significant difference (Supplementary Fig. 5c). Given the lower levels of both antisense and sense piRNA in Group B, these results suggest that primary processing and subsequent PIWI loading may preferentially promote corresponding piRNA generation in Group A. However, we cannot exclude the possibility that a fraction of antisense and sense piRNAs in Group A also arises from Ping-Pong amplification.\nCollectively, these analyses indicate that the reduced Ping-Pong signature observed in certain L1 subfamilies is likely due to their lower transcript abundance and limited antisense and sense piRNA generation through primary processing.\nThese new analyses have been incorporated into the revised Supplementary Fig. 5 and are described in the corresponding Results section (Line 270-280) of the revised manuscript, which reads as:\n“Notably, approximately two-thirds of L1 subfamilies exhibited strong Ping-Pong activity (Group A) and consistently showed higher 1U and 10A proportions, greater antisense and sense piRNA abundance, and elevated transcript levels compared with those with low signals (Group B) (Supplementary Fig. 5a, b). Interestingly, genomic analysis of the top 20 piRNA clusters revealed that Group B L1 subfamilies carried more antisense insertions than Group A, while sense insertions were comparable (Supplementary Fig. 5c). Given the lower abundance of both antisense and sense piRNAs for Group B, these results suggest that primary processing and PIWI loading preferentially enhance piRNA biogenesis from Group A L1 subfamilies. These findings indicate that weak Ping-Pong activity in certain L1 subfamilies likely reflects their lower transcript abundance and inefficient primary piRNA processing.”\n7. Figure 6a-d (line 311-313): The authors claim a strong bias toward L1 and ERV fragments within the top 20 clusters. However, 11 out of the top 20 clusters show L1 antisense bias, which does not seem statistically strong. Could these data imply that piRNA clusters are differentially expressed at each stage of oocyte development? Do piRNA cluster share common promoters/transcription factor binding sites in a stage-specific manner, for example?\nResponse: We thank the reviewer for these insightful comments. As the reviewer noted, although L1 showed higher antisense piRNA production than sense piRNA production in the top 20 piRNA clusters, particularly within the top 10 clusters, they do not display a clear antisense insertion bias when directly compared with their sense insertions. We have therefore revised the description to more accurately reflect this point, as now stated:\n“Interestingly, although L1 elements did not show a pronounced difference in antisense versus sense insertion frequencies, the top 10 piRNA clusters still produced markedly more antisense than sense piRNAs (Fig. 6d and Supplementary Fig. 12a). This suggests selective primary processing, in which L1 antisense piRNAs may be preferentially generated or more efficiently loaded onto PIWI proteins.” (Lines 376–381).\nAs the reviewer noted, piRNA clusters are indeed differentially expressed during oocyte development, and their contributions to the TE-related piRNA population vary by stage. For example, the cluster chr4:78106722–78396420|+ contributed the second-largest number of L1-related antisense short-piRNAs at the PrF/PF stages but became the predominant source at the GV/MII stages. Similarly, the cluster chr12:114437233–114625671|− ranked second for ERV-related antisense short-piRNA production at the PrF stage but became the leading contributor after the PF stage. These examples have been added to the revised manuscript to further support our conclusion that piRNA clusters are differentially expressed and dynamically regulated across oocyte developmental stages, as now stated:\n“These clusters displayed dynamic expression patterns during oocyte development, with their contributions to TE-related piRNA population varying by stage (Fig. 6a, b, and Supplementary Data 9). For example, the cluster chr4:78106722–78396420|+ ranked second for L1 antisense short-piRNA generation at the PrF and PF stages, but became the predominant source at the GV/MII stages. Similarly, the cluster chr12:114437233–114625671|− produced the second-largest number of ERV antisense short piRNAs at the PrF stage and emerged as the primary contributor after the PF stage.” (Lines 356-363).\npiRNA clusters are generally defined as genomic regions enriched for piRNAs, although their primary transcripts are often difficult to delineate due to extensive precursor processing. Previous work (5) demonstrated that bidirectional piRNA clusters are particularly informative for identifying promoters and associated transcription factors. That study further showed that A-MYB acts as a key transcriptional regulator of pachytene piRNA biogenesis and PIWIL1 expression during mouse spermatogenesis, with a conserved role in other species.\nIn our analysis of the top 20 piRNA clusters, we identified 2 pairs of bidirectional clusters (chr5:73361914–73441831|– vs chr5:73438508–73512499|+, and chr6:39934078–40135208|– vs chr6:40129348–40169244|+) (Supplementary Fig. 11). Notably, one of these four loci (chr6:40129348–40169244|+) was ranked outside the top 20 based on relatively lower unique piRNA generation. Expanding the analysis to the top 100 piRNA clusters revealed two additional pairs of bidirectional clusters. These findings suggest that, unlike pachytene piRNAs in mouse testis, only a small subset of human oocyte piRNA clusters originates from bidirectionally transcribed loci, consistent with previous observations in golden hamster oocytes (6). Among the four bidirectional pairs identified within the top 100 piRNA clusters, two appear to arise from back-to-back transcribed primary transcripts (chr5:73361914–73441831|– vs chr5:73438508–73512499|+; chr6:39934078–40135208|– vs chr6:40129348–40169244|+), which may share common promoter or transcription factor binding regions. We analyzed these shared regions to predict potential motifs or regulatory factors; however, no significantly enriched motifs were detected, likely due to the limited sequence regions available for analysis.\nWe have added the genome browser views for these two back-to-back transcribed piRNA clusters in revised Supplementary Fig. 11, and a description of this observation into the revised manuscript:\n“Notably, only two pairs of bidirectionally transcribed clusters were identified among the top 20 piRNA clusters (Supplementary Fig. 11), and expanding the analysis to the top 100 clusters revealed two additional pairs. This indicates that, in contrast to pachytene piRNAs in mouse testis (5), only a small fraction of human oocyte piRNA clusters arises from bidirectionally transcribed loci, consistent with observations in golden hamster oocytes (6).” (Lines 363–369).\nSupplementary Figure 11. UCSC genome browser views of two back-to-back transcribed piRNA clusters\n(a) Genome browser view of back-to-back transcribed piRNA cluster pair (chr5:73361914 - 73441831|– vs chr5:73438508 - 73512499|+), showing the surrounding gene annotations, density of TE superfamily insertions, long RNA-seq signal, and expression levels of uniquely and total mapped short- and long-piRNAs.\n(b) Genome browser view of back-to-back transcribed piRNA cluster pair (chr6:39934078–40135208|– vs chr6:40129348–40169244|+), showing the surrounding gene annotations, density of TE superfamily insertions, long RNA-seq signal, and expression levels of uniquely and total mapped short- and long-piRNAs.\nStrand-specific piRNA read coverage on opposing strands reveals convergent transcription and reciprocal processing, characteristic of bidirectional piRNA cluster organization. Long RNA-seq signal is strand-unspecific.\n8. Figure 6e: How do the authors know that the piRNA cluster is unidirectionally transcribed? Are short-piRNAs produced by the ping-pong cycle derived from unidirectionally transcribed cluster transcripts?\nResponse: In the current version of the manuscript, we did not specifically discuss the transcriptional orientation of piRNA clusters following their identification. Our identification strategy was based on detecting regions enriched for piRNAs mapped to the same strand, and regions containing more than five unique piRNAs were defined as piRNA clusters.\nAs described in our response to Question 7, we identified 3 bidirectional piRNA clusters among the top 20 clusters, with the remaining 17 classified as unidirectional. We next examined the origin of short-piRNAs involved in the Ping-Pong amplification cycle relative to bidirectional and unidirectional cluster types. Our analysis showed that over 75% of TE-related short-piRNAs exhibiting a characteristic 10-nt overlap were derived from unidirectional clusters, whereas fewer than 25% originated from the bidirectional clusters (Rebuttal Figure. 4a). A similar distribution was observed for long-piRNAs, indicating that both short- and long-piRNAs with Ping-Pong signatures are predominantly generated from unidirectional clusters (Rebuttal Figure. 4b). It is important to note that many piRNAs with 10-nt overlaps are multi-mapped; therefore, all mapped piRNAs were included in this analysis, which means may lead to partial overlap between piRNAs assigned to unidirectional versus bidirectional clusters.\nRebuttal Figure 4. Proportion of TE-related short- (a) and long-piRNAs (b) with 10-nt overlap signatures derived from the top 20 unidirectional and bidirectional piRNA clusters. All piRNAs (uniquely and multiply mapped) were included in the analysis.\nReviewer #2 (Remarks to the Author):\nIntegrated small and long RNA sequencing in single oocytes reveals piRNA-mediated transposon repression during human oogenesis\nThe manuscript describes the expression, by nucleotide sequencing, of small non coding RNAs, mRNAs and transposon elements in human oocytes of 4 different developmental stages: primordial follicle, primary follicle, germinal vesicle and metaphase II. It is concluded that short piwi interacting RNAs (piRNAs) act as primary suppressors of transposon activity, while long piRNAS coordinate silencing of specific ERV transposons.\nIn general, the manuscript provides a nice descriptive dataset of expressed piRNA and transposon sequences at different stage of human oocytes. In that respect it is of value to the scientific community.\nResponse: We sincerely thank the reviewer for the thoughtful and constructive evaluation of our manuscript. We are encouraged that the reviewer considers our dataset valuable to the scientific community and appreciates the comprehensive characterization of piRNAs and transposon dynamics across human oocyte development. We also thank the reviewer for the insightful comments and suggestions regarding additional analyses and clarifications that could enhance the biological significance of our study. We have carefully revised the manuscript in response to these points. Below, we provide a point-by-point response, highlighting the changes made to the text, figures, and supplementary materials where appropriate.\nFunctional data that demonstrate the activity of certain piRNAs, or how possible transposon activity would affect the quality of the oocyte or a possible developing embryo are lacking however. This limits the conclusions that can be drawn from the data, and as such it is merely a descriptive manuscript. There may be correlation, but there is no demonstration of a causation.\nResponse: We thank the reviewer for highlighting the important distinction between correlation and causation, and for pointing out the limitations arising from the absence of direct functional perturbation experiments. We agree that functional assays would provide definitive evidence for the mechanistic roles of short- and long-piRNAs in transposon regulation. However, such experiments remain technically and ethically challenging in human oocytes. Given these constraints, our study aims to establish a foundational, high-resolution atlas of piRNA populations and TE expression, and to propose a potential regulatory model supported by multiple, independent lines of correlative evidence. We have now revised the Discussion to explicitly acknowledge this limitation and to emphasize the need for further functional, clinical, or genetic studies to further elucidate the consequences of PIWI/piRNA pathway activity on TE regulation, as stated:\n“Second, definitive functional validation will be required to establish the causal roles of piRNAs in human oocyte development. However, such experiments remain largely infeasible due to the extremely limited availability of human oocytes and strict ethical constraints. Even so, evidence from model organisms shows that disruption of piRNA biogenesis leads to defective oocyte development and early embryonic arrest (7-11), underscoring the central importance of this pathway. Future clinical and genetic studies will therefore be critical to elucidate the functional consequences of PIWI/piRNA dysfunction in human reproductive biology.” (Lines 475–482).\nMajor comments.\nThe analysis of piRNA, PIWI protein genes and transposon expression of human oocytes at different developmental stages is solid, well described and well presented in clear figures. The sequencing, though solid and extensive, is not novel. Similar data have been published previously eg PMID: 34340992; PMID: 33622565; PMID: 25818294; PMID: 31358756.\nResponse:\nWe thank the reviewer for the positive comments on our data analysis and for noting that previous studies have profiled small RNAs during human oogenesis and early embryonic development. However, we believe that our study provides at least two important advances beyond these earlier efforts.\n1. Integrated profiling of small and long RNAs simultaneously at the single-oocyte level. Our dataset simultaneously profiles both small and long RNAs from the same individual human oocytes, allowing integrated analyses of small regulatory RNAs, including piRNAs and miRNAs, together with long transcripts, such as protein-coding mRNAs and transposable element transcripts. This approach enables a more accurate and comprehensive understanding of the relationships between small RNAs and their targets at single-oocyte resolution.\n2. Comprehensive developmental coverage, particularly follicular stages. Our dataset encompasses four key stages of human oogenesis—primordial follicle, primary follicle, germinal vesicle, and metaphase II—providing a broader temporal view of both small and long RNA dynamics. As the reviewer may note, small RNA datasets for human oocytes at the follicular stages have been largely unavailable. Characterizing small RNAs during these early developmental stages is crucial for understanding transcriptional and post-transcriptional regulation underlying follicle activation and growth, offering valuable insights into follicle recruitment and the restoration of oocyte competence. The lack of such datasets likely reflects the extreme scarcity of early-stage human oocytes and the low abundance of small RNAs at these stages, which require highly sensitive library construction methods. To overcome these challenges, in this study, we optimized our small RNA library construction protocol to make it suitable for early-stage human oocytes.\nTogether, these features make our dataset a unique and valuable resource for studying human oogenesis, providing integrated, stage-resolved, single-oocyte–level insights into small RNA functions and their potential regulatory roles in gene expression and transposable element control. We have revised the Discussion section to emphasize the novelty and uniqueness of our dataset, highlighting its advantage compared with previously published human oocyte datasets. The added text now reads:\n“Unlike previous studies that examined only small RNAs during human oogenesis and early embryonic development (2, 12, 13), our approach integrates matched small- and long-RNA sequencing from the same individual oocytes, including the previously underrepresented PrF and PF stages. This design provides a unified view of RNA regulation—from small RNA biogenesis and TE transcripts or mRNA expression dynamics—providing unprecedented resolution for investigating stage-specific small RNA-mediated regulatory processes in human oocytes.” (Line 420-426).\nFunctional data are lacking in the manuscript. In that sense it is descriptive, and conclusions drawn are rather speculative. Different interpretations of the data are possible, and how presence of absence of piRNAs would affect the quality and functionality of oocytes remains unknown.\nResponse: We agree that direct functional validation would provide the strongest evidence for the role of piRNAs in human oogenesis. However, such experiments are currently infeasible due to the extremely limited availability of human oocytes, particularly those with piRNA deficiency. Notably, studies in the model organisms have provided compelling evidence that disruption of the piRNA pathway leads to abnormal oocyte development and early embryonic arrest. We and two other groups independently established golden hamsters as a mammalian model for piRNA-related research and demonstrated that oocytes lacking piRNAs fail to support embryogenesis (4, 8, 14). Specifically, disruption of PIWIL1-associated long-piRNAs resulted in developmental arrest at the 2-cell stage, whereas loss of PIWIL3-associated short-piRNAs caused partial arrest between the 2- and 8-cell stages. More recently, a study from Dr. Haifan Lin’s group, using rabbits as an animal model, showed that depletion of PIWIL3-associated short-piRNAs leads to complete infertility, accompanied by severe defects in both oocyte maturation and early embryonic development (11). Together, these findings support a conserved and essential role for piRNAs in mammalian oogenesis and embryogenesis.\nAlthough direct functional evidence in human oocytes remains unavailable, these results from model organisms strongly suggest that piRNAs play a crucial role in maintaining oocyte quality and developmental competence. We have revised the Discussion section to acknowledge this limitation and clarify that our study provides correlative rather than functional evidence for the role of piRNAs in human oogenesis. The revised text reads:\n“Second, definitive functional validation will be required to establish the causal roles of piRNAs in human oocyte development. However, such experiments remain largely infeasible due to the extremely limited availability of human oocytes and strict ethical constraints. Even so, evidence from model organisms shows that disruption of piRNA biogenesis leads to defective oocyte development and early embryonic arrest (7-11), underscoring the central importance of this pathway. Future clinical and genetic studies will therefore be critical to elucidate the functional consequences of PIWI/piRNA dysfunction in human reproductive biology.” (Lines 475–482).\nAlthough in the introduction it is mentioned that in mice oocytes transposon activity is thought to be regulated via miRNAs and a specific Dicer isoform, crucial information here that the mouse genome does not harbour Piwil3 should be mentioned.\nResponse: We thank the reviewer for this helpful suggestion. We have revised the Introduction section to explicitly note that the mouse genome lacks the Piwil3 gene. The revised text now reads:\n“Interestingly, the Piwil3 gene is absent from the mouse genome and mouse oocytes primarily express endo-siRNAs rather than piRNAs” (lines 83-85)\nIn the results section it is stated that 30 oocytes were collected. The value of this number is limited, as only 23 oocytes were used. The numbers is Fig1A ae misleading: PrF n=7: only 5 were used, PF n=12: only 9 were used, MII n=7: only 5 were used.\nResponse: We thank the reviewer for pointing out this issue. In total, 30 oocytes were collected and sequenced, including 7 primordial follicle (PrF), 12 primary follicle (PF), 4 germinal vesicle (GV), and 7 metaphase II (MII) stage oocytes. However, simultaneous profiling of small and long RNAs from the same single oocyte remains technically challenging, and not all paired libraries passed our quality control (QC) criteria. To ensure high-quality data for downstream analysis, we applied stringent QC filters to both the small RNA and long RNA datasets.\nFor small RNA libraries, we excluded samples with (i) a useful read rate below 60% after 3′ adaptor trimming and short-read removal, (ii) a mapped read rate below 40%, and (iii) rRNA/tRNA-derived reads exceeding 30%. For long RNA libraries, we excluded samples with (i) a mapped read rate below 80%, (ii) mitochondrial-derived reads exceeding 10%, or (iii) rRNA-derived reads exceeding 10%. After applying these QC criteria, 23 oocytes passed QC for both RNA types and were retained for downstream analyses. We have added these QC details to the revised manuscript to clarify why 23 oocytes remained, as stated:\n“To ensure high data quality for downstream analysis, we applied stringent quality-control criteria to both the small RNA and long RNA datasets. For small RNA data, we removed samples with a useful read rate below 60% after 3’ adaptor trimming and short-read filtering, a mapped read rate below 40%, or rRNA/tRNA-derived reads exceeding 30%. For long RNA data, we excluded samples with a mapped read rate below 80%, mitochondrial-derived reads exceeding 10%, or rRNA-derived reads exceeding 10%.” (Lines 108-113).\nAlthough the total number of oocytes is limited, we incorporated statistical testing across all relevant analyses following the reviewer’s suggestion. Nearly all results supporting our main conclusions reached statistical significance, indicating that the number of oocytes analyzed is sufficient to ensure the robustness of our key findings. To avoid further misunderstanding, we have updated Fig. 1a and added the number of oocytes used for library construction as well as those retained after quality filtering in the figure legend.\nFigure 1. Global downregulation of TEs occurs after the primary follicle stages during human oocyte development and maturation\n(a) Schematic of the experimental strategy for simultaneous profiling of small and long RNA transcriptomes from single human oocytes at four developmental stages: primordial follicle (PrF), primary follicle (PF), germinal vesicle (GV), and metaphase II (MII). A total of 7 (Prf), 12 (PF), 4 (GV), and 7 (MII) oocytes were collected at each respective stage, of which 5 (Prf), 9 (PF), 4 (GV), and 5 (MII) were retained after quality filtering.\nIt is not clear whether the somatic cells surrounding the oocytes (granulosa and cumulus) were removed. For the GV and MII oocytes it is evident that these were removed but for the oocytes from primordial follicles this is less clear. How certain are the authors that all somatic cells were indeed removed? Indeed, in figure 1a, the oocytes are presented with somatic cells present.\nResponse: We thank the reviewer for raising this important point. All oocytes used for sequencing were thoroughly denuded of surrounding somatic cells, regardless of their developmental stage. We acknowledge that the previous description in the Methods section may have inadvertently suggested that this procedure was applied only to GV and MII oocytes. To avoid potential misunderstanding, we have now revised the text to explicitly state that the same denudation procedure was used for all follicular stages:\n“For all oocytes at all developmental stages, surrounding granulosa and cumulus cells were removed via incubation with 1 mg/mL hyaluronidase.” (Line 504-505).\nIn addition, Fig. 1a has been updated to clearly indicate that granulosa and cumulus cells were removed from all oocytes prior to sequencing.\nTo further confirm the absence of somatic cell contamination, we now provide two lines of evidence:\n1. Morphological validation: We provided representative bright-field images of oocytes before and after removal of surrounding somatic cells (Rebuttal Figure 5).\nRebuttal Figure 5. Morphology of human oocytes before and after somatic cell removal.\nRepresentative bright-field images of human oocytes before (left), during (middle), and after (right) removal of granulosa and cumulus cells. These images demonstrate the clean isolation of denuded oocytes used for downstream small RNA and long RNA sequencing. In the middle panel, black arrows indicate removed somatic cells, and white arrows indicate the denuded oocytes. Scale bars, 50 μm.\n2. Molecular validation: We additionally performed small RNA sequencing on (i) cumulus cell (CC) samples, and (ii) GV-stage oocytes surrounded by cumulus cells (GV+CC samples), which were obtained by gently needle flushing to remove most, but not all cumulus cells to avoid oocyte damage. We firstly compared the miRNA profiles between denuded GV oocytes and CC samples and identified both GV oocyte-specific and CC-specific expressed miRNAs (Rebuttal Figure. 6a). GV oocyte-specific expressed miRNAs such as miR-372, miR-516a-5p, and miR-5683 were detected in both denuded GV oocytes and GV+CC samples but absent in CC-only samples (Rebuttal Figure. 6b). Conversely, cumulus cell-specific miRNAs (miR-144-3p, miR-379-5p and miR-382-5p) were detected in CC-only and GV+CC samples but not in denuded GV oocytes and were also absent across all other oocyte stages analyzed in this study (Rebuttal Figure. 6b).\nTogether, these results confirm that somatic cell contamination in the oocytes used for sequencing was minimal.\nRebuttal Figure 6. Molecular validation of oocyte denudation.\n(a) Differential expressed analysis for miRNAs between cumulus cells (CC) and denuded GV stage oocytes.\n(b) Expression levels of cumulus cell–specific miRNAs (miR-144-3p, miR-379-5p, and miR-382-5p) in CC samples, GV oocytes surrounded by cumulus cells (GV+CC), denuded GV oocytes, and other stage oocytes (PrF, PF, and MII).\n(c) Expression levels of oocyte-enriched miRNAs (miR-372, miR-516a-5p, miR-5683) across CC samples, GV+CC samples, denuded GV oocytes, and oocytes at other developmental stages (PrF, PF, and MII). Significance was calculated by t-test and is shown in each panel.\nFig 1c, and 2d, the colored areas are a bit misleading, just the colored dots would be more informative and realistic.\nResponse: We have updated all PCA plots in Fig. 1b-c, Fig. 2b, Fig. 2f, Fig. 4a, and Fig. 5a by removing the colored background areas while retaining the elliptical outlines to improve visual clarity and facilitate interpretation. The ellipses were generated using the stat_ellipse function in the ggplot2 R package, which provides a statistically robust representation of data dispersion.\nRebuttal Figure 7. Updated PCA plots in Fig. 1b-c, Fig. 2b, Fig. 2f, Fig. 4a and Fig. 5a\nEllipses were generated using the stat_ellipse function in the ggplot2 R package, which provides a statistically supported visualization of group dispersion.\nIn the different violin plots such as Fig 1g, 2c, 2e, 4b, 4d,5f, 5g indications of which data sets differ statistically from others is essential.\nResponse: We thank the reviewer for this suggestion. Statistical details have been incorporated for all relevant figure panels.\nThe use of the word ‘repression’ is suggestive, or misleading, as it indicates and active mechanism. For instance in line 150 where it is stated that Alu and SVA elements underwent repression with fold changes of 4.71 and 5.36. What is the indication that these elements were repressed, rather than having reduced expression? From these data there is no indication of repression. Similar for ‘robust silencing’ in line 149. There are no experimental data indicating that the expression is actively silenced, instead of expression simply being at a lower level.\nResponse: We thank the reviewer for this comment and agree that the terms “repression” and “robust silencing” may not be fully supported by our current assays. Accordingly, we have replaced them with “reduced expression” or “downregulation,” which more accurately describe the observed decreases in TE expression levels without implying direct transcriptional or post-transcriptional silencing. We have carefully reviewed and revised the manuscript to ensure the consistent use of accurate and neutral terminology.\nThe manuscript is concentrated on indicating that piRNAs are important for reducing transposon activity. At the same time, expression of miRNAs, even at low levels, is suggested to ‘ contribute to oocyte development’ line 188. Although it is understandable that piRNAs are linked with suppression of TE activity, the data presented in this manuscript do not exclude that piRNA, similar to miRNAs, have another role besides repression of TE activity that is beneficial for acquiring oocyte developmental competence.\nResponse: We thank the reviewer for this insightful comment. We agree that, although TE silencing remains the best-characterized function of the piRNA pathway, growing evidence suggests that piRNAs and PIWI proteins also participate in gene regulation (15-18). In our analysis, we found that approximately 35-45% of piRNAs originate from repetitive elements, whereas ~5–10% map to annotated coding or exonic non-TE regions (Supplementary Fig. 4b). Notably, these region-aligned piRNAs lacked the characteristic 10A bias, suggesting potential roles in gene regulation through Ping-Pong–independent mechanisms (Supplementary Fig. 4c).\nIn the revised manuscript, we have now highlighted the presence and characteristics of these gene-derived piRNAs in the Result section, as stated:\n“Although the majority of piRNAs mapped to various TE families, a notable minority (5-10%) originated from coding or other exonic non-TE genomic regions (Supplementary Fig. 4b, c). Unlike TE-related piRNAs, these piRNAs lacked the characteristic 10A bias (Supplementary Fig. 4c), suggesting they may participate in gene regulation through Ping-Pong–independent mechanisms.” (Line 283-287).\nWe also expanded the Discussion to address the potential regulatory roles pf these gene-aligned piRNAs during oocyte maturation and early embryogenesis, which now reads:\n“This study has two important limitations. First, although TE silencing is the well-established function of the piRNA pathway (19), growing evidence indicates that piRNAs and PIWI proteins also participate in broader gene regulatory processes (15-18). Here, we did not assess whether piRNAs engage in additional regulatory mechanisms that may influence oocyte competence. Nevertheless, the single-oocyte co-sequencing dataset generated in this study provides a valuable foundation for uncovering potential noncanonical piRNA-mediated regulatory mechanisms during human oogenesis.” (Line 468-474).\nOn the same note, expression of several TE increase during oocyte development (line 148; fig 1f). What is the proposed function of this increased expression. Could it be important for acquisition of developmental competence?\nResponse: We sincerely thank the reviewer for this insightful question. To address it, we performed a comparative analysis of upregulated and downregulated TE subfamilies, examining the expression levels of sense and antisense piRNAs, the strength of Ping-Pong signatures for both short- and long-piRNAs, overall 1U/10A ratios, and TE transcript levels (Supplementary Fig. 10a). Upregulated TE subfamilies displayed significantly higher transcript levels than downregulated TE subfamilies. However, apart from their elevated transcript abundance, all other piRNA-related features were markedly lower in upregulated TEs—most notably, the abundance of corresponding antisense short- and long-piRNAs. Consistently, these up-regulated TE subfamilies exhibited weaker Ping-Pong signatures, reflected by reduced 10-nt overlap ratios and diminished 1U/10A biases.\nSupplementary Figure 10. Comparative characteristics of up- and down-regulated TE subfamilies.\n(a) Comparison of up- and down-regulated TE subfamilies based on related piRNA abundance, Ping-Pong Z-scores, overall 1U and 10A ratios, and TE transcript levels.\n(b) Boxplots showing the density of upregulated and downregulated TE insertions within the top 20 piRNA clusters, separated by strand orientation (Antisense vs. Sense). For both (a) and (b), statistical significance was calculated using a t-test.\nWe next analyzed the genomic distribution of upregulated and downregulated TE subfamilies within the top 20 piRNA clusters, which collectively generate over 90% of uniquely mapped piRNAs. Although upregulated TEs were present within these clusters, their antisense insertions were significantly less frequent than those of the downregulated TE subfamilies, consistent with their reduced antisense short- and long-piRNA expression levels (Supplementary Fig. 10b).\nNotably, seven out of the eight upregulated TE subfamilies are human-specific elements. Given their relatively high transcript abundance and reduced antisense piRNA expression, we propose that these upregulated elements likely represent evolutionarily younger TE subfamilies for which the piRNA/PIWI pathway has not yet fully adapted to achieve efficient repression—potentially due to insufficient antisense insertions within highly productive piRNA clusters.\nNevertheless, we cannot exclude the possibility that increased expression of these TE subfamilies may contribute functionally to oocyte competence, rather than merely reflecting incomplete silencing. Indeed, several studies have shown that specific TE families can contribute to early embryogenesis (20, 21). For instance, full-length MERVL transcripts are required for proper regulation of the host transcriptome and chromatin organization during early embryogenesis, and repression of MERVL results in embryonic lethality due to defects in differentiation and genomic stability (20).\nTogether, these new analyses are presented in the revised Supplementary Fig. 10 and described in the updated Results section:\n“Compared with TE subfamilies that were downregulated, these upregulated TE subfamilies showed higher transcript abundance but produced markedly fewer corresponding antisense and sense short- and long-piRNAs, accompanied by weaker Ping-Pong signatures and reduced 1U and 10A biases (Supplementary Fig. 10a). Genomic distribution analyses within the top 20 piRNA clusters further revealed substantially fewer antisense insertions for these upregulated TE subfamilies relative to downregulated TE subfamilies, consistent with their reduced antisense piRNA levels (Supplementary Fig. 10b).” (Lines 333-341).\nWe have also added a brief discussion highlighting the potential biological relevance of increased TE subfamily expression during oocyte development:\n“Notably, several TEs showed increased expression during oocyte development (Fig. 1f and Fig. 5e). Given their elevated transcript levels and the concomitant reduction in antisense piRNAs (Supplementary Fig. 10), these elements may represent evolutionarily young TEs that the PIWI-piRNA system has not yet fully adapted to silence, potentially due to insufficient antisense insertions within highly productive piRNA clusters. An alternative possibility is that their upregulation may have functional relevance for oocyte maturation. Indeed, specific TEs have been shown to play essential roles in early embryogenesis (22, 23). For example, full-length MERVL transcripts are required for proper transcriptome regulation and chromatin organization, and loss of MERVL results in embryonic lethality due to impaired differentiation and genomic instability (23). Whether the increased expression of specific TEs contributes to the acquisition of oocyte developmental competence, rather than merely reflecting incomplete silencing, therefore warrants further investigation.” (Lines 455–467).\nFig 2g is difficult to read, individual lines/arrows cannot be distinguished. As such the figure is not informative.\nResponse: We have revised Fig. 2g (now Fig. 2i) into bubble charts to more clearly illustrate whether different miRNA families potentially target the same pathways. In the updated figure, bubble size represents the gene ratio, defined as the proportion of predicted target genes of a given miRNA family that fall within an enriched pathway relative to the total number of genes annotated to that pathway. Bubble color indicates the p-value, with more intense colors reflecting stronger pathway enrichment.\nFigure 2. piRNAs are the predominant small non-coding RNAs during human oogenesis.\n(i) Bubble plots display the top enriched Reactome pathways for representative miRNA families that are highly expressed in earlier (PrF and PF, panel a) and late (GV and MII, panel b) stages of human oocytes. Bubble size reflects the gene ratio (the proportion of predicted target genes within each pathway), and bubble color represents the p-value, reflecting the significance of pathway enrichment.\nFigure 3c, ‘ Monkey’ is a general term, it should be further specified; eg Macaca fascicularis, or whichever species was used here.\nResponse: The monkey used for oocyte collection in this study was cynomolgus monkey (Macaca fascicularis), and we have now specified this information wherever relevant in the text and figures of the revised manuscript.\nLine 267 The Ping-Pong signature was either undetectable or 267 substantially diminished for most TE-aligned long-piRNAs. How does this vibe with other publications such as PMID: 25818294?\nResponse: We thank the reviewer for this insightful suggestion. In the study by Elke F. Roovers et al. (12), the authors examined PIWI protein expression and piRNA populations in bovine and human ovaries. To evaluate whether the strength of the Ping-Pong signature observed in their dataset is consistent with ours, we reanalyzed their publicly available small RNA-seq data from adult human ovaries. Whereas Roovers et al. analyzed whole ovarian tissues containing a mixture of cell types, our study focused specifically on single oocytes at defined developmental stages. Nevertheless, because adult ovaries primarily contain oocytes at the primordial follicle stage, their dataset provides a meaningful parallel comparison with ours.\nOur reanalysis revealed that long-piRNAs demonstrated Ping-Pong signatures for SVA and Alu elements, whereas L1 and ERV elements exhibited very low or undetectable signals in the human ovarian samples from Roovers et al (Supplementary Figure 7). Moreover, even for SVA and Alu, the Ping-Pong signatures of long-piRNAs were weaker than those of short-piRNAs. These observations are consistent with our own results, particularly for early follicular-stage oocytes (PrF- and PF-stage), reinforcing the reliability of our data and conclusions.\nWe have incorporated this analysis into Supplementary Figure 7 and added a brief description of the results in the revised manuscript, which now reads:\n“Analysis of TE-related piRNAs using the adult human ovary dataset from a previous study (12) further confirmed that long-piRNAs exhibited weaker Ping-Pong signatures than short-piRNAs for SVA and Alu elements, whereas L1 and ERV showed very low or undetectable Ping-Pong signatures (Supplementary Figure 7). (Lines 304-308).\nSupplementary Figure 7. Ping-Pong signatures of adult human ovarian piRNAs from Roovers et al. (12) Ping-Pong signatures for LINE/L1, SINE/Alu, LTR/ERV, and SVA-related short (17-24 nt) and long-piRNAs (25-35) in adult human ovaries. Z-scores for 10-nt overlaps are shown in each panel.\nThe conclusions are too strong. For instance line 352, a role for piRAs has not been provided, but suggested. By examining piRNA and TE expression levels, it is by no means demonstrated that these affect each other. Line 352, the results do not ‘emphasize the suppressive effects of piRNAs on TE’, but instead provide a correlation between piRNAs and TEs.\nResponse: We agree that our data demonstrate correlative rather than causal relationships between piRNAs and their corresponding TEs. In the revised manuscript, we have carefully refined the wording in the Results section to more accurately reflect this interpretation. Specifically, in line 352 and related passages, we replaced statements such as “emphasize the suppressive effects of piRNAs on TEs” with “providing a comprehensive view of the relationship between short- and long-piRNAs and TE downregulation.” (Line 418-419).\nLine 364-371: It should be made clear that (control of) TE activity in germ cells is different from TE activity in somatic cells.\nResponse: We thank the reviewer for this insightful comment. We agree that TE regulation differs fundamentally between somatic and germline tissues. In somatic contexts, TE repression operates largely independently of small RNA pathways. By contrast, in germ cells, TE activity is tightly controlled by the PIWI–piRNA pathway to preserve genomic integrity and ensure proper gametogenesis.\nTo clarify this distinction, we have revised the Discussion section to explicitly state the difference between TE control in somatic versus germ cells, which now reads:\n“In somatic cells, TE repression is maintained largely independent of small RNA pathways (24-26). By contrast, in germ cells, TE activity is tightly regulated by the PIWI-piRNA pathway, which safeguards genomic stability throughout gametogenesis (19, 27-29).” (Lines 414-417).\nIn addition, to avoid potential misunderstanding, we have removed the previous description suggesting ERV activity in somatic cells in the revised manuscript.\nMinor\nLine 34 ‘(long-piRNAs)’ can be removed.\nResponse: We thank the reviewer for this suggestion. We have removed “(long-piRNAs)” as suggested.\nReviewer #3 (Remarks to the Author):\nManuscript No.: NCOMMS-25-40496 “Integrated small and long RNA sequencing in single oocytes reveals piRNA-mediated transposon repression during human oogenesis” Fengjuan Zhang, Hongdao Zhang, Yali Xiao, Miao Liu, Aimin Ren, Suying Liu2, Ligang Wu Although the piRNA pathway has been well characterized in model organisms, little is known about its functions in human gametogenesis. Large scale exome sequencing projects have identified mutations in piRNA pathway proteins that cause human male infertility, but much less is understood about piRNAs in human oocytes. In this manuscript, the authors profile the piRNAs and long RNAs in four stages of developing human oocytes. The analyses are of high quality, the data clearly presented, and the authors interpretations well supported by data and clearly described. The manuscript will be of broad interest to those studying piRNAs, gametogenesis, mammalian fertility, and TE silencing. My few concerns should be easily addressed by adding a few additional analyses and revising the text.\nResponse: We sincerely thank the reviewer for the positive evaluation of our study and for the constructive comments and suggestions. We have carefully addressed all points raised, adding additional analyses and revising the text accordingly to further improve the clarity and completeness of the manuscript.\nMajor Concerns 1.\nI don’t understand the logic of these conclusions, “Collectively, these results suggest that long-piRNAs exhibit a conserved, preferential targeting of ERVs across species, indicating their specialized role in ERV repression. In contrast, short-piRNAs exhibit broader TE targeting, likely reflecting their more generalized silencing function across multiple TE families” (lines 211-214). I agree that antisense piRNAs are expected to target sense TE transcripts, but why would sense piRNAs have any function at all, other than to help amplify antisense piRNAs via the Ping-Pong cycle? Are antisense TE transcripts detected in the long RNA-seq data? Other than piRNA precursors, are there transcripts with extensive complementarity to sense piRNAs?\nResponse: We thank the reviewer for raising this question. To clarify this point, we have revised the manuscript to explicitly include strand-specific information and to distinguish the regulatory patterns of TE-related sense and antisense piRNAs. The revised text now states:\n“Collectively, these results suggest that TE-related antisense long-piRNAs exhibit a conserved, preferential targeting of ERVs across species, indicating their specialized role in ERV repression. In contrast, TE-related antisense short-piRNAs exhibit broader TE targeting, likely reflecting their more generalized function across multiple TE families.” (Line 226-230)\nIn this study, we primarily focused on the roles of short- and long-piRNAs in TE regulation and observed distinct targeting preferences between these two piRNA populations across different TE superfamilies. However, we also recognize that the piRNA may contribute to regulatory processes beyond canonical TE silencing, particularly when TE-derived sense or antisense sequences are embedded with host protein-coding or noncoding genes. Indeed, Watanabe et al. demonstrated that in mouse spermatocytes, many mRNAs and lncRNAs harboring retrotransposon-derived sequences are repressed by the PIWIL1-piRNA pathway, and a pseudogene (Stambp-ps1) generates piRNAs that repress its parental gene (Stambp) via direct base-pairing and PIWIL1 slicer activity (30). These findings suggest that piRNAs may have broader regulatory roles beyond their canonical function in TE repression.\nBecause Smart-seq2 used in this study is not strand-specific, our long-RNA data cannot determine transcript orientation. Therefore, we assigned TE or gene transcript orientation based on reference annotation (TE annotation downloaded from UCSC genome browser and RepBase; gene annotation downloaded from GECODE V19). TE-related sense and antisense piRNAs were analyzed separately, and their distribution across exonic and intronic regions of different gene categories was examined (Rebuttal Figure 8). Both TE-related sense and antisense piRNAs were dispersed across multiple gene types, suggesting potential targeting activities beyond canonical TE-derived transcripts.\nWithin exonic regions, we observed that TE-related sense piRNAs more frequently aligned to the sense strand of protein-coding genes (>5%) compared with antisense piRNAs (~2.5%) (Rebuttal Figure 8a). Conversely, TE-related antisense piRNAs showed greater enrichment on the antisense strand of lncRNA, antisense transcripts, and pseudogenes. These observations indicate potential strand-specific regulatory interactions between piRNAs and exonic regions across diverse gene categories, mediated by the presence of TE-derived sequences.\nWithin intronic regions, a greater number of TE-related piRNAs were distributed across multiple gene categories and displayed distinct patterns compared with exonic regions. In addition to protein-coding genes, lincRNAs, antisense genes, and pseudogenes, TE-related piRNAs frequently aligned to intronic regions of processed transcripts and sense-overlapping loci (Rebuttal Figure 8b). Notably, TE-related sense piRNAs preferentially aligned to the antisense strand of sense-overlapping genes, whereas TE-related antisense piRNAs were more frequently associated with the sense strand of processed transcripts. These patterns raise the possibility that TE-related piRNAs may modulate RNA processing, transcript maturation, or co-transcriptional regulation through recognition of intronic TE-derived sequences.\nAlthough these observations require further validation and the underlying mechanisms remain to be elucidated, our single-oocyte co-sequencing dataset provides a valuable resource for future studies investigating potential piRNA-mediated regulation beyond classical TE repression.\nRebuttal Figure 8. Distribution of TE-related sense and antisense piRNAs across gene regions.\n(a) Strand-specific alignment of TE-related sense and antisense piRNAs to exonic regions across gene categories.\n(b) Strand-specific alignment of TE-related sense and antisense piRNAs to intronic regions across gene categories. TE-related piRNAs from all oocyte stages were aggregated for this analysis. Proportions were calculated relative to total TE-related sense or antisense piRNAs.\n2. Was phased piRNA biogenesis detected among short or long piRNAs? Please add plots and Z-scores for the 3′ to 5′ and 5′ to 5′ distances of piRNA ends on the same genomic strand.\nResponse: To determine whether phased piRNAs occur in developing human oocytes, we analyzed the 3’-to-5’ and 5’-to-5’ distances between piRNA ends mapped on the same genomic strand and calculated Z-scores for phasing enrichment, following the analytical framework described by Bo W Han et al. (31). To validate our analytical pipeline, we first reanalyzed the published Tdrkh⁺/⁻ and Tdrkh⁻/⁻ mouse testis datasets, where strong phasing signatures have been reported (31, 32). Consistent with previous studies, we observed a prominent 1-nt phasing signal in the 3’-to-5’ distance distribution and broad peaks around 30-40 nt in the 5’-to-5’ distance distribution, hallmark features of phased piRNA processing (Rebuttal Figure 9a, b). Nucleotide composition analysis of downstream phased piRNA further showed a pronounced 1U bias at their 5’ ends (Rebuttal Figure 9c, d).\nRebuttal Figure 9. Phasing signatures of piRNAs in TDRKH-/- and TDRKH+/- mouse testes.\n(a) Probability distribution of 3′-to-5′ distances between piRNAs mapped on the same genomic strand in Tdrkh-/- and Tdrkh+/- mouse testes. Z-scores for position 1 (phasing) are indicated.\n(b) Probability distribution of 5′-to-5′ distances between piRNAs mapped on the same genomic strand in Tdrkh-/- and Tdrkh+/- mouse testes.\n(c, d) Nucleotide composition of the last nucleotide of upstream piRNAs and the first nucleotide of downstream piRNAs for pairs separated by a 1-nt 3′-to-5′ distance in Tdrkh-/- (c) and Tdrkh+/- (d) mouse testes.\nIn developing human oocytes, the 3’-to-5’ distance analysis revealed significant enrichment at a 1-nt separation between adjacent piRNAs, with Z-scores ranging from 3.6 to 6.8 (Rebuttal Figure 10a). Short-piRNAs exhibited slightly stronger phasing signals than long-piRNAs. Nucleotide composition analysis of downstream phased piRNA showed a pronounced 1U bias at their 5’ ends, particularly among short-piRNAs (Rebuttal Figure 10c-d), consistent with canonical phased piRNA processing signatures observed in other mammals. No obvious peaks were detected in the 5′-to-5′ distance analysis (Rebuttal Figure 10b).\nInterestingly, we also observed even stronger enrichment at distance = 0 than at distance = 1 (phasing) in the 3’-to-5’ distance analysis (Rebuttal Figure 10a), indicating frequent 1-nt overlaps between adjacent piRNAs. These 1-nt-overlapping downstream piRNAs similarly exhibited a strong 5′-terminal U bias, while their upstream counterparts exhibited marked enrichment for 3’-terminal U (Rebuttal Figure 10e-f). We speculate that this distinctive overlap pattern may arise from the high frequency of 3’ non-templated U tailing. However, this interpretation cannot be definitively tested with the current data, as terminal U residues cannot be reliably distinguished from tailing events. Future studies—such as perturbation of enzymes involved in piRNA 3’-end tailing—will be necessary to evaluate this model.\nRebuttal Figure 10. Phasing signatures for short- and long-piRNAs in human developing oocytes.\n(a) Probability distribution of 3’-to-5’ distances between piRNAs mapped on the same genomic strand across human oocyte stages (PrF, PF, GV, and MII stages). Z-scores for position 1 (phasing) and position 0 are indicated.\n(b) Probability distribution of 5’-to-5’ distances between piRNAs mapped on the same genomic strand across human oocyte stages (PrF, PF, GV, and MII stages).\n(c) Nucleotide composition of the last nucleotide of upstream piRNAs and the first nucleotide of downstream piRNAs for short-piRNA pairs with a 1-nt 3’-to-5’ distance.\n(d) Nucleotide composition of the last nucleotide of upstream piRNAs and the first nucleotide of downstream piRNAs for long-piRNA pairs with a 1-nt 3’-to-5’ distance.\n(e) Nucleotide composition of the last nucleotide of upstream piRNAs and the first nucleotide of downstream piRNAs for short-piRNA pairs with a 0-nt 3’-to-5’ distance.\n(f) Nucleotide composition of the last nucleotide of upstream piRNAs and the first nucleotide of downstream piRNA for long-piRNA pairs with a 0-nt 3’-to-5’ distance.\nShort- and long-piRNAs were analyzed separately. Data are from 4-9 representative biological samples, excluding piRNAs with non-templated 3’ nucleotides.\n3. What fraction of short piRNAs were prefixes of long piRNAs (i.e., same 5′ ends)? Are the short piRNAs and long piRNAs produced from the same precursors, then trimmed to the length characteristic of the bound PIWI protein (i.e., they are prefixes) or produced by separate mechanisms?\nResponse: To determine whether short-piRNAs are produced from the same precursors as long-piRNAs, we analyzed the positional relationships between their 5′ and 3′ ends. Specifically, we compared the genomic coordinates of short- and long-piRNAs and calculated the probabilities that their 5′ or 3′ ends coincide. In developing human oocytes, long-piRNAs were substantially more likely to share 5′ ends with short-piRNAs than 3′ ends (0.51 vs. 0.02 in PrF; 0.68 vs. 0.01 in PF; 0.66 vs. 0.01 in GV; and 0.63 vs. 0.01 in MII). These results indicate that approximately half of long-piRNAs arise from the same 5′-cleaved precursors that generate short-piRNAs. This pattern agrees with previous findings in golden hamster MII oocytes (4) and supports a model in which short- and long-piRNAs are generated from common precursors but undergo differential 3’-end processing, likely determined by their associated PIWI proteins. Thus, the observed size differences between the two piRNA populations most likely reflect PIWI-dependent 3′-end maturation rather than distinct biogenesis pathways.\nThese new analyses are presented in revised Supplementary Fig. 2 and described in the corresponding Results section (Line 181-185) of the revised manuscript, which reads as:\n“Genomic coordinate analysis showed that short- and long-piRNAs were far more likely to share 5’ ends than 3’ ends across all stages (PrF: 0.51 vs. 0.02; PF: 0.68 vs. 0.01; GV: 0.66 vs. 0.01; MII: 0.63 vs. 0.01) (Supplementary Fig. 2), suggesting that both populations are derived from the same 5’-cleaved precursors, consistent with findings in golden hamster (4).”\nSupplementary Figure 2. 5’- and 3’-end overlap between short and long-piRNAs\nProbability distributions of distance between the 5’ ends and 3’ ends of short and long-piRNAs. Numbers indicate the total frequency of coincident ends. Data are based on 4-9 representative biological samples.\nMinor Points\n1. I very much enjoyed reading the Introduction. It was well written and should make the paper more accessible to non-specialists.\nResponse: We sincerely thank the reviewer for the positive feedback.\n2. Figs 1f and 5d: It is hard to estimate fold-change based on circle size. Perhaps vertical bars could be used instead, with a scale on the right? It is far easier to compare the length of two objects than their areas.\nResponse: We thank the reviewer for the suggestion. We have changed the circle-based plot into bar-based plot in Fig. 1f and 5d instead.\n3. Significant figures: it seems unlikely that the average increase or decrease in expression of an RNA can be measured to an accuracy of four or five significant 1 figures. For example, it would be better to write 35-fold instead of 34.93-fold, 3.3-fold instead of 3.39-fold, etc., using just two significant figures throughout.\nResponse: We thank the reviewer for this suggestion. All fold-change values have been rounded to two significant figures throughout the manuscript and figures.\n4. Violin plots of small data sets are not really useful, since there is too little data to reliably estimate the probability distribution of values. The choice of kernel is more likely to influence the shape of the plot than the actual data. Typically a minimum of ten data points is recommended for violin plots. Figures such as Fig. 2c and 4b should just show the individual data points and the median (dot plots).\nResponse: In the revised manuscript, we have replaced all bar graphs based on small data sets with dot plots and added statistical significance indicators to the updated figures.\n5. Lines 166-167: “3’-O-methylation” should be “3′ terminal 2′-O-methylation”.\nResponse: We thank the reviewer for pointing out this error, and have made the correction.\n6. Please show the oxidized piRNA sequencing separately from the non-oxidized data in Figure 2.\nResponse: We thank the reviewer for this valuable suggestion. In the original manuscript, small RNA oxidation sequencing was not performed because our previous studies had already demonstrated that short-piRNAs in GV and MII oocytes lack 2′-O-methylation, whereas long-piRNAs possess this modification (1, 2). Furthermore, oocytes at early follicular stages yield extremely limited amounts of RNA, and the scarcity of human oocytes makes the construction of oxidation sequencing libraries technically challenging. Nevertheless, as suggested, we collected oocytes at the primordial and primary follicular stages and performed oxidation sequencing. Comparison between oxidized and non-oxidized libraries revealed that short-piRNAs were largely eliminated after oxidation, whereas long-piRNAs were retained, particularly in PF-stage oocytes (Fig. 2e). These findings confirm that in oocytes at the primordial and primary follicular stages, PIWIL3-associated short-piRNAs lack 2’-O-methylation at their 3’ ends, while long-piRNAs carry this modification, consistent with our observations in GV and MII oocytes. The incomplete depletion of short-piRNAs is likely due to the limited input RNAs (approximately 2.5 oocytes per reaction), which may reduce the efficiency of the oxidation treatment.\nThese new data have been incorporated into revised Fig. 2e, and the corresponding descriptions have been added to the Results (Lines 173-179) and Materials and Methods (Lines 541–552) sections of the revised manuscript.\nFigure 2. piRNAs are the predominant small non-coding RNAs during human oogenesis.\n(e) Length profiles of small RNAs before and after oxidation treatment in primordial follicle oocytes. Short-piRNAs (17-24 nt) were largely depleted upon oxidation, whereas long-piRNAs (25-35 nt) were largely retained, confirming the absence of 2’-O-methylation in short-piRNAs and its presence in long-piRNAs. Approximately 2.5 oocytes were pooled per sample for each oxidation assay.\n7. Is the heatmap in Fig. 2f on a linear or log2 scale?\nResponse: For each miRNA, we used RPM values as input and performed row-wise normalization by calculating Z-scores. The heatmap, therefore, displays normalized expression values represented as row-wise Z-scores. We chose Z-score normalization because it emphasizes relative expression patterns across samples, facilitating the comparison of expression dynamics for each miRNA during human oocyte development. We have now added this information to the figure legend.\n8. Throughout the text, please clarify whether the piRNA or long RNA data being discussed is restricted to uniquely mapping reads or includes multiply mapping reads, and, if it includes multi-mappers, whether they were apportioned among all loci or randomly assigned to individual loci.\nResponse: For multiply mapped piRNAs or long reads, each read was assigned to all possible mapped loci and categorized by TE superfamily in the following order: LINE/L1 → SINE/Alu → LTR/ERV → LTR/others → SVA → others. Because individual TE subfamilies (e.g., L1PA1 and L1PA2) share high sequence similarity, we applied a greedy mapping strategy to capture all potential alignments. If a read or piRNA species mapped to multiple TE subfamilies within the same superfamily, it was counted in each subfamily. However, when summarizing results at the superfamily level (e.g., LINE/L1), each piRNA or read was counted only once. For instance, a piRNA mapping to both L1PA1 and L1PA2 was counted once in the LINE-1 superfamily summary but included in both subfamilies when L1PA1- and L1PA2-derived piRNAs were analyzed separately. This strategy was applied consistently across all TE-related analyses, including TE transcript level, TE-associated short- and long-piRNAs, and strand-specific subsets, to ensure uniform quantification and accurate cross-dataset comparison.\nWe have carefully revised the text to clarify which long RNA reads or piRNA categories (uniquely mapped, multiply mapped, or both) were used in each analysis and have explicitly specified these details in the Materials and Methods section (Lines 565-573) and relevant figure legends.\n9. In addition to the heat map of Z-scores, it would be helpful to see the overall plot of Ping-Pong for all piRNAs (i.e., 5′ to 5′ distance on opposite genomic strands). (with different subfamily of TEs)\nResponse: We analyzed the Ping-Pong signatures of all piRNAs associated with different TE superfamilies. Consistent with our observations for piRNAs related to individual TE subfamilies shown in Fig. 4f, all piRNAs related to Alu and SVA elements exhibited the strongest Ping-Pong activity, followed by moderate activity for L1-related piRNAs, and low activity for ERV-related piRNAs. During human oocyte development, Ping-Pong signals for Alu- and SVA-related piRNAs increased markedly, whereas those for L1-related piRNAs showed only a minor increase, and ERV-related piRNAs remained largely stable or slightly decreased relative to the PrF stage.\nThese new results are presented in the revised Supplementary Fig. 4d and described in the corresponding Results section (Line 262-264) of the revised manuscript, which reads as\n“Among TE superfamilies, Alu- and SVA-associated piRNAs exhibited the strongest Ping-Pong signatures, followed by moderate activity for L1 and weak activity for ERV (Supplementary Fig. 4d).”\nSupplementary Figure 4. Differential upregulation patterns of TE superfamily-aligned short-piRNAs and presence of Ping-Pong signatures during human oogenesis\n(d) Ping-Pong signatures of all piRNAs mapped to LINE, SINE/Alu, LTR/ERV, and SVA superfamilies. Z-scores for 10-nt overlaps are shown in each stage of each panel.\n10. Please add Z-scores to the Ping-Pong plots in Fig. 5f and 5g.\nResponse: We have updated Fig. 5f and 5g by adding Z-scores, as well as all other relevant figures that require showing the strength of the Ping-Pong signature.\n11. It would be helpful to show long RNA coverage by genomic strand in Fig. 6e and 6f.\nResponse: We have updated Fig. 6e, Fig. 6f, and Supplementary Fig. 12 and 14 to include long RNA coverage by genomic strand.\n12. Is APOBEC expressed in the developing human oocyte? If it is, do TE transcripts show base changes consistent with deamination of C to U?\nResponse: We examined the expression of all APOBEC family members, including APOBEC1, APOBEC2, APOBEC3A-H, and APOBEC4, in our dataset. Nearly all APOBEC genes were undetectable during human oogenesis, except for APOBEC3C, which exhibited very low expression in MII-stage oocytes.\nRebuttal Figure 11. Expression of APOBEC family members in developing human oocytes.\nLog transformed expression level (FPKM) was shown across developmental stages."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made considerable efforts to improve the manuscript and incorporate the suggestions raised in the previous submission. The story is cleaner and less fragmented, and they have carefully addressed each concern that was raised. There are no further technical concerns about the content and presentation of the manuscript in its current form.", "reviewer_2": "All the questions and remarks I provided for the original manuscript have been properly addressed and the manuscript has been adapted following those questions. In that the authors did an excellent job, and I now consider this manuscript suitable for publication in Nature communications.\nThere are only 2 typos that I identified in the figures:\nFigure 2 panel E, title of the graph: beforre should be before\nSupplementary Figure 10, panel b: Down-reguated should be Down-regulated, Up-reguated should be Up-regulated\nIn addition in the figures both Ping-Pong, ping-pong and Ping-pong are used, but I do not expect this would lead to any confusion.", "reviewer_3": "The authors have satisfied all my concerns, with one remaining issue: many of the reviewers’ comments, including my own, are addressed in the authors’ rebuttal and “rebuttal figures,” but the new data and analyses are not always incorporated into the manuscript itself. For example, the phasing analysis (rebuttal figures 9 and 10) needs to be included and discussed in the manuscript itself. Otherwise, the study is ready to be published in Nature Communications.\nReviewer: Phil Zamore"}, "author_response": "We appreciate the reviewer’s thorough reassessment and positive comments.\nResponse: We thank the reviewer for the positive evaluation and support for publication.\nResponse: We thank the reviewer for identifying the typographical errors. These have been corrected in Figure 2E and Supplementary Figure 10b (now Supplementary Figure 11b).\nResponse: We thank the reviewer for noting this point. For consistency, the terminology has been unified to “Ping-Pong” throughout all figures.\nResponse: We thank the reviewer for this important comment. The phasing analysis previously presented in rebuttal Figures 9 and 10 has now been incorporated into the manuscript as Supplementary Figure 3, with corresponding description and discussion added to the Results section. Relevant text has been updated accordingly (Line 178-201)."}]

Decision letter