Esr1 in MPOAVgat+Esr1+ neurons governs adolescent maturation of sexual behaviors

a, Schematic timeline for behavioral experiments testing the role of Esr1 in MPOAVgat+Esr1+ during adolescent maturation of mating behaviors in Slc32a1Flp::Esr1lox/lox mice (Vgat-Esr1KO), refers to panels b-h.

b, Representative images of viral-reporter and Esr1 expression from Control (left) and Cre (right) groups. Quantification of Esr1 in the MPOA.

c-e, Quantitative comparisons of female mice: first vaginal opening (VO) age (c), estrous age (d, left), receptive age (d, right), number of intromissions received (e, left), and receptivity (e, right).

f-h, Quantitative comparisons of male mice: balanopreputial separation (BPS) age (f), first mount age (g), number of mounts (h, left) and thrusts (h, right).

i, Schematic timeline for behavioral experiments testing the role of Esr1 in MPOAVglut2+Esr1+ during adolescent maturation of mating behaviors in Slc17a6Flp::Esr1lox/lox mice (Vglut2-Esr1KO), refers to panels j-p.

j, Representative images of viral-reporter and Esr1 expression from Control (left) and Cre (right) groups. Quantification of Esr1 in the MPOA.

k-m, Quantitative comparisons of female mice: first vaginal opening (VO) age (k), estrous age (l, left), receptive age (l, right), number of intromissions received (m, left), and receptivity (m, right).

n-p, Quantitative comparisons of male mice: balanopreputial separation (BPS) age (n), first mount age (o), number of mounts (p, left) and thrusts (p, right).

Line plots are shown in mean ± S.E.M and analyzed with a two-way repeated measures ANOVA followed by multiple comparisons. Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with unpaired t-test. ***p < 0.001, **p < 0.01. *p < 0.05. Statistical details in Methods.

Single cell RNA sequencing identifies neuronal cell types across adolescence in the MPOA

a, Schematic illustrating MPOA scRNAseq experiment methods.

b, Joint clustering of 24,627 neurons from all groups (P23, P35, P50, GDX). Total identified clusters = 32. UMAP plot is color coded by neuronal cluster type, listed in legend.

c, Visualization of neuron cluster transcriptional associations via PAGA analysis. Neuron cluster similarity is represented by the thickness of each connecting line. Vgat+Esr1+ clusters (Vgat 2, 4, & 16) demonstrate high connectedness and are highlighted within the dotted circle. Clusters: Pink = Vglut2+, Blue = Vgat+Esr1-, Purple = Vgat+Esr1+, Grey = Mixed.

d, Heatmap showing the scaled sum of DEG log fold changes for each Vgat+ cluster in comparison to GDX samples.

e, Linear regression analysis of percentage of Esr1 expressing cells (x-axis) and the number of P50M-enriched DEGs in comparison with GDX samples (y-axis), where each dot is a Vgat+ cluster.

f, Left: Scatter plots showing aR2 values of hormone receptor genes (dots) in females (left) and males (right) comparing the percentage of percent expressing cells to the number of DEGs in comparison with GDX samples at each Vgat+ cluster for P50-enriched (top) and P35-enriched (bottom) DEGs. Right: SCENIC analysis-computed ranked sums of TFs associated with HA-DEGs and their co-expression scores in females (top, 952 TFs) and males (bottom, 915 TFs). HA-DEGs in both sexes show highest co-expression with Esr1. Each TF is plotted along the x-axis in descending rank order.

g, Quantitative comparison of aggregate HA-DEG expression (represented by AUCell score) within each Vgat+ cluster across P23, P35, P50, and GDX groups.

h, Top: Schematic illustrating integrative analysis, establishing correspondence between MERFISH20 (Moffitt et al., 2018) and scRNAseq datasets. The MERFISH label transfer algorithm (details in Methods) identified behaviorally-relevant cells within the defined scRNAseq clusters. Bottom: The dot plot graph represents the percent of cells, indicated by dot size, identified as behaviorally-relevant (selected female behaviors: mating & parental; selected male behaviors: mating & aggression) in each Vgat+ cluster. Data integration revealed the top two marker genes related to mating as Esr1 and Apoc3 (female, left) and Esr1 and Lamp5 (male, right) (see panel i). The dot plot graph also represents the percent of cells, indicated by dot size, expressing the marker gene in each Vgat+ cluster, in addition to its z-scored gene expression indicated by dot color intensity.

