Introduction

Spermatogonial cells (SCs) are the initiators and supporting cellular foundation of spermatogenesis in testis in many species including mammals. In mice, SCs become active one to two days after birth, when they exit mitotic arrest and start dividing to populate the basement membrane of seminiferous tubules. During the first week of postnatal life, a population of SCs continues to proliferate to give rise to undifferentiated A_single (A_s), A_paired (A_pr) and A_aligned (A_al) cells. The remaining SCs differentiate to form chains of daughter cells that become primary and secondary spermatocytes around postnatal day (PND) 10 to 12. Spermatocytes then undergo meiosis and give rise to haploid spermatids that develop into spermatozoa. Spermatozoa are then released in the lumen of seminiferous tubules and continue to mature in the epididymis until becoming capable of fertilization by PND 42-48 (De Rooij, 2017; Kubota and Brinster, 2018; Oatley and Griswold, 2017).

Recent work showed that SCs in early postnatal life have distinct transcriptional signatures (Green et al., 2018; Hammoud et al., 2014; Hermann et al., 2018; Law et al., 2019). During the first week of postnatal development, SCs display unique features necessary for the rapid establishment and expansion of the cell population along the basement membrane. This includes high expression of genes involved in cell cycle regulation, stem cell proliferation, transcription and RNA processing (Grive et al., 2019). In comparison, genes expressed in adult SCs are involved in the maintenance of a balance between proliferation and differentiation, and help constitute a steady cell population that ensures sperm formation across life. Previous transcriptome analyses have revealed that adult SCs prioritize pathways related to paracrine signaling and niche communication, as well as mitochondrial functions and oxidative phosphorylation (Grive et al., 2019; Hermann et al., 2018). Epigenetic changes such as histone tail modifications and DNA methylation have also been reported in SCs during postnatal development (Hammoud et al., 2014; Hammoud et al., 2015). Today however, the relationship between the transcriptome dynamics and chromatin landscape in SCs during the transition from early postnatal to adult stage has not been fully characterized.

We characterized the transcriptome and the profile of chromatin accessibility in SCs during the transition from early postnatal to adult stages. The results reveal extensive changes in transcription and in chromatin accessibility in particular, an increase in chromatin accessibility at enhancer regions in adult SCs compared to postnatal SCs. Regions with changes in chromatin accessibility are enriched in binding motifs of several TFs that are differentially expressed between postnatal and adult stages, suggesting a possible role in changes in chromatin accessibility. Analyses of chromatin accessibility at transposable elements (TEs) identified previously uncharacterized changes at long terminal repeats (LTR) and LINE L1 subtypes between developing and adult SCs. Together, these findings suggest a functional link between transcriptional dynamics and chromatin accessibility in SCs during development, and underscore the plasticity of genome organization in germ cells.

Results

FACS enriches SCs collected from postnatal and adult mouse testis

We collected testes from 8- and 15-day old pups (PND8 and 15) and adult males and prepared cell suspensions from each animal by enzymatic digestion. The preparations were enriched for SCs by fluorescence-activated cell sorting (FACS) using specific surface markers (Kubota et al., 2004) (Fig. 1A, B). The purity of sorted cells was evaluated by immunocytochemistry using PLZF (ZBTB16), a well-established marker for undifferentiated SCs (Costoya et al., 2004). FACS enriched cell populations from 3-6% PLZF+ (before FACS) to 85-95% PLZF+ (Fig. 1B), suggesting a high SCs enrichment.

Figure 1.

To validate the molecular identity of enriched SCs populations, we profiled their transcriptome at PND8, PND15 and adulthood by total RNA sequencing (RNA-seq) (n=6 for each group) and examined the expression of known spermatogonial stem cell and somatic markers (Leydig and Sertoli cells) (Fig. 1C and Fig. S1A, B). Classical SCs markers such as c-kit (Schrans-Stassen et al., 1999), Id4 (Helsel et al., 2017; Sun et al., 2015), Lin28a (Chakraborty et al., 2014; Wang et al., 2020), Zbtb16 (Costoya et al., 2004; Song et al., 2020) and Gfra1 (He et al., 2007) were robustly expressed in all SCs samples at each developmental stage (Fig. 1C and Fig. S1A, B), indicating high purity of sorted cells. In contrast, low to negligible expression of somatic cells markers such as Vim (Bernardino et al., 2018) and Tspan17 (Gewiss et al., 2021) for Sertoli cells, and Fabp3 and Hsd3b1 for Leydig cells (Sararols et al., 2021) (Fig. 1C) was detected. To further validate the samples purity, we manually curated a list of spermatogonial and somatic cell markers derived from recent single-cell RNA-seq datasets (Cao et al., 2021; Sararols et al., 2021) and determined their expression level in our datasets. Again, we detected robust and consistent expression of all reported spermatogonial markers and low expression of somatic cells markers in all samples at each developmental stage (Fig. S1A). These results confirm the efficiency of our FACS-based enrichment method.

Dynamic transcriptomic states characterize postnatal and adult SCs

We examined the transcriptome dynamics of SCs from postnatal to adult stages by identifying differentially expressed genes (DEGs) in the RNA-seq datasets. We used a stringent cut-off (adjusted-P≤0.05, abs Log₂ fold change (FC)≥1) on 17,000 genes expressed during at least one developmental stage (See Methods). A total of 663 DEGs were identified between PND8 and PND15 (146 down-regulated and 517 up-regulated) and 2483 DEGs between PND15 and adult stage (914 down-regulated and 1579 up-regulated) (Fig. 1D-E and Table S1). Consistent with previous reports (Hammoud et al., 2015), we observed a dynamic regulation of germ cell factors, transcription factors (TF) involved in core pluripotency pathways and signaling molecules important for self-renewal (Fig. S1B). For instance, Pou5f1 (Oct4), a TF necessary for pluripotency (Nichols et al., 1998), is significantly down-regulated in adult SCs while the TFs Klf4 and Sox2, also needed for pluripotency (An et al., 2019), are expressed similarly in postnatal and adult stages although at different level (i.e Klf4 is expressed more than Sox2) (Fig. S1B).

