Chromatin and Single-Cell RNA-Seq Profiling Reveal Dynamic Signaling and Metabolic Transitions during Human Spermatogonial Stem Cell Development

Abstract

Human adult spermatogonial stem cells (hSSCs) must balance self-renewal and differentiation. To understand how this is achieved, we profiled DNA methylation and open chromatin (ATAC-seq) in SSEA4⁺ hSSCs, analyzed bulk and single-cell RNA transcriptomes (RNA-seq) in SSEA4⁺ hSSCs and differentiating c-KIT⁺ spermatogonia, and performed validation studies via immunofluorescence. First, DNA hypomethylation at embryonic developmental genes supports their epigenetic "poising" in hSSCs for future/embryonic expression, while core pluripotency genes (OCT4 and NANOG) were transcriptionally and epigenetically repressed. Interestingly, open chromatin in hSSCs was strikingly enriched in binding sites for pioneer factors (NFYA/B, DMRT1, and hormone receptors). Remarkably, single-cell RNA-seq clustering analysis identified four cellular/developmental states during hSSC differentiation, involving major transitions in cell-cycle and transcriptional regulators, splicing and signaling factors, and glucose/mitochondria regulators. Overall, our results outline the dynamic chromatin/transcription landscape operating in hSSCs and identify crucial molecular pathways that accompany the transition from quiescence to proliferation and differentiation.

Keywords: DNA methylation; hormone receptors; human spermatogonial stem cells; metabolism; open chromatin; pluripotency; single-cell RNA-seq; spermatogenesis.

Graphical abstract

Figure 1

Genomic Profiling of Human Spermatogonial Stem Cells (A) Schematic illustration of human adult male germline development and niche, depicting a small section of the seminiferous tubule. (B) Experimental workflow in this study. sc, single cell; WGBS, whole-genome bisulfite sequencing. (C) Expression profiles of selected key genes following the enrichment procedure with SSEA4 (blue) or c-KIT (red). Browser snapshots of DDX4 (germ cell marker), FGFR3 (hSSC marker), KIT and SYCP3 (differentiating spermatogonia marker), GATA4 (Sertoli cell marker), and LHCGR (Leydig cell marker). The intron/exon (box) genomic structure of each gene is shown in black. (D) Distribution of DNAme in human PGCs, hSSCs, sperm, egg, ICMs (inner cell mass), ESCs, FC (frontal cortex), and liver. Human PGC and liver methylation data are from Guo et al. (2015); ICM and FC methylation data are from Guo et al. (2014a); egg methylation data are from Okae et al. (2014); ESC methylation data are from Gifford et al. (2013). (E) Hierarchical clustering of correlation of global DNAme in human PGCs, hSSCs, sperm, egg, ICMs, ESCs, FC, and liver. See also Figures S1 and S2.

Figure 2

Unique Chromatin Landscape in hSSCs Revealed by ATAC-Seq (A) Heatmap of k-means clustering (n = 4) showing ATAC-seq signals at ESC and hSSC peaks and motifs enriched in each cluster. (B) Distance between NFY sites, DMRT1 sites, and HRE sites. (C) Expression of genes adjacent (within 10 kb) to DMRT1 sites, NFY sites, and HRE sites are specifically upregulated in hSSCs. See also Figures S3 and S4.

Figure 3

Chromatin Underlying Poised Pluripotency in hSSCs (A) Principal component analysis (PCA) of the transcriptome of SSEA4⁺ hSSCs, c-KIT⁺ spermatogonia, PGCs, and ESCs. Human PGC and ESC RNA-seq data are from Gkountela et al. (2015). (B) Hierarchical clustering of pluripotency-related factors from ESCs, PGCs, SSEA4⁺ hSSCs, and c-KIT⁺ spermatogonia. (C) Browser snapshots of ATAC-seq signals and DNAme at selected key genes. Note: POU5F1, NANOG, and SOX2 encode core pluripotency factors; SALL4, TCF3, KLF4, KLF2, STAT3, and MBD3 encode ancillary pluripotency factors; DDX4 and DAZL are germ cell-specific markers. See also Figure S5.

