Revealing cellular and molecular transitions in neonatal germ cell differentiation using single cell RNA sequencing

In mammals, male fertility is sustained by male germline stem cells (GSCs) known as spermatogonial stem cells (SSCs) that either self-renew or differentiate to progenitors to initiate spermatogenesis (Oatley and Brinster, 2008). The current model suggests that the initial pool of SSCs is derived from the transition of postnatal gonocytes arising from primordial germ cells (PGCs) during embryogenesis, which enter a quiescent state after embryonic day 13.5 (Kluin and Derooij, 1981). Shortly after birth, gonocytes resume mitotic activity and begin migrating to the basal lamina of testicular cords. Once they reside at the basement membrane, they become spermatogonia, including self-renewing SSCs (Culty, 2009; Manku and Culty, 2015; Saitou and Yamaji, 2012). Violating this cell fate decision process can lead to serious clinical consequences (Kamisawa et al., 2012; Rajpert-de Meyts and Hoei-Hansen, 2007). SSC self-renewal and differentiation must also be balanced to maintain normal spermatogenesis and avoid tumorigenesis (Ishii et al., 2012).

Although much attention has been dedicated to the characterization of adult SSCs, the molecular relationship between the different germ cell populations at the neonatal stage has not been sufficiently addressed. The main challenge is that the germ cells in neonatal testes are not developmentally synchronized, resulting in the co-existence of germ cells at different developmental stages (Nagano et al., 2000). Moreover, there is significant heterogeneity regarding marker expression among the germ cell populations, where different markers are restricted to distinct or overlapping subsets of cells (Hermann et al., 2015; Nakagawa et al., 2007; Suzuki et al., 2009; Yoshida et al., 2006). Recently, another layer of heterogeneity in marker expression also received substantial attention: SSC markers can display a gradient of expression level and define different cellular states. For example, adult SSCs are present in populations of cells expressing a high level of ID4, whereas progenitor populations express relatively lower levels of ID4 (Helsel et al., 2017; Komai et al., 2014). To understand this further, it will be important to clarify whether expression gradients can be applied to other markers. OCT4 (POU5F1) is of particular relevance as it can exert a dose-dependent effect in mouse embryonic stem cells, whereby different expression levels of OCT4 define differentiation, de-differentiation or self-renewal (Zeineddine et al., 2014).

Single cell RNA sequencing (RNA-Seq) technology is increasingly used to deconstruct cellular subpopulations and transcriptional heterogeneity. Here, we profiled single cell transcriptomes of germ cells from mouse testes on postnatal day (PND) 5.5. We provide a snapshot of the cellular composition of undifferentiated germ cells, reveal global patterns of gene expression heterogeneity, and present potential novel transcriptional regulators during differentiation. Further analysis reveals a subpopulation within the most primitive spermatogonia labeled by CD87 (also known as uPAR and PLAUR), as well as a differentiation state marked by higher Oct4 expression and the potential role of OCT4 in regulation of the MAPK and cell adhesion pathway.

To resolve cellular heterogeneity in early germ cell development, we focused on PND5.5 when germ cells consist of gonocytes and primitive spermatogonia, including SSCs and differentiating spermatogonia. We found that a portion of Oct4-GFP⁺ cells resided within the center of testicular cords or were positioned adjacent to the basement membrane with protrusions towards the periphery, which confirmed the existence of transitioning gonocytes yet to finish the migratory phase at PND5.5 (Fig. 1). Oct4-GFP⁺ cells exhibit significant overlap with PLZF (ZBTB16), a marker for the undifferentiated state (Fig. 1A,D), suggesting that Oct4-GFP⁺ cells indeed label most, if not all, germ cells at this stage. Co-expression of Oct4-GFP with the self-renewal factor GFRA1 was also observed in a subset of Oct4-GFP⁺ cells (Fig. 1B,D). In addition, many Oct4-GFP⁺ cells begin to express the differentiation marker KIT, which suggests that PND5.5 testis contains differentiating spermatogonia (Fig. 1C). Altogether, Oct4-GFP⁺-labeled germ cells at this single time point should represent a mixture of germ cells at different cellular states and single cell transcriptome profiling of these cells would allow reconstruction of dynamic gonocyte-to-spermatogonia transition (GST) and subsequent differentiation in vivo (Fig. S1A).

To maximize the representation of the most undifferentiated germ cells, we isolated Oct4-GFP⁺/KIT⁻ cells from PND5.5 testis by fluorescence-activated cell sorting (FACS) (Fig. S1B-F). We used isotype control to confirm specific sorting (Fig. S1B). qRT-PCR analysis of Kit mRNA expression also showed a more than 10-fold increase in the KIT⁺ sorting fraction, which confirms the specific enrichment of Oct4-GFP⁺/KIT⁻ cells (Fig. S1F). The sorted Oct4-GFP⁺/KIT⁻ cells were subjected to Fluidigm C1 single cell RNA-Seq (Fig. 1E, Fig. S2, Tables S1 and S2). The germ cell markers Dazl and Ddx4 were abundantly expressed in all cells, indicating no contamination from somatic cells (Fig. S3A). A significant portion of germ cell-related genes, such as Bcl6b, Egr2, Sohlh1 and Csf3r, showed substantial variability, which revealed the significant gene expression heterogeneity in early male germ cell populations (Fig. S3A,B).