i, Scatter plots of enriched genes in mating-related scRNAseq clusters in females (top) and males (bottom).

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. ***p < 0.001. Statistical details in Methods.

aR2: adjusted R squared; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Identification of adolescent transcriptional trajectories in MPOAVgat+Esr1+ neurons

a-b, UMAP visualization of Vgat+Esr1+ cells and their transcriptional trajectories depicted by a solid black line. Vgat+Esr1+ cells are color coded by group (left) and pseudotime (right), where progression of time is delineated from dark to bright coloring. Dashed arrows indicate the direction of transcriptional progression. Kinetics plots show the relative expression of HA-DEGs Pgr (shared), Tead1 (female-related), and Napb (male-related) in Vgat+Esr1+ cells for each group across adolescent pseudotime (x-axis). a: females; b: males.

c, e, Box plots showing scaled gene expression of Vgat+Esr1+ branch enriched genes in female (c) and male (e) mice.

d, f, Box plots showing Vgat+Esr1+ cell placement along pseudotime for each group.

g, i, UMAP visualization of Vgat+HRLow cells and their transcriptional trajectories depicted by a solid black line. Vgat+HRLow cells are color coded by group in females (g) and males (i).

h, j, Cumulative distributions of decoding accuracy by SVM classification between mature groups (P50, P35) and immature groups (P23, GDX) using expression data from Vgat+Esr1+ (black, real line), Vgat+ HRLow (grey) or shuffled data (dashed line) in females (h) and males (j).

k, m, Volcano plots comparing gene expression between late and early trajectories in female (k) and male (m). Sex-specific and shared gene programs are highlighted (shared: yellow; female-specific: salmon; male-specific: blue).

l, n, Box and line plots showing scaled expression of sex-specific and shared enriched genes from early to late trajectories in females (l) and males (n).

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Kruskal-Wallis H test followed by multiple comparisons test with p-values Bonferroni corrected (c-f) or Wilcoxon rank-sum test (l, n). Cumulative distributions were analyzed one-way ANOVAs followed by multiple comparisons. ***p < 0.001.

OVX: ovariectomy; ORX: orchiectomy; GDX: gonadectomy; HRLow: hormone receptor-low; SVM: support vector machine.

HM-HCR FISH reveals spatial transcriptional trajectories of MPOAVgat+Esr1+ neurons across adolescence and hormonal state

a, Schematic timeline for highly multiplexed-HCR FISH (HM-HCR FISH) experiments with tissue collected at different combinations of ages and hormone manipulation.

b, Schematic illustrating HM-HCR FISH procedure.

c, Schematic of a representative mouse coronal brain section for analyzing MPOA. A total of 41,549 cells from MPOA were analyzed.

d, e, Representative images of reconstructed cells in MPOA, color coded by scaled expression of genes at P23, GDX, P28 with hormone supplementation, and P50. d: females; e: males. Quantitative analysis of each gene is reported in Supplementary Fig. 9a-b, i-j.

f, g, Pseudotime spatial visualization of Vgat+Esr1+cells across all 6 groups in females (f) and males (g), where progression of time is delineated from dark to bright coloring. Pseudotime score was computed using HM-HCR FISH gene expression data, described in detail in Methods. Quantitative comparisons of this data is reported in Supplementary Fig. 9r-u.

T: testosterone; E: estradiol; AC: anterior commissure; MPOA: medial preoptic area; MPN: medial preoptic nucleus; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Identification of sexually dimorphic genes in MPOA

a, Heatmap representing the scaled sum of log fold changes of sexually dimorphic genes for each Vgat+ cluster comparing females and males (bidirectionally) in hormonally-intact P50 mice or GDX mice.

b, Linear regression analysis comparing the percentage of Esr1 expressing cells and the number of female dimorphic genes within each Vgat+ cluster (individual dots).

c, Scatter plots showing adjusted R2 values of hormone receptor genes (dots) in females (top) and males (bottom) as a result of linear regression analysis comparing the percentage of hormone receptor gene expressing cells and number of sexually dimorphic genes across Vgat+ clusters.