We conducted Gene Ontology (GO) enrichment analyses to identify the biological processes which DEGs are involved in. GO analyses showed that DEGs between PND8 and PN15 are involved in cell differentiation, cilium movement, sperm motility and transcription and include for instance, TFs such as Junb, Hoxb6, Cebpa and Pparg (Fig. 1F and Table S2). Further, genes coding for histone proteins such as histone H2B Hist2h3b and epigenetic modifiers like the methyltransferase Pygo1 and histone acetyl-transferase Ncoa1 are also differentially expressed between PND8 and 15. Interestingly, genes down- or up-regulated between postnatal stages are involved in different processes. While down-regulated genes are implicated in cell differentiation, cell migration and regulation of proliferation for instance, Cdkn1a that is down-regulated at PND15 regulates genetic diversity during spermatogenesis (Kanatsu-Shinohara et al., 2022), up-regulated genes are rather involved in germ cell development, cell signalling, insulin secretion and transcription (Fig. S1C and Table S2). Further, DEGs between PND15 and adulthood have roles in reproduction, spermatogenesis, cell differentiation, DNA replication, glycolysis and extracellular matrix (ECM)-receptor interaction pathways (Fig. 1G, Fig. S1D and Table S2). For example, the DNA methyltransferase Dnmt1, necessary for genome regulation (Edwards et al., 2017) and Col4a2, a subunit of type IV collagen and component of the basal membrane (Reissig et al., 2019), are specifically down-regulated in adult SCs (Fig. S1C). Interestingly, expression of collagen has been associated to a high proliferative potential and the ability to form germ cell colonies in SCs (He et al., 2005), suggesting that regulation of collagen genes in adult SCs may decrease germ-stem cell potential. Genes up-regulated in adult SCs are involved in cilium movement, germ cell development and glycolysis (Table S2).

We then examined the differences and similarities of transcriptional profiles across the three developmental stages. To identify unique or common changes in gene expression during the transition from early postnatal to adult stages, we compared the list of DEGs at PND8 versus PND15 and at PND15 versus adulthood. Remarkably, the vast majority of DEGs show significant stage-specific changes in transcription (Fig. S2A). For instance, 75% (495/663) of DEGs between PND8 and PND15 are not statistically changed in the transition to adult SCs (Fig. S2B and Table S1 and S2). Similarly, 93% (2325/2493) of DEGs in adult SCs are not statistically changed when compared to PND8 versus PND15 SCs (Fig. 1H and Table S1 and S2). GO enrichment analyses showed that adult-specific down-regulated genes are involved in cell differentiation, cell migration, regulation of signal transduction and Wnt signalling, regulation of DNA replication and transcription as well as collagen biosynthesis and laminin interactions while adult-specific up-regulated genes are involved in sexual reproduction, male gamete generation, germ cell development, cilium movement and glycolysis (Fig. 1H and Table S1 and S2).

Interestingly, a small fraction of all DEGs (4.9%) are detected as significantly changed in the same direction and magnitude (adjusted-P≤0.05, absLog₂FC≥1) across the three developmental stages (Fig. S2C and Table S1 and S2). For instance, 121 genes are consistently up-regulated from PND8 to adult stage, while just 26 are consistently down-regulated from PND8 to adult stage (Fig. S2C). Such DEGs are involved in different biological pathways like chromatin organization (Table S2). In particular, the histone gene clusters displayed significant down-regulation across postnatal development and in adulthood (Table S2).

Landscape of chromatin accessibility in SCs during postnatal development

Cellular differentiation is generally accompanied by changes in chromatin accessibility at regulatory elements (Atlasi et al., 2017). We examined how chromatin accessibility in SCs is modified during development using omni-ATAC-seq (Corces et al., 2017). We focused on SCs at PND15 and adulthood, stages that showed the highest changes in gene expression (Fig. 1). Accessible regions were identified by peak-calling on merged nucleosome-free fragments (NFF) as proxy for genomic regions with potential regulatory activity. We identified 158,977 peaks with clear ATAC-seq signal compared with surrounding genomic regions (Fig. S3A and Table S3). Most accessible regions were located at distal intergenic regions (38%), introns (28%) and putative promoter regions (23%) (Fig. S3B, C) and encompassed sequences enriched for histone PTMs associated with active chromatin such as H3K4 mono-, di- and tri-methylation (Fig. S2D). Notably, signal enrichment was higher for H3K4 methylation than for other histone PTMs (Fig. S3D). These results are consistent with previous observations that ATAC-seq peaks are identified at regulatory elements such as enhancers and promoters and are preferentially located at genomic regions with nucleosomes carrying H3K4me (Henikoff et al., 2020).

We then examined differences in chromatin accessibility between PND15 and adult SCs by differential accessibility analysis. 3212 differentially accessible regions (DARs) were identified (adjusted-P≤0.05, absLog₂FC≥1) with a total of 760 regions with decreased chromatin accessibility (DARs-down) and 2452 regions with increased accessibility (DARs-up) in adult SCs (Figure 2A, B and Table S3). DARs were predominantly localized in distal intergenic regions (54% DARs-down, 37% DARs-up) and introns (33% DARs-down, 36% DARs-up) with a minority located in putative promoter regions (8% DARs-down, 12% DARs-up), consistent with the genomic distribution of all detected ATAC-seq peaks (Fig. 2C). GO analyses of the closest genes assigned to each DAR showed that these genes are involved in male gonad development, cell adhesion, sex differentiation and regulation of cell communication among others (Fig. 2D, E).

Figure 2.

Epigenomic annotation of DARs using high-quality publicly available ChIP-seq datasets for postnatal SCs (Cheng et al., 2020) revealed that DARs are predominantly located in regions enriched for histone PTMs associated with enhancers and promoters. 51% of regions with increased accessibility are located at genomic loci that carry histone PTMs since early postnatal stages, and about half (48%) overlap with regions significantly enriched in H3K4me1, a histone PTM associated to enhancers. In contrast, 80% of regions with decreased accessibility are located in regions with H3K4me1 and 1/3 (33%) also overlap with the repressive histone PTM, H3K27me3 (Fig. 2F). Interestingly, while 16% of regions with increased accessibility are located at potential active regulatory elements that carry H3K27ac, none of the regions with decreased accessibility overlap with H3K27ac (Fig. 2F). Therefore, our data indicate that the transition from postnatal to adult SCs is accompanied by discrete, yet robust, changes in chromatin accessibility at potential regulatory elements, suggesting their involvement in the control of gene transcription.

DARs are associated with binding sites for distinct families of TFs

We examined the relationship between differences in chromatin accessibility and TF binding by TF motif analysis using DARs as input. We observed that regions with decreased accessibility in adult SCs are enriched for binding motifs of members of specific TF families such as KLF (KLF2, KLF5, KLF11, KLF12, KLF15, KLF16), SP (SP1-5, SP9), FOXO (FOXO1, FOXO3), ETS (ETV1-6) among others (Fig. 3A and

Figure 3.

Table S4). Expression analyses showed that many of these TFs are differentially expressed in postnatal or adult SCs (adjusted-P<0.05, 10 down-regulated and 11 up-regulated) (Fig. 3B, C). For example, KLF11,12 and 15, TFs that regulate stem cell maintenance and development (Bialkowska et al., 2017), are upregulated in adult SCs compared to PND15 (Fig. 3C). These TFs and their motif match top enriched motifs detected in DARs in adult SCs (Fig. 3A). The SP family of TFs also show differential expression and enrichment of binding motifs in regions with decreased accessibility. SP5, which has its lowest expression in adult SCs (Fig. 3B, C), promotes self-renewal of mESCs by directly regulating Nanog (Tang et al., 2017) .