We next sought to interrogate the heterogeneous cellular composition. Principal component analysis (PCA) was used to identify the top 500 most variable genes that accounted for the majority of the variance in expression (Table S3) (Pollen et al., 2014). Unbiased clustering with this reduced gene list revealed three clusters of cells (Fig. 2A): one main cluster containing 56 cells, and two smaller clusters (I and IV: nine cells and six cells, respectively). It is notable that Cluster IV cells lack expression of several self-renewal genes (e.g. Gfra1, Etv5 and Ccnd2) (Fig. 2A), which were named as differentiation-primed spermatogonia. The aforementioned characterization (Fig. 1) implies that the majority of cells captured at PND5.5 should be undifferentiated spermatogonia and may be linked to the main cluster. This cluster can be further divided into two subclusters (Cluster II and Cluster III) based on the dendrogram (Fig. 2A). To reconstruct the transitional dynamics and ontogeny relationship of the cell subpopulations, we used the single cell analysis tool Monocle to order cells along a hypothetical timeline (pseudotime) of development (Trapnell et al., 2014). The resulting trajectory path began with Cluster I, passing through undifferentiated spermatogonia (Cluster II), then Cluster III, and ending at Cluster IV (Fig. 2B, upper panel). This indicated that Cluster II might sit at an earlier developmental stage than Cluster III in the pseudotime. Consistently, our PCA results implied that the first component represents the developmental transition from less differentiated cells (Cluster II) to the differentiation-primed spermatogonia (Cluster IV) on the x-axis (Fig. 2C). The PC scores of the top 500 genes revealed SSC self-renewal genes (Etv5 and Ccnd2) were expressed at higher levels in Cluster II cells, whereas differentiation regulators (Sohlh1 and Trib3) were expressed at higher levels in Cluster III cells (Fig. 2D). The majority of cells with high Id4 expression fell into Cluster II {15 out of 23 cells, log2[transcripts per million (TPM)>10]}, consistent with the observation that cells with a high level of Id4 represent true SSCs (Helsel et al., 2017). Erbb3 and Pax7 were also preferentially expressed in Cluster II (Fig. 2E). Together with the Monocle ordering result, these data indicate that Cluster II cells are likely to be the more stem cell-like spermatogonia (hereafter referred to as ‘primitive stem cells’) and Cluster III cells are more progenitor-like spermatogonia (hereafter referred to as ‘progenitors’). Cluster I corresponds to gonocytes, based on their earliest pseudotime stage and the fact that the existence of a gonocyte population should be reflected in single cell transcriptomes of the captured cells (Fig. 2B, upper panel). It is notable that Cd81, which mediates migratory gonocyte movement (Tres and Kierszenbaum, 2005), showed higher expression in Cluster I (mean TPM: Cluster I=62, cluster II=18, Cluster III=22, Cluster IV=10) (Fig. 2E). Monocle analysis also revealed three major cell states, which corresponded to GST (State 1), SSC self-renewal (State 2) and differentiation priming (State 3) (Fig. 2B, lower panel).

To validate the single cell RNA-Seq result, we sorted the second batch of PND5.5 Oct4-GFP⁺/KIT⁻ germ cells and measured the expression of a panel of genes that were preferentially expressed in each subpopulation using single cell RT-qPCR (Fig. S4A-C, Table S4). In line with the single cell RNA-Seq results, four clusters were recovered from hierarchical clustering (Fig. S4D) and the PCA plot showed a similar structure with the transition from stem cell subpopulation (green) to progenitor (red) and then to differentiation-primed cells (orange) (Fig. S4E). Similarly, the gene expression data generally recapitulated patterns we identified in the RNA-Seq experiment (Fig. S4F).

To investigate the transcriptional signature of the gonocyte population, we first focused on the cells in GST state (State 1) and compared gene expression between the first half (nine cells from Cluster I) and the second half (nine cells from Cluster II) in the trajectory (Fig. 3A). Top-ranked differentially expressed genes identified were postulated to regulate cell migration. For example, this analysis revealed Cd81 as one of the top genes (Table S5). Cluster I also expressed higher Crk, which encodes the adapter protein of PDGFR and might be involved in PDGF-mediated gonocyte survival and migration (Fig. 3A) (Basciani et al., 2008). Moreover, integrin-binding associated genes, such as Nf2, Itgb1, Utrn, Mfge8 and P4hb, are also upregulated in the gonocytes (Fig. 3A). In contrast, the gonocyte population expressed a lower level of Tnfaip1 (also known as Bacurd2) (Fig. 2D, Table S5), which influences the directional cellular migration through the regulation of RHOA protein turnover (Sailland et al., 2014). Interestingly, top-ranked genes enriched in migratory gonocytes included Trpm7, indicating its possible role in gonocyte migration through calcium-mediated cell movement (Chen et al., 2010).

Gene set enrichment analysis (GSEA) further revealed that gene sets associated with cell cycle and protein synthesis (aminoacyl-tRNA biosynthesis) were under-represented in the gonocyte population (Fig. 3B). Interestingly, gonocytes were shown to express lower levels of basal transcriptional factors (TFIIB) such as Taf4b, which encodes a chromatin-associated factor and has been demonstrated to regulate early germ stem cell specification (Falender et al., 2005). In line with the presence of 12.6% (9/71) cells (Cluster I) within the total cells captured in the single cell RNA-seq experiment, we observed a similar fraction that was negative for Ki-67 (Mki67) expression (Fig. 3C). Immunostaining also revealed that the centrally located cells contained more Ki-67-negative cells than peripherally located cells (Fig. 3D,E). These results implied a resting phase of these gonocytes prior to their transition to spermatogonia.