d, SCENIC analysis-computed ranked sums of TFs associated with sexually dimorphic genes and their co-expression scores in females (top) and males (bottom). Sexually dimorphic genes in both sexes show highest co-expression with Esr1. Each TF is plotted along the x-axis in descending rank order.

e, Volcano plot comparing P50M and P50F gene expression in Vgat+Esr1+ cells. Dimorphic genes are numbered and color coded (P50F-rich: salmon; P50M-rich: blue).

f, g, Scatter plots showing fold change differences of sexually dimorphic genes (f) or adolescence-related genes (g) between P50F and P23F (x-axis) and P50M and P23M (y-axis). Female-rich genes plotted on the left and male-rich genes plotted on the right. Adolescent genes and dimorphic genes were only partially overlapping. Adolescent genes were partially shared between sexes.

h, Motif-enrichment analysis of male-rich dimorphic genes reveals deconstructed Ar-regulons. The co-expression score between Ar and a regulon gene is represented by the thickness of their connecting line in the visualization (left) and via a bar graph (right).

i, Box plots comparing P50M to P50F expression of male-rich dimorphic Ar-regulon genes at each Vgat+ cluster via AUCell analysis.

j, Linear regression analysis between the percentage of Ar (left) or Esr1 (right) expressing cells and AUCell score of male-rich dimorphic Ar-regulon genes at each Vgat+ cluster (each dot).

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Wilcoxon rank-sum test. p-values were Bonferroni corrected. ***p < 0.001. Statistical details in Methods.

aR2: adjusted R squared; TF: transcription factor; FC: fold change; Ar: androgen.

Adolescent dynamics of sexually dimorphic genes in MPOAVgat+Esr1+ neurons

a, b, Bar plot of adult dimorphic genes in females (a) and males (b) in each group. Box plots show scaled gene expression of sexually dimorphic genes (left: all; middle: adolescent-unrelated; right: adolescent-related) within Vgat+Esr1+ cells in females (a) and males (b).

c, d, Cumulative distributions of decoding accuracy by SVM classification for adolescent-unrelated (left) and -related (right) dimorphic genes between females and males of matched groups using female (c) or male (d) dimorphic genes. Shuffled data depicted with dashed lines.

e, UMAP transcriptional trajectory visualization of combined female and male Vgat+Esr1+ cells (dots), color coded by group (left) or pseudotime (right) where progression of time is delineated from dark to bright coloring.

f, Box plots showing pseudotime assignment of Vgat+Esr1+ cells across groups in females (left) and males (right).

g, h, Box plots comparing AUCell aggregate expression of female dimorphic (g) or male dimorphic (h) genes (top: adolescent-unrelated; bottom: adolescent-related) across groups.

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. Cumulative line graphs were analyzed with one-way ANOVAs followed by multiple comparisons. ***p < 0.001. Statistical details in Methods.

GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy; SVM: support vector machine.

Regulation of adolescent transcriptional dynamics via Esr1 activation in MPOAVgat+Esr1+ neurons

a, Schematic illustrating cell type and site-specific deletion of Esr1 in MPOA using viral vector AAV-frtFlex-Cre in Slc32a1Flp::Esr1lox/lox mice, followed by scRNAseq.

b, Heatmaps show the scaled average expression (z-score) of Vgat+ HA-DEGs for Esr1+ and hormoneRLow cells across Esr1KO, GDX, and P50 groups for females (left) and males (right).

c, d, UMAP visualization of Vgat+Esr1+ cells and their transcriptional trajectories depicted by a solid black line in females (c) and males (d). Vgat+Esr1+ cells are color coded by group (left) and pseudotime (right), where progression of time is delineated from dark to bright coloring.

e, g, Box plots showing scaled gene expression of Vgat+Esr1+ branch enriched genes for each group in female (e) and male (g) mice.

f, h, Box plots showing Vgat+Esr1+ cell placement along pseudotime for each group.

i, k, Volcano plot comparing gene expression between late and early trajectories in female (i) and male

(k). Sex-specific and shared gene programs are highlighted (shared: yellow; female-specific: salmon; male-specific: blue).

j, l, Box and line plots showing scaled expression of sex-specific and shared enriched genes from early to late trajectories in females (j) and males (l).

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. ***p < 0.001. Statistical details in Methods.

