Single-Cell RNA Sequencing of Human, Macaque, and Mouse Testes Uncovers Conserved and Divergent Features of Mammalian Spermatogenesis

Abstract

Spermatogenesis is a highly regulated process that produces sperm to transmit genetic information to the next generation. Although extensively studied in mice, our current understanding of primate spermatogenesis is limited to populations defined by state-specific markers from rodent data. As between-species differences have been reported in the duration and differentiation hierarchy of this process, it remains unclear how molecular markers and cell states are conserved or have diverged from mice to man. To address this challenge, we employ single-cell RNA sequencing to identify transcriptional signatures of major germ and somatic cell types of the testes in human, macaque, and mice. This approach reveals similarities and differences in expression throughout spermatogenesis, including the stem/progenitor pool of spermatogonia, markers of differentiation, potential regulators of meiosis, RNA turnover during spermatid differentiation, and germ cell-soma communication. These datasets provide a rich foundation for future targeted mechanistic studies of primate germ cell development and in vitro gametogenesis.

Keywords: evolution; meiosis; signaling; single-cell RNA-seq; soma; spermatogenesis; spermatogonial stem cell.

Figure 1.. Overview of major cell types and cellular attributes inferred from single-cell RNA-Seq analyses of human and macaque testes

A. Workflow overview of data collection and analysis. B. Visualization of major testis cell types from global clustering of all 13,837 human (left) and 21,574 macaque (right) cells in UMAP (top). Inset: Focused clustering of 3,722 human (left) and 2,098 macaque (right) somatic cells identifies 7 distinct cell types (bottom). C. Heatmap of marker gene expression in the 11 major cell types for human (left) and macaque (right). Shown are each gene’s mean expression in each cell type (i.e., the centroid), standardized over the 11 centroids. Representative markers are listed. D. Distribution profiles of per-cell attributes compared across the 11 cell types, including mouse data from Green et. al. 2018. From top to bottom: nGene, total number of detected genes per cell; nUMI, total number of Unique Molecular Identifiers (UMI) per cell; %ChrX, percent of total UMIs from genes on the X chromosome. E. Putative MSCI escapee genes, shown as red dots, in context of other highly variable genes (in gray). X and Y axis show the mean expression level in Scyte and SPG, respectively. The vertical line demarcates genes with an average expression >= 0.5 log normalized UMI. The diagonal solid line highlights genes with a fold increase >2 from spermatogonia to spermatocytes. F. Intron FISH localizes RIBC1 transcripts (magenta) to the XY-body (H2AX, green) in macaque. Arrows indicate co-localization. Scale bars 10um. See also Figure S1 and Table S1.

Figure 2:. Merged analysis of mouse, human and macaque spermatogonia identifies six molecular states (SPG1–6).

A. Re-clustering of 1688, 747, 2174 SPG cells from human, macaque and mouse, respectively identifies six molecular states, shown in t-SNE space. B. Distribution of SPG cells across the six states, colored by species. C. Biological annotation of SPG1–6 based on expression patterns of established conserved markers for different stages of spermatogonia. Shown are expression levels in individual cells without distinguishing by species. D-E. Expression patterns of several newly identified noncanonical markers of spermatogonia states. Each column is for a gene, showing the expression heatmap over all cells without distinguishing by species (top), and the corresponding violin plots for each species (bottom). Note that VCX genes are primate-specific but lack annotated 1–1 orthologs in macaque. F. Similarity among SPG states within and across species, shown as heatmaps of average Jaccard index among cell clusters. G. Graphic summary of between-species comparison of state-dependent expression patterns for select marker genes, colored by species. See also Figure S2 and Table S3.

Figure 3.. Localization and functional characterization of primate undifferentiated spermatogonia.

A-D. Immunofluorescence co-staining analysis of SPG1–2 markers, TCF3 (A), MAGEB2 (B), CDK17 (C), and MORC1 (D), with undifferentiated spermatogonia marker UCHL1 (top) or differentiating spermatogonia marker cKIT (bottom) in human seminiferous epithelium cross-sections. Dot plots show the number of single positive (left, right), or double positive (middle) cells per tubule. Error bars represent standard deviation. E-H. Co-localization analysis of SPG markers PIWIL4 (E), MORC1 (F) and TCF3 (G) with human spermatogonia subtypes: Adark, Apale, B-type, and Spermatocytes. Dot plots summarize the number of marker positive cells per cross section for each SPG subtype. H demonstrates background staining. Error bars represent standard deviation. I. Workflow of human to mouse xenotransplantation experiments. J. FACS gating strategy to enrich for TSPAN33 cells. K. Number of colonies obtained from TSPAN33-negative and TSPAN33-positive cells, per 10⁵ cells transplanted. At least 10 testes were counted per fraction from 3 different donors. Error bars represent standard deviation. L. Normalization of colonies from xenotransplantation of FACS sorted cell fractions by total fraction size as estimated by FACS. Error bars represent standard deviation. All scale bars = 50 μm. See also Figure S3.

