Differential RA responsiveness directs formation of functionally distinct spermatogonial populations at the initiation of spermatogenesis in the mouse

Morphogens induce cellular differentiation based on their concentrations across a tissue, the competence of cells within that tissue to respond to morphogens, or a combination of both (Turing, 1952; Waddington, 1940). Retinoic acid (RA) is a classic morphogen that, as a small lipophilic molecule, predominantly signals in a paracrine fashion in localized concentration gradients. RA signaling is then initiated in target cells following binding to a monomeric RA receptor (RAR), of which there are three identified subtypes (RARA, RARB and RARG). This is followed by heterodimerization with a monomer of one of three identified retinoid X receptor subtypes (RXRA, RXRB and RXRG) (Chambon, 1996; Mark et al., 2015). In the mammalian testis, RA acts within germ cells and somatic cells to regulate gene expression at both the transcriptional and post-transcriptional levels (Chambon, 1996; Mark et al., 2015; Busada et al., 2015b; Busada and Geyer, 2016; Busada et al., 2014; Chung et al., 2004; Chung et al., 2005; Chung et al., 2009; Gely-Pernot et al., 2012; Gely-Pernot et al., 2015; Hasegawa and Saga, 2012; Zhou et al., 2008a). Ultimately, RA signaling has pleiotropic roles during development (reviewed by Niederreither and Dolle, 2008; Gudas and Wagner, 2011; Janesick et al., 2015), and multiple studies have unequivocally demonstrated that RA is required for male fertility (reviewed by Busada and Geyer, 2016).

RA levels are generally tightly controlled by the combined action of enzymes regulating its synthesis (retinol dehydrogenases and aldehyde dehydrogenases) and catabolism (CYP26 enzymes) (reviewed by Mark et al., 2015; Griswold, 2016). RA is synthesized in target tissues by converting circulating or stored retinol via a two-step reaction, and the proposed source of RA in the developing testis is Sertoli cells; Sertoli cell expression of the synthesizing enzymes aldehyde dehydrogenase family 1 (ALDH1A1-3) and retinol dehydrogenase 10 (RDH10) is essential for initiating the first wave of spermatogenesis (Kent et al., 2016; Raverdeau et al., 2012; Tong et al., 2013). From the site of biosynthesis, RA then diffuses and can act either through binding its cognate receptors or by being converted into inert polar metabolites (Niederreither et al., 2002) by the intracellular cytochrome P450 enzymes CYP26A1, CYP26B1 and CYP26C1 (reviewed by Nelson et al., 2013). In terms of substrate specificity, CYP26C1 has greater specificity for 9-cis RA, whereas CYP26A1 and CYP26B1 predominantly catabolize all-trans RA (hereafter referred to as RA) (Thatcher and Isoherranen, 2009). Cyp26a1 and Cyp26b1-knockout (KO) mice exhibit embryonic or perinatal lethality, whereas Cyp26c1-KO mice are viable and fertile, with no apparent malformations (Abu-Abed et al., 2001; MacLean et al., 2001; Taimi et al., 2004; Uehara et al., 2007; Yashiro et al., 2004).

Spermatogenesis begins in the neonatal mouse testis with the segregation, starting at postnatal days (P)3-4, of precursor prospermatogonia into either undifferentiated or differentiating spermatogonial populations (Busada and Geyer, 2016; Griswold, 2016; Geyer, 2017; Niedenberger et al., 2015). A subpopulation of nascent undifferentiated spermatogonia are spermatogonial stem cells (SSCs), which retain perpetual self-renewal capacity and lack the capacity to directly differentiate. The remaining undifferentiated progenitor spermatogonia lack self-renewal capacity and undergo transit-amplifying cell divisions before responding to RA and dividing a finite number of additional times before entering meiosis. This process is continuously maintained during steady-state spermatogenesis in the adult testis. For over 60 years, research has shown that most rodent undifferentiated spermatogonia are progenitors that differentiate at stage VIII of the seminiferous epithelial cycle (Clermont and Leblond, 1953). This coincides with a transient increase in RA levels (Hogarth et al., 2014), apparently produced by spermatocytes (Griswold, 2016; Raverdeau et al., 2012). SSCs within those same tubules somehow avoid RA-induced differentiation, whereas their nonself-renewing/transit-amplifying progenitor progeny will expand in subsequent stages without responding to RA, and then differentiate 8.6 days later at the next stage VIII (Clermont, 1972). Therefore, identifying the underlying mechanisms that control response to RA is central to understanding how spermatogonial fate is regulated.

It is currently unclear how the spermatogonial response to RA is regulated at the molecular level. We envision two general scenarios, which are not mutually exclusive, for RA responsiveness: (1) all spermatogonia are exposed to RA, but only some can respond (intrinsic or germ cell-autonomous control); or (2) all spermatogonia are equivalently primed to respond, but exposure to RA is restricted by a specialized niche microenvironment (extrinsic control). Most previous studies investigating the RA response of neonatal germ cells have utilized genetic models, in which regulators of RA binding and sequestration, metabolism (synthesis versus catabolism) or reception have been ablated using whole-body and conditional KO approaches (reviewed by Mark et al., 2015; Busada and Geyer, 2016). In most cases, loss of single or even multiple proteins (e.g. the cellular RA-binding proteins CRABP1/2 or RARA/B/G) predicted to have crucial roles in RA metabolism during spermatogenesis resulted in spermatogenic defects that are mild in comparison to the gold-standard, vitamin A deficiency, in which spermatogonial differentiation is completely halted (Gely-Pernot et al., 2012, 2015; Hogarth et al., 2015; Li et al., 2009). It is evident from these studies that multiple proteins, including unidentified contributors, manage the normal spermatogonial RA response.

Here, we used the neonatal mouse as a model system to determine whether RA responsiveness is regulated to form the heterogeneous populations of stem, progenitor and differentiating spermatogonia. Our results, gathered from both in vitro and in vivo pharmacological approaches, show that, whereas most undifferentiated progenitor spermatogonia are poised to respond to RA, CYP26-mediated catabolism prevents the response of this population. Therefore, CYP26 enzyme activity contributes to patterning the germ cell population in the neonatal testis by creating one or more functional RA concentration gradients. We also report that another subset of undifferentiated spermatogonia (the presumptive foundational SSCs) fail to respond to high RA levels either in vitro or in vivo. Lastly, we report for the first time that these functionally distinct populations exist when the germline is quiescent, at the prospermatogonial stage before initiation of the first wave of spermatogenesis. Consequently, we propose a model whereby the distinct spermatogonial populations are determined before the initiation of spermatogenesis with differential intrinsic capacity to respond to RA.

