Adolescence is a secondary critical period, following the perinatal steroid surge, during which juvenile animals undergo extensive physiological maturation of both the body and nervous system as they transition into adulthood ^1,2. This maturation is primarily orchestrated by the hypothalamus-pituitary-gonad (HPG) axis, which activates gonadal cells in the testes and ovaries to release testosterone and estrogen, initiated during adolescence ^1,3,4. These circulating sex steroids are essential not only for reproductive maturation but also for guiding the development and execution of sex-specific behaviors through their direct actions on the central nervous system ^5–7. Despite the well-established role of these hormones in shaping behavior, the molecular mechanisms underlying their influence on brain development during adolescence remain poorly understood.

The medial preoptic area (MPOA) is a critical region in the brain, known for its involvement in mating and other social behaviors ^7–14. Within the MPOA, specific molecularly defined subpopulations of neurons express various neuropeptides and hormone receptors that are closely tied to reproductive behavior ^14–18. Notably, the MPOA is enriched with the expression of steroid hormone receptor genes such as Esr1, or estrogen receptor 1 ^19,20, leading us to hypothesize that gonadal steroid signaling during adolescence facilitates the transcriptional maturation of these genetically defined MPOA neurons, which are crucial for the development of mating behaviors. While recent advances have shed light on the molecular diversity of neurons in sexually dimorphic brain regions during the perinatal period ²¹ and adulthood ^20,22,23, the transcriptional dynamics that occur during adolescence, particularly at the single-cell level, remain largely unexplored. This gap in our understanding represents a critical frontier in unraveling how hormonal signaling shapes the neural circuits that govern sex-specific behaviors.

In this study, we used single-cell RNA sequencing (scRNAseq) and hybridized chain reaction fluorescent in situ hybridization (HCR-FISH) to investigate the transcriptional dynamics of the MPOA throughout adolescence in both male and female mice. We show that Esr1 plays a critical role in the transcriptional maturation of GABAergic neurons (Vgat+) within the MPOA, which is essential for the development of mating behaviors. Deleting Esr1 in Vgat+ neurons disrupted the development of mating behaviors in both sexes, while deletion in glutamatergic neurons (Vglut2+) had no observable effect. Our findings highlight the importance of Esr1 in guiding sex-specific transcriptional programs during adolescence, revealing distinct hormone-dependent gene expression patterns and gene-regulatory networks crucial for reproductive behaviors.

Results

Selective knockout of Esr1 in GABAergic or glutamatergic MPOA neurons at pre-adolescence

Previous studies have shown that knocking out Esr1, either globally or in specific cell types, reduces mating behavior ^7,24. However, the role of Esr1 in genetically defined MPOA neurons during adolescent development remains unclear, as these studies could not easily differentiate between effects occurring during perinatal and adolescent periods. To address this gap and directly test whether the emergence of mating behavior requires steroid hormone signaling in molecularly defined MPOA cell types, we developed an AAV-frtFlex-Cre virus. We injected this virus into the MPOA of Slc32a1^Flp::Esr1^lox/lox (MPOA^Vgat-Esr1KO) or Slc17a6^Flp::Esr1^lox/lox (MPOA^{Vglut2-Esr1KO}) male and female mice at postnatal day (P) 14-18. This viral strategy allowed for Esr1 deletion in MPOA GABAergic or glutamatergic neurons, depending on the transgenic mouse line, throughout adolescence (Fig. 1a-b,i-j and Supplementary Fig. 1a-b). The onset of genital changes, the maturation of the HPG axis (as determined by the first day of estrus in female mice²⁵ and first day of initiation of sexual behaviors in male mice²⁶), and locomotor behaviors were not altered in MPOA^Vgat-Esr1KO and MPOA^{Vglut2-Esr1KO} mice (Fig. 1c-d,f,k-l,n and Supplementary Fig. 1e-f). The establishment of mating behavior was abolished in MPOA^Vgat-Esr1KO mice in both sexes (Fig. 1d-e,g-h and Supplementary Fig. 1a-b,g), but unaffected in MPOA^{Vglut2-Esr1KO} mice (Fig. 1l-m,o-p and Supplementary Fig. 1h). These results extend prior work to show that Esr1 signaling in MPOA GABAergic neurons is required for the emergence of mating behavior during adolescence ^7,12,13.

Fig. 1:

Transcriptional dynamics of MPOA cell types during adolescent development

To characterize the transcriptional dynamics occurring in the MPOA specific to the adolescent period, we collected tissue from male and female wildtype mice at pre- (P23), mid- (P35) and post- (P50) adolescence timepoints, as well as from mice gonadectomized (GDX) prior to adolescence onset. We used droplet-based scRNAseq to recover the transcriptomes of 58,921 cells and combined the data from all samples across the 8 experimental groups (Fig. 2a and Supplementary Fig. 2a-e and Supplementary Table 1). Using known neuronal marker genes Stmn2 and Thy1, we performed high-resolution clustering to focus our subsequent analyses on neuronal cell data (24,627 cells). Thirty-five neuron specific clusters were resolved with 434 marker genes in total (UMIs per cell: 5,395; median genes per cell: 2,592), resulting in 20 GABAergic clusters (Vgat+ clusters, 13,334 cells; 54.4% of neurons), 12 glutamatergic clusters (Vglut2+ clusters, 7,658 cells; 31.1% of neurons), and 3 clusters with mixed glutamatergic and GABAergic markers (3,635 cells; 14.8% of neurons) (Fig. 2b and Supplementary Figs. 2f-l, 3a,j and Supplementary Table 2), largely consistent with previously published data ^20,27. To further quantify transcriptional associations between MPOA neuronal cell types, we performed partition-based graph abstraction (PAGA) analysis ²⁸ to identify cell cluster similarities through a graphical representation of cell differentiation (Fig. 2c). Consistent with their distinct roles, Vgat+ and Vglut2+ neurons were readily separated in the PAGA graph. Additionally, Vgat+ clusters 2, 4 and 16, all of which expressed Esr1, showed higher connectivity compared to other Vgat+ cell type clusters (Fig. 2c and Supplementary Figs. 2m, 3h-i), suggesting that Vgat+Esr1+ clusters form a distinct transcriptional unit.

Fig. 2:

Irrespective of age, sex, and hormonal states, we found that each identified cell type was represented across all 8 experimental groups (Supplementary Figs. 2d,h-i, 3a), indicating that new cell types do not emerge during adolescence. However, when assessing gene expression differences across timepoints and sex, we found that hormone-associated differentially expressed genes (HA-DEGs) showed substantial variability across neuronal cell types, ranging from 0 to 88 HA-DEGs per cell type (Fig. 2d-f and Supplementary Fig. 3a-g and Supplementary Tables 3-5). We then used these HA-DEGs in regression, co-expression, and AUCell analyses to identify the neuronal types that are transcriptionally sensitive at different hormonal states ²⁹. Regression analysis revealed a strong enrichment of HA-DEGs in Vgat+Esr1+ neurons (Fig. 2e-f). Consistent with this result, co-expression analysis demonstrated that out of more than 1,000 transcription factors (TFs), Esr1 was one of the most co-expressed in both sexes (Fig. 2f and Supplementary Table 6). Furthermore, in predominantly Vgat+Esr1+ clusters, the aggregate expression of HA-DEGs (AUCell) displayed patterns sensitive to both adolescence and hormonal changes, with the highest expression observed in adult mice (P50) and notably lower expression at P23 and during hormonal deprivation (GDX) (Fig. 2g). Together, these analyses demonstrate that sex-hormone signaling during adolescence development significantly influences the transcriptional states of specific MPOA neuronal cell types, particularly cells belonging to Vgat+Esr1+ clusters.

Having identified the distinct transcriptional sensitivity of Vgat+Esr1+ neurons to hormone signaling during adolescence, we next examined whether these neurons might play a direct role during behavior. Given that mating behavior was abolished in MPOA^Vgat-Esr1KO mice in both sexes (Fig. 1), this finding, along with our scRNAseq analyses, suggests that Vgat+Esr1+ neurons are critically involved during mating behaviors. To test this hypothesis, we conducted an integrative analysis using our scRNAseq dataset and a publicly available multiplexed error-robust fluorescence in situ hybridization (MERFISH) dataset ²⁰ which identified social behavior-activated POA cells via mating, parenting, and fighting behavior assays using Fos (Fig. 2h). We clustered the neuronal cells in the MERFISH data and computed the enrichment of Fos in the derived clusters. Consistent with previous observations, only a subset of MERFISH Gad1+ inhibitory clusters showed Fos enrichment (Supplementary Fig. 4a,b). After integrating our scRNAseq clusters with the MERFISH clusters (Details in Methods. Seurat V3 Integrative Analysis: reference-based integration; label transfers ³⁰) we were able to identify mating, parenting, and fighting-relevant neuronal clusters within our dataset (Fig. 2h and Supplementary Fig. 4c-d and Supplementary Table 7). We found that Esr1 was the most enriched and selective gene in the mating-relevant scRNAseq Vgat+ clusters (Fig. 2h-i and Supplementary Table 7). This observation is consistent with a previous study demonstrating that optogenetic stimulation of Vgat+Esr1+ neurons in the MPOA is sufficient to induce sexual behaviors in male mice ¹⁶. Collectively, these findings suggest that a subset of Vgat+ clusters expressing Esr1, which exhibit transcriptional changes during adolescence, are also functionally linked to mating behaviors in both males and females.

While these differential gene expression (DE) analyses quantify changes in individual gene expression at the pseudo-bulk level (Fig. 2d and Supplementary Fig. 3b-c), it does not capture the transitions in transcriptional states of individual cells as they progress through distinct biological stages. These transitions, driven by complex interactions between multiple genes, are key to understanding how single cells change over time ³¹. To resolve these dynamic changes, we applied pseudotime analysis (Fig. 3 and Supplementary Figs. 5-6), a method that infers the transcriptional progression of individual cells through a biological process, such as adolescence, by constructing a principal graph based on combinatorial gene expression patterns ^32,33.

Fig. 3:

Recognizing the potential for sex-specific differences in transcriptional progression, we created separate pseudotime manifold models for males and females to capture nuanced dynamics during adolescence ³⁴. In both sexes, Vgat+Esr1+ trajectories revealed that transcriptional states at P35 and P50 branch and diverge from preadolescence (Fig. 3a-f and Supplementary Fig. 5a-f), suggesting an acceleration of transcriptional dynamics between these timepoints. Neurons from GDX mice, on the other hand, showed arrested transcriptional progression, closely resembling preadolescent states, demonstrating the importance of circulating sex hormones in the maturation of these cells during adolescence. In contrast, Vgat+ cells lacking steroid hormone receptor gene expression (Vgat+ hormone R^Low) exhibited minimal transcriptional changes in response to hormone changes (Fig. 3g,i and Supplementary Fig. 6a-c). Furthermore, all Vgat+Esr1+ clusters were found to co-express the androgen receptor gene (Ar) to a high degree (68.3 – 88.3 %, Supplementary Fig. 6a). However, Vgat+ clusters expressing Ar but not Esr1 (Vgat+ Esr1-AR+) showed only a moderate transcriptional progression in males and minimal changes in females (Supplementary Fig. 6b-d), suggesting that Esr1 is more influential in driving transcriptional changes during adolescence than Ar alone. Lastly, we analyzed Vglut2+ clusters enriched with Esr1, and again found only moderate transcriptional progression in both sexes (Supplementary Fig. 7), altogether indicating that cell types outside of Vgat+Esr1+ clusters are less transcriptionally dynamic during this period.

Two major branches of transcriptional trajectories emerged from subpopulations of Vgat+Esr1+ cells in both sexes. In females, Branch1 was dominated by Vgat+ cluster 4, while Branch2 was occupied by clusters 2 and 16. In males, Branch1 was composed of Vgat+ clusters 4 and 16, while Branch 2 was primarily formed by Vgat+ cluster 2 (Supplementary Fig. 5g-j). The analysis of DEGs along these branches revealed distinct gene enrichment, indicating that specific subpopulations of Vgat+Esr1+ cells undergo unique transcriptional progressions during adolescence (Supplementary Fig. 5c-d and Supplementary Table 8). When we combined the trajectories for both sexes into a joint manifold, the model failed to capture hormone-dependent dynamics (Supplementary Fig. 6e-g), emphasizing that separate models are necessary for accurate analysis of sex-specific and sex-shared transcriptional progressions (Fig. 3k-n).

To strengthen our understanding of the developmental trajectory and transcriptional maturity of Vgat+Esr1+ clusters, we employed a support vector machine (SVM) classifier to predict the developmental state of single cells based solely on HA-DEGs ³⁵. The SVM accurately classified Vgat+Esr1+ single cells as either transcriptionally mature (intact adolescent or adult) or immature (preadolescent or hormonally deprived) with high accuracy (median prediction accuracy 90.6 – 93.6 %). In contrast, Vgat+ hormone R^Low cells had significantly lower prediction accuracy (median 66.7 - 73.0 %) (Fig. 3h,j), further underscoring the importance of Vgat+Esr1+ clusters in the development of MPOA transcriptional states.

Finally, we validated our pseudotime trajectories (Monocle V3) using Manifold Enhancement Latent Dimension (MELD) analysis, which measures continuous transcriptional progression in response to experimental conditions (P23, P35, P50, GDX) ^36,37. MELD analysis corroborated our findings from pseudotime, showing that intact adolescents and adults were spatially segregated in transcriptional space and exhibited a higher likelihood of reaching mature transcriptional states compared to preadolescent and hormonally deprived groups (Supplementary Fig. 6h-m). These results, visualized in PHATE (Potential of Heat-diffusion for Affinity-based Trajectory Embedding) space ³⁶, further highlight that while MPOA neuronal cell types are terminally diversified by preadolescence (P23) (Supplementary Figs. 2h, 3e), Vgat+Esr1+ neurons continue to undergo hormonally dependent transcriptional refinement from adolescence into adulthood.

Spatial phenotyping of MPOA cell types during adolescence

To assess whether adolescent hormonal changes affect spatially resolved gene expression in the MPOA and to cross-validate our scRNAseq trajectory analyses, we performed highly multiplexed hybridization chain reaction fluorescent in situ hybridization (HM-HCR FISH, V3)^38,39. This technique allowed us to detect transcripts of ∼12 genes across 41,549 MPOA cells at single-molecule resolution (Fig. 4a-c and Supplementary Figs. 8-9). The rationale for selecting this gene panel is related to Supplementary Fig. 6k,m, with details provided in the Methods. In addition to the experimental groups used in the scRNAseq experiments, we prepared tissue from hormonally supplemented mice. These mice received testosterone (males) or estrogen (females) from P23-P27 (see timeline in Fig. 4a) to investigate whether early sex steroid supplementation accelerated adolescent transcriptional trajectories.

