The non-canonical bivalent gene Wfdc15a controls spermatogenic protease and immune homeostasis Free

Infertility affects 8-12% of couples, among which 50% is caused by male factors as primary or contributing causes (Agarwal et al., 2021). Although the highest frequency of genetic factors is found in azoospermia, ∼20% of azoospermia cases remain idiopathic, largely owing to the complexity of spermatogenesis, which involves at least 2000 genes (Krausz and Riera-Escamilla, 2018). Spermatogenesis is a multi-stage process to produce spermatozoa (haploid) involving sequential cell differentiation from spermatogonia (diploid cells), including spermatogonial stem cells (SSCs), to spermatocytes (meiosis division I) and finally to haploid spermatids (meiosis division II) (Handel and Schimenti, 2010).

In addition to genetic causes, evidence suggests the importance of epigenetic regulation in spermatogenesis and male infertility. SSC differentiation requires DNA methylation factors and histone modifiers (Dura et al., 2022; Kuroki et al., 2013; Lambrot et al., 2015; Lin et al., 2022; Shirakawa et al., 2013; Tomizawa et al., 2018). Meiotic progression from spermatocytes to spermatids is compromised by the depletion of major epigenetic factors, such as Prmt1, Np95 (also known as Uhrf1), Dnmt1, and Dnmt3a (Barau et al., 2016; Bourc'his and Bestor, 2004; Hata et al., 2006; Hirota et al., 2018; Kaneda et al., 2004; Tachibana et al., 2007; Takada et al., 2021; Waseem et al., 2021). Histone modifications also exert a regulatory influence on differentiation across an extended temporal range: SSCs and other spermatogenic cells possess bivalent and monovalent histone modifications that pre-program genes required at later processes such as spermatid maturation and embryonic development (Hasegawa et al., 2015; Lesch et al., 2013; Tomizawa et al., 2018). Canonical bivalent chromatin is poised for gene expression and typically carries KMT2B-mediated H3K4me3 and Polycomb-mediated H3K27me3 at a promoter, which is later resolved for activation or silencing upon differentiation (Bernstein et al., 2006; Denissov et al., 2014; Hanna et al., 2018; Tomizawa et al., 2018). The crucial roles of one of the KMT2B-programmed canonical bivalent genes in spermatids have been validated by a knockout (KO) study (Kobayashi et al., 2021). Monovalent genes found in the male germline are also H3K4me3-modified by KMT2B but lack H3K27me3 in SSCs (Tomizawa et al., 2018). They are generally repressed in SSCs by unknown mechanisms, but many are later activated in spermatids, indicating potential functions in spermiogenesis. However, despite spermatids expressing a high proportion of known infertility genes (Murat et al., 2023), the nature and importance of over 200 monovalent genes expressed in spermatids remain largely elusive (Tomizawa et al., 2018). Several other forms of bivalent chromatin have been reported. Among them, evidence suggests that H3K4me3/H3K9me3 non-canonical (nc)-bivalent genes play important roles for cellular differentiation in various types of tissues (Du et al., 2021; Khromov et al., 2011; Matsumura et al., 2015; Rugg-Gunn et al., 2010). Thus, a subset of monovalent H3K4me3 genes might carry H3K9me3 as a repressive mark. However, the role of H3K4me3/H3K9me3 nc-bivalent modification in the male germline on epigenetic characteristics and physiological significance have not been well explored compared with the classical canonical bivalent genes.

Postnatal spermatogenic cells including spermatids are physically and nutritionally supported by Sertoli cells. Distinct combinations of different germ cells are associated with each Sertoli cell, which can be classified into stages I to XII in the mouse (Meistrich and Hess, 2013). The release of spermatozoa into the tubule lumen is tightly regulated by hormones at the hypothalamic-pituitary-testicular axis, and takes place at late stage VIII (Cheng and Mruk, 2015). Until then, spermatids adhere to Sertoli cells with apical ectoplasmic specialization (ES), which consists of the adherens junction and the actin-based anchoring junction (Cheng and Mruk, 2015). In Sertoli cells at late stage VIII, neural Wiskott-Aldrich syndrome protein (N-WASP; Wasl) and cortactin activate actin-related protein 2/3 (Arp2/3; Actr2/3) and convert microfilaments from the bundled to de-bundled configuration. These dynamic changes are followed by apical-ES degeneration, disruption of the tubulobulbar complex (TBC) and induction of matrix metallopeptidase 2 (MMP2), which finally leads to spermiation (O'Donnell et al., 2011). Failure of apical ES formation by genetic disturbance or anti-androgenic reagents, as well as abnormality in spermatids themselves, could lead to premature release of spermatids (exfoliation) and may cause infertility (Li et al., 2017; Wong et al., 2008; Xu et al., 2021).

Meiotic germ cells, including spermatids that emerge after the establishment of immune competence, are under a unique immunological regulation: they are protected from the autoimmune attack by the basal ES structure, the tight junction and gap junction-based adhesion between two opposing Sertoli cells (Cheng et al., 2011b). This structure makes testes an immuno-tolerant organ and functions as barriers from auto-antigen and exogenous invasions. The basal ES is dynamic, as remodeling by MMPs increases the permeability of the tight junction (Yao et al., 2010). Disruption of basal ES is associated with bacterial/viral infection, spermatic antigens and low-grade inflammatory conditions such as obesity (Fan et al., 2015; Govero et al., 2016). These extrinsic and intrinsic molecules are recognized via pattern recognition receptors (PRRs) such as Toll-like receptors (TLRs) and nucleotide-binding oligomerization domain-like receptors (NLRs) (Fitzgerald and Kagan, 2020). Activation of PRRs often induces adhesion molecules accelerating inflammation, migration of effector cells of the innate immunity and induces the release of pro-inflammatory cytokines such as TNFα, IL1, IL6 and IL10. These adverse innate immunity reactions in the testes can interfere with spermatogenesis.

To gain insights into the epigenetic regulation of spermatogenic genes by the H3K4me3/H3K9me3 nc-bivalent mechanism, we investigated a subset of monovalent genes that carry H3K9me3 in spermatogonia. We identify Wfdc15a, which belongs to the whey acidic protein four-disulfide core (WFDC) serine-type protease inhibitor family that is linked to innate immunity and inflammation (Bingle and Vyakarnam, 2008; Clauss et al., 2005; Nakajima et al., 2019; Scott et al., 2011; Small et al., 2017), as a novel H3K4me3/H3K9me3 nc-bivalent gene with a potential function in spermiogenesis. Wfdc15a loses H3K9me3 during meiosis and becomes highly expressed in spermatids. Depletion of Wfdc15a in mice not only disturbs protease homeostasis but also causes orchitis-like inflammation associated with TNFα expression in round spermatids. As a result, the KO males show epididymal azoospermia due to a complete loss of spermatids after seminiferous epithelium stage IX. Notably, we find that the testis is a unique organ in that it expresses many WFDC protease inhibitors, suggesting their family-wide involvement in reproduction. These results add a new layer of epigenetic regulation for correct spermatogenesis through modulating protease and immune homeostasis.

To identify previously unreported nc-bivalent genes in spermatogenesis, we re-analyzed the 292 KMT2B-modified H3K4me3⁺/H3K27me3⁻ genes in spermatogonia (Tomizawa et al., 2018). Our re-analysis using germline stem cell (GSC) ChIP-seq data (Tomizawa et al., 2018) identified a small number of nc-bivalent genes (nine genes) carrying both H3K4me3 and H3K9me3 at the promoter in GSCs (Fig. 1A,H; Fig. S1A). The promoters of these genes lost H3K4me3 upon Kmt2b depletion in GSCs (Fig. 1B). We hereafter call these ‘spermatogenic nc-bivalent genes’ (snc-bivalent genes). H3K9me3 is generally a hallmark of constitutive heterochromatin along with DNA methylation and is frequently found at repetitive elements including long terminal repeat (LTR) elements, telomeres and centromeres (Kato et al., 2018; Padeken et al., 2022). ChIP-seq analysis showed that the snc-bivalent genes have moderate levels of H3K9me3 enrichment at the transcription start sites (TSSs), which contrasts to the rest of the monovalent genes having no H3K9me3 enrichment (Fig. 1C). By contrast, entire LTR regions had high levels of H3K9me3. No difference was observed in the H3K4me3 enrichment between snc-bivalent and monovalent genes (Fig. 1D). In contrast to LTRs, DNA methylation on CpG sites of the snc-bivalent gene promoters in undifferentiated (KIT⁻) and differentiating (KIT⁺) spermatogonia (Kubo et al., 2015) was low (mean: 83.3% versus 8.1%; Fig. 1E). Comparison of CpG density showed a lower density in the snc-bivalent TSSs than the canonical bivalent TSSs carrying H3K27me3 (Fig. 1F). Thus, the promoters of the snc-bivalent genes can be regarded as facultative H3K9me3 heterochromatin in spermatogonia. Consistent with this, bulk RNA-seq data (Gan et al., 2013) showed that five of the snc-bivalent genes (Cypt1, Ccdc110, Gm6632, Wfdc15a and 4930469K13Rik) were among the H3K4me3-marked monovalent genes (no H3K27me3) that are activated during spermiogenesis (called Mono-2 genes) (Tomizawa et al., 2018) (Fig. 1G). This included an uncharacterized gene Wfdc15a, which belongs to the WFDC protease-inhibitor family linked to innate immunity and inflammation. As inflammation is a key feature often associated with male infertility (Hasan et al., 2022; Hedger, 2011), we next focused on the regulation of Wfdc15a. Analysis of ChIP-seq and CpG methylation data in GSCs and spermatogonia confirmed the existence of H3K4me3 and H3K9me3 and the absence of H3K27me3 and CpG methylation at the Wfdc15a promoter (Fig. 1H). Furthermore, in vivo ChIP-seq data (Adams et al., 2018; Liu et al., 2019; Maezawa et al., 2020) showed that Wfdc15a promoter loses H3K9me3 during germline differentiation, and inversely acquires H3K27ac as differentiation proceeds (Fig. S1B). These findings strongly suggest that an snc-bivalent gene Wfdc15a is epigenetically regulated in late spermatogenesis and prompted us to investigate the role of Wfdc15a in vivo further.

Gene expression patterns reflect tissue-specific functions. We therefore investigated tissue-specificity of Wfdc15a by comparing it with all proteases and protease inhibitors across 30 mouse tissues using the ENCODE datasets to assess the potential role for spermatogenesis. Unexpectedly, K-means clustering of 842 genes associated with protease or protease inhibitor ontologies revealed that the testis specifically expresses 57% of the WFDC family genes (12 out of 21 genes), which includes Wfdc15a (cluster 1 of Fig. 2A). This is consistent with the evidence suggesting co-evolution and similarity of WFDC family genes, among which 16 are clustered in a 566-kb region of chromosome 2 in the mouse (Clauss et al., 2005; Hurle et al., 2007). Classification showed that cluster 1 contains 100 protease genes and 34 protease inhibitor genes, and revealed that WFDC is the largest family of testis-specific protease inhibitors (Fig. 2A; Fig. S2A-C; Table S1). Of the WFDC family genes, Wfdc15a was the most testis specific in its expression (Fig. 2B). We then asked what cells express Wfdc15a in the testis. Single-cell RNA-seq (scRNA-seq) analysis of adult whole testis (Hermann et al., 2018) demonstrated that Wfdc15a is primarily expressed in spermatocytes and early spermatids [spermatids (ST) 4-6], as well as weakly in Sertoli cells (Fig. 2C). This profile partly overlapped with the spermatocyte marker Sycp3 and was distinct from elongating spermatid (Prm1) or SSC (Zbtb16) markers (Fig. 2C). This is consistent with the bulk RNA-seq data showing strong expression of Wfdc15a in pachytene spermatocytes and round spermatids (Fig. 1G). Analysis of other WFDC family genes showed that Wfdc2, Wfdc3, Wfdc10, Wfdc12 and Eppin are expressed in germ cells, whereas the rest of the genes are primarily expressed in somatic cells, including Sertoli cells and interstitial somatic cells (Fig. S3A). This different expression profile of Wfdc15a relative to other WFDC family members suggests a unique role for this gene in the testis. H3K4me3 levels in GSCs showed that four of the germ cell-expressed genes (Wfdc2, Wfdc10, Wfdc12 and Wfdc15a) are affected in the Kmt2b KO GSCs, suggesting their potential roles in the germline are in a KMT2B-dependent manner (Fig. S3B).

We next asked whether the expression pattern of WFDC family genes is conserved in humans. Sequence analysis using the Genotype-Tissue Expression (GTEx) database revealed that the mRNA of WFDC15D (also known as RP1-300I2.2 and ENSG00000290777 in UCSC and GTEx), which was previously annotated as a pseudogene and one of the human orthologs of Wfdc15a (Clauss et al., 2005), contains an open reading frame (ORF) having a WAP-like domain (Fig. S4A,B). Moreover, the WFDC15D ORF was almost exclusively expressed in the testis (Fig. S4C), suggesting that WFDC15D is not a pseudogene and may encode a functional protein. Further profiling of 19 WFDC family genes in 54 human tissues on the GTEx database revealed that the testis expresses the most WFDC genes relative to other organs (Fig. S4D). These results demonstrate a conserved expression pattern of WFDC family genes between mouse and human testis.

Thus, our data suggest that the WFDC genes constitute the largest family of testicular protease inhibitors in mice and humans and that mouse Wfdc15a is a unique snc-bivalent gene epigenetically programmed for expression in germ cells during late spermatogenesis.

To investigate the functional role of Wfdc15a for spermatogenesis, we performed a CRISPR-Cas9 KO targeted to the protein-coding sequence of the first exon of mouse Wfdc15a. We obtained two independent KO lines of Wfdc15a (#1 and #2) carrying a mutation at the guide RNA (gRNA) target site (Fig. S5A,B). KO #1 contained a stop codon generated by the mutation, whereas #2 carried a 1-bp deletion causing a frameshift mutation. By crossing these mutants with wild-type (WT) animals, we obtained heterozygous (Wfdc15a^+/−) and homozygous (Wfdc15a^−/−) males for downstream analysis. We designed PCR primers that can distinguish WT and the two KO alleles according to the method published previously (Fig. S5A) (Kobayashi et al., 2021). The PCR results confirmed the mutations carried by the KO mice (Fig. S5C).

Crossing of Wfdc15a^−/− males with WT females revealed that both lines of Wfdc15a KO males are infertile (no pups obtained), whereas heterozygous males were not significantly different from WT males (Fig. 3A). Both lines of Wfdc15a mutants exhibited no gross abnormality in adult males and showed a similar body weight to their heterozygous counterparts or WT males (Fig. S5D), but the testis weight declined significantly (47.3% smaller on average compared with heterozygous samples) (Fig. 3B). The testis of the Wfdc15a KO males was also visibly smaller and the epididymal duct of the KO mice was not visible from outside, which is suggestive of a reduced amount of sperm contained (Fig. 3C).

To investigate the cause of infertility, we examined the testicular and epididymal tissues using hematoxylin and eosin (H&E) staining and immunofluorescence (IF). Seminiferous tubule epithelium sections can be classified into 12 stages according to the combinations of cell types contained, and spermatid development is further classified in spermiogenic steps (Fig. S5E) (Meistrich and Hess, 2013). H&E-staining showed that the tubule lumens lack elongating spermatids of approximately step 9 and later, which can be found in the heterozygous seminiferous tubules (Fig. 3D; summarized in Fig. S5F). Notably, stage VIII-XII seminiferous tubules in KO contained abnormal round spermatids clustered in the lumen, which is suggestive of dying cells (Fig. 3D). From this analysis, we observed no abnormal features in all stages of heterozygous seminiferous tubules when compared with WT (Fig. S5G). Consistent with the germ cell degeneration in the KO testis, H&E-staining for epididymis showed that both caput and cauda epididymis from the KO males contained abnormal cells and lacked spermatozoa having thin heads (Fig. S6A,B).

To investigate the defective spermatid development more precisely, we performed IF using markers for different cell types: PLZF (ZBTB16; undifferentiated spermatogonia), KIT (differentiating spermatogonia), SYCP3 (spermatocytes) and TSGA8 (round spermatids) (Fig. S6C-E). This experiment confirmed strong and specific expression of these markers in expected cell types. Moreover, the expected shapes of round spermatid acrosomes stained with peanut agglutinin (PNA) suggested that acrosome formation in round spermatids up to step 8 was not affected (Fig. S6E). Quantification of TSGA8-expressing round spermatids in stage VII-VIII seminiferous tubules showed that the round spermatid number is not significantly different between control and KO (Fig. S6F). Collectively, these results suggest that germ cells can differentiate up to stage VII-VIII round spermatids in the KO testis.

We next performed transmission electron microscopy (TEM) to clarify morphological defects of spermatids. As observed in the IF experiments, round spermatids up to steps 7-8 have normal morphology of nuclei and acrosomes in the KO (Fig. S7A,B). However, many spermatids at approximately steps 8-9 and later contained a large vacuole in the nuclei and were highly deformed (Fig. S7A,B). In the KO testis, we found a very small number of elongating spermatids, which showed dense and malformed nuclei with branched appearances, suggesting that KO spermatids cannot undergo extensive morphological transition after step 9 (Fig. 3E,F; Fig. S7A,B).

To examine whether increased apoptosis is associated with the reduction of elongating spermatids in the KO seminiferous tubules, we performed the terminal deoxynucleotidyl transferase (TdT)-mediated d-UTP nick end-labeling (TUNEL) assay. We counted the number of TUNEL-positive cells inside seminiferous tubules in control and KO testis from adult mice and normalized for the size of the testis. We observed a statistically significant increase in the number of apoptotic cells per tubule in the KO testis compared with the control testis (Fig. S7C,D). Together, these results demonstrate that Wfdc15a regulates germ cell development during the round-to-elongating transition of spermatids at step 9.

Disturbance of apical ES can cause spermatid exfoliation or premature spermiation. We investigated whether the apical ES structure supported by Sertoli cell F-actin is preserved in the KO testis. In the Sertoli cells of control mice, strong actin bundle formation labeled with Phalloidin was observed around spermatid nuclei (Fig. 3G,H). However, in the KO mice, dispersed spreading of F-actin throughout the Sertoli cell cytoplasm was observed at stage IX and later, where developing spermatids were lost (Fig. 3I,J; Fig. S7E). Notably, there was no accumulation of actin bundles around the remaining elongating spermatids. Likewise, diffusion of espin, an actin-bundling protein, was also observed in the KO Sertoli cells (Fig. 3K-N). Meanwhile, preservation of the blood-testis-barrier (BTB), an ultrastructure linked to infertility (Cheng et al., 2011a; Jiang et al., 2014), was confirmed by the co-staining of a BTB-tight junction component CLDN11, with F-actin (Fig. 3G-J; Fig. S7E). These data suggest that the Wfdc15a deficiency in spermatids, and potentially in Sertoli cells, induced changes, particularly in the apical structural organization of Sertoli cells, which may, in turn, promote the loss of anchoring spermatids (Fig. 3O).

To investigate the mechanism of spermatid deficiency, we performed fluorescent activated cell sorting (FACS) to collect round spermatids from control and KO testes (Fig. 4A; Fig. S8A). DNA staining and sorting enriched haploid round spermatids to the average purity of 87.2% (Fig. S8B). We performed RNA-seq using the sorted cells and analyzed gene expression changes. We identified 1634 upregulated and 423 downregulated genes in the KO spermatids relative to control spermatids (P<0.05, log2 fold change >1; Fig. 4B). Gene ontology analysis revealed significant enrichment of upregulated genes in immune and inflammatory pathways (P-value, 3.0e-16 and 1.1e-14, respectively; Fig. 4C). Notably, genes known to be closely associated with an inflammatory response including Tnf, Fos, Atf3, Ctss (Klinngam et al., 2018; Kwon et al., 2015; Ray et al., 2006; van Loo and Bertrand, 2023; Zhou et al., 2022) were significantly upregulated (Fig. 4B).

Transcriptome analysis is a powerful method for capturing protease-antiprotease imbalance (Harbig et al., 2020; Johansen et al., 2022; Verbovšek et al., 2014). Interestingly, our analysis revealed that both upregulated and downregulated terms of biological process and molecular function contained protease-related functions such as protease binding, peptidase activity, peptidase inhibitor activity and zymogen activation (Fig. 4C; Fig. S8C). In total, 81 upregulated and 18 downregulated genes were associated with protease-related functions. To understand the effect of Wfdc15a KO on testicular protease homeostasis, we sought to investigate what kind of proteases and related genes are dysregulated in spermatids. Because little is known about the function and cell-type expression of a large number of proteases and protease inhibitors in the testis, we analyzed whole testis scRNA-seq data from WT mice (Hermann et al., 2018) and mapped the expression patterns of all protease-related genes. We clustered the genes according to the expression patterns in different cell types (Fig. 4D). This analysis revealed highly cell type-dependent expression of 348 testicular protease-related genes expressed in the adult testis and identified six clusters according to the cell-type specificity, suggesting their different physiological roles in maintaining homeostasis. Annotation of Wfdc15a KO-upregulated genes to the clusters showed that 96% (23 genes) are not spermatid-expressed in the normal testis (Fig. 4E), suggesting that these genes were ectopically upregulated in the KO spermatids. Among these, 88% (21 genes) were specific for either spermatogonia or somatic cells including Sertoli cells and peritubular macrophages in the WT testis. Classification of the upregulated testicular protease inhibitors showed that 56% of them (19 genes) are serine-type (Fig. 4E). In contrast, among downregulated, all genes annotated to the single-cell clusters belonged to the kallikrein family of serine-type proteases expressed in spermatids (Klk1, Klk1b8, Klk1b26 and Klk1b3) (Fig. 4D). Of these, Klk1b8, Klk1b26, and Klk1b3 were predominantly expressed in the testis (Fig. S8D). These results unexpectedly demonstrated that deletion of a single protease inhibitor Wfdc15a causes dysregulation of the protease and protease inhibitor network.

Testicular inflammation is a common condition associated with male infertility, involving various cytokines (Hasan et al., 2022; Hedger, 2011). In particular, TNFα, encoded by the Tnf gene, is a cytokine that plays a central role in inflammation (van Loo and Bertrand, 2023). Protein-protein interaction analysis of the top 400 upregulated genes in the KO spermatids demonstrated co-upregulation of Tnf and its interacting partners (Fig. S9A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Kanehisa and Goto, 2000) demonstrated that crucial genes involved in the TNF signaling pathway such as p38 (Mapk14), AP-1 (Jun), IL1B and CXCL family chemokine genes were also upregulated (Fig. S9B). To validate the upregulation of the inflammatory genes at the protein level, we performed IF using an antibody specific for TNFα. This experiment confirmed significant upregulation of TNFα in round spermatids inside seminiferous tubules (Fig. 4F). Staging revealed tubules of stages IX to XII and I, which coincides with the timing of spermatid impairment (Fig. S5G), contained increased numbers of TNFα-positive round spermatids in the KO samples (Fig. 4G; Fig. S9C). Thus, these results suggest that Wfdc15a maintains protease homeostasis in spermatids and highlight its crucial role in innate immunity.

Protease imbalance and inflammatory cytokines from spermatids can act both in autocrine and paracrine manners and may induce broad alteration across the whole testicular tissue, which may in turn cause testicular inflammation (Huleihel and Lunenfeld, 2004; Yan et al., 2008). To investigate this possibility, we analyzed gene expression changes using adult whole testis from Wfdc15a^−/−, Wfdc15a^+/− and WT mice (Fig. S10A). Principal component analysis (PCA) of RNA-seq samples confirmed similarity of WT and heterozygous testes and their clear separation from the KO samples (Fig. S10A). Hierarchical clustering of differentially expressed genes (DEGs) confirmed the consistency of the two KO lines (Fig. S10B). Comparison between control (+/+ and +/−) and KO (−/−) data identified 1126 upregulated and 624 downregulated genes (P<0.05, log2 fold change >1; Fig. S10C). Gene ontology analysis of the upregulated genes revealed enrichment for immune and inflammatory terms including ‘response to cytokine’, ‘immune system process’ and ‘inflammatory response’ (Fig. 5A; Fig. S10D,E). When performing KEGG pathway analysis, we found that upregulated genes were also linked to inflammatory related pathways such as cytokine receptor interactions and infectious conditions (Fig. S10F). Comparison with the spermatid DEGs showed that 65.3% of the genes associated with immune or inflammatory functions (47 genes) were uniquely upregulated in the whole testis (Fig. 5B), suggesting that inflammatory response also occurred outside spermatids. Meanwhile, the downregulated genes were enriched for ‘multicellular organism development’, ‘cell differentiation’ and ‘sperm chromatin condensation’, consistent with the disturbance of spermatid development (Fig. S5G). Other terms included ‘calcium-ion binding’, ‘ion channel activity’ and ‘signaling receptors’, which are linked to male fertility (Fig. 5C; Fig. S10G,H) (Chen et al., 2016). Expression analysis of male germ cell markers for different cell types showed a moderate reduction of early spermatid genes and a more severe reduction of late spermatid genes in the Wfdc15a KO testis (Fig. 5D), which is in line with the histological observations showing a defect at the round-to-elongating spermatid transition at around step 9 (Fig. S5G). Further analysis using scRNA-seq data from whole testis (Hermann et al., 2018) showed that 87% of the downregulated genes are late-spermatid genes (Fig. S10I). These data suggest that inflammatory response also occurred outside spermatids, whereas most of the gene downregulation reflects spermatid defects.

Notably, both upregulated and downregulated terms included ‘peptidase’, ‘peptidase inhibitor’ and ‘zymogen activation’ (Fig. 5A,C; Fig. S10D-H). We asked whether the protease-related genes were affected broadly in the testicular tissue. Comparison of protease-related DEGs between spermatids and whole testis showed that 70.9% of the protease-related DEGs (56 genes) were specifically altered in the whole testis (Fig. 5E), suggesting that these genes were induced outside spermatids. Characterization of the whole testis DEGs revealed changes in 54 genes with protease or associated functions (43, upregulated; 11, downregulated) and 25 genes with protease inhibitory or associated functions (23, upregulated; 2, downregulated) (Fig. 5F,G). Interestingly, KO induced upregulation of three testis non-specific WFDC genes (Wfdc1, Wfdc17, Wfdc18). Classification of the genes indicated that 80.0% of Wfdc15a KO-upregulated protease inhibitors are serine-type (16 genes) (Fig. 5H).

Thus, the RNA-seq comparison of control and KO samples from the adult samples captured events underlying adverse conditions in steady-state spermatogenesis. However, the lack of elongating spermatids only in the KO samples may have affected DEG identification. To validate our finding while excluding any influence of elongating spermatids, we analyzed postnatal day (P) 25 testicular samples which would contain equivalent spermatid cell types between control and KO, as elongating spermatids are absent (Ernst et al., 2019). Our IF analysis using the P25 samples for TSGA8 confirmed that WT, heterozygous and KO seminiferous tubules contained similar types of spermatids: there were abundant round spermatids but elongating spermatids were missing (Fig. S11A). Using these samples, we collected round spermatids and whole testicular tissues and performed RNA-seq. For the round spermatid collection, we used the same FACS-based method applied for adult samples (Fig. 4A) and confirmed the average of 86.8% of purity by microscopy validation (Fig. S11B). Consistent with the adult spermatid RNA-seq data, comparison of control and KO round spermatids from P25 mice showed that upregulated genes were enriched for terms associated with inflammatory response such as ‘apoptotic signaling pathway’ and ‘immune system process’ (Fig. S11C). Likewise, whole testis RNA-seq from P25 mice showed enrichment of upregulated genes in inflammation-related terms such as ‘innate immune system’ and ‘inflammatory response’ (Fig. S11D). To assess the consistency of the protease-related DEGs identified in adult spermatids (Fig. 4D), we analyzed the expression of the same genes using the P25 round spermatid data. Although changes of the downregulated genes were not significant, we observed similar trends of expression changes in P25 spermatids (Fig. S11E). The protease-related DEGs in adult whole testis (Fig. 5F,G) also showed similar trends of changes in P25 KO testes (Fig. S11F). Thus, these results demonstrate that disturbance of protease homeostasis and inflammatory response are occurring in round spermatids from the first wave of spermatogenesis in the absence of Wfdc15a.

To further validate the changes of protease-related genes at the protein level, we performed IF for tissue inhibitor of metalloproteinase 1 (TIMP1), the mRNA of which showed an 8.4-fold increase in the adult KO testis (Fig. 5F,G). We found that the TIMP1 protein was significantly upregulated in seminiferous tubules of stages VII to XII and I (Fig. 5I,J; Fig. S12A), which coincided with the onset of spermatid defects in the KO testis (Fig. S5G). This analysis also showed increased signal intensity of TIMP1 outside seminiferous tubules of the KO samples (Fig. S12A), suggesting that the influence of Wfdc15a loss extends to interstitial cells. Co-staining for TIMP1 and Sertoli cell F-actin with high-resolution microscopy confirmed that TIMP1 upregulation in the KO occurred in the Sertoli cell cytoplasm (Fig. S12B). These results suggest that the protease-antiprotease disturbance at the protein level also occurred in Sertoli cells.

Testicular inflammation has been studied using a well-established experimental autoimmune orchitis (EAO) model (Tung and Teusher, 1995) and a lipopolysaccharide (LPS)-induced inflammation model in rodents (Hedger et al., 2005; Winnall et al., 2011). We investigated whether the conditions we observed in the Wfdc15a KO testis are common with the ones seen in testicular inflammation. We confirmed upregulation of most inflammatory genes including Tnf and Cd14, as well as hypoxia genes such as Angptl4, Egr1 and Ier3, that are activated in mouse EAO (Nicolas et al., 2017) and/or rat LPS (Palladino et al., 2018) in the KO testis (Fig. 5K). Additionally, scRNA-seq data detected the expression of major receptor genes for IL1 and Tnf in Sertoli cells, Leydig cells, and peritubular macrophages in a physiological condition (Fig. S12C). Together, these results suggest that spermatid abnormality caused by the loss of Wfdc15a disturbs protease-antiprotease homeostasis broadly in the seminiferous epithelium and induces a response reminiscent of experimental models of testicular inflammation.

The etiology of ∼40% of male infertility is unknown and the search for causative genes linked to idiopathic infertility remains a challenge (Krausz and Riera-Escamilla, 2018). The genes pre-programmed by KMT2B are generally expressed at low levels in SSCs and activated later during differentiation (Tomizawa et al., 2018). Our epigenome analysis revealed that H3K9me3 represses a subset of these genes in SSCs by non-canonical bivalent chromatin and their activation accompanies H3K9me3 removal and H3K27ac acquisition. Although H3K9me3 in lineage-committed cells typically co-exists with highly methylated CpGs to form constitutive heterochromatin at centromeres, telomeres and retroelements (Nicetto and Zaret, 2019; Nicetto et al., 2019), the H3K9me3 at the snc-bivalent genes forms facultative heterochromatin. A report suggests that facultative H3K9me3 is found at silent protein-coding genes (Saksouk et al., 2015), but their promoter sequence features and DNA methylation have not been well characterized. In this study, we showed that the snc-bivalent genes have CpG-poor and hypomethylated promoters, which contrast with the Polycomb-targeted canonical H3K4me3/H3K27me3 bivalent promoters showing CpG-dense characteristics (Schübeler, 2015). Thus, the choice of whether to use H3K27me3 or H3K9me3 to repress genes in SSCs likely depends on CpG density. We previously found that a new infertility gene, Tsga8, which has the canonical bivalent chromatin in SSCs, is specifically and highly expressed in round spermatids (Kobayashi et al., 2021; Tomizawa et al., 2018). Because the snc-bivalent gene promoters are unmethylated on cytosines, one role of KMT2B may be to prevent the formation of constitutive heterochromatin over a set of promoters by blocking the access of DNA methyltransferases through H3K4me3 (Jia et al., 2007; Ooi et al., 2007), even in the presence of other repressive modifications (Fig. 6A). Taken together, our study suggests that KMT2B modulates a facultative heterochromatin state of genes, particularly for spermatogenesis and development, and disruption of these epigenetic programs could lead to a profound impact on newly generated cells, which may link to male infertility. Thus, our epigenetic approach adds a new path to infertility gene identification.

We showed that Wfdc15a KO results in dysregulation of the protease network and upregulation of inflammatory associated pathways not only in the steady-state spermatogenesis but also in the first wave of spermatogenesis. Because Wfdc15a does not encode a transcription factor, the mRNA changes we observed are likely the indirect effect of the knockout. Wfdc15a may affect various signaling cascades as observed in other proteases and protease inhibitors affecting transcription through MAPK or RAS signaling (Heuberger and Schuepbach, 2019; Sharony et al., 2010). Additionally, Wfdc15a KO may have induced a compensatory mechanism, through feedback loops or genetic compensation, to activate transcription of related genes (El-Brolosy et al., 2019; Sztal and Stainier, 2020), given the activation of serine-type protease inhibitors in KO. The effect of Wfdc15a depletion is consistent with other WFDC family members showing anti-protease and innate immunity properties in other tissues such as lung, skin and gastrointestinal tract (Clauss et al., 2002; Glasgow et al., 2015; Nakajima et al., 2019; Nugteren and Samsom, 2021; Scott et al., 2011; Small et al., 2017). In Wfdc15a-deficient round spermatids where apical-ES is absent and immune homeostasis fails, TLRs including Tlr1 and Tlr2 were upregulated. Also, the pro-inflammatory cytokines/chemokines IL1, Tnf, Ccl2 and Cxcl2, downstream molecules of TLRs, and effector transcription factors AP-1 and Fos were upregulated in the KO round spermatids. These observations suggest that immunological pathways inside the defective round spermatids are substantially activated, which then affects surrounding cells. Additionally, the loss of the weak expression of Wfdc15a mRNA in Sertoli cells or potential secretion of WFDC15A (Bingle and Vyakarnam, 2008) may have influenced the surrounding cells and tissues. TNFα induces the expression of CCL2 in Sertoli cells, Leydig cells and testicular macrophages (Wang et al., 2021). CCL2 may recruit leukocytes and causes a negative impact on spermatogenesis and induces acute epididymo-orchitis (Wu et al., 2021). This could also explain why the inflammatory response that primarily occurred in round spermatids by the loss of Wfdc15a extended immune and inflammatory responses to the entire testicular tissue. In support of this, the receptor genes for IL1 and Tnf were expressed in Sertoli cells, Leydig cells and peritubular macrophages. It is also notable that the KO phenotype was similar to that of the EAO model or LPS-exposed testis. Both are characterized by increased expression of cytokine genes such as Tnf and Ccl12 as well as TLRs that are closely linked to cytokine production, including Tlr1 and Tlr2 (Nicolas et al., 2017; Palladino et al., 2018). The activation of inflammatory pathways may be linked to potential premature spermiation due to the loss of F-actin accumulation at apical ES after stage IX in the Wfdc15a-depleted Sertoli cells. In addition to the TNFα expression, we observed a significant increase of apoptosis in the KO mice. Therefore, the reduced Sertoli-spermatid adhesion in the Wfdc15a KO mice through inflammation and apoptosis could be complicit in the loss of spermatids.

Using a combinatorial analysis of ENCODE expression data, single-cell transcriptome data as well as the mouse KO models, we revealed specialized expression and functions of the WFDC family of protease inhibitors in the testis. Our KO mouse analysis revealed a specific role of a KMT2B-preprogrammed Wfdc15a gene in regulating testicular inflammation during spermatid development (Fig. 6B). Depletion of the Wfdc15a gene alone caused substantial dysregulation of proteases and protease inhibitors as well as activation of inflammatory pathway genes in the testis. Previous studies have shown that genes with similar functions and expression, such as immune genes and reproductive genes, tend to cluster together in the genome (Borgmann et al., 2016; Miller et al., 2004; Williams and Hurst, 2002). WFDC genes are coded in two subclusters in mice and humans (Clauss et al., 2005; Small et al., 2017), suggesting that they co-evolved to acquire similar functions. Indeed, of the 21 WFDC genes in the mouse, 12 showed specific expression in the testis. Interestingly, one of the WFDC subclusters contains six seminal vesicle-secreted protein (SVS) genes, which also play a role in reproduction (Clauss et al., 2005; Noda and Ikawa, 2019). Thus, the WFDC loci can be regarded as large reproductive clusters and may cooperatively play a role in inflammatory regulation for successful reproduction.

Currently, the mechanistic details of how Wfdc15a regulates immune homeostasis and how it is linked to defective spermiogenesis are unclear. However, the specific effect of Wfdc15a KO on four kallikrein genes in round spermatids suggests that Wfdc15a potentially regulates inflammation through the kallikrein-kinin system (Bekassy et al., 2022). The positive outcome of asthenozoospermia and/or oligozoospermia treatment with oral KLK1 administration is compatible with our results (Aydin et al., 1995; Giovenco et al., 2009; Schill and Miska, 2009). Wfdc15a may also be involved in various cellular processes of protease-regulated spermatogenesis, such as sperm maturation, spermiation and fertilization (Kiyozumi and Ikawa, 2022; Mruk and Cheng, 2004; Scovell et al., 2021; Wong et al., 2008). WFDC15D, a human ortholog of Wfdc15a, was previously annotated as a pseudogene (Clauss et al., 2005). Pseudogenes are thought to be nonfunctional due to a disrupted ORF or lack of expression (Bruford, 2014). However, our analysis revealed that WFDC15D contains an ORF and is expressed specifically in the testis. This finding opens the possibility that this gene is not a pseudogene and may play similar roles in humans. An implication of our study is that an imbalance of a complex network of protease-antiprotease involving WFDC genes may underlie inflammatory-associated male infertility in human patients.

All mouse experiments were approved by the Committee for Animal Care and Use at Yokohama City University (approval ID: F-A-20-001, F-A-23-015) and carried out according to their guidelines.

KO mouse generation was performed as previously reported (Aoto et al., 2021). Firstly, three gRNAs designed for Wfdc15a were cloned into pSpCas9 (BB)-2A-GFP (px458, a gift from Feng Zhang, Broad Institute of MIT and Harvard; Addgene plasmid #48138) with U6 promoter-driven empty guide RNA, Streptococcus pyogenes Cas9 and checked for cleavage efficiency in the Neuro2a mouse neuroblastoma cell line. Cells were collected and lysed in 50 mM NaOH/0.1 M Tris-HCl (pH 7.5) and analyzed with the high-resolution melting (HRM) analysis using an Eco 48 Real-Time PCR system (PCRmax). The gRNAs showing the highest efficiency of cleavage were selected for the KO mouse generation. The sequence used was as follows: GCTCTAGGCTGCGCCATGCTG.

For KO mouse generation, fertilized one-cell eggs collected from ICR females (Japan SLC) were incubated in EmbryoMax^® KSOM Mouse Embryo Medium (MR-121-D, Merck) until electroporation. CRISPR ribonucleoprotein (RNP) containing 0.1 μg/μl Cas9 protein, 0.1 μg/μl crRNA, 0.1 μg/μl tracrRNA (Integrated DNA Technologies) and 0.2 μg/μl single-stranded oligodeoxynucleotides (ssODN; Macrogen) within 130 bp homologies were mixed in Opti-MEM solution (31985062, Thermo Fisher Scientific). Collected eggs (10-15 cells) with 5 μl CRISPR RNP/ssODN solution were placed in a 1 mm gap electrode (CUY501P1-1.5, NEPAGENE), and electroporated (NEPA21 electroporator, NEPAGENE). The conditions were: poring pulse (40 V, pulse duration 3.5 ms, pulse interval 50 ms, four pulses) and transfer pulse (5 V, 50 ms, 50 ms, five pulses). The eggs were cultured for 1 day and 55 two-cell-stage embryos were transferred to the ampulla of the oviduct of pseudopregnant ICR recipients under the use of the following anesthesia: Domitor/Midazolam/Betorufare (Nakakita Yakuhin); Antisedan, a recovery drug (Nakakita Yakuhin). The sequence used as ssODN was: ATGAAGCCAAGCAGCCTCCTACTGTTCACAACAACCATCCTCCTTTGCCTCAGCtaGGCGCAGCCTAGAGCAACTAGGAAAGGAG. Of the eight pups obtained by embryo transfer, two were confirmed to have a mutation by genomic DNA sequencing. Detailed analysis of the KO phenotype was performed using the progenies of mutant #1 obtained by crossing heterozygous mutants unless otherwise indicated, because no obvious difference was observed in histological and RNA-seq data between #1 and #2.

PCR primers were designed to distinguish the small differences generated by the KO as previously reported (Kobayashi et al., 2021). This method is applied from the detection of single nucleotide substitutions (Mamotte, 2006) whereby the 3′ end of one of the two primers used for a single PCR is designed to overlap with the mutation site. Such primers amplify DNA only when the 3′ end of the primer is complementary to the mutation site. The primer pairs were: WT forward, ATCCTCCTTTGCCTCAGCAT; KO #1 specific forward, CATCCTCCTTTGCCTCAGCT; KO #2 specific forward, ATCCTCCTTTGCCTCAGCC; common reverse, GGTTCCACAGAATTCTCAGACG. The PCR condition was 30 cycles of 94°C for 30 s, 60°C for 30 s and 72°C for 30 s, followed by a 7 min extension at 72°C with GoTaq DNA polymerase (Promega), which lacks the 5′ to 3′ exonuclease activity.

Each adult male mouse from the Wfdc15a lines was paired with two WT ICR females for 1 week, and the number of pups obtained from each female was counted and compared between control and KO.

Preparation of testicular frozen blocks and staining of frozen sections were performed as described previously (Ohmura et al., 2004; Shirakawa et al., 2013). Perfusion-fixation of 10- to 15-week-old mice was performed with 4% paraformaldehyde (PFA) and the testes were incubated in 4% PFA for 1 h per 25 mg of testis. The tissues were embedded in OCT compound (Sakura Finetek) and frozen sections were made for staining. Samples were incubated in 1% bovine serum albumin in PBS and incubated for either 2 h at room temperature or 4°C overnight with an appropriate primary antibody. After washing, a diluted secondary antibody was applied and incubated for 1 h at room temperature.

The following antibodies were used: goat anti-PLZF (1:200; R&D Systems, AF2944, lot# VUG0319041), rat anti-KIT (1:1000; Invitrogen, MA5-17836, lot# TE2578616), rabbit anti-SYCP3 (1:750; NovusBio, NB300-231, lot# F-3), rabbit anti-TSGA8 (Uchida et al., 2000; 1:800), FITC-conjugated anti-PNA (1:100; EY Laboratories, F2301, lot# 030917), rabbit anti-CLDN11 (1:100; Thermo Fisher Scientific, 36-4500, lot# XI362819), rabbit anti-espin (1:100; Abcam, ab254774, lot# 1056683-1), goat anti-TIMP1 (1:500; R&D Systems, AF980-SP, lot# CNL0421041), mouse anti-TNFα (1:500; Santa Cruz Biotechnology, sc-52746, lot# F0414), donkey anti-rabbit IgG Alexa555 (1:200; Thermo Fisher Scientific, A31572, lot# 1818686), biotin-SP donkey anti-rat IgG (1:1000; Jackson ImmunoResearch, 712-066-153, lot# 104189), donkey anti-rabbit IgG Alexa488 (1:200; Thermo Fisher Scientific, A21206, lot# 2072687) and Alexa546 Phalloidin (1 unit/ml; Thermo Fisher Scientific, A22283, lot# 727953). The TSA Biotin System (Perkin Elmer) was used to amplify signals of KIT, TIMP1 and TNFα. HistoVT One (Nacalai Tesque) was used for antigen retrieval in the CLDN11 staining experiments. DAPI (Sigma-Aldrich) was used for nuclear staining. The sections were mounted with ProLong Diamond (Thermo Fisher Scientific) and observed using confocal laser microscopy (Olympus FV-1000). Staging of seminiferous epithelium was performed as previously reported (Meistrich and Hess, 2013). We also used an acrosome marker PNA, which helps to determine precise stages by the spermatid morphology. In the Wfdc15a KO testis, stages IX-XII were grouped together due to the difficulty of staging in the absence of elongating spermatids.

The TUNEL assay was performed using ApopTag Plus Fluorescein In Situ Apoptosis Detection Kit (Chemicon International, S7111) according to the manufacturer's instructions. TUNEL-positive signals of >10 μm² within seminiferous tubules were detected and counted by using the National Institutes of Health ImageJ software.

For H&E staining, testes collected from 10-week-old mice were fixed in Bouin's solution at 4°C overnight followed by a transfer to 70% ethanol and embedding in paraffin. The blocks were sliced at 6 µm with a Microm HM355 (Thermo Fisher Scientific) and imaged with Keyence BZ-X800. Staging of Wfdc15a control tubules was performed as previously reported (Meistrich and Hess, 2013). Tubules of Wfdc15a KO mice were grouped into stages I-VIII or VIII-XII according to the morphology of germ cells up to round spermatids.

Sample preparation for transmission electron microscopy (TEM) and image analysis were performed as previously described (Nakajima et al., 2019). Briefly, 10-week-old males were perfusion-fixed using phosphate buffer (0.1 M, pH 7.4) containing 2% glutaraldehyde and 2% paraformaldehyde. Small pieces of tissues were gathered from the fixed testes and placed in the same fixative for 2 h at room temperature. Tissues were then washed with 0.1 M phosphate buffer and fixed again in 1% OsO₄ for 1 h at 4°C. Fixed tissues were washed and incubated in 4% uranyl acetate for 1 h at room temperature. After dehydration in graded concentrations of ethanol, the samples were placed in propylene oxide, embedded in Epon812 resin (TAAB Laboratories Equipment) and incubated for 48 h at 60°C. Ultra-thin sections were made using Reichert Ultracut N Ultramicrotome (Leica Microsystems) and stained with 2% uranyl acetate in 70% ethanol and 0.4% lead citrate. TEM imaging was performed using JEM1400 flash (JEOL).

Round spermatids from 15-week-old adults and P25 mice were collected according to the method previously reported (Lesch et al., 2019) with modifications as follows. Dissociated seminiferous tubules from each adult male at 9-10 weeks of age were treated with 1 mg/ml of collagenase (Sigma-Aldrich) and 1 unit/ml of DNase I (Promega) in 1 ml of PBS at 32°C for 10 min using a shaker incubator. The reaction was quenched with 10 ml of Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% fetal calf serum (FCS) and the cells were passed through a 100 μm strainer (Falcon). After centrifugation (400 g for 5 min), cells were re-suspended to the concentration of ∼5×10⁶ cells/ml in 5% FCS/PBS. Cells were incubated with Vybrant DyeCycle Green (Thermo Fisher Scientific; 1:500) at 37°C for 30 min. Dead cells were labeled with Propidium Iodide (Sigma-Aldrich) at a concentration of 2 μg/ml. Cells with 1C DNA content and large size (round spermatids) were collected using the SONY MA900 cell sorter. The purity of the sorted cell population was verified after DAPI staining with a fluorescent microscope (Olympus BX51).

Whole testes were collected from 15-week-old adult and P25 mice. Tunica albuginea was removed and the whole testicular tissues were homogenized to isolate total RNA using Isogen (Nippon Gene) according to the manufacturer's instructions. About 2 ng of the purified total RNA was used for RNA-seq library preparation using the Smart-seq2 protocol (Picelli et al., 2013) with 11 cycles of PCR. For the generation of libraries from sorted round spermatids, 100-200 cells were lysed in 2.5 μl of RLT plus buffer (Qiagen) supplemented with 0.06 μl SUPERase-In RNase inhibitor (Thermo Fisher Scientific) per sample. The samples were directly processed according to the Smart-seq2 protocol with 15 cycles of PCR. The amplified cDNA (200 pg) was used for the Tagmentation reaction using Nextera XT DNA Sample Preparation Kit (Illumina) with dual index primers from Nextera XT Index Kit v2 (Illumina) according to the manufacturer's instructions. The barcoded libraries were pooled and sequenced using either Illumina HiSeq 2500 or NovaSeq 6000.

For the tissue-specific expression analysis of proteases and related genes, we obtained fastq files of RNA-seq from 30 mouse tissues. Accession codes for each tissue are summarized in Table S2. To check for consistency, we used three independent RNA-seq datasets from adult mouse testis including the one we published previously (Kobayashi et al., 2021). The list of the genes associated with protease and related functions was obtained from the Mouse Genome Informatics (MGI) gene ontology database using ‘peptidase’ and ‘peptidase inhibitor’ as keywords. The heat map of human WFDC family gene expression was generated using the GTEx Multi Gene Query. For the analysis of a potential Wfdc15a ortholog in humans, tissue expression data for WFDC15D (RP1-300I2.2) were obtained from the GTEx Analysis Release V8 (dbGaP Accession phs000424.v8.p2; https://gtexportal.org/home/). The corresponding cDNA sequence for this gene (ENST0000s0461248.1) was obtained from Ensembl (GRCh38.p13 dataset; http://www.ensembl.org) (Cunningham et al., 2022) and NCBI GenBank (EST AI222267; https://www.ncbi.nlm.nih.gov/).

For gene expression analysis from fastq files, RNA-seq reads were quality-checked with FastQC (version 0.11.8, Babraham Bioinformatics), trimmed using Trim Galore! (version 0.6.6, Babraham Bioinformatics) and mapped with Hisat2 (version 2.1.0) (Kim et al., 2015) on the mm10 mouse genome using default parameters. Read count was performed by using SeqMonk (version 1.48.1, Babraham Bioinformatics), and differential gene expression analysis was performed using edgeR (version 3.36.0) (Robinson et al., 2010). Volcano plots were generated using a Bioconductor package EnhancedVolcano (version 1.22.0). ChIP-seq reads were mapped onto mouse genome mm10 using Bowtie (version 1.1.1) (Langmead et al., 2009) with the parameters -n 2 -l 30 --best --strata -m 1. ChIP-seq peaks were identified using MACS2 (Zhang et al., 2008). Histone modification enrichment analysis was performed using ngs.plot (Shen et al., 2014). Gene ontology analysis was performed using DAVID 2021 (Huang et al., 2009a,b). For the DAVID analysis, genes showing more than a 2-fold change in expression with P-values <0.05 in edgeR were used unless otherwise indicated. DNA methylation analysis of P7 spermatogonia was performed using published whole genome bisulfite sequencing data and the methods reported (DRA002477 and DRA002402) (Kubo et al., 2015).

scRNA-seq data from adult whole testis (GSE109033) (Hermann et al., 2018) were downloaded and analyzed using Seurat version 4.0.5 (Hao et al., 2021). After creating a Seurat object with min.cells=3, min.features=200, quality filtering was performed by counts of RNA and mitochondrial reads. Cluster identification was performed with UMAP.

The number of samples used for experiments are indicated in figure legends. Data are mean±s.e.m. One-sided Wilcoxon rank-sum test was used to calculate statistical significance. RNA-seq samples having excessive numbers of reads derived from PCR duplication were excluded from analysis. No randomization was used, and investigators were aware of sample group allocation.

We thank Dr Atsushi Toyoda for technical support in sequencing experiments and Dr Bluma Lesch for valuable advice on round spermatid sorting experiments.