It is unknown how the initial populations of mammalian stem, progenitor and differentiating spermatogonial populations are formed at the initiation of spermatogenesis. We and others have shown that administration of exogenous RA initiates protein expression of both stimulated by RA gene 8 (STRA8) and kit proto-oncogene receptor tyrosine kinase (KIT) in most spermatogonia in the developing testis (Busada et al., 2015b; Busada et al., 2014; Zhou et al., 2008a; Niedenberger et al., 2015; Zhou et al., 2008b). To confirm and quantify this, we assessed the extent of this differential responsiveness to RA in wild-type (WT) mice in vivo. We used mice at P5-6, ages when spermatogonia (stem, progenitor and differentiating) are the only germ cells present (Bellve et al., 1977; Drumond et al., 2011; Kluin et al., 1982; Kluin et al., 1984). In vehicle (dimethyl sulfoxide; DMSO)-treated controls, only 20% of spermatogonia were STRA8+ (Fig. 1A,C), similar to our previous observations in untreated WT testes (Busada et al., 2015b; Niedenberger et al., 2015; Hermann et al., 2015). We defined these as ‘RA-competent/RA-responsive’ spermatogonia. In mice treated with a bolus of 100 μg RA, 85% of spermatogonia became STRA8+, revealing that an additional 65% of spermatogonia were capable of responding to RA, and we defined these as ‘RA-competent/RA-nonresponsive’ spermatogonia because they required elevated RA levels to become STRA8+ (Fig. 1B,C, P<0.05). This >4-fold increase in the number of STRA8+ spermatogonia following RA injection revealed that approximately two-thirds of spermatogonia are capable of responding to RA but, for whatever reason, do not normally do so in vivo. Therefore, a mechanism must exist to prevent these RA-competent/RA nonresponsive germ cells from responding to RA in the developing testis until the appropriate time. Moreover, 15% of spermatogonia did not respond to RA, which we defined as ‘RA-incompetent/RA-nonresponsive’, which could represent the nascent SSC pool that must remain insensitive to RA to ensure ongoing spermatogenesis (Fig. 1B,C, P<0.05).

We next addressed whether unequal availability of RA to germ cells following subcutaneous injection at a site distant from the testis in vivo underlies the differential ability of spermatogonia to respond to RA. This could be mediated by several factors, including differential RA concentrations throughout the testis or the effect of distinct spermatogonial niche microenvironments, as suggested recently (Lord et al., 2018). We tested the intrinsic, or germ cell-autonomous, responsiveness of neonatal germ cells in vitro by culturing freshly isolated testicular single cell suspensions ±RA for 24 h to maximize the uniformity of RA exposure. Only 20% of P5 spermatogonia were STRA8+ (RA-competent/RA-responsive) in vehicle-treated controls (Fig. 1D,F). Twenty-four hours after administration of RA to the culture medium, 93% of spermatogonia became STRA8+ (Fig. 1E,F), demonstrating that 93% were RA-competent/RA-nonresponsive and 7% were ‘RA-incompetent/RA-nonresponsive’.

We then examined the intrinsic RA responsiveness of isolated spermatogonia maintained in vitro without somatic cells. We isolated Enhanced Green Fluorescent Protein (EGFP)+ spermatogonia from inhibitor of DNA binding 4 (Id4)-eGfp mice, in which epifluorescence marks spermatogonia with differing regenerative capacity following transplantation (Hermann et al., 2015; Chan et al., 2014; Helsel et al., 2017; Mutoji et al., 2016). In this model, ID4-EGFP^bright spermatogonia are highly enriched for regenerative SSCs, whereas ID4-EGFP^dim and ID4-EGFP– spermatogonia are depleted of regenerative capacity and represent progenitor and differentiating spermatogonia, respectively. In vehicle-treated controls, few ID4-EGFP^bright spermatogonia were STRA8+ in the absence of RA (3% at 6 h and 0% at 12 h, Fig. S1A,C,E). In response to RA, there was a ∼3-fold increase in the percentage of STRA8+ (RA-competent/RA-nonresponsive) ID4-EGFP^bright spermatogonia (22% at 6 h and 10% at 12 h, Fig. S1B,D,E). These results corroborate previous in vivo results (Figs 1, 2), revealing that a limited proportion of postnatal ID4-EGFP^bright spermatogonia are RA competent (with and without RA response). A similar inability of ID4-EGFP^bright spermatogonia to respond to RA in vitro and in vivo revealed that differential RA responsiveness is an intrinsic feature of developing male germ cells, with minimal reliance upon an in vivo niche microenvironment. Importantly, these results also strongly support the existence of three functionally distinct spermatogonial groups (RA-competent/RA-responsive, RA-competent/RA-nonresponsive and RA-incompetent/RA-nonresponsive).

We predicted that most neonatal spermatogonia that can respond to RA but do not normally do so in vivo represent undifferentiated progenitors that are poised to differentiate. At P5, only 11% of ID4-EGFP^bright spermatogonia were STRA8+ (RA-competent/RA-responsive) (Fig. 2A,E), which suggests that this small subset of ID4-EGFP^bright spermatogonia are not SSCs, but early transit-amplifying progenitors. We next determined what percentage of the ID4-EGFP^bright spermatogonial population would respond to RA (by becoming STRA8+) in response to exogenous RA. By 6 h after RA injection, 33% of the ID4-EGFP^bright population had become STRA8+ (RA-competent/RA-nonresponsive) (Fig. 2B,E), a 3-fold increase over the 11% of ID4-EGFP^bright spermatogonia that became STRA8+ in vehicle-treated mice over this interval (P=0.009). Twelve hours after RA injection, only 8% of ID4-EGFP^bright spermatogonia were STRA8+, and this dropped to 0% at 24 h (Fig. 2C-E). We previously showed that administration of RA does not increase germ cell apoptosis (Busada et al., 2014; Niedenberger et al., 2015); therefore, these results revealed that approximately one-third of neonatal ID4-EGFP^bright spermatogonia were RA-competent/RA-nonresponsive, which could mean that they subsequently differentiate and lose EGFP expression. Conversely, 67% of the ID4-EGFP^bright population were ‘RA-incompetent/RA nonresponsive’ (Fig. 2E), which is a requisite functional feature of SSCs in the developing and adult testis; therefore, these cells probably represent the bona fide SSC population in vivo.