Lastly, gene ontology (GO) analysis identified a significant upregulation of genes associated with the death receptor pathway in Cluster I cells, such as Fadd, Tnfrsf1b, Icam1, Madd, Stx4a and Skil. We also found elevated expression of genes involved in intrinsic apoptosis signaling (Itm2b) (Fig. 3F, Table S5) (Fleischer et al., 2002). Genes associated with the DNA damage/repair pathway were shown to be repressed in gonocytes, including Ercc1, Trp53 and Bre (Fig. 3F, Table S5). These results are in line with the notion that apoptosis plays an integral role in gonocyte cell fate determination (Wang et al., 1998).

Collectively, our analyses revealed transcriptional signatures that might be involved in cell fate decision and migration of gonocytes during GST.

We next investigated the underlying molecular features of the progression from the stem-like stage to the progenitor state of spermatogonia. We identified 511 genes differentially expressed when comparing Cluster II with Cluster III (Table S6). GO analysis revealed significant enrichment for signatures associated with integrin complex, Rap1 signaling and lipid metabolism in primitive stem cell pools (Fig. S5A, Table S6), whereas TGF signaling genes (e.g. Smad3 and Creb1) were over-represented in progenitor pool (Fig. S5B, Table S6). GSEA analysis further revealed that ‘JAK/STAT signaling’, ‘cell adhesion molecules (CAMs)’,’ neuroactive ligand-receptor interaction’, and ‘RIG-I-like receptor signaling pathway’ were enriched in the stem cell cluster, whereas the ribosome biogenesis pathway was upregulated in the progenitor subset (Fig. 4A, Fig. S6A).

The transition of SSCs into a differentiating state requires exposure to retinoic acid (RA), leading to the onset of KIT expression. The differentiation-primed population (Cluster IV) lacks both Gfra1 and Kit expression (Gfra1⁻/Kit⁻), suggesting that these spermatogonia are exiting self-renewal but uncommitted to differentiation. In support of this, GSEA analysis showed a strong repression of glial cell line-derived neurotrophic factor (GDNF)-regulated genes in Gfra1⁻/Kit⁻ cells (Fig. 4B) (Oatley et al., 2006). Along with this, previously identified pathways involved in the maintenance of SSCs, such as the Ras signaling pathway, cell-cell communication and the gonadotropin-releasing hormone pathway, are also inhibited in Gfra1⁻/Kit⁻ cells (Fig. S7, Table S7) (Lee et al., 2009). Moreover, Gfra1⁻/Kit⁻ cells express an increased level of differentiation regulator Rarg, indicating that they begin to acquire responsiveness to RA (Fig. S4F, Table S7) (Ikami et al., 2015). To understand better the molecular signatures associated with the priming state, we searched for genes with increased expression in Gfra1⁻/Kit⁻ cells. We identified several transcriptional regulators, such as Cfp1 (Cxxc1), Ctcf, Ehmt2 (also known as G9a), Tdrd5 and Trib3 (Fig. 4C, Table S7).

An advantage of sequencing single cells is that it allows reconstruction of the gene expression dynamics along the pseudotime. We applied Monocle analysis and identified 1223 genes with significant variation in expression levels along the differentiation trajectory (P<0.05). These genes are clustered into four groups according to their expression trend (Fig. 4D). During the first stage in the trajectory, the expression of known self-renewal genes, such as Ccnd2 and Etv5, as well as genes involved in MAPK/Erk1 and VEGF signaling began to decrease (Group 1 genes). Concurrently, a series of genes highly enriched for cell cycle and cell division was activated (Group 2 genes). Second, the cells passed through an intermediate phase and transiently expressed an elevated level of Nanos3 and genes related to the DNA-repair pathway (Bre and Trp53) (Group 3 genes). Third, the priming phase featured complete repression of SSC self-renewal genes and elevation of differentiation genes, such as Sohlh1 (Group 4 genes). Genes involved in the mTOR pathway (Eif4ebp1, Pdcd4, Cycs and Prkab2) were also upregulated, which is in line with enhanced mTOR pathway activity in adult spermatogonia differentiation (Fig. 4E) (Hobbs et al., 2010).

To identify novel regulatory networks associated with SSC self-renewal and differentiation, we applied the weighted gene co-expression network analysis (WGCNA) algorithm on transcriptomes of spermatogonia excluding gonocytes (Cluster I cells) (Langfelder and Horvath, 2008). In the search for modules that correlate with our differentiation model, we identified the pink module, which was enriched in primitive stem cells and decreased across pseudotime (Fig. S6B, left). GO terms of this module are associated with lipoprotein particle receptor activity, cell differentiation and development (Table S8). This unbiased approach highlighted key regulators as central to this network, including Etv5, Ccnd2 and Gfra1, and GDNF signaling regulation (19 out of 107 genes are GDNF regulated), which confirmed the current knowledge of SSC self-renewal (He et al., 2009). The second module was mostly repressed in differentiation-primed cells and was enriched for genes involved in mitotic cell cycle and RNA binding (Fig. S6B, right, Table S8). Hub-gene-network analysis indicated Tex13 and Rbm35a (Esrp1) are the most highly connected genes, which might exert potential functional roles in spermatogonia differentiation.