Fig. 4:

To distinguish Vgat+Esr1+ from Vgat+ hormone R^Low cells, we measured Slc32a1, Esr1, Ar, and 8-9 DEGs during adolescence (Fig. 4d-e and Supplementary Figs. 8-9. Details in Methods). Consistent with the scRNAseq data, we observed a transcriptional progression from the preadolescent state (P23) to a matured state (P35 and P50), where the dynamics were bidirectionally influenced by circulating steroids– accelerated by early testosterone or estrogen supplementation (P28 T or E2) and delayed in a hormonally-deprived state (GDX) (Fig. 4f-g and Supplementary Fig. 9a-c,f,i-q). Interestingly, the transcriptional trajectories of female and male adolescents differed spatially, but both were regulated by sex hormones (Supplementary Fig. 9r-u). This suggests that neurons defined transcriptionally as adult female or male occupy partially overlapping spatial distributions in the MPOA. SVM classification revealed that the expression of ≥10 genes in Vgat+Esr1+ cells from the HM-HCR FISH data was sufficient to accurately classify individual cells as transcriptionally mature or immature, with high accuracy (84.4– 85.5 %). However, the same set of genes significantly underperformed in classifying Vgat+ hormone R^Low cells (63.6–64.6 %) (Supplementary Fig. 9d,g). Classification accuracy further decreased when individual genes were iteratively removed (Supplementary Fig. 9e,h), further indicating that a relatively small combinatorial set of HA-DEGs is sufficient for accurately identifying age- and sex-specific transcriptional states.

Molecular profiling of the MPOA in adult mice has established the presence of sexual dimorphism ^40,41. Our scRNAseq and in situ analyses further show hormone-dependent transcriptional progression during adolescence in both females and males. However, critical questions remain about whether (1) sexual dimorphism is evident within specific MPOA cell types, (2) sexually dimorphic genes overlap with adolescent gene sets, and (3) the degree of sexual dimorphism shifts during adolescence. To address these questions, we first identified sexually dimorphic genes across Vgat+ clusters (Supplementary Fig. 10a and 11). Like HA-DEGs, sexually dimorphic genes were enriched within Vgat+Esr1+ clusters. The number of these dimorphic genes correlated with Esr1 expression, where sexually dimorphic genes co-expressed most often with Esr1 in both sexes (Fig. 5a-e and Supplementary Fig. 10a and Supplementary Table 6). Notably, only subsets of sexually dimorphic genes were also classified as adolescent genes (Figs. 5f, 6a-b) and some of these adolescent genes were shared across sexes (Fig. 5g).

Fig. 5:

Fig. 6:

Previous research has shown that Esr1 plays a role in regulating sexual dimorphisms within sexually dimorphic brain nuclei ²¹. However, androgen receptor (AR) activation may also be essential for the expression of sexually dimorphic genes. To determine whether AR regulates male dimorphic genes, we performed single-cell regulatory network inference (SCENIC) analysis ²⁹. SCENIC identified a regulon of 33 AR-regulated genes, which were highly co-expressed with Ar and enriched with Ar’s consensus DNA regulatory element within their gene loci, in the male dimorphic gene set (Fig. 5h). AUCell analysis, however, showed that male dimorphic AR-regulon genes were enriched within male Vgat+Esr1+ clusters but not in Vgat+ Esr1-Ar+ clusters (specifically Vgat+ 5, 8, and 13) (Fig. 5i). Consistent with this observation, regression analysis showed that Esr1 expression, rather than Ar, predicted the expression of male dimorphic AR-regulon genes (Fig. 5j).

As previously noted, a subset of sexually dimorphic genes overlaps with adolescent-related genes (Fig. 5f). The number of sexually dimorphic genes increases during adolescence but reverts following gonadectomy (Fig. 6a-b). Among 45 female dimorphic genes, 12 showed increased expression during adolescence (adolescence-related dimorphic genes), while 7 were not associated with adolescence (Fig. 6a, right). Of the 110 dimorphic genes in males, 52 were adolescence-related dimorphic genes and 10 were not (Fig. 6b, right). The expression levels of adolescence-related dimorphic genes were the highest at P50 and were again reduced by gonadectomy in both sexes (Fig. 6a-b). SVM analysis revealed that adolescence-related dimorphic genes decoded sex with the highest accuracy in intact P50 males and females, followed by P35, and were less accurate in P23 and GDX groups. In contrast, adolescence-unrelated dimorphic genes successfully predict the sexes of single cells with high accuracy (>90%) irrespective of age and conditions (Fig. 6c-d), emphasizing the influence of adolescence-related dimorphic genes in shaping MPOA sex-specific transcriptional profiles. AUCell and trajectory analyses were then performed to quantify the transcriptional dynamics of sexual dimorphic genes during adolescence. These analyses consistently demonstrated that transcriptional states were the most sexually dimorphic at P50 driven by adolescence-related dimorphic genes and were the least dimorphic at P23 and GDX (Fig. 6e-h). Thus, combined trajectory analysis of all Vgat+Esr1+ neurons from both sexes, along with AUCell analysis, revealed that P50F and P50M cells were the most sexually dimorphic from each other. However, these transcriptional states were bridged via a common immature state observed largely prior to adolescence onset or in hormonally depleted conditions, suggesting that sexually dimorphic MPOA gene expression is bimodal in adults but largely continuous between males and females prior to adolescence. These analyses (Figs. 5 and 6) indicate that: 1) sexual dimorphism is most pronounced in Vgat+Esr1+ cells; 2) a subset of sexually dimorphic genes are linked to adolescence, with a fraction of adolescence-related genes shared between sexes; and 3) sexual dimorphism in the MPOA expands during adolescence as sex steroid hormone levels rise.

Esr1 knockout at preadolescence alters the transcriptional dynamics of Vgat+ MPOA neurons

Co-expression of Esr1 with adolescence genes and dimorphic genes, along with partially distinct transcriptional dynamics in Vgat+Esr1+ cells between adolescent females and males, suggests that Esr1 uniquely regulates gene expression in each sex. To further test this, we virally deleted Esr1 from Vgat+ MPOA neurons in female and male mice prior to the onset of adolescence and conducted scRNAseq at P50, comparing Esr1KO to Esr1-intact controls (Fig. 7a and Supplementary Fig. 11. Details in Methods). Esr1KO in Vgat+ MPOA neurons reduced the DEGs in males and females (Esr1-DEGs) to 600 and 824, respectively, while only reducing 5.3% and 3.3% of male and female Esr1-DEGs in Vgat+ hormone R^Low cells (Fig. 7b and Supplementary Fig. 10c and Supplementary Table 9). Previously identified HA-DEGs (Fig. 2) were largely represented in this list of Esr1-DEGs, suggesting that Esr1 deletion recapitulates hormone deprivation effects on the transcriptomes of Vgat+Esr1+ cells (males: 56.5%, 78/138 genes; females: 73.3%, 124/169 genes; Fig. 7b). Again, because of the complex nature of this relationship between dimorphic and adolescent genes, we generated a separate manifold for each sex to assess MPOA adolescent transcriptional dynamics. Pseudotime trajectory analysis indicated that Esr1KO near-completely prevented adolescent transcriptional maturation in both sexes (Fig. 7c-h and Supplementary Fig. 12). DE analysis also identified sex-shared and -specific processes (Fig. 7i-l and Supplementary Table 10). In addition, dimorphic genes in the Esr1-DEGs of Vgat+Esr1+ cells were sufficient to predict the sex in P50 females and males with the highest accuracy (81.4 ± 0.11 %), followed by P35, and lowest accuracy in P23 (Supplementary Fig. 10d). Like our previous observation (Supplementary Fig. 6), we observed two major branches of transcriptional trajectories with branch specific gene expression profiles, primarily originating from one to two subpopulations of Vgat+Esr1+ clusters (Supplementary Fig. 12c-d and Supplementary Table 10) highlighting Esr1 as a key regulator of adolescent transcriptional dynamics in the MPOA of both sexes.

Fig. 7:

Gene regulatory networks (GRNs) are complex systems of genes, transcription factors, and other molecules that interact to control gene expression within a cell. To identify GRNs that are shared between sexes or specific to one sex, we used functional-SCENIC to examine if Esr1 cis-regulates its own differentially expressed genes (Esr1-DEGs) and/or transcription factors (Esr1-TFs) (Fig. 8a). Our analysis revealed that a significant portion of Esr1-DEGs contained Esr1 binding sites–14.9% in females (123 out of 824 genes) and 16.3% in males (98 out of 600 genes) (Fig. 8b-c)–indicating that Esr1 directly regulates many target genes. Additionally, 9 female and 6 male Esr1-DEGs encoded transcription factors that function as central regulatory hubs within the network (Fig. 8d-e), suggesting that Esr1’s influence extends indirectly through these secondary transcription factors. 54.1% of female and 28.3% of male Esr1-DEGs were cis-regulated by Esr1-TFs, and collectively, 56.1 % female and 36.7% male Esr1-DEGs were cis-regulated by Esr1, Esr1-TFs, or a combination of the two (Esr1-GRNs; Fig. 8d-e and Supplementary Fig. 13).

Fig. 8:

To validate these Esr1-GRNs in situ, we virally deleted Esr1 from MPOA neurons in female and male mice prior to the onset of adolescence and conducted HM-HCR FISH on brain sections at P50, comparing Esr1KO to Esr1-intact controls. The gene panel included genes identified via SCENIC (Details in Methods)–3 Esr1-TFs: Zmat4 (sex-shared), Pou6f2 (male-specific), and Creb3l1 (female-specific), and 6 Esr1-DEGs regulated by Esr1 or Esr1-TFs: Gria3, Pdzrn4, Ctnna2, Nefl, Tead1, and Esr1. Consistent with our scRNAseq analysis (Fig. 8b-e), sex-shared Esr1-TF, Zmat4, was reduced by Esr1KO in both sexes, while male-specific Esr1-TF, Pou6f2, was decreased only in males. Female-specific Esr1-TF, Creb3l1, was reduced in females, but also in males to a lesser degree, perhaps because of high sensitivity in HM-HCR FISH assays. The expression of sex-shared Esr1-DEGs: Gria3, Pdzrn4, and Ctnna2, were decreased in Esr1KO females and males, and female-specific Esr1-DEG, Tead1, was reduced only in females (Fig. 8f-g). These results further suggest that Esr1 orchestrates sex-shared and -specific transcriptional states via an Esr1-specific gene regulatory network. These dynamic gene regulatory networks may in turn dictate the sex-shared and -specific adolescent changes that facilitate the maturation of sexually dimorphic neuronal circuits for mating.

Discussion

Adolescence represents a critical period for the maturation of behaviors and brain functions necessary for reproduction and social interactions. Unlike homeostatic processes such as eating and breathing, many behaviors that define adulthood emerge through hormonally driven neural circuit refinements during adolescence ^1,3,4. Our study provides a detailed analysis of transcriptional dynamics in MPOA Esr1+ neurons, revealing their essential role in the adolescent maturation of mating behaviors in both sexes. Using a combination of single-cell RNA sequencing, advanced in situ hybridization techniques, and functional genomic manipulations, we show that Esr1 is a central regulator of these adolescent transcriptional changes. We show that MPOA Esr1+ neurons are transcriptionally dynamic during adolescence, progressing from an immature state in preadolescence to mature adult-like states in a hormone-dependent manner. This maturation process was arrested in mice lacking Esr1 in MPOA Vgat+ neurons during adolescence, leading to profound deficits in mating behaviors. Importantly, these findings highlight that transcriptional programs driving sexual maturation are cell-type specific and temporally constrained, dependent on precise hormonal signaling during adolescence. Adolescent transcriptional dynamics in MPOA Esr1+ neurons appear to be continuous, with trajectories shaped by circulating sex hormones. This hormone-dependent maturation was confirmed by pseudotime analysis, which revealed that hormonal deprivation arrests transcriptional progression, while hormone supplementation accelerates it. These findings align with the broader view that the adolescent brain undergoes extensive molecular and circuit-level adaptations to support the emergence of sexually dimorphic behaviors ^3,14,42.

Integration of MPOA sexual dimorphism and adolescence-related transcriptional dynamics

One of the major contributions of this study is the disentanglement of sexual dimorphism and adolescent dynamics in Esr1+ neurons. We identified subsets of genes that are both sexually dimorphic and adolescently dynamic, as well as genes that are unique to each process. This complex interplay underscores the importance of analyzing these processes separately in a sex-specific fashion to uncover their contributions to neural development. Our findings suggest that Esr1 not only regulates common gene programs shared by males and females but can also orchestrate sexually distinct transcriptional networks. The male-specific enrichment of androgen receptor-related genes and the identification of sexually dimorphic Esr1-dependent transcription factors (e.g., Pou6f2 in males, Creb3l1 in females) further highlight how Esr1 integrates hormonal signals to shape sex-specific neural circuits ^40,41.

Functional implications for neural circuitry and behavior

The adolescent transcriptional changes observed in MPOA Esr1+ neurons likely translate to functional adaptations in mating-relevant neural circuits. Previous studies have shown that MPOA neurons are critical for a wide range of sexually dimorphic behaviors, including mounting, vocalization, and parental care ^10,14,18,43. Our study provides molecular evidence linking the transcriptional maturation of MPOA Esr1+ neurons to these behaviors. This link is especially compelling considering our integration with published MERFISH datasets, which further validated the association of Vgat+Esr1+ neurons with mating behaviors in both sexes ^16,20. Dynamic changes in the circuit properties of MPOA neurons during adolescence may involve alterations in cellular excitability, synaptic connectivity, and interactions with downstream targets. Ontology analysis of Esr1-dependent genes revealed enrichment for pathways related to synaptic transmission and axonal structure, supporting the hypothesis that transcriptional maturation directly contributes to circuit remodeling ^27,44. Future studies combining transcriptional profiling and CRISPR gene editing with in vivo imaging and circuit manipulation will be essential to dissect these connections. Our study complements previous work on perinatal transcriptional dynamics in sexually dimorphic brain regions, such as the BNST ²¹. While the perinatal period is critical for establishing primary sexual differentiation, our findings emphasize adolescence as a secondary critical period where sexually dimorphic transcriptional programs are refined. Notably, both periods rely heavily on Esr1-mediated regulation, suggesting a conserved role for estrogen signaling across developmental stages.

Limitations of the study

