Decision letter: The single-cell chromatin accessibility landscape in mouse perinatal testis development

Abstract

Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Spermatogenesis depends on an orchestrated series of developing events in germ cells and full maturation of the somatic microenvironment. To date, the majority of efforts to study cellular heterogeneity in testis has been focused on single-cell gene expression rather than the chromatin landscape shaping gene expression. To advance our understanding of the regulatory programs underlying testicular cell types, we analyzed single-cell chromatin accessibility profiles in more than 25,000 cells from mouse developing testis. We showed that single-cell sequencing assay for transposase-accessible chromatin (scATAC-Seq) allowed us to deconvolve distinct cell populations and identify cis-regulatory elements (CREs) underlying cell-type specification. We identified sets of transcription factors associated with cell type-specific accessibility, revealing novel regulators of cell fate specification and maintenance. Pseudotime reconstruction revealed detailed regulatory dynamics coordinating the sequential developmental progressions of germ cells and somatic cells. This high-resolution dataset also unveiled previously unreported subpopulations within both the Sertoli and Leydig cell groups. Further, we defined candidate target cell types and genes of several genome-wide association study (GWAS) signals, including those associated with testosterone levels and coronary artery disease. Collectively, our data provide a blueprint of the ‘regulon’ of the mouse male germline and supporting somatic cells. Editor's evaluation This manuscript by Liao et al. aims to understand the genetic networks that underlie or modulate gonadogenesis and germ cell maturation during the fetal to neonatal transition. This goal was achieved by performing scATACseq on multiple timepoints (E18.5 and Postnatal days 1,2,5). Clustering of thousands of cells has striking cellular diversity and convincingly led to the identification of additional novel populations, of both germ cell and somatic origins. This is an important paper with far-reaching implications in reproductive biology, but additional validation would be needed to confirm the correlative observations and the functionality of newly identified testis cells. https://doi.org/10.7554/eLife.75624.sa0 Decision letter Reviews on Sciety eLife's review process Introduction Mammalian testis consists of germ cells and distinct somatic cell types that coordinately underpin the maintenance of spermatogenesis and fertility. These testicular cells display extensive developmental dynamics during the perinatal period. Primordial germ cells give rise to M-prospermatogonia at about embryonic day (E) 12, which enter G0 mitotic arrest at about E14 to form T1-prospermatogonia (T1-ProSG) (McCarrey, 2013). Shortly after birth, T1-ProSG resume mitotic activity and begin migrating from the center of the testis cords to the basal lamina of testicular cords, and become T2-prospermatogonia (T2-ProSG). Once resident at the basement membrane, T2-ProSG generate spermatogonia including self-renewing spermatogonial stem cells (SSCs) or directly transition into differentiating spermatogonia that participate in the first round of spermatogenesis (Kluin and de Rooij, 1981; Manku and Culty, 2015). Specialized somatic cells play a pivotal role in maintaining normal germ cell development and spermatogenesis. SSCs and their initial progenies reside in a niche on the basement membrane and are surrounded by Sertoli cells, which nourish SSCs. In mice, Sertoli cells actively proliferate during the neonatal period for 2 weeks (Vergouwen et al., 1991). Outside the seminiferous tubules, the interstitial compartment of testis mainly contains stroma, peritubular myoid cells (PTMs), Leydig cells, macrophages, and vascular cells. PTMs are smooth muscle cells that distribute over the peripheral surface of the basement membrane (Maekawa et al., 1996). They are mainly involved in tubule contractions to facilitate the movement of sperm to the epididymis and secrete extracellular matrix materials (Chen et al., 2014). Leydig cells are responsible for steroidogenesis and provide critical support for spermatogenesis. In mammals, there are two types of Leydig cells, fetal Leydig cells (FLCs) and adult Leydig cells (ALCs), which develop sequentially in the testis (Shima, 2019). While FLCs start to degenerate after birth, they are replaced by stem Leydig cells (SLCs), which are the progenitors of ALCs (Su et al., 2018). One of the goals of developmental biology is to identify transcriptional networks that regulate cell differentiation. Recent advances in single-cell transcriptomic methods have enabled an unbiased identification of cell types and gene expression networks in testis (Green et al., 2018; Tan et al., 2020). A key remaining question is how the gene expression is dynamically regulated during the development of distinct cell types. Since all cells share the same genetic information, cell-type development must be regulated by differential chromatin accessibility in a cell type-specific manner. Thus, the understanding of the testis development, especially the cell type-specific transcription factors (TFs), is of paramount importance. scRNA-Seq provides limited information of TFs, which are usually lowly expressed. Although sequencing methods such as ATAC-Seq and DNase-Seq have been developed for profiling chromatin accessibility landscapes across samples and classification of regulatory elements in the genome, the nature of bulk measurements masks the cellular and regulatory heterogeneity in subpopulations within a given cell type (Shema et al., 2019). Moreover, uncovering regulatory elements in testicular cell types has been particularly challenging as samples are limited and heterogeneous. Notably, only one previous study examined the genome-wide chromatin accessibility in Sertoli cells, using the Sox9 transgenic line followed by DNase-Seq (Maatouk et al., 2017). By profiling the genome-wide regulatory landscapes at a single-cell level, recent single-cell sequencing assay for transposase-accessible chromatin (scATAC-Seq) studies have demonstrated the potential to discover complex cell populations, link regulatory elements to their target genes, and map regulatory dynamics during complex cellular differentiation processes. We reasoned that the advent of new single-cell chromatin accessibility sequencing methods, combined with single-cell transcriptomic data, would be instrumental in advancing our understanding of gene regulatory networks in mammalian testis development. Here, we applied scATAC-Seq to deconvolve cell populations and identify cell type-specific epigenetic regulatory circuits during perinatal testis development. The dataset led to identification of key cell type-specific TFs, defined the cellular differentiation trajectory, and characterized regulatory dynamics of distinct cell types. Furthermore, our results shed light on the identification of target cell types for genetic variants. To enable public access to our data, we constructed the mouse testis epigenetic regulatory atlas website at http://testisatlas.s3-website-us-west-2.amazonaws.com/. Results Single-cell ATAC-Seq captures developmental and cell type-specific heterogeneity in the testis To delineate the dynamic changes on cellular populations in a developing testis, we profiled the chromatin accessibility landscapes of mouse perinatal testis across E18.5 and postnatal stages (P0, P2, and P5) by scATAC-Seq (Figure 1A). These time points were chosen to represent the diversity of cell-type compositions involved in the key developmental events in the testis (Figure 1B). Altogether, we profiled chromatin accessibility in 25,613 individual cells after stringent quality control filtration and heterotypic doublet removal (Figure 1—figure supplement 1). These samples showed no clustering based on covariates such as transcription start site (TSS) enrichment and number of fragments detected (Figure 1—figure supplement 2A). Figure 1 with 2 supplements see all Download asset Open asset Classification and identification of germ cells and somatic cells during perinatal testicular development. (A) Experimental design. The workflow of testis collection and single-cell sequencing assay for transposase-accessible chromatin (scATAC-Seq) to measure single nuclei accessibility on BioRad SureCell ATAC-Seq platform. (B) Illustration of the testicular microenvironment. GC: germ cell; SC: Sertoli cell; LC: Leydig cell; BV: blood vessel; BM: basement membrane; ST: seminiferous tubule; PTM: peritubular myoid cell. (C) Uniform manifold approximation and projection (UMAP) representations with cells colored by the gene score of marker genes for each cell type. (D) UMAP representation of cells captured from four time points. Cells are colored by predicted groups. (E) Bar chart showing the distribution of cells in each cluster for different time points. (F) Heatmap of 12,250 marker genes across cell types (FDR ≤ 0.05, Log2FC ≥ 0.2). Several clusters showed developmental stage specificity, which were made up almost entirely of cells from a single time point (Figure 1—figure supplement 2B). To improve cell-type annotation, we used Harmony to integrate datasets of four time points and project cells onto a shared embedding in which cells were grouped by cell type rather than developmental stage (Korsunsky et al., 2019). Unbiased iterative clustering of these single cells after integration identified 11 distinct clusters (Figure 1—figure supplement 2B). Some clusters could be assigned to known testicular cell types based on gene activity scores of key marker genes compiled from chromatin accessibility signals within the gene body and promoter (Figure 1C; Tan et al., 2020). While this approach provided broad classifications for cell-type annotation, an unbiased method is needed for more accurate classification. Therefore, we leveraged a previously published scRNA-Seq dataset of perinatal testis samples to predict cell types in scATAC-Seq data (Tan et al., 2020). We first re-analyzed scRNA-Seq data to determine the cellular composition and annotate cells based on their transcriptional profiles. Prediction of cell types in scATAC-Seq was then performed by directly aligning cells from scATAC-Seq with cells from scRNA-Seq through comparing the ‘query’ gene activity scores matrix with ‘reference’ gene expression matrix based on the top variable genes in the scRNA-Seq dataset (Supplementary file 1). The results showed that the vast majority of cells had a high prediction score and were confidently assigned to a single cell type (Figure 1D, Figure 1—figure supplement 2C). Cell-type proportions varied across time points, such as the expansion of germ cells during the early neonatal period (Figure 1E). We further validated the cluster assignment by gene score and chromatin accessibility profiles of marker genes (Figure 1F). Taken together, scATAC-Seq allowed the detection and assignment of cell identities in the developing testis. Chromatin accessibility defines cell types in developing testis Cell types can be distinguished based on whether differentially accessible chromatin regions (DARs) are ‘open’ or ‘closed’. After identifying 214,890 accessible chromatin regions in the scATAC-Seq library (Supplementary file 2), we investigated cell type-specific chromatin accessibility profiles. We compared differences in chromatin accessibility among cell types directly using Wilcoxon testing to identify DARs while accounting for TSS enrichment and the number of unique fragments per cell (Figure 2A, Supplementary files 3 and 4). Deconvolution of chromatin accessibility by cell types revealed that accessible sites are primarily located in the distal and intron region (>3 kb from TSS), suggesting an enrichment of gene regulatory elements (Figure 2—figure supplement 1A). Figure 2 with 1 supplement see all Download asset Open asset Characterization of differentially accessible regions and identification of cell type-specific transcription factors (TFs). (A) Heatmap of 51,937 differentially upregulated accessible peaks (FDR ≤ 0.01, Log2FC ≥ 2) across cell types. (B) Aggregated single-cell sequencing assay for transposase-accessible chromatin (scATAC-Seq) profiles of selected markers. (C) Heatmap of enriched motifs (FDR ≤ 0.1, Log2FC ≥ 0.5) across cell types. (D) TF footprints (average ATAC-Seq signal around predicted binding sites) for selected TFs. (E) Schematic of identifying positive TF regulators through correlating gene score (scATAC-Seq data)/gene expression (integrating scATAC-Seq and scRNA-Seq data) with TF motif activity (scATAC-Seq data). (F) Scatter plot of positive TF regulators (correlation >0.5, adjusted p-value <0.01). (G) Heatmaps of differential TF motif activity (left) and gene expression (right) of positive TF regulators in F. (H) TF overlay on scATAC uniform manifold approximation and projection (UMAP) of TF chromVAR deviations (top) and gene expression (bottom). We found cluster-specific DARs were associated with cell type-specific marker genes identified from scRNA-Seq (Figure 2B). For example, Amh is a marker gene in Sertoli cells, and it showed increases in both number and amplitude of ATAC peaks within its promoter and gene body. We further compared DARs to a previously published DNase-Seq experiments in bulk Sertoli cells and found that DNase I hypersensitive sites were clearly enriched in the Sertoli cell population in our scATAC-Seq (Figure 2—figure supplement 1B; Maatouk et al., 2017). These data confirm that scATAC-Seq is a robust method for the detection of cell type-specific chromatin accessibility. Chromatin accessibility is associated with cell type-specific TF activity Currently, the identities of cell type-specific TFs involved in testis development are poorly defined. Accessibility at regulatory sites is driven by TF binding and histone modifications of local chromatin (Cui et al., 2013). To characterize the determinants of chromatin accessibility variation among cell types, we predicted TF ‘activity’ for individual cell types based on the presence of binding motifs within DARs. Assessment of enriched TFs and their cognate motifs identified several known cell type-specific regulators – including the nuclear receptors (NR4A1 and NR5A1) in Sertoli cells and Leydig cells, ESR2 in Leydig cells, MYOG in PTMs, and previously uncharacterized TFs as potential cell type-specific regulators (Figure 2C). For example, we found that ZEB1 and SNAI2 motifs were enriched in germ cells, indicating they may undergo mesenchymal-like transition in perinatal development (Hammoud et al., 2015; Liao et al., 2020). DNA bound by TFs is protected from transposition by Tn5, which can be visualized by plotting the ‘footprint’ pattern of each TF as the local chromatin accessibility surrounding the motif midpoint. Examining the footprint validated the cell type-dependent differential footprint occupancy of identified TFs (Figure 2D). Although motif enrichment for DAR can be informative, this measurement is not calculated on a per-cell basis and they do not take into account the insertion sequence bias of Tn5. Therefore, in the second analysis approach, we used chromVAR to infer TF motif activity, which can reflect the enrichment level of the TF motif in accessible regulatory elements in an scATAC-Seq dataset. We first identified deviant TF motifs by stratifying motifs based on the degree of variation observed across clusters. Since TFs from the same family often share a similar motif, this makes it challenging to identify the specific TFs that actually drive the observed changes in chromatin accessibility. This is illustrated by SOX9, which exhibited enriched TF activity in both Sertoli and germ cell clusters, despite being known to be expressed only in Sertoli cells. This could be attributed to the expression of other SOX family genes, such as SOX3, which share similar DNA binding motifs, in germ cells (Raverot et al., 2005). However, SOX9 exhibited concordant enrichment of TF activity and gene expression only in Sertoli cells and can be considered a strong candidate in Sertoli cells, but not in germ cells. To reduce false discovery, we systematically identified putative positive regulators determined from the correlation between the gene expression based on scRNA-Seq dataset (or inferred gene activity based on scATAC-Seq dataset) and the chromVAR motif activity score, reasoning that expression of high-confidence TFs is correlated with their motif accessibility (Figure 2E). Clustering analysis of positive regulators showed that diverse combinatorial TF motif landscapes were apparent across cell types and closely mirrored gene accessibility profiles of respective TFs (Figure 2F–H). There was an increased GATA1 TF ‘activity’ (motif activity) in the Sertoli cell cluster, in addition to increased chromatin accessibility in Gata1 (gene activity) and increased Gata1 transcription (gene expression) (Figure 2G). It has been shown that mutation of GATA1 causes human cryptorchidism (Nichols et al., 2000) and its expression in Sertoli cells is conserved between human and mouse (Yomogida et al., 1994). A similar pattern was seen for DMRT1 in germ cells, which is one of the top motifs enriched in human SSC-specific ATAC-Seq peaks (Guo et al., 2017). This analysis also revealed shared and unique regulatory programs across cell types. For example, PTM and stromal cells shared similar regulators, but PTM demonstrated higher activity of AR and TCF21 (Figure 2G). Similarly, NR5A1 and GATA4 were more active in both Leydig cell and Sertoli cell populations. However, only Sertoli cells showed increased activity in the SOX family. Importantly, we observed that individual cell types can be defined by TF ‘activity’, suggesting that cell type-specific TFs likely regulate chromatin accessibility. Collectively, these results are indicative of robust inference of TF activity at the level of single cells and reveal TF dynamics central to cis-regulatory specification of diverse cell states. Chromatin accessibility is associated with cell type-specific chromatin interaction networks As enhancers play a critical role in establishing tissue-specific gene expression patterns during development, we predicted that active enhancers would be enriched around lineage-specific genes. To test this, we used an analytical framework to link distal peaks to genes in cis, based on the coordination of chromatin accessibility and gene expression levels across cells (Figure 3A). We identified 35,245 peak-to-gene links by correlating accessibility changes of ATAC peaks within 250 kb of the gene promoter with the mRNA expression of the gene from scRNA-Seq. Some of these peak-to-gene links are likely to be promoter-enhancer regulatory units, as 3262 regions overlapped with previously identified testis enhancers (Gao and Qian, 2020; Figure 3—figure supplement 1A). Figure 3 with 2 supplements see all Download asset Open asset Chromatin interaction networks in different cell types. (A) Schematic for identifying significant peak-to-gene links by correlating accessible peaks (single-cell sequencing assay for transposase-accessible chromatin [scATAC-Seq] data) to gene expression (integrating scATAC-Seq data and scRNA-Seq data). (B) Heatmaps of peak accessibility (left) and gene expression (right) of 22,545 peak-to-gene linkages across cell types. (C) Aggregated scATAC-Seq profiles showing peak-to-gene links to the Sox9 locus overlapped with known enhancer regions. (D) Aggregated scATAC-Seq profiles showing peak-to-gene links to the Wt1 locus overlapped with known enhancer regions. (E) Aggregated scATAC-Seq profiles showing peak-to-gene links to the Dlk1 locus overlapped with SNP. (F) Sorting strategy for isolation of DLK1- and DLK1+ cells from P6 whole testis. The majority of DLK1+ cells are located in P1 (upper left). The DLK1-/+ population was gated using Red-X-labeled sample compared with unstained control (lower left). RT-PCR analysis (right) of relative expression of peritubular myoid cell (PTM) marker (Acta2), Leydig cell marker (Cyp11a1), and stromal cell marker (Igf1) of DLK1-/+ cells compared with whole testis sample (p<0.001, n=3, one-way ANOVA). Gapdh was used as endogenous control. Error bars are plotted with SD. (G) Representative confocal images of testis sections from Oct4-GFP transgenic mice at P6. Stromal cells (asterisks) and some PTMs (arrowheads) are positive for DLK1 (red). Oct4-GFP indicated germ cells. Cell nuclei were stained with DAPI. Scale bar = 50 μm. To identify putative cis-regulatory elements (CREs) specific to each cell type, we performed clustering analysis, with each cluster containing peak-to-gene links enriched in one or two specific cell types (Figure 3B). Gene Ontology (GO) analysis of the targets of peak-to-gene links in each cluster confirmed that they were highly enriched in terms related to regulations of each cell type (Figure 3—figure supplement 1B). The full list of peak-to-gene links in each cluster can be found in Supplementary file 5. We next examined whether we can use this information to link DARs to known cell type-specific enhancers. During male sex determination, Sry activates male-specific transcription of Sox9 in the male genital ridge via the testis-specific enhancer core element (TESCO) enhancer (Sekido and Lovell-Badge, 2008). Comparison of the genomic region around Sox9 among all cell types revealed a region 13 kb upstream formed a peak-to-gene link with the Sox9 TSS, which is unique to the Sertoli cell population and overlapped with the 3 kb TES enhancer including the TESCO elements (Figure 3C). Enhancer activity was previously narrowed to a subregion of TES: the 1.3 kb TESCO element (Sekido and Lovell-Badge, 2008). Besides known enhancers, our data linked an additional three elements located 9 kb 5′, 21 kb 3’, and 68 kb 3’ to Sox9 representing novel regulatory elements of Sox9 gene regulation. Our analysis also successfully revealed a previously identified functional enhancer as a novel candidate to regulate Sertoli cell marker Wt1 (Figure 3D). After confirming that our approach can be used to reveal functionally relevant regulatory regions, we next identified putative regulatory elements for important cell type-specific regulators in each cell type, including Nanos2 and Uchl1 in germ cells, Cyp11a1 and Nr5a1 in Leydig cells, Tpm1 and Socs3 in PTMs (Figure 3—figure supplement 2), and Dlk1 in stromal cells and PTMs (Figure 3E). Although DLK1 is considered as a marker for immature Leydig cells in human (Lottrup et al., 2014), the Dlk1-Gtl2 locus demonstrated preferential accessibility in stromal cells and PTMs. This coincided with the high number of peak-to-gene links including the largest mammalian miRNA mega-cluster located approximately 150 kb downstream (Seitz et al., 2004). Therefore, we examined the expression of Dlk1 in mouse testicular cells. We sorted DLK1+ mouse testicular cells using fluorescence-activated cell sorting and observed higher expression of Igf1 mRNA compared to whole testis samples. This suggests that the sorted cells were enriched for stromal cells, since Igf1 is commonly used as a marker for this cell type (Figure 3F). Immunostaining of neonatal testis tissue demonstrated that stromal cells and PTMs were positive for DLK1 (Figure 3G). In conclusion, our results highlight the occurrence of diverse cell type-specific regulatory configurations among CREs and their target genes in the testis. Stage-specific TF regulators and chromatin co-accessibility during gonocyte to spermatogonia transition Next, we analyzed the chromatin accessibility characteristics of the individual germ cell subsets in our datasets. Since the goal was to reveal developmental dynamics, we did not perform ‘harmony’ integration since we didn’t want to remove the contribution of developmental stage-of-origin from the embedding. Re-clustering of germ cells from the E18.5, P0, P2, and P5 testicular datasets revealed seven-cell clusters (Figure 4A and Figure 4—figure supplement 1A). Notably, germ cells from E18.5 and P0 are largely clustered together and occupy clusters GC1 and GC3, indicating a minimal change of chromatin accessibility before and after birth. In contrast, P2 cells are present in GC2 and P5 cells occupy the remaining clusters. Figure 4 with 4 supplements see all Download asset Open asset Identification of germ cell clusters during the perinatal period. (A) Uniform manifold approximation and projection (UMAP) representation of germ cells. Cells are colored by time points (left) and clustering based on constrained integration with scRNA-Seq data (right). (B) Heatmaps of differential transcription factor (TF) motif activity (left) and gene activity (right) of positive TF regulators across cell clusters (correlation >0.5, adjusted p-value <0.01). (C) TF overlay on scATAC UMAP of gene expression (top) and TF chromVAR deviations (bottom) for positive TF regulator examples in B. (D) Representative confocal images of immunostaining on sections from P6 testis demonstrate that a subset of germ cells (TRA98+) express the Sertoli cell marker NR5A1 (arrowhead), while the majority of germ cells are NR5A1-negative (arrow). Scale bar = 50 μm. (E) Gene activity of Ngn3 shown in UMAP. Deconvoluting cell states using scATAC-Seq measurements alone is difficult within a single cell type. Therefore, we integrated the germ cell subsets based on scATAC-Seq data with the published perinatal testis scRNA-Seq dataset (Figure 4—figure supplement 1B). The prediction scores of individual cells were overall high, indicating the cluster identity assignment was reliable (Figure 4—figure supplement 1C). This prediction revealed four clusters of known developmental stages, T1-ProSG (T1), T2-ProSG (T2), undifferentiated spermatogonia (Undiff), and differentiation-primed spermatogonia (Diff), together with two clusters with unknown identity (Figure 4A). We first identified the TFs important for each cluster, revealing 57 putative positive TF regulators in germ cell development (Figure 4B and C). For instance, consistent with previous studies, Foxo1 and Dmrt1 exhibited increased TF motif activity and gene expression in undifferentiated spermatogonia (Goertz et al., 2011; Matson et al., 2010). E2f4 shows highest activity in differentiation-primed spermatogonia, and is known to be critical for the development of the male reproductive system (Danielian et al., 2016). Interestingly, this analysis identified Nr5a1 as the positive TF regulators in the unknown clusters (Figure 4B and C, Figure 4—figure supplement 1D). Since Nr5a1 is widely considered as a somatic cell marker in testis (Luo et al., 1994), we performed immunostaining to examine NR5A1’s expression in germ cells (Figure 4D). Indeed, a subset of germ cells expressed NR5A1, which ruled out the possible contamination from somatic cells. Interestingly, previous scRNA-Seq of neonatal pro-spermatogonia identified a spermatogonial signature cluster showed high levels of mRNAs characteristic of Sertoli cells, including Nr5a1, Sox9, and Wt1 (Hermann et al., 2015). Additionally, the expression of the somatic cell marker WT1 has been observed in some germ cells through immunostaining (Wen et al., 2021). These observations reinforced our findings that cells with germ cell identity can express somatic cell genes (Figure 4—figure supplement 1D). Our results also revealed several TF candidates regulating T1-ProSG, which have previously been difficult to identify due to technical challenges in isolating this cell population, such as Ybx2, Smad3, and Msc. In line with its role in testicular development, Gata6 is upregulated in T1- and T2-ProSG (Padua et al., 2015). We then aimed to reconstruct the differentiation trajectories by ordering the clusters with developmental stages predicted by scRNA-Seq integration. Two possible trajectories were observed, as germ cell differentiation appeared to diverge at P0 via two distinct branches. The first trajectory represents the differentiation fit into the conventional model as it charts a trajectory from gonocyte to undifferentiated and then differentiating spermatogonia in P5. The second path bypassed the undifferentiated state but passed through the unknown populations and directly reached the differentiating state by P5. It has been reported that the first round of spermatogonia arise from a unique neurogenin-3 (Ngn3) negative pool of ProSG that transitions directly into A1 spermatogonia (Yoshida et al., 2006). Interestingly, the unknown population (Unknown-2) displayed the lowest level of Ngn3 gene activity, which raised the possibility that the second trajectory represents the origin of the first wave of spermatogenesis (Figure 4E). To determine the key genes driving the spermatogonial development in the first trajectory, we generated a pseudotime trajectory and uncovered a list of genes with dynamic changes (Figure 4—figure supplement 1E and F). The pattern of TF dynamics suggested a model of differentiation as a transition between two phases involving progressive loss of gonocyte-specific TF activities and gradual increase of TF activities relevant to spermatogonia. We observed that the motif binding activity of Id4 was initially high but then declined after birth, while that of ETS and Sp/KLF family members increased in spermatogonia (Figure 4—figure supplement 1E). To further reveal TFs that drive the germ cell development, we pruned the data by correlating gen