Dynamics of transcriptional programs and chromatin accessibility in mouse spermatogonial cells from early postnatal to adult life

In mammals, spermatogonial cells (SPGs) are undifferentiated male germ cells in testis that are quiescent until birth and then self-renew and differentiate to produce spermatogenic cells and functional sperm from early postnatal life throughout adulthood. The transcriptome of SPGs is highly dynamic and timely regulated during postnatal development. We examined if such dynamics involves changes in chromatin organization by profiling the transcriptome and chromatin accessibility of SPGs from early postnatal stages to adulthood in mice using deep RNA-seq, ATAC-seq and computational deconvolution analyses. By integrating transcriptomic and epigenomic features, we show that SPGs undergo massive chromatin remodeling during postnatal development that partially correlates with distinct gene expression profiles and transcription factors (TF) motif enrichment. We identify genomic regions with significantly different chromatin accessibility in adult SPGs that are marked by histone modifications associated with enhancers and promoters. Some of the regions with increased accessibility correspond to transposable element subtypes enriched in multiple TFs motifs and close to differentially expressed genes. Our results underscore the dynamics of chromatin organization in developing germ cells and complement existing datasets on SPGs by providing maps of the regulatory genome at high resolution from the same cell populations at early postnatal, late postnatal and adult stages collected from single individuals.