RARG is the predominant RAR subtype expressed in spermatogonia (Gely-Pernot et al., 2012; Vernet et al., 2006), and was recently shown to be sufficient to drive precocious spermatogonial differentiation in the mixed SSC/progenitor population among adult glial cell-derived neurotrophic factor family receptor alpha 1 (GFRA1)+ undifferentiated spermatogonia (Ikami et al., 2015). The clear implication is that SSCs remain undifferentiated because they lack RARG. However, a recent report challenged this concept by reporting that RARG was detectable in ID4-EGFP+ spermatogonia in the developing testis (Lord et al., 2018). As previously shown in the adult and developing testis (Helsel et al., 2017), the spermatogonial population labeled by the Id4-Egfp transgene is a mixture of those with the highest regenerative capacity (ID4-EGFP^bright; thus, enriched for SSCs) as well as progenitor spermatogonia (ID4-EGFP^dim). Therefore, we examined the association of RARG with ID4-EGFP intensity along with established markers of spermatogonial fate in vivo. We performed co-immunostaining for GFRA1+RARG and KIT+RARG in P6 Id4-eGfp testes and found that nearly all GFRA1+ spermatogonia were also ID4-EGFP+, and nearly all KIT+ spermatogonia were ID4-EGFP–, as expected (Fig. 3A,B). Almost none of the spermatogonia were ID4-EGFP^bright/RARG+ (7 out of 1235 cells counted, 0.6%). By contrast, 33% of spermatogonia were ID4-EGFP^dim/RARG+, and 10% were ID4-EGFP–/RARG+ (Fig. 3C). Therefore, we conclude that, in the developing testis, RARG status correlates with RA responsiveness in progenitor and differentiating spermatogonia.

We next wanted to determine the temporal origin of the differential RA responsiveness we linked to the cellular fate of newly formed spermatogonia. Mammalian spermatogonia arise from a quiescent pool of precursor prospermatogonia, which have long been presumed to be equipotent with regards to the types of spermatogonial fate they can assume, because they exhibit no known features of molecular or functional heterogeneity (Niedenberger et al., 2015; Drumond et al., 2011; Kluin and de Rooij, 1981). Identification of heterogeneity in the precursor prospermatogonial population would provide a significant advance by suggesting that spermatogonial fates are either predetermined or established before their formation.

We explored whether prospermatogonia at P1-2 also exhibit the differential RA responsiveness characteristic of newly formed spermatogonia at P5. Until ∼P2, the only germ cells present in the mouse testis are prospermatogonia, which have not been exposed to RA, because the first STRA8+/KIT+ differentiating spermatogonia do not appear until ∼P3-4 (Busada and Geyer, 2016; Zhou et al., 2008a; Niedenberger et al., 2015; Yoshida et al., 2006). We treated P1 pups with vehicle alone (DMSO) or 100 μg of RA and quantified the numbers of prospermatogonia able to respond to RA in vivo 24 h later. As expected, we observed no STRA8 staining among DEAD (Asp-Glu-Ala-Asp) box polypeptide 4 (DDX4)+ prospermatogonia in testes from mice treated with vehicle alone (Fig. 4A,C). By contrast, RA treatment induced STRA8 in 57% of prospermatogonia (Fig. 4B,C). Therefore, more than half of the precursor prospermatogonia population at P1 were RA-competent/RA-nonresponsive (Busada and Geyer, 2016; Zhou et al., 2008a; Niedenberger et al., 2015; Yoshida et al., 2006).

We then assessed whether the differential responsiveness of precursor prospermatogonia to RA is an intrinsic feature of germ cells or patterned by putative niche microenvironments within the testis cords. We cultured freshly isolated P1 testicular single cell suspensions with or without RA for 24 h to ensure uniform RA exposure among prospermatogonia. As expected, no STRA8+ cells were detected in vehicle-treated controls (Fig. 4D,F). In response to RA, 67% of prospermatogonia became STRA8+ (Fig. 4E,F), mimicking results gleaned in vivo (Fig. 4A-C).

We used single cell RNA-seq to profile the transcriptomes of individual testis cells and identify potential biological mechanisms driving the differential response to RA in vivo. We selected the P1 testis for examination, because the only germ cells present are precursor prospermatogonia that have not yet responded to RA (Fig. 4); ∼60% are able to respond to RA (RA-competent/RA-nonresponsive), whereas ∼40% are not able to respond to RA (RA-incompetent/RA-nonresponsive). Testis cells were isolated 12 h after P1 pups received injections of DMSO or RA, and were then subjected to 10x Genomics analysis (Zheng et al., 2017). Results from a total of 34,153 cells were highly consistent among duplicate cell preparations from independent litters treated with DMSO (15,417 cells) or RA (18,736 cells) (Fig. S2A). Unsupervised clustering of the aggregated single cell transcriptomes projected onto t-distributed stochastic neighbor embedding (t-SNE) plots revealed multiple heterogeneous cell clusters (Fig. 5A) that represented the complete complement of testicular somatic cells and prospermatogonia that we identified based on expression of key cell type-specific gene markers (Fig. 5B, Fig. S2B, Table S1). Prospermatogonia were identified in cell cluster 12 based on expression of Ddx4 (Fig. 5C) and other germ cell markers [e.g. deleted in azoospermia-like (Dazl), Id4, reproductive homeobox 10 (Rhox10) and spermatogenesis and oogenesis-specific basic helix-loop-helix 1 (Sohlh1)] (Fig. S2B) and limited detection of the somatic marker GATA-binding protein 4 (Gata4, Fig. 5C), along with an extensive panel of somatic cell type-specific genes (Fig. S2B). When examined in isolation, this cluster of unselected prospermatogonia exhibited further heterogeneity. This allowed resolution of five distinct Ddx4-expressing germ cell clusters (Fig. S2C), one of which was contaminated with Gata4-expressing somatic cells (Fig. S2C). Further interrogation of the remaining four prospermatogonial clusters revealed no cluster bias between cells arising from DMSO- or RA-treated testes (Fig. 5D), but identified 30 genes that were differentially expressed between cells from the two treatment groups (Fig. 5E, Table S1). This list included several genes shown previously to be RA responsive in P1 prospermatogonia in vitro (Zhou et al., 2008b) (Fig. 5F). As expected, these results revealed considerable heterogeneity in the response to RA (Fig. 5E,F).

To investigate the differential response to RA, we examined prospermatogonia from RA-treated testes in isolation and identified four cell clusters (Fig. 5G) characterized by panels of significantly differentially expressed genes (DEGs) (Fig. 5H). Surprisingly, these clusters did not clearly separate prospermatogonia based on RA response (Fig. 5I). To provide enhanced resolution of this distinction, we identified 213 putative RA-responsive prospermatogonia [Stra8+ or superoxide dismutase 3, extracellular (Sod3)+, 57% of all RA-treated prospermatogonia] and compared their transcriptomes to those of 159 prospermatogonia negative for both genes (Fig. 5J). In total, we identified 65 genes that were detected at significantly higher levels in RA-responsive prospermatogonia (Fig. 5K), including those involved in cell cycle regulation, oxidative phosphorylation and signaling pathways active in cancer (Fig. 5L). Ten genes were detected at higher levels in RA-nonresponsive prospermatogonia (Fig. 5K), and are involved in translational initiation {e.g. EIF2 signaling [ribophorin 1 (Rpn1) and ribosomal protein L6 (Rpl6)]}, EEF2-dependent translational elongation (diphthamide biosynthesis) and inhibition of gene expression by demethylation of histone H3 K4 [lysine (K)-specific demethylase 1B (Kdm1b)] (Fig. 5K,L).