Figure 4

Single-Cell Transcriptome Analysis by t-SNE and Monocle (A) t-SNE analysis plot of single-cell transcriptome. t-SNE, t-distributed stochastic neighbor embedding. (B) Expression profiles of selected key genes in SSEA4-enriched or c-KIT-enriched single cells projected on the t-SNE plot. (C) Monocle analysis plot of scRNA-seq data, in which gene expression in multi-dimensional space is compressed to two dimensions/components. Most cells were positioned along a central branch, with two small branches emanating at the transition between SSEA4⁺ and c-KIT⁺ cells. The states assignment involved subsequent hierarchical clustering shown in Figure 6B. (D) Expression of selected key genes along pseudotime development. SSEA4⁺ (blue) or c-KIT⁺ (red) cells are projected along pseudotime. Genes associated with self-renewal are depicted on the left column, and genes associated with c-KIT⁺ differentiating cells along the right column. Note: depicted as compressed (log10) transformed expression data, and only 30%–70% of single cells typically provide non-zero expression of individual genes. See also Figure S5.

Figure 5

Signaling Pathways Differentially Expressed in SSEA4⁺ hSSCs or c-KIT⁺ Spermatogonia (A) Expression levels of different cell signaling pathway components and other key genes along pseudotime. (B) Schematic summary of signaling pathways singled-out by our RNA-seq analysis, with the ligands currently used in mouse SSCs cultures (outset box). Note: this schematic is based on detection of RNA transcripts and potential signaling activity, not flux measurements. See also Figure S6.

Figure 6

Monocle and Clustering Analyses Reveal Four Cell States For a Figure360 author presentation of Figure 6, see the figure legend at http://dx.doi.org/10.1016/j.stem.2017.09.003. (A) K-means clustering (n = 4) of genes exhibiting differential expression in SSEA4⁺ hSSCs versus c-KIT⁺ spermatogonia. Note: each row represents a gene, and each column represents a single cell, with columns/cells placed in pseudotime order (as defined in Figure 4C). Gene expression levels utilize a Z score, which depicts variance from the mean. (B) Hierarchical clustering of state 2 and state 3 cells from Figure 6A. Note: columns (cells) were re-ordered using hierarchical clustering, while genes (rows) were kept in the same order as Figure 6A. These state assignments were then used to refine the identity and trajectory of the minor branches highlighted on the Monocle plot in Figure 4C. (C) Summary schematic of the combinatorial distribution of the four gene expression clusters combine to define four distinct cellular states and proposed dynamic ordering model based on gene identities and pathways. (D) Violin plots of representative genes from each gene cluster, and their relative expression levels in each cellular state. y axis represents Z score of expression levels. Here, the mean levels for each state are linked by lines to depict the developmental trajectory. Sig., signaling; TF, transcription factor; Glu. inh., glucose inhibition; Diff., differentiation.

Figure 7

Validation of scRNA-Seq by Immunostaining of Human Seminiferous Tubules Immunolocalization of FGFR3 (cluster A marker used as a surrogate marker to SSEA4, in green), c-KIT (cluster D marker, in red), and 11 different cluster-specific antigens (in blue) on formalin-fixed paraffin embedded (FFPE) sections of human seminiferous tubules. Each antigen (name in blue on the left side) is represented by three panels (left, co-staining with FGFR3; middle, co-staining with c-KIT; right, triple co-IF staining). For clusters A and B, the cluster-specific antigens (blue) are expressed in FGFR3⁺ c-KIT^– cells, while, in clusters C and D, the cluster-specific antigens (blue) are expressed in FGFR3^– c-KIT⁺ cells. The bottom right tryptic represents the negative (no primary) controls. All pictures are at the same magnification, and the white bar in the top left panel is 10 μm. See also Figure S7.