Esr1-dependent signaling and transcriptional maturation in the medial preoptic area of the hypothalamus shape the development of mating behavior during adolescence

Previous studies have shown that knocking out Esr1, either globally or in specific cell types, reduces mating behavior (Wu and Tollkuhn, 2017; Ogawa et al., 1997). 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 days (P) 14–18. This viral strategy allowed for Esr1 deletion in MPOA GABAergic or glutamatergic neurons, depending on the transgenic mouse line, throughout adolescence (Figure 1a–b and i–j and Figure 1—figure supplement 1a–b). The onset of genital changes, the maturation of the HPG axis (as determined by the first day of estrus in female mice Cheong et al., 2015 and first day of initiation of sexual behaviors in male mice Karigo and Deutsch, 2022), and locomotor behaviors were not altered in MPOA^Vgat-Esr1KO and MPOA^{Vglut2-Esr1KO} mice (Figure 1c–d, f, k–l and n and Figure 1—figure supplement 1e–f). The establishment of mating behavior was abolished in MPOA^Vgat-Esr1KO mice in both sexes (Figure 1d–e and g–h and Figure 1—figure supplement 1a–b and g), but unaffected in MPOA^{Vglut2-Esr1KO} mice (Figure 1l–m and o–p and Figure 1—figure supplement 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 (Wu and Tollkuhn, 2017; Sano et al., 2016; Sano et al., 2013).

Figure 1 with 1 supplement see all

Download asset Open asset

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 eight experimental groups (Figure 2a, Figure 2—figure supplement 1a–e, and Supplementary file 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: 5395; median genes per cell: 2592), resulting in 20 GABAergic clusters (Vgat +clusters, 13,334 cells; 54.4% of neurons), 12 glutamatergic clusters (Vglut2 +clusters, 7658 cells; 31.1% of neurons), and 3 clusters with mixed glutamatergic and GABAergic markers (3635 cells; 14.8% of neurons; Figure 2b, Figure 2—figure supplements 1f-l and 2a,j, and Supplementary file 2), largely consistent with previously published data (Moffitt et al., 2018; Hrvatin et al., 2020). To further quantify transcriptional associations between MPOA neuronal cell types, we performed partition-based graph abstraction (PAGA) analysis (Wolf et al., 2019a) to identify cell cluster similarities through a graphical representation of cell differentiation (Figure 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 (Figure 2c and Figure 2—figure supplements 1m and 2h–i), suggesting that Vgat +Esr1+clusters form a distinct transcriptional unit.

Figure 2 with 3 supplements see all

Download asset Open asset

Irrespective of age, sex, and hormonal states, we found that each identified cell type was represented across all 8 experimental groups (Figure 2—figure supplements 1d,h-i and 2a), 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 (Figure 2d-f, Figure 2—figure supplement 2a-g and Supplementary files 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 (Aibar et al., 2017). Regression analysis revealed a strong enrichment of HA-DEGs in Vgat +Esr1+neurons (Figure 2e–f). Consistent with this result, co-expression analysis demonstrated that out of more than 1000 transcription factors (TFs), Esr1 was one of the most co-expressed TFs in both sexes (Figure 2f and Supplementary file 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; Figure 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 (Figure 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 (Moffitt et al., 2018) which identified social behavior-activated POA cells via mating, parenting, and fighting behavior assays using Fos (Figure 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 (Figure 2—figure supplement 3a and b). After integrating our scRNAseq clusters with the MERFISH clusters (Details in Materials and methods. Seurat V3 Integrative Analysis: reference-based integration; label transfers Stuart et al., 2019), we were able to identify mating, parenting, and fighting-relevant neuronal clusters within our dataset (Figure 2h, Figure 2—figure supplement 3c–d, and Supplementary file 7). We found that Esr1 was the most enriched and selective gene in the mating-relevant scRNAseq Vgat +clusters (Figure 2h–i and Supplementary file 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 (Karigo et al., 2021). 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 (Figure 2d and Figure 2—figure supplement 2b–c), they do 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 (Heumos et al., 2023). To resolve these dynamic changes, we applied pseudotime analysis (Figure 3 and Figure 3—figure supplements 1–2), 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 (Cao et al., 2019; Qiu et al., 2017).

Figure 3 with 3 supplements see all

Download asset Open asset

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 (Pape et al., 2024). In both sexes, Vgat +Esr1+trajectories revealed that transcriptional states at P35 and P50 branch and diverge from preadolescence (Figure 3a-f, Figure 3—figure supplement 1a-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 (Figure 3g,i, Figure 3—figure supplement 2a-c). Pairwise differentially expressed gene (DEG) analysis consistently showed a larger number of DEGs between P35 and P23 in Vgat +Esr1+ (male: 146 genes; female: 162 genes) than in Vgat +hormone R^Low (male: 26 genes; female: 1 gene). Furthermore, all Vgat +Esr1+clusters were found to co-express the androgen receptor (AR) gene (Ar) to a high degree (68.3–88.3%, Figure 3—figure supplement 2a). However, Vgat +clusters expressing Ar but not Esr1 (Vgat +Esr1- Ar+) showed only a moderate transcriptional progression in males and minimal changes in females (Figure 3—figure supplement 2b–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 (Figure 3—figure supplement 3), 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, Branch 1 was dominated by Vgat +cluster 4, while Branch 2 was occupied by clusters 2 and 16. In males, Branch 1 was composed of Vgat +clusters 4 and 16, while Branch 2 was primarily formed by Vgat +cluster 2 (Figure 3—figure supplement 1g–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 (Figure 3—figure supplement 1c–d and Supplementary file 8). When we combined the trajectories for both sexes into a joint manifold, the model failed to capture hormone-dependent dynamics (Figure 3—figure supplement 2e–g), emphasizing that separate models are necessary for accurate analysis of sex-specific and sex-shared transcriptional progressions (Figure 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 (Weinreb et al., 2020). 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%; Figure 3h and 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; Burkhardt et al., 2021; Moon et al., 2019). 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 (Figure 3—figure supplement 2h–m). These results, visualized in Potential of Heat-diffusion for Affinity-based Trajectory Embedding (PHATE) space (Burkhardt et al., 2021), further highlight that while MPOA neuronal cell types are terminally diversified by preadolescence (P23; Figure 2—figure supplements 1h and 2e), Vgat +Esr1+neurons continue to undergo hormonally dependent transcriptional refinement from adolescence into adulthood.

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; von Buchholtz et al., 2021; Choi et al., 2018). This technique allowed us to detect transcripts of ~12 genes across 41,549 MPOA cells at single-molecule resolution (Figure 4a–c and Figure 4—figure supplements 1–2). The rationale for selecting this gene panel is related to Figure 3—figure supplement 2k and m with details provided in the Materials and 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 Figure 4a) to investigate whether early sex steroid supplementation accelerated adolescent transcriptional trajectories.

Figure 4 with 2 supplements see all

Download asset Open asset

To distinguish Vgat +Esr1+from Vgat +hormone R^Low cells, we measured Slc32a1, Esr1, Ar, and eight to nine DEGs during adolescence (Figure 4d–e and Figure 4—figure supplements 1–2. Details in Materials and 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; Figure 4f–g and Figure 4—figure supplement 2a–c, f and i–q). Interestingly, the transcriptional trajectories of female and male adolescents differed spatially, but both were regulated by sex hormones (Figure 4—figure supplement 2r–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%; Figure 4—figure supplement 2d and g). Classification accuracy further decreased when individual genes were iteratively removed (Figure 4—figure supplement 2e and 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 (Knoedler et al., 2022; Xu et al., 2012). 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 (Figure 5—figure supplement 1a and Figure 7—figure supplement 1). 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 (Figure 5a–e, Figure 5—figure supplement 1a, and Supplementary file 6). Notably, only subsets of sexually dimorphic genes were also classified as adolescent genes (Figures 5f and 6a,b) and some of these adolescent genes were shared across sexes (Figure 5g).

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 6

Download asset Open asset

Previous research has shown that Esr1 plays a role in regulating sexual dimorphisms within sexually dimorphic brain nuclei (Gegenhuber et al., 2022). However, 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 (Aibar et al., 2017). 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 (Figure 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; Figure 5i). Consistent with this observation, regression analysis showed that Esr1 expression, rather than Ar, predicted the expression of male dimorphic AR-regulon genes (Figure 5j).

As previously noted, a subset of sexually dimorphic genes overlaps with adolescent-related genes (Figure 5f). The number of sexually dimorphic genes increases during adolescence but reverts following gonadectomy (Figure 6a–b). Among 45 female dimorphic genes, 12 showed increased expression during adolescence (adolescence-related dimorphic genes), while 7 were not associated with adolescence (Figure 6a, right). Of the 110 dimorphic genes in males, 52 were adolescence-related dimorphic genes and 10 were not (Figure 6b, right). The expression levels of adolescence-related dimorphic genes were the highest at P50 and were again reduced by gonadectomy in both sexes (Figure 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 (Figure 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 (Figure 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 (Figures 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.

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 (Figure 7a and Figure 7—figure supplement 1. Details in Materials and 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 (Figure 7b, Figure 5—figure supplement 1c, and Supplementary file 9). Previously identified HA-DEGs (Figure 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; Figure 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 (Figure 7c–h and Figure 7—figure supplement 2). DE analysis also identified sex-shared and -specific processes (Figure 7i–l and Supplementary file 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 (Figure 5—figure supplement 1d). Like our previous observation (Figure 3—figure supplement 2), 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 (Figure 7—figure supplement 2c–d and Supplementary file 10) highlighting Esr1 as a key regulator of adolescent transcriptional dynamics in the MPOA of both sexes (a comprehensive list of sexually dimorphic genes in the Vgat +Esr1+population is provided in Supplementary file 11).

Figure 7 with 2 supplements see all

Download asset Open asset

Gene regulatory networks (GRNs) are complex systems of genes, TFs, 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 TFs (Esr1-TFs; Figure 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; Figure 8b–c) – indicating that Esr1 directly regulates many target genes. Additionally, nine female and six male Esr1-DEGs encoded TFs that function as central regulatory hubs within the network (Figure 8d–e), suggesting that Esr1’s influence extends indirectly through these secondary TFs. 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; Figure 8d–e and Figure 8—figure supplement 1).

Figure 8 with 1 supplement see all

Download asset Open asset

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 Materials and methods) – three Esr1-TFs: Zmat4 (sex-shared), Pou6f2 (male-specific), and Creb3l1 (female-specific), and six Esr1-DEGs regulated by Esr1 or Esr1-TFs: Gria3, Pdzrn4, Ctnna2, Nefl, Tead1, and Esr1. Consistent with our scRNAseq analysis (Figure 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 (Figure 8f–g). These results further suggest that Esr1 orchestrates sex-shared and -specific transcriptional states via an Esr1-specific GRN. These dynamic GRNs may in turn dictate the sex-shared and -specific adolescent changes that facilitate the maturation of sexually dimorphic neuronal circuits for mating.

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 AR-related genes and the identification of sexually dimorphic Esr1-dependent TFs (e.g. Pou6f2 in males, Creb3l1 in females) further highlight how Esr1 integrates hormonal signals to shape sex-specific neural circuits (Knoedler et al., 2022; Xu et al., 2012).

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 (Oomura et al., 1988; Wei et al., 2018; Wu et al., 2014; Ishii et al., 2024). 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 (Karigo et al., 2021; Moffitt et al., 2018). 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 (Hrvatin et al., 2020; Shimogori et al., 2010). 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 (Gegenhuber et al., 2022). 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.

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.

1. Although HM-HCR experiments showed the bidirectional control of transcriptional progression during adolescence, it is unclear if the facilitation in males by testosterone supplement is via activation of AR or Esr1 or both because testosterone will likely be converted to estrogen in the brain. Future studies using dihydrotestosterone (DHT) and estrogen in males may address this issue.

2. Although we have identified hormone/Esr1-dependent transcriptional trajectories during adolescence, the relations and interplay with genetically determined perinatal events, which are earlier and robust, are unclear. Some sex differences during adolescence might be an extension of perinatally established sex differences, while others might be unique adolescent changes.

3. While we have observed a robust effect of Esr1 KO in the scRNAseq experiment, which was further validated with the FISH experiment, it is possible that there are further heterogeneous Vgat-Esr1 populations in the MPOA which might be differentially targeted in each virally injected sample. To mitigate this, three to four mice were pooled for each sample in the scRNAseq experiment and HCR-FISH experiment. In addition to confirming recombinase RNA expression within the MPOA, we included samples with robust Esr1 deletion in the MPOA. Interestingly, due to the technical challenge, Esr1 deletion tends to be more robust than weakly detected recombinase RNA expression (data not shown).

4. While we have observed robust transcriptional progression in Vgat +Esr1+neurons during adolescence, we observed more mild alterations in Vglut2 +neurons. Although the scale of our study is comparable to or exceeds prior scRNAseq studies in MPOA (Moffitt et al., 2018; Hrvatin et al., 2020), future larger studies may have more sensitivity to detect adolescent transcriptional dynamics in Vglut2 +neurons.

5. Although we demonstrated adolescent transcriptional changes were observed as early as P35, and either hormonal deprivation or Esr1 KO prior to adolescence prevented the transcriptional progression (arrested transcriptional state even at adult), given the viral incubation time and permanent deletion of Esr1 after viral injection, it is challenging to disambiguate the role of Esr1 during adolescence and adulthood. Future studies injecting the virus in adults may provide additional insights on the similarity and difference between transcriptional changes during puberty and maintained transcriptional states in adulthood.

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 (Fang et al., 2018; Brock et al., 2015; Romanov et al., 2020). 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 (Blakemore et al., 2010; McCarthy and Arnold, 2011). While our study focused on Esr1, it is likely that other TFs and hormone receptors synergistically contribute to adolescent transcriptional dynamics. For example, the AR may play a role in males, while additional regulatory networks may be governed by sex chromosome-linked genes (Xu et al., 2012). Advanced genomic approaches, such as perturb-seq, could help unravel the complexity of these gene regulatory networks in a cell-type-specific manner (Aibar et al., 2017). 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 (Gegenhuber and Tollkuhn, 2020; McCarthy et al., 2012; Shansky, 2019).

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/6 J 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 hr 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.

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 s after every 60 nL infusion, and the capillaries were left in place for 10 min 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 SW foster dams. SW dams effectively nursed surgically treated pups, whereas C57BL/6 dams frequently committed infanticide following such procedures.

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- to 8-week-old C57BL/6 J female mice were ovariectomized using the same procedure. After a 3-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 hr and 5 µg at 24 hr before the assay) and Progesterone (Sigma, 500 µg) administered 6 hr prior to the RI assay.

The split-Cre design used in this study was previously described (Chen et al., 2018; Heymann et al., 2020). 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 (Heymann et al., 2020). 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 hr later. Purification of AAV1 capsids was performed using gradient centrifugation. Viral titer was estimated at approximately 2x10¹² particles/mL based on densitometry analysis following gel electrophoresis, using a known standard for reference.

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; Figure 1—figure supplement 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 (Figure 1—figure supplement 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 Figure 1 and Figure 1—figure supplement 1.

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 (Hashikawa et al., 2020). 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% O₂ 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 min. MPOA tissues were dissected under a light microscope (Leica MZFL) using a micro scalpel (Feather) from three to four animals per group (10–15 slices pooled per group). The tissue was enzymatically digested with 1 mg/mL pronase (Sigma-Aldrich) for 30 min 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–1000 cells/µL. Typically, 30,000–80,000 cells were collected from three to four 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 (10 x 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 10 x Genomics Cell Ranger pipeline (V3) to generate cell-by-gene count matrices for downstream analysis.

HM-HCR FISH assays were conducted to validate trajectory inference from scRNAseq data (HCR1; Figure 4 and Figure 4—figure supplements 1–2) and to cross-validate the Esr1-GRN analysis (HCR2; Figure 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 P23 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; Figure 1c and f) and reports from the pituitary (P25–30; Mayer et al., 2010). 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 unilaterally injected into the MPOA of Esr1^lox/lox mice of both sexes at P30–35, and brain tissues were harvested 3 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 von Buchholtz et al., 2021. Tissue sections were fixed in pre-chilled 4% paraformaldehyde (PFA) for 30 min on ice, rinsed twice in 1 x PBS at room temperature, dehydrated in ethanol (50%, 70%, 100%, 100%; 5 min each), and air-dried for 5 min. Sections were then permeabilized with protease IV (ACD, 322336) for 5 min at room temperature, followed by rinses in 1 x PBS and 2 x 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 min 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 5 x SSCT (75%, 50%, 25%; 15 min each), then in 5 x SSCT at 37 °C (15 min) and at room temperature (5 min).

For signal amplification, tissues were equilibrated in amplification buffer for 30 min 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 5 x SSCT at room temperature, followed by two 30-min washes in fresh 5 x SSCT and rinses in 2 x SSC. Autofluorescence was minimized with a quenching kit (Vector Laboratories, SP-8400–15). Sections were treated with quenching solution for 2 min, rinsed in 2 x SSC, counterstained with DAPI (ACD, 320858) for 30 s, 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 hr.

Probes and amplifiers were stripped between rounds to enable multiplexing. Coverslips were floated off in 2 x SSC for 30 min at room temperature, and tissue was incubated in DNase I (250 U/mL in 1 x DNase I buffer, Roche, 04716728001) for 1.5 hr at room temperature, followed by six washes in 2 x SSC (5 min 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 (full probe list provided in Supplementary file 12).

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 (Figure 3—figure supplement 2k,m). For HCR2, Esr1-TFs and Esr1-DEGs regulated by Esr1-TFs or Esr1 were selected from the Esr1-GRN (Figure 8). To disambiguate the MPOA and adjacent brain regions, quantitative analysis is restricted to Vgat +Esr1+neurons and is devoid of posterior BNST. For HCR2, AAV was injected unilaterally so that successful targeting of the MPOA with AAV-Cre-YFP (detection of recombinase RNA within the MPOA) and the deletion of Esr1 were confirmed for inclusion of samples.

AAV-frtFlex-Cre (300 nl) was unilaterally injected into the MPOA of VgatFlp or Vglut2Flp mice (n=2 per experiment). After 2 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 20 x 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 (Figure 1 and Figure 1—figure supplement 1).

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 BPS in males or VO in females. For female subjects, vaginal smears were also collected daily following VO to monitor estrous cycles. Behavioral experiments, including the 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-min 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 min. To avoid ejaculation, females were separated from the male after intromission began or approximately 5 s 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 (Cheong et al., 2015), the male subjects in this study, with Esr1 knockouts limited to MPOA Vgat +neurons during adolescence, showed no aggressive behaviors, consistent with previous findings (Sano et al., 2013).

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-min habituation period without stimuli, a social stimulus and an object were introduced, and the subject mouse was allowed to explore for 10 min. 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 min. 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 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-min 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.

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% PFA in PBS. Brains were extracted and post-fixed in 4% PFA in PBS on ice for 5–7 hr 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 IHC. Sections were first washed in PBS (3×5 min) and blocked with 15% normal donkey serum (NDS) in PBST (0.3% Triton X-100) for 2 hr 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 hr at 4 °C. After incubation, sections were washed in PBST (0.3% Triton X-100, 3×30 min) 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 hr at room temperature. Sections were then washed in PBST (2×15 min), incubated with DAPI for 2 min, rinsed in PBS (2×15 min), mounted on slides, and coverslipped with mounting medium. Tiled images of the MPOA were acquired using a Zeiss ApoTome2 fluorescence microscope with a 20 x 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). 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.

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 :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 GRN deconstruction.

Clustering and differential gene expression analyses were performed using the Seurat V3 package (Stuart et al., 2019), while LISI (Korsunsky et al., 2019) was applied to evaluate integration performance across conditions. Lineage inference, representing transcriptional connectivity of neuronal clusters, was conducted with PAGA (Wolf et al., 2019a). Single-cell gene set activity was quantified using AUCell (Aibar et al., 2017). GRNs were inferred and deconstructed using the SCENIC package, which integrates GENIE3 (Aibar et al., 2017; Huynh-Thu et al., 2010) and RcisTarget. These analyses depicted cis-regulatory relationships between TFs and DEGs. Gene co-expression networks were identified using scWGCNA (Morabito et al., 2021). To examine transcriptional state progressions by age and hormonal state, Monocle V3 (Cao et al., 2019) and MELD (Burkhardt et al., 2021) were applied. Ontology analysis was conducted using the online Enrichr platform (Chen et al., 2013; Kuleshov et al., 2016).

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.

Preprocessing followed the procedures described in our previously published study (Hashikawa et al., 2020). 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; DePasquale et al., 2019b) with default settings, applying a doublet rate of 5.6 ± 1.0%.

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 (Stuart et al., 2019). This method combines canonical correlation analysis (CCA; Butler et al., 2018) and mutual nearest neighbor analysis (Haghverdi et al., 2018). 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 (Hashikawa et al., 2020; Tasic et al., 2018; Zeisel et al., 2018). Neuronal and non-neuronal cell types were combined for visualization (Figure 2—figure supplement 1a–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.

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 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; Korsunsky et al., 2019) was computed in UMAP space, with LISI scores reflecting the effective number of distinct groups in the local neighborhood of each cell (Figure 2—figure supplement 1i). Clustering stability was further evaluated by examining lineage relationships between clusters at varying resolutions (0.2–2.0) using the clustree package (Zappia and Oshlack, 2018). Clustering tree analysis showed the emergence of heterogeneous Esr1 clusters at resolution 0.8 or higher (Figure 2—figure supplement 1k). Specific marker genes in neuronal clusters (Figure 2—figure supplement 2a) and Vgat +Esr1+subclusters (Figure 2—figure supplement 2h) further validated clustering.

To explore transcriptional connectivity between neuronal clusters, PAGA (Wolf et al., 2019a) 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 (Figure 2c).

Seurat V3 was used to jointly analyze publicly available POA MERFISH data (Moffitt et al., 2018) and our scRNAseq data. To establish correspondence between MERFISH and scRNAseq clusters related to reproductive behaviors in adult mice, cells from behaviorally naive 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 naive 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 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.

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 HA-DEGs (Figure 2g), Ar-regulon in male-dimorphic genes (Figure 5h–j), and adolescence-related versus unrelated dimorphic genes (Figure 6g–h).

Monocle 3 was used to quantify adolescent transcriptional progression in scRNAseq and HM-HCR FISH datasets (Cao et al., 2019). 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 HA-DEGs (using the top 10–15 principal components [PCs]). 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 (Figure 3—figure supplement 1 and Figure 7—figure supplement 2). Branch-specific DEG analysis identified genes enriched in each branch and those shared across branches (Figure 3—figure supplement 1 and Figure 7—figure supplement 2).

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; Figure 3—figure supplement 2k, 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 (Figure 4—figure supplement 2). 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.

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; Burkhardt et al., 2021; Moon et al., 2019). 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. 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; Figure 3—figure supplement 2h–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 (Figure 3—figure supplement 2h–i, j,l).

Data from Esr1KO groups were independently processed for quality control, including the removal of low-quality cells and doublets (DePasquale et al., 2019b). Seurat V3 was used to perform clustering and identify neuronal clusters (detailed in Integrative Clustering and Differential Gene Expression Analysis). The dataset included 6727 cells from Esr1KOF (female) and 6464 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 (Figure 7—figure supplement 1; Esr1KOF: 2182 neuronal cells; Esr1KOM: 4761 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 (Figure 7—figure supplement 1). 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 (Figure 7—figure supplement 1). 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. PCs 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).

To infer the cis-regulation of genes by TFs, we used the SCENIC computational framework (Aibar et al., 2017; Hashikawa et al., 2020; Davie et al., 2018). SCENIC was combined with cell-type-specific gene manipulation to causally infer the 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 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, nine (female) and six (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, Figure 8d–e, and Figure 8—figure supplement 1). 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 (scWGCNA) was employed to identify sex-specific co-expression networks in Vgat +Esr1+ cells (Morabito et al., 2021). 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 PC of gene expression for each module, were computed using the MEs 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 (Chen et al., 2013; Kuleshov et al., 2016). 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 was performed using Enrichr to examine the enrichment of functionally related genes within Esr1-DEGs (Figure 8—figure supplement 1). 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 Figure 8—figure supplement 1. 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.

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 SVM classification using the scikit-learn library. GridSearch with 10-fold cross-validation was applied to optimize the classification as described in previous studies (Otis et al., 2017; Rossi et al., 2019). SVM models were tested using both linear and radial basis function (rbf) kernels. Hyperparameters were optimized across a range of values for γ(10^–3, 10^–2, 10^–1, 10¹, 10², 10³) and C (10^–3, 10^–2, 10^–1, 10¹, 10², 10³). Input features for classification included 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).

Statistical analyses for each experiment are described in the Materials and 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.

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

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

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

Figure 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.

Figure 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.

Figure 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.

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

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

Figure 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.

Figure 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.

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

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

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

Figure 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.

Figure 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).

Figure 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).

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

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

Figure 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).

Figure 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).

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

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

Figure 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.

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

male:

• 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

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

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

Figure 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.

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

Figure 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.

Figure 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.

Figure 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.

Figure 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, Pp<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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

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

Figure 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.

Figure 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.

Figure 1—figure supplement 1a: paired t-test. t(11)=-0.523, p=0.612.

Figure 1—figure supplement 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.

Figure 1—figure supplement 1d left: unpaired t-test. t(24) = –2.966, p=0.0067.

Figure 1—figure supplement 1d right: unpaired t-test. t(15.466)=–3.7502, p=0.0018.

Figure 1—figure supplement 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.

Figure 1—figure supplement 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.

Figure 1—figure supplement 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.

Figure 1—figure supplement 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.

Figure 1—figure supplement 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.

Figure 1—figure supplement 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.

Figure 2—figure supplement 1i: unpaired Wilcoxon test. W=396364535, p<0.0001.

Figure 2—figure supplement 2d 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.176780, 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.

Figure 2—figure supplement 2d 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.

Figure 2—figure supplement 2e 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.

Figure 2—figure supplement 2e 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: ANOVA F (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, Pp<0.0001. P35 vs P50, t(55) = 7.715, p<0.0001.

Figure 2—figure supplement 2f: Linear regression analysis. Adjusted R-squared=0.43, p=0.013.

Figure 2—figure supplement 2g: Linear regression analysis. Adjusted R-squared for each hormone receptor gene was reported.

Figure 2—figure supplement 3a, 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.

Figure 3—figure supplement 1c 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.

Figure 3—figure supplement 1c 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).

Figure 3—figure supplement 1c 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.

Figure 3—figure supplement 1c 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).

Figure 3—figure supplement 1c 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.

Figure 3—figure supplement 1c 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.

Figure 3—figure supplement 1d 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.

Figure 3—figure supplement 1d 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).

Figure 3—figure supplement 1d 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.

Figure 3—figure supplement 1d 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).

Figure 3—figure supplement 1d 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.

Figure 3—figure supplement 1d 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.

Figure 3—figure supplement 1e 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.

Figure 3—figure supplement 1e 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.

Figure 3—figure supplement 1f 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.

Figure 3—figure supplement 1f 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.

Figure 3—figure supplement 1g branch1: unpaired t-test. T(152.39)=216.35, p<0.0001.

Figure 3—figure supplement 1g branch2: unpaired t-test. T(147.5)=153.52, p<0.0001.

Figure 3—figure supplement 1h branch1: unpaired t-test. T(130.85)=244.88, p<0.0001.

Figure 3—figure supplement 1h branch2: unpaired t-test. T(124.65)=204.9, p<0.0001.

Figure 3—figure supplement 1i branch1 AUCell score: unpaired Wilcoxon test. p<0.0001.

Figure 3—figure supplement 1i branch2 AUCell score: unpaired Wilcoxon test. <0.0001.

Figure 3—figure supplement 1i branch1 SVM: unpaired t-test. T(134.62)=142.72, p<0.0001.

Figure 3—figure supplement 1i branch2 SVM: unpaired t-test. T(176.75)=66.219, p<0.0001.

Figure 3—figure supplement 1j branch1 AUCell score: unpaired Wilcoxon test. p<0.0001.

Figure 3—figure supplement 1j branch2 AUCell score: unpaired Wilcoxon test. p=0.0050.

Figure 3—figure supplement 1j branch1 SVM: unpaired t-test. T(164.27)=83.34, p<0.0001.

Figure 3—figure supplement 1j branch2 SVM: unpaired t-test. T(190.61)=71.154, p<0.0001.

Figure 3—figure supplement 2d: 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.

Figure 3—figure supplement 2f: 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.

Figure 3—figure supplement 2j: Kruskal-Wallis rank sum test followed by multiple comparisons. Kruskal-Wallis chi-squared=570.12, df = 3, p<0.0001.

Figure 3—figure supplement 2k: unpaired t-test. T(198) = 61.044, p<0.0001.

Figure 3—figure supplement 2l: 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.

Figure 3—figure supplement 2m: unpaired t-test. T(198) = 99.879, p<0.0001.

Figure 3—figure supplement 2o: 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.

Figure 3—figure supplement 2q: 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.

Figure 3—figure supplement 3d: 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.

Figure 3—figure supplement 3e: 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.

Figure 4—figure supplement 2a: 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.

Figure 4—figure supplement 2b: 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.

Figure 4—figure supplement 2d: 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.

Figure 4—figure supplement 2g: 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.

Figure 4—figure supplement 2n: unpaired t-test. T(198)=91.01, p<0.0001.

Figure 4—figure supplement 2q: unpaired t-test. T(198)=157.2, p<0.0001.

Figure 4—figure supplement 2r: 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.

Figure 4—figure supplement 2s: 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.

Figure 4—figure supplement 2u: 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.

Figure 5—figure supplement 1d: 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.

Figure 7—figure supplement 2c 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.

Figure 7—figure supplement 2c 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).

Figure 7—figure supplement 2c 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.

Figure 7—figure supplement 2c 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).

Figure 7—figure supplement 2c 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.

Figure 7—figure supplement 2c 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.

Figure 7—figure supplement 2d 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.

Figure 7—figure supplement 2d 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).

Figure 7—figure supplement 2d 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.

Figure 7—figure supplement 2d 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).

Figure 7—figure supplement 2d 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.

Figure 7—figure supplement 2d 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.

Figure 7—figure supplement 2e 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.

Figure 7—figure supplement 2e 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.

Figure 7—figure supplement 2f 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.

Figure 7—figure supplement 2f 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.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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 10 x Genomics library generation. 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.

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, protocol 4450-01) of the University of North Carolina and University of Washington.

1. Preprint posted: February 26, 2025
2. Sent for peer review: February 26, 2025
3. Reviewed Preprint version 1: May 1, 2025
4. Reviewed Preprint version 2: April 7, 2026
5. Version of Record published: June 19, 2026

You can cite all versions using the DOI https://doi.org/10.7554/eLife.106347. This DOI represents all versions, and will always resolve to the latest one.

© 2025, Hashikawa, Hashikawa et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.