These data also revealed considerable heterogeneity among prospermatogonia binned as RA-responsive and RA-nonresponsive (Fig. 5J,K), suggesting poor separation of biologically meaningful groups. Thus, we then performed unbiased dynamic cell trajectory analyses using all RA-treated prospermatogonia, which placed each germ cell in pseudotime order on the basis of their transcriptomes (Fig. 5M). This trajectory analysis yielded four discrete cell states along a developmental continuum that contained two predicted nodes with minor branching points (Fig. 5M,N). Most prospermatogonia were found in State 1 or 3, which existed at the beginning and end of pseudotime, respectively (Fig. 5M,N). Although levels of pan-prospermatogonial genes, such as Dazl and Ddx4, were essentially invariant across pseudotime, levels of Stra8 increased over pseudotime (Fig. 5O), reflecting the activation of the RA response. We found that mRNA levels encoding the RA catabolic enzyme CYP26B1 were highest at the beginning of pseudotime and declined as pseudotime proceeded (Fig. 5O). Indeed, for cells in which either Stra8 or Cyp26b1 was detected at levels of >1 count per cell, expression of these genes was weakly inversely correlated (Pearson's r=−0.19) in individual RA-treated prospermatogonia (Fig. 5O). Importantly, levels of Cyp26b1 were also variable among DMSO-treated P1.5 prospermatogonia (data not shown). Having established Cyp26b1 expression as a potential intrinsic regulator of the differential prospermatogonial RA response, we sought to identify transcriptome-wide gene expression patterns that could elucidate additional contributors to this process. We identified six clusters of genes that exhibited differential expression among RA-treated spermatogonia arranged in pseudotime (Fig. 5P), with Stra8 appearing in gene cluster 2 alongside a panel of prototypical spermatogonial genes [Gfra1, Id4, nanos C2HC-type zinc finger 3 (Nanos3), ret proto-oncogene (Ret), spalt like transcription factor 4 (Sall4), and zinc finger and BTB domain containing 16 (Zbtb16)] that exhibit elevated levels towards the end of pseudotime. Reciprocally, Cyp26b1 was found in gene cluster 1, which contained genes involved in several signaling pathways, including phosphatase and tensin homolog (PTEN), Netrin and G-protein-coupled receptor signaling (Fig. 5P). We also identified two other clusters of genes that exhibited maximal expression among prospermatogonia preceding the RA response in pseudotime (gene clusters 3 and 4), which could point to even more intrinsic regulators of RA responsiveness (Fig. 5P).

We previously reported single cell transcriptomes enriched for, or depleted of, transplantable SSCs from spermatogonial populations of P6 mice (Hermann et al., 2018). We compared these transcriptomes to those from prospermatogonia from P1.5 RA-treated testes to determine whether subsets of P1.5 prospermatogonia were more similar to SSC-enriched or -depleted fractions at P6. Two clusters comprising primarily P1.5 prospermatogonia (clusters 26 and 30) were located apart from P6 spermatogonia in principal component space (Fig. 6A,B). However, using unsupervised hierarchical clustering, we found that the P1.5 prospermatogonial clusters (clusters 26 and 30) were located immediately between clusters primarily comprising P6 SSC-enriched ID4-EGFP^bright spermatogonia and clusters primarily comprising P6 SSC-depleted ID4-EGFP^dim spermatogonia (Fig. 6C). Cluster 30 was adjacent to ID4-EGFP^bright and cluster 26 was adjacent to ID4-EGFP^dim, suggesting that these clusters represent precursors of the initial SSC and progenitor spermatogonia populations, respectively, seen at P5-6 (Fig. 6C,D). Differential gene expression analysis identified 450 genes that were expressed at significantly higher levels in cluster 30 than in cluster 26 (≥log2 fold-change; Fig. 6D, Table S2), including known SSC-related genes, such as Id4, Gfra1, Etv5, Ret and several others, whereas, reciprocally, there were 161 genes that were expressed at higher levels in cluster 26 than in cluster 30 (≥log2 fold-change; Fig. 6D, Table S2), including the progenitor marker Ndrg4 (Hermann et al., 2018). Pathway analysis demonstrated that P1 prospermatogonia in SSC-like cluster 30 expressed higher levels of proliferative genes, whereas prospermatogonia in cluster 26 (non-SSC-like) expressed higher levels of developmental genes (Fig. 6E). Together, these data indicate that a subset of P1.5 prospermatogonia (cluster 30) are proliferative and SSC-like, whereas the other subset, which are more RA-responsive, do not exhibit SSC-like characteristics.

Based on the results from the single cell RNA-seq analyses, we tested the hypothesis that CYP26-mediated catabolism prevented RARG+ undifferentiated progenitor spermatogonia from responding to RA. Multiple studies have shown a role for CYP26B1 in preventing the fetal male germline from responding to RA (Hogarth et al., 2015; Li et al., 2009; Bowles et al., 2006; MacLean et al., 2007). Results from KO studies suggest that CYP26B1 acts within multiple cell types in the developing testis (Hogarth et al., 2015; Li et al., 2009; MacLean et al., 2007), which is evident from immunolocalization (Fig. S3A-D). Therefore, we chose a pharmacological approach to rapidly and specifically inhibit CYP26 enzyme activity in all cell types both in vitro and in vivo. We utilized RA metabolism-blocking agents (RAMBAs) with either broad-spectrum CYP inhibitory activity (ketoconazole) or specific CYP26 inhibitory activity (talarozole) [reviewed in (Nelson et al., 2013)]. We first treated freshly isolated testicular cell suspensions from P5 mice with vehicle alone or a CYP inhibitor and then, 2 h later, added vehicle alone or RA for another 10 h (Fig. 7A). Similar to results shown in Fig. 1, RA increased the percentage of STRA8+ spermatogonia in primary cultures receiving no inhibitor pretreatment from 20% (vehicle alone, already STRA8+ when isolated) to 85%, similar to results when CYP26 inhibitors were added with RA (ketoconazole+RA=85%, talarozole+RA=86%, Fig. 7B,C,F). Primary cultures pretreated with either ketoconazole or talarozole exhibited a statistically significant 2-fold increase in the percentage of STRA8+ spermatogonia even without the addition of exogenous RA (Fig. 7D-F, Fig. S4A,B). These results reveal that, even in vitro, there was sufficient latent RA to activate STRA8 in ∼20% of previously STRA8– spermatogonia when CYP26 catabolism was inhibited. In addition, these results revealed that, even in the absence of cognate niches in vitro, ∼15% of the spermatogonial population remained RA-nonresponsive in a CYP26-independent manner (Fig. 7F).

The in vitro results provided proof of concept that CYP26 catabolism mediated RA responsiveness in freshly isolated spermatogonia in the absence of niches. However, the full extent of activity was difficult to determine because talarozole-treated samples in vitro lacked a consistent source of RA (because RA is labile, and retinol was not added to those samples). Therefore, we assessed the germ cell RA response in vivo 24 h after feeding P5 Id4-eGfp mice the CYP26 inhibitor talarozole (Fig. 7G). CYP26 inhibition induced a significant 6.5-fold increase in the numbers of STRA8+ spermatogonia (from 14 to 91%, Fig. 7H-J), and this increase was similar to treatment with exogenous RA (Fig. 1A-C). These results revealed that ∼60-70% of spermatogonia in the developing testis are poised to respond to RA, but are prevented from doing so by CYP26 catabolism.

Finally, we tested whether CYP26 activity prevented precursor prospermatogonia from responding to RA. We fed P1 pups talarozole or vehicle-alone for 6 h, and then examined testes 18 h later (Fig. 8A). In testes from control mice fed vehicle alone, <1% of prospermatogonia were STRA8+, as expected (Fig. 8B,F). In response to talarozole, only 3% of prospermatogonia became STRA8+ (Fig. 8C,F), which suggests that insufficient RA was available to prospermatogonia at P1-2 to direct differentiation. To assess this, we injected pups with vehicle or RA after the 6 h vehicle or talarozole feeding, and then euthanized mice 18 h later. When pups were fed vehicle and then injected with RA, 58% of prospermatogonia became STRA8+, and this increased to 81% when mice were fed talarozole (Fig. 8D-F). These results revealed that, although RA levels are minimal in the P1-2 testis, most neonatal prospermatogonia can respond to RA, and CYP26-mediated catabolism prevents precocious differentiation.

Here, using the mouse as a model system, we examined how postnatal male germ cells in the developing testis differ in their responsiveness to the differentiation signal provided by RA. This differential responsiveness allows a proportion of spermatogonia to retain stem cell potential (as SSCs) and proliferate as undifferentiated progenitor spermatogonia before differentiating in response to RA. Using approaches in vivo and in vitro with freshly isolated spermatogonia, we quantified the RA-responsive and RA-nonresponsive populations and defined the underlying gene expression signatures of each population. We report that the ability of germ cells to respond is correlated with the presence of RARG expression and absence of Cyp26b1 expression. Utilizing a pharmacological approach, we determined that CYP26-mediated catabolism prevents most postnatal male germ cells from responding to RA, except for SSCs, which remain RA-nonresponsive in a CYP26-independent manner. In addition, our study provides the first evidence of functional heterogeneity in the quiescent precursor prospermatogonia population, which reveals that spermatogonial fate might be determined before their emergence, in the form of differential responsiveness to RA.

The initial population of mouse spermatogonia is heterogeneous with regards to cellular fate, and already contains populations of stem, progenitor and differentiating spermatogonia (Niedenberger et al., 2015; Hermann et al., 2015; Yoshida et al., 2006). Here, we provided evidence that establishment of three functionally distinct populations of spermatogonia at P6 is based in part on their capacity to respond to RA and the presence and absence of an RA response: RA-competent/RA-responsive characterized differentiating spermatogonia; RA-competent/RA-nonresponsive characterized progenitors; and RA-incompetent/RA-nonresponsive characterized SSCs. As shown here and in other reports, RA signaling normally initiates in a subset (∼20%) of spermatogonia beginning at P3-4 (Busada et al., 2015b; Busada et al., 2014; Niedenberger et al., 2015; Snyder et al., 2010). However, the presence of exogenous RA drove most spermatogonia (85% at P5-6) to activate STRA8, proving their competence to respond to RA both in vitro and in vivo. Underlying these distinct spermatogonial populations are molecular controls that prevent the RA response in the RA-competent/RA-nonresponsive spermatogonia. Indeed, CYP26-mediated catabolism prevented a response among RARG+ undifferentiated progenitor spermatogonia. Such a catabolic buffer is irrelevant to putative SSCs, which are RA-nonresponsive and also RA-incompetent, because they lack detectable RARG expression in vivo. The striking similarity between our in vitro and in vivo results reveals that RA responsiveness is a cell-autonomous feature independent of the niche microenvironment. This conclusion contrasts with that of a recent study (Lord et al., 2018); however, in that study RA responsiveness was not measured in vivo in the presence of niches.

Here, we propose that the RA-incompetent/RA-nonresponsive spermatogonia are SSCs because all of these cells were ID4-EGFP^bright spermatogonia, a population that was previously shown to be highly enriched for regenerative activity by transplantation (Helsel et al., 2017). Our results indicate that most (>99%) ID4-EGFP^bright spermatogonia lacked detectable expression of RARG protein. These data agree with previous reports of RARG expression in progenitor and differentiating spermatogonia in the adult testis (Gely-Pernot et al., 2012) and prior transgenic mouse results in which forced expression of RARG in adult undifferentiated spermatogonia was sufficient for precocious RA-induced spermatogonial differentiation (Ikami et al., 2015). However, despite the congruency in those observations, a recent report speculated that RARG is present in SSCs; therefore, RARG expression did not distinguish competent and noncompetent cells (Lord et al., 2018). One explanation for this discrepancy is that the RARG measurements were made almost exclusively in a mixed ID4-EGFP+ spermatogonial population (Lord et al., 2018). As noted earlier, the ID4-EGFP+ population contains both SSC-enriched ID4-EGFP^bright, which we showed to be RARG–, as well as ID4-EGFP^dim progenitor spermatogonia, which we demonstrated are RARG+. Taken together, the compendium of evidence supports a model whereby RARG expression normally confers RA-response (differentiation) competency upon the progenitor subset of spermatogonia, whereas RARG expression is absent among SSCs and allows those spermatogonia to remain incompetent for an RA response. Therefore, a key regulatory step in establishing which SSC progeny self-renew and retain stemness and which attain differentiation competency is the transcriptional or post-transcriptional mechanisms controlling activation of RARG expression. However, it is curious that deletion of Rarg (or other RAR receptors singly or in combination) do not recapitulate lack of RA (VAD), but rather, lead to comparatively minor reproductive defects (Gely-Pernot et al., 2012). One possible explanation is that RARG normally directs spermatogonia to differentiate, but in the absence of RARG, alternative, compensatory mechanisms are activated to permit acquisition of RA competency and RA response to facilitate at least some spermatogonial differentiation (e.g. RARA compensates, as mentioned earlier). Moreover, because spermatogonial differentiation provides the pool of germ cells that enter meiosis, which is vital for male fertility, it is perhaps not surprising that acquisition of a crucial parameter, such as RA responsiveness, might be supported by multiple redundant receptor proteins.

Within the first day after birth in the mouse, quiescent precursor prospermatogonia in rodents reside in the center of the testis cords; besides subjective evidence of histomorphological heterogeneity (Drumond et al., 2011; Kluin and de Rooij, 1981), they have been widely considered to be a homogeneous population. Our observations of differential RA responsiveness provide the first clear evidence of functional heterogeneity among neonatal prospermatogonia as well as the full (transcriptome-wide) extent of molecular heterogeneity. This is a significant departure from fetal prospermatogonia, which are also not normally exposed to RA because of CYP26B1-mediated catabolism; however, fetal prospermatogonia from Cyp26b1-KO or ketoconazole-WT mice uniformly respond to RA by prematurely activating meiotic markers and dying by apoptosis before birth (Hogarth et al., 2015; Li et al., 2009; Yoshida et al., 2006; MacLean et al., 2007).

In humans, the transition from quiescent prospermatogonia occurs with the formation of two populations of spermatogonia (A_dark and A_pale) by 2-3 months after birth (Paniagua and Nistal, 1984). In the rodent testis, where this has been best studied, the initial spermatogonial population forms by P3-4 and is heterogeneous, containing distinct populations of undifferentiated SSCs and progenitor spermatogonia as well as differentiating spermatogonia that have already responded to RA (Busada et al., 2014; Zhou et al., 2008a; Niedenberger et al., 2015; Yoshid et al., 2006; Snyder et al., 2010). The mechanisms underlying the formation of the spermatogonial pool are unknown for any mammal, although two general scenarios can be envisioned, and these might not be mutually exclusive. In the first scenario, precursor prospermatogonia are equipotent and could adopt any spermatogonial fate (SSC, progenitor or differentiating) based on selection by extrinsic factors provided by the soma. There is currently little evidence supporting this selection mechanism in mammals. Developing and adult mammalian testes are an open stem cell niche environment in which SSCs appear to be randomly positioned throughout the seminiferous tubules. This is in stark contrast to the closed or polarized environment in Drosophila or Caenorhabditis elegans testes, in which spermatogenic cells differentiate as they move in a linear fashion away from a bona fide stem cell niche (Stine and Matunis, 2013; Yoshida, 2016). Therefore, if a niche microenvironment regulates SSC fate, then somatic cells surrounding that niche should provide differential signals; however, to the best of our knowledge, this heterogeneity has not been shown experimentally. Reciprocally, spermatogonia do differentially express receptors facilitating differential responses to widely distributed ligands, such as glial cell-derived neurotrophic factor (GDNF), fibroblast growth factor (FGF), Kit ligand (KITL), and RA, which regulate spermatogonial fate (Gely-Pernot et al., 2012; Gely-Pernot et al., 2015; Ikami et al., 2015; Blume-Jensen et al., 2000; Fouchecourt et al., 2006; Grasso et al., 2012; Ishii et al., 2012; Naughton et al., 2006; Packer et al., 1995; Schrans-Stassen et al., 1999; Sorrentino et al., 1991). In this study, RA responsiveness was examined in vivo (where male germ cells occupy specific niches) and in vitro (where niches are disrupted). The striking similarity in RA responsiveness in these two conditions supports the concept that germ cell RA responsiveness is a germ cell-autonomous and niche-independent feature. In the second scenario, the fates of the emergent spermatogonial subtypes at P3-4 are determined at some point earlier during development. Accordingly, precursor prospermatogonia would exhibit phenotypic and functional heterogeneity before P3-4 that predisposes them to acquire specific spermatogonial fates. Here, we provide the first evidence of transcriptome-wide phenotypic and functional heterogeneity in quiescent prospermatogonia at P1. Precursor postnatal prospermatogonia differed in their responsiveness to RA, revealing that molecular changes are already in place to modulate their behavior, or fate, as spermatogonia. We posit that the subset of prospermatogonia and newly formed spermatogonia that are RA-incompetent/RA-nonresponsive represent the foundational SSC pool, whereas those that are RA-competent/RA-nonresponsive are progenitors poised to differentiate during the first wave of spermatogenesis and RA-competent/RA-responsive spermatogonia are the first differentiating spermatogonia.

Our single cell RNA-seq results revealed that quiescent precursor prospermatogonia have differential expression of Cyp26b1. However, in developing testes, the highest detectable levels of CYP26B1 are in peritubular myoid cells (PTMs), which physically surround the testis cords. This suggests that a major role for CYP26B1 in the postnatal testis is akin to insulation around electrical wires, which prevents electricity from escaping the conducting metal. Using that analogy, CYP26B1 expression in PTMs would prevent RA produced by Sertoli cells within the testis cords (Tong et al., 2013) from escaping to the interstitium. It would also prevent RA produced in the interstitium, potentially by peritubular macrophages (which express RA-synthesizing enzymes), from indiscriminately reaching the developing spermatogonia (DeFalco et al., 2015). This system could be envisioned to provide intracord regulation of RA levels by preventing outside influences, and this concept could be tested by generating PTM-specific Cyp26b1-KO mice. It is unclear how discrete regions of testis cords could contain differentiating STRA8+/KIT+ spermatogonia (exposed to RA), whereas others have only undifferentiated spermatogonia (not exposed to RA) (Fig. 9). One possibility is that differential synthesis of RA by Sertoli cells occurs in a regional manner in cords, although no studies to date have reported heterogeneity of gene expression within the developing Sertoli cell population.

We found that 33% of ID4-EGFP^bright spermatogonia responded to RA to become STRA8+ at 6 h, and that the ID4-EGFP^bright population was reduced by one-third 24 h later. We interpret those results to mean that, at any one time in the neonatal testis, one-third of the ID4-EGFP^bright spermatogonia are responsive to RA (as progenitors), whereas the remaining two-thirds are true RA-nonresponsive SSCs; these results support those from previous studies from the Oatley et al. laboratory that have indicated that the ID4-EGFP^bright population is highly enriched for SSCs (Chan et al., 2014; Helsel et al., 2017; Oatley et al., 2011). The ∼30% of spermatogonia that did respond to RA by becoming STRA8 and losing ID4-EGFP^bright status probably represent those that are normally transitioning away from the self-renewing SSC state to become ID4-EGFP^dim progenitors, and might express RARG below the level of detection by IIF. A recent study from the Griswold laboratory reported the numbers of ID4-EGFP^bright spermatogonia and the numbers of functional SSCs in neonatal testes via transplantation of a suspension of Percoll-enriched germ cells in response to the following modulations of RA levels: (1) P2 pups were treated for 24 h with vehicle alone or RA; and (2) P2 pups were treated for 6 days with WIN 18,446 (which effectively blocks RA synthesis). With regards to the numbers of SSCs capable of colonizing recipient testes, there was no difference in response to exogenous RA [(1) above], and a very small increase in response to blocking RA synthesis [(2) above] (Agrimson et al., 2016). These results agree with our conclusions here, that the true SSC population is RA-nonresponsive.

Here, we have shown that the differential responsiveness of neonatal testicular germ cells to RA is established as early as P1, providing the first evidence to support the functional heterogeneity of the prospermatogonial population. We predict that differential responsiveness to RA is an intrinsic or cell-autonomous feature defining differentiating spermatogonia (that are responding), undifferentiated progenitor spermatogonia poised to differentiate (that can respond, but do not yet) and undifferentiated spermatogonia destined to become foundational SSCs (that cannot respond). This characteristic is then maintained throughout the remainder of steady-state spermatogenesis in the adult, because all spermatogonia must make the decision to differentiate in response to RA or remain an undifferentiated SSC at each stage VIII of the seminiferous epithelial cycle.

All animal procedures were performed in accordance with the National Research Council Guide for the Care and Use of Laboratory Animals and approved by the Institutional Animal Care and Use Committees of East Carolina University (Assurance A3469-01) and the University of Texas at San Antonio (Assurance A3592-01). All animals were maintained under conditions of ad libitum water and food with constant light-dark cycles. Analyses were done using male CD-1 mice (Charles River Laboratories) or Id4-eGfp mice on a C57Bl/6 background (Chan et al., 2014). The administration of exogenous RA was done as previously described (Busada et al., 2015b; Busada et al., 2014). Briefly, at least three neonatal mice received one subcutaneous injection of 100 μg all-trans RA (#R2625, Sigma-Aldrich) dissolved in 10 μl vehicle (DMSO) or vehicle alone at P5, which is after distinct populations of undifferentiated and differentiating spermatogonia have formed (Niedenberger et al., 2015), and were euthanized 6-24 h later. For oral administration, talarozole (Medchem Express) was also dissolved in DMSO and then diluted to a final concentration of 1 μg/μl in a solution of 5% polyethylene glycol 400 (Sigma-Aldrich) and 5% polysorbate-80 (Sigma-Aldrich). Talarozole was fed to P5 pups using a 24-gauge feeding needle, and mice were euthanized 24 h later on P6. Using dosages similar to a previous study (Puighermanal et al., 2009), we fed Id4-eGfp mice talarozole at 5 μg/g body weight. We determined that this was the lowest dosage that led to Stra8 activation without adverse systemic health problems. Control mice received an equivalent volume of vehicle alone, and each treatment experiment was repeated in triplicate.

Testes of P5 mice were used to generate single cell suspensions according to previously reported methods (Chappell et al., 2013). Briefly, testes were removed from neonatal mice and placed in Krebs-Ringer Bicarbonate Buffer (EKRB). Testes from at least five pups from a given litter were detunicated and transferred to a solution containing 4.5 ml 0.25% trypsin (Gibco) and 0.5 ml DNase1 (7 mg/ml, Sigma-Aldrich) and incubated in a 37°C water bath for 3 min. An additional 1 ml DNase1 was added and the mixture was pipetted up and down to further break up the testes, followed by a 1-min incubation in a 37°C water bath. Then, 1 ml fetal bovine serum (FBS, ATCC) was added to deactivate the trypsin, and the mixture was triturated to break apart any remaining clumps of intact testis. The cell mixture was filtered through a 40-μm sieve and centrifuged at 500×g for 7 min. The resulting cell pellet was suspended in DMEM media (Gibco) containing 1X penicillin-streptomycin (Thermo-Fisher Scientific). Cells were plated at a density of 1.4×10⁶ cells/well into a 12-well plate containing glass coverslips. The cells were allowed to adhere for 2 h at 37°C with 5% CO₂ before treatment. Cells were treated with either vehicle control (DMSO) or the CYP inhibitors ketoconazole (5 μM, Sigma-Aldrich) or talarozole (25 nM) for 2 h, followed by the addition of either vehicle alone or RA (1 μM) for 8 h. These concentrations were calculated based on published IC₅₀ values and titered to achieve effective inhibition of the CYP enzymes at the minimum dosage (Thatcher and Isoherranen, 2009). These cells were in culture for a total of 12 h. Each cell culture experiment was repeated in triplicate.

Indirect immunofluorescence (IIF) was performed as previously described (Busada et al., 2015b; Busada et al., 2014; Niedenberger et al., 2015). Briefly, testes from at least three different mice were fixed in 4% PFA pH 7.2 for 4.5 h at 4°C, washed in 1× PBS for at least 8 h, and incubated in 30% sucrose overnight at 4°C. Testes were embedded in optimal cutting temperature compound (OCT), frozen and then cut into 5-μm sections. Sections were incubated with primary antibodies diluted in 1× PBS (see Table S3) for 1 h at room temperature. Primary antibody was omitted in the negative controls. Sections were washed in 1× PBS+Triton X-100 and then incubated with either Alexa Fluor anti-rabbit or anti-goat secondary antibody (1:500, Invitrogen) and phalloidin-635 (1:500, Invitrogen) for 1 h at room temperature. Coverslips were mounted with Vectastain containing DAPI to visualize any nuclei (Vector Laboratories). Sections were imaged using a Fluoview FV1000 confocal laser-scanning microscope (Olympus America). Each IIF was repeated in triplicate.

Determination of spermatogonia expressing a given protein in immunostained frozen tissue sections was carried out as previously described (Busada et al., 2015a, b). Briefly, germ cells within 21-30 testis cords were counted from at least three different animals. Cells were identified as positive for a marker by the threshold tool in Image J (US National Institutes of Health) using the default algorithm. Intensity thresholds were determined by averaging the thresholds of the control images in which all cells were negative. Threshold ranges were as follows: DDX4=100-255, STRA8=63-255 and EGFP=126 (EGFP is measured as epifluorescence rather than being visualized by immunostaining, and assigned a single threshold to maintain the boundary between EGFP^bright and EGFP^dim cells). The channel-splitting tool in ImageJ was used to count ID4-EGFP/RARG spermatogonia. Cells that were RARG+ (red) and ID4-EGFP– (not green) were counted. The same method was used for counting ID4-EGFP/KIT. Quantification of testicular single cell suspensions followed a similar process. Cells from 5-10 frames per replicate were identified as positive, as described earlier. Cells were quantitated based on the number of STRA8+ cells divided by the number of DDX4+ spermatogonia, which marks all germ cells in the neonatal testis (Niedenberger et al., 2015). Photomicrographs were captured with an Axio Observer A1 microscope (Carl Zeiss Microscopy, LLC) equipped with an XL16C digital camera and Exponent version 1.3 software (Dage-MTI).

Cell suspensions were loaded into Chromium microfluidic chips with 3′ v2 chemistry and used to generate single cell gelbead emulsions (GEMs) using the Chromium controller (10x Genomics) essentially as described (Zheng et al., 2017). Suspensions containing ∼17,400 cells were loaded on the instrument with the expectation of collecting up to 10,000 GEMs containing single cells. GEM-RT was performed in a T100 Thermal cycler (Bio-Rad) and all subsequent steps to generate single cell libraries were performed according to manufacturer recommendations. Libraries were sequenced at the Genome Sequencing Facility (GSF) at Greehey Children's Cancer Research Institute in UT Health San Antonio (UTHSA) on a HiSeq3000 instrument (Illumina). Trimmed FASTQ files (26 bp Cell barcode and UMI Read1, 8bp i7 index, and 100 bp Read2) were generated using the CellRanger mkfastq command (a 10x Genomics wrapper around BCL2Fastq). Primary data analysis (alignment, filtering and UMI counting) to determine gene transcript counts per cell (producing a gene-barcode matrix), quality control (QC), clustering and statistical analysis were performed using CellRanger count (10x Genomics) and the mm10 genome assembly/annotation reference. Outputs from replicate independent samples of single cells were combined using CellRanger aggr (10x Genomics) based on mapped read counts to normalize sequencing depth and produce aggregated gene × cell barcode matrices and clustering models. The Loupe Cell Browser v2.0.0 (10x Genomics) was used to visualize results. QC metrics from these data are detailed in Table S4. All newly generated genomics data reported in this study were deposited in the Gene Expression Omnibus (GEO) and Sequence Read Archive (SRA) databases under the accession number GSE117437. To facilitate the ease of interrogation of these data by the scientific community, analyzed data sets in Loupe Cell Browser format were deposited in the Mendeley Data repository (data.mendeley.com) under doi:10.17632/5d3h9dwpzs.1. Users should download the appropriate Loupe Cell Browser version (see support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest) to facilitate queries of these data. These data are also available in analyzed format through the ReproGenomics Viewer (rgv.genouest.org/).

Raw count single cell gene expression matrices were imported into the R package Seurat v1.4.03 (Satija et al., 2015), filtered [cells expressing ≥200 detected genes and ≥2500 unique molecular identifiers (UMIs), genes expressed in ≥3 cells, and cells with <10% transcripts corresponding to mitochondrial genes] and gene expression values were log normalized and scaled. Unsupervised cell clustering and dimensionality reduction with t-SNE were performed in Seurat based on the statistically significant principal components. The top-ten differentially expressed genes (marker genes) of each cell cluster were determined by log fold change ≥0.25 using a default Wilcox test. After identification of cell clusters, raw count matrices without cluster subsets (e.g. without testicular somatic cells) or sample subsets (e.g. RA only) were re-analyzed with Seurat and subsequently imported to Monocle2 (Qiu et al., 2017; Trapnell et al., 2014) for additional combined t-SNE, unsupervised density peaks clustering and trajectory analyses (genes expressed in ≥10 cells, minimum expression of 0.1). The top 1000 differentially expressed genes (DEGs; q value <0.01) over pseudotime were used to perform trajectory analysis, which ordered cells in pseudotime by cell state. A pseudotime heatmap was plotted using these DEGs using the ward.D2 method (Hermann et al., 2018). Lists of differentially expressed genes were analyzed by Ingenuity Pathway Analysis [Qiagen, Build 481437M, content version 44691306 (6/2018)] to identify biological pathways that were significantly over-represented among the genes in each list.

The P1.5 RA-treated data sets were merged with previously published single cell transcriptomes from sorted P6 ID4-EGFP^bright and ID4-EGFP^dim spermatogonia (Hermann et al., 2018), comprising a total of 31,830 cells; data were analyzed in Seurat with the same QC and scaling as noted earlier. Linear dimensional reduction was performed using the RunPCA function based on 2195 variable genes, principal components (PC) 1-10, and a resolution of 0.6. This resulted in 19 clusters, among which clusters 0, 3, 5 and 16 were identified as germ cell clusters, using expression of germ cell markers (Ddx4 and Gata4) visualized as violin plots. This resulted in a log-normalized gene expression matrix of 9922 germ cells and 21,324 genes for subsequent analysis. Following preprocessing by Seurat, the germ cell clusters were analyzed in Pagoda2 (Fan et al., 2016) using the default QC function. Gene variance normalization was performed, followed by dimensional reduction by PC analysis (PCA; nPcs=50). Subsequently, a k-nearest neighbor (KNN) graph identified 50 cell clusters in PCA space and embeddings were then generated using largeVis based on PCA reduction. Differential gene expression of each cluster (≥log2 fold change) and pathway overdispersion were determined to generate an unsupervised hierarchical clustering dendrogram. Differentially expressed gene lists between clusters 30 and 26 were determined within the PAGODA2 web app using a Mann–Whitney U test and were used for gene set enrichment analysis in Metascape (http://metascape.org/).

Statistical analyses of cell counts were performed using the Student's t-test, and the level of significance was set at P≤0.05.

The authors thank Dr Douglas Weidner (ECU) for technical assistance with flow cytometry and Dr Jon Oatley (Washington State University) for providing the Id4-eGfp mice. The authors also appreciate computational support from the UTSA Research Computing Support Group and the high-performance computing cluster ‘Shamu’ operated by the UTSA Office of Information Technology. Single cell transcriptome samples were generated in the UTSA Genomics Core (supported by the UTSA NIH grant G12-MD007591 and NSF grant DBI-1337513). Sequencing data were generated in the Genome Sequencing Facility [supported by UTHSA, NIH/NCI P30-CA054174 (Cancer Center at UTHSA), NIH Shared Instrument grant S10-OD021805 and a CPRIT Core Facility Award RP160732].