Members of the FOXO family of TFs also have an enrichment of motifs in regions with decreased accessibility. While FOXO3 has increased expression from PND15 onwards, FOXO1 is decreased in adult compared to postnatal SCs (Fig. 3B, C). Interestingly, FOXO1 and FOXO3 regulate spermatogonial stem cell function and maintenance in mouse (Goertz et al., 2011). Finally, ETS-type of TFs also show differential expression during SCs development and their motif is enriched in regions with decreased accessibility. ELK4 is up-regulated in adult SCs and can act as a transcriptional repressor via recruitment of the HDAC Sirt7 and deacetylation of H3K18ac (Fig. 3B-D) (Barber et al., 2012). EKL4 can also act as transcriptional activator of immediate early genes such as c-Fos (Dalton and Treisman, 1992). Interestingly, c-Fos itself codes for a TF important for the regulation of proliferation (He et al., 2008) and shows a trend towards up-regulation in adult SCs (Table S1). In contrast, GABPA is downregulated in adult SCs compared to PND15, and has been involved in the regulation of proliferation in mESCs (Ueda et al., 2017).

Regions with decreased chromatin accessibility also show an enrichment for the binding motif of the TF DMRT1. Dmrt1 is progressively repressed during SCs postnatal development and is the lowest in adult SCs (Fig 3B, C and E). DMRT1 can act either as a repressor or activator and controls testis development and male germ cell proliferation (Zhang et al., 2016). DMRT1 can inhibit meiosis and promote mitosis in SCs by repressing Stra8 (Matson et al., 2010). Consistently, we observed transcriptional repression of Stra8 in adult SCs (Table S1).

Next, we identified motifs overrepresented in regions with increased chromatin accessibility in adult SCs. The identified motifs correspond to binding motifs of members of specific TF families such as POU (POU2F2, POU1F1, POU5F1), RFX (RXF1-7), DMRT (DMRT1, DMRT3, DMRTA1-2, DMRTC2), SOX (SOX2, SOX3, SOX5, SOX13), NFY (NFYA, NFYB, NFYC) and AP-1 (JUN, JUNB, JUND, FOS, ATF3, BATF) (Fig. 3F and Table S4). A subset of them is also differentially expressed (adjusted-P<0.05, 13 down-regulated and 5 up-regulated). The most overrepresented motif is similar to the binding motif of the POU family of TFs, which are critical regulators of stem cells (Nichols et al., 1998). Pou1f1 and Pou5f1 are transcriptionally repressed in adult SCs while Pou2f2 is maximally expressed in PND15 SCs and down-regulated in adult cells (Fig. 3G, H). We also identified motifs for members of the RFX family of TFs, which are master regulators of ciliogenesis (Choksi et al., 2014) implicated in regulation of neural stem cells (Kawase et al., 2014). RFX2 is robustly expressed in adult SCs (Fig. 3G, H and I) and has been reported to induce the expression of ciliary genes in association with the TF FOXJ, which has a trend towards up-regulation in adult SCs (Fig. S4) (Quigley and Kintner, 2017). Interestingly, ciliary genes are among the top genes specifically up-regulated in adult SCs (Fig. 1H), suggesting a regulatory relationship between RFX2 and ciliary genes expression in adult SCs.

Other binding motifs enriched at regions with increased accessibility correspond to members of the NF-Y complex, NF-YA, NF-YB and NF-YC. In mESCs, NF-Y TFs facilitate a permissive chromatin conformation and play an important role in the expression of core ESC pluripotency genes (Oldfield et al., 2014). NF-YA/B motif enrichment has also been found in regions of open chromatin in human SCs (Guo et al., 2017). Interestingly, the expression of NF-YA and NF-YC is progressively down-regulated starting at PND15 (Fig. 3G, H). We also detected an enrichment for the binding motifs of members of the AP-1 family of TFs, which are involved in many processes, from regulation of cell proliferation to differentiation and acute responses to environmental clues (Bejjani et al., 2019; Vierbuchen et al., 2017). Interestingly, except for Fosl2, other TFs do not show any significant change in transcription in postnatal or adult SCs and are either constitutively expressed or repressed. For instance, c-FOS and JUN are robustly expressed in postnatal and adult SCs (Fig. 3J) consistent with their role in promoting the proliferative potential of spermatogonial stem cells (He et al., 2008; Wang et al., 2018). Together, these data reveal that a shift in TFs repertoire accompanies chromatin reorganization occurring in SCs during postnatal to adult development.

A subset of DARs is associated with differential gene expression

Thus far, our results show that in SCs, the transition from postnatal to the adult stages is accompanied by changes in gene expression and chromatin accessibility. In some instances, differences in gene expression can be correlated with differences in chromatin accessibility of regulatory elements. Since DARs in postnatal and adult SCs overlap with putative enhancer and promoter elements, we examined if these changes in chromatin accessibility correlate with changes in gene expression. First, we tested if promoter regions of DEGs (+/-2.5kb from TSS) have differential chromatin accessibility by identifying all possible ATAC-seq peaks located in promoters and calculating an average accessibility signal. No significant difference was detected for down-regulated genes in adult SCs, despite a trend for decreased accessibility around the TSS (Fig. 4A). However, down-regulated genes had a significant increase in chromatin accessibility in their body (Fig. 4B). Up-regulated genes had a significant increase in chromatin accessibility both at the promoter region and gene body (Fig. 4A, B). Then, we conducted motif enrichment analyses and identified families of TFs with overrepresented binding motifs at the promoter of DEGs. Binding motifs for members of KLF and SP family of TFs were overrepresented in the promoter of both down- and up-regulated DEGs, while motifs for members of TF families SOX and NFY were detected only in the promoter of down-regulated genes. Promoters of up-regulated genes were specifically enriched in motifs for RFX and AP-1 (Fig. 4C). These observations suggest that changes in chromatin accessibility, in particular at the promoter of up-regulated DEGs may be associated with differential binding of RFX and AP-1 TFs.

Figure 4.

Next, we examined the overlap between DARs and the promoter of DEGs. Unexpectedly, only 3.2% of all promoters of DEGs overlap with at least one DAR (Fig. 4D). The majority of overlapping DARs (71%) is associated with promoters of up-regulated genes in adult SCs (Fig. 4D). For example, Gpx4, a peroxidase involved in the control of stemness (Hu et al., 2021), is transcriptionally active in adult SCs and has an open chromatin at its promoter region (Fig. 4E). Since the majority of DARs are located at distal intergenic regions and have histone PTMs associated with active regulatory elements, we extended our search of their target genes by assigning a large region (+/-100kb) around each DEG. We observed that a third of DEGs were associated with at least one DAR (Fig. 4F), with the majority of DARs (69%) overlapping regulatory elements of up-regulated genes. For instance, Braf is transcriptionally up-regulated in adult SCs, which correlates with increased chromatin accessibility at distal regulatory elements (Fig. 4G). Our results overall show that only a fraction of DEGs is associated with DARs, suggesting that changes in gene expression are not always mirrored at the level of chromatin accessibility, consistent with recent observations (Kiani et al., 2022).

Chromatin accessibility at transposable elements (TEs) undergoes remodeling during the transition from early postnatal to adulthood in SCs