While these results suggest that male and female adolescent transcriptional maturation has unique and overlapping components, some technical limitations of the results are worth noting. First, UMIs per cell are low due to older scRNAseq technology being used. We used HM-HCR FISH to corroborate our findings because this technique is very sensitive for detecting targeted genes. However, it is likely that the transcriptional state changes we observed are blunted due to sparse gene expression detection in 10X V2 kits compared to now-standard versions of these reagents. Additionally, while KO of Esr1 has dramatic impacts on mating behavior in both sexes, due to the massive number of genes that are regulated by Esr1, this is not surprising. Now that we have a list of targeted genes that are likely directly under the influence of Esr1, future studies are needed to resolve how individual or groups of Esr1-associated DEGs are critical for specific aspects of mating behavior.

Conclusions and future directions

These findings have implications beyond mating behaviors. The role of Esr1+ neurons in the MPOA likely extends to other hormonally sensitive behaviors and physiological processes, such as thermoregulation, metabolic control, and social bonding ^15,19,45. Additionally, understanding how disruptions in these transcriptional programs contribute to developmental disorders—such as delayed puberty, hypogonadism, or neuropsychiatric conditions with sex-biased prevalence—represents an important avenue for future research ^1,2. While our study focused on Esr1, it is likely that other transcription factors and hormone receptors synergistically contribute to adolescent transcriptional dynamics. For example, the androgen receptor may play a role in males, while additional regulatory networks may be governed by sex chromosome-linked genes ⁴¹. Advanced genomic approaches, such as perturb-seq, could help unravel the complexity of these gene regulatory networks in a cell-type-specific manner ²⁹. Although we focused on adolescent transcriptional trajectories of Vgat+Esr1+ neurons, this dataset provides a rich resource for studying transcriptional states across other neuronal, glial, and stromal MPOA cell types during adolescent brain development. Understanding steroid-induced changes in transcriptional trajectories during two major developmental periods (perinatal and adolescence periods) has implications for the healthy and maladaptive development of human brains, and provides additional insight into neurobiological mechanisms that underlie sex differences ^46–48.

Materials and Methods

Mice

Several lines of male and female mice were used for scRNAseq, HM-HCR FISH, immunohistochemistry (IHC), and behavioral experiments, with specific ages and experimental conditions described in their respective sections. All mice used in these experiments were on a C57BL/6 background. Wild-type C57BL/6J mice were originally obtained from Jackson Laboratory (JAX) and bred in-house. Transgenic Slc32a1^Flp::Esr1^lox/lox mice were generated by crossing Slc32a1^Flp mice (Jackson Laboratory stock no. 029591) with Esr1^lox/lox mice (kindly provided by Kenneth Korach; Jackson Laboratory stock no. 032173), and then breeding the resulting progeny with Slc32a1^Flp::Esr1^lox/+ mice. Similarly, Slc17a6^Flp::Esr1^lox/lox mice were generated by crossing Slc17a6^Flp mice (Jackson Laboratory stock no. 030212) with Esr1^lox/lox mice and then crossing the progeny with Slc17a6^Flp::Esr1^lox/+ mice. For both lines, only mice heterozygous for Flp and homozygous for Esr1^lox/lox were used in experiments. For scRNAseq and HCR experiments, mice were group-housed except for the 1–2 days prior to tissue isolation, during which they were singly housed. Behavioral experiment mice were singly housed for the duration of the experiments. Wild-type Swiss Webster (SW) mice, originally obtained from Taconic and bred in-house, were used as lactating foster dams to rear pups following brain surgery, as C57BL/6 dams frequently committed infanticide after such procedures. All mice were maintained with ad libitum access to food and water under a reverse 12-hour light-dark cycle. All experiments were conducted in accordance with the National Institutes of Health’s Guide for the Care and Use of Laboratory Animals and approved by the Institutional Animal Care and Use Committees (IACUC) of the University of North Carolina and University of Washington.

Stereotaxic viral injection

Juvenile mice aged P14–18 were separated from their dams, anesthetized with isoflurane, and positioned in a stereotaxic frame (Kopf Instruments). Isoflurane anesthesia was maintained at less than 0.8% throughout the procedure. Following a craniotomy above the target region, viral injections were performed using glass capillaries at a controlled rate of 60 nl/min with a Nanoject system (Drummond). To ensure accurate delivery and minimize backflow, injections were paused for 30 seconds after every 60 nl infusion, and the capillaries were left in place for 10 minutes before withdrawal. The stereotaxic coordinates used were 0.2 mm AP and 0.23 mm ML from Bregma; -4.95 mm DV from the brain surface. After recovery from anesthesia, mice were transferred to the care of lactating Swiss Webster (SW) foster dams. SW dams effectively nursed surgically treated pups, whereas C57BL/6 dams frequently committed infanticide following such procedures.

Gonadectomy

Gonadectomized mice were prepared for scRNAseq and HCR experiments. At P22–23, mice were anesthetized with isoflurane and maintained under anesthesia at less than 1.0%. After shaving the hair around the flank area for females or the scrotum for males, local anesthesia was administered, and a small incision was made through the skin and muscle above the target gonads. The gonads were identified, gently exteriorized with forceps, and removed using a heat cautery pen (Bovie Medical). The remaining tissues were returned to the abdominal cavity (females) or scrotum (males). The incision site was closed with sutures, and tissue adhesive (Vetbond) was applied to ensure secure skin closure. Brain tissues were collected at P50 ± 2, ensuring the mice were deprived of sex hormones throughout adolescence. Gonadectomized adult female mice were also prepared as stimulus animals for Resident Intruder (RI) and three-chamber sociability assays. For these experiments, 7–8-week-old C57BL/6J female mice were ovariectomized using the same procedure. After a three-week recovery period, a standard hormonal priming protocol was used to induce estrus. This involved subcutaneous injections of β-Estradiol 3-Benzoate (Sigma, 10 µg at 48 hours and 5 µg at 24 hours before the assay) and Progesterone (Sigma, 500 µg) administered six hours prior to the RI assay.

Generation of AAV-frtFlex-Cre virus

The split-Cre design used in this study was previously described ^49,50. In this system, an intron is introduced into the middle of the Cre sequence, with the second exon flanked by frt sites. The action of FLP recombinase inverts the second exon, enabling splicing to generate functional Cre. For these experiments, the construct was incorporated into an AAV vector. The frtFlex-Cre-EGFP cassette was subcloned into the pAAV-CAG-WPRE vector between AscI and XhoI restriction sites. AAV1 virus carrying frtFlex-Cre-EGFP was produced as described previously ⁵⁰. Briefly, pAAV-CAG-frtFlex-Cre-EGFP-WPRE was transiently transfected into HEK293T/17 cells (ATCC) along with the pDG1 packaging plasmid. Twenty-four hours post-transfection, cells were transferred to serum-free media, and viral capsids were harvested 48 hours later. Purification of AAV1 capsids was performed using gradient centrifugation. Viral titer was estimated at approximately 2 × 10¹² particles/mL based on densitometry analysis following gel electrophoresis, using a known standard for reference.

Validation of efficiency and specificity of AAV-frtFlex-split-Cre virus

To confirm that AAV-frtFlex-Cre does not produce functional Cre in the absence of FLP recombinase, 300 nl of AAV-frtFlex-Cre was unilaterally injected into the MPOA of Esr1^lox/lox mice (n = 3 males and 3 females; data combined across sexes) (Supplementary Fig. 1a). To validate the specificity of AAV-frtFlex-Cre, 300 nl of the virus was unilaterally injected into the MPOA of Slc32a1^Flp and Slc17a6^Flp mice. Two weeks post-transduction, brains were collected and analyzed using FISH, as described in the RNAscope section (Supplementary Fig. 1b). Six weeks after the injection, brain tissue was collected and analyzed via IHC, as detailed in the Histology section. To assess the efficiency of AAV-frtFlex-Cre, 300 nl of the virus was bilaterally injected into the MPOA of Slc32a1^Flp::Esr1^lox/lox and Slc17a6^Flp::Esr1^lox/lox mice. These mice were first used for behavioral experiments, as described in the Behavioral Experiment section, and efficiency was quantified post-behavioral testing. Quantitative results for each mouse line and sex are presented in Fig. 1 and Supplementary Fig. 1.

Single-cell preparation and cDNA library construction for scRNAseq

Male and female wild-type mice (n = 3–4 animals pooled per group) were prepared at P23 ± 1 (preadolescence, hormonally intact), P35 ± 1 (mid-adolescence, hormonally intact), and P50 ± 2 (hormonally intact early adulthood or gonadectomy groups). For gonadectomy groups, mice were gonadectomized at P22–23. Additionally, male and female Slc32a1^Flp::Esr1^lox/lox mice (n = 3–4 per group), which had received bilateral MPOA injections of 300 nl AAV-frtFlex-Cre at P14–18, were collected at P50 ± 2 (early adulthood). All mice were singly housed 1–2 days before tissue dissection.

Single-cell preparation followed previously published procedures ⁵¹. Immediately after removal from their home cages, mice were deeply anesthetized with an intraperitoneal injection of 0.2 mL sodium pentobarbital (39 mg/mL) and phenytoin sodium (5 mg/mL), followed by transcardial perfusion with ice-cold NMDG-aCSF containing transcription and translation inhibitors to minimize transcriptional events induced by anesthesia, perfusion, or brain extraction. NMDG-aCSF consisted of 96 mM NMDG, 2.5 mM KCl, 1.35 mM NaH2PO4, 30 mM NaHCO3, 20 mM HEPES, 25 mM glucose, 2 mM thiourea, 5 mM Na+ ascorbate, 3 mM Na+ pyruvate, 0.6 mM glutathione-ethyl-ester, 2 mM N-acetyl-cysteine, 0.5 mM CaCl2, and 10 mM MgSO4 (pH 7.35–7.40, 300–305 mOsm, oxygenated with 95% O2 and 5% CO2). Inhibitor cocktails included 500 nM TTX, 10 μM APV, 10 μM DNQX, 5 μM actinomycin, and 37.7 μM anisomycin. Brains were extracted, and coronal sections containing the MPOA were sliced at 300 µm using a vibratome (Leica VT1200). Slices (3–4 per animal) were recovered on ice in a chamber for 30 minutes. MPOA tissues were dissected under a light microscope (Leica MZFL) using a micro scalpel (Feather) from 3–4 animals per group (10–15 slices pooled per group). The tissue was enzymatically digested with 1 mg/mL pronase (Sigma-Aldrich) for 30 minutes at room temperature, followed by mechanical trituration using fire-polished glass capillaries (tip diameter 200–300 µm). Cell suspensions were filtered twice through 40 µm strainers to remove aggregates, and dead cells were eliminated using a dead cell removal kit (Miltenyi Biotec). After centrifugation, the supernatant was removed, and cells were resuspended (20– 30 mL). A fraction (∼5 mL) was mixed with trypan blue for viability assessment and cell concentration measurement using a hemocytometer. Samples with >80% viability were used for cDNA library preparation, and final cell concentrations were adjusted to 800–1,000 cells/µL. Typically, 30,000–80,000 cells were collected from 3–4 mice per group, exceeding the number required for cDNA library construction and reducing biological variability between subjects.

cDNA libraries were prepared using the Chromium Single Cell 3’ Reagent Kits V2 according to the manufacturer’s instructions (10x Genomics). Approximately 17,000 dissociated cells were mixed with reverse transcription mix and loaded onto a chip to recover up to 10,000 single cells. mRNAs from single cells were captured by barcoded beads in droplets using a Chromium Controller. Reverse-transcribed cDNAs were PCR amplified, fragmented, and ligated with adapters, followed by sample indexing via PCR. cDNA libraries were sequenced on an Illumina NextSeq 500 (v2.5) or HiSeq system. Sequencing reads were aligned to the mouse genome using the 10x Genomics Cell Ranger pipeline (V3) to generate cell-by-gene count matrices for downstream analysis.

Highly Multiplexed HCR FISH

HM-HCR FISH assays were conducted to validate trajectory inference from scRNAseq data (HCR1; Fig. 4 and Supplementary Fig. 8–9) and to cross-validate the Esr1-GRN analysis (HCR2; Fig. 8).

For HCR1, we analyzed groups included in the scRNAseq experiments, including preadolescence (P23 ± 1), adolescence (P35 ± 1), intact adult (P50 ± 2), and gonadectomy groups (P50 ± 2, gonadectomy performed at P22–23). Additionally, a hormone-supplemented group received daily injections of sex hormones (females: 5 µg β-Estradiol 3-Benzoate in 50 µl sesame oil; males: 200 µg Testosterone Propionate in 50 µl sesame oil) from P22 to P27, with tissue collected at P28. A control group for hormone supplementation received daily injections of sesame oil (50 µl) during the same period, and brains were harvested at P28. The timing of hormone supplementation aligned with the onset of adolescence in peripheral tissues (P28–33; Fig. 1c, f) and reports from the pituitary (P25–30; ⁵²). Adolescence onset was verified in hormone-supplemented animals by assessing male balanopreputial separation (BPS) and female vaginal opening (VO). All hormone-supplemented animals displayed BPS or VO by P28, whereas none of the controls reached adolescence by this age. All mice in the 12 groups were singly housed 1–2 days prior to tissue collection.

For HCR2, AAV-Cre-YFP or AAV-Flp-YFP was injected into the MPOA of Esr1^lox/lox mice of both sexes at P30–35, and brain tissues were harvested three weeks later. For both HCR1 and HCR2, mice (n = 2 per group) were deeply anesthetized with isoflurane, and brains were rapidly extracted and frozen on dry ice. Coronal sections (20 µm) were prepared using a cryostat (Leica) and stored at -80°C until use.

Sequential HCR FISH followed a modified protocol based on ³⁸. Tissue sections were fixed in pre-chilled 4% paraformaldehyde (PFA) for 30 minutes on ice, rinsed twice in 1x PBS at room temperature, dehydrated in ethanol (50%, 70%, 100%, 100%; 5 minutes each), and air-dried for 5 minutes. Sections were then permeabilized with protease IV (ACD, 322336) for 5 minutes at room temperature, followed by rinses in 1x PBS and 2x SSC. Probes, amplifiers, and buffers (Molecular Instruments) were prepared as per the manufacturer’s guidelines (Choi et al., 2018). Tissue was equilibrated in probe hybridization buffer for 10 minutes at 37°C in a humidified chamber, then incubated overnight with probe mixtures (final concentration 4–10 nM) at 37°C. Coverslips were removed in 100% pre-warmed wash buffer at 37°C, and tissues were sequentially washed in dilutions of wash buffer in 5x SSCT (75%, 50%, 25%; 15 minutes each), then in 5x SSCT at 37°C (15 minutes) and at room temperature (5 minutes).