We also identified putative markers or regulators that were specific to the earlier, potentially self-renewing spermatogonia (Cluster II) (Fig. S8A). We found CD87 (PLAUR), the uPA receptor (uPAR), particularly interesting as it has been implicated in hematopoietic cell migration and homing and has been used as a marker for a subset of hematopoietic stem/progenitor cells (HSPCs) (Tjwa et al., 2009). Immunostaining of CD87 protein confirmed that expression is restricted to a subset of Oct4-GFP⁺ cells (Fig. 5A). In addition, RNA-Seq analysis of FACS-sorted cells showed that CD87⁺ cells expressed a higher level of self-renewal genes (Id4 and Gfra1) and a reduced level of early differentiating genes (Crabp1, Akt1) (Fig. 5B,C, Table S9). Importantly, in the CD87⁺ subpopulation, genes involved in cell migration were more abundantly expressed (Mmp9, Hes1, Chd1, Lrp1 and Procr). We then explored whether CD87 participated in germ-cell migration using a matrix adhesion assay. A portion of sorted cells seeded on Matrigel formed cellular protrusions, which represents an essential step for cell migration (Ladoux and Nicolas, 2012). Interestingly, CD87⁺ cells displayed cellular protrusions more frequently than CD87⁻ cells, suggesting a greater migration potential (Fig. 5D, Fig. S8B). The pseudotime analysis revealed a potential interaction between the cell cycle genes (Group II genes, Fig. 4D,E) and differentiation. Therefore, we investigated whether CD87⁺ cells, as a proxy to primitive stem cells, differed from CD87⁻ cells with regard to cell cycle activity. FACS analysis showed a higher frequency of cells in the G2/M phase in CD87⁺ versus CD87⁻ cells, likely reflecting a faster proliferation rate (Fig. 5E). Indeed, we observed a significantly higher growth rate for germ cell cultures derived from CD87⁺ cells (Fig. 5F). Interestingly, all Oct4-GFP⁺ cells expressed CD87 at early stages and the proportion gradually decreased from PND4.5 to PND7.5, and only a very small fraction of cells expressed CD87 at PND13.5 (Fig. 5G). Collectively, these results indicate that CD87⁺ cells in neonatal testis mark a unique cell state with higher expression of self-renewal genes and greater migration potential.

Oct4 has been shown to be crucial for SSC maintenance as Oct4 is downregulated by RA-induced differentiation and knockdown of Oct4 results in differentiation (Dann et al., 2008). Surprisingly, differentiation-primed cells expressed a notably higher level of Oct4 and Rarg (Fig. 6A). We confirmed that Oct4-GFP fluorescence intensity positively correlates with endogenous levels of both Oct4 mRNA and protein in single cells (Figs S4B, S9A). We then sorted the Oct4-GFP⁺/KIT⁻ cells into fractions with low, medium and high GFP expression levels and examined their gene expression profile by bulk RNA-Seq. Notably, the gene expression signature of Oct4-GFP-high cells strongly correlates with the differentiation-primed population, which is characterized by decreased expression of SSC self-renewal genes (e.g. Gfra1, Id4 and Etv5) and increased expression of genes associated with differentiation [e.g. Ngn3 (Neurog3), Rarg and Sohlh1] as well as the aforementioned candidates upregulated in differentiation-primed cells (e.g. Ehmt2, Trib3 and Tdrd5) (Fig. 6B). We validated this observation by FACS and qRT-PCR analyses on cells from PND7.5 testes, which did not contain gonocytes, and found consistently higher expression of Rarg, Sohlh1 and Ngn3 in Oct4-GFP-high cells (Fig. 6C). We further associated Oct4-GFP intensity to GFRA1 protein level by FACS analysis and found that Oct4-GFP-high cells contained significantly fewer GFRA1⁺ cells (Fig. 6D, Fig. S9B,C).

Having established that Oct4-GFP-high cells represent differentiation-primed cells in neonatal testis, we investigated whether this population differed in their potential for RA responsiveness. To test this, primary cultures derived from PND5.5 Oct4-GFP⁺/KIT⁻ cells were treated with all-trans RA or vehicle. Control cultures without RA exposure showed signs of spontaneous differentiation and Oct4-GFP-high cells contained a higher proportion of KIT⁺ cells (Fig. 6E). RA treatment did not alter the intensity of Oct4-GFP or the proportions of cell subsets with different Oct4-GFP intensities. This suggested that Oct4-GFP level was stable in this experimental setup and that cells of the same intensity in control or treatment groups could be directly compared. Indeed, there were significantly more Oct4-GFP-high cells that became KIT⁺ compared with Oct4-GFP-low and Oct4-GFP-mid cells (Fig. 6E).

Collectively, these data show that Oct4-GFP-high cells display molecular and functional properties associated with cells primed to differentiation.