TEs are tightly regulated in the germline by coordinated epigenetic mechanisms involving DNA methylation, chromatin silencing and PIWI proteins-piRNA pathway (Deniz et al., 2019). SCs were recently shown to have a unique landscape of chromatin accessibility at long-terminal repeats (LTRs) retrotransposons, compared with other stages of germ cells in testis (Sakashita et al., 2020). We examined chromatin accessibility at TEs in PND15 and adult SCs by quantifying ATAC-seq reads overlapping TEs defined by UCSC RepeatMasker then analyzing differential accessibility at subtypes level (see Methods section). These analyses showed that SCs transitioning from PND15 to adult stages have significantly different chromatin accessibility at 135 TEs subtypes (adjusted-P≤0.05, absLog₂FC≥1) (Fig. 5A, B and Table S5). Most TEs subtypes have decreased chromatin accessibility (93/135, 68,9%) (Fig. 5A) and 42 have increased accessibility (Fig. 5B) in adult SCs. Differentially accessible TEs loci are predominantly located in intergenic (68%) and intronic (25%) regions and 6% re close to a gene (+/-1kb from TSS) (Fig. 5C).

Figure 5.

LTR retrotransposons were the most abundant TEs with altered chromatin accessibility, specifically ERVK and ERV1 subtypes (Fig. 5A, B). Exemplary ERVK subtypes with less accessible chromatin included RLTR17, RLTR9A3, RLTR12B and RMER17B (Table S5). Enrichment of RLTR17 and RLTR9 repeats was previously reported in mESCs, specifically at TFs important for pluripotency maintenance such as Oct4 and Nanog (Fort et al., 2014). Interestingly, we identified the promoter region of the long non-coding RNA (lncRNA) Lncenc1, an important regulator of pluripotency in mESCs (Fort et al., 2014; Sun et al., 2018), to harbor several LTR loci with decreased accessibility in adult SCs. One of these loci, RLTR17, overlap with the TSS of Lncenc1, and its decreased accessibility correlates with markedly lower Lncenc1 expression in adult SCs (Fig. 5D). Lncenc1 (also known as Platr18) is part of the pluripotency-associated transcript (Platr) family of lncRNAs recently identified as potential regulators of pluripotency-associated genes Oct4, Nanog and

Zfp42 in mESCs (Bergmann et al., 2015). We also identified several other Platr genes such as Platr27 and Platr14, whose TSS overlap with LTRs with reduced accessibility, RLTR17 and RLTR16B_MM, respectively (Fig. 5D and Table S5). These two pluripotency-associated transcripts tended to be downregulated in adult SCs but were unchanged at PND8 and PND15 (Fig. 5D and Table S1). The other LTR subtypes with decreased accessibility in adult SCs belong to ERV1, ERVL and MaLR families (Fig. 5A). A few non-LTR TEs also had decreased chromatin accessibility particularly 7 DNA element subtypes, 2 satellite subtypes and 1 LINE subtype, respectively (Table S5).

TE subtypes with increased accessibility mostly included members of ERVK, ERVL and ERV1 families (24/42, 57.1%) (Table S5). Interestingly, many LINE L1 subtypes had increased chromatin accessibility in adult SCs (Table S5). When parsing the data for more accessible loci within L1 subtypes, we found several L1 loci situated less than +/-5kb from the TSS of olfactory (Olfr) genes, particularly Olfr gene clusters on chromosome 2, 7 and 11 (Table S5).

TEs are known to provide tissue-specific substrates for TF binding (Sundaram and Wysocka, 2020). We examined the regulatory potential of differentially accessible LTR subtypes by determining the enrichment of TF motifs in these regions. We grouped LTR subtypes by family (EVK, ERV1, ERVL and ERVL-MaLR families) and observed that among families with less accessible chromatin, ERVKs had the highest number of enriched TF motifs in adult SCs (Table S5). Top hits included TFs with known regulatory functions in cell proliferation and differentiation such as FOXL1 and FOXQ1, stem cell maintenance factors ELF1, EBF1 and THAP11 and TFs important for spermatogenesis PBX3, ZNF143 and NFYA/B (Table S5). ERVLs displayed motif enrichment for very few TFs, among which ETV2, recently reported to be the spermatogonial stem cell factor ZBTB7A and the testis-specific CTCF paralog CTCFL (Table S5) (Green et al., 2018).

More accessible LINE L1 subtypes were highly enriched in TF motifs, particularly in multiple members of ETS, E2F and FOX families (Table S5). The most significant motifs belonged to spermatogonial stem cell maintenance and stem cell potential regulators FOXO1 and ZEB1, as well as TFs recently associated with active enhancers of the stem cell-enriched population of SCs such as ZBTB17 and KLF5 (Table S5) (Cheng et al., 2020). More accessible ERV1s were also enriched in several TF binding sites, including spermatogenesis-related TFs (PBX3, PRDM1, NFYA/B), hypoxia inducible HIF1A and cytokine regulators STAT5A/B, suggesting different metabolic demand between postnatal and adult stage (Table S5).

Discussion

This study examines the transcriptomic landscape and chromatin accessibility of mouse SCs from early postnatal development to adulthood. Using FACS-enriched SCs populations at PND8, PND15 and adulthood, it shows that the transcriptome of SCs is dynamically regulated during development and that many factors including pluripotency-related TFs and signaling molecules are differentially expressed. This correlates with changes in chromatin accessibility although not all DEGs have altered accessibility. TEs such as LTR retrotransposons have modified chromatin accessibility, which in some cases affects pluripotency-associated lncRNAs. TF motif enrichment in differentially accessible TEs suggested various regulatory roles in SCs.

Different transcriptional states characterize the postnatal development of SCs

SCs are known to undergo genome-wide transcriptional changes during development (Grive et al., 2019; Hammoud et al., 2015; Hermann et al., 2018). The present data complement previous findings by showing that postnatal maturation of SCs is accompanied by an overall decrease in the expression of pluripotency markers and differential expression of signaling molecules and regulators of cell cycle and spermatogenesis. The high quality and depth of our RNA-seq datasets allowed us to uncover specific transitional states of SCs postnatal maturation that were not known before. For instance, we observed that while the transition from early to mid-postnatal stages (PND8 to PND15) is accompanied by changes in gene expression, the transition from mid-postnatal to adult stages is more extensive, both in number of regulated genes and magnitude of changes. We also observed a clear trend towards up-regulation of gene expression as postnatal maturation proceeds, that may reflect gradual transcriptional increase for some genes and possibly sharp changes at certain developmental stages for others.

Consistently, many transcriptional regulators and chromatin modifiers are differentially regulated during SCs development. Different families of TFs including factors with a specific role in SCs like RFX and DMRT, and factors with more widespread functions across cell types such as members of the AP-1 family are modulated, suggesting the contribution of multiple classes of TFs to SCs development. Notably, several transcriptional repressors and epigenetic silencers are differentially down-regulated which may explain the marked trend towards transcriptional up-regulation in our datasets (Fig. 1D). Interestingly, we also observed a global down-regulation of histone genes starting at PND15 (Table S1), which may modify the composition and structure of chromatin and be associated with chromatin openness (Prado et al., 2017). Overall, our transcriptional data highlight the yet unappreciated dynamic regulation of TFs and epigenetic/chromatin modifiers in SCs during postnatal development.