For signal amplification, tissues were equilibrated in amplification buffer for 30 minutes at room temperature, then incubated overnight with snap-cooled hairpins conjugated to Alexa488, 546/594, and 647 (final concentration 60 nM). Coverslips were removed in 5x SSCT at room temperature, followed by two 30-minute washes in fresh 5x SSCT and rinses in 2x SSC. Autofluorescence was minimized with a quenching kit (Vector Laboratories, SP-8400-15). Sections were treated with quenching solution for 2 minutes, rinsed in 2x SSC, counterstained with DAPI (ACD, 320858) for 30 seconds, mounted in Vectashield (Vector Laboratories, H-1700-10), and stored at 4°C. DAPI staining was performed only in the first round, and all images were acquired within 24 hours.

Probes and amplifiers were stripped between rounds to enable multiplexing. Coverslips were floated off in 2x SSC for 30 minutes at room temperature, and tissue was incubated in DNase I (250 U/mL in 1x DNase I buffer, Roche, 04716728001) for 1.5 hours at room temperature, followed by six washes in 2x SSC (5 minutes each). Tissue was then equilibrated in pre-hybridization buffer for subsequent rounds. This process was repeated for up to four rounds, allowing detection of up to 12 different mRNAs.

Images were acquired using an Axio Imager M2 fluorescence microscope (Zeiss) equipped with Zen software. Channels for green (Alexa488), red (Alexa546/594), and far-red (Alexa647) fluorescence were captured, along with brightfield images to aid registration. DAPI signals from the first round defined nuclear regions of interest (ROIs), which were expanded to include cytoplasmic transcripts. Image registration was performed using HiPlex software (ACD), with brightfield images from subsequent rounds aligned to the first round using transformation matrices. Up to 17 images were overlaid, and overlapping regions were cropped to generate a single 12-plex image.

For quantification, HCR FISH images were analyzed in ImageJ. DAPI-defined ROIs were transferred to the 12-plex image to measure transcript numbers for each gene. For HCR1, genes enriched in Vgat+Esr1+ neurons and HA-DEGs from scRNAseq data were analyzed, with validation by SVM (Supplementary Fig. 6k, m). For HCR2, Esr1-TFs and Esr1-DEGs regulated by Esr1-TFs or Esr1 were selected from the Esr1-GRN (Fig. 8).

RNAscope

AAV-frtFlex-Cre (300 nl) was unilaterally injected into the MPOA of VgatFlp or Vglut2Flp mice (n = 2 per experiment). After two weeks of viral incubation, mice were deeply anesthetized with isoflurane, and brains were rapidly extracted and frozen on dry ice. Coronal sections (20 µm) were prepared using a cryostat (Leica) and stored at -80°C until further use. Probe hybridization and signal detection were performed according to the manufacturer’s instructions (ACDbio). In VgatFlp mice, probes targeted Slc32a1 and eGFP, while in Vglut2Flp mice, probes targeted Slc17a6 and eGFP. Tiled images of the MPOA were acquired using a Zeiss ApoTome2 fluorescence microscope with a 20x objective and Zen software (Zeiss). Image acquisition settings were consistent across all experiments. The resulting CZI files were analyzed using Fiji. The number of cells expressing Slc32a1 or Slc17a6 and eGFP, as well as double-positive cells expressing both markers, were quantified for each mouse line (Fig.1 and Supplementary Fig. 1).

Behavioral experiments

Male and female Slc32a1^Flp::Esr1^lox/lox or Slc17a6^Flp::Esr1^lox/lox mice were bilaterally injected with 300 nl of AAV-frtFlex-Cre (Cre group; virus generated in-house, detailed in the Generation of AAV-frtFlex-Cre virus section) or AAV-fDIO-eYFP (control group; UNC Vector Core) at P14–18. Starting at P25, sexual organ development was inspected daily to determine the age of first balano-preputial separation (BPS) in males or vaginal opening (VO) in females. For female subjects, vaginal smears were also collected daily following vaginal opening to monitor estrous cycles. Behavioral experiments, including the Resident Intruder (RI) assay, locomotion tests, sociability assays, and the Elevated Plus Maze (EPM), were conducted during the second half of the dark cycle under red-light illumination. Body weights were recorded daily from P30 to P54. Upon completion of the behavioral experiments, subjects were deeply anesthetized, and brain tissues were collected for histological analysis (detailed in the Histology section). Mice that failed to gain weight or received mistargeted viral injections were excluded from the analysis.

To assess the adolescent maturation of sexual behaviors, male subjects underwent RI assays from P34 to P54 every other day. In this assay, a hormonally primed adult female mouse was introduced into the male subject’s home cage for a 15-minute interaction (Cre groups: Slc32a1^Flp::Esr1^lox/lox, n = 12; Slc17a6^Flp::Esr1^lox/lox, n = 9; control groups: Slc32a1^Flp::Esr1^lox/lox, n = 13; Slc17a6^Flp::Esr1^lox/lox, n = 9). To prevent ejaculation, male subjects were separated from the female intruder as soon as thrusting began. The interactions were video-recorded, and the number of mounting and thrusting behaviors was manually counted.

Female sexual receptivity was assessed in RI assays conducted from P35 to P54 for Slc32a1^Flp::Esr1^lox/lox mice (Cre group: n = 16; control group: n = 10) and from P40 to P60 for Slc17a6^Flp::Esr1^lox/lox mice (Cre group: n = 9; control group: n = 11). Assay timing was guided by pilot experiments and conducted only when vaginal smears indicated the female was in proestrus or estrus. Female subjects were introduced into the home cage of a sexually experienced adult male mouse and allowed to interact freely for 10 minutes. To avoid ejaculation, females were separated from the male after intromission began or approximately 5 seconds after an unsuccessful mounting attempt. These interactions were recorded using an IR camera controlled by Ethovision (Noldus), and behaviors such as being mounted, intromitted, escaping, or displaying submissive postures were manually counted. Receptivity was calculated as the proportion of intromission bouts out of total mounting attempts. Stress levels were minimized and equalized across groups by selecting highly sexually experienced but minimally aggressive male residents. In rare instances of aggression, male residents were promptly removed. Unlike mice with brain-wide Esr1 deletions ²⁵, the male subjects in this study, with Esr1 knockouts limited to MPOA Vgat+ neurons during adolescence, showed no aggressive behaviors, consistent with previous findings ¹³.

Sociability tests were conducted after P45 on days when RI assays were not performed. Subjects were placed in a standard three-chamber choice arena, where one side contained a caged social stimulus and the other an object. Stimulus mice were either adult males or hormonally primed females. After a 5-minute habituation period without stimuli, a social stimulus and an object were introduced, and the subject mouse was allowed to explore for 10 minutes. The stimuli were then exchanged for a new mouse of the opposite sex and a novel object, and the subject mouse explored for an additional 10 minutes. Movements and locations were tracked using an IR camera controlled by Ethovision (Noldus). Total distance traveled during the habituation period was used to assess locomotion, while the time spent in the chamber with the social stimulus was used for statistical comparisons.

The Elevated Plus Maze (EPM) test was performed after P45 on days when neither RI assays nor sociability tests were conducted. Subjects were placed in a standard EPM arena and allowed to explore freely. After a 5-minute habituation period, movements were tracked using an IR camera controlled by Ethovision (Noldus). Time spent in open and closed arms of the maze was recorded for statistical analysis.

Histology

Histological experiments were conducted to: (1) examine viral infection in the MPOA of mice used in behavioral experiments (Slc32a1^Flp::Esr1^lox/lox or Slc17a6^Flp::Esr1^lox/lox), (2) validate the specificity of AAV-frtFlex-Cre in Esr1^lox/lox mice, and (3) confirm the specificity of reporter gene expression in Slc32a1^Flp::Esr1-Cre::RC-FLTG mice. Mice were deeply anesthetized with pentobarbital and transcardially perfused with 40 mL of 4% paraformaldehyde (PFA) in PBS. Brains were extracted and post-fixed in 4% PFA in PBS on ice for 5–7 hours before being transferred to PBS with 0.05% sodium azide (Sigma) for storage at 4°C until sectioning. Free-floating coronal brain sections (50 µm thick) were prepared using a vibratome (Leica) and stored in PBS with 0.05% sodium azide until immunohistochemistry (IHC). Sections were first washed in PBS (3 × 5 minutes) and blocked with 15% normal donkey serum (NDS) in PBST (0.3% Triton X-100) for 2 hours at room temperature. Following blocking, sections were incubated with an anti-Esr1 primary antibody (1:500, Santa Cruz, sc-542) diluted in 15% NDS in PBST for 72 hours at 4°C. After incubation, sections were washed in PBST (0.3% Triton X-100, 3 × 30 minutes) and incubated with a secondary antibody (1:500, Life Technologies donkey anti-rabbit 568 or Jackson ImmunoResearch donkey anti-rabbit 647) diluted in 15% NDS in PBST for 2 hours at room temperature. Sections were then washed in PBST (2 × 15 minutes), incubated with DAPI for 2 minutes, rinsed in PBS (2 × 15 minutes), mounted on slides, and coverslipped with mounting medium. Tiled images of the MPOA were acquired using a Zeiss ApoTome2 fluorescence microscope with a 20x objective and Zen software (Zeiss). Image acquisition settings were consistent across all experiments. The resulting CZI files were analyzed with HALO software using the ISH-IF version 1.2 module (Indica Labs). Regions of interest (ROIs) corresponding to the MPOA were manually outlined in both hemispheres across three adjacent sections from each brain. GFP- or YFP-positive cells, Esr1-positive cells, and double-positive cells were automatically detected and quantified using specific threshold settings to define phenotypes.

Data analysis for scRNAseq, MERFISH, and HM-HCR data

Analysis of scRNAseq, MERFISH, and HM-HCR FISH data utilized several R and Python packages, which were adapted and modified as needed. Previously published MERFISH data, deposited in DRYAD, was acquired from https://datadryad.org/stash/dataset/doi:10.5061/dryad.8t8s248. The data analysis workflow encompassed clustering, differential gene expression analysis, integration of differential conditions and modalities, lineage inference, gene set activity measurements, and gene regulatory network (GRN) deconstruction.

Clustering and differential gene expression analyses were performed using the Seurat V3 package ³⁰, while LISI ⁵³ was applied to evaluate integration performance across conditions. Lineage inference, representing transcriptional connectivity of neuronal clusters, was conducted with PAGA (partition-based graph abstraction; Wolf et al., 2019). Single-cell gene set activity was quantified using AUCell ²⁹. GRNs were inferred and deconstructed using the SCENIC package, which integrates GENIE3 ^29,54 and RcisTarget. These analyses depicted cis-regulatory relationships between transcription factors and differentially expressed genes. Gene co-expression networks were identified using scWGCNA ⁵⁵. To examine transcriptional state progressions by age and hormonal state, Monocle V3 ³² and MELD ³⁶ were applied. Ontology analysis was conducted using the online Enrichr platform ^56,57.

The data analysis workflow consisted of four major steps: (1) clustering and identifying cell types that exhibit transcriptional dynamics during adolescence and are relevant for reproductive behaviors; (2) trajectory analysis to quantify adolescent transcriptional progression; (3) transcriptional analysis in spatial contexts; and (4) GRN deconstruction of adolescent transcriptional dynamics. Each step is detailed in subsequent sections.

Data preprocessing and doublet removal for scRNAseq data

Preprocessing followed the procedures described in our previously published study ⁵¹. Low-expression genes, defined as those detected in fewer than three cells, and low-quality cells (total UMI < 700, total UMI > 15,000, or >20% mitochondrial gene expression) were excluded from downstream analysis. Suspected doublets were computationally removed using the DoubletDecon package (Version 1.1.5; ⁵⁸) with default settings, applying a doublet rate of 5.6 ± 1.0%.

Integrative clustering and differential gene expression analysis

To minimize the effects of batch variability and differences due to sex, age, and hormonal state while preserving the global similarity of transcriptional states within cell types, we applied the Seurat V3 integrative approach ³⁰. This method combines canonical correlation analysis (CCA; ⁵⁹) and mutual nearest neighbor analysis ⁶⁰. After preprocessing the data, which included 58,921 cells in total, gene counts were normalized using the NormalizeData function to scale counts by total UMI with a constant scale factor (10,000), followed by natural-log transformation (log1p). The FindVariableFeatures function was used to select 2,000 highly variable genes from each sample based on variance stabilizing transformation. Integration was performed pairwise using the FindIntegrationAnchors function (CCA1-40) to identify anchors and assign scores, followed by the IntegrateData function to compute an integrated gene expression matrix by iteratively constructing and subtracting transformation matrices from the original data.

The integrated expression data was used for clustering. Expression matrices were scaled, centered, and reduced using principal component analysis (PCA). A nearest-neighbor graph was constructed in PCA space (FindNeighbors; default settings), followed by Louvain clustering at a resolution of 0.8 (FindClusters). Uniform Manifold Approximation and Projection (UMAP) was then generated for visualization (RunUMAP). Initial clustering identified 32 clusters, 13 of which were neuronal, expressing canonical markers such as Thy1 or Stmn2. Cell types were determined based on marker expression (Neuron: Stmn2, Thy1; Astrocyte: Ntsr2; OPC: Gpr17; Oligodendrocyte: Mog; Microglia: C1qc; Mural cell: Tagln; Endothelial cell: Flt1; Intermediate cell: B2m; Ependymal cell: Foxj1) as reported in previous studies ^51,61,62. Neuronal and non-neuronal cell types were combined for visualization (Supplementary Fig. 2a-d).

To identify conserved markers, we used the FindConservedMarkers function to compare gene expression in each cluster against the rest of the dataset using the Wilcoxon rank-sum test. P-values were adjusted for the number of genes tested, and genes with an adjusted p-value < 0.05 were considered significantly enriched. Clustering robustness was assessed by sub-sampling 10–100% of cells and re-clustering using the same pipeline, repeated 10 times at each sampling rate.

For higher-resolution analysis, the 24,831 cells in the 13 neuronal clusters were extracted and re-clustered using the same integrative clustering approach. This resulted in 36 clusters, of which one (204 cells) was excluded due to low expression of neuronal markers. The remaining 24,627 cells in 35 neuronal clusters were analyzed using the FindConservedMarkers function to identify cluster-specific marker genes. These neuronal clusters included cells from all groups across different ages, sexes, and hormonal states.