HA-DEG: hormone-associated differentially expressed gene; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Deconstructing sex-specific gene regulatory networks underlying adolescent transcriptional dynamics in MPOAVgat+Esr1+ neurons

a, Schematic demonstrating the deconstruction of causal GRNs through the combination of functional genomics and SCENIC analysis. (Left to right) Functional-SCENIC analysis involves identification of pTFs that co-express previously identified HA-DEGs, then those pTF-DEGs are analyzed for regulons and sorted into primary and secondary cis-regulated genes. With that data, a GRN can be reconstructed for each pTF.

b, Heatmap for Vgat+Esr1+ cells scaled expression (z-score) of Esr1 and Esr1-regulated TFs in Esr1KO, P23, GDX, P35, and P50 groups for females (left) and males (right).

c, Representative images of reconstructed cells color coded by scaled expression of Creb3l1 (female, top) or Pou6f2 (male, bottom) in the MPOA at ages P23 (left) and P50 (right). Violin plots show the number of Creb3l1 (female, top) or Pou6f2 (male, bottom) transcripts per cell at ages P23 and P50.

d, e, Motif-enrichment analysis of Esr1-DEG deconstructed GRNs in females (d) and males (e). The visual representation of the GRNs show Esr1 and Esr1-regulated TFs cis-regulate Esr1-DEGs. TFs are depicted by numbered boxes and circles denote Esr1-DEGs color coded by ontology category.

f, g, Violin plots comparing Esr1KO to control Vgat+ cell gene expression (measured as number of transcripts per cell) in MPOA of females (f) and males (g).

Violin plots are outlined with a distribution line, individual dots represent each cell, and depict a red dot indicating the median. Violin plots are analyzed with Wilcoxon rank-sum test. ***p < 0.001, **p < 0.01. *p < 0.05. Statistical details in Methods.

GRN: gene regulatory network; pTF: primary transcription factor; HA-DEG: hormone-associated differentially expressed gene; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Supporting data for cell-type specific deletion of Esr1 and behavioral experiments, related to Fig. 1

a, Schematic illustrating unilateral injection of AAV-frtFlex-cre virus in the MPOA of Esr1lox/lox mice. Representative image of Esr1 expression in MPOA after a left hemisphere injection. Quantification of Esr1 in the MPOA. Note that this is a control experiment to test viral specificity, where mice did not express Flp-recombinase.

b, Schematic illustrating AAV-frtFlex-Cre injection in the MPOA of Slc32a1Flp or Slc17a6Flp mice. Representative images showing expression of GFP viral-reporter and Slc32a1 (top) or Slc17a6 (bottom) in the MPOA of Slc32a1Flp (top) or Slc17a6Flp (bottom) mice.

c-d, Quantitative comparisons of female mice: number of mounts received (c), escape attempts (d, left) and submissive postures (d, right) in Slc32a1Flp::Esr1lox/lox mice.

e-f, From left to right, females (top) to males (bottom): Quantifications of body weight, locomotion, time investigating a male conspecific, time investigating a female conspecific, social preference in three chamber test, and time spent in the open arm of EPM of Slc32a1Flp::Esr1lox/lox (e) or Slc17a6Flp::Esr1lox/lox (f) mice.

g-h, Linear regression analysis between number of Esr1 cells and mating-related behavioral measurements in female (left) and male (right) Slc32a1Flp::Esr1lox/lox (g) or Slc17a6Flp::Esr1lox/lox (h) mice.

Line plots are shown in mean ± S.E.M and analyzed with a two-way repeated measures ANOVA followed by multiple comparisons. Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with unpaired t-test. **p < 0.01. *p < 0.05. Statistical details in Methods.

aR2: adjusted R squared; EPM: elevated plus maze.

Supporting data for scRNAseq clustering, related to Fig. 2

a, Joint clustering of 58,921 cells from all groups (P23, P35, P50, GDX). UMAP plot is color coded by 9 cell category types, listed in legend.

b, Violin plots showing the gene (top) and UMI (bottom) distributions in each cell category type.

c, Violin plots of marker gene expression in each cell category type.

d, Proportion of cells represented in each cell category type from each group.

e, Graph representation of cluster number(s) from MPOA cell sub-sampling.

f, Violin plots showing the gene (top) and UMI (bottom) distributions in each neuronal cell type cluster.