Chromatin accessibility in SCs is modified at enhancers during postnatal development

Our results suggest that the transition of SCs from early postnatal to adult life is accompanied by changes in chromatin accessibility. In particular, accessibility is increased in regions enriched in H3K4me1, a histone PTM typically associated with primed enhancers, but also to a lesser extent, in regions enriched in H3K27ac, a mark associated with active regulatory elements (Gasperini et al., 2020; Shlyueva et al., 2014). This suggests that the transition of SCs from postnatal to adult stages involves changes in the regulatory genome, in particular at primed enhancers. Enhancer priming has recently been proposed as an important mechanism for SCs differentiation (McCarrey, 2023). Differences in chromatin accessibility may be partially explained by differential binding of TFs to regulatory elements (Klemm et al., 2019). This may involve at least three non-mutually exclusive events: 1) a change in the abundance of TFs, 2) a change in the epigenetic status of TF binding sites due to different amount and/or activity of epigenetic modifiers, and 3) a change in nucleosome positioning around TF binding sites due to different amount and/or activity of chromatin remodelers. For instance, CpG methylation and nucleosome positioning can influence TF binding at regulatory elements and affect chromatin accessibility (Barisic et al., 2019; Kaluscha et al., 2022). The fact that Dnmt1 is downregulated in adult SCs concomitantly with differential expression of TFs and chromatin modifiers may explain changes in chromatin accessibility at potential enhancers.

Modest correspondence between transcription and chromatin accessibility in SCs

Despite robust changes in gene expression in adult SCs, only a subset of DARs could be assigned to DEGs. This suggests that transcriptional regulation and chromatin accessibility at regulatory elements are not always linked and can be dissociated, as previously reported (Chereji et al., 2019; Kiani et al., 2022). This implies that chromatin accessibility may be modified without any effects on transcription and conversely, transcription can be activated or repressed without requiring any change in chromatin accessibility. In turn, it is possible that a repressor can be exchanged with an activator at a regulatory element but have no detectable consequences. It is also possible that changes in chromatin accessibility have no immediate effects on transcription but reflect an intermediate state that primes a gene or genomic locus for later activation or repression such as observed with developmentally- or experience-dependent primed enhancers (Marco et al., 2020; Rada-Iglesias et al., 2011). Consistently, the majority of DARs with increased accessibility overlap with regions enriched in H3K4me1, a histone PTM characteristic of primed enhancers thought to control future transcriptional responses in different cell types (Rada-Iglesias et al., 2011). These regions may be important for responses to external cues in adult SCs. The apparent lack of correspondence between transcription and accessibility may also be technical and due to the stringency of filters we used to assign genomic loci to DARs. We cannot exclude the possibility that modest changes in accessibility (<log2FC1) correspond to changes in TF occupancy. Chromatin accessibility has been suggested to not be the primary determinant of chromatin-mediated gene regulation (Chereji et al., 2019). Such possibility could be addressed by examining genome-wide maps of TFs, cofactors and RNA-Pol II occupancy that would provide more accurate information about the regulatory genome and its relationship with observed gene expression programs.

Chromatin accessibility at TEs during SCs maturation

Our results reveal differences in chromatin accessibility and TFs motif landscape at different TEs subtypes between PND15 and adult SCs. ERVK and ERV1 subtypes were the most abundant categories of TEs that become less accessible in adult SCs, whilst LINE L1 subtypes gained accessibility. Although the majority of these TEs reside in intergenic and intronic regions, we could detect specific loci belonging to differentially accessible ERVK and LINE L1 subtypes localized nearby TSS of distinct gene families. Notably, ERVKs are known to play an important role in the regulation of mRNA and lncRNAs transcription during spermatogenesis (Davis et al., 2017; Sakashita et al., 2020). The landscape of chromatin accessibility at LTRs loci in SCs is known to be unique compared to other mitotic germ cells and meiotic gametes in testis (Sakashita et al., 2020).

RLTR17, one of the LTR subtypes with decreased chromatin accessibility in adult SCs overlaps with the TSS of several downregulated lncRNAs from the Platr family. Platr genes, including Lncenc1 and Platr14 identified in our study, are LTR-associated lncRNAs important for gene expression programs in ESCs (Bergmann et al., 2015). Interestingly, RLTR17 has also previously been linked to pluripotency maintenance. In mESCs, it is highly expressed and enriched in open chromatin regions and has been shown to provide binding sites for the pluripotency factors Oct4 and Nanog (Fort et al., 2014). Therefore, we suggest that RLTR17 chromatin organization may play a role in regulating pluripotency programs between early postnatal and adult SCs.

In contrast to decreased accessibility at LTRs loci, LINE L1 subtype had increased accessibility in adult SCs. Some of these L1 loci were located close to olfactory receptor genes with upregulated mRNA expression. Recent findings in mouse and human ESCs suggest a non-random genomic localization of L1 elements, specifically at genes encoding proteins with specialized functions (Lu et al., 2020). Among them, the Olfr gene family was the most enriched in L1 elements (Lu et al., 2020). Although their role in SCs is currently not established, Olfr proteins have been implicated in swimming behaviour of sperm cells (Fukuda and Touhara, 2005). Given their dynamic regulation in SCs, we speculate that Olfr genes play additional roles in spermatogenesis, other than in sperm physiology. Together with the high number of enriched TF motifs identified at differentially accessible ERVKs and LINE L1 elements, these data highlight a regulatory role for chromatin organization of TEs in SCs during the transition from development to adult stage (Sundaram et al., 2014; Sundaram and Wysocka, 2020).

Limitations

One limitation of our study is the partial purity of SCs populations obtained by FACS, due to the fact that it does not fully remove other testicular cells from the samples. This therefore does not exclude the possible influence of contaminating cells on the data and their interpretation. Such influence may explain differences between our datasets and datasets in the literature. Nevertheless, by comparing chromatin landscape in developing and adult SCs, the study reveals an age-dependent dynamic reorganization of chromatin accessibility in these cells. When integrated with our transcriptomic datasets and published histone PTMs profiles, the data provide novel insight into transcriptome-chromatin dynamics of mouse SCs from early postnatal life to adulthood. This represents an important resource for future studies on mouse germ cells.

Methods

Mouse husbandry

Male C57Bl/6J mice were purchased from Janvier Laboratories (France) and bred in-house to C57Bl/6J primiparous females to generate males for the experiments. All animals were kept on a reversed 12-h light/12-h dark cycle in a temperature- and humidity-controlled facility with food (M/R Haltung Extrudat, Provimi Kliba SA, Switzerland) and water provided ad libitum. Cages were changed once weekly. Animals from 2 independent cohorts generated separately were used for the experiments.