To enable a more systematic analysis of the gene expression program accompanied by the Oct4-GFP intensity, we conducted Fuzzy clustering using all the differentially expressed genes among the three groups (Fig. 7A) (Kumar and Futschik, 2007). Interestingly, over 90% of the differentially expressed genes either show a gradually decreasing (Clusters 5 and 6) or increasing trend (Clusters 1, 3 and 4) from Oct4-GFP-low to Oct4-GFP-high transition (Fig. 7A, Table S10). Notably, Clusters 5 and 6 contain self-renewal genes, such as Gfra1, Fgfr1, Ret and Ccnd2. Further GO analyses revealed over-representation of the MAPK cascade (e.g. Dusp4, Dusp5, Nefl, Psme3) (Fig. 7B). These genes are co-enriched in several pathways related to SSC self-renewal, such as RET signaling and signaling to RAS and ERKs (Table S10). Surprisingly, stronger Oct4-GFP intensity also correlated with higher expression of mesenchymal markers (e.g. Vim, Dab2 and Emp3) along with genes significantly enriched in extracellular matrix organization, cell migration and cellular response to growth factor stimulus (Figs 6B and 7A,B). These results indicate that SSCs may modify their cellular adhesion property when they are primed for differentiation. Taken together, these results reveal a gradient of Oct4 expression coupled with self-renewal and differentiation transcriptional programs.

We hypothesized that higher Oct4 expression is important for the SSC differentiating potential and intended to reduce Oct4 expression closer to the ‘ground level’. Mild knockdown of Oct4 using siRNA showed that there is a significant reduction but not complete loss of Oct4 mRNA and protein level (Fig. S9D). The expression of key stem cell markers (i.e. Gfra1 and Id4) remained unchanged, while 590 differentially transcribed genes (FDR<0.05, fold change>1.5, Fig. 7C) were identified. Interestingly, the majority of differentially expressed genes were upregulated (n=500, 85%). The top genes included Etv4, a downstream factor of GDNF, and Egr4, which is important for male fertility (Fig. S10) (Lu et al., 2009; Tourtellotte et al., 1999). GO analysis clearly revealed significant enrichment for focal adhesion, which may be important for stem cell-niche interaction (e.g. Rhou, Ctnnb1 and Itgb1) (Table S11, Fig. S10). Of particular interest is Itgb1 (β1-integrin), which has been implicated as an essential adhesion receptor required for both gonocyte migration and SSC homing (Kanatsu-Shinohara et al., 2008; Tres and Kierszenbaum, 2005). The upregulated genes are also highly enriched in FGF signaling pathway components (e.g. Fgfr1, Etv4, Mapk1, Ctnnb1 and Fst) and MAPK protein kinases (e.g. Mapk1, Mapk12, Wee1 and Map3k2), which are downstream regulators of GDNF and FGF signaling (Fig. 7D) (Hasegawa et al., 2013). Only 90 downregulated genes were identified as a consequence of Oct4 knockdown (Fig. 7C). Among these, Pten is a direct target of OCT4 in cancer cells and in SSC maintenance (Huang et al., 2011; Tang et al., 2015). Dll3 is a ligand of the Notch signaling pathway (Fig. S10). To characterize further the functional consequence of the altered gene expression, we analyzed the ability of cultured germ cells to respond to short-pulse RA-induced differentiation with decreased Oct4 expression. Knockdown of Oct4 led to a 30% reduction of RA-induced Kit mRNA expression (Fig. 7E). These results suggest Oct4 may contribute to the differentiation-primed state of SSCs through transcriptional regulation.

Unraveling the heterogeneity of SSC by high-throughput single cell analysis has recently drawn attention. Studies on neonatal germ cells using single cell qRT-PCR analysis or Id4-GFP⁺ cells using single cell RNA-Seq has provided a better understanding of their heterogeneity. However, these studies focused on a panel of genes or a restricted population (Id4-GFP⁺) (Hermann et al., 2018, 2015; Song et al., 2016). Here, we extended these analyses to whole-transcriptome profiling on a broader undifferentiated germ cell population in PND5.5 testes.

To aid in the interpretation of these data, we first discuss challenges and limitations. Comprehensive characterization in an unbiased manner requires isolation of all subtypes during development. The Oct4-GFP⁺ population harbors GFRA1⁺ and PLZF⁺ cells, making it a good surrogate for undifferentiated germ cells (Fig. 1). However, we cannot rule out the possibility that rare or unknown subpopulations (i.e. Oct4-GFP⁻, or GFRA1⁻ or PLZF⁻) may not fit into the current classification scheme and yet contribute to the establishment of spermatogenesis. It is also possible that rare subpopulations within the Oct4-GFP⁺ compartment may not have been captured in our 71 samples during profiling. In this sense, our profiling might not be completely exhaustive and the existence of other subtypes remains a possibility. Nevertheless, our study presents an outlook on the heterogeneity of the main germ cell populations in neonatal testis.

Our study provides a global view of the gene expression signature in transitional gonocytes. First, we discovered several molecules that might be involved in gonocyte migration. Second, the gonocyte and spermatogonia population is distinguished by global transcription activity and change of cell cycle status. Third, our results showed that the expression of apoptosis-related genes is actively regulated during GST. Particularly, transitional gonocytes seem to have a lower threshold for initiating apoptosis (elevated death receptor pathway) and increased susceptibility to DNA damage-induced apoptosis (reduced DNA-repair pathway). This raises the possibility that transitional gonocytes possess poised apoptotic machinery, which renders them more vulnerable to apoptotic cell death. Further studies using our resources will help to shed light on the mechanistic nature of such transcriptional regulation in GST.