g, Graph representation of cluster number(s) from MPOA neuron sub-sampling.

h, Proportion of cells represented in each neuron cluster from each group. All neuronal cell types were present irrespective of age, sex, and hormonal states.

i, UMAP clustering without integration (left) and then joint clustering (right) of neuronal cells from all groups. UMAP plot is color coded by experimental group. Boxplot of LISI (Local Inverse Simpson’s Index, or cell “mixing” index) was computed to assess integration performance. Note that the theoretical maximum of LISI score is 8.

j, UMAP plot with color code for major neuronal types: Salmon = Vglut2+, Blue = Vgat+Esr1-, Purple = Vgat+Esr1+, Grey = Mixed.

k, Cluster tree illustrating the lineage of clusters at each clustering resolution. Esr1 expression and the number of cells in each cluster are represented by the color and the size of the dot, respectively. Number of cells diverging from a node is represented by the color of the connecting line.

l, Violin plots showing normalized expression values of Slc32a1, Slc17a6, several steroid hormone receptor genes, and canonical marker genes in the MPOA at each cluster.

m, In-silico analysis generated UMAP plots, using simulated subsets of data. UMAP plots are color coded by neuronal types. Vgat+Esr1+ clusters are highlighted.

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed with Wilcoxon rank-sum test. ***p < 0.001. Statistical details in Methods.

UMIs: unique molecular identifiers; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Supporting data for scRNAseq experiments, related to Fig. 2

a, Dot plot illustrating scaled expression (z-score, color intensity) and percentage of expressing cells (dot size) for a list of marker genes (y-axis) in each cluster (x-axis). On the left y-axis each clusters gene expression score is presented for each group: P23, P35, P50, and GDX (females and males).

b, c, Bar plot to show the difference in number of DEGs between P50 and GDX groups (females: salmon; males: blue) in each Vgat+ cluster (b) and Vglut2+ cluster (c).