Germ cells isolation

Germ cells were isolated at postnatal day (PND) 8 or 15 for RNA-seq and ATAC-seq, and at postnatal week 20 (adult) for ATAC-seq. Testicular single-cell suspensions were prepared as previously described with slight modifications (Kubota et al., 2004a, 2004b). For PND8 and PND15 cells, testes from 2 animals were pooled for each sample while for adults, testes from individual males were used. Pup testes were collected in sterile HBSS on ice. Tunica albuginea was gently removed from each testis, making sure to keep seminiferous tubules as intact as possible. Tubules were enzymatically digested in 0.25% trypsin-EDTA (ThermoFisher Scientific) and 7mg/ml DNase I (Sigma-Aldrich) solution for 5 min at 37oC. The suspension was vigorously pipetted up and down 10 times and incubated again for 3 min at 37oC. The digestion was stopped by adding 10% fetal bovine serum (ThermoFisher Scientific) and cells were passed through a 20μm-pore-size cell strainer (Miltenyi Biotec) and pelleted by centrifugation at 600g for 7 min at 4oC. Cells were resuspended in PBS-S (PBS with 1% PBS, 10 mM HEPES, 1 mM pyruvate, 1mg/ml glucose, 50 units/ml penicillin and 50 μg/ml streptomycin) and used for sorting. Adult testes were collected and the tunica was removed. Seminiferous tubules were digested in 1 mg/mL collagenase type IV (Sigma Aldrich) and then in 0.25 % Trypsin-EDTA (Gibco), both times for 8-12 minutes and in the presence of DNase I solution (Sigma Aldrich). FBS (Cytiva HyClone) was added to a concentration of 10 % and the cell suspension was filtered through a 40μm-pore-size cell strainer (Corning) and washed with DPBS-S (1 % FBS, 10 mM HEPES (Gibco), 1mM sodium pyruvate (Gibco), 1mg/mL Glucose (Gibco) and 1X Penicillin-Streptomycin (Gibco) in DPBS (Gibco). 5 mL of the single cell suspension was then overlayed on 2 mL 30 % Percoll (Sigma) and centrifuged with disabled break at 600g for 8 minutes at 4 °C. The supernatant was removed and cell suspension was washed twice in DPBS-S and used for sorting.

SCs enrichment by FACS

For pup testis, dissociated cells were stained with BV421-conjugated anti-β2M, biotin-conjugated anti-THY1 (53-2.1) and PE-conjugated anti-αv-integrin (RMV-7) antibodies. THY1 was detected by staining with Alexa Fluor 488-Sav. For adult testes, cells were stained with anti-α6-integrin (CD49f; GoH3), BV421-conjugated anti-β2 microglobulin (β2M; S19.8) and R-phycoerythrin (PE)-conjugated anti-THY1 (CD90.2; 30H-12) antibodies. α6-Integrin was detected by Alexa Fluor 488-SAv after staining with biotin-conjugated rat anti-mouse IgG1/2a (G28-5) antibody. Prior to FACS, 1 μg/ml propidium iodide (Sigma) was added to cell suspensions to discriminate dead cells. All antibody incubations were performed in PBS-S for at least 30 min at 4oC followed by washing in PBS-S. Antibodies were obtained from BD Biosciences (San Jose, United States) unless otherwise stated. Cell sorting was performed at 4oC on a FACS Aria III 5L using an 85μm nozzle at the Cytometry Facility of University of Zurich. For RNA-seq at PND8 and PND15, cells were collected in 1.5 ml Eppendorf tubes in 500 μL PBS-S, immediately pelleted by centrifugation and snap frozen in liquid N2. Cell pellets were stored at -80oC until RNA extraction. For RNA-seq of adults, 1000 spermatogonial cells (MHC I negative, alpha-6-integrin and Thy1 positive) per male were sorted into PBS. For Omni-ATAC at PND15, 25’000 cells were collected in a separate tube, pelleted by centrifugation and immediately processed using a library preparation protocol (Corces et al., 2017). For Omni-ATAC in adults, 5000 cells from each animal were collected in a separate tube and processed using the same protocol.

Immunocytochemistry

The protocol used for assessing SCs enrichment after sorting was kindly provided by Jon Oatley at Washington State University, Pullman, USA (Yang et al., 2013). Briefly, 30,000-50,000 cells were adhered to poly-L-Lysine coated coverslips (Corning Life Sciences) in 24-well plates for 1 h. Cells were fixed in freshly prepared 4% PFA for 10 min at room temperature then washed in PBS with 0.1% Triton X-100 (PBS-T). Non-specific antibody binding was blocked by incubation with 10% normal goat serum for 1 h at room temperature. Cells were incubated overnight at 4oC with mouse anti-PLZF (0.2 μg/ml, Active Motif) primary antibody. Alexa488 goat anti-mouse IgG (1 µg/mL, ThermoFisher Scientific) was used for secondary labelling at 4oC for 1 h. Coverslips were washed 3 times, mounted onto glass slides with VectaShield mounting medium containing DAPI (Vector Laboratories) and examined by fluorescence microscopy. SCs enrichment was determined by counting PLZF+ cells in 10 random fields of view from each coverslip and dividing by the total number of cells present in the same fields of view (DAPI-stained nuclei). The number of PLZF+ and PLZF-cells from the 10 different fields of view was averaged.

RNA extraction and RNA-seq library preparation

For RNA-seq at PND8 and PND15, total RNA was extracted from sorted cells using AllPrep RNA/DNA Micro kit (Qiagen). RNA quality was assessed using a Bioanalyzer 2100 (Agilent Technologies). Samples were quantified using Qubit RNA HS Assay (ThermoFisher Scientific). RNA sequencing libraries were prepared using SMARTer Stranded Total RNA-Seq Kit v3 (Takara Bio USA, Inc) following the recommended protocol with minor adjustments. For PND8 and PND15 SCs libraries, RNA was fragmented for 4 minutes at 94 °C. After reverse transcription, samples were barcoded using SMARTer Unique Dual Index Kit (Takara Bio USA, Inc). Five PCR cycles were run and PCR products were purified using AMPure XP reagent (Beckman Coulter Life Sciences; bead:sample ratio 0.8X). Ribosomal cDNA was depleted following the protocol and final library amplification was performed with 12 PCR cycles. Libraries were purified twice using AMPure XP reagent (bead:sample ratio 1X) and eluted in Tris buffer. Adult SCs libraries were prepared using the option to start directly from cells (1000 cells as input). Cells were incubated for 6 minutes at 85 °C and processed as indicated for PND8 and PND15 SCs. Samples were pooled at equal molarity, as determined in a 100-800 bp window on the Bioanalyzer.

Omni-ATAC library preparation and sequencing

Chromatin accessibility was profiled from PND15 and adult SCs. Libraries were prepared according to Omni-ATAC protocol, starting from 25,000 PND15 and 5000 adult sorted SCs (Corces et al., 2017). Briefly, sorted cells were lysed in cold lysis buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% NP40, 0.1% Tween-20, and 0.01% digitonin) and nuclei were pelleted and transposed using Nextera Tn5 (Illumina) for 30 min at 37oC in a thermomixer with shaking at 1000 rpm. Transposed fragments were purified using MinElute Reaction Cleanup Kit (Qiagen). Following purification, libraries were generated by PCR amplification using NEBNext High-Fidelity 2X PCR Master Mix (New England Biolabs) and purified using Agencourt AMPure XP magnetic beads (Beckman Coulter) to remove primer dimers (78bp) and fragments of 1000-10,000 bp. Library quality was assessed on an Agilent High Sensitivity DNA chip using the Bioanalyzer 2100 (Agilent Technologies). 6 samples were sequenced from PND15 and 5 from adult SCs.