In corroboration of a recent ‘hierarchical model’ of adult SSCs, we showed that neonatal primitive stem cell and progenitor populations are distinguished by expression gradients of SSC maintenance genes (Id4, Gfra1, Etv5 and Tspan8) and differentiation-promoting genes (Sohlh1 and Rarg) (de Rooij, 2017). The ‘ultimate’ stem cells proposed in this model might reside in the primitive stem cell population we identified. CD87 is a glycosylphosphatidylinositol-anchored protein and acts as the receptor of urokinase-type plasminogen activator (uPA), which induces cell adhesion and migration (Ploug et al., 1991). We showed that CD87⁺ cells are most abundant within the first week after birth and have greater potential for migration in vitro. Interestingly, CD87 is also a ligand for several integrins including ITGB1 (β1-integrin). Thus, it is tempting to speculate that CD87 may play an important role in the initial establishment of the SSC pool, potentially through the regulation of migration.

Our results also provided an opportunity to reflect on whether the regulations of neonatal and adult spermatogonia are distinct. The stem cell-progenitor-differentiation-primed spermatogonia transition is reminiscent of the adult differentiation trajectory reported recently (Green et al., 2018; Hermann et al., 2018). Interestingly, Cluster I cells (including stem cells) expressed a higher level of a set of genes that overlapped with the FGF target genes (e.g. Spry4, Dusp6, Lhx1 and Etv4) as reported recently (Kitadate et al., 2019). Kitadate et al. demonstrated the role of FGFs in the cell fate decision of adult SSCs to either self-renew or differentiate. Echoing this, we found that syndecan genes (Sdc1 and Sdc2, regulators of FGF signaling) are more abundantly expressed in Cluster I primitive stem cells (Fig. S6A), suggesting that Cluster I cells may be better furnished with the reception machinery for FGFs. Therefore, it is possible that the mechanism of mitogen (FGF) competition in adult SSCs is also shared by neonatal spermatogonia.

At the same time, this high-resolution data set provides a more comprehensive characterization of transcriptome-wide features and differences between stem cells and progenitors. Our analysis showed that genes related to lipid metabolism are highly enriched in the stem cell pool, which extends previous findings on the role of fatty acid metabolism in stem cell maintenance (Ito and Suda, 2014). Our results also indicated that ribosome biogenesis increased in progenitors compared with stem cell spermatogonia, implying an enhanced need for global protein synthesis during differentiation. Our results further demonstrated that the SSC-to-progenitor transition is intertwined with cell cycle states.

Although the action of RA signaling on the neonatal spermatogonia differentiation process is known, the transcriptional processes that regulate the first steps in differentiation are not well defined. Surprisingly, our study uncovered a differentiation-primed population with a higher level of Rarg. Higher Rarg expression may enhance the sensitivity of these cells to RA exposure, which closely resembles the recently described differentiation model in adult spermatogenesis (Ikami et al., 2015). These observations revealed the earliest event of SSC differentiation in which cells turn off the self-renewal program before committing to differentiation, which can only be resolved by single cell analysis.

Differentiation-primed cells notably display a surprisingly higher level of Oct4 (Fig. 6A, Fig. S4D,F). The association between Oct4 and SSC differentiation has been overlooked but it has been recently demonstrated that Oct4-GFP⁻ and Oct4-GFP⁺ undifferentiated spermatogonia population in adult testis represent stem cells and progenitors, respectively (La et al., 2018). Our study also showed that higher Oct4-GFP was associated with a more differentiated state and that the gradient of Oct4 induction was paralleled by the gradual loss of self-renewal genes and gain of differentiation-promoting genes. Functionally, we showed that Oct4-GFP-high cells displayed stronger RA responsiveness, which shared the functional property of Id4-GFP-dim progenitors as reported previously (Lord et al., 2018).

A prominent unanswered question is whether and how OCT4 contributes to the enhanced RA response of SSCs. In our mild knockdown experiment, the expression of key self-renewal markers (i.e. Gfra1 and Id4) does not appear to be affected, suggesting that there was little change in the relative abundance of cell subpopulations. It is also unlikely that OCT4 directly controls the level of Rarg and positively regulates the RA pathway because we did not identify a significant change in the expression of differentiation-related genes (i.e. Rarg, Kit and Stra8). Instead, OCT4 appears to modulate downstream signaling of the SSC self-renewal pathway, such as the MAPK pathway. Our results also raised the possibility that OCT4might modulate the stem cell-niche interaction in SSCs. Oct4-GFP-high cells express a higher level of mesenchymal markers and enrichment of genes associated with cell migration. In addition, many OCT4-regulated genes were associated with focal adhesion (e.g. Itgb1). Thus, it is possible that OCT4 could modulate the exposure of SSCs to niche-secreted short-range signals and cell contact-dependent signals. Taken together, we speculate that different levels of OCT4 expression may act as a molecular rheostat to regulate self-renewal versus differentiation of SSCs in neonatal testis. Elevation of Oct4 in the earliest stage of SSC differentiation leads to the repression of several downstream components of the GDNF/FGF pathway, as well as the exposure of SSCs to niche signals, facilitating the response to the differentiation cue. Our work thus provides insights into the important issue of the coupling between cell adhesion, stem cell maintenance and differentiation.

In conclusion, we envisage that knowledge of transcriptional dynamics and cell fate decisions during early germ cell establishment should provide hints to aid further dissection of the heterogeneity in adult human and mouse spermatogonia. It will also help the study of SSCs in the context of aging and disease.

