Transcriptome dynamics between early postnatal and adult SPGs.

(A) Schematic representation of the experimental strategy to isolate and analyse SPGs from postnatal or adult male mice.

(B) Upper panel: Representative dot plots of the sorting strategy for enrichment of postnatal and adult SPGs. Gating based on side scatter/forward scatter (SSC-A/FSC-A) and forward scatter (height/forward scatter) (FSC-H/FSC-A) area was conducted to exclude cell debris and clumps. Lower panel: Representative immunocytochemistry on unsorted and sorted PND15 cells with anti-PLZF antibody and DAPI (nuclei) showing enrichment of PLZF+ cells after FACS.

(C) Heatmap of expression profile of selected markers of SPGs and different testicular somatic cells extracted from total RNA-seq data on PND8, PND15 and adult samples (n=6 for each group). Key genes for stem cell potential, stem and progenitor spermatogonia and Leydig and Sertoli cells were chosen to assess the enrichment of SPGs in the sorted cell populations. Each row in the heatmap represents a biological replicate from two experimental batches (B1 and B2). At the bottom are shown Integrative Genomics Viewer (IGV) tracks for aggregated RNA-seq signal for PND8, PND15 and adult over key genes used as markers. Gene expression is represented as Log2CPM (counts per million).

(D) Volcano plot of differentially expressed genes (DEGs) (adjusted-P≤0.05 and absolute Log2FC≥1) between PND8 and PND15 (left) and PND15 and adult (right). Figures in grey boxes indicate the number of down- and up-regulated DEGs in each comparison.

(E) Genomic IGV snapshots of exemplary DEGs showing aggregated RNA-seq signal for PND8, PND15 and adult SPGs.

(F) Left: Bar-plot of GO categories enriched in DEGs between PND8 and PND15 (adjusted-P≤0.05). Dotted line indicates threshold value for significance of 0.05. Right: Heatmap of DEGs between PND8 and PND15 belonging to the GO category “Transcription”. Each row represents a biological replicate from two experimental batches. Shown are Log2FC with respect to average at PND8.

(G) Left: Bar-plot of GO categories enriched in DEGs between PND15 and adult (adjusted-P≤0.05). Dotted line indicates threshold value for significance of 0.05. Right: Heatmap of DEGs between PND8 and PND15 belonging to the GO category “Spermatogenesis”. Each row represents a biological replicate from two experimental batches. Shown are Log2FC with respect to average at PND8.

(H) Left: Heatmap of all DEGs specific to adult SPGs. Each row represents a biological replicate from two experimental batches (B1 and B2). Shown are Log2FC with respect to average at PND8. Right, Bar-plots of GO categories enriched in up-regulated (adjusted-P≤0.05) (top) or down-regulated (adjusted-P≤0.05) (bottom) genes in adult SPGs. Dotted line indicates threshold value for significance of adjusted-P=0.05.

Dynamics of chromatin accessibility of SPGs during the transition from postnatal to adult stage (A) Volcano plot of DARs obtained by ATAC-seq (adjusted-P≤0.05 and absolute Log2FC≥1) between PND15 and adult SPGs. Figures in grey boxes indicate the number of DARs with decreased (Down) or increased (Up) chromatin accessibility in PND15 compared to adult SPGs.

(B) Heatmaps showing normalized ATAC-seq signal for all identified DARs comparing PND15 and adult stages. Each row represents a 2kb genomic region extended 1kb down- and upstream from the centre of each identified DAR (as shown in A). Rows are ordered in a decreasing level of mean accessibility.

(C) Bar plot illustrating the genomic distribution of DARs between PND15 and adult SPGs.

(D) Bar plot of GO categories for associated genes assigned by proximity to all DARs. Dotted line indicates threshold value for significance of adjusted-P=0.05.

(E) IGV tracks for ATAC-seq signal for PND15 and adult SPGs showing DARs (red box) located at promoters and intergenic regions.

(F) Heatmaps showing normalized ChIP-seq signal for different histone marks in adult SPGs (public data derived from Cheng et al., 2020) at all identified ATAC-seq peaks from PND15 and adult stages. Each row in the heatmap represents a 2kb genomic region extended 1kb down- and up-stream from the centre of each identified peak. Shown data correspond to non-regenerative SPGs as stated in Cheng et al., 2020.

(G) Pie charts showing the overlap of all identified DARs with genomic regions significantly enriched for different histone marks in adult SPGs (public data derived from Cheng et al., 2020).

TF binding motifs at DARs and their transcriptional dynamics during the transition from postnatal to adult SPGs

(A) Table of top five enriched motifs in genomic regions with decrease in chromatin accessibility in adult SPGs and TF matching motifs. A full list of TF motifs is provided in Table S4.