d, e, Heatmaps with hierarchical clustering dendrograms (top) and line graphs (bottom) to demonstrate gene expression patterns within a cluster relative to other clusters for each group: P23, P35, P50, GDX in females (left) and males (right). Clusters are denoted by color. Cluster numbers that are bolded and marked with a hashtag (#) show higher expression at P35 and P50 than P23 and GDX. Z-scored gene expression is denoted by a separate color intensity scale. Line graph data shows normalized expression of each cluster relative to P23 expression. Vgat+ clusters (d); Vglut2+ clusters (e).

f, Linear regression analysis between percentage of Esr1-expressing cells and number of P50M-enriched DEGs compared to GDX for each Vglut2+ cluster (dot).

g, Scatter plots showing aR2 values of hormone receptor genes (dots) in females (left) and males (right) comparing the percentage of percent expressing cells to the number of DEGs in comparison with GDX samples at each Vglut2+ cluster for P50-enriched (top) and P35-enriched (bottom) DEGs.

h, Violin plots showing the normalized expression values of Slc32a1, Slc17a6, Esr1, and several subtype-specific genes in the MPOA for each Vgat+Esr1+ cluster (Vgat 2, 4, 16).

i, Heatmap illustrating the Pearson correlation coefficient between Vgat+ clusters including or excluding Esr1 expression data. Near identical clusters were detected irrespective of Esr1 expression.

j, Bar plot showing the percent of identical cell barcodes within each Vgat+ cluster comparing inclusion and exclusion of Esr1 expression data. Near identical clusters were detected irrespective of Esr1 expression.

Line graphs include standard error and were analyzed using one-way repeated-measures ANOVA followed by multiple comparisons. Statistical details in Methods.

GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy; aR2: adjusted R squared; DEG: differentially expressed gene.

Supporting data for scRNAseq and MERFISH analysis, related to Fig. 2

a, b, Percentage of Fos positive cells associated with female parental and mating behaviors (a) and male aggression and mating behaviors (b) from MERFISH cell clusters (Moffitt et al., 2018). Dashed lines indicate the overall mean percentage. Each cluster was compared to the mean using a Fisher’s exact test. **p < 0.01, ***p < 0.001. Statistical details in Methods.

c, Alluvial plots showing the correspondence between MERFISH and scRNAseq clusters color coded by behavior.

d, Heatmaps and dot plots showing both Pearson correlation coefficients of MERFISH + scRNAseq GABAergic (mGad; Vgat) clusters and their expression of selected genes: Slc32a1, Slc17a6, Esr1, Ar, or Gad1 in the scRNAseq dataset (top) and the MERFISH dataset (right). Increased color intensity corresponds to greater values. Dot size represents percentage of cells expressing the corresponding gene. Left: females; right: males.

Supporting data for transcriptional trajectory branches from scRNAseq data, related to Fig. 3

a, b, UMAP visualization of Vgat+Esr1+ cells and their transcriptional trajectories depicted by a solid black line. Vgat+Esr1+ cells are color coded by group (left) and pseudotime (right), where progression of time is delineated from dark to bright coloring. Dashed arrows indicate the direction of transcriptional progression. Smaller UMAPs on the right visually isolate the trunk and two identified branch trajectories. a: females; b: males.

c, d, Box plots show scaled gene expression in Vgat+Esr1+ cells from Branch 1 (left column) and Branch 2 (right column) trajectories for Branch 1-specific genes (top row), Branch 2-specific genes (middle row), and Branch-common genes (bottom row) across all groups.

e, f, Box plots show pseudotime score for Branch 1 (left column) and Branch 2 (right column) trajectories across groups in females (e) and males (f).

g, h, Bar plots represent the percentage of Vgat+Esr1+ cells that belong to cluster Vgat 2, 4, and 16 in the trunk and branches. Cumulative distributions for decoding accuracy between more mature (P35, P50) and immature (P23, GDX) groups using Branch 1 (yellow) or Branch 2 (grey) -specific genes compared to shuffled data (dashed line). g: females; h: males.

i, j, Box plots of HA-DEG aggregate expression (AUCell score) separated by branch in P35 and P50. Cumulative distributions for decoding accuracy between P35 and using Branch 1 (yellow) or Branch 2 (grey) -specific genes compared to shuffled data (dashed line). i: females; j: males.

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed using Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. Cumulative distributions were analyzed with unpaired t-tests for each branch. **p < 0.01, ***p < 0.001. Statistical details in Methods.

HA-DEG: hormone-associated differentially expressed gene; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Supporting data for scRNAseq trajectory analysis, related to Fig. 3

a, Scatter plot showing the percentage of Esr1 and/or Ar expressing cells in each Vgat+ cluster (dot).

b,c, UMAP visualization of transcriptional trajectories (black line) color coded by experimental group in Vgat+hormoneRLow (left), Vgat+Esr1-Ar+ (middle), and Vgat+Esr1+ (right) in females (b) and males (c).

d, Cumulative distributions of decoding accuracy by SVM classification between mature groups (P50, P35) and immature groups (P23, GDX) using expression data from Vgat+Esr1-Ar+ (salmon: female, blue: male, solid line), or shuffled data (dashed line).

e, UMAP visualization of Vgat+Esr1+ cells combining males and females, depicting their transcriptional trajectory denoted by a solid black line. Vgat+Esr1+ cells are color coded by group (right) and pseudotime (left), where progression of time is delineated from dark to bright coloring.

f, Box plot of pseudotime score from combined male and female UMAP (e) separated by group.

g, Bar plot showing the numbers of sex-specific and shared genes enriched in later trajectory, which were identified using a separate manifold for each sex or a joint manifold. Majority of sex-specific gene programs were not identified in joint manifold analysis.

h,i, PHATE visualization of transcriptional progression in Vgat+Esr1+ cells (dot) color coded by group (left) and relative mature likelihood (right) in females (h) and males (i).

j, l, Box plot showing MELD score (mature likelihood) of Vgat+Esr1+ cells in females (j) and males (l). k,m, Cumulative distributions of decoding accuracy between mature groups (P50, P35) and immature groups (P23, GDX) using highest expressing (yellow, used in the HCR experiments) or lowest expressing (grey) HA-DEGs in females (k) and males (m).

n,p, In silico simulation using a 50% subset of the data. UMAP visualization of transcriptional trajectory (black line) and cells (dots) color coded by group (right) and pseudotime (left) in female (n) and male (p) Vgat+Esr1+ cells. Dashed arrows indicate the direction of transcriptional progression.

o, q, In silico simulation analysis using a 50% subset of the data. Box plots show pseudotime values assigned to each Vgat+Esr1+ cell across groups in females (o) and males (q).

o, UMAP visualization of transcriptional trajectory (black line) and cells (dots) color-coded by group (right) or pseudotime (left) in Vgat+Esr1+populations, including both sexes.

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed using Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. Cumulative distributions were analyzed with unpaired t-tests. ***p < 0.001. Statistical details in Methods.

SVM: support vector machine.

Supporting information for trajectory analysis on scRNAseq data, related to Fig. 3

a, Violin plots show normalized expression values of Slc32a1, Slc17a6, several steroid hormone receptor genes, and canonical marker genes in the MPOA at each Vglut2+ cluster. Vglut2+Esr1+ clusters are highlighted.

b, c, UMAP visualization of transcriptional trajectory (black line) and cells (dots) color coded by pseudotime (left) or group (right) in female (b) and male (c) Vglut2+Esr1+ populations.

d, e, Box plot showing pseudotime of Vglut2+Esr1+ cells in females (d) and males (e).

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed using Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. ***p < 0.001. Statistical details in Methods.

Supporting data for HM-HCR FISH, related to Fig. 4

a, b, Representative images showing detected genes in 4 iterative rounds of HM-HCR FISH from the MPOA at P50 (left) and P23 (right). a: males; b: females. Scale bars: 200 µm.

Supporting data for HM-HCR FISH experiments related to Fig. 4.

a, b, Violin plots showing the expression values of individual genes across experimental groups in MPOA Vgat+Esr1+ cells. a: females; b: males.

c, f, Heatmaps showing the percentage of Vgat+Esr1+ cells (top) and scaled gene expression (bottom) across pseudotime. c: females; f: males.

d, g, Cumulative distributions of decoding accuracy between mature groups (P50, P35, hormone-supplemented P28) and immature groups (P23, GDX, P28 control) using expression data from Vgat+Esr1+ (black, solid line), Vgat+hormoneRLow (grey, solid line), and shuffled data (dashed lines). d: females; g: males.

e, h, Cumulative distributions of decoding accuracy between mature (P50, P35, hormone-supplemented P28) and immature (P23, GDX, P28 control) groups using subsets of Vgat+Esr1+ gene expression data, color coded by number of genes used for decoding. e: females; h: males.

i, j, Correlograms of Pearson correlation coefficients between all pairs of DEGs in MPOA Vgat+Esr1+ cells. i: females; j: males.

k, Schematic illustrating integrative analysis between HM-HCR FISH and scRNAseq datasets to predict scRNAseq gene expression.

l, o, Scatter plots showing the correlated expressions of real and predicted scRNAseq data. l: females; o: males.

m, p, Heatmaps showing the percentage of Vgat+Esr1+ cells (top) and scaled predicted gene expression (bottom) across pseudotime. m: females; p: males.

n, p, Cumulative distributions of decoding accuracy between mature groups (P50, P35) and immature groups (P23, GDX) using predicted expression data in females (n) and males (p).

r-u, Quantitative analyses of pseudotime in space, corresponding to Fig. 4f, g.

r, s, Boxplots showing pseudotime of Vgat+Esr1+ cells in the medial and lateral MPOA of females (r) and males (s).

t, Schematic illustrating the quantification of sex differences across the transcriptional progression rate. Neighborhood analysis was performed in physical space to compute the maturation index of cells from P28 control groups.

u, Box plot of maturation index data from Vgat+Esr1+ cells in the medial (M) and lateral (L) MPOA of P28 female and male controls. At P28, male mice showed an early onset of transcriptional progression in the lateral part of the MPOA.

Violin plots are outlined with a distribution line, individual dots represent each cell, and red dot indicates the median. Violin plots are analyzed with Wilcoxon rank-sum test. Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed using Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. Cumulative distributions were analyzed with unpaired t-tests. ***p < 0.001. Statistical details in Methods.

R: Pearson correlation coefficient; M: medial MPOA; L: lateral MPOA; P28FC: P28 female control; P28FE: P28 females with estrogen treatment; P28MC: P28 male control; P28MT: P28 males with testosterone treatment; GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy.

Supporting data for sexual dimorphism analysis, related to Fig. 5-7

a, Bar plot showing the number of sexually dimorphic genes at each Vgat+ cluster comparing female (salmon) and male (blue) groups at P50.

b, Volcano plot comparing P50 female (salmon) and male (blue) gene expression in Vgat+Esr1+ cells (left). Dimorphic genes are highlighted and indicated numerically. This is the same panel shown in Fig. 5e and reused here to facilitate the comparison with GDX data shown on the right. Volcano plot comparing OVX (salmon) and ORX (blue) gene expression in Vgat+Esr1+ cells. Dimorphic genes are highlighted and indicated numerically.

c, Volcano plots comparing intact P50 to Esr1KO P50 gene expression in Vgat+Esr1+ cells. The percent of sexually dimorphic genes (P50F-enriched genes: salmon, P50M-enriched genes: blue) present in intact P50 and Esr1KOs are quantified in the bar plot (top). Left: females; right: males.

d, Cumulative distributions of decoding accuracy between males and females (P50: black, P35: dark grey, P23: light grey, shuffled: dashed).

Cumulative distributions were analyzed with one-way ANOVA followed by multiple comparisons. Details in Methods. ***p < 0.001. Statistical details in Methods.

GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy; LogFC: Log fold change.

Supporting data for trajectory and gene regulatory network analysis followed by cell type specific deletion of Esr1, related to Fig. 7

a, b, UMAP visualization of neuronal clusters in Esr1KOF (a, left) and Esr1KOM (b, left). Feature plots showing normalized expression of Slc32a1 and Slc17a6 in Esr1KOF (a, right) and Esr1KOM (b, right). c, d, Violin plots showing gene (top) and UMI (middle) distributions in each neuronal cell type. Bar graphs show the fraction representation of each neuronal cluster in Esr1KOF (c, bottom) and Esr1KOM (d, bottom).

e, f, Heatmaps illustrating Pearson correlation coefficient (red) between Vgat+ clusters of Esr1KO and intact groups. Dot plots are attached illustrating scaled expression (color intensity) and the percentage of expressing cells (dot size) of Slc32a1, Slc17a6, Esr1, and Ar in intact (top) and Esr1KOs (right side) in females (e) and males (f).

g, Schematic illustrating integrative analysis to establish correspondence between cell types of Esr1KO data and the reference scRNAseq data.

h, i, Bar graphs showing the correspondence between Esr1KO and the reference Vgat+Esr1+ clusters in females (h) and males (i) using the label transfer technique (details in Methods).

Violin plots are outlined with a distribution line and red dot indicates the median.

Supporting data for trajectory analysis followed by cell type specific deletion of Esr1, related to Fig. 7

a,b, UMAP visualization of transcriptional trajectories (black line) and cells (dots) color coded by experimental group (left) and pseudotime (inset) in Vgat+Esr1+ cells. UMAPs on the right visualize where cells from the branches and trunk lie on the trajectories. a: females; b: males.

c, d, Box plots show scaled gene expression in Vgat+Esr1+ cells from Branch 1 (left column) and Branch 2 (right column) trajectories for Branch 1-specific genes (top row), Branch 2-specific genes (middle row), and Branch-common genes (bottom row) across all groups. c: females; d: males.

e, f, Box plots show pseudotime score for Branch 1 (left column) and Branch 2 (right column) trajectories across groups in females (e) and males (f).

g,h, Bar plot of Vgat+Esr1+ cluster (eVgat for Esr1KOs) proportion of cells that reside in the trunk, branch1, or branch2. g: females; h: males.

Box plots are shown with box (25%, median line, and 75%) and whiskers and analyzed using Kruskal-Wallis H test followed by multiple comparisons test. p-values were Bonferroni corrected. ***p < 0.001, **p < 0.01. Statistical details in Methods.

GDX: gonadectomy; OVX: ovariectomy; ORX: orchiectomy; KO: knockout.

Supporting information for gene regulatory network analysis followed by cell-type specific deletion of Esr1, related to Fig. 8

a, b, Enriched Esr1-DEG gene ontology terms are placed alongside bar plots corresponding to gene expression percentage. Top 25 GO terms are shown and color-coded by gene category. Bars are fractionated into direct Esr1-regulated genes (yellow), secondary TF-regulated genes (black), and other genes regulated by Esr1-TFs and other TFs (grey). GO Biological Process, GO Molecular Function, and GO Cellular Component are referenced. a: females; b: males.