All animal procedures were performed according to protocols approved by the Animal Experimentation Ethics Committee (AEEC) of the Chinese University of Hong Kong and following the Animals (Control of Experiments) Ordinance (Cap. 340) licensed from the Hong Kong Government Department of Health. Oct4-GFP transgenic mice [B6; CBA-Tg(Pou5f1-EGFP)2Mnn/J, Stock no.: 004654] were acquired from The Jackson Laboratory and maintained in the CUHK Laboratory Animal Services Centre.

Neonatal testes were harvested and fixed in 4% paraformaldehyde (PFA) solution for 2 h followed by dehydration in 30% sucrose. Then, the testes were embedded in OCT (Sakura) and snap-frozen. Cryosections of 10 μm thickness were collected on positively charged slides (Thermo Fisher Scientific). Sections were blocked with PBS supplemented with 0.2% Triton X-100 and 10% normal donkey serum (Jackson ImmunoResearch) for 1 h at room temperature, followed by overnight incubation with antibodies against KIT, GFRA1, PLZF, APC-Ki-67 and CD87 at 4°C. After washing in PBS, sections were incubated with the corresponding secondary antibody (except for Ki-67 staining). Antibodies used are listed in Table S12. Finally, tissue slides were washed with PBS and mounted with mounting medium (Thermo Fisher Scientific). A minimum of three different animals were tested and negative controls lacking primary antibody were used. Fluorescent microscopy and brightfield images were acquired using a Leica SP2 confocal microscope. All images were captured using Leica Application Suite Software. For quantification of cells with overlapping expression of markers, different localization or Ki-67 expression, each staining was performed in triplicate and images were taken for at least three sections for each testis. Cells were counted in ten randomly selected testis cords in each testis section.

Undifferentiated germ cells from testes of Oct4-GFP transgenic mice at PND5.5 or PND7.5 were purified using the method described previously (Garcia and Hofmann, 2012). For the Fluidigm C1 single cell capture experiment, testes of five pups at PND5.5 from one litter were pooled before FACS sorting. In the Biomark qRT-PCR experiment, testes of four pups at PND5.5 from another litter were pooled before FACS sorting. Isolated seminiferous tubules were digested with 1 mg/ml type 2 collagenase (Worthington), 1 mg/ml hyaluronidase (Sigma-Aldrich) and 5 μg/ml DNase I (Sigma-Aldrich) at 37°C for 20 min with occasional shaking. The suspension was passed through a 40-μm strainer cap (BD Falcon) to yield a uniform single cell suspension. After incubation with antibodies against KIT or/and CD87 at 4°C for 30 min, cell fractions were collected with a FACS Aria II cell sorter (BD Biosciences). In order to obtain pure populations, sorting gates were set within the boundaries defined by to ‘fluorescence minus one’ (FMO) controls. Indirect labeling was used for GFRA1 and OCT4 analysis. Cells were stained for 30 min on ice with APC-H7-conjugated anti-KIT antibody and washed with PBS. Cells were then fixed in 2% PFA, permeabilized in Perm buffer (0.5% Tween 20, 2% fetal bovine serum in PBS). Washed cells were incubated with anti-GFRA1 or anti-OCT4 antibodies at 4°C for 30 min. After extensive washing, cells were incubated with Alexa Fluor 647 secondary antibody at 4°C for 30 min. Positive GFRA1 or OCT4 antibody labeling was determined by comparison with staining with Alexa Fluor 647 antibody only. Age- and strain-matched non-reporter mice, i.e. C57BL/6 mice, were used for cell purification as a negative control devoid of GFP expression for the cell-sorting experiment. The antibodies used in this study can be found in Table S12.

Quantitative real-time PCR (qRT-PCR) was performed using the Fast SYBR Green Master Mix (Applied Biosystems) on a ViiA 7 Real-Time PCR System (Applied Biosystems) according to the manufacturer's instructions. qRT-PCR was carried out with at least three biological replicates per group, and each biological replicate was performed in duplicate or triplicate. Primers used are listed in Table S13.

Single cells were captured on a medium-sized (10-17 μm cell diameter) microfluidic RNA-Seq chip (Fluidigm) using the Fluidigm C1 Auto Prep system. Cells were loaded onto the chip at a concentration of 300-500 cells/μl and stained for viability (LIVE/DEAD cell viability assay; Molecular Probes, Life Technologies). The Fluidigm C1 IFC Chip after cell capture was inspected under the microscope and dead cells and doublets captured sites were eliminated. cDNAs were prepared using the SMARTer Ultra Low RNA Kit for Illumina (Clontech). Libraries were generated using the Nextera XT (Illumina) protocol and were sequenced as 50-bp single-end reads on Illumina HiSeq 2000.

Raw reads were pre-processed with trim galore for trimming Illumina adaptor sequences and mapped to mouse mm9 genome. We estimated the expression level for all RefSeq genes using Salmon (Patro et al., 2017) and selected the 500 top-ranked PCA genes for clustering. Pseudotime ordering was performed with Monocle (v1.2.0) as previously described (Trapnell et al., 2014). Differential gene expression analysis was performed with DESeq2 method (Love et al., 2014). R package ‘WGCNA’ was used to construct the weighted gene co-expression network. Further details can be found in supplementary Materials and Methods.