Figure 4.. Three-species comparison of gene expression dynamics across the germ cell differentiation trajectory.

A. Pseudotime heterochrony between pairs of species, shown as correlation matrices between 200 ordered centroids in each species, defined from initial pseudotime assignments in each species. SC, spermatocyte. Elong, elongating spermatid. B. Construction of a consensus pseudotime by principal curve (black line). The spermatogonia centroid is shown as an anchor (triangle) to indicate the direction of differentiation. Unaligned pseudo-time bins at the end of macaque trajectory are indicated as crosses. Size of symbol indicates the number of cells contained by each bin. C. Rank correlation for pairs of species using 20 pseudotime bins after aligning to the consensus pseudotime. D. Relative cell abundance across 20 aligned pseudotime bins (GC1–20), preceded by 6 spermatogonia bins representing consensus SPG 1–6. Size of points indicate number of cells contained by each bin. Mouse data used only the cells from non-enriched experiments. E. Expression patterns of established marker genes across the 6 SPG states and 20 germ cell states allow mapping between the molecular states and the major cell types/processes, showing on the top of 4D. See also Figure S4 and Table S4.

Figure 5.. Three-species comparison of dynamically expressed genes in the germ cells.

A. Comparison of dynamically expressed genes between pairs of species. Rows indicate genes, ordered by six gene clusters; and columns indicate the 20 aligned bins for germ cells, plus one SPG bin. Grey lines link orthologous pairs. Genes lacking connecting lines have orthologs that are not highly expressed or variable in the other species. B. Number of genes with “shifted” phase between species, with the direction of change shown in arrows, labeled with the main biological processes (gene ontology) enriched for. C. Between-species similarity of transcriptome across the 6 SPG and 20 GC states. Values are rank correlations of gene expression centroid between human and macaque (blue) or human and mouse (pink), for either all 11,023 1-1-1 orthologs (solid line) or 835 transcription factor orthologs (dotted line). Line indicates loess best fit. TF, transcription factor. D. Values are rank correlations of gene expression centroid between human and macaque (blue) or human and mouse (pink) for 835 transcription factor orthologs (dash line) or 835 randomly subsampled orthologs (solid line) with matched distribution of expression levels as the TFs. E. Group-level gene expression changes across the 26 germ cell molecular states, for several mutually exclusive groups of orthologs. Venn diagram (top left) shows the groups of genes that are unique to one species (α’s), with orthologs in two species (β’s), or in three species (γ). Group δ are 1-1-1 orthologs. Graphs show the percentage of total expression accounted for by 1-1-1 orthologs (top right), species-specific genes (bottom left), and primate-specific genes (bottom right). See also Figure S5 and Table S5.

Figure 6.. Somatic cells show divergence of transcriptome and signaling relationships.

A. Expression of a representative marker for each of the 7 somatic cell types for human (top) and macaque (bottom). B. Schematic of somatic cell localization patterns in primate testis. C. Hematoxylin and eosin staining of human, macaque, and mouse seminiferous tubules. Arrowheads indicate peritubular myoid cell nuclei; red dotted lines indicate layers of myoid cells. D. Cross-species comparison of somatic transcriptome by a joint PCA of cell type centroids from all 3 species, using 673 1-1-1 ortholog genes that are highly variable in each species. Symbol shape represents the species. Symbol color represents the cell type. Ellipsoids group cell types with similar transcriptome across species. Elliptic lines depict 50% confidence interval for cell types observed across species. E. Cross-species comparison of somatic cell types using rank correlation of gene expression centroids for 7 human, 7 macaque, and 6 mouse cell types. Correlations are calculated using two gene sets: 673 orthologs described in Figure 6D (upper right triangle) or 311 1-1-1 ortholog genes encoding signaling ligands (lower left triangle). See also Figure S6 and Table S6.

Figure 7.. Potential soma-germline signaling relationships by ligand-receptor analysis.

A. Summary of the extent of potential interactions between ligands expressed in somatic cells (left) and receptors expressed in germ cells (right), counting top 5% of all interactions for each species. Symbol size indicates the number of receptor-ligand interactions contributed by a cell type; and line width shows the number of interactions between the two cell types. B. Pattern of expression for selected ligand-receptor pairs. Symbol size indicates the level of expression; line width shows the expected signaling strength. Predicted interactions are shown as arrows (details in Methods). See also Figure S7 and Table S7.