RNA-seq analysis

Quality control and alignment: 100bp single end sequencing was performed at the Functional Genomics Center Zurich (FGCZ) using a Novaseq SP flowcell on the Novaseq 6000 platform. Quality assessment of FASTQ files was done using FastQC (Andrews et al., 2012) (version 0.11.8). TrimGalore (Krueger, 2015) (version 0.6.2) was used to trim adapters and low-quality ends from reads with Phred score less than 30 (-q 30) and to discard trimmed reads shorter than 30bp (--length 30). Trimmed reads were pseudo-aligned using Salmon (Patro et al., 2017) (version 0.9.1) with automatic detection of the library type (-l A), correcting for sequence-specific bias (--seqBias) and correcting for fragment GC bias correction (--gcBias) on a transcript index prepared for the Mouse genome (GRCm38) from GENCODE (version M18) (Harrow et al., 2012), with additional piRNA precursors and TEs (concatenated by family) from Repeat Masker as in (Gapp et al., 2018).

Downstream analysis

Data analyses and plotting were conducted with R (R Core Team, 2020) (version 3.6.2) using packages from The Comprehensive R Archive Network (CRAN) (https://cran.r-project.org) and Bioconductor (Huber et al., 2015). Pre-filtering of genes was done using the filterByExpr function from edgeR (Robinson et al., 2009) (version 3.28.1) with a design matrix requiring at least 15 counts (min.counts=15). Normalization factors were obtained using TMM normalization (Robinson and Oshlack, 2010) from edgeR package and differential gene expression analyses were done using limma-voom (Law et al., 2014) pipeline from limma (Ritchie et al., 2015) (version 3.42.2). GO enrichment analyses were performed using g:Profiler (https://biit.cs.ut.ee/gprofiler/gost) querying against Mus musculus database.

ATAC-seq analysis

Quality control, alignment, and peak calling: Paired-end (PE) sequencing was performed on PND15 and adult SCs on Illumina HiSeq2500 platform at FGCZ. FASTQ files were assessed for quality using FastQC (Andrews et al., 2012) (version 0.11.8). QC was performed using TrimGalore (Krueger, 2015) (version 0.6.2) in PE mode (--paired), trimming adapters, low-quality ends (-q 30) and discarding reads shorter than 30bp after trimming (--length 30). Alignment on GRCm38 genome was performed using Bowtie2 (Langmead and Salzberg, 2012) (version 2.3.5) with the following parameters: allowing fragments up to 2kb to align (-X 2000), entire read alignment (--end-to-end), suppressing unpaired alignments for paired reads (--no- mixed), suppressing discordant alignments for paired reads (--no-discordant) and minimum acceptable alignment score with respect to read length (--score-min L,- 0.4,-0.4). Using alignmentSieve (version 3.3.1) from deepTools (Ramirez et al., 2016) (version 3.4.3), aligned data (BAM files) were adjusted for read start sites to represent the center of the transposon cutting event (--ATACshift), and filtered for reads with a high mapping quality (--minMappingQuality 30). Reads mapping to the mitochondrial chromosome and ENCODE blacklisted regions were filtered out. To call nucleosome-free regions, all aligned files were merged within groups (PND15 and adult), sorted and indexed using SAMtools (Li et al., 2009) (version 0.1.19) and nucleosome-free fragments (NFFs) were obtained by selecting alignments with a template length between 40 and 140bp in length. Peak calling on NFFs was performed using MACS2 (Zhang et al., 2008) (version 2.2.7.1) with mouse genome size (-g 2744254612) and PE BAM file format (-f BAMPE).

Differential accessibility analysis

These analyses were conducted in R (version 3.6.2) using packages from CRAN (https://cran.r-project.org) and Bioconductor (Huber et al., 2015). Peaks were annotated based on overlap with GENCODE (version M18) (Harrow et al., 2012) transcript and/or distance to the nearest TSS (https://github.com/mansuylab/SC_postnatal_adult/bin/annoPeaks.R). The number of extended reads overlapping in peak regions was calculated using csaw package (Lun and Smyth, 2015) (version 1.20.0). Peak regions that did not have at least 15 reads in at least 40% of the samples were filtered out. Normalization factors were obtained on filtered peak regions using TMM normalization method (Robinson and Oshlack, 2010) and differential analysis on peaks (PND15 versus adult) was performed using Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood (glmQLFit) Tests from edgeR package (Robinson et al., 2009) (version 3.28.1). Peak regions with an absLog2FC≥1 and adjusted-P≤0.05 were categorized as DARs.

Downstream analysis

Heatmaps of normalized ATAC-seq signal were created using deepTools with default parameters. Genomic annotation of ATAC-seq peak and DARs and identification of closest gene DARs were performed using ChIPpeakAnno (Zhu et al., 2010). GO enrichment analysis of genes associated to DARs was performed using g:Profiler (Raudvere et al., 2019) querying against the Mus musculus database. For the epigenomic annotation of DARs, publicly available signal files and peak files for all histone marks were used (Cheng et al., 2020). Heatmaps of ChIP-seq signal for histone PTMs around ATAC-seq peaks and DARs were generated with deepTools. Overlap between ChIP-seq peaks and DARs was generated using Intervene (Khan et al., 2017). TF motif analysis was performed using MEME-ChIP (Ma et al., 2014) and motif identification was done querying identified motifs against JASPAR database (Castro-Mondragon et al., 2021). Quantification and identification of differences in chromatin accessibility between promoter regions of DEGs was performed with deepStats (Richard 2019). To quantify differences in chromatin accessibility between DEGs, bamscale cov (Pongor et al., 2020) was employed using ATAC-seq BAM files as input and genomic coordinates of up-regulated and down-regulated genes. Kruskal-Wallis test was used to test for significant differences in chromatin accessibility between each category of DEGs at postnatal and adult stage.

Differential accessibility analysis at TEs

TE gene transfer format (GTF) file was obtained from http://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/mm10_rmsk_TE.gtf.gz on 03.02.2020. This file provides hierarchical information about TEs: Class (level 1, eg. LTR), family (level 2, eg. LTR - >L1), subtype (level 3, eg. LTR-> L1->L1_Rod), and locus (level 4, eg. LTR->L1 -> L1_Rod -> L1_Rod_dup1). TE loci were annotated based on overlap with GENCODE (version M18) as described above for ATAC-seq peaks. Filtered BAM files (without reads mapping to blacklisted or mitochondrial regions) were used to analyze TEs. Mapped reads were assigned to TEs using featureCounts from R package Rsubread (Liao et al., 2019) (version 2.0.1) and were summarized to Subtypes (level 3), allowing for multi-overlap with fractional counts, while ignoring duplicates. The number of extended reads overlapping at TE loci were obtained using csaw package (Lun and Smyth, 2015) (version 1.20.0). Subtypes without at least 15 reads, and loci without at least 5 reads in at least 40% of samples were filtered out. Normalization and differential accessibility analysis were performed as described above. Subtypes which had an absolute Log2FC≥0.5 and adjusted-P≤0.05 were categorized as differentially accessible subtypes and loci with an absLog2FC≥1 and adjusted-P≤0.05 were categorized as differentially accessible loci. For further downstream data analysis, only differentially accessible loci or differentially accessible subtypes were considered. TF motif enrichment analysis was performed using the marge package (Amezquita, 2018) (version 0.0.4.9999), which is a wrapper around the Homer tool (Heinz et al., 2010) (version 4.11.1).

Authors contribution

ILC and IMM designed and conceived the study. ILC performed all RNA-seq, ICC and ATAC-seq experiments with support of LCS. DKT and RGA-M analyzed the RNA-seq, ATAC-seq, ChIP-seq with significant support from PLG. OUF and MC performed chromatin accessibility analysis of DEGs. ILC, DKT and RGA-M prepared figures. ILC and RGA-M interpreted the data with significant input from DKT, PLG and IMM. ILC, DKT, RGA-M and IMM wrote the manuscript, with significant input from PLG. All authors read and accepted the final version of the manuscript.

Acknowledgements

We thank Martin Roszkowski and Francesca Manuella for taking care of animals breeding, Yvonne Zipfel for animal care, Silvia Schelbert and Alberto Corcoba for taking care of animal licenses and lab organization. We thank Catherine Aquino and Emilio Yángüez from Functional Genomics Center Zurich for support and advice with libraries preparation and sequencing. We are very grateful to Jon Oatley, Melissa Oatley, Tessa Lord and Nathan Law for conceptual advice, hands-on training, and for providing detailed protocols for testis dissection and preparation, and immunocytochemistry of SCs. We thank Zuguang Gu for support with heatmaps. We thank Service and Support for Science IT (www.s3it.uzh.ch) for computational infrastructure.

Competing interest

The authors declare that they have no competing interests.

Funding

The Mansuy lab is funded by the University Zürich, the ETH Zürich, the Swiss National Science Foundation grant number 31003A_175742/1, the National Centre of Competence in Research (NCCR) RNA&Disease funded by the Swiss National Science Foundation (grant number 182880/Phase 2 and 205601/Phase 3), ETH grants (ETH-10 15-2 and ETH-17 13-2), European Union Horizon 2020 Research Innovation Program Grant number 848158, European Union projects FAMILY and HappyMums funded by the Swiss State Secretariat for Education, Research and Innovation (SERI), FreeNovation grant from Novartis Forschungsstiftung, and the Escher Family Fund. RGA-M received an ETH Postdoctoral Fellow grant number 20-1 FEL-28. DKT received a Swiss Government Excellence Scholarship.

Data and material availability

Sequencing data is available via the EBI BioStudies database within the ArrayExpress collection (RNA-seq: E-MTAB-12721; ATAC-seq: E-MTAB-12722). The code employed for data analysis is available from https://github.com/mansuylab/SC_postnatal_adult/. Additional information is available upon request.

Supplementary Figure 1. Expression of SCs markers and their dynamics between early postnatal and adult stage.

(A) Heatmap of the expression profile of an extended list of selected markers of spermatogonial and different testicular somatic cells extracted from total RNA-seq data on PND8, PND15 and adult samples (n=6 for each group). Each row in the heatmap represents a biological replicate from two experimental batches. Gene names in bold correspond to cellular markers presented in Figure 1C. Gene expression is represented as Log2CPM (counts per million).

(B) Left, Heatmap of the expression profile of genes involved in germ-cell maintenance, pluripotency and signalling as reported by Hammoud, 2015. Each row in the heatmap represents a biological replicate from two experimental batches. Unsupervised hierarchical clustering was applied to each of the rows of the heatmaps and a dendrogram indicating similarity of expression profiles among genes in each biological category is shown with each heatmap. Right, line-plots of the average expression of each gene displayed in the heatmaps showing the dynamics and breadth of gene expression across SCs postnatal development.

(C) Genomic snapshots from IVG of exemplary DEGs showing aggregated RNA-seq signal from PND8, PND15 and adult stages.

(D) Left, bar-plot of GO KEGG pathways categories enriched in DEGs between PND15 and adult SCs (adjusted-P≤0.05). Dotted line indicates a threshold value for significance of 0.05. Right, Heatmap of DEGs between PND8 and PND15 belonging to the GO category “ECM-receptor interaction”. Each row in the heatmap represents a biological replicate from two experimental batches. Shown are Log2FC with respect to the average of PND8.

Supplementary Figure 2. Specific transcriptional programs between postnatal and adult SCs.

(A) Four-way Venn diagram of DEGs detected between PND8-PND15 and PND15-adult comparisons.

(B) Heatmap of all DEGs specific to PND8 stage. Each row in the heatmap represents a biological replicate from two experimental batches. Shown are the Log2 FC with respect to the average of PND8.

(C) Heatmap of all DEGs with significant changes across the 3 developmental time points. Each row in the heatmap represents a biological replicate from two experimental batches. Shown are the Log2 FC with respect to the average of PND8.

Supplementary Figure 3. Chromatin accessibility landscape of SCs from PND15 to adult stage.

(A) Heatmaps showing normalized ATAC-seq signal for all identified ATAC-seq peaks from PND15 and adult stages. Each row in the heatmap represents a 2kb genomic region extended 1kb down- and upstream from the center of each identified peak. Each row is ordered in a decreasing level of average accessibility. On top of each heatmap is shown a plot of the signal profile over all ATAC-seq peaks and their extended genomic region.

(B) Bar plot of the genomic distribution of all identified ATAC-seq peaks from PND15 and adult stages.

(C) IGV tracks for ATAC-seq signal for PND15 and adult SCs showing ATAC-seq peaks located at promoters and intragenic regions of genes with important functions in SCs. Highlighted genomic segments correspond to ATAC-seq peaks.

(D) Heatmaps showing normalized ChIP-seq signal for different histone marks in adult Scs (public data derived from Cheng et al., 2020) at DARs from PND15 and adult stages. Each row in the heatmap represents a 2kb genomic region extended 1kb down- and upstream from the centre of each identified peak. Shown data corresponds to non-regenerative SCs as stated in Cheng et al., 2020.