(B) Heatmap of expression profile of TFs with motifs in genomic regions with decreased chromatin accessibility shown in (A). Each row represents a biological replicate from two experimental batches. Unsupervised hierarchical clustering was applied to each row and a dendrogram indicating similarity of expression profiles among genes is shown. Shown are Log2FC with respect to average of PND8 SPGs.

(C) Line-plots of average expression of each gene displayed in the heatmap in (B) showing the dynamics of gene expression across PND8, PND15 and adult SPGs.

(D, E) Genomic snapshot from IGV showing aggregated RNA-seq signal from PND8, PND15 and adult SPGs for the TFs D) Dmrt1 and E) Elk4.

(F) Table of top five enriched motifs in genomic regions with increase in chromatin accessibility in adult SPGs and TF matching motifs. A full list of TF motifs is provided in Table S4.

(G) Heatmap of expression profile of TFs with motifs in genomic regions with increase in chromatin accessibility shown in (F). Shown are Log2FC with respect to the average of PND8 SPGs.

(H) Line-plots of average expression of each gene displayed in the heatmap in (G) showing the dynamics of gene expression across SPGs development.

(I, J) Genomic snapshot from IGV showing the aggregated RNA-seq signal from PND8, PND15 and adult SPGs for the TFs I) Jun and J) Rfx2.

Chromatin accessibility dynamics at DEGs during the transition from postnatal to adult SPGs

(A) Top, profile plots of ATAC-seq signal from PND15 and adult SPGs over all ATAC-seq peaks located within a genomic region of 5kb surrounding the TSS of DEGs. Mean (line) and bootstrap confidence interval (area, shadow lines) computed over all non-overlapping 50bp genomic regions. n.s, non-significant differences at TSS region. *, significant difference in signal between PND15 and adult at p-value<0.05. PC, peak center from ATAC-seq data. Profile plots for all ATAC-seq peaks surrounding the TSS of down-regulated (top left) or up-regulated (top right) DEGs. Middle, profile plots of ATAC-seq signal from PND15 and adult SPGs over all ATAC-seq peaks located within a genomic region of 500bp surrounding the TSS of DEGs. Bottom, heatmaps showing normalized ATAC-seq signal from PND15 and adult for all ATAC-seq peaks surrounding the TSS of all down-regulated (left) or up-regulated (right) DEGs. Each row in the heatmap represents a 2kb genomic region extended 1kb down- and up-stream from the center of each identified peak. Row Are ordered in a decreasing level of mean accessibility taking PND15 as reference.

(B) Violin plots of average ATAC-seq signal for PND15 and adult SPGs over gene body of down-regulated (left) or up-regulated (right) DEGs. ****p<0.0001.

(C) TF motif enrichment analysis for all ATAC-seq peaks surrounding the TSS of all down-regulated (left) or up-regulated (right) DEGs.

(D) Left, pie chart of percentage of DARs overlapping the extended promoter of DEGs (+/- 5kb from TSS). Right, pie chart of percentage of DEGs associated with a DAR.

(E) IGV tracks for an example of DAR associated to the extended promoter (+/- 5kb from TSS) of a DEG. Shown are tracks for RNA-seq and ATAC-seq signal for PND15 and adult SPGs generated in this study and signal track for H3K4me3 from SPGs derived from public ChIP-seq data (Cheng et al., 2020).

(F) Left, pie chart of percentage of DARs overlapping potential regulatory sequences of DEGs (+/-100kb from start and end of gene). Right, pie chart of the percentage of DEGs associated with a DAR.

(G) IGV tracks for an example of DAR associated with potential regulatory sequences (+/-100kb from start and end of gene) of a DEG.

Differential chromatin accessibility at TEs in PND15 and adult SPGs.

(A) Heatmap of LTR and LINE subtypes with decreased accessibility between PND15 and adult SPGs extracted from Omni-ATAC data (adjustedP≤0.05 and Log2FC≥0.5)

(B) Heatmap of LTR and LINE subtypes with increased accessibility between PND15 and adult SPGs extracted from Omni-ATAC data (adjustedP≤0.05 and Log2FC≥0.5). Expression changes of these subtypes between adult and PND15 SPGs RNA-seq from literature are represented as Log2FC.

(C) Bar plot illustrating the genomic distribution of DARs between PND15 and adult SPGs.

(D) Genomic snapshot from IVG showing aggregated ATAC-seq and RNA-seq signal from PND8, PND15 and adult SPGs for Lncenc1 and Platr14. Below signal tracks is shown the RepeatMasker track for the mouse genome, and TEs with changes in chromatin accessibility are highlighted.