Differentially expressed genes (DEGs) between groups were calculated using the FindMarkers function. Pairwise comparisons were performed within aggregate Vgat+ or Vglut2+ clusters for each sex (criteria: >10% expression, logFC > 0, adjusted p-value < 0.05). Hierarchical clustering of DEGs generated dendrogram trees, with a threshold (h = 3.15) used to identify DEG clusters. This analysis revealed DEG clusters with higher expression in P50 and P35 groups compared to P23 and gonadectomy (GDX) groups in both Vgat+ and Vglut2+ cells. Additional analyses identified hormonally associated DEGs (HA-DEGs) by comparing gene expression between P50 or P35 and GDX groups in individual or aggregate clusters (e.g., Vgat+Esr1+, Vgat+Esr1-Ar+, Vgat+ hormone RLow; where + indicates >50% positive cells and - indicates <10% positive cells).

To assess integration performance, Local Inverse Simpson’s Index (LISI; ⁵³) was computed in UMAP space, with LISI scores reflecting the effective number of distinct groups in the local neighborhood of each cell (Supplementary Fig. 2i). Clustering stability was further evaluated by examining lineage relationships between clusters at varying resolutions (0.2–2.0) using the clustree package ⁶³. Clustering tree analysis showed the emergence of heterogeneous Esr1 clusters at resolution 0.8 or higher (Supplementary Fig. 2k). Specific marker genes in neuronal clusters (Supplementary Fig. 3a) and Vgat+Esr1+ subclusters (Supplementary Fig. 3h) further validated clustering.

To explore transcriptional connectivity between neuronal clusters, PAGA ²⁸ was applied. In the resulting PAGA graph, nodes represented clusters, and edge thickness corresponded to connectivity strength. Globally, the PAGA graph revealed distinct separation between Vgat+ and Vglut2+ clusters. Locally, Vgat+Esr1+ clusters (Vgat2, Vgat4, and Vgat16) exhibited high connectivity and were distinct from other clusters (Fig. 2c).

Integrative analysis of MERFISH and scRNAseq data

Seurat V3 was used to jointly analyze publicly available POA MERFISH data ²⁰ and our scRNAseq data. To establish correspondence between MERFISH and scRNAseq clusters related to reproductive behaviors in adult mice, cells from behaviorally naïve subjects (lacking Fos expression data), lactating females, or individuals exhibiting aggression toward pups were excluded. Cells categorized as “Inhibitory” or “Excitatory” in the MERFISH metadata were extracted and normalized using the NormalizeData function. Highly variable genes were identified with FindVariableFeatures, and the top 60 variable genes (excluding Fos) were used for clustering, following the procedures described earlier.

Seurat clustering identified 19 neuronal clusters within the MERFISH dataset, including 10 GABAergic clusters enriched with Gad1 and 6 excitatory clusters expressing Slc17a6. Because Fos data from behaviorally naïve animals was unavailable in the MERFISH dataset, Fos-enriched clusters were defined using the following approach: For each behavioral category and sex, a threshold for Fos enrichment was set at the one-sided 95th percentile expression level. The proportion of cells above this threshold was calculated for each cell type, and Fisher’s exact test was performed to identify Fos-enriched clusters (p < 0.05 after multiple comparison correction). As Fos-enriched clusters were only observed in a subset of GABAergic clusters, subsequent analyses focused exclusively on the correspondence of GABAergic clusters between modalities.

To map MERFISH clusters onto scRNAseq data, reference (MERFISH) and query (scRNAseq, P50) datasets were aligned for each sex using the FindTransferAnchors (CCA1-30) and TransferData functions. This generated a weights matrix (W), which was used to transfer cell-type labels through the equation P = LW^T, where L is a binary classification matrix and P represents label predictions. The imputed MERFISH cluster labels on scRNAseq cells enabled the inference of scRNAseq clusters associated with social behaviors. Genes selectively enriched in scRNAseq clusters linked to sexual behaviors were identified using the FindMarkers function.

The integrative analysis of MERFISH and scRNAseq data, along with enrichment analysis of hormonally associated DEGs (HA-DEGs), consistently demonstrated that Vgat+Esr1+clusters (defined as clusters where >50% of cells express Esr1) were transcriptionally dynamic and strongly linked to sexual behaviors. Consequently, downstream analyses primarily focused on Vgat+Esr1+clusters, while differences with other Vgat+ populations (e.g., Vgat+ Esr1-Ar+ or Vgat+ hormone RLow) were highlighted when relevant.

Gene set activity analysis of scRNAseq data

To assess the aggregate expression of specific gene sets in single cells, we performed AUCell analysis (Aibar et al., 2017). AUCell calculates the Area Under the Curve (AUC) to determine whether a given gene set is enriched among the expressed genes in each cell. Using this approach, we computed gene set activity for hormonally associated DEGs (HA-DEGs; Fig. 2g), Ar-regulon in male-dimorphic genes (Fig. 5h–j), and adolescence-related versus unrelated dimorphic genes (Fig. 6g–h).

Trajectory analysis of scRNAseq and HM-HCR FISH data

Monocle 3 was used to quantify adolescent transcriptional progression in scRNAseq and HM-HCR FISH datasets ³². For scRNAseq data, Seurat clustering results were used to separately analyze Vgat+Esr1+, Vgat+Esr1-Ar+, Vgat+hormoneR^Low, and Vglut2+Esr1+ clusters in each sex. The preprocessing step involved normalizing and scaling the data using the preprocess_cds function, followed by PCA dimensional reduction based on hormonally associated DEGs (HA-DEGs; using the top 10–15 principal components). The reduce_dimension function was applied to generate UMAP embeddings, and a principal graph in UMAP space was constructed using the learn_graph function with default settings. The root node of the trajectory was assigned to cells from P23 or GDX groups, and pseudotime, representing transcriptional progression from the root state, was computed as the geodesic distance of each cell from the root node in the UMAP.

To test whether hormonally dependent transcriptional trajectories observed in Vgat+Esr1+ cells were also present in other cell types, similar Monocle trajectory analyses were conducted for Vgat+Esr1-Ar+, Vgat+hormoneR^Low, and Vglut2+Esr1+ cells. Additionally, branches of transcriptional trajectories were analyzed due to the presence of two major branches originating from one or two subpopulations of Vgat+Esr1+ cells in both sexes (Supplementary Fig. 5, 12). Branch-specific DEG analysis identified genes enriched in each branch and those shared across branches (Supplementary Fig. 5, 12).

In the HM-HCR experiments, seven genes were selected from the top 50 HA-DEGs based on fold change and selectivity (ratio of pct1 to pct2), and two genes (male) or one gene (female) were selected from the top 10 genes enriched in P23 over P50. To confirm that these selected DEGs represented pubertal and adolescent transcriptional states, decoding accuracy of experimental conditions (age and hormonal states; detailed in Decoding Conditions from scRNAseq or HM-HCR FISH) was computed and compared between HCR DEGs and control DEGs (seven from the bottom 50 DEGs and two or one from the bottom 10 genes enriched in P23 over P50; Supplementary Fig. 6k, m).

Trajectory analysis of HM-HCR data was performed similarly to scRNAseq data using Monocle 3, focusing on Vgat+Esr1+ cells. The primary difference was that log1p values (not normalized by total UMI) were used, as most detected genes were HA-DEGs, and equivalent transcript numbers were not assumed for single cells. To quantify spatial-temporal dynamics of pubertal and adolescent transcriptional progression across ages, hormonal states, and sexes, the MPOA size was normalized across groups. Pseudotime values for the medial and lateral MPOA were computed to identify spatially specific patterns. To examine sex differences in the transcriptional onset of puberty and adolescence, Vgat+Esr1+ cells from P23, P28 (non-hormonally treated controls), and P50 groups were projected into the normalized MPOA space. A transcriptional maturation index was calculated for each P28 cell based on the pseudotime of its neighborhood P23 and P50 cells.

To jointly compute transcriptional trajectories from scRNAseq and HM-HCR data, scRNAseq data was imputed using HM-HCR data via Seurat V3 (Supplementary Fig. 9). In Vgat+Esr1+ cells, Slc32a1 and one HM-HCR gene were excluded from the scRNAseq dataset. The FindTransferAnchors (CCA) and TransferData functions were used to identify anchors between the HM-HCR reference and the scRNAseq query, generating a weights matrix (W). Gene expression features were transferred as P = FW^T, where F is the gene expression matrix and P is the predicted expression matrix. This process was iterated for all HM-HCR genes to generate predicted scRNAseq data (11–12 genes). The predicted scRNAseq data was validated against the real dataset using Pearson correlation coefficient analysis. Transcriptional trajectories were then learned from the predicted data using the same Monocle pipeline as described above.

MELD analysis

To computationally validate the monocle trajectory analysis of scRNAseq data, we used the MELD Python package, which learns the transition of transcriptional states as the relative likelihood of conditions for each cell in the PHATE space (Potential of Heat-diffusion for Affinity-based Trajectory Embedding) ^36,37. PHATE preserves local and distal relationships of transcriptional states, outperforming other embedding techniques in maintaining the structure of single-cell data. The input data consisted of log-normalized expression values and clustering metadata for Vgat+Esr1+ cells, generated from Seurat V3 analysis. The analysis proceeded in two main steps. First, data was embedded into low-dimensional PHATE space. Principal components (PCs) were computed using scprep.reduce.pca (number of components = 8–15), and PHATE embeddings were generated with phate.PHATE and phate_op.fit_transform using default settings. Second, the kernel density for each condition was estimated to quantify the likelihood of cells belonging to specific states. Conditions were categorized as mature (P50, P35) or immature (P23, GDX) (Supplementary Fig. 6h-i, j, l). Kernel density estimation was performed using meld.MELD (beta = 67, knn = 7) and meld_op.fit_transform with default settings. Relative likelihoods for each condition (the probability of observing a cell in each condition) were computed using the replicate_normalize_densities function, which applies L1 normalization across samples. These MELD analyses were used to measure hormonally associated adolescent trajectories, providing computational cross-validation for the monocle trajectory analysis (Supplementary Fig. 6h-i, j, l).

DEG and trajectory analysis of scRNAseq in Esr1KO experiments

Data from Esr1KO groups were independently processed for quality control, including the removal of low- quality cells and doublets ⁵⁸. Seurat V3 was used to perform clustering and identify neuronal clusters (detailed in Integrative Clustering and Differential Gene Expression Analysis). The dataset included 6,727 cells from Esr1KOF (female) and 6,464 cells from Esr1KOM (male). Neuronal cells were re-clustered to identify Vgat+ clusters. In females, six Vgat+ clusters (eVgat1–6), seven Vglut2+ clusters (eVglut1–7), and one mixed cluster (eMix1) were identified, while in males, 12 Vgat+ clusters (eVgat1–12), six Vglut2+ clusters (eVglut1–6), and two mixed clusters (eMix1–2) were identified (Supplementary Fig. 11; Esr1KOF: 2,182 neuronal cells; Esr1KOM: 4,761 neuronal cells).

To identify Vgat+ clusters in the Esr1KO group corresponding to specific clusters in the intact (viral-free) group, a Pearson correlation coefficient analysis was performed for each sex. This correlation matrix determined which clusters in the Esr1KO dataset corresponded to Vgat+Esr1+ clusters (e.g., Vgat2, 4, 16) or Vgat+hormoneR^Low clusters (e.g., Vgat14, 17, 20) in the intact group (Supplementary Fig. 11). To address potential experimental variations introduced by viral manipulations, this correspondence was cross-validated using anchor-based integrative analysis. Using the FindTransferAnchors (CCA1-30) and TransferData functions, anchors were identified between the intact group (reference) and Esr1KO group (query), constructing a weights matrix (W). Labels (cell-type classifications) were transferred using the equation P = LW^T, where L is a binary classification matrix and P is the label prediction. Label imputation confirmed that eVgat1 and eVgat3 (female), and eVgat3 and eVgat4 (male), corresponded to the Vgat+Esr1+ population of intact groups. Similarly, eVgat5 (female) and eVgat9 (male) were closely related to the Vgat+hormoneR^Low clusters of P50 groups (Supplementary Fig. 11). Differentially expressed genes (DEGs) specific to Esr1 were computed using Seurat’s FindMarkers function by comparing the gene expression of Vgat+Esr1+ clusters between intact P50 and Esr1KO groups. For trajectory analysis, Vgat+Esr1+ cells from all groups were combined and analyzed using Monocle 3. Principal components were computed based on the same HA-DEGs as described previously, and these were used to construct UMAP embeddings. In UMAP space, a principal graph was learned, root nodes were defined, pseudotime was computed, and branch-specific analysis was performed (detailed in Trajectory Analysis of scRNAseq and HM-HCR FISH Data).

Deconstruction of gene regulatory network

To infer the cis-regulation of genes by transcription factors (TFs), we used the SCENIC computational framework ^29,51,64. SCENIC was combined with cell-type-specific gene manipulation to causally infer the gene regulatory network (GRN). First, co-expression scores between TFs and target genes were computed using GENIE3 (runGenie3) based on scRNAseq expression data from Vgat+ cells at P50. TFs were ranked by the sum of their weights for hormonally associated DEGs (HA-DEGs; logFC > 0.25). Among all 952 TFs (female) and 915 TFs (male) in the dataset, Esr1 achieved the highest sum score in both sexes. These results prompted further scRNAseq analysis of subjects with Esr1 deleted in Vgat+ cells in the MPOA (detailed in Trajectory Analysis of scRNAseq in Esr1KO Experiments). After identifying Vgat+Esr1+ clusters, genes reduced by Esr1 deletion (Esr1-DEGs) were identified using the FindMarkers function in Seurat V3. To determine which Esr1-DEGs were cis-regulated by Esr1, we applied the RcisTarget package. The RcisTarget database (mm9-tss-centered-10kb-7species.mc9nr.feather) was used to score motifs within the transcription start site (TSS) ±10 kb region and annotate associated TFs. Overrepresentation of each motif was assessed using the calcAUC function, and significantly enriched motifs were identified using the addMotifAnnotation function (normalized enrichment score threshold: NES ≥ 2). This analysis identified 123 genes in females and 98 genes in males that were significantly enriched with Esr1 motifs. Of these, 9 (female) and 6 (male) genes were TFs (Esr1-TFs). To investigate whether Esr1-TFs could cis-regulate Esr1-DEGs, we repeated RcisTarget analysis for pairs of Esr1-TFs and Esr1-DEGs. These analyses allowed us to infer genes cis-regulated by Esr1 or each Esr1-TF, collectively defining regulons for each TF. Using these regulons, we constructed the Esr1-GRN with the igraph package. In the Esr1-GRN, regulons were categorized and highlighted based on gene ontology (detailed in Ontology Analysis, Fig. 8d-e, and Supplementary Fig. 13). This approach provided a causal framework to understand the regulatory influence of Esr1 and its downstream TFs on hormonally associated transcriptional networks.

Single-cell consensus weighted gene co-expression network analysis