Introduction
Spermatogonial cells are the initiators and supporting foundation of mammalian spermatogenesis. In mice, they become active at postnatal day (PND) 1 to 2 when they exit mitotic arrest and start cycling, in order to populate the basement membrane of the seminiferous tubules. During the first week after birth, a population of spermatogonial cells continue to proliferate and give rise to the undifferentiated Asingle (As), Apaired (Apr) and Aaligned (Aal) spermatogonial pool. The remaining spermatogonial cells proceed into differentiation and form chains of daughter cells that become primary and secondary spermatocytes around PND 10 to 12. A second round of meiosis gives rise to haploid spermatids, which then mature to form adult spermatozoa capable of fertilization by PND 35 to 37 Brinster, 2006, 2012;De Rooij, 2017). Recent work has shown distinct epigenetic and transcriptional profiles in spermatogonial cell populations in postnatal compared to adult testis (Green et al., 2018;Hammoud et al., 2014Hammoud et al., , 2015Hermann et al., 2018;Law et al., 2019). The pool of spermatogonial cells in early postnatal testis displays unique characteristics necessary for rapid establishment and expansion. In contrast, in adulthood the focus lies on the maintenance of a steady spermatogonial population which balances proliferation and differentiation capabilities to ensure adult sperm formation across life (Law and Oatley, 2020  OmniATAC-seq  revealed an extensive change, in particular increase in accessibility, in adult spermatogonia compared to PND15 cells. Furthermore, motif enrichment analysis identified candidate transcription factors (TFs) specifically enriched in the regions of differential accessibility, suggesting distinct regulatory mechanisms contributing to the differential chromatin organization between pup and adult spermatogonial cells. By combining our datasets with previously published chromatin immunoprecipitation sequencing (ChIP-seq) and whole genome bisulfite sequencing (WGBS) data from postnatal and adult spermatogonia (Hammoud et al., 2014(Hammoud et al., , 2015, we were able to identify chromatin regions for which the change in accessibility was marked by distinct histone modifications and DNA methylation profiles, providing novel insight into the epigenetic programming associated with specific classes of differentially expressed genes. Lastly, by performing chromatin accessibility analysis specifically on repeat elements (REs), we revealed changes in accessibility at all main families of repeats in adult spermatogonial cells, including in promoters of genes with differential expression, suggesting a potential role for REs in regulating gene expression between postnatal and adult stages.

Fluorescence activated cell sorting highly enriches for spermatogonial cells from postnatal and adult mouse testis
To enrich for spermatogonial cells from postnatal and adult mice testis, we performed surface marker-based fluorescence activated-cell sorting (FACS) following the protocol previously established by Kubota et al (Kubota et al., 2003, 2004a, 2004b. To validate the enrichment of spermatogonial cells in our sorted populations, we performed immunocytochemistry using PLZF, a well-established marker of undifferentiated spermatogonia (Costoya et al., 2004). Following staining of unsorted and sorted cell fractions, we identified between 85-99% PLZF + cells in the sorted fractions compared to 3-6% PLZF + cells in the unsorted samples (Fig. S1A). To further confirm the high enrichment of spermatogonial cells in our sorted samples, we also compared expression of known pluripotency and spermatogonial stem cell markers between our sorted fractions and previously published datasets from Oct4-GFP and Id4-GFP transgenic mice (Garcia and Hofmann, 2012;Helsel et al., 2017) and found similar expression levels for key stem cell and undifferentiated spermatogonial markers (Fig. S1B). Noteworthy, our sorted populations still showed low expression of genes expressed in differentiating spermatogonial cells and scarce expression of somatic contaminants. Additionally, to exclude biological variability from our inferential analysis and detect statistically significant differences, we have included between 5 and 9 biological replicates for each individual timepoint in both RNA-seq and ATAC-seq analyses.

Dynamic transcriptome states characterize spermatogonial cells between early postnatal development and adult stage
Extensive work in rodents has underlined the unique characteristics of the first wave of spermatogenesis, including distinct transcriptome features of early postnatal spermatogonial cells (Law et al., 2019;Yoshida et al., 2006). Moreover, discernable differences between gene expression profiles of spermatogonia during the first two weeks of postnatal development have also been documented (Hammoud et al., 2015). To capture the early postnatal transcriptome transitions, and to obtain a comprehensive view of the gene expression programs in the developing testis, we performed RNA-seq on FACS sorted spermatogonial cells from 2 distinct postnatal stages (PND8 and PND15) and from adult 20 weeks old (PNW20) mouse testis. Pairwise comparisons revealed a high correlation between biological replicates at each of the three timepoints (Fig. S2), while principal component analysis (PCA) pointed to an age-dependent grouping of the samples, with postnatal spermatogonial cells more closely related and distanced from adult spermatogonia (Fig. 1A).
Initially, we performed gene set enrichment analysis (GSEA) to identify distinct classes of over / under represented GO terms between PND8 and PND15 spermatogonial cells. Our analysis revealed a downregulation of genes involved in cell cycle and cell proliferation pathways, consistent with the increased proliferative and colonization capacities of spermatogonial cells specific to the first week of postnatal development, and essential for successful colonization of the testis basement membrane (Hammoud et al., 2015). Among other candidates, the most downregulated GO terms were RNA splicing, DNA replication, and ribosome biogenesis (Fig. S3A). GSEA between PND15 and adult spermatogonial cells revealed an upregulation of genes related to mitochondrial function, mRNA metabolism and RNA splicing (Fig. S3B). To better understand the relationship between the two transitions, from PND8 to PND15 and from PND15 to adult stage, we continued by performing unsupervised k-means clustering of all differentially expressed genes across the three stages, which led to five clusters with distinct gene expression profiles (Fig. 1B). Overrepresentation analysis (ORA) revealed distinct sets of Gene Ontology (GO) terms enriched within each gene cluster (Fig. 1C). Genes with a marked downregulation in adult spermatogonia compared to postnatal stages were enriched for DNA packaging, protein-DNA complex assembly and nucleosome assembly (Cluster 1). A similar expression trajectory was detected for genes involved in GTPase-mediated signal transduction and activity, and histone and chromatin modifications (Cluster 2). Genes involved in synapse organization, synaptic vesicle exocytosis, axon guidance and modulation of chemical synaptic transmission were maximally expressed at PND15, followed by a downregulation in adult spermatogonia (Cluster 3). In contrast, mitochondrial function, genes involved in ribosome biogenesis and ncRNA processing (Cluster 4) displayed an upregulation, reaching peak expression in adult spermatogonial cells. Cluster 5, which showed the most elevated expression levels in adults, consisted of genes involved in the electron transport chain, spermatogenesis and ribonucleoprotein complex (Fig. 1C). Collectively, our RNA-seq data reveal the main transcriptional changes accompanying the transition from postnatal to adult spermatogonia, and support similar findings from recent single cell profiling studies of rodent and human spermatogonia (Green et al., 2018;Guo et al., 2017;Hermann et al., 2018;Law et al., 2019;Wang et al., 2018)

Spermatogonia in postnatal and adult testes have different expression profiles of signaling factors important for self-renewal and differentiation
Next, we focused on the expression of factors important for spermatogonial cell fate.
Several key transcription factors involved in stem cell maintenance and proliferation exhibited highly dynamic gene expression profiles between postnatal and adult stages. The majority of genes known to be involved in spermatogonial stem cell maintenance and selfrenewal (Lhx1,Etv5,Id4,Cxcr4,Bcl6b,Pou5f1,Pax7) showed a gradual downregulation from PND8 to PND15 to adult stage (Fig. 1D). Signaling factors important for spermatogonial cell maintenance in vitro, such as LIF receptors (Lifr), showed a transient and moderate upregulation between PND8 and 15, followed by downregulation in adulthood, whereas FGF receptor genes (Fgfr1 and Fgfr2) were maximally expressed at PND8, followed by a gradual decrease at PND15, and an even more pronounced decrease in adult spermatogonial cells (Fig. 1D). Self-renewal factors belonging to the TGFb family, Tgfb1, Tgfb2 and Tgfb3 were also significantly reduced in adult spermatogonial cells (Fig1. D).
Altogether these findings indicate a switch in the stem cell maintenance and renewal strategies employed by spermatogonial cells between postnatal and adult stages, reinforcing previous findings suggesting a change from autocrine to mainly paracrine (relying on the niche) signaling in adult spermatogonial cells (Hammoud et al., 2015). One other class of genes with important roles in spermatogonial cell maintenance and proliferation is the testisspecific melanoma-associated antigen (MAGE) family, specifically the Mage-a and Mage-b subclasses, which were recently found to protect the adult germline against environmental stressors (Fon Tacer et al., 2019;Lee et al., 2020). Notably, our RNA-seq data revealed high levels of MAGE genes in adult spermatogonial cells, with all MAGE genes going through a significant upregulation in adult spermatogonial cells compared to postnatal stages (Fig. 1D).
These findings may suggest that gene programs meant to ensure steady-state sperm production and protect against potentially damaging environmental fluctuations are preferentially upregulated and utilized in adult spermatogonial cells, compared to early postnatal development, when the emphasis falls on rapid and successful expansion and proliferation.

Chromatin accessibility transitions to a more open state from postnatal to adult spermatogonia
To investigate if the transcriptome changes we found between postnatal and adult spermatogonial cells are accompanied by changes in chromatin organization, we profiled chromatin accessibility of 6 biological replicates from PND15 and 5 biological replicates from adult spermatogonial cells using Omni-ATAC-seq. Omni-ATAC-seq assesses chromatin accessibility with a higher signal-to-noise ratio compared to the original ATACseq protocol, starting from as low as few hundred cells . Accessible regions were identified by peak-calling on the merged nucleosome-free fragments (NFF).
Following the removal of lowly enriched regions, we included 158,978 peaks in downstream analyses (see Methods section for details). Most of the Tn5-accesible regions were intergenic (38%), in gene bodies (33%) and in proximity of gene transcription starting sites (+/-1 kbp from TSS, 28%) (Fig. S4). Differential accessibility analysis revealed 3212 differentially accessible regions between PND15 and adult spermatogonia, with the majority of regions showing a gain in accessibility at adult stage ( Fig. 2A). Accessibility changes were localized predominantly to intergenic regions and gene bodies, whilst 15% of all differentially accessible regions resided +/-1kbp from the TSS of a gene (Fig. 2B).
Next, to investigate the functional significance of the regions with differential accessibility between PND15 and adult spermatogonia, we performed GO analysis using the Genomic Regions Enrichment of Annotations Tool (GREAT) (McLean et al., 2010). Overall, regions of increased accessibility in adults were associated with cell fate and stem cell population maintenance, protein metabolism and RNA metabolic processes (Fig. 2C).
Interestingly, when performing GREAT analysis at specific genomic locations where regions with increased accessibility where situated (intergenic, gene bodies, +/-1kbp from TSS) we found that regions located in gene bodies (exons and introns) were enriched for GO terms related mainly to spermatogenesis, protein metabolism and stem cell population maintenance, whilst regions close to, or overlapping with TSSs, were strongly enriched for developmental processes and tissue morphogenesis (Fig. 2D). Regions with decreased accessibility in adult spermatogonia featured an overabundance of terms related to embryonic development ( Fig.   2C), with the majority of these regions located in intergenic spaces (Table S1). Collectively, our differential accessibility analysis and functional analysis of differentially accessible regions suggest that chromatin undergoes extensive remodeling from postnatal to adult stage, and that increased and decreased accessibility regions associate with distinct gene functions, depending on their genomic localization.

Transcription factor motif enrichment reveals potential regulators of chromatin accessibility from postnatal to adult spermatogonial cells
Previous work has underlined the importance of transcription factors (TFs) in establishing and maintaining distinct transcriptional signatures across the developmental trajectory of a cell population or subpopulation (Fushan et al., 2015;Shavlakadze et al., 2019). To predict which TF interactions may mediate the chromatin accessibility remodeling we detected between postnatal and adult spermatogonia, we performed motif enrichment analysis using the Hypergeometric Optimization of Motif EnRichment (HOMER) tool (Heinz et al., 2010a). In regions that showed an increased chromatin accessibility at adult stage, we identified 41 enriched TF motifs (q-value < 0.05) (Fig. 3A). The top candidate list (q-value<0.0001) was dominated by members of the Fos/ Jun family (FOS, FOSB, FOSL1 and FOSL2, JUN, JUNB and JUND) (Fig. 3B). Notably, almost all TFs with enriched motifs displayed stage-specific differences in mRNA levels as well: Fos, Jun and Junb were upregulated in adult spermatogonial cells, while the other FOS/JUN family TFs displayed a decrease in mRNA levels (Fig. 3A). JUN, FOS and CREB are all part of the AP-1 (activating protein-1) superfamily, and play an important role in regulating cell proliferation and death, by mediating the senescence-associated chromatin and transcriptional landscape (Martínez-Zamudio et al., 2020;Shaulian and Karin, 2002). Specific members of this complex, such as JUND and c-FOS have also been studied in the context of spermatogonial stem cell proliferation (He et al., 2008;Wang et al., 2018). Noteworthy, we identified binding sites for JUND and c-FOS in proximity of genes essential for stem-cell self-renewal (Gata2 and Psmd2), spermatogonial cell differentiation (Kit), chromatin dynamics (Chd2) and growth factor adaptor proteins (Frs2) ( Fig. 3E and Table S2). These genes also displayed differential expression between postnatal and adult stage, with Psmd2 showing an increased expression, whilst Gata2, Kit, Chd2 and Frs2 had a decreased expression in adult spermatogonial cells ( Fig. 3D). TF motif analysis using HOMER also revealed enriched binding sites for retinoic acid receptors such as RXRA and RAR ( Fig 3A). Notably, among the regions of increased accessibility, we observed motifs for RXRA and RAR situated in proximity of chromatin remodeling factors Satb1, Hdac4 and Eed, with all 3 genes displaying differential expression between postnatal and adult spermatogonial cells. (Fig 3D-E). These results suggest a role for nuclear receptors in the stage-specific chromatin remodeling that we found in adult spermatogonial cells, potentially via transcriptional regulation of HDAC4 and Polycomb Repressive Complex 2 (PRC2) core component EED. We also performed motif enrichment analysis for the more accessible chromatin regions situated in gene bodies, intergenic regions and +/-1kbp from TSS, and identified several TF motifs specifically enriched in intergenic regions and not present in gene bodies such as the ubiquitously expressed NF-Y complex, NF-YA, NF-YB and NF-YC (Fig. 3C). Interestingly, in distal regions, NF-Y has been shown to facilitate a permissive chromatin conformation, and play an important role in the expression of core embryonic stem cell (ESC) pluripotency genes (Oldfield et al., 2014). Its enrichment in intergenic regions of increased chromatin accessibility could suggest a similar role for NF-Y trimers in adult spermatogonial cells. Furthermore, NF-YA/B motif enrichment has also been found in regions of open chromatin in human spermatogonial cells (Guo et al., 2017), suggesting a potentially similar role for both rodent and human spermatogonia.
Although regions of increased chromatin accessibility predominate in adult spermatogonia compared to PND15, we have also identified TFs with enriched motifs in the regions of decreased accessibility (Fig. 3A). Interestingly, almost all of these TFs were not found in the regions of increased chromatin accessibility. Top hits included members of the FOX family (FOXO1, FOXO3, FOXP2, FOXK1, FOXA2) and members of the ETS and ETS-related families (ETS1, GABPA, ETV4, ELF1, ELF3) (Fig. 3B). The gene expression levels of most of these TFs were also decreased in adult spermatogonial cells (Fig. 3A).
FOXO1 is a pivotal regulator in the self-renewal and differentiation of spermatogonial stem cells in both pup and adult testis, acting via the PI3K-Akt signaling pathway (Chan et al., 2014;Goertz et al., 2011). The roles of the various ETS-related TFs in spermatogonial cells have not been clarified, however recently published data found ETV4 in the stem-cell

Chromatin accessibility differences between postnatal and adult spermatogonia associate with distinct histone modifications and gene programs
Next, we investigated how the differentially accessible chromatin regions associate with the transcriptome changes that we have identified from our RNA-seq data. To aid the functional interpretation of our findings, we divided differentially accessible regions in distal -situated more than 2.5 kb from the TSS of a gene, and proximal -situated less than 2.5 kb from the TSS of a gene. Additionally, we made use of published ChIP-seq and WGBS data from Hammoud et al (Hammoud et al., 2014(Hammoud et al., , 2015 and investigated histone mark  Table S3). Genes important in cell cycle regulation and stem cell maintenance were also present in Category 1, including DAZL-mediated Psmd1 and Psmd2 (Mikedis et al., 2020), and GDNF signaling receptor Gfra2 (Hammoud et al., 2015) (Table S3). Category 2 included proximal regions with increased chromatin accessibility and decreased expression of nearby genes in adult spermatogonia, indicative of active repression taking place (Fig. 4A). Notably, regions in Category 2 were enriched in binding motifs for FOS and JUND, suggesting a role for these TFs in mediating the transcriptional repression of the nearby genes (Fig. S5D). Additionally, a subset of Category 2 regions, which mainly associated with developmental genes, were marked by H3K27me3 together with H3K4me3, suggesting a poised state of these genes, which ensures silencing in the germline, but a competent transcriptional state in future embryos (Hammoud et al., 2014) (Fig. 4B). Examples include developmental factors Satb1, Hmx1, Hoxb5 and Gata2 ( Fig. 4C and Table S3). GO analysis using GREAT confirmed that regions in Category 2 were mainly associated with genes involved in regulation of developmental processes (Fig. S5C). Not surprisingly, DNAme levels in both Category 1 and 2 regions were inversely correlated with the enrichment for histone marks. Regions in Category 3 displayed decreased chromatin accessibility and a downregulation of nearby genes in adult spermatogonia and were relatively few in number compared to Categories 1 and 2 (Fig. 4A). Notably, most of the regions in this category did not show ChIP-signal for any of the 3 histone marks investigated. Category 4 comprised an even smaller number of proximal regions with decreased accessibility (<20 regions) at genes which were upregulated in adult SCs (Fig. 4A). Categories 5 and 6 consisted of regions with an increased and decreased chromatin accessibility respectively, for which nearby gene expression was not detected in PND15 and adult spermatogonial cells (Fig. S5A). Notably, DNAme profiles from postnatal to adult stages did not show drastic changes across any of the 4 categories of proximal regions, suggesting a relatively stable DNAme profile in the transition from postnatal to adult stage. Aside from proximal regions, we also identified regions of differential chromatin accessibility between PND15 and adult spermatogonia located distal from TSSs of genes (Fig. 5A). The majority of these regions did not show enrichment for any of the evaluated histone modifications (Fig. 5B). Nevertheless, in distal regions of increased accessibility in adult spermatogonia for which ChIP-seq signal was present, most regions were enriched in H3K4m3, H3K27ac and coupled H3K4me3/H3K27ac, indicative of a potential distal enahncer role of these regions (Fig. 5B). Furthermore, TF motif analysis revealed enrichment of binding sites for nuclear receptors RXRA, RFX3, NR4A1 and RFX2 in these regions (Fig. 5C). In distal regions which were less accessible in adult spermatogonia, a subset of regions displayed H3K27me3 alone or together with H3K4me3 ( Fig. 5B). Similar to proximal regions, DNAme levels were decreased in the presence of histone marks, and did not display major changes in the transition from postnatal to adult stage ( Fig. 5B). Altogether, our data integration reveals distinct epigenetic landscapes associated with differentially expressed genes in postnatal and adult spermatogonial cells, which may be required to coordinate specific signaling pathways important for spermatogonial cell maintenance and spermatogenic potential.

Chromatin accessibility dynamics at repeat elements in adult spermatogonial cells compared to postnatal development
Repeat elements (REs) are under tight control in the germline, achieved through coordinated epigenetic mechanisms involving DNA methylation, chromatin silencing and the PIWI proteins-piRNA pathway (Deniz et al., 2019;Thompson et al., 2016a). To explore potential differences in REs regulation between postnatal and steady-state adult spermatogonial cells, we performed differential accessibility analysis of REs using our ATAC-seq data. Interestingly, we found that the transition from postnatal to adult stage spermatogonia is accompanied by chromatin accessibility changes at members of all major RE families including LTRs, LINEs and SINEs ( Fig. 6A and Table S4). Most of the REs harboring changes in chromatin accessibility were situated in intergenic and intronic regions, and around 10% were located in proximity of a gene (+/-1kbp from a TSS) (Fig. 6B). Next, we investigated if the distribution of REs harboring changes in chromatin accessibility may relate to distinct gene functions. GO enrichment using GREAT revealed that differentially accessible LINE elements were mainly associated with genes involved in embryogenesis and nervous system development and synaptic function, differentially accessible LTRs associated with genes important for immune responses, while SINEs were strongly linked to genes involved in chromatin organization and transcriptional regulation (Fig. 6C). Moreover, for differentially accessible SINEs inserted in proximity of a gene, GREAT analysis identified a strong enrichment for RNA-related functions including transport, processing and rRNA metabolism (Fig. S6A). Similar gene functions were recently found for promoter-situated SINEs in mESCs, suggesting a conserved regulatory role for this repeat family in different types of stem cells (Lu et al., 2020). LTRs integrated in promoters may also contribute to the transcriptional regulation of nearby genes (Lowe et al., 2007;Thompson et al., 2016b). 25% of the differentially accessible LTR elements in our data resided in proximity of genes (Fig.   S6B). For example, the promoter region of the lncRNA Lncenc1, an important regulator of pluripotency in mESCs (Fort et al., 2014;Sun et al., 2018), harbored two LTR repeats with decreased accessibility in adult spermatogonial cells, which correlated with a decreased expression of Lncenc1 in adult spermatogonia (Fig. 6D). Other examples of promoterintegrated LTRs which displayed changes in accessibility in adult spermatogonia included miR-21, an important regulator of stem-cell renewal capacity of spermatogonial cells (Niu et al., 2011), and miR-34, and miR-146, small RNAs preferentially expressed in spermatogonial stem cell-enriched populations compared to somatic cells (Niu et al., 2011) (Table S4).
To further understand the regulatory potential of differentially accessible REs, we performed TF motif enrichment analysis on LINEs, SINEs and LTRs which harbored changes in chromatin accessibility, and localized in proximity of a gene (+/-1kbp from TSS).
Differentially accessible LINE elements were highly enriched in binding sites for transcription initiating factors KLF3/4/5, E2F family, spermatogonial stem cell maintenance factor SALL4 and FOX family members FOXJ2 and FOXP2 (Fig S7). LTRs displayed high enrichment for several TF motifs mainly important for embryonic development such as HXB8, HXA9, ALX1 and ZN143 (Fig S7). Differentially accessible SINEs comprised the highest number of enriched TF motifs, showing specific enrichment for DMRT1/B binding sites ( Fig S7). By analyzing chromatin accessibility and TF motif enrichment specifically at REs, we were able to reveal extensive chromatin transitions taking place at numerous SINEs, LINEs and LTRs between postnatal and adult spermatogonial cells. Furthermore, the nonrandom distribution in the vicinity of genes with distinct functions, and their enrichment in TF binding sites suggests an important contribution of REs to shaping gene regulatory networks in spermatogonial cells during postnatal life.

Discussion
Mammalian spermatogenesis is critical for ensuring reproductive potential across life. Lastly, by investigating chromatin accessibility specifically at REs, we were able to show that REs undergo changes in chromatin accessibility between postnatal and adult stages. We initially revealed a non-random distribution of the differentially accessible repeats, with LINEs and LTRs more likely to associate with genes important for specialized functions such as nervous system development and immune processes respectively, whereas SINEs were mainly associated with genes important for general cellular functions and were particularly enriched in promoters of genes responsible for RNA processing and metabolism.
Furthermore, we identified promoters of genes important for spermatogonial stem cell renewal enriched in differentially accessible LTRs, supporting a novel role for repeat elements in stage-specific regulation of genes with distinct functions. Similar roles for SINEs and L1 elements have been found in gene and chromatin regulation in mESCs (Lu et al., 2020;Thompson et al., 2016c). Additionally, the abundant enrichment in TF motifs we identified in the differentially accessible repeats points towards an important contribution of REs to the regulation of gene expression dynamics between postnatal and adult spermatogonial cells (Sundaram and Wysocka, 2020;Sundaram et al., 2014).
To our knowledge, this is the first time chromatin accessibility has been profiled in adult mouse spermatogonial cells. However, we cannot entirely exclude the influence of rare contaminating cells on the transcriptome and chromatin accessibility data interpretation.
Therefore, we acknowledge that our findings are representative for the broader spermatogonial cell population in the testis, and that further refinement of our investigation will clarify how these stage-specific mechanisms are fine-tuned in certain spermatogonial subpopulations. Collectively, our cross-functional analyses reveal novel informational content of the spermatogonial chromatin states at different postnatal stages, integrate it with specific transcriptome changes, and provide a vast resource for further analysis of mouse postnatal spermatogonial cell programming.

Mouse husbandry
Male C57Bl/6J mice were purchased from Janvier Laboratories (France) and bred inhouse to generate male mice used for 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 breedings were used for the experiments.

Germ cells isolation
Germ cells were isolated from male mice at postnatal day (PND) 8 or 15, and 20 weeks of age (PNW20). Testicular single-cell suspensions were prepared as previously described with slight modifications (Kubota et al., 2003(Kubota et al., , 2004a(Kubota et al., , 2004b. For preparations using PND8 and PND15 pups, testes from 2 animals were pooled for each sample. Pup testes were collected in sterile HBSS on ice. Tunica albuginea was gently removed from each testis, making sure to keep the 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 37 o C. The suspension was vigorously pipetted up and down 10 times and incubated again for 3 min at 37 o C. The digestion was stopped by adding 10% fetal bovine serum (ThermoFisher Scientific) and the cells were passed through a 20µm-pore-size cell strainer (Miltenyi Biotec) and pelleted by centrifugation at 600g for 7 min at 4 o C. 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. For preparations from adults, one adult male was used for each sample. The tunica was removed from testes, and seminiferous tubules were digested in 2 steps. The first consisted in an incubation in 1mg/ml collagenase type IV (Sigma-Aldrich) for 5 min at 37 o C and vigorous swirling until the tubules were completely separated. Then tubules were placed on ice for 5 min to sediment, the supernatant removed and washed with HBSS.
Washing/sedimentation steps were repeated 3 times and were necessary to remove interstitial cells. After the last washing step, sedimented tubule fragments were digested again with 0.25% trypsin-EDTA and 7mg/ml DNase I solution, and the digestion was stopped by adding 10% FBS. The resulting single-cell suspension was filtered through a 40µm strainer (Corning Life Sciences) and washed with HBSS. After centrifugation at 600g for 7 min at 4 o C, the cells were resuspended in PBS-S, layered on a 30% Percoll solution (Sigma-Aldrich) and centrifuged at 600g for 8 min at 4 o C without braking. The top 2 layers (HBSS and Percoll) were removed and the cell pellets resuspended in PBS-S and used for sorting.

Spermatogonial cell enrichment by FACS
For pup testis, dissociated cells were stained with BV421-conjugated anti-β2M, biotin-conjugated anti-Thy-1 (53-2.1), and PE-conjugated anti-αv-integrin (RMV-7) antibodies. Thy-1 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-Thy-1 (CD90.2; 30H-12) antibodies. α6-Integrin was detected by Alexa Fluor 488-SAv after staining with biotinconjugated rat anti-mouse IgG1/2a (G28-5) antibody. Prior to FACS, 1 µg/ml propidium iodide (Sigma) was added to the cell suspensions to discriminate dead cells. All antibody incubations were performed in PBS-S for at least 30 min at 4 o C followed by washing in PBS-S excess. Antibodies were obtained from BD Biosciences (San Jose, United States) unless otherwise stated. Cell sorting was performed on a FACS Aria III 5L, at 4 o C and using a 85µm nozzle, at the Cytometry Facility of University of Zurich. For RNA-seq on PND8 and PND15 spermatogonia, cells were collected in 1.5 ml Eppendorf tubes in 500 µL PBS-S, immediately pelleted by centrifugation and snap frozen in liquid N2. Cells pellets were stored at -80 o C until DNA/RNA extraction. For RNA-seq on adult spermatogonia, 1000 cells from each animal were collected for processing using the SmartSeq2 protocol (Picelli et al., 2014). For OmniATAC-seq on PND15 spermatogonia, 25'000 cells were collected in a separate tube, pelleted by centrifugation and immediately processed using OmniATAC-seq library preparation protocol . For OmniATAC-seq on adult spermatogonia, 5000 cells from each animal were collected in a separate tube and further processed.

Immunocytochemistry
The protocol used for assessing spermatogonial cell enrichment after sorting was kindly provided by the Oatley Lab 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 4 o C with mouse anti-PLZF (0.2 µg/ml, Active Motif, clone 2A9) primary antibody. Alexa488 goat anti-mouse IgG (1 µg/mL, ThermoFisher Scientific) was used for secondary labelling at 4 o C for 1 h. Coverslips were washed 3x and mounted onto glass slides with VectaShield mounting medium containing DAPI (Vector Laboratories) and examined by fluorescence microscopy. Stem cell 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 field of view (DAPI-stained nuclei).

RNA extraction and RNA-seq library preparation
For RNA-seq at PND8 and PND15, total RNA was extracted from sorted cells using the 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). 10 ng of total RNA from each sample were used to prepare long RNA sequencing libraries using the SMARTer® Stranded Total RNA-Seq Kit v2 -Pico Input Mammalian (Takara Bio USA, Inc.) according to the manufacturer's instructions. For adult spermatogonia, 5000 cells from each sample were lysed in 10 uL buffer RLT Plus (Qiagen) immediately after sorting and processed following the SmartSeq2 protocol (Picelli et al., 2014). The number of biological replicates sequenced at each timepoint was 8 for PND8, 9 for PND15 and 6 for adult spermatogonia.

OmniATAC-seq library preparation and sequencing
Chromatin accessibility was profiled from PND15 and adult spermatogonial cells.
The libraries were prepared according to the OmniATAC-seq protocol, starting from 25 000 PND15 and 5000 adult sorted spermatogonia, respectively . 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 the nuclei were pelleted and transposed using the Nextera Tn5 (Illumina) for 30 min at 37 o C in a thermomixer with shaking at 1000 rpm. The transposed fragments were purified using the MinElute Reaction Cleanup Kit (Qiagen). Following purification, libraries were generated by PCR amplification for 12 cycles using the NEBNext High-Fidelity 2X PCR Master Mix (New England Biolabs) and purified using Agencourt AMPure XP magnetic beads (Beckman Coulter) in order to remove primer dimers (78bp) and large fragments of 1000-10,000 bp in length. Library quality was assessed on an Agilent High Sensitivity DNA chip using the Bioanalyzer 2100 (Agilent Technologies). 6 biological replicates were sequenced from PND15 and 5 biological replicates from adult spermatogonial cells.

High-throughput sequencing data analysis
Data availability: The datasets used in this study are available from the following GEO accessions: GSE_____, GSE49621, and GSE49623. An overview of the datasets included in the study is shown in the following table:

RNA-seq
Quality control and alignment: Single-end (SE) sequencing was performed on the PND8, PND15, and adult samples on the Illumina HiSeq4000, Illumina HiSeq4000, and NovaSeq platform, respectively. PND8 raw data (FASTQ files) was merged from two runs for the RNA-Seq samples. The FASTQ files were assessed for the quality using FastQC (Andrews et al., 2012) (version 0.11.8). TrimGalore (Krueger, 2012) (version 0.6.2) was used to trim adapters, low-quality ends from reads with Phred score less than 30 (-q 30), and discarding trimmed reads shorter than 30 bp (--length 30). Trimmed reads were then 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 the GENCODE (Harrow et al., 2012) annotation (version M18) with additional piRNA precursors, and transposable elements from repeat masker (concatenated by family), as in (Gapp et al., 2018) Downstream analysis: The downstream analysis was performed using 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). Genes that had less than 15 mapped reads in at least 40% of the samples were filtered out. Surrogate variable analysis was performed using thesva package (Leek, 2014) (version 3.32.1) on variance-stabilized data from the DESeq2 package (Love et al., 2014) (version 1.24.0), identifying 3 variables which were then included in the model for differential expression. Normalization factors were obtained using TMM normalization (Robinson and Oshlack, 2010) from the edgeR package (Robinson et al., 2009) (version 3.28.1) and differential gene expression (DGE) analysis was performed using the limma-voom (Law et al., 2014) pipeline from the limma (Ritchie et al., 2015) (version 3.42.2), adjusting for the three surrogate variables.

ATAC-seq
Quality control, alignment, and peak calling: Paired-end (PE) sequencing was performed on the PND15 and adult samples on the Illumina HiSeq2500 platform. The FASTQ files were assessed for the quality using FastQC (Andrews et al., 2012) (version 0.11.8). Quality control (QC) was performed using TrimGalore (Krueger, 2012) (version 0.6.2) in PE mode (-paired), trimming adapters, low-quality ends (-q 30) and discarding reads that become shorter than 30 bp after trimming (--length 30). Alignment of the QC data was performed using

Downstream analysis:
The downstream analysis was performed in R (R Core Team, 2020) (version 3.6.2), using packages from CRAN (https://cran.r-project.org) and Bioconductor (Huber et al., 2015). The peaks were annotated based on overlap with GENCODE (Harrow et al., 2012) (version M18) transcript and/or the distance to the nearest transcription start site (https://gist.github.com/dktanwar/f873af6c4acd8c38ee928aa6f3d5b81d). The number of extended reads overlapping in the peaks regions was calculated using the csaw package (Lun and Smyth, 2015) (version 1.20.0). The peak regions that had less than 15 mapping reads in at least 40% of the samples were filtered out. Normalization factors were obtained on the filtered peak regions using the TMM normalization method (Robinson and Oshlack, 2010) and differential analysis on the peaks was performed using the Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood (glmQLFit) Tests from the edgeR package (Robinson et al., 2009)

Figures
All the figures in this study were generated using the ggplot2 (Wickham, 2016) ggpubr (Kassambara, 2020) EnrichedHeatmap (Gu et al., 2018) ggforce (Pedersen, 2020) and ComplexHeatmap (Gu et al., 2016) packages, and using base plotting functions in R (R Core Team, 2020).  (B) K-means clustering (k = 5) of genes with differential expression in postnatal versus adult spermatogonial cells (cluster 1: n = 325 genes, cluster 2: n = 2064 genes, cluster 3: n = 5131 genes, cluster 4: n = 3122 genes, cluster 5: n = 1386 genes). Each row depicts a gene and each column depicts a sample. Gene expression fold-change is calculated by subtracting log2 counts per million (CPM) of PND15 and adults from PND8. Samples are clustered using Ward's method and genes are ordered by "MDS_angle method" using seriation (R package); (C) Dot plot of top 10 Gene Ontology (GO) Biological Processes (BP) terms enriched for each gene cluster selected based on the adjusted P. The dot color is representative of the number of genes annotated for each GO term; (D) Line plots depicting expression patterns of key genes at PND8, PND15 and in adult spermatogonial cells.

Figure 2. Chromatin accessibility changes between PND15 and adult spermatogonial cells
(A) Volcano plot of differentially accessible ATAC regions (adjusted P ≤ 0.05 and Log2FC ≥ 1) between PND15 and adult spermatogonial cells. More acc. in Adult = regions which gain in accessibility between PND15 and adult spermatogonial cells, Less acc. in Adult = regions which loose in accessibility between PND15 and adult spermatogonial cells; (B) Bar plot illustrating the genomic distribution of differentially accessible ATAC regions in adult spermatogonia compared to PND15 cells; (C) Dot plots of top enriched GO BP terms for regions with increased and decreased chromatin accessibility in adult spermatogonia. GO analysis was performed using GREAT. The size of the dot is representative of the number of genes annotated in the GO term, and the color corresponds to the adjusted P; (D) Dot plots of top enriched GO BP for more accessible regions located in gene body and in proximity of TSS (+/-1kb) in adult spermatogonia. GO analysis was performed using GREAT. The size of the dot is representative of the number of genes annotated in the GO term, and the color corresponds to the adjusted P. (C) Dot plots of top transcription factor motifs enriched in differentially accessible chromatin regions situated in gene bodies or in intergenic areas of the genome. Each dot corresponds to a motif. The differential gene expression of the corresponding transcription factor between PND15 and adult spermatogonia was extracted from the RNA-seq data, and is shown as color coded Log2 fold-change. The size of the dot indicates the expression adjusted P of each motif; (D) Heatmap of exemplary differentially expressed genes between PND15 and adult spermatogonia displaying motifs of JUND/cFOS and RFX1, RXRa and RXRg in their proximity. Each row depicts a gene and each column depicts a sample. Gene expression fold-change was calculated by subtracting log2CPM of Adult from PND15. Samples are clustered using Ward's method and genes are ordered by "PCA method" using seriation (R package); (E) Genomic snapshots of exemplary genes Gata2 and Eed showing relative abundance of: transcripts from RNA-seq, chromatin accessibility from ATAC-seq, 3 and TF motifs from HOMER analysis, in PND15 and adult spermatogonial cells. with H3K27ac (n = 32), Category 3 regions without any H3 mark (n = 25). Each line represents a peak region and the regions are ordered by the ATAC signal. Mid-x-axis corresponds to the middle of a peak region and is extended to +/-2.5 kbp. The colorkey of the ATAC and ChIP heatmaps represent ATAC and ChIP-seq signal respectively. For WGBS the color key shows the percentage of methylation. For RNA-seq, gene expression fold-change is calculated by subtracting log2 CPM of PND15 and Adults from PND8; (C) Genomic snapshots of exemplary genes from Category 1 (Mael) and Category 2 (Gata2 and Hmx1) showing relative abundance of: transcripts from RNA-seq, chromatin accessibility from ATAC-seq, 3 different histone marks (H3K27ac, H3K4me3 and H3K27me3) from ChIP-seq and TF motifs from HOMER analysis, in PND15 and adult spermatogonial cells.   AHR  ATF3  CREB1  DLX3  E2F3  E2F4  E2F6  EGR1  EGR2  ELF1  ELF2  ELK1  ELK4  ERG  ETS1  ETS2  ETV4  ETV6  FEV  FLI1  FOS  FOSB  FOSL1  FOSL2  FOXA2  FOXC1  FOXJ2  FOXJ3  FOXK1  FOXM1  FOXO1  FOXO3  FOXO4  FOXP2  FOXQ1  GABPA  ISL1  JUN  JUNB  JUND  LHX2  LHX6  MITF  NFIA  NFYA  NOBOX  NR1D2  NR5A2  OVOL2  PBX3  PDX1  PRRX2  RARA  RFX1  RFX2  RFX3  RXRA  RXRG  SOX3  SOX4  SP3  SP4  SPI1  TBX20  TFE3  USF1  USF2  VSX2  WT1 −5 0 5 Log 2 foldchange -Log10 adjusted P 2.5 5.0 7.5 10.0 12.5

More acc. in Adult Less acc. in Adult
Gene body Intergenic Gene body Intergenic