Oct4-GFP⁺/KIT⁻ cells were FACS-sorted from PND5.5 testes based on different GFP intensity levels. Germ cell culture in the Oct4 knockdown experiments were FACS-sorted to remove MEF feeder cells. Oct4-GFP⁺/KIT⁻/CD87⁺ and Oct4-GFP⁺/KIT⁻/CD87⁻ cells were also sorted by FACS. Total RNA was isolated from sorted cells using the AllPrep DNA/RNA Mini kit (Qiagen) and subjected to low input RNA-Seq library preparation. Further details of methods and bioinformatic analysis can be found in supplementary Materials and Methods.

Single cell gene expression analysis was performed using the BioMark 48.48 Dynamic Array platform (Fluidigm) and custom-designed gene-specific primer pairs (Table S13). We sorted Oct4-GFP⁺/KIT⁻ single cells using FACS into 96-well PCR plates containing 5 μl of CellsDirect reaction mix and these cells were immediately stored at −80°C. Upon thawing, a mix containing 2.5 μl gene-specific primers, 1.2 μl CellsDirect RT/Taq mix and 0.3 μl TE buffer were added to each well. RT-PCR pre-amplification cycling conditions were: 50°C, 15 min; 95°C, 2 min; 22× (95°C, 15 s; 60°C, 4 min). Samples were diluted 1:5 in TE buffer and a 2.25 μl aliquot of each cDNA was mixed with 2.5 μl of 2× SsoFast EvaGreen Supermix with Low ROX (Bio-Rad) and 0.25 μl of 20× DNA Binding Dye Sample Loading Reagent (Fluidigm). The sample mix and custom-designed primer pairs were loaded separately into wells of 48.48 Gene Expression Dynamic Arrays (Fluidigm) in the presence of appropriate loading reagents. The arrays were read using a Biomark analysis system (Fluidigm). Data were analyzed using Fluidigm Real-Time PCR Analysis software (v3.0.2). Cts were recovered from BioMark. The quality threshold was set to 0.65 and a linear derivative as baseline correction was used. The limit of detection was set to a Ct of 28; Cts for qRT-PCR reactions that failed the quality threshold were converted to 28. Cts were converted to log2 expression by subtracting Cts from 28. Expression matrix of the genes analyzed can be found in Table S4.

For RA-induced differentiation, 5 mM all-trans-RA (Sigma-Aldrich) stock in DMSO was diluted to 1 μM in germ cell culture medium before applying to cells; vehicle (DMSO) alone was added in the control group. For Fig. 6E, the cells were collected for FACS analysis after 16 h. For Fig. 7E, the cells were collected by FACS 4 h after treatment and sorted cells were subjected to RNA extraction and qRT-PCR analysis.

Cells were stained for 30 min on ice with APC-H7-conjugated anti-KIT and APC-conjugated anti-CD87 antibodies and washed with PBS. Cells were then fixed in 2% PFA and permeabilized in Perm buffer. Cells were washed and stained for 45 min on ice with APC-conjugated anti-Ki-67 antibody. The cells were washed again and stained with Hoechst 33342 (Sigma-Aldrich) at 5 μg/ml final concentration in complete medium for 30 min at 37°C. The stained cells were re-suspended in cold PBS and run on the FACS analyzer. Data were analyzed using FlowJo software (Tree Star).

Culture plates were coated with Matrigel (BD Biosciences) by the following procedure. A Matrigel bottle was thawed on ice in a 4°C refrigerator overnight until the Matrigel liquified. Matrigel was divided into 300 μl aliquots and stored at −20°C until use. Working Matrigel solution was prepared by diluting 300 μl of Matrigel with 29 ml of DMEM/F12 medium (Invitrogen) and thorough mixing. This solution was added to 24-well plates (0.5 ml per well) to cover the whole surface of the wells. The plates were allowed to sit for 1 h at room temperature. Excess Matrigel solution was then removed, and the plates were washed once with DMEM/F12. An equal number of the FACS-sorted Oct4-GFP⁺/KIT⁻/CD87⁺ or Oct4-GFP⁺/KIT⁻/CD87⁻ cells were plated on three Matrigel-coated wells. After 16 h incubation, images were taken from four random fields of each well and the proportion of cells with protrusion was calculated as the number of cells with protrusions divided by the total number of cells in all 12 images in each biological replicate.

Small interference RNAs (siRNAs) against Oct4 were synthesized by GenePharma (Shanghai, China). The sequences of siRNAs are: sense, 5′-AAGGAUGUGGUUCGAGUAUGG-3′; antisense, 5′-CCATACTCGAACCACATCCTT-3′. The negative control sequence has no homology with known mouse genes. The siRNA was transfected into SSC culture cells using Lipofectamine 2000. Briefly, a 100 μl OptiMEM mixture containing 2.5 μl Lipofectamine and 10 pmol siRNA was prepared according to the manufacturer's instructions. The mixture was added to 400 μl of SSC culture medium containing 1.5×10⁵ trypsin-dissociated SSCs and incubated at 37°C for 30 min. The SSC-Lipofectamine-siRNA mixture was then added to a 24-well of pre-seeded MEF and incubated for 6-8 h after which the medium was replaced with normal SSC culture medium. The SSCs were allowed to grow for another 48 h and subjected to direct gene expression analysis or RA-induced differentiation.

Experimental data are shown as mean±s.e.m. All statistical analyses were either conducted with JMP or as specified in relevant sections.