Single-cell consensus-weighted gene co-expression network analysis (scWGCNA) was employed to identify sex-specific co-expression networks in Vgat+Esr1+ cells ⁵⁵. The analysis aimed to detect modules (groups of coexpressed genes) highly expressed in Vgat+Esr1+ clusters. To begin, metacells were constructed using the MetacellsByGroups function. Co-expression networks were then generated with the ConstructNetwork function, using the optimal soft-power threshold identified through TestSoftPowers. Module eigengenes (MEs), which represent the principal component of gene expression for each module, were computed using the ModuleEigengenes function. Module connectivity, calculated with the ModuleConnectivity function, was used to identify hub genes—those highly connected within each module. The top 25 hub genes and their connectivity were visualized using the ModuleNetworkPlot function. This analysis revealed nine modules in females and six modules in males that were highly expressed in Vgat+Esr1+ cells (Vgat2, 4, 16). Functional enrichment for each module was assessed using Enrichr via the RunEnrichr function ^56,57. To directly compare modules between sexes, Jaccard similarity was computed across all module pairs to quantify sex-specific similarities and differences in gene co-expression networks.

Ontology analysis

Ontology analysis was performed using Enrichr to examine the enrichment of functionally related genes within Esr1-DEGs (Supplementary Fig. 13). Gene ontology (GO) terms were identified from the GO Biological Process, GO Molecular Function, and Cellular Component databases. Detailed GO terms and their associated genes are reported in Supplementary Fig. 13. To simplify interpretation, GO terms were grouped into four categories: Gene Expression, Synapse/Excitability, Morphology, and Energy. Genes associated with multiple categories were manually assigned to the most relevant category to streamline classification and analysis.

Decoding conditions from scRNAseq or HM-HCR FISH

To evaluate whether the selected features from trajectory analyses of scRNAseq and HM-HCR FISH experiments were sufficient to decode cellular conditions (e.g., age, sex, hormonal states, transcriptional maturity), we performed support vector machine (SVM) classification using the scikit-learn library. GridSearch with 10-fold cross-validation was applied to optimize the classification as described in previous studies ^65,66. SVM models were tested using both linear and radial basis function (rbf) kernels. Hyperparameters were optimized across a range of values for γ(10⁻³, 10⁻², 10⁻¹, 10¹, 10², 10³) and C (10⁻³, 10⁻², 10⁻¹, 10¹, 10², 10³). Input features for classification included differentially expressed genes (DEGs) used in constructing trajectories for scRNAseq and HM-HCR FISH. Decoding accuracies were compared against baseline accuracies calculated using randomized features (see Detailed Statistical Procedures).

Detailed statistical procedures

Statistical analyses for each experiment are described in the methods and summarized below. Behavioral data and SVM classification analyses were performed using GraphPad Prism v9.0.2. Non-parametric tests were used to compare gene expression level distributions in scRNAseq and RNAscope data, conducted in R. No formal statistical tests were used to predetermine sample sizes.

Fig. 1b: Female: unpaired t-test. t(19)=5.7941, p < 0.001., Male: unpaired t-test. t(18)=3.8128, p=0.0029.

Fig. 1c: unpaired t-test. t(24)=0.7317, p=0.4714.

Fig. 1d left: unpaired t-test. t(24)=0.1302, p=0.8975.

Fig. 1d right: unpaired t-test. Subjects, which did not become receptive by P55, were given the value 60. T(24)=4.138, p=0.0004.

Fig. 1e left: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (1.92778, 46.2667) = 21.8041, p<0.0001), main effect of group (F (1, 24) = 36.9300, p<0.0001) and interaction between group and age (F (2, 48) = 10.8824, p=0.0010). Holm-Sidak multi-comparison test was conducted. *p < 0.05, ***p < 0.001.

Fig. 1e right: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (1.88961, 45.3507) = 23.1245, p<0.0001), main effect of group (F (1, 24) = 75.2342, p<0.0001) and interaction between group and age (F (2, 48) = 10.3659, p=0.0002). Holm-Sidak multi-comparison test was conducted. *p < 0.05, **p < 0.01, ***p < 0.001.

Fig. 1f: unpaired t-test. t(23)=0.09601, p=0.9243.

Fig. 1g: unpaired t-test. t(23)=4.534, p=0.0001.

Fig. 1h left: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (3.94868, 90.8196) = 22.7095, p<0.0001), main effect of group (F (1, 23) = 66.9792, p<0.0001) and interaction between group and age (F (10, 230) = 18.7448, p<0.0001). Holm-Sidak multi- comparison test was conducted. *p < 0.05, **p < 0.01, ***p < 0.001.

Fig. 1h right: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (4.27919, 98.4214) = 13.4311, p<0.0001), main effect of group (F (1, 23) = 19.1087, p=0.0002) and interaction between group and age (F (10, 230) = 8.95218, p<0.0001). Holm-Sidak multi-comparison test was conducted. *p < 0.05, **p < 0.01.

Fig. 1j: Female: unpaired t-test. t=2.523, p=0.025. Male: unpaired t-test. t= 4.7744, p=0.00024.

Fig. 1k: First VO age: unpaired t-test. t(18)=1.285, p=0.2152.

Fig. 1l left: First estrous age: unpaired t-test. t(18)=0.5142, p=0.6133.

Fig. 1l right: First receptive age: unpaired t-test. Subjects, which did not become receptive by P55, were given the value 65. t(18)=1.664, p=0.1134.

Fig. 1m left: The number of received intromissions: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (1.933, 34.79) = 10.92, p=0.0002), no effect of group (F (1, 18) = 2.471, p=0.1334) and no interaction between group and age (F (2, 36) = 0.5656, p=0.5730).

Fig. 1m right: Receptivity: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (1.990, 35.82) = 13.80, p<0.0001), no effect of group (F (1, 18) = 2.221, p=0.1534) and no interaction between group and age (F (2, 36) = 0.1445, p=0.8660).

Fig. 1n: First BPS age: unpaired t-test. t(16)=1.000, p=0.3322.

Fig. 1o: First mount age: unpaired t-test. t(16)=0.6030, p=0.5549.

Fig. 1p left: The number of mounts: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (3.494, 55.90) = 39.72, p<0.0001), no effect of group (F (1, 16) = 0.001556, p=0.9690) and no interaction between group and age (F (10, 160) = 0.5442, p=0.8566).

Fig. 1p right: The number of thrusts: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (3.831, 61.29) = 24.70, p<0.0001), no effect of group (F (1, 16) = 0.4016, p=0.5352) and no interaction between group and age (F (10, 160) = 0.3325, p=0.9713).

Fig. 2e: Linear regression analysis. Adjusted R-squared=0.51, p=0.00025.

Fig. 2f: Linear regression analysis. Adjusted R-squared for each hormone receptor gene was reported.

Fig. 2g: Kruskal-Wallis rank sum test followed by multiple comparisons at each Vgat+ cluster. Female:

Vgat2: Kruskal-Wallis chi-squared= 235.66, df=3, p<0.0001.

Vgat4: Kruskal-Wallis chi-squared= 273.36, df=3, p<0.0001.

Vgat16: Kruskal-Wallis chi-squared= 95.768, df=3, p<0.0001.

Vgat6: Kruskal-Wallis chi-squared= 96.91, df=3, p<0.0001.

Vgat13: Kruskal-Wallis chi-squared= 8.2627, df=3, p= 0.04088.

Vgat8: Kruskal-Wallis chi-squared= 5.4462, df=3, p= 0.1419.

Vgat1: Kruskal-Wallis chi-squared= 124.49, df=3, p<0.0001.

Vgat10: Kruskal-Wallis chi-squared= 11.429, df=3, p= 0.009619.

Vgat18: Kruskal-Wallis chi-squared= 8.5876, df=3, p= 0.03531.

Vgat7: Kruskal-Wallis chi-squared= 2.5441, df=3, p= 0.4674.

Vgat3: Kruskal-Wallis chi-squared= 5.9714, df=3, p= 0.113.

Vgat9: Kruskal-Wallis chi-squared= 6.7542, df=3, p= 0.08016.

Vgat17: Kruskal-Wallis chi-squared= 7.4935, df=3, p= 0.05773.

Vgat15: Kruskal-Wallis chi-squared= 3.2132, df=3, p= 0.3599.

Vgat12: Kruskal-Wallis chi-squared= 5.316, df=3, p= 0.1501.

Vgat19: Kruskal-Wallis chi-squared= 1.7054, df=3, p= 0.6357.

Vgat5: Kruskal-Wallis chi-squared= 3.6455, df=3, p= 0.3024.

Vgat14: Kruskal-Wallis chi-squared= 1.9379, df=3, p= 0.5854.

Vgat11: Kruskal-Wallis chi-squared= 12.301, df=3, p=0.006419.

male:

Vgat20: Kruskal-Wallis chi-squared= 5.3732, df=3, p = 0.1464.

Vgat2: Kruskal-Wallis chi-squared= 309.73, df=3, p<0.0001.

Vgat4: Kruskal-Wallis chi-squared= 235.14, df=3, p<0.0001.

Vgat16: Kruskal-Wallis chi-squared= 147.85, df=3, p<0.0001.

Vgat6: Kruskal-Wallis chi-squared= 99.699, df=3, p<0.0001.

Vgat13: Kruskal-Wallis chi-squared= 69.569, df=3, p<0.0001.

Vgat8: Kruskal-Wallis chi-squared= 45.965, df=3, p<0.0001.

Vgat1: Kruskal-Wallis chi-squared= 120.17, df=3, p<0.0001.

Vgat10: Kruskal-Wallis chi-squared= 78.183, df=3, p<0.0001.

Vgat18: Kruskal-Wallis chi-squared= 54.178, df=3, p<0.0001.

Vgat7: Kruskal-Wallis chi-squared= 34.197, df=3, p<0.0001.

Vgat3: Kruskal-Wallis chi-squared= 100.71, df=3, p<0.0001.

Vgat9: Kruskal-Wallis chi-squared= 84.388, df=3, p<0.0001.

Vgat17: Kruskal-Wallis chi-squared= 24.424, df=3, p<0.0001.

Vgat15: Kruskal-Wallis chi-squared= 88.496, df=3, p<0.0001.

Vgat12: Kruskal-Wallis chi-squared= 116.7, df=3, p<0.0001.

Vgat19: Kruskal-Wallis chi-squared= 55.483, df=3, p<0.0001.

Vgat5: Kruskal-Wallis chi-squared= 163.73, df=3, p<0.0001.

Vgat14: Kruskal-Wallis chi-squared= 17.716, df=3, p=0.00050.

Vgat11: Kruskal-Wallis chi-squared= 61.638, df=3, p<0.0001.

Vgat20: Kruskal-Wallis chi-squared= 8.5476, df=3, p = 0.03595.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23, GDX) vs (P35, P50) when (P35, P50) was higher: ***p < 0.001.

Fig. 3c: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 473.01, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Fig. 3d: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared=448.36, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Fig. 3e: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1361.5, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Fig. 3f: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared=712.23, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Fig. 3h: One-way ANOVA followed by multiple comparisons. ANOVA, F (3, 396) = 13016.7483, p<0.0001. Tukey’s multi-comparison test: real Vgat+Esr1+vs shuffle Vgat+Esr1+Vgat⁺, p< 0.0001. real Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. real hormoneR^Low vs shuffle hormoneR^Low, p< 0.0001.

Fig. 3j: One-way ANOVA followed by multiple comparisons. ANOVA, F (3, 396) = 21372.43634, p<0.0001. Tukey’s multi-comparison test: real Vgat+Esr1+vs shuffle Vgat+Esr1+Vgat⁺, p< 0.0001. real Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. real hormoneR^Low vs shuffle hormoneR^Low, p< 0.0001.

Fig. 3l: Wilcoxon test. Female only genes: W=0, p<0.0001; Male only genes: W=583, p=0.23; shared genes: W=0, p<0.0001.

Fig. 3n: Wilcoxon test. Female only genes: W=1414, p=0.057; Male only genes: W=0, p<0.0001; shared genes: W=0, p<0.0001.

Fig. 5b: Linear regression analysis. Adjusted R-squared=0.51, p=0.00025.

Fig. 5c: Linear regression analysis. Adjusted R-squared for each hormone receptor gene was reported.

Fig. 5i: unpaired Wilcoxon test. From left to right, p: <0.0001, <0.0001, <0.0001, =0.49, =1, =1, <0.0001, =0.08, =1, =0.40, =0.40, =1, =1, =1, =1, =1, =1, =1, =1, =1

Fig. 5j: Linear regression analysis. (Left) Adjusted R-squared=0.12, p=0.073. (right) Adjusted R-squared=0.53, p=0.00017.

Fig. 6a:

All dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 38.68, df=3, p<0.0001.

Adolescence related dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 37.977, df=3, p<0.0001.

Adolescence unrelated dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 0.83392, df=3, p=0.8413.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Fig. 6b:

All dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 214.46, df=3, p<0.0001.

Adolescence related dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 161.79, df=3, p<0.0001.

Adolescence unrelated dimorphic genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 2.9488, df=3, p=0.3996.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Fig. 6c: One-way ANOVA followed by multiple comparisons.

Adolescence unrelated dimorphic genes: ANOVA, F (7, 792) = 34970, p<0.0001. Tukey’s multi-comparison test: real vs random shuffle, p<0.0001. real P23 vs real OVX, p=1. Real P35 vs real OVX, p=0.35. real P50 vs real OVX, p=0.44. real P35 vs real P23, p=0.49. real P50 vs real P23, p=0.59. real P50 vs real P35, p=1.

Adolescence related dimorphic genes: ANOVA, F (7, 792) = 5095, p<0.0001. Tukey’s multi-comparison test: real vs random shuffle, p<0.0001. real P23 vs real OVX, p<0.0001. Real P35 vs real OVX, p<0.0001. real P50 vs real OVX, p<0.0001. real P35 vs real P23, p<0.0001. real P50 vs real P23, p<0.0001. real P50 vs real P35, p<0.0001.

Post-hoc statistical results of Real P50 vs real others were highlighted when P50 was statistically larger than other real groups.

Fig. 6d: One-way ANOVA followed by multiple comparisons.

Adolescence unrelated dimorphic genes: ANOVA, F (7, 792) = 24080, p<0.0001. Tukey’s multi-comparison test: real vs random shuffle, p<0.0001. real P23 vs real OVX, p<0.0001. Real P35 vs real OVX, p<0.0001. real P50 vs real OVX, p<0.0001. real P35 vs real P23, p<0.0001. real P50 vs real P23, p<0.0001. real P50 vs real P35, p=0.186.

