Sequencing and QC statistics for single-nucleus multiome libraries of E11.5-E13.5 XX and XY gonads.

Single-nucleus multiomics sequencing of E11.5-E13.5 XX and XY PGCs.

(a) Workflow of single-nucleus multiomics analyses. (b-c) Dimensional reduction (UMAP) and cell-type specific clustering of snRNA-seq data collected from E11.5-E13.5 XX and XY PGCs. Each dot represents an individual cell color coded by either embryonic stage and sex (a) or cluster number (c). (d) Cell type composition of snRNA-seq clusters. The embryonic stage and sex of the PGC population is indicated on the x-axis and the unbiased clustering number is represented on the y-axis. (e-f) UMAP of snATAC-seq data collected from E11.5-E13.5 XX and XY PGCs color coded by either embryonic stage and sex (e) or cluster number (f). (g) Cell type composition of snATAC-seq clusters. (h-i) UMAP visualization of E11.5-E13.5 XX and XY PGCs using a joint, or weighted, measurement of snRNA- and snATAC-seq modalities. (h) Trajectory of transcriptomic changes (also known as pseudo-time analysis in Monocle) among PGCs overlaid on the joint UMAP. Arrows indicate the path of the developmental trajectory. Numbers on the pseudotime trajectory indicate the following: ‘0’ represents the starting cell population; ‘1’ indicates the branching point in the trajectory path; and ‘2’ represents the terminal cell populations of the trajectory paths. Each PGC is color coded by embryonic stage and sex. (i) Unbiased clustering of PGCs using combined snRNA- and snATAC-seq data. Each PGC is color coded by cluster number. (j) Cell type composition of joint UMAP clusters.

Molecular characterization of genes known for their roles in male PGC differentiation.

(a) Dotplot of the average expression and percentage of cells expressing Rb1, Rbl2, Cdkn1b, Cdkn2b, Bnc2, Cnot1, Dnd1, Nanos2, Nanos3, and Zkscan5. The DNA binding motif for ZKSCAN5 is indicated above. The color scale represents the average expression level, and the size of the dot represents the percentage of cells expressing the gene. (b) Left: Coverage plot of the normalized snATAC-seq signal at the Bnc2 locus. Right: Violin plot of Bnc2 expression in E11.5-E13.5 XX and XY PGCs. (c) Joint UMAP showing Bnc2 expression levels for E11.5-E13.5 XX and XY PGCs. Individual cells are color coded by the level of Bnc2 expression. (d) Candidate factors (x-axis) regulating Bnc2 expression based on their Cistrome DB regulatory potential score (y-axis). (e) Dotplot of the average expression and percentage of cells expressing the candidate factors that potentially regulate Bnc2 expression. The color scale represents the average expression level, and the size of the dot represents the percentage of cells.

Molecular characterization of genes known for their roles in female PGC sex determination.

(a) Dotplot of the average expression and percentage of cells expressing Rec8, Sycp2, Sycp3, Ctnnb1, Meioc, Rnf2, Stra8, Ythdc2, Zglp1, and Dmrt1. The DNA binding motif for DMRT1 is indicated. The color scale represents the average expression level, and the size of the dot represents the percentage of cells expressing the gene. (b) Left: Coverage plot of the normalized snATAC-seq signal at the Stra8 locus. Right: Violin plot of Stra8 expression in E11.5-E13.5 XX and XY PGCs. (c) Joint UMAP showing Stra8 expression levels for E11.5-E13.5 XX and XY PGCs. Individual cells are color coded by the level of Stra8 expression. (d) Candidate factors (x-axis) regulating Stra8 expression based on their Cistrome DB regulatory potential score (y-axis). (e) Dotplot of the average expression and percentage of cells expressing the candidate factors that potentially regulate Stra8 expression. The color scale represents the average expression level, and the size of the dot represents the percentage of cells.

Identification of differentially expressed genes and differentially accessible peaks in E11.5-E13.5 XX and XY PGCs.

(a) Bar graph showing the number of differentially expressed genes (DEGs) between XX and XY PGCs at E11.5-E13.5. (b) Bar graph showing the number of differentially accessible chromatin peaks (DAPs) between XX and XY PGCs at E11.5-E13.5. (c) Left: Coverage plot of the normalized snATAC-seq signal at the DEG Porcn locus. Peak-to-gene linkages are indicated by the ‘links’ line. Right: Violin plot of Porcn expression in E11.5-E13.5 XX and XY PGCs. (d) Left: Coverage plot of the normalized snATAC-seq signal at the DEG Rimbp1 locus. Peak-to-gene linkages are indicated by the ‘links’ line. Right: Violin plot of Rimbp1 expression in E11.5-E13.5 XX and XY PGCs.

Analysis of gene networks in E11.5-E13.5 XX and XY PGCs.

(a) Illustration of the analytical flow to infer gene networks (i.e. transcription factors (TFs) and their predicted target genes) involved in PGC sex determination. (b-c; i-j) Line plots showing the expression levels and in silico chromatin binding, termed ‘motif activity’, of Tfap2c (b), Tcfl5 (c), Foxk2 (i), and Pou6f2 (j) in E11.5-E13.5 XX (pink) and XY (teal) PGCs. The p-value of motif enrichment and the TF binding motif are indicated. (d and k) Enrichment analysis of biological process gene ontology (GO) of the predicted target genes for TFAP2c (d) and FOXK2 (k). The color scale represents the p-value of GO term enrichment, and the size of the dot indicates the gene ratio for each GO term. (e and l) The number of predicted target genes for enriched TFs in XX PGCs (e) and XY PGCs (l). (f) Circos plot showing whether the motifs for TFAP2C, TCFL5, ZFX, MGA, or NR6A1 (indicated in black) are detected in the regulatory loci linked to Tbx4, Nr6a1, Mga, Tcfl5, or Tfap2c (indicated in grey). (m) Venn diagram showing the overlap of predicted target genes for POU6F1/2, FOXK2, and FOXO1. (g-h) Icon array showing the presence (pink/purple) or absence (light grey) of the motifs for TFAP2c, TCFL5, ZFX, MGA, NR6A1, TBX4, and GATA2 in peaks linked to genes related to meiosis (g) or WNT signaling (h). (n-o) Icon array showing the presence (teal) or absence (light grey) of the motifs for FOXO1, FOXK2, and POU6F1/2 in peaks linked to genes related to mitosis (n) or signal transduction (o).

Sexually dimorphic interactions between PGCs and gonadal supporting cells.

(a) Illustration of cell communication between PGCs and gonadal supporting cells. (b) Number of significant interactions or ligand-receptor pairs (p-value < 0.05) between supporting cells and PGCs in E11.5-E13.5 XX and XY gonads. (c) Number of significantly enriched ligand-receptor pairs, or interactions, per signaling pathways induced by either WNTs, inhibins/activins, or BMPs in E13.5 XX and XY gonads. (d, f, & h) Venn diagrams showing the overlap of significant interactions between supporting cells and PGCs in XX and XY gonads at E11.5 (d), E12.5 (f), and E13.5 (h). (e, g, & i) Dotplots showing the mean expression for all interacting partners in representative ligand-receptor pairs in XX and XY gonads at E11.5 (e), E12.5 (g), and E13.5 (i). Ligand-receptor pairs in pink represent XX-specific interactions, ligand-receptor pairs in blue represent XY-specific interactions, and ligand-receptor pairs in grey represent interactions found in both XX and XY gonads. The size of the dot represents the mean expression value.