Adolescence related dimorphic genes: ANOVA, F (7, 792) = 5040, p<0.0001. Tukey’s multi-comparison test: real vs random shuffle, p<0.0001. real P23 vs real OVX, p<0.0001. Real P35 vs real OVX, p<0.0001. real P50 vs real OVX, p<0.0001. real P35 vs real P23, p<0.0001. real P50 vs real P23, p<0.0001. real P50 vs real P35, p<0.0001.

Post-hoc statistical results of Real P50 vs real others were highlighted when P50 was statistically larger than other real groups.

Fig. 6f: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 767.39, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F) and (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Fig. 6g: Kruskal-Wallis rank sum test followed by multiple comparisons.

(top) Kruskal-Wallis chi-squared= 1283.8, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P35F, P50F) vs others. ***p < 0.001. (bottom) Kruskal-Wallis chi-squared= 446.72, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P35F, P50F) vs others. ***p < 0.001.

Fig. 6h: Kruskal-Wallis rank sum test followed by multiple comparisons.

(top) Kruskal-Wallis chi-squared= 600.41, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P35M, P50M) vs others. ***p < 0.001. (bottom) Kruskal-Wallis chi-squared= 484.69, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P35M, P50M) vs others. ***p < 0.001.

Fig. 7e: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1583.4, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX, Esr1KO) vs (P35F, P50F): ***p < 0.001.

Fig. 7f: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1194.6, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX, Esr1KO) vs (P35F, P50F): ***p < 0.001.

Fig. 7g: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 2990.2, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX, Esr1KO) vs (P35M, P50M): ***p < 0.001.

Fig. 7h: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared=1341.1, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX, Esr1KO) vs (P35M, P50M): ***p < 0.001.

Fig. 7j: Wilcoxon test. Female only genes: W=0, p<0.0001; Male only genes: W=9327, p=0.17; shared genes: W=0, p<0.0001.

Fig. 7l: Wilcoxon test. Female only genes: W=2672, p=0.65; Male only genes: W=0, p<0.0001; shared genes: W=0, p<0.0001.

Fig. 8c: Female: unpaired Wilcoxon test. W=197164, p<0.0001. Male: unpaired Wilcoxon test. W=506476, p<0.0001.

Fig. 8f: unpaired Wilcoxon test. (Top to bottom) from left to right: W= 653278, 490664, 411540, 534038, 731136, 484492, 711128, 666599, 550366; P: <0.0001, <0.0001, =0.5057, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001

In the figure, p values were reported as the number of stars when values in the cre group were lower than those in control.

Fig. 8g: unpaired Wilcoxon test. (Top to bottom) from left to right: W= 154110, 96298, 154052, 94662, 25612, 116958, 133079, 139160, 83886; P: <0.0001, =0.0026, <0.0001, =0.0102, <0.0001, <0.0001, <0.0001, <0.0001, =0.73

In the figure, p values were reported as the number of stars when values in the cre group were lower than those in control.

Supplementary Fig. 1a: paired t-test. t(11)=-0.523, p=0.612.

Supplementary Fig. 1c: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (2, 48) = 11.4416728, p<0.0001), main effect of group (F (1, 24) = 16.8147242, p<0.001) and interaction between group and age (F (2, 48) = 0.5084058, p=0.60). Tukey multi-comparison test was conducted. *p < 0.05.

Supplementary Fig. 1d left: unpaired t-test. t(24)= -2.966, p=0.0067

Supplementary Fig. 1d right: unpaired t-test. t(15.466)= -3.7502, p= 0.0018

Supplementary Fig. 1e top: Body weight: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (3.062, 73.49) = 51.24, p<0.0001), no effect of group (F (1, 24) = 1.366, p=0.2539) and no interaction between group and age (F (12, 288) = 0.2786, p=0.9922). Locomotion: unpaired t-test. t(24)=1.134, p=0.2682. Time investigating male: unpaired t-test. t(24)=0.2027, p=0.8411. Time investigating female: unpaired t-test. t(24)=1.499, p=0.1469. Social Preference: unpaired t-test. t(24)=0.8857, p=0.3846. Time in open arm: unpaired t-test. t(14)=0.2707, p=0.7906.

Supplementary Fig. 1e bottom: Body weight: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (3.182, 73.19) = 202.5, p<0.0001), no effect of group (F (1, 23) = 1.788, p=0.1942) and interaction between group and age (F (12, 276) = 3.375, p=0.0001). Holm-Sidak multi-comparison test was conducted. No significance was observed at any age between groups. Locomotion: unpaired t-test. t(23)=0.06150, p=0.9515. Time investigating male: unpaired t-test. t(23)=0.5738, p=0.5717. Time investigating female: unpaired t-test. t(23)=2.048, p=0.0521. Social Preference: unpaired t-test. t(23)=1.563, p=0.1317. Time in open arm: unpaired t-test. t(14)=0.4188, p=0.6817.

Supplementary Fig. 1f top: Body weight: Two-way repeated measures ANOVA. ANOVA revealed main effect of age (F (4.942, 88.95) = 104.5, p<0.0001), no effect of group (F (1, 18) = 2.849, p=0.1087) and no interaction between group and age (F (12, 216) = 1.006, p=0.4442). Locomotion: unpaired t-test. t(18)=0.7283, p=0.4758.Time investigating male: unpaired t-test. t(18)=1.060, p=0.3030. Time investigating female: unpaired t-test. t(18)=0.5965, p=0.5583. Social Preference: unpaired t-test. t(18)=1.684, p=0.1094. Time in open arm: unpaired t-test. t(18)=0.2933, p=0.7727.

Supplementary Fig. 1f bottom: Body weight: Two-way repeated measures ANOVA followed by multiple comparisons. ANOVA revealed main effect of age (F (4.126, 66.02) = 197.4, p<0.0001), main effect of group (F (1, 16) = 8.571, p=0.0099) and interaction between group and age (F (12, 192) = 9.574, p<0.0001). Holm-Sidak multi-comparison test was conducted. *p < 0.05. Locomotion: unpaired t-test. t(16)=1.826, p=0.0865. Time investigating male: unpaired t-test. t(16)=0.9004, p=0.3813. Time investigating female: unpaired t-test. t(16)=0.1816, p=0.8582. Social Preference: unpaired t-test. t(16)=0.5261, p=0.6061. Time in open arm: unpaired t-test. t(10)=0.8410, p=0.4200.

Supplementary Fig. 1g: #Esr1 and # intromission in female: Linear regression analysis. Adjusted R-squared=0.52, p=0.00014. #Esr1 and # receptivity in female: Linear regression analysis. Adjusted R-squared=0.48, p=0.00028. #Esr1 and # mounts in male: Linear regression analysis. Adjusted R-squared=0.29, p=0.0086. #Esr1 and # thrusts in male: Linear regression analysis. Adjusted R-squared=0.20, p=0.027.

Supplementary Fig. 1h: #Esr1 and # intromission in female: Linear regression analysis. Adjusted R-squared=-0.031, p=0.47. #Esr1 and # receptivity in female: Linear regression analysis. Adjusted R-squared=-0.033, p=0.48. #Esr1 and # mounts in male: Linear regression analysis. Adjusted R-squared= -0.056, p=0.76. #Esr1 and # thrusts in male: Linear regression analysis. Adjusted R-squared=0.078, p=0.14.

Supplementary Fig. 2i: unpaired Wilcoxon test. W=396364535, p<0.0001.

Supplementary Fig. 3d left: one-way repeated measures ANOVA followed by multiple comparisons. Cluster1: ANOVA, F (135, 405) = 2.9438259, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(135)= 0.073054964, p= 0.9419. P23 vs P35, t(135)= 6.6450382, p <0.0001. P23 vs P50, t(135)= 12.347378, p <0.0001. P35 vs OVX, t(135)= 4.0952302, p< 0.0001. P50 vs OVX, t(135)= 13.269141, p <0.0001. P35 vs P50, t(135)= 11.17 6780, p <0.0001.

Cluster2: ANOVA, F (282, 846) = 4.6794342, P<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(282)= 22.482676, p <0.0001. P23 vs P35, t(282)= 29.289408, p <0.0001. P23 vs P50, t(282)= 31.215903, p <0.0001. P35 vs OVX, t(282)= 5.2700226, p< 0.0001. P50 vs OVX, t(282)= 12.850838, p <0.0001. P35 vs P50, t(282)= 9.6808179, p <0.0001.

Cluster3: ANOVA, F (354, 1062) = 3.9834674, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(354)= 51.772642, p <0.0001. P23 vs P35, t(354)= 20.469338, p <0.0001. P23 vs P50, t(354)= 40.905873, p <0.0001. P35 vs OVX, t(354)= 28.730343, p< 0.0001. P50 vs OVX, t(354)= 11.831025, p <0.0001. P35 vs P50, t(354)= 20.702854, p <0.0001.

Cluster4: ANOVA, F (356, 1068) = 4.6642539, P<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(356)= 2.6050178, p <0.0096. P23 vs P35, t(356)= 23.110498, p <0.0001. P23 vs P50, t(356)= 6.1529826, p <0.0001. P35 vs OVX, t(356)= 33.395442, p< 0.0001. P50 vs OVX, t(356)= 11.259059, p <0.0001. P35 vs P50, t(356)= 30.387666, p <0.0001.

Cluster5: ANOVA, F (29, 87) = 6.4594382, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(29)= 3.9693469, p = 0.0017. P23 vs P35, t(29)= 1.4257659, p = 0.3021. P23 vs P50, t(29)= 0.040574193, p = 0.9679. P35 vs OVX, t(29)= 6.7801762, p< 0.0001. P50 vs OVX, t(29)= 7.3425963, p <0.0001. P35 vs P50, t(29)= 3.2891624, p = 0.0079.

Cluster6: ANOVA, F (339, 1017) = 3.7357024, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(339)= 27.925103, p < 0.0001. P23 vs P35, t(339)= 1.1969343, p = 0.2322. P23 vs P50, t(339)= 39.356609, p < 0.0001. P35 vs OVX, t(339)= 26.966204, p< 0.0001. P50 vs OVX, t(339)= 24.477635, p <0.0001. P35 vs P50, t(339)= 41.561976, p < 0.0001.

Supplementary Fig. 3d right: one-way repeated measures ANOVA followed by multiple comparisons. Cluster1: ANOVA, F (770, 2310) = 4.8206939, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(770)= 12.113758, p < 0.0001. P23 vs P35, t(770)= 9.6124052, p < 0.0001. P23 vs P50, t(770)= 26.901991, p < 0.0001. P35 vs ORX, t(770)= 28.711090, p< 0.0001. P50 vs ORX, t(770)= 54.978801, p <0.0001. P35 vs P50, t(770)= 31.862981, p < 0.0001.

Cluster2: ANOVA, F (1089, 3267) = 5.6620007, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(1089)= 77.978449, p < 0.0001. P23 vs P35, t(1089)= 48.804692, p < 0.0001. P23 vs P50, t(1089)= 55.618808, p < 0.0001. P35 vs ORX, t(1089)= 36.331287, p< 0.0001. P50 vs ORX, t(1089)= 19.946091, p <0.0001. P35 vs P50, t(1089)= 10.943465, p < 0.0001.Cluster3: ANOVA, F (189, 567) = 3.4899665, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(189)= 5.6579721, p < 0.0001. P23 vs P35, t(189)= 11.699347, p < 0.0001. P23 vs P50, t(189)= 8.0593044, p < 0.0001. P35 vs ORX, t(189)= 14.879007, p< 0.0001. P50 vs ORX, t(189)= 12.488563, p <0.0001. P35 vs P50, t(189)= 14.346162, p < 0.0001.

Cluster4: ANOVA, F (1030, 3090) = 6.0976581, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(1030)= 30.413779, p < 0.0001. P23 vs P35, t(1030)= 37.654898, p < 0.0001. P23 vs P50, t(1030)= 38.951160, p < 0.0001. P35 vs ORX, t(1030)= 11.231385, p< 0.0001. P50 vs ORX, t(1030)= 13.001623, p <0.0001. P35 vs P50, t(1030)= 5.5005486, p < 0.0001.

Cluster5: ANOVA, F (118, 354) = 3.3165807, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(118)= 6.3281474, p < 0.0001. P23 vs P35, t(118)= 7.1294177, p < 0.0001. P23 vs P50, t(118)= 6.5731275, p < 0.0001. P35 vs ORX, t(118)= 2.0573177, p= 0.0419. P50 vs ORX, t(118)= 11.363377, p <0.0001. P35 vs P50, t(118)= 14.575549, p < 0.0001.

Supplementary Fig. 3e left: one-way repeated measures ANOVA followed by multiple comparisons. Cluster1: ANOVA, F (106, 318) = 1.784, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(106)= 17.25, p<0.0001. P23 vs P35, t(106)= 5.938, p<0.0001. P23 vs P50, t(106)= 3.823, p =0.0004. P35 vs OVX, t(106)= 17.54, p< 0.0001. P50 vs OVX, t(106)= 9.249, p <0.0001. P35 vs P50, t(106)= 0.9242, p =0.3575.

Cluster2: ANOVA, F (87, 261) = 6.329, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(87)= 10.19, p<0.0001. P23 vs P35, t(87)= 8.772, p<0.0001. P23 vs P50, t(87)= 10.93, p<0.0001. P35 vs OVX, t(87)= 0.6321, p=0.5291. P50 vs OVX, t(87)= 7.370, p <0.0001. P35 vs P50, t(87)= 6.112, p<0.0001.

Cluster3: ANOVA, F (110, 330) = 2.172, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(110)= 4.682, p<0.0001. P23 vs P35, t(110)= 16.20, p<0.0001. P23 vs P50, t(110)= 3.059, p=0.0056. P35 vs OVX, t(110)= 11.90, p <0.0001. P50 vs OVX, t(110)= 1.708, p =0.0904. P35 vs P50, t(110)= 15.13, p<0.0001.

Cluster4: ANOVA, F (11, 33) = 6.390, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(11)= 3.374, p=0.0185. P23 vs P35, t(11)= 4.954, p=0.0026. P23 vs P50, t(11)= 3.743, p=0.0129. P35 vs OVX, t(11)= 4.328, p =0.0060. P50 vs OVX, t(11)= 1.324, p =0.2122. P35 vs P50, t(11)= 2.866, p=0.0305.

Cluster5: ANOVA, F (201, 603) = 4.559, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(201)= 29.20, p<0.0001. P23 vs P35, t(201)= 14.46, p<0.0001. P23 vs P50, t(201)= 32.37, p<0.0001. P35 vs OVX, t(201)= 26.25, p <0.0001. P50 vs OVX, t(201)= 2.931, p =0.0038. P35 vs P50, t(201)= 17.71, p<0.0001.

Supplementary Fig. 3e right: one-way repeated measures ANOVA followed by multiple comparisons. Cluster1: F (706, 2118) = 5.635, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(706)= 52.31, p < 0.0001. P23 vs P35, t(706)= 28.25, p < 0.0001. P23 vs P50, t(706)= 33.39, p < 0.0001. P35 vs ORX, t(706)= 39.25, p< 0.0001. P50 vs ORX, t(706)= 19.27, p <0.0001. P35 vs P50, t(706)= 10.93, p < 0.0001.

Cluster2: ANOVA, F (520, 1560) = 5.850, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(520)= 13.78, p < 0.0001. P23 vs P35, t(520)= 21.64, p < 0.0001. P23 vs P50, t(520)= 22.69, p < 0.0001. P35 vs ORX, t(520)= 9.028, p< 0.0001. P50 vs ORX, t(520)= 10.51, p <0.0001. P35 vs P50, t(520)= 5.259, p < 0.0001.

Cluster3: ANOVA, F (176, 528) = 3.522, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(176)= 15.68, p < 0.0001. P23 vs P35, t(176)= 4.169, p < 0.0001. P23 vs P50, t(176)= 8.638, p < 0.0001. P35 vs ORX, t(176)= 17.51, p< 0.0001. P50 vs ORX, t(176)= 26.15, p <0.0001. P35 vs P50, t(176)= 13.70, p < 0.0001.

Cluster4: ANOVA, F (164, 492) = 5.971, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(164)= 29.20, p < 0.0001. P23 vs P35, t(164)= 33.47, p < 0.0001. P23 vs P50, t(164)= 31.91, p < 0.0001. P35 vs ORX, t(164)= 5.824, p< 0.0001. P50 vs ORX, t(164)= 9.851, p <0.0001. P35 vs P50, t(164)= 4.902, p < 0.0001.

Cluster5: ANOVAF (55, 165) = 3.797, p<0.0001. Holm-Sidak multi-comparison test: P23 vs ORX, t(55)= 6.306, p < 0.0001. P23 vs P35, t(55)= 6.165, p < 0.0001. P23 vs P50, t(55)= 1.518, p = 0.1348. P35 vs ORX, t(55)= 3.378, p= 0.0027. P50 vs ORX, t(55)= 11.23, p <0.0001. P35 vs P50, t(55)= 7.715, p < 0.0001.

Supplementary Fig. 3f: Linear regression analysis. Adjusted R-squared=0.43, p=0.013.

Supplementary Fig. 3g: Linear regression analysis. Adjusted R-squared for each hormone receptor gene was reported.

Supplementary Fig. 4a, b: Fisher’s exact test was conducted to compare each cluster with all clusters. P-values were Bonferroni corrected. **p < 0.01, ***p < 0.001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5c left branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 142.61, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5c left branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 133.17, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F).

Supplementary Fig. 5c left branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 128.99, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5c right branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 27.447, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F).

Supplementary Fig. 5c right branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 510.85, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5c right branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 117.43, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5d left branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 278.57, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5d left branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 192.7, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M).

Supplementary Fig. 5d left branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 137.74, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5d right branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 72.355, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M).

Supplementary Fig. 5d right branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 646.97, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5d right branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared=139.53, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5e branch1: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 435.41, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5e branch2: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 271.89, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 5f branch1: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 565.82, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5f branch2: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 540.99, df=3, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 5g branch1: unpaired t-test. T(152.39)= 216.35, p<0.0001.

Supplementary Fig. 5g branch2: unpaired t-test. T(147.5)= 153.52, p<0.0001.

Supplementary Fig. 5h branch1: unpaired t-test. T(130.85)= 244.88, p<0.0001.

Supplementary Fig. 5h branch2: unpaired t-test. T(124.65)= 204.9, p<0.0001.

Supplementary Fig. 5i branch1 AUCell score: unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5i branch2 AUCell score: unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5i branch1 SVM: unpaired t-test. T(134.62)= 142.72, p<0.0001.

Supplementary Fig. 5i branch2 SVM: unpaired t-test. T(176.75)= 66.219, p<0.0001.

Supplementary Fig. 5j branch1 AUCell score: unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5j branch2 AUCell score: unpaired Wilcoxon test. P=0.0050.

Supplementary Fig. 5j branch1 SVM: unpaired t-test. T(164.27)= 83.34, p<0.0001.

Supplementary Fig. 5j branch2 SVM: unpaired t-test. T(190.61)= 71.154, p<0.0001.

Supplementary Fig. 6d: One-way ANOVA followed by multiple comparisons. ANOVA, F (3, 396) = 15075, p<0.0001. Tukey’s multi-comparison test: real male vs shuffle male, p< 0.0001. real male vs real female, p< 0.0001. real male vs shuffle female, p< 0.0001. shuffle male vs real female, p< 0.0001. shuffle male vs shuffle female, p= 0.6720. real female vs shuffle female, p< 0.0001.

Supplementary Fig. 6f: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 830.08, df=7, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23, GDX) vs (P35, P50): ***p < 0.001.

Supplementary Fig. 6j: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 570.12, df=3, p<0.0001.

Supplementary Fig. 6k: unpaired t-test. T(198)= 61.044, p<0.0001.

Supplementary Fig. 6l: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 643.82, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 6m: unpaired t-test. T(198)= 99.879, p<0.0001.

Supplementary Fig. 6o: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 228.06, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 6q: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 310.29, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 7d: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 146.73, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig. 7e: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 515.86, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig. 9a: Kruskal-Wallis rank sum test followed by multiple comparisons.

Pdzrn4: Kruskal-Wallis chi-squared=1998, df=5, p<0.001. Creb3l1: Kruskal-Wallis chi-squared=507.5, df=5, p<0.001. Sccpdh: Kruskal-Wallis chi-squared=2517.7, df=5, p<0.001. Esr1: Kruskal-Wallis chi-squared=166.37, df=5, p<0.001. Socs2: Kruskal-Wallis chi-squared=1609.2, df=5, p<0.001. Pgr: Kruskal-Wallis chi-squared=1928.7, df=5, p<0.001. Tmem35a: Kruskal-Wallis chi-squared=1450.7, df=5, p<0.001. Nts: Kruskal-Wallis chi-squared=631.36, df=5, p<0.001. Slc32a1: Kruskal-Wallis chi-squared=133.17, df=5, p<0.001. Alcam: Kruskal-Wallis chi-squared=657.22, df=5, p<0.001. Ar: Kruskal-Wallis chi-squared=546.52, df=5, p<0.001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23F, P28C, OVX) vs (P28E, P35F, P50F): *p < 0.05, **p < 0.01, ***p < 0.001.

Supplementary Fig. 9b: Kruskal-Wallis rank sum test followed by multiple comparisons.

Sytl4: Kruskal-Wallis chi-squared=1532.5, df=5, p<0.001. Acvr1c: Kruskal-Wallis chi-squared=936.37, df=5, p<0.001. Greb1: Kruskal-Wallis chi-squared=512.3, df=5, p<0.001. Nrip: Kruskal-Wallis chi-squared=759.41, df=5, p<0.001. Lamp5: Kruskal-Wallis chi-squared=1052.2, df=5, p<0.001. Apoc3: Kruskal-Wallis chi-squared=938.64, df=5, p<0.001. Pgr: Kruskal-Wallis chi-squared=606.25, df=5, p<0.001. Esr1: Kruskal-Wallis chi-squared=87.281, df=5, p<0.001. Slc32a1: Kruskal-Wallis chi-squared=192.37, df=5, p<0.001. Npy2r: Kruskal-Wallis chi-squared=528.2, df=5, p<0.001. Pgr151: Kruskal-Wallis chi-squared=1920.3, df=5, p<0.001. Ar: Kruskal-Wallis chi-squared=235.83, df=5, p<0.001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (P23M, P28C, ORX) vs (P28T, P35M, P50M): *p < 0.05, **p < 0.01, ***p < 0.001.

Supplementary Fig. 9d: One-way ANOVA followed by multiple comparisons. ANOVA, F (3, 396) = 43666.7034, p<0.0001. Tukey’s multi-comparison test: real Vgat+Esr1+vs shuffle Vgat+Esr1+Vgat⁺, p< 0.0001. real Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. real hormoneR^Low vs shuffle hormoneR^Low, p< 0.0001.

Supplementary Fig. 9g: One-way ANOVA followed by multiple comparisons. ANOVA, F (3, 396) = 36690.76438, p<0.0001. Tukey’s multi-comparison test: real Vgat+Esr1+vs shuffle Vgat+Esr1+Vgat⁺, p< 0.0001. real Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR^Low, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR^Low, p< 0.0001. real hormoneR^Low vs shuffle hormoneR^Low, p< 0.0001.

Supplementary Fig. 9n: unpaired t-test. T(198)=91.01, p<0.0001.

Supplementary Fig. 9q: unpaired t-test. T(198)=157.2, p<0.0001.

Supplementary Fig. 9r: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 2815.7, df=11, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparisons between (P23F, P28FC, OVX) vs (P28FE, P35F, P50F) in M or L and for comparisons between M and L within each group: **p < 0.01, ***p < 0.001.

Supplementary Fig. 9s: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 2847.3, df=1, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparisons between (P23M, P28MC, ORX) vs (P28MT, P35M, P50M) in M or L and for comparisons between M and L within each group: **p < 0.01, ***p < 0.001.

Supplementary Fig. 9u: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 73.89, df=3, p<0.0001.

p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between groups. ***p < 0.001.

Supplementary Figure 10d: One-way ANOVA followed by multiple comparisons. ANOVA, F (5, 594) = 5120.211650, p<0.0001. Tukey’s multi-comparison test: real P50 vs shuffle P50, p< 0.0001. real P50 vs real P35, p< 0.0001. real P50 vs shuffle P35, p< 0.0001. real P50 vs real P23, p< 0.0001. real P50 vs shuffle P23, p< 0.0001. shuffle P50 vs real P35, p< 0.0001. shuffle P50 vs shuffle P35, p< 0.0001. shuffle P50 vs real P23, p< 0.0001. shuffle P50 vs shuffle P23, p= 0.9968. real P35 vs shuffle P35, p< 0.0001. real P35 vs real P23, p= 0.0004. real P35 vs shuffle P23, p< 0.0001. real P23 vs shuffle P35, p< 0.0001. shuffle P35 vs shuffle P23, p< 0.0001. real P23 vs shuffle P23, p< 0.0001.

Supplementary Fig.12c left branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 330.26, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F): ***p < 0.001, **p < 0.01.

Supplementary Fig.12c left branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 26.497, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F).

Supplementary Fig.12c left branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 333.74, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig.12c right branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 84.947, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F).

Supplementary Fig.12c right branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 319.6, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig.12c right branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 380.94, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF, P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig.12d left branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1136.2, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig.12d left branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 109.38, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M).

Supplementary Fig.12d left branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 613.41, df=4, p<0.0001. p values for multi-

comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig.12d right branch1 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 153.86, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M).

Supplementary Fig.12d right branch2 genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 447.68, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig.12d right branch common genes: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 596.01, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig.12e branch1: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1025.1, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF,P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig.12e branch2: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1174.6, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOF,P23F, OVX) vs (P35F, P50F): ***p < 0.001.

Supplementary Fig.12f branch1: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1014.8, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Supplementary Fig.12f branch2: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared= 1113.1, df=4, p<0.0001. p values for multi-comparison tests (Bonferroni corrected) were reported for comparison between (Esr1KOM, P23M, ORX) vs (P35M, P50M): ***p < 0.001.

Resource and Reagent Table

Data and materials availability

The NCBI Gene Expression Omnibus accession number for the scRNAseq data reported in this paper for the reviewers is GEO: GSE172177 All the codes used to analyze scRNAseq, and HM-HCR FISH are available at a Github repository affiliated with Stuber Laboratory group and this manuscript title (http://www.github.com/stuberlab/).

Acknowledgements

We thank M. Borrego for critical comments on the manuscript. We thank C. Trapnell for helpful discussions. We thank Y. Tao and Z. Hu (University of North Carolina) for assistance with initial 10x Genomics library generation.

Additional information

Funding

National Institute of Health grant DA038168 (GDS)

National Institute of Health grant DA032750 (GDS)

National Institute of Health grant DA054317 (GDS)

The Foundation of Hope (GDS, JM)

Brain and Behavior Research Foundation (KH and MAR)

National Institute of Health grant DK121883 (MAR)

National Institute of Health grant P30DA048736

Author contributions

Conceptualization: GDS, KH, YH, KI, BB, JM, DR

Methodology: KH, YH, KI, BB, MLB, RDP, LSZ

Investigation: KH, YH, KI, BB, YL, MAR, JYC, ORA, RVM, NLJ, RCS, JCS, JRS,

Data analysis: KH, YH, KI, BB, GDS, JCS, JRS

Manuscript writing: GDS, KH, YH, BB, KI

Manuscript editing: all authors

Supervision: GDS

Additional files

Supplementary Table 1. Data table showing quality of sequenced libraries and corresponding scRNAseq data from samples used in this study.

Supplementary Table 2. Data table showing the marker genes of neuronal clusters, related to Fig. 2.

Supplementary Table 3. Data table showing the DEGs in the pseudobulk populations between groups, related to Supplementary Fig. 3.

Supplementary Table 4. Data table showing the HA-DEGs in each Vgat+ cluster, related to Supplementary Fig. 3.

Supplementary Table 5. Data table showing the HA-DEGs in the Vgat+Esr1+ population, related to Supplementary Fig. 3.

Supplementary Table 6. Data table showing co-expression scores between transcription factors and adolescent gene programs or sexually dimorphic programs related to Fig. 2 and 5.

Supplementary Table 7. Data table showing the marker genes in the clusters associated with mating, related to Fig. 2 and Supplementary Fig. 4.

Supplementary Table 8. Data table showing the genes enriched in later trajectories or branch-DEGs in the Vgat+Esr1+ population, related to Fig. 3 and Supplementary Fig. 5.

Supplementary Table 9. Data table showing the Esr1-DEGs in the Vgat+Esr1+ or Vgat+/HormoneR^low populations, related to Fig. 7.

Supplementary Table 10. Data table showing the genes enriched in later trajectories or branch-DEGs in the Vgat+Esr1+ population, related to Fig. 7 and Supplementary Fig. 12.

Supplementary Table 11. Data table showing the sexually dimorphic genes in the Vgat+Esr1+ population, related to Fig. 7.

Supplementary Table 12. In situ hybridization probe list.