Adolescence is a secondary critical period, following the perinatal steroid surge, during which juvenile animals undergo extensive physiological maturation of both the body and nervous system as they transition into adulthood


1

,

2


. This maturation is primarily orchestrated by the hypothalamus-pituitary-gonad (HPG) axis, which activates gonadal cells in the testes and ovaries to release testosterone and estrogen, initiated during adolescence


<sup>1

,

3

,

4</sup>


. These circulating sex steroids are essential not only for reproductive maturation but also for guiding the development and execution of sex-specific behaviors through their direct actions on the central nervous system


<sup>5

–

7</sup>


. Despite the well-established role of these hormones in shaping behavior, the molecular mechanisms underlying their influence on brain development during adolescence are still limited to brain-region level (bulk)


⁸


in humans and model organisms and adolescent transcriptional dynamics at single cell resolution in the brain remain poorly understood (but see a pioneering study in the human testis


⁹


).

The medial preoptic area (MPOA) is a critical region in the brain, known for its involvement in mating and other social behaviors


<sup>7

,

10

–

16</sup>


. Within the MPOA, specific molecularly defined subpopulations of neurons express various neuropeptides and hormone receptors that are closely tied to reproductive behavior


<sup>16

–

20</sup>


. Notably, the MPOA is enriched with the expression of steroid hormone receptor genes such as

Esr1,

or estrogen receptor 1


<sup>21

,

22</sup>


, leading us to hypothesize that gonadal steroid signaling during adolescence facilitates the transcriptional maturation of these genetically defined MPOA neurons, which are crucial for the development of mating behaviors. While recent advances have shed light on the molecular diversity of neurons in sexually dimorphic brain regions during the perinatal period


²³


and adulthood


<sup>22

,

24

,

25</sup>


, the transcriptional dynamics that occur during adolescence, particularly at the single-cell level, remain largely unexplored. This gap in our understanding represents a critical frontier in unraveling how hormonal signaling shapes the neural circuits that govern sex-specific behaviors.

In this study, we used single-cell RNA sequencing (scRNAseq) and hybridized chain reaction fluorescent

in situ

hybridization (HCR-FISH) to investigate the transcriptional dynamics of the MPOA throughout adolescence in both male and female mice. We show that

Esr1

plays a critical role in the transcriptional maturation of GABAergic neurons (Vgat+) within the MPOA, which is essential for the development of mating behaviors. Deleting

Esr1

in Vgat+ neurons disrupted the development of mating behaviors in both sexes, while deletion in glutamatergic neurons (Vglut2+) had no observable effect. Our findings highlight the importance of

Esr1

in guiding sex-specific transcriptional programs during adolescence, revealing distinct hormone-dependent gene expression patterns and gene-regulatory networks crucial for reproductive behaviors.

Results

Selective knockout of

Esr1

in GABAergic or glutamatergic MPOA neurons at pre-adolescence

Previous studies have shown that knocking out

Esr1

, either globally or in specific cell types, reduces mating behavior


<sup>7

,

26</sup>


. However, the role of

Esr1

in genetically defined MPOA neurons during adolescent development remains unclear, as these studies could not easily differentiate between effects occurring during perinatal and adolescent periods. To address this gap and directly test whether the emergence of mating behavior requires steroid hormone signaling in molecularly defined MPOA cell types, we developed an AAV-frtFlex-Cre virus. We injected this virus into the MPOA of

Slc32a1


^Flp

::

Esr1


^lox/lox

(MPOA

^Vgat-Esr1KO

) or

Slc17a6


^Flp

::

Esr1


^lox/lox

(MPOA

^{Vglut2-Esr1KO}

) male and female mice at postnatal day (P) 14-18. This viral strategy allowed for

Esr1

deletion in MPOA GABAergic or glutamatergic neurons, depending on the transgenic mouse line, throughout adolescence (

Fig. 1a-b,i-j

and

Supplementary Fig. 1a-b

). The onset of genital changes, the maturation of the HPG axis (as determined by the first day of estrus in female mice


²⁷


and first day of initiation of sexual behaviors in male mice


²⁸


), and locomotor behaviors were not altered in MPOA

^Vgat-Esr1KO

and MPOA

^{Vglut2-Esr1KO}

mice (

Fig. 1c-d,f,k-l,n

and

Supplementary Fig. 1e-f

). The establishment of mating behavior was abolished in MPOA

^Vgat-Esr1KO

mice in both sexes (

Fig. 1d-e,g-h

and

Supplementary Fig. 1a-b,g

), but unaffected in MPOA

^{Vglut2-Esr1KO}

mice (

Fig. 1l-m,o-p

and

Supplementary Fig. 1h

). These results extend prior work to show that

Esr1

signaling in MPOA GABAergic neurons is required for the emergence of mating behavior during adolescence

<sup>7,

14

,

15</sup>


.

Transcriptional dynamics of MPOA cell types during adolescent development

To characterize the transcriptional dynamics occurring in the MPOA specific to the adolescent period, we collected tissue from male and female wildtype mice at pre- (P23), mid- (P35) and post- (P50) adolescence timepoints, as well as from mice gonadectomized (GDX) prior to adolescence onset. We used droplet-based scRNAseq to recover the transcriptomes of 58,921 cells and combined the data from all samples across the 8 experimental groups (

Fig. 2a

and

Supplementary Fig. 2a-e

and Supplementary Table 1). Using known neuronal marker genes

Stmn2

and

Thy1

, we performed high-resolution clustering to focus our subsequent analyses on neuronal cell data (24,627 cells). Thirty-five neuron specific clusters were resolved with 434 marker genes in total (UMIs per cell: 5,395; median genes per cell: 2,592), resulting in 20 GABAergic clusters (Vgat+ clusters, 13,334 cells; 54.4% of neurons), 12 glutamatergic clusters (Vglut2+ clusters, 7,658 cells; 31.1% of neurons), and 3 clusters with mixed glutamatergic and GABAergic markers (3,635 cells; 14.8% of neurons) (

Fig. 2b

and

Supplementary Figs. 2f-l, 3a,j

and Supplementary Table 2), largely consistent with previously published data


<sup>22

,

29</sup>


. To further quantify transcriptional associations between MPOA neuronal cell types, we performed partition-based graph abstraction (PAGA) analysis


³⁰


to identify cell cluster similarities through a graphical representation of cell differentiation (

Fig. 2c

). Consistent with their distinct roles, Vgat+ and Vglut2+ neurons were readily separated in the PAGA graph. Additionally, Vgat+ clusters 2, 4 and 16, all of which expressed

Esr1

, showed higher connectivity compared to other Vgat+ cell type clusters (

Fig. 2c

and

Supplementary Figs. 2m

,

3h-i

), suggesting that Vgat+Esr1+ clusters form a distinct transcriptional unit.

Irrespective of age, sex, and hormonal states, we found that each identified cell type was represented across all 8 experimental groups (

Supplementary Figs. 2d,h-i

,

3a

), indicating that new cell types do not emerge during adolescence. However, when assessing gene expression differences across timepoints and sex, we found that hormone-associated differentially expressed genes (HA-DEGs) showed substantial variability across neuronal cell types, ranging from 0 to 88 HA-DEGs per cell type (

Fig. 2d-f

and

Supplementary Fig. 3a-g

and Supplementary Tables 3-5). We then used these HA-DEGs in regression, co-expression, and AUCell analyses to identify the neuronal types that are transcriptionally sensitive at different hormonal states


³¹


. Regression analysis revealed a strong enrichment of HA-DEGs in Vgat+Esr1+ neurons (

Fig. 2e-f

). Consistent with this result, co-expression analysis demonstrated that out of more than 1,000 transcription factors (TFs),

Esr1

was one of the most co-expressed in both sexes (

Fig. 2f

and Supplementary Table 6). Furthermore, in predominantly Vgat+Esr1+ clusters, the aggregate expression of HA-DEGs (AUCell) displayed patterns sensitive to both adolescence and hormonal changes, with the highest expression observed in adult mice (P50) and notably lower expression at P23 and during hormonal deprivation (GDX) (

Fig. 2g

). Together, these analyses demonstrate that sex-hormone signaling during adolescence development significantly influences the transcriptional states of specific MPOA neuronal cell types, particularly cells belonging to Vgat+Esr1+ clusters.

Having identified the distinct transcriptional sensitivity of Vgat+Esr1+ neurons to hormone signaling during adolescence, we next examined whether these neurons might play a direct role during behavior. Given that mating behavior was abolished in MPOA

^Vgat-Esr1KO

mice in both sexes (

Fig. 1

), this finding, along with our scRNAseq analyses, suggests that Vgat+Esr1+ neurons are critically involved during mating behaviors. To test this hypothesis, we conducted an integrative analysis using our scRNAseq dataset and a publicly available multiplexed error-robust fluorescence

in situ

hybridization (MERFISH) dataset


²²


which identified social behavior-activated POA cells via mating, parenting, and fighting behavior assays using

Fos

(

Fig. 2h

). We clustered the neuronal cells in the MERFISH data and computed the enrichment of

Fos

in the derived clusters. Consistent with previous observations, only a subset of MERFISH Gad1

+

inhibitory clusters showed

Fos

enrichment (

Supplementary Fig. 4a,b

). After integrating our scRNAseq clusters with the MERFISH clusters (Details in Methods. Seurat V3 Integrative Analysis: reference-based integration; label transfers


³²


) we were able to identify mating, parenting, and fighting-relevant neuronal clusters within our dataset (

Fig. 2h

and

Supplementary Fig. 4c-d

and Supplementary Table 7). We found that

Esr1

was the most enriched and selective gene in the mating-relevant scRNAseq Vgat+ clusters (

Fig. 2h-i

and Supplementary Table 7). This observation is consistent with a previous study demonstrating that optogenetic stimulation of Vgat+Esr1+ neurons in the MPOA is sufficient to induce sexual behaviors in male mice


¹⁸


. Collectively, these findings suggest that a subset of Vgat+ clusters expressing

Esr1

, which exhibit transcriptional changes during adolescence, are also functionally linked to mating behaviors in both males and females.

While these differential gene expression (DE) analyses quantify changes in individual gene expression at the pseudo-bulk level (

Fig. 2d

and

Supplementary Fig. 3b-c

), it does not capture the transitions in transcriptional states of individual cells as they progress through distinct biological stages. These transitions, driven by complex interactions between multiple genes, are key to understanding how single cells change over time


³³


. To resolve these dynamic changes, we applied pseudotime analysis (

Fig. 3

and

Supplementary Figs. 5-6

), a method that infers the transcriptional progression of individual cells through a biological process, such as adolescence, by constructing a principal graph based on combinatorial gene expression patterns


<sup>34

,

35</sup>


.

Recognizing the potential for sex-specific differences in transcriptional progression, we created separate pseudotime manifold models for males and females to capture nuanced dynamics during adolescence


³⁶


. In both sexes, Vgat+Esr1+ trajectories revealed that transcriptional states at P35 and P50 branch and diverge from preadolescence (

Fig. 3a-f

and

Supplementary Fig. 5a-f

), suggesting an acceleration of transcriptional dynamics between these timepoints. Neurons from GDX mice, on the other hand, showed arrested transcriptional progression, closely resembling preadolescent states, demonstrating the importance of circulating sex hormones in the maturation of these cells during adolescence. In contrast, Vgat+ cells lacking steroid hormone receptor gene expression (Vgat+ hormone R

^Low

) exhibited minimal transcriptional changes in response to hormone changes (

Fig. 3g,i

and

Supplementary Fig. 6a-c

). Pair-wise DEG analysis consistently showed that larger number of DEGs between P35 and P23 in Vgat+Esr1+ (male: 146 genes; female: 162 genes) than Vgat+ hormone R

^Low

(male: 26 genes; female: 1 gene). Furthermore, all Vgat+Esr1+ clusters were found to co-express the androgen receptor gene (

Ar

) to a high degree (68.3 – 88.3 %,

Supplementary Fig. 6a

). However, Vgat+ clusters expressing

Ar

but not

Esr1

(Vgat+ Esr1- AR+) showed only a moderate transcriptional progression in males and minimal changes in females (

Supplementary Fig. 6b-d

), suggesting that

Esr1

is more influential in driving transcriptional changes during adolescence than

Ar

alone. Lastly, we analyzed Vglut2+ clusters enriched with

Esr1

, and again found only moderate transcriptional progression in both sexes (

Supplementary Fig. 7

), altogether indicating that cell types outside of Vgat+Esr1+ clusters are less transcriptionally dynamic during this period.

Two major branches of transcriptional trajectories emerged from subpopulations of Vgat+Esr1+ cells in both sexes. In females, Branch1 was dominated by Vgat+ cluster 4, while Branch2 was occupied by clusters 2 and 16. In males, Branch1 was composed of Vgat+ clusters 4 and 16, while Branch 2 was primarily formed by Vgat+ cluster 2 (

Supplementary Fig. 5g-j

). The analysis of DEGs along these branches revealed distinct gene enrichment, indicating that specific subpopulations of Vgat+Esr1+ cells undergo unique transcriptional progressions during adolescence (

Supplementary Fig. 5c-d

and Supplementary Table 8). When we combined the trajectories for both sexes into a joint manifold, the model failed to capture hormone-dependent dynamics (

Supplementary Fig. 6e-g

), emphasizing that separate models are necessary for accurate analysis of sex-specific and sex-shared transcriptional progressions (

Fig. 3k-n

).

To strengthen our understanding of the developmental trajectory and transcriptional maturity of Vgat+Esr1+ clusters, we employed a support vector machine (SVM) classifier to predict the developmental state of single cells based solely on HA-DEGs


³⁷


. The SVM accurately classified Vgat+Esr1+ single cells as either transcriptionally mature (intact adolescent or adult) or immature (preadolescent or hormonally deprived) with high accuracy (median prediction accuracy 90.6 – 93.6 %). In contrast, Vgat+ hormone R

^Low

cells had significantly lower prediction accuracy (median 66.7 - 73.0 %) (

Fig. 3h,j

), further underscoring the importance of Vgat+Esr1+ clusters in the development of MPOA transcriptional states.

Finally, we validated our pseudotime trajectories (Monocle V3) using Manifold Enhancement Latent Dimension (MELD) analysis, which measures continuous transcriptional progression in response to experimental conditions (P23, P35, P50, GDX)


<sup>38

,

39</sup>


. MELD analysis corroborated our findings from pseudotime, showing that intact adolescents and adults were spatially segregated in transcriptional space and exhibited a higher likelihood of reaching mature transcriptional states compared to preadolescent and hormonally deprived groups (

Supplementary Fig. 6h-m

). These results, visualized in Potential of Heat-diffusion for Affinity-based Trajectory Embedding (PHATE) space


³⁸


, further highlight that while MPOA neuronal cell types are terminally diversified by preadolescence (P23) (

Supplementary Figs. 2h

,

3e

), Vgat+Esr1+ neurons continue to undergo hormonally dependent transcriptional refinement from adolescence into adulthood.

Spatial phenotyping of MPOA cell types during adolescence

To assess whether adolescent hormonal changes affect spatially resolved gene expression in the MPOA and to cross-validate our scRNAseq trajectory analyses, we performed highly multiplexed hybridization chain reaction fluorescent

in situ

hybridization (HM-HCR FISH, V3)


<sup>40

,

41</sup>


. This technique allowed us to detect transcripts of ∼12 genes across 41,549 MPOA cells at single-molecule resolution (

Fig. 4a-c

and

Supplementary Figs. 8-9

). The rationale for selecting this gene panel is related to

Supplementary Fig. 6k,m

, with details provided in the Methods. In addition to the experimental groups used in the scRNAseq experiments, we prepared tissue from hormonally supplemented mice. These mice received testosterone (males) or estrogen (females) from P23-P27 (see timeline in

Fig. 4a

) to investigate whether early sex steroid supplementation accelerated adolescent transcriptional trajectories.

To distinguish Vgat+Esr1+ from Vgat+ hormone R

^Low

cells, we measured

Slc32a1

,

Esr1, Ar,

and 8-9 DEGs during adolescence (

Fig. 4d-e

and

Supplementary Figs. 8-9

. Details in Methods). Consistent with the scRNAseq data, we observed a transcriptional progression from the preadolescent state (P23) to a matured state (P35 and P50), where the dynamics were bidirectionally influenced by circulating steroids– accelerated by early testosterone or estrogen supplementation (P28 T or E2) and delayed in a hormonally-deprived state (GDX) (

Fig. 4f-g

and

Supplementary Fig. 9a-c,f,i-q

). Interestingly, the transcriptional trajectories of female and male adolescents differed spatially, but both were regulated by sex hormones (

Supplementary Fig. 9r-u

). This suggests that neurons defined transcriptionally as adult female or male occupy partially overlapping spatial distributions in the MPOA. SVM classification revealed that the expression of ≥10 genes in Vgat+Esr1+ cells from the HM-HCR FISH data was sufficient to accurately classify individual cells as transcriptionally mature or immature, with high accuracy (84.4–85.5 %). However, the same set of genes significantly underperformed in classifying Vgat+ hormone R

^Low

cells (63.6–64.6 %) (

Supplementary Fig. 9d,g

). Classification accuracy further decreased when individual genes were iteratively removed (

Supplementary Fig. 9e,h

), further indicating that a relatively small combinatorial set of HA-DEGs is sufficient for accurately identifying age- and sex-specific transcriptional states.

Molecular profiling of the MPOA in adult mice has established the presence of sexual dimorphism


<sup>42

,

43</sup>


. Our scRNAseq and

in situ

analyses further show hormone-dependent transcriptional progression during adolescence in both females and males. However, critical questions remain about whether (1) sexual dimorphism is evident within specific MPOA cell types, (2) sexually dimorphic genes overlap with adolescent gene sets, and (3) the degree of sexual dimorphism shifts during adolescence. To address these questions, we first identified sexually dimorphic genes across Vgat+ clusters (

Supplementary Fig. 10a and 11

). Like HA-DEGs, sexually dimorphic genes were enriched within Vgat+Esr1+ clusters. The number of these dimorphic genes correlated with

Esr1

expression, where sexually dimorphic genes co-expressed most often with

Esr1

in both sexes (

Fig. 5a-e

and

Supplementary Fig. 10a

and Supplementary Table 6). Notably, only subsets of sexually dimorphic genes were also classified as adolescent genes (

Figs. 5f

,

6a-b

) and some of these adolescent genes were shared across sexes (

Fig. 5g

).

Previous research has shown that

Esr1

plays a role in regulating sexual dimorphisms within sexually dimorphic brain nuclei


²³


. However, androgen receptor (AR) activation may also be essential for the expression of sexually dimorphic genes. To determine whether AR regulates male dimorphic genes, we performed single-cell regulatory network inference (SCENIC) analysis


³¹


. SCENIC identified a regulon of 33 AR-regulated genes, which were highly co-expressed with

Ar

and enriched with

Ar

’s consensus DNA regulatory element within their gene loci, in the male dimorphic gene set (

Fig. 5h

). AUCell analysis, however, showed that male dimorphic AR-regulon genes were enriched within male Vgat+Esr1+ clusters but not in Vgat+ Esr1- Ar+ clusters (specifically Vgat+ 5, 8, and 13) (

Fig. 5i

). Consistent with this observation, regression analysis showed that

Esr1

expression, rather than

Ar

, predicted the expression of male dimorphic AR-regulon genes (

Fig. 5j

).

As previously noted, a subset of sexually dimorphic genes overlaps with adolescent-related genes (

Fig. 5f

). The number of sexually dimorphic genes increases during adolescence but reverts following gonadectomy (

Fig. 6a-b

). Among 45 female dimorphic genes, 12 showed increased expression during adolescence (adolescence-related dimorphic genes), while 7 were not associated with adolescence (

Fig. 6a, right

). Of the 110 dimorphic genes in males, 52 were adolescence-related dimorphic genes and 10 were not (

Fig. 6b, right

). The expression levels of adolescence-related dimorphic genes were the highest at P50 and were again reduced by gonadectomy in both sexes (

Fig. 6a-b

). SVM analysis revealed that adolescence-related dimorphic genes decoded sex with the highest accuracy in intact P50 males and females, followed by P35, and were less accurate in P23 and GDX groups. In contrast, adolescence-unrelated dimorphic genes successfully predict the sexes of single cells with high accuracy (>90%) irrespective of age and conditions (

Fig. 6c-d

), emphasizing the influence of adolescence-related dimorphic genes in shaping MPOA sex-specific transcriptional profiles. AUCell and trajectory analyses were then performed to quantify the transcriptional dynamics of sexual dimorphic genes during adolescence. These analyses consistently demonstrated that transcriptional states were the most sexually dimorphic at P50 driven by adolescence-related dimorphic genes and were the least dimorphic at P23 and GDX (

Fig. 6e-h

). Thus, combined trajectory analysis of all Vgat+Esr1+ neurons from both sexes, along with AUCell analysis, revealed that P50F and P50M cells were the most sexually dimorphic from each other. However, these transcriptional states were bridged via a common immature state observed largely prior to adolescence onset or in hormonally depleted conditions, suggesting that sexually dimorphic MPOA gene expression is bimodal in adults but largely continuous between males and females prior to adolescence. These analyses (

Figs. 5

and

6

) indicate that: 1) sexual dimorphism is most pronounced in Vgat+Esr1+ cells; 2) a subset of sexually dimorphic genes are linked to adolescence, with a fraction of adolescence-related genes shared between sexes; and 3) sexual dimorphism in the MPOA expands during adolescence as sex steroid hormone levels rise.

Esr1

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

Co-expression of

Esr1

with adolescence genes and dimorphic genes, along with partially distinct transcriptional dynamics in Vgat+Esr1+ cells between adolescent females and males, suggests that

Esr1

uniquely regulates gene expression in each sex. To further test this, we virally deleted

Esr1

from Vgat+ MPOA neurons in female and male mice prior to the onset of adolescence and conducted scRNAseq at P50, comparing Esr1KO to Esr1-intact controls (

Fig. 7a

and

Supplementary Fig. 11.

Details in Methods). Esr1KO in Vgat+ MPOA neurons reduced the DEGs in males and females (Esr1-DEGs) to 600 and 824, respectively, while only reducing 5.3% and 3.3% of male and female Esr1-DEGs in Vgat+ hormone R

^Low

cells (

Fig. 7b

and

Supplementary Fig. 10c

and Supplementary Table 9). Previously identified HA-DEGs (

Fig. 2

) were largely represented in this list of Esr1-DEGs, suggesting that

Esr1

deletion recapitulates hormone deprivation effects on the transcriptomes of Vgat+Esr1+ cells (males: 56.5%, 78/138 genes; females: 73.3%, 124/169 genes;

Fig. 7b

). Again, because of the complex nature of this relationship between dimorphic and adolescent genes, we generated a separate manifold for each sex to assess MPOA adolescent transcriptional dynamics. Pseudotime trajectory analysis indicated that Esr1KO near-completely prevented adolescent transcriptional maturation in both sexes (

Fig. 7c-h

and

Supplementary Fig. 12

). DE analysis also identified sex-shared and -specific processes (

Fig. 7i-l

and Supplementary Table 10). In addition, dimorphic genes in the Esr1-DEGs of Vgat+Esr1+ cells were sufficient to predict the sex in P50 females and males with the highest accuracy (81.4 ± 0.11 %), followed by P35, and lowest accuracy in P23 (

Supplementary Fig. 10d

). Like our previous observation (

Supplementary Fig. 6

), we observed two major branches of transcriptional trajectories with branch specific gene expression profiles, primarily originating from one to two subpopulations of Vgat+Esr1+ clusters (

Supplementary Fig. 12c-d

and Supplementary Table 10) highlighting

Esr1

as a key regulator of adolescent transcriptional dynamics in the MPOA of both sexes.

Gene regulatory networks (GRNs) are complex systems of genes, transcription factors, and other molecules that interact to control gene expression within a cell. To identify GRNs that are shared between sexes or specific to one sex, we used functional-SCENIC to examine if

Esr1

cis-regulates its own differentially expressed genes (Esr1-DEGs) and/or transcription factors (Esr1-TFs) (

Fig. 8a

). Our analysis revealed that a significant portion of Esr1-DEGs contained

Esr1

binding sites–14.9% in females (123 out of 824 genes) and 16.3% in males (98 out of 600 genes) (

Fig. 8b-c

)–indicating that

Esr1

directly regulates many target genes. Additionally, 9 female and 6 male Esr1-DEGs encoded transcription factors that function as central regulatory hubs within the network (

Fig. 8d-e

), suggesting that Esr1’s influence extends indirectly through these secondary transcription factors. 54.1% of female and 28.3% of male Esr1-DEGs were cis-regulated by Esr1-TFs, and collectively, 56.1 % female and 36.7% male Esr1-DEGs were cis-regulated by

Esr1

, Esr1-TFs, or a combination of the two (Esr1-GRNs;

Fig. 8d-e

and

Supplementary Fig. 13

).

To validate these Esr1-GRNs

in situ

, we virally deleted

Esr1

from MPOA neurons in female and male mice prior to the onset of adolescence and conducted HM-HCR FISH on brain sections at P50, comparing Esr1KO to Esr1-intact controls. The gene panel included genes identified via SCENIC (Details in Methods)–3 Esr1-TFs:

Zmat4

(sex-shared),

Pou6f2

(male-specific), and

Creb3l1

(female-specific), and 6 Esr1-DEGs regulated by

Esr1

or Esr1-TFs:

Gria3

,

Pdzrn4

,

Ctnna2

,

Nefl

,

Tead1

, and

Esr1

. Consistent with our scRNAseq analysis (

Fig. 8b-e

), sex-shared Esr1-TF,

Zmat4

, was reduced by Esr1KO in both sexes, while male-specific Esr1-TF,

Pou6f2

, was decreased only in males. Female-specific Esr1-TF,

Creb3l1

, was reduced in females, but also in males to a lesser degree, perhaps because of high sensitivity in HM-HCR FISH assays. The expression of sex-shared Esr1-DEGs:

Gria3

,

Pdzrn4

, and

Ctnna2

, were decreased in Esr1KO females and males, and female-specific Esr1-DEG,

Tead1

, was reduced only in females (

Fig. 8f-g

). These results further suggest that

Esr1

orchestrates sex-shared and -specific transcriptional states via an

Esr1

-specific gene regulatory network. These dynamic gene regulatory networks may in turn dictate the sex-shared and -specific adolescent changes that facilitate the maturation of sexually dimorphic neuronal circuits for mating.

Discussion

Adolescence represents a critical period for the maturation of behaviors and brain functions necessary for reproduction and social interactions. Unlike homeostatic processes such as eating and breathing, many behaviors that define adulthood emerge through hormonally driven neural circuit refinements during adolescence

<sup>1,

3

,

4</sup>


. Our study provides a detailed analysis of transcriptional dynamics in MPOA Esr1+ neurons, revealing their essential role in the adolescent maturation of mating behaviors in both sexes. Using a combination of single-cell RNA sequencing, advanced

in situ

hybridization techniques, and functional genomic manipulations, we show that Esr1 is a central regulator of these adolescent transcriptional changes. We show that MPOA Esr1+ neurons are transcriptionally dynamic during adolescence, progressing from an immature state in preadolescence to mature adult-like states in a hormone-dependent manner. This maturation process was arrested in mice lacking

Esr1

in MPOA Vgat+ neurons during adolescence, leading to profound deficits in mating behaviors. Importantly, these findings highlight that transcriptional programs driving sexual maturation are cell-type specific and temporally constrained, dependent on precise hormonal signaling during adolescence. Adolescent transcriptional dynamics in MPOA Esr1+ neurons appear to be continuous, with trajectories shaped by circulating sex hormones. This hormone-dependent maturation was confirmed by pseudotime analysis, which revealed that hormonal deprivation arrests transcriptional progression, while hormone supplementation accelerates it. These findings align with the broader view that the adolescent brain undergoes extensive molecular and circuit-level adaptations to support the emergence of sexually dimorphic behaviors


<sup>3

,

16

,

44</sup>


.

Integration of MPOA sexual dimorphism and adolescence-related transcriptional dynamics

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

Pou6f2

in males,

Creb3l1

in females) further highlight how Esr1 integrates hormonal signals to shape sex-specific neural circuits


<sup>42

,

43</sup>


.

Functional implications for neural circuitry and behavior

The adolescent transcriptional changes observed in MPOA Esr1+ neurons likely translate to functional adaptations in mating-relevant neural circuits. Previous studies have shown that MPOA neurons are critical for a wide range of sexually dimorphic behaviors, including mounting, vocalization, and parental care


<sup>12

,

16

,

20

,

45</sup>


. 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


<sup>18

,

22</sup>


. 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


<sup>29

,

46</sup>


. Future studies combining transcriptional profiling and CRISPR gene editing with

in vivo

imaging and circuit manipulation will be essential to dissect these connections. Our study complements previous work on perinatal transcriptional dynamics in sexually dimorphic brain regions, such as the BNST


²³


. While the perinatal period is critical for establishing primary sexual differentiation, our findings emphasize adolescence as a secondary critical period where sexually dimorphic transcriptional programs are refined. Notably, both periods rely heavily on Esr1-mediated regulation, suggesting a conserved role for estrogen signaling across developmental stages.

Limitations of the study

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

Technical limitations of the study

1. 
   

   
   Although HM-HCR experiments showed the bidirectional control of transcriptional progression during adolescence, it is unclear if the facilitation in male 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 to 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 event, which is 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 robust effect of Esr1-KO in scRNAseq experiment which was further validated with 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, 3-4 mice were pooled for each sample in scRNAseq experiment and in 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 alternations in VgluT2
   
   ⁺
   
   neurons. Although the scale of our study is comparable or exceeds prior scRNAseq studies in MPOA
   
   
   <sup>22
   
   ,
   
   29</sup>
   
   
   , 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 in 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 adult. Future studies injecting the virus at adult may provide additional insights on the similarity and difference between transcriptional changes during puberty and maintained transcriptional states at adult.
   

   
   

Conclusions and future directions

These findings have implications beyond mating behaviors. The role of Esr1+ neurons in the MPOA likely extends to other hormonally sensitive behaviors and physiological processes, such as thermoregulation, metabolic control, and social bonding


<sup>17

,

21

,

47</sup>


. 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


<sup>1

,

2</sup>


. While our study focused on Esr1, it is likely that other transcription factors and hormone receptors synergistically contribute to adolescent transcriptional dynamics. For example, the androgen receptor may play a role in males, while additional regulatory networks may be governed by sex chromosome-linked genes


⁴³


. Advanced genomic approaches, such as perturb-seq, could help unravel the complexity of these gene regulatory networks in a cell-type-specific manner


³¹


. Although we focused on adolescent transcriptional trajectories of Vgat+Esr1+ neurons, this dataset provides a rich resource for studying transcriptional states across other neuronal, glial, and stromal MPOA cell types during adolescent brain development. Understanding steroid-induced changes in transcriptional trajectories during two major developmental periods (perinatal and adolescence periods) has implications for the healthy and maladaptive development of human brains, and provides additional insight into neurobiological mechanisms that underlie sex differences


<sup>48

–

50</sup>


.

Materials and Methods

Mice

Several lines of male and female mice were used for scRNAseq, HM-HCR FISH, immunohistochemistry (IHC), and behavioral experiments, with specific ages and experimental conditions described in their respective sections. All mice used in these experiments were on a C57BL/6 background. Wild-type C57BL/6J mice were originally obtained from Jackson Laboratory (JAX) and bred in-house. Transgenic

Slc32a1


^Flp

::

Esr1


^lox/lox

mice were generated by crossing

Slc32a1


^Flp

mice (Jackson Laboratory stock no. 029591) with

Esr1


^lox/lox

mice (kindly provided by Kenneth Korach; Jackson Laboratory stock no. 032173), and then breeding the resulting progeny with

Slc32a1


^Flp

::

Esr1


^lox/+

mice. Similarly,

Slc17a6


^Flp

::

Esr1


^lox/lox

mice were generated by crossing

Slc17a6


^Flp

mice (Jackson Laboratory stock no. 030212) with

Esr1


^lox/lox

mice and then crossing the progeny with

Slc17a6


^Flp

::

Esr1


^lox/+

mice. For both lines, only mice heterozygous for Flp and homozygous for Esr1

^lox/lox

were used in experiments. For scRNAseq and HCR experiments, mice were group-housed except for the 1–2 days prior to tissue isolation, during which they were singly housed. Behavioral experiment mice were singly housed for the duration of the experiments. Wild-type Swiss Webster (SW) mice, originally obtained from Taconic and bred in-house, were used as lactating foster dams to rear pups following brain surgery, as C57BL/6 dams frequently committed infanticide after such procedures. All mice were maintained with

ad libitum

access to food and water under a reverse 12-hour light-dark cycle. All experiments were conducted in accordance with the National Institutes of Health’s

Guide for the Care and Use of Laboratory Animals

and approved by the Institutional Animal Care and Use Committees (IACUC) of the University of North Carolina and University of Washington.

Stereotaxic viral injection

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

Gonadectomy

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

Generation of AAV-frtFlex-Cre virus

The split-Cre design used in this study was previously described


<sup>51

,

52</sup>


. In this system, an intron is introduced into the middle of the Cre sequence, with the second exon flanked by frt sites. The action of FLP recombinase inverts the second exon, enabling splicing to generate functional Cre. For these experiments, the construct was incorporated into an AAV vector. The frtFlex-Cre-EGFP cassette was subcloned into the pAAV-CAG-WPRE vector between AscI and XhoI restriction sites. AAV1 virus carrying frtFlex-Cre-EGFP was produced as described previously


⁵²


. Briefly, pAAV-CAG-frtFlex-Cre-EGFP-WPRE was transiently transfected into HEK293T/17 cells (ATCC) along with the pDG1 packaging plasmid. Twenty-four hours post-transfection, cells were transferred to serum-free media, and viral capsids were harvested 48 hours later. Purification of AAV1 capsids was performed using gradient centrifugation. Viral titer was estimated at approximately 2 x 10

¹²

particles/mL based on densitometry analysis following gel electrophoresis, using a known standard for reference.

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

To confirm that AAV-frtFlex-Cre does not produce functional Cre in the absence of FLP recombinase, 300 nl of AAV-frtFlex-Cre was unilaterally injected into the MPOA of

Esr1


^lox/lox

mice (n = 3 males and 3 females; data combined across sexes) (

Supplementary Fig. 1a

). To validate the specificity of AAV-frtFlex-Cre, 300 nl of the virus was unilaterally injected into the MPOA of

Slc32a1


^Flp

and

Slc17a6


^Flp

mice. Two weeks post-transduction, brains were collected and analyzed using FISH, as described in the RNAscope section (

Supplementary Fig. 1b

). Six weeks after the injection, brain tissue was collected and analyzed via IHC, as detailed in the Histology section. To assess the efficiency of AAV-frtFlex-Cre, 300 nl of the virus was bilaterally injected into the MPOA of

Slc32a1


^Flp

::

Esr1


^lox/lox

and

Slc17a6


^Flp

::

Esr1


^lox/lox

mice. These mice were first used for behavioral experiments, as described in the Behavioral Experiment section, and efficiency was quantified post-behavioral testing. Quantitative results for each mouse line and sex are presented in

Fig. 1

and

Supplementary Fig. 1

.

Single-cell preparation and cDNA library construction for scRNAseq

Male and female wild-type mice (n = 3–4 animals pooled per group) were prepared at P23 ± 1 (preadolescence, hormonally intact), P35 ± 1 (mid-adolescence, hormonally intact), and P50 ± 2 (hormonally intact early adulthood or gonadectomy groups). For gonadectomy groups, mice were gonadectomized at P22–23. Additionally, male and female

Slc32a1


^Flp

::

Esr1


^lox/lox

mice (n = 3–4 per group), which had received bilateral MPOA injections of 300 nl AAV-frtFlex-Cre at P14–18, were collected at P50 ± 2 (early adulthood). All mice were singly housed 1–2 days before tissue dissection.

Single-cell preparation followed previously published procedures


⁵³


. Immediately after removal from their home cages, mice were deeply anesthetized with an intraperitoneal injection of 0.2 mL sodium pentobarbital (39 mg/mL) and phenytoin sodium (5 mg/mL), followed by transcardial perfusion with ice-cold NMDG-aCSF containing transcription and translation inhibitors to minimize transcriptional events induced by anesthesia, perfusion, or brain extraction. NMDG-aCSF consisted of 96 mM NMDG, 2.5 mM KCl, 1.35 mM NaH2PO4, 30 mM NaHCO3, 20 mM HEPES, 25 mM glucose, 2 mM thiourea, 5 mM Na+ ascorbate, 3 mM Na+ pyruvate, 0.6 mM glutathione-ethyl-ester, 2 mM N-acetyl-cysteine, 0.5 mM CaCl2, and 10 mM MgSO4 (pH 7.35–7.40, 300–305 mOsm, oxygenated with 95% O2 and 5% CO2). Inhibitor cocktails included 500 nM TTX, 10 μM APV, 10 μM DNQX, 5 μM actinomycin, and 37.7 μM anisomycin. Brains were extracted, and coronal sections containing the MPOA were sliced at 300 µm using a vibratome (Leica VT1200). Slices (3–4 per animal) were recovered on ice in a chamber for 30 minutes.

MPOA tissues were dissected under a light microscope (Leica MZFL) using a micro scalpel (Feather) from 3–4 animals per group (10–15 slices pooled per group). The tissue was enzymatically digested with 1 mg/mL pronase (Sigma-Aldrich) for 30 minutes at room temperature, followed by mechanical trituration using fire-polished glass capillaries (tip diameter 200–300 µm). Cell suspensions were filtered twice through 40 µm strainers to remove aggregates, and dead cells were eliminated using a dead cell removal kit (Miltenyi Biotec). After centrifugation, the supernatant was removed, and cells were resuspended (20–30 mL). A fraction (∼5 mL) was mixed with trypan blue for viability assessment and cell concentration measurement using a hemocytometer. Samples with >80% viability were used for cDNA library preparation, and final cell concentrations were adjusted to 800–1,000 cells/µL. Typically, 30,000–80,000 cells were collected from 3–4 mice per group, exceeding the number required for cDNA library construction and reducing biological variability between subjects.

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

Highly Multiplexed HCR FISH

HM-HCR FISH assays were conducted to validate trajectory inference from scRNAseq data (HCR1;

Fig. 4

and

Supplementary Fig. 8–9

) and to cross-validate the Esr1-GRN analysis (HCR2;

Fig. 8

).

For HCR1, we analyzed groups included in the scRNAseq experiments, including preadolescence (P23 ± 1), adolescence (P35 ± 1), intact adult (P50 ± 2), and gonadectomy groups (P50 ± 2, gonadectomy performed at P22–23). Additionally, a hormone-supplemented group received daily injections of sex hormones (females: 5 µg β-Estradiol 3-Benzoate in 50 µl sesame oil; males: 200 µg Testosterone Propionate in 50 µl sesame oil) from P22 to P27, with tissue collected at P28. A control group for hormone supplementation received daily injections of sesame oil (50 µl) during the same period, and brains were harvested at P28. The timing of hormone supplementation aligned with the onset of adolescence in peripheral tissues (P28–33;

Fig. 1c, f

) and reports from the pituitary (P25–30;


⁵⁴


). Adolescence onset was verified in hormone-supplemented animals by assessing male balanopreputial separation (BPS) and female vaginal opening (VO). All hormone-supplemented animals displayed BPS or VO by P28, whereas none of the controls reached adolescence by this age. All mice in the 12 groups were singly housed 1–2 days prior to tissue collection.

For HCR2, AAV-Cre-YFP or AAV-Flp-YFP was unilaterally injected into the MPOA of Esr1

^lox/lox

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

Sequential HCR FISH followed a modified protocol based on


⁴⁰


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

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

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

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

For quantification, HCR FISH images were analyzed in ImageJ. DAPI-defined ROIs were transferred to the 12-plex image to measure transcript numbers for each gene. For HCR1, genes enriched in Vgat+Esr1+ neurons and HA-DEGs from scRNAseq data were analyzed, with validation by SVM (

Supplementary Fig. 6k,m

). For HCR2, Esr1-TFs and Esr1-DEGs regulated by Esr1-TFs or Esr1 were selected from the Esr1-GRN (

Fig. 8

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

RNAscope

AAV-frtFlex-Cre (300 nl) was unilaterally injected into the MPOA of VgatFlp or Vglut2Flp mice (n = 2 per experiment). After two weeks of viral incubation, mice were deeply anesthetized with isoflurane, and brains were rapidly extracted and frozen on dry ice. Coronal sections (20 µm) were prepared using a cryostat (Leica) and stored at -80°C until further use. Probe hybridization and signal detection were performed according to the manufacturer’s instructions (ACDbio). In VgatFlp mice, probes targeted

Slc32a1

and

eGFP

, while in Vglut2Flp mice, probes targeted

Slc17a6

and

eGFP

. Tiled images of the MPOA were acquired using a Zeiss ApoTome2 fluorescence microscope with a 20x objective and Zen software (Zeiss). Image acquisition settings were consistent across all experiments. The resulting CZI files were analyzed using Fiji. The number of cells expressing

Slc32a1

or

Slc17a6

and

eGFP

, as well as double-positive cells expressing both markers, were quantified for each mouse line (

Fig.1

and

Supplementary Fig. 1

).

Behavioral experiments

Male and female

Slc32a1


^Flp


::Esr1


^lox/lox

or

Slc17a6


^Flp


::Esr1


^lox/lox

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

To assess the adolescent maturation of sexual behaviors, male subjects underwent RI assays from P34 to P54 every other day. In this assay, a hormonally primed adult female mouse was introduced into the male subject’s home cage for a 15-minute interaction (Cre groups:

Slc32a1


^Flp


::Esr1


^lox/lox

, n = 12;

Slc17a6


^Flp


::Esr1


^lox/lox

, n = 9; control groups:

Slc32a1


^Flp


::Esr1


^lox/lox

, n = 13;

Slc17a6


^Flp


::Esr1


^lox/lox

, n = 9). To prevent ejaculation, male subjects were separated from the female intruder as soon as thrusting began. The interactions were video-recorded, and the number of mounting and thrusting behaviors was manually counted.

Female sexual receptivity was assessed in RI assays conducted from P35 to P54 for

Slc32a1


^Flp


::Esr1


^lox/lox

mice (Cre group: n = 16; control group: n = 10) and from P40 to P60 for

Slc17a6


^Flp


::Esr1


^lox/lox

mice (Cre group: n = 9; control group: n = 11). Assay timing was guided by pilot experiments and conducted only when vaginal smears indicated the female was in proestrus or estrus. Female subjects were introduced into the home cage of a sexually experienced adult male mouse and allowed to interact freely for 10 minutes. To avoid ejaculation, females were separated from the male after intromission began or approximately 5 seconds after an unsuccessful mounting attempt. These interactions were recorded using an IR camera controlled by Ethovision (Noldus), and behaviors such as being mounted, intromitted, escaping, or displaying submissive postures were manually counted. Receptivity was calculated as the proportion of intromission bouts out of total mounting attempts. Stress levels were minimized and equalized across groups by selecting highly sexually experienced but minimally aggressive male residents. In rare instances of aggression, male residents were promptly removed. Unlike mice with brain-wide

Esr1

deletions


²⁷


, the male subjects in this study, with

Esr1

knockouts limited to MPOA Vgat+ neurons during adolescence, showed no aggressive behaviors, consistent with previous findings


¹⁵


.

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

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

Histology

Histological experiments were conducted to: (1) examine viral infection in the MPOA of mice used in behavioral experiments (

Slc32a1


^Flp


::Esr1


^lox/lox

or

Slc17a6


^Flp


::Esr1


^lox/lox

), (2) validate the specificity of AAV-frtFlex-Cre in

Esr1


^lox/lox

mice, and (3) confirm the specificity of reporter gene expression in

Slc32a1


^Flp


::Esr1-Cre::RC-FLTG

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

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

Analysis of scRNAseq, MERFISH, and HM-HCR FISH data utilized several R and Python packages, which were adapted and modified as needed. Previously published MERFISH data, deposited in DRYAD, was acquired from

. The data analysis workflow encompassed clustering, differential gene expression analysis, integration of differential conditions and modalities, lineage inference, gene set activity measurements, and gene regulatory network (GRN) deconstruction.

Clustering and differential gene expression analyses were performed using the Seurat V3 package


³²


, while LISI


⁵⁵


was applied to evaluate integration performance across conditions. Lineage inference, representing transcriptional connectivity of neuronal clusters, was conducted with PAGA (partition-based graph abstraction; Wolf et al., 2019). Single-cell gene set activity was quantified using AUCell


³¹


. GRNs were inferred and deconstructed using the SCENIC package, which integrates GENIE3


<sup>31

,

56</sup>


and RcisTarget. These analyses depicted cis-regulatory relationships between transcription factors and differentially expressed genes. Gene co-expression networks were identified using scWGCNA


⁵⁷


. To examine transcriptional state progressions by age and hormonal state, Monocle V3


³⁴


and MELD


³⁸


were applied. Ontology analysis was conducted using the online Enrichr platform


<sup>58

,

59</sup>


.

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

Data preprocessing and doublet removal for scRNAseq data

Preprocessing followed the procedures described in our previously published study


⁵³


. Low-expression genes, defined as those detected in fewer than three cells, and low-quality cells (total UMI < 700, total UMI > 15,000, or >20% mitochondrial gene expression) were excluded from downstream analysis. Suspected doublets were computationally removed using the DoubletDecon package (Version 1.1.5;


⁶⁰


) with default settings, applying a doublet rate of 5.6 ± 1.0%.

Integrative clustering and differential gene expression analysis

To minimize the effects of batch variability and differences due to sex, age, and hormonal state while preserving the global similarity of transcriptional states within cell types, we applied the Seurat V3 integrative approach


³²


. This method combines canonical correlation analysis (CCA;


⁶¹


) and mutual nearest neighbor analysis


⁶²


. After preprocessing the data, which included 58,921 cells in total, gene counts were normalized using the NormalizeData function to scale counts by total UMI with a constant scale factor (10,000), followed by natural-log transformation (log1p). The FindVariableFeatures function was used to select 2,000 highly variable genes from each sample based on variance stabilizing transformation. Integration was performed pairwise using the FindIntegrationAnchors function (CCA1-40) to identify anchors and assign scores, followed by the IntegrateData function to compute an integrated gene expression matrix by iteratively constructing and subtracting transformation matrices from the original data.

The integrated expression data was used for clustering. Expression matrices were scaled, centered, and reduced using principal component analysis (PCA). A nearest-neighbor graph was constructed in PCA space (FindNeighbors; default settings), followed by Louvain clustering at a resolution of 0.8 (FindClusters). Uniform Manifold Approximation and Projection (UMAP) was then generated for visualization (RunUMAP). Initial clustering identified 32 clusters, 13 of which were neuronal, expressing canonical markers such as

Thy1

or

Stmn2

. Cell types were determined based on marker expression (Neuron:

Stmn2, Thy1;

Astrocyte:

Ntsr2;

OPC:

Gpr17;

Oligodendrocyte:

Mog;

Microglia:

C1qc;

Mural cell:

Tagln;

Endothelial cell:

Flt1;

Intermediate cell:

B2m;

Ependymal cell:

Foxj1

) as reported in previous studies


<sup>53

,

63

,

64</sup>


. Neuronal and non-neuronal cell types were combined for visualization (

Supplementary Fig. 2a-d

).

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

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

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

To assess integration performance, Local Inverse Simpson’s Index (LISI;


⁵⁵


) was computed in UMAP space, with LISI scores reflecting the effective number of distinct groups in the local neighborhood of each cell (

Supplementary Fig. 2i

). Clustering stability was further evaluated by examining lineage relationships between clusters at varying resolutions (0.2–2.0) using the clustree package


⁶⁵


. Clustering tree analysis showed the emergence of heterogeneous

Esr1

clusters at resolution 0.8 or higher (

Supplementary Fig. 2k

). Specific marker genes in neuronal clusters (

Supplementary Fig. 3a

) and Vgat+Esr1+ subclusters (

Supplementary Fig. 3h

) further validated clustering.

To explore transcriptional connectivity between neuronal clusters, PAGA


³⁰


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

Fig. 2c

).

Integrative analysis of MERFISH and scRNAseq data

Seurat V3 was used to jointly analyze publicly available POA MERFISH data


²²


and our scRNAseq data. To establish correspondence between MERFISH and scRNAseq clusters related to reproductive behaviors in adult mice, cells from behaviorally naïve subjects (lacking

Fos

expression data), lactating females, or individuals exhibiting aggression toward pups were excluded. Cells categorized as “Inhibitory” or “Excitatory” in the MERFISH metadata were extracted and normalized using the NormalizeData function. Highly variable genes were identified with FindVariableFeatures, and the top 60 variable genes (excluding

Fos

) were used for clustering, following the procedures described earlier.

Seurat clustering identified 19 neuronal clusters within the MERFISH dataset, including 10 GABAergic clusters enriched with

Gad1

and 6 excitatory clusters expressing

Slc17a6

. Because

Fos

data from behaviorally naïve animals was unavailable in the MERFISH dataset,

Fos

-enriched clusters were defined using the following approach: For each behavioral category and sex, a threshold for

Fos

enrichment was set at the one-sided 95th percentile expression level. The proportion of cells above this threshold was calculated for each cell type, and Fisher’s exact test was performed to identify

Fos

-enriched clusters (

p

< 0.05 after multiple comparison correction). As

Fos

-enriched clusters were only observed in a subset of GABAergic clusters, subsequent analyses focused exclusively on the correspondence of GABAergic clusters between modalities.

To map MERFISH clusters onto scRNAseq data, reference (MERFISH) and query (scRNAseq, P50) datasets were aligned for each sex using the FindTransferAnchors (CCA1-30) and TransferData functions. This generated a weights matrix (

W

), which was used to transfer cell-type labels through the equation P = LW

^T

, where L is a binary classification matrix and P represents label predictions. The imputed MERFISH cluster labels on scRNAseq cells enabled the inference of scRNAseq clusters associated with social behaviors. Genes selectively enriched in scRNAseq clusters linked to sexual behaviors were identified using the FindMarkers function.

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

Gene set activity analysis of scRNAseq data

To assess the aggregate expression of specific gene sets in single cells, we performed AUCell analysis (Aibar et al., 2017). AUCell calculates the Area Under the Curve (AUC) to determine whether a given gene set is enriched among the expressed genes in each cell. Using this approach, we computed gene set activity for hormonally associated DEGs (HA-DEGs;

Fig. 2g

),

Ar

-regulon in male-dimorphic genes (

Fig. 5h–j

), and adolescence-related versus unrelated dimorphic genes (

Fig. 6g–h

).

Trajectory analysis of scRNAseq and HM-HCR FISH data

Monocle 3 was used to quantify adolescent transcriptional progression in scRNAseq and HM-HCR FISH datasets


³⁴


. For scRNAseq data, Seurat clustering results were used to separately analyze Vgat+Esr1+, Vgat+Esr1-Ar+, Vgat+hormoneR

^Low

, and Vglut2+Esr1+ clusters in each sex. The preprocessing step involved normalizing and scaling the data using the preprocess_cds function, followed by PCA dimensional reduction based on hormonally associated DEGs (HA-DEGs; using the top 10–15 principal components). The reduce_dimension function was applied to generate UMAP embeddings, and a principal graph in UMAP space was constructed using the learn_graph function with default settings. The root node of the trajectory was assigned to cells from P23 or GDX groups, and pseudotime, representing transcriptional progression from the root state, was computed as the geodesic distance of each cell from the root node in the UMAP.

To test whether hormonally dependent transcriptional trajectories observed in Vgat+Esr1+ cells were also present in other cell types, similar Monocle trajectory analyses were conducted for Vgat+Esr1-Ar+, Vgat+hormoneR

^Low

, and Vglut2+Esr1+ cells. Additionally, branches of transcriptional trajectories were analyzed due to the presence of two major branches originating from one or two subpopulations of Vgat+Esr1+ cells in both sexes (

Supplementary Fig. 5

,

12

). Branch-specific DEG analysis identified genes enriched in each branch and those shared across branches (

Supplementary Fig. 5

,

12

).

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

Supplementary Fig. 6k,m

).

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

To jointly compute transcriptional trajectories from scRNAseq and HM-HCR data, scRNAseq data was imputed using HM-HCR data via Seurat V3 (

Supplementary Fig. 9

). In Vgat+Esr1+ cells,

Slc32a1

and one HM-HCR gene were excluded from the scRNAseq dataset. The FindTransferAnchors (CCA) and TransferData functions were used to identify anchors between the HM-HCR reference and the scRNAseq query, generating a weights matrix (

W

). Gene expression features were transferred as P = FW

^T

, where F is the gene expression matrix and P is the predicted expression matrix. This process was iterated for all HM-HCR genes to generate predicted scRNAseq data (11–12 genes). The predicted scRNAseq data was validated against the real dataset using Pearson correlation coefficient analysis. Transcriptional trajectories were then learned from the predicted data using the same Monocle pipeline as described above.

MELD analysis

To computationally validate the monocle trajectory analysis of scRNAseq data, we used the MELD Python package, which learns the transition of transcriptional states as the relative likelihood of conditions for each cell in the PHATE space (Potential of Heat-diffusion for Affinity-based Trajectory Embedding)


<sup>38

,

39</sup>


. PHATE preserves local and distal relationships of transcriptional states, outperforming other embedding techniques in maintaining the structure of single-cell data. The input data consisted of log-normalized expression values and clustering metadata for Vgat+Esr1+ cells, generated from Seurat V3 analysis. The analysis proceeded in two main steps. First, data was embedded into low-dimensional PHATE space. Principal components (PCs) were computed using scprep.reduce.pca (number of components = 8–15), and PHATE embeddings were generated with phate.PHATE and phate_op.fit_transform using default settings. Second, the kernel density for each condition was estimated to quantify the likelihood of cells belonging to specific states. Conditions were categorized as mature (P50, P35) or immature (P23, GDX) (

Supplementary Fig. 6h-i, j, l

). Kernel density estimation was performed using meld.MELD (beta = 67, knn = 7) and meld_op.fit_transform with default settings. Relative likelihoods for each condition (the probability of observing a cell in each condition) were computed using the replicate_normalize_densities function, which applies L1 normalization across samples. These MELD analyses were used to measure hormonally associated adolescent trajectories, providing computational cross-validation for the monocle trajectory analysis (

Supplementary Fig. 6h-i, j, l

).

DEG and trajectory analysis of scRNAseq in Esr1KO experiments

Data from

Esr1

KO groups were independently processed for quality control, including the removal of low-quality cells and doublets


⁶⁰


. Seurat V3 was used to perform clustering and identify neuronal clusters (detailed in Integrative Clustering and Differential Gene Expression Analysis). The dataset included 6,727 cells from

Esr1

KOF (female) and 6,464 cells from

Esr1

KOM (male). Neuronal cells were re-clustered to identify Vgat+ clusters. In females, six Vgat+ clusters (eVgat1–6), seven Vglut2+ clusters (eVglut1–7), and one mixed cluster (eMix1) were identified, while in males, 12 Vgat+ clusters (eVgat1–12), six Vglut2+ clusters (eVglut1–6), and two mixed clusters (eMix1–2) were identified (

Supplementary Fig. 11

;

Esr1

KOF: 2,182 neuronal cells;

Esr1

KOM: 4,761 neuronal cells).

To identify Vgat+ clusters in the

Esr1

KO 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

Esr1

KO dataset corresponded to Vgat+Esr1+ clusters (e.g., Vgat2, 4, 16) or Vgat+hormoneR

^Low

clusters (e.g., Vgat14, 17, 20) in the intact group (

Supplementary Fig. 11

). To address potential experimental variations introduced by viral manipulations, this correspondence was cross-validated using anchor-based integrative analysis. Using the FindTransferAnchors (CCA1-30) and TransferData functions, anchors were identified between the intact group (reference) and

Esr1

KO group (query), constructing a weights matrix (

W

). Labels (cell-type classifications) were transferred using the equation P = LW

^T

, where L is a binary classification matrix and P is the label prediction. Label imputation confirmed that eVgat1 and eVgat3 (female), and eVgat3 and eVgat4 (male), corresponded to the Vgat+Esr1+ population of intact groups. Similarly, eVgat5 (female) and eVgat9 (male) were closely related to the Vgat+hormoneR

^Low

clusters of P50 groups (

Supplementary Fig. 11

). Differentially expressed genes (DEGs) specific to

Esr1

were computed using Seurat’s FindMarkers function by comparing the gene expression of Vgat+Esr1+ clusters between intact P50 and

Esr1

KO groups. For trajectory analysis, Vgat+Esr1+ cells from all groups were combined and analyzed using Monocle 3. Principal components were computed based on the same HA-DEGs as described previously, and these were used to construct UMAP embeddings. In UMAP space, a principal graph was learned, root nodes were defined, pseudotime was computed, and branch-specific analysis was performed (detailed in Trajectory Analysis of scRNAseq and HM-HCR FISH Data).

Deconstruction of gene regulatory network

To infer the cis-regulation of genes by transcription factors (TFs), we used the SCENIC computational framework


<sup>31

,

53

,

66</sup>


. SCENIC was combined with cell-type-specific gene manipulation to causally infer the gene regulatory network (GRN). First, co-expression scores between TFs and target genes were computed using GENIE3 (runGenie3) based on scRNAseq expression data from Vgat+ cells at P50. TFs were ranked by the sum of their weights for hormonally associated DEGs (HA-DEGs; logFC > 0.25). Among all 952 TFs (female) and 915 TFs (male) in the dataset,

Esr1

achieved the highest sum score in both sexes. These results prompted further scRNAseq analysis of subjects with

Esr1

deleted in Vgat+ cells in the MPOA (detailed in Trajectory Analysis of scRNAseq in

Esr1

KO Experiments). After identifying Vgat+Esr1+ clusters, genes reduced by

Esr1

deletion (

Esr1

-DEGs) were identified using the FindMarkers function in Seurat V3. To determine which

Esr1

-DEGs were cis-regulated by

Esr1

, we applied the RcisTarget package. The RcisTarget database (mm9-tss-centered-10kb-7species.mc9nr.feather) was used to score motifs within the transcription start site (TSS) ±10 kb region and annotate associated TFs. Overrepresentation of each motif was assessed using the calcAUC function, and significantly enriched motifs were identified using the addMotifAnnotation function (normalized enrichment score threshold: NES ≥ 2). This analysis identified 123 genes in females and 98 genes in males that were significantly enriched with

Esr1

motifs. Of these, 9 (female) and 6 (male) genes were TFs (

Esr1

-TFs). To investigate whether

Esr1

-TFs could cis-regulate

Esr1

-DEGs, we repeated RcisTarget analysis for pairs of

Esr1

-TFs and

Esr1

-DEGs. These analyses allowed us to infer genes cis-regulated by

Esr1

or each

Esr1

-TF, collectively defining regulons for each TF. Using these regulons, we constructed the

Esr1

-GRN with the igraph package. In the

Esr1

-GRN, regulons were categorized and highlighted based on gene ontology (detailed in Ontology Analysis,

Fig. 8d-e

, and

Supplementary Fig. 13

). This approach provided a causal framework to understand the regulatory influence of

Esr1

and its downstream TFs on hormonally associated transcriptional networks.

Single-cell consensus weighted gene co-expression network analysis

Single-cell consensus-weighted gene co-expression network analysis (scWGCNA) was employed to identify sex-specific co-expression networks in Vgat+Esr1+ cells


⁵⁷


. The analysis aimed to detect modules (groups of coexpressed genes) highly expressed in Vgat+Esr1+ clusters. To begin, metacells were constructed using the MetacellsByGroups function. Co-expression networks were then generated with the ConstructNetwork function, using the optimal soft-power threshold identified through TestSoftPowers. Module eigengenes (MEs), which represent the principal component of gene expression for each module, were computed using the ModuleEigengenes function. Module connectivity, calculated with the ModuleConnectivity function, was used to identify hub genes—those highly connected within each module. The top 25 hub genes and their connectivity were visualized using the ModuleNetworkPlot function. This analysis revealed nine modules in females and six modules in males that were highly expressed in Vgat+Esr1+ cells (Vgat2, 4, 16). Functional enrichment for each module was assessed using Enrichr via the RunEnrichr function


<sup>58

,

59</sup>


. To directly compare modules between sexes, Jaccard similarity was computed across all module pairs to quantify sex-specific similarities and differences in gene co-expression networks.

Ontology analysis

Ontology analysis was performed using Enrichr to examine the enrichment of functionally related genes within

Esr1

-DEGs (

Supplementary Fig. 13

). Gene ontology (GO) terms were identified from the GO Biological Process, GO Molecular Function, and Cellular Component databases. Detailed GO terms and their associated genes are reported in

Supplementary Fig. 13

. To simplify interpretation, GO terms were grouped into four categories: Gene Expression, Synapse/Excitability, Morphology, and Energy. Genes associated with multiple categories were manually assigned to the most relevant category to streamline classification and analysis.

Decoding conditions from scRNAseq or HM-HCR FISH

To evaluate whether the selected features from trajectory analyses of scRNAseq and HM-HCR FISH experiments were sufficient to decode cellular conditions (e.g., age, sex, hormonal states, transcriptional maturity), we performed support vector machine (SVM) classification using the scikit-learn library. GridSearch with 10-fold cross-validation was applied to optimize the classification as described in previous studies


<sup>67

,

68</sup>


. 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 differentially expressed genes (DEGs) used in constructing trajectories for scRNAseq and HM-HCR FISH. Decoding accuracies were compared against baseline accuracies calculated using randomized features (see Detailed Statistical Procedures).

Detailed statistical procedures

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

Fig. 1b

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

Fig. 1c

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

Fig. 1d left:

unpaired t-test. t(24)=0.1302, p=0.8975.

Fig. 1d right:

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

Fig. 1e left:

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

Fig. 1e right:

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

Fig. 1f

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

Fig. 1g

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

Fig. 1h left:

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

Fig. 1h right:

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

Fig. 1j

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

Fig. 1k

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

Fig. 1l left:

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

Fig. 1l right:

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

Fig. 1m left:

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

Fig. 1m right:

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

Fig. 1n

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

Fig. 1o

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

Fig. 1p left:

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

Fig. 1p right:

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

Fig. 2e

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

Fig. 2f

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

Fig. 2g

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

Vgat2: Kruskal-Wallis chi-squared= 235.66, df=3, p<0.0001. Vgat4: Kruskal-Wallis chi-squared= 273.36, df=3, p<0.0001. Vgat16: Kruskal-Wallis chi-squared= 95.768, df=3, p<0.0001. Vgat6: Kruskal-Wallis chi-squared= 96.91, df=3, p<0.0001.

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

Vgat1: Kruskal-Wallis chi-squared= 124.49, df=3, p<0.0001. Vgat10: Kruskal-Wallis chi-squared= 11.429, df=3, p= 0.009619. Vgat18: Kruskal-Wallis chi-squared= 8.5876, df=3, p= 0.03531. Vgat7: Kruskal-Wallis chi-squared= 2.5441, df=3, p= 0.4674.

Vgat3: Kruskal-Wallis chi-squared= 5.9714, df=3, p= 0.113. Vgat9: Kruskal-Wallis chi-squared= 6.7542, df=3, p= 0.08016. Vgat17: Kruskal-Wallis chi-squared= 7.4935, df=3, p= 0.05773. Vgat15: Kruskal-Wallis chi-squared= 3.2132, df=3, p= 0.3599. Vgat12: Kruskal-Wallis chi-squared= 5.316, df=3, p= 0.1501.

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

male:

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.

Vgat2: Kruskal-Wallis chi-squared= 309.73, df=3, p<0.0001. Vgat4: Kruskal-Wallis chi-squared= 235.14, df=3, p<0.0001. Vgat16: Kruskal-Wallis chi-squared= 147.85, df=3, p<0.0001. Vgat6: Kruskal-Wallis chi-squared= 99.699, df=3, p<0.0001. Vgat13: Kruskal-Wallis chi-squared= 69.569, df=3, p<0.0001. Vgat8: Kruskal-Wallis chi-squared= 45.965, df=3, p<0.0001. Vgat1: Kruskal-Wallis chi-squared= 120.17, df=3, p<0.0001. Vgat10: Kruskal-Wallis chi-squared= 78.183, df=3, p<0.0001. Vgat18: Kruskal-Wallis chi-squared= 54.178, df=3, p<0.0001. Vgat7: Kruskal-Wallis chi-squared= 34.197, df=3, p<0.0001. Vgat3: Kruskal-Wallis chi-squared= 100.71, df=3, p<0.0001. Vgat9: Kruskal-Wallis chi-squared= 84.388, df=3, p<0.0001. Vgat17: Kruskal-Wallis chi-squared= 24.424, df=3, p<0.0001. Vgat15: Kruskal-Wallis chi-squared= 88.496, df=3, p<0.0001. Vgat12: Kruskal-Wallis chi-squared= 116.7, df=3, p<0.0001. Vgat19: Kruskal-Wallis chi-squared= 55.483, df=3, p<0.0001. Vgat5: Kruskal-Wallis chi-squared= 163.73, df=3, p<0.0001. Vgat14: Kruskal-Wallis chi-squared= 17.716, df=3, p=0.00050. Vgat11: Kruskal-Wallis chi-squared= 61.638, df=3, p<0.0001. Vgat20: Kruskal-Wallis chi-squared= 8.5476, df=3, p = 0.03595.

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

Fig. 3c

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

Fig. 3d

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

Fig. 3e

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

Fig. 3f

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

Fig. 3h

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

⁺

, p< 0.0001. real Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. real hormoneR

^Low

vs shuffle hormoneR

^Low

, p< 0.0001.

Fig. 3j

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

⁺

, p< 0.0001. real Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. real hormoneR

^Low

vs shuffle hormoneR

^Low

, p< 0.0001.

Fig. 3l

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

Fig. 3n

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

Fig. 5b

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

Fig. 5c

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

Fig. 5i

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

Fig. 5j

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

Fig. 6a

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

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

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

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

Fig. 6b

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

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

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

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

Fig. 6c

: One-way ANOVA followed by multiple comparisons.

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

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

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

Fig. 6d

: One-way ANOVA followed by multiple comparisons.

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

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

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

Fig. 6f

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

Fig. 6g

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

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

Fig. 6h

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

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

Fig. 7e

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

Fig. 7f

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

Fig. 7g

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

Fig. 7h

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

Fig. 7j

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

Fig. 7l

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

Fig. 8c

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

Fig. 8f

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

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

Fig. 8g

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

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

Supplementary Fig. 1a

: paired t-test. t(11)=-0.523 , p=0.612.

Supplementary Fig. 1c

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

Supplementary Fig. 1d left

: unpaired t-test. t(24)= -2.966, p=0.0067

Supplementary Fig. 1d right

: unpaired t-test. t(15.466)= -3.7502, p= 0.0018

Supplementary Fig. 1e top

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

Supplementary Fig. 1e bottom

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

Supplementary Fig. 1f top

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

Supplementary Fig. 1f bottom

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

Supplementary Fig. 1g

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

Supplementary Fig. 1h

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

Supplementary Fig. 2i

: unpaired Wilcoxon test. W=396364535, p<0.0001.

Supplementary Fig. 3d left

: one-way repeated measures ANOVA followed by multiple comparisons.

Cluster1: ANOVA, F (135, 405) = 2.9438259, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(135)= 0.073054964, p= 0.9419. P23 vs P35, t(135)= 6.6450382, p <0.0001. P23 vs P50, t(135)= 12.347378, p <0.0001. P35 vs OVX, t(135)= 4.0952302, p< 0.0001. P50 vs OVX, t(135)= 13.269141, p <0.0001. P35 vs P50, t(135)= 11.17 6780, p <0.0001.

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

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

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

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

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

Supplementary Fig. 3d right

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

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

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

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

Supplementary Fig. 3e left

: one-way repeated measures ANOVA followed by multiple comparisons.

Cluster1: ANOVA, F (106, 318) = 1.784, p<0.0001. Holm-Sidak multi-comparison test: P23 vs OVX, t(106)= 17.25, p<0.0001. P23 vs P35, t(106)= 5.938, p<0.0001. P23 vs P50, t(106)= 3.823, p =0.0004. P35 vs OVX, t(106)= 17.54, p< 0.0001. P50 vs OVX, t(106)= 9.249, p <0.0001. P35 vs P50, t(106)= 0.9242, p =0.3575.

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

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

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

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

Supplementary Fig. 3e right:

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

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

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

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

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

Supplementary Fig. 3f:

Linear regression analysis. Adjusted R-squared=0.43, p=0.013.

Supplementary Fig. 3g:

Linear regression analysis. Adjusted R-squared for each hormone receptor gene was reported.

Supplementary Fig. 4a, b:

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

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

Supplementary Fig. 5c left branch1 genes:

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

Supplementary Fig. 5c left branch2 genes:

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

Supplementary Fig. 5c left branch common genes:

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

Supplementary Fig. 5c right branch1 genes:

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

Supplementary Fig. 5c right branch2 genes:

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

Supplementary Fig. 5c right branch common genes:

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

Supplementary Fig. 5d left branch1 genes:

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

Supplementary Fig. 5d left branch2 genes:

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

Supplementary Fig. 5d left branch common genes:

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

Supplementary Fig. 5d right branch1 genes:

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

Supplementary Fig. 5d right branch2 genes:

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

Supplementary Fig. 5d right branch common genes:

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

Supplementary Fig. 5e branch1:

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

Supplementary Fig. 5e branch2:

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

Supplementary Fig. 5f branch1:

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

Supplementary Fig. 5f branch2:

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

Supplementary Fig. 5g branch1:

unpaired t-test. T(152.39)= 216.35, p<0.0001.

Supplementary Fig. 5g branch2:

unpaired t-test. T(147.5)= 153.52, p<0.0001.

Supplementary Fig. 5h branch1:

unpaired t-test. T(130.85)= 244.88, p<0.0001.

Supplementary Fig. 5h branch2:

unpaired t-test. T(124.65)= 204.9, p<0.0001.

Supplementary Fig. 5i branch1 AUCell score:

unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5i branch2 AUCell score:

unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5i branch1 SVM:

unpaired t-test. T(134.62)= 142.72, p<0.0001.

Supplementary Fig. 5i branch2 SVM:

unpaired t-test. T(176.75)= 66.219, p<0.0001.

Supplementary Fig. 5j branch1 AUCell score:

unpaired Wilcoxon test. p<0.0001.

Supplementary Fig. 5j branch2 AUCell score:

unpaired Wilcoxon test. P=0.0050.

Supplementary Fig. 5j branch1 SVM:

unpaired t-test. T(164.27)= 83.34, p<0.0001.

Supplementary Fig. 5j branch2 SVM:

unpaired t-test. T(190.61)= 71.154, p<0.0001.

Supplementary Fig. 6d:

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

Supplementary Fig. 6f:

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

Supplementary Fig. 6j:

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

Supplementary Fig. 6k:

unpaired t-test. T(198)= 61.044, p<0.0001.

Supplementary Fig. 6l:

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

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

Supplementary Fig. 6m:

unpaired t-test. T(198)= 99.879, p<0.0001.

Supplementary Fig. 6o:

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

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

Supplementary Fig. 6q:

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

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

Supplementary Fig. 7d:

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

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

Supplementary Fig. 7e:

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

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

Supplementary Fig. 9a:

Kruskal-Wallis rank sum test followed by multiple comparisons.

Pdzrn4

: Kruskal-Wallis chi-squared=1998, df=5, p<0.001.

Creb3l1

: Kruskal-Wallis chi-squared=507.5, df=5, p<0.001.

Sccpdh

: Kruskal-Wallis chi-squared=2517.7, df=5, p<0.001.

Esr1

: Kruskal-Wallis chi-squared=166.37, df=5, p<0.001.

Socs2

: Kruskal-Wallis chi-squared=1609.2, df=5, p<0.001.

Pgr

: Kruskal-Wallis chi-squared=1928.7, df=5, p<0.001.

Tmem35a

: Kruskal-Wallis chi-squared=1450.7, df=5, p<0.001.

Nts

: Kruskal-Wallis chi-squared=631.36, df=5, p<0.001.

Slc32a1

: Kruskal-Wallis chi-squared=133.17, df=5, p<0.001.

Alcam

: Kruskal-Wallis chi-squared=657.22, df=5, p<0.001.

Ar

: Kruskal-Wallis chi-squared=546.52, df=5, p<0.001.

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

Supplementary Fig. 9b:

Kruskal-Wallis rank sum test followed by multiple comparisons.

Sytl4

: Kruskal-Wallis chi-squared=1532.5, df=5, p<0.001.

Acvr1c

: Kruskal-Wallis chi-squared=936.37, df=5, p<0.001.

Greb1

: Kruskal-Wallis chi-squared=512.3, df=5, p<0.001.

Nrip

: Kruskal-Wallis chi-squared=759.41, df=5, p<0.001.

Lamp5

: Kruskal-Wallis chi-squared=1052.2, df=5, p<0.001.

Apoc3

: Kruskal-Wallis chi-squared=938.64, df=5, p<0.001.

Pgr

: Kruskal-Wallis chi-squared=606.25, df=5, p<0.001.

Esr1

: Kruskal-Wallis chi-squared=87.281, df=5, p<0.001.

Slc32a1

: Kruskal-Wallis chi-squared=192.37, df=5, p<0.001.

Npy2r

: Kruskal-Wallis chi-squared=528.2, df=5, p<0.001.

Pgr151

: Kruskal-Wallis chi-squared=1920.3, df=5, p<0.001.

Ar

: Kruskal-Wallis chi-squared=235.83, df=5, p<0.001.

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

Supplementary Fig. 9d:

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

⁺

, p< 0.0001. real Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. real hormoneR

^Low

vs shuffle hormoneR

^Low

, p< 0.0001.

Supplementary Fig. 9g:

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

⁺

, p< 0.0001. real Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. real Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs real hormoneR

^Low

, p< 0.0001. shuffle Vgat+Esr1+vs shuffle hormoneR

^Low

, p< 0.0001. real hormoneR

^Low

vs shuffle hormoneR

^Low

, p< 0.0001.

Supplementary Fig. 9n:

unpaired t-test. T(198)=91.01, p<0.0001.

Supplementary Fig. 9q:

unpaired t-test. T(198)=157.2, p<0.0001.

Supplementary Fig. 9r:

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

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

Supplementary Fig. 9s:

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

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

Supplementary Fig. 9u:

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

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

Supplementary Figure 10d:

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

Supplementary Fig.12c left branch1 genes:

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

Supplementary Fig.12c left branch2 genes:

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

Supplementary Fig.12c left branch common genes:

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

Supplementary Fig.12c right branch1 genes:

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

Supplementary Fig.12c right branch2 genes:

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

Supplementary Fig.12c right branch common genes:

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

Supplementary Fig.12d left branch1 genes:

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

Supplementary Fig.12d left branch2 genes:

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

Supplementary Fig.12d left branch common genes:

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

Supplementary Fig.12d right branch1 genes:

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

Supplementary Fig.12d right branch2 genes:

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

Supplementary Fig.12d right branch common genes:

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

Supplementary Fig.12e branch1:

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

Supplementary Fig.12e branch2:

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

Supplementary Fig.12f branch1:

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

Supplementary Fig.12f branch2:

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

Data availability

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

).

Acknowledgements

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

Additional information

Funding

National Institute of Health grant DA038168 (GDS) National Institute of Health grant DA032750 (GDS) National Institute of Health grant DA054317 (GDS) The Foundation of Hope (GDS, JM)

Brain and Behavior Research Foundation (KH and MAR) National Institute of Health grant DK121883 (MAR) National Institute of Health grant P30DA048736

Author contributions

Conceptualization: GDS, KH, YH, KI, BB, JM, DR Methodology: KH, YH, KI, BB, MLB, RDP, LSZ

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

Data analysis: KH, YH, KI, BB, GDS, JCS, JRS Manuscript writing: GDS, KH, YH, BB, KI Manuscript editing: all authors

Supervision: GDS

Abbreviations

medial preoptic area (MPOA)

single-cell RNA sequencing (scRNAseq) estrogen receptor 1 (Esr1)

GABAergic neurons (Vgat+) glutamatergic neurons (Vglut2+) hybridized chain reaction fluorescent

in situ

hybridization (HCR-FISH) gonadectomized (GDX) partition-based graph abstraction (PAGA) hormone-associated differentially expressed genes (HA-DEGs) multiplexed error-robust fluorescence

in situ

hybridization (MERFISH) differential gene expression (DE) differentially expressed genes (DEGs) support vector machine (SVM) Manifold Enhancement Latent Dimension (MELD) Potential of Heat-diffusion for Affinity-based Trajectory Embedding (PHATE) androgen receptor (AR) single-cell regulatory network inference (SCENIC)

References

• 
  
  (1)
  
  
  
  
  1. 
     Blakemore S.-J.
     
  
  
  2. 
     Burnett S.
     
  
  
  3. 
     Dahl R. E
     
  
  
  
  
  2010
  
  
  
  The Role of Puberty in the Developing Adolescent Brain
  
  
  
  
  Human Brain Mapping
  
  
  31
  
  :926–933
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (2)
  
  
  
  
  1. 
     McCarthy M. M.
     
  
  
  2. 
     Arnold A. P
     
  
  
  
  
  2011
  
  
  
  Reframing Sexual Differentiation of the Brain
  
  
  
  
  Nat Neurosci
  
  
  14
  
  :677–683
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (3)
  
  
  
  
  1. 
     Romeo R. D.
     
  
  
  2. 
     Richardson H. N.
     
  
  
  3. 
     Sisk C. L
     
  
  
  
  
  2002
  
  
  
  Puberty and the Maturation of the Male Brain and Sexual Behavior: Recasting a Behavioral Potential
  
  
  
  
  Neuroscience & Biobehavioral Reviews
  
  
  26
  
  :381–391
  
  
  
  https://doi.org/10.1016/S0149-7634(02)00009-X
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (4)
  
  
  
  
  1. 
     Schulz K. M.
     
  
  
  2. 
     Sisk C. L
     
  
  
  
  
  2006
  
  
  
  Pubertal Hormones, the Adolescent Brain, and the Maturation of Social Behaviors: Lessons from the Syrian Hamster
  
  
  
  
  Molecular and Cellular Endocrinology
  
  
  254–255
  
  :120–126
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (5)
  
  
  
  
  1. 
     Juntti S. A.
     
  
  
  2. 
     Tollkuhn J.
     
  
  
  3. 
     Wu M. V.
     
  
  
  4. 
     Fraser E. J.
     
  
  
  5. 
     Soderborg T.
     
  
  
  6. 
     Tan S.
     
  
  
  7. 
     Honda S.-I.
     
  
  
  8. 
     Harada N.
     
  
  
  9. 
     Shah N. M
     
  
  
  
  
  2010
  
  
  
  The Androgen Receptor Governs the Execution, but Not Programming, of Male Sexual and Territorial Behaviors
  
  
  
  
  Neuron
  
  
  66
  
  :260–272
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (6)
  
  
  
  
  1. 
     Wu M. V.
     
  
  
  2. 
     Manoli D. S.
     
  
  
  3. 
     Fraser E. J.
     
  
  
  4. 
     Coats J. K.
     
  
  
  5. 
     Tollkuhn J.
     
  
  
  6. 
     Honda S.-I.
     
  
  
  7. 
     Harada N.
     
  
  
  8. 
     Shah N. M.
     
  
  
  
  
  2009
  
  
  
  Estrogen Masculinizes Neural Pathways and Sex-Specific Behaviors
  
  
  
  
  Cell
  
  
  139
  
  :61–72
  
  
  
  
  Google Scholar
  
  
  


• 
  
  (7)
  
  
  
  
  1. 
     Wu M. V.
     
  
  
  2. 
     Tollkuhn J
     
  
  
  
  
  2017
  
  
  
  Estrogen Receptor Alpha Is Required in GABAergic, but Not Glutamatergic, Neurons to Masculinize Behavior
  
  
  
  
  Hormones and Behavior
  
  
  95
  
  :3–12
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (8)
  
  
  
  
  1. 
     Hou H.
     
  
  
  2. 
     Uusküla-Reimand L.
     
  
  
  3. 
     Makarem M.
     
  
  
  4. 
     Corre C.
     
  
  
  5. 
     Saleh S.
     
  
  
  6. 
     Metcalf A.
     
  
  
  7. 
     Goldenberg A.
     
  
  
  8. 
     Palmert M. R.
     
  
  
  9. 
     Wilson M. D
     
  
  
  
  
  2017
  
  
  
  Gene Expression Profiling of Puberty-Associated Genes Reveals Abundant Tissue and Sex-Specific Changes across Postnatal Development
  
  
  
  
  Hum Mol Genet
  
  
  26
  
  :3585–3599
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (9)
  
  
  
  
  1. 
     Guo J.
     
  
  
  2. 
     Nie X.
     
  
  
  3. 
     Giebler M.
     
  
  
  4. 
     Mlcochova H.
     
  
  
  5. 
     Wang Y.
     
  
  
  6. 
     Grow E. J.
     
  
  
  7. 
     Kim R.
     
  
  
  8. 
     Tharmalingam M.
     
  
  
  9. 
     Matilionyte G.
     
  
  
  10. 
      Lindskog C.
      
  
  
  11. 
      Carrell D. T.
      
  
  
  12. 
      Mitchell R. T.
      
  
  
  13. 
      Goriely A.
      
  
  
  14. 
      Hotaling J. M.
      
  
  
  15. 
      Cairns B. R
      
  
  
  
  
  2020
  
  
  
  The Dynamic Transcriptional Cell Atlas of Testis Development during Human Puberty
  
  
  
  
  Cell Stem Cell
  
  
  26
  
  :262–276
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (10)
  
  
  
  
  1. 
     Clemens L. G.
     
  
  
  2. 
     Yang L.-Y
     
  
  
  
  
  2000
  
  
  
  MPOA Lesions Affect Female Pacing of Copulation in Rats
  
  
  
  
  Behavioral Neuroscience
  
  
  114
  
  :1191–1202
  
  
  
  
  Google Scholar
  
  
  


• 
  
  (11)
  
  
  
  
  1. 
     Heimer L.
     
  
  
  2. 
     Larsson K
     
  
  
  
  
  1967
  
  
  
  Impairment of Mating Behavior in Male Rats Following Lesions in the Preoptic-Anterior Hypothalamic Continuum
  
  
  
  
  Brain Research
  
  
  3
  
  :248–263
  
  
  
  https://doi.org/10.1016/0006-8993(67)90076-5
  
  
  Google Scholar
  
  
  


• 
  
  (12)
  
  
  
  
  1. 
     Oomura Y.
     
  
  
  2. 
     Aou S.
     
  
  
  3. 
     Koyama Y.
     
  
  
  4. 
     Fujita I.
     
  
  
  5. 
     Yoshimatsu H
     
  
  
  
  
  1988
  
  
  
  Central Control of Sexual Behavior
  
  
  
  
  Brain Research Bulletin
  
  
  20
  
  :863–870
  
  
  
  https://doi.org/10.1016/0361-9230(88)90103-7
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (13)
  
  
  
  
  1. 
     Rosenblatt J. S.
     
  
  
  2. 
     Ceus K
     
  
  
  
  
  1998
  
  
  
  Estrogen Implants in the Medial Preoptic Area Stimulate Maternal Behavior in Male Rats
  
  
  
  
  Horm Behav
  
  
  33
  
  :23–30
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (14)
  
  
  
  
  1. 
     Sano K.
     
  
  
  2. 
     Nakata M.
     
  
  
  3. 
     Musatov S.
     
  
  
  4. 
     Morishita M.
     
  
  
  5. 
     Sakamoto T.
     
  
  
  6. 
     Tsukahara S.
     
  
  
  7. 
     Ogawa S
     
  
  
  
  
  2016
  
  
  
  Pubertal Activation of Estrogen Receptor α in the Medial Amygdala Is Essential for the Full Expression of Male Social Behavior in Mice
  
  
  
  
  PNAS
  
  
  113
  
  :7632–7637
  
  
  
  
  Google Scholar
  
  
  


• 
  
  (15)
  
  
  
  
  1. 
     Sano K.
     
  
  
  2. 
     Tsuda M. C.
     
  
  
  3. 
     Musatov S.
     
  
  
  4. 
     Sakamoto T.
     
  
  
  5. 
     Ogawa S
     
  
  
  
  
  2013
  
  
  
  Differential Effects of Site-Specific Knockdown of Estrogen Receptor α in the Medial Amygdala, Medial Pre-Optic Area, and Ventromedial Nucleus of the Hypothalamus on Sexual and Aggressive Behavior of Male Mice
  
  
  
  
  European Journal of Neuroscience
  
  
  37
  
  :1308–1319
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (16)
  
  
  
  
  1. 
     Wei Y.-C.
     
  
  
  2. 
     Wang S.-R.
     
  
  
  3. 
     Jiao Z.-L.
     
  
  
  4. 
     Zhang W.
     
  
  
  5. 
     Lin J.-K.
     
  
  
  6. 
     Li X.-Y.
     
  
  
  7. 
     Li S.-S.
     
  
  
  8. 
     Zhang X.
     
  
  
  9. 
     Xu X.-H
     
  
  
  
  
  2018
  
  
  
  Medial Preoptic Area in Mice Is Capable of Mediating Sexually Dimorphic Behaviors Regardless of Gender
  
  
  
  
  Nature Communications
  
  
  9
  
  :1–15
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (17)
  
  
  
  
  1. 
     Fang Y.-Y.
     
  
  
  2. 
     Yamaguchi T.
     
  
  
  3. 
     Song S. C.
     
  
  
  4. 
     Tritsch N. X.
     
  
  
  5. 
     Lin D
     
  
  
  
  
  2018
  
  
  
  A Hypothalamic Midbrain Pathway Essential for Driving Maternal Behaviors
  
  
  
  
  Neuron
  
  
  98
  
  :192–207
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (18)
  
  
  
  
  1. 
     Karigo T.
     
  
  
  2. 
     Kennedy A.
     
  
  
  3. 
     Yang B.
     
  
  
  4. 
     Liu M.
     
  
  
  5. 
     Tai D.
     
  
  
  6. 
     Wahle I. A.
     
  
  
  7. 
     Anderson D. J
     
  
  
  
  
  2020
  
  
  
  Distinct Hypothalamic Control of Same- and Opposite-Sex Mounting Behaviour in Mice
  
  
  
  
  Nature
  
  :1–6
  
  
  
  
  Google Scholar
  
  
  


• 
  
  (19)
  
  
  
  
  1. 
     McHenry J. A.
     
  
  
  2. 
     Otis J. M.
     
  
  
  3. 
     Rossi M. A.
     
  
  
  4. 
     Robinson J. E.
     
  
  
  5. 
     Kosyk O.
     
  
  
  6. 
     Miller N. W.
     
  
  
  7. 
     McElligott Z. A.
     
  
  
  8. 
     Budygin E. A.
     
  
  
  9. 
     Rubinow D. R.
     
  
  
  10. 
      Stuber G. D
      
  
  
  
  
  2017
  
  
  
  Hormonal Gain Control of a Medial Preoptic Area Social Reward Circuit
  
  
  
  
  Nature Neuroscience
  
  
  20
  
  :449–458
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (20)
  
  
  
  
  1. 
     Wu Z.
     
  
  
  2. 
     Autry A. E.
     
  
  
  3. 
     Bergan J. F.
     
  
  
  4. 
     Watabe-Uchida M.
     
  
  
  5. 
     Dulac C. G
     
  
  
  
  
  2014
  
  
  
  Galanin Neurons in the Medial Preoptic Area Govern Parental Behaviour
  
  
  
  
  Nature
  
  
  509
  
  :325–330
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (21)
  
  
  
  
  1. 
     Brock O.
     
  
  
  2. 
     Mees C. D.
     
  
  
  3. 
     Bakker J
     
  
  
  
  
  2015
  
  
  
  Hypothalamic Expression of Oestrogen Receptor α and Androgen Receptor Is Sex-, Age- and Region-Dependent in Mice
  
  
  
  
  Journal of Neuroendocrinology
  
  
  27
  
  :264–276
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (22)
  
  
  
  
  1. 
     Moffitt J. R.
     
  
  
  2. 
     Bambah-Mukku D.
     
  
  
  3. 
     Eichhorn S. W.
     
  
  
  4. 
     Vaughn E.
     
  
  
  5. 
     Shekhar K.
     
  
  
  6. 
     Perez J. D.
     
  
  
  7. 
     Rubinstein N. D.
     
  
  
  8. 
     Hao J.
     
  
  
  9. 
     Regev A.
     
  
  
  10. 
      Dulac C.
      
  
  
  11. 
      Zhuang X
      
  
  
  
  
  2018
  
  
  
  Molecular, Spatial, and Functional Single-Cell Profiling of the Hypothalamic Preoptic Region
  
  
  
  
  Science
  
  
  362
  
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (23)
  
  
  
  
  1. 
     Gegenhuber B.
     
  
  
  2. 
     Wu M. V.
     
  
  
  3. 
     Bronstein R.
     
  
  
  4. 
     Tollkuhn J
     
  
  
  
  
  2022
  
  
  
  Gene Regulation by Gonadal Hormone Receptors Underlies Brain Sex Differences
  
  
  
  
  Nature
  
  
  606
  
  :153–159
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (24)
  
  
  
  
  1. 
     Kim D.-W.
     
  
  
  2. 
     Yao Z.
     
  
  
  3. 
     Graybuck L. T.
     
  
  
  4. 
     Kim T. K.
     
  
  
  5. 
     Nguyen T. N.
     
  
  
  6. 
     Smith K. A.
     
  
  
  7. 
     Fong O.
     
  
  
  8. 
     Yi L.
     
  
  
  9. 
     Koulena N.
     
  
  
  10. 
      Pierson N.
      
  
  
  11. 
      Shah S.
      
  
  
  12. 
      Lo L.
      
  
  
  13. 
      Pool A.-H.
      
  
  
  14. 
      Oka Y.
      
  
  
  15. 
      Pachter L.
      
  
  
  16. 
      Cai L.
      
  
  
  17. 
      Tasic B.
      
  
  
  18. 
      Zeng H.
      
  
  
  19. 
      Anderson D. J
      
  
  
  
  
  2019
  
  
  
  Multimodal Analysis of Cell Types in a Hypothalamic Node Controlling Social Behavior
  
  
  
  
  Cell
  
  
  179
  
  :713–728
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (25)
  
  
  
  
  1. 
     Welch J. D.
     
  
  
  2. 
     Kozareva V.
     
  
  
  3. 
     Ferreira A.
     
  
  
  4. 
     Vanderburg C.
     
  
  
  5. 
     Martin C.
     
  
  
  6. 
     Macosko E. Z
     
  
  
  
  
  2019
  
  
  
  Single-Cell Multi-Omic Integration Compares and Contrasts Features of Brain Cell Identity
  
  
  
  
  Cell
  
  
  177
  
  :1873–1887
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (26)
  
  
  
  
  1. 
     Ogawa S.
     
  
  
  2. 
     Lubahn D. B.
     
  
  
  3. 
     Korach K. S.
     
  
  
  4. 
     Pfaff D. W
     
  
  
  
  
  1997
  
  
  
  Behavioral Effects of Estrogen Receptor Gene Disruption in Male Mice
  
  
  
  
  PNAS
  
  
  94
  
  :1476–1481
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (27)
  
  
  
  
  1. 
     Cheong R. Y.
     
  
  
  2. 
     Czieselsky K.
     
  
  
  3. 
     Porteous R.
     
  
  
  4. 
     Herbison A. E
     
  
  
  
  
  2015
  
  
  
  Expression of ESR1 in Glutamatergic and GABAergic Neurons Is Essential for Normal Puberty Onset, Estrogen Feedback, and Fertility in Female Mice
  
  
  
  
  J. Neurosci
  
  
  35
  
  :14533–14543
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (28)
  
  
  
  
  1. 
     Karigo T.
     
  
  
  2. 
     Deutsch D
     
  
  
  
  
  2022
  
  
  
  Flexibility of Neural Circuits Regulating Mating Behaviors in Mice and Flies
  
  
  
  
  Front. Neural Circuits
  
  
  16
  
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (29)
  
  
  
  
  1. 
     Hrvatin S.
     
  
  
  2. 
     Sun S.
     
  
  
  3. 
     Wilcox O. F.
     
  
  
  4. 
     Yao H.
     
  
  
  5. 
     Lavin-Peter A. J.
     
  
  
  6. 
     Cicconet M.
     
  
  
  7. 
     Assad E. G.
     
  
  
  8. 
     Palmer M. E.
     
  
  
  9. 
     Aronson S.
     
  
  
  10. 
      Banks A. S.
      
  
  
  11. 
      Griffith E. C.
      
  
  
  12. 
      Greenberg M. E
      
  
  
  
  
  2020
  
  
  
  Neurons That Regulate Mouse Torpor
  
  
  
  
  Nature
  
  
  583
  
  :115–121
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (30)
  
  
  
  
  1. 
     Wolf F. A.
     
  
  
  2. 
     Hamey F. K.
     
  
  
  3. 
     Plass M.
     
  
  
  4. 
     Solana J.
     
  
  
  5. 
     Dahlin J. S.
     
  
  
  6. 
     Göttgens B.
     
  
  
  7. 
     Rajewsky N.
     
  
  
  8. 
     Simon L.
     
  
  
  9. 
     Theis F. J
     
  
  
  
  
  2019
  
  
  
  PAGA: Graph Abstraction Reconciles Clustering with Trajectory Inference through a Topology Preserving Map of Single Cells
  
  
  
  
  Genome Biology
  
  
  20
  
  :59
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (31)
  
  
  
  
  1. 
     Aibar S.
     
  
  
  2. 
     González-Blas C. B.
     
  
  
  3. 
     Moerman T.
     
  
  
  4. 
     Huynh-Thu V. A.
     
  
  
  5. 
     Imrichova H.
     
  
  
  6. 
     Hulselmans G.
     
  
  
  7. 
     Rambow F.
     
  
  
  8. 
     Marine J.-C.
     
  
  
  9. 
     Geurts P.
     
  
  
  10. 
      Aerts J.
      
  
  
  11. 
      Oord J.
      
  
  
  12. 
      Atak Z. K.
      
  
  
  13. 
      Wouters J.
      
  
  
  14. 
      Aerts S
      
  
  
  
  
  2017
  
  
  
  SCENIC: Single-Cell Regulatory Network Inference and Clustering
  
  
  
  
  Nat Methods
  
  
  14
  
  :1083–1086
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (32)
  
  
  
  
  1. 
     Stuart T.
     
  
  
  2. 
     Butler A.
     
  
  
  3. 
     Hoffman P.
     
  
  
  4. 
     Hafemeister C.
     
  
  
  5. 
     Papalexi E.
     
  
  
  6. 
     Mauck W. M.
     
  
  
  7. 
     Hao Y.
     
  
  
  8. 
     Stoeckius M.
     
  
  
  9. 
     Smibert P.
     
  
  
  10. 
      Satija R
      
  
  
  
  
  2019
  
  
  
  Comprehensive Integration of Single-Cell Data
  
  
  
  
  Cell
  
  
  177
  
  :1888–1902
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (33)
  
  
  
  
  1. 
     Heumos L.
     
  
  
  2. 
     Schaar A. C.
     
  
  
  3. 
     Lance C.
     
  
  
  4. 
     Litinetskaya A.
     
  
  
  5. 
     Drost F.
     
  
  
  6. 
     Zappia L.
     
  
  
  7. 
     Lücken M. D.
     
  
  
  8. 
     Strobl D. C.
     
  
  
  9. 
     Henao J.
     
  
  
  10. 
      Curion F.
      
  
  
  11. 
      Schiller H. B.
      
  
  
  12. 
      Theis F. J
      
  
  
  
  
  2023
  
  
  
  Best Practices for Single-Cell Analysis across Modalities
  
  
  
  
  Nat Rev Genet
  
  
  24
  
  :550–572
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (34)
  
  
  
  
  1. 
     Cao J.
     
  
  
  2. 
     Spielmann M.
     
  
  
  3. 
     Qiu X.
     
  
  
  4. 
     Huang X.
     
  
  
  5. 
     Ibrahim D. M.
     
  
  
  6. 
     Hill A. J.
     
  
  
  7. 
     Zhang F.
     
  
  
  8. 
     Mundlos S.
     
  
  
  9. 
     Christiansen L.
     
  
  
  10. 
      Steemers F. J.
      
  
  
  11. 
      Trapnell C.
      
  
  
  12. 
      Shendure J
      
  
  
  
  
  2019
  
  
  
  The Single-Cell Transcriptional Landscape of Mammalian Organogenesis
  
  
  
  
  Nature
  
  
  566
  
  :496–502
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (35)
  
  
  
  
  1. 
     Qiu X.
     
  
  
  2. 
     Mao Q.
     
  
  
  3. 
     Tang Y.
     
  
  
  4. 
     Wang L.
     
  
  
  5. 
     Chawla R.
     
  
  
  6. 
     Pliner H. A.
     
  
  
  7. 
     Trapnell C
     
  
  
  
  
  2017
  
  
  
  Reversed Graph Embedding Resolves Complex Single-Cell Trajectories
  
  
  
  
  Nat. Methods
  
  
  14
  
  :979–982
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (36)
  
  
  
  
  1. 
     Pape M.
     
  
  
  2. 
     Miyagi M.
     
  
  
  3. 
     Ritz S. A.
     
  
  
  4. 
     Boulicault M.
     
  
  
  5. 
     Richardson S. S.
     
  
  
  6. 
     Maney D. L
     
  
  
  
  
  2024
  
  
  
  Sex Contextualism in Laboratory Research: Enhancing Rigor and Precision in the Study of Sex-Related Variables
  
  
  
  
  Cell
  
  
  187
  
  :1316–1326
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (37)
  
  
  
  
  1. 
     Weinreb C.
     
  
  
  2. 
     Rodriguez-Fraticelli A.
     
  
  
  3. 
     Camargo F. D.
     
  
  
  4. 
     Klein A. M
     
  
  
  
  
  2020
  
  
  
  Lineage Tracing on Transcriptional Landscapes Links State to Fate during Differentiation
  
  
  
  
  Science
  
  
  367
  
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (38)
  
  
  
  
  1. 
     Burkhardt D. B.
     
  
  
  2. 
     Stanley J. S.
     
  
  
  3. 
     Tong A.
     
  
  
  4. 
     Perdigoto A. L.
     
  
  
  5. 
     Gigante S. A.
     
  
  
  6. 
     Herold K. C.
     
  
  
  7. 
     Wolf G.
     
  
  
  8. 
     Giraldez A. J.
     
  
  
  9. 
     Dijk D.
     
  
  
  10. 
      Krishnaswamy S
      
  
  
  
  
  2021
  
  
  
  Quantifying the Effect of Experimental Perturbations at Single-Cell Resolution
  
  
  
  
  Nature Biotechnology
  
  :1–11
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (39)
  
  
  
  
  1. 
     Moon K. R.
     
  
  
  2. 
     Dijk D.
     
  
  
  3. 
     Wang Z.
     
  
  
  4. 
     Gigante S.
     
  
  
  5. 
     Burkhardt D. B.
     
  
  
  6. 
     Chen W. S.
     
  
  
  7. 
     Yim K.
     
  
  
  8. 
     Elzen A.
     
  
  
  9. 
     Hirn M. J.
     
  
  
  10. 
      Coifman R. R.
      
  
  
  11. 
      Ivanova N. B.
      
  
  
  12. 
      Wolf G.
      
  
  
  13. 
      Krishnaswamy S
      
  
  
  
  
  2019
  
  
  
  Visualizing Structure and Transitions in High-Dimensional Biological Data
  
  
  
  
  Nature Biotechnology
  
  
  37
  
  :1482–1492
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (40)
  
  
  
  
  1. 
     Buchholtz L. J.
     
  
  
  2. 
     Ghitani N.
     
  
  
  3. 
     Lam R. M.
     
  
  
  4. 
     Licholai J. A.
     
  
  
  5. 
     Chesler A. T.
     
  
  
  6. 
     Ryba N. J. P
     
  
  
  
  
  2021
  
  
  
  Decoding Cellular Mechanisms for Mechanosensory Discrimination
  
  
  
  
  Neuron
  
  
  109
  
  :285–298
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (41)
  
  
  
  
  1. 
     Choi H. M. T.
     
  
  
  2. 
     Schwarzkopf M.
     
  
  
  3. 
     Fornace M. E.
     
  
  
  4. 
     Acharya A.
     
  
  
  5. 
     Artavanis G.
     
  
  
  6. 
     Stegmaier J.
     
  
  
  7. 
     Cunha A.
     
  
  
  8. 
     Pierce N. A
     
  
  
  
  
  2018
  
  
  
  Third-Generation in Situ Hybridization Chain Reaction: Multiplexed, Quantitative, Sensitive, Versatile, Robust
  
  
  
  
  Development
  
  
  145
  
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (42)
  
  
  
  
  1. 
     Knoedler J. R.
     
  
  
  2. 
     Inoue S.
     
  
  
  3. 
     Bayless D. W.
     
  
  
  4. 
     Yang T.
     
  
  
  5. 
     Tantry A.
     
  
  
  6. 
     Davis C.
     
  
  
  7. 
     Leung N. Y.
     
  
  
  8. 
     Parthasarathy S.
     
  
  
  9. 
     Wang G.
     
  
  
  10. 
      Alvarado M.
      
  
  
  11. 
      Rizvi A. H.
      
  
  
  12. 
      Fenno L. E.
      
  
  
  13. 
      Ramakrishnan C.
      
  
  
  14. 
      Deisseroth K.
      
  
  
  15. 
      Shah N. M
      
  
  
  
  
  2022
  
  
  
  A Functional Cellular Framework for Sex and Estrous Cycle-Dependent Gene Expression and Behavior
  
  
  
  
  Cell
  
  
  185
  
  :654–671
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (43)
  
  
  
  
  1. 
     Xu X.
     
  
  
  2. 
     Coats J. K.
     
  
  
  3. 
     Yang C. F.
     
  
  
  4. 
     Wang A.
     
  
  
  5. 
     Ahmed O. M.
     
  
  
  6. 
     Alvarado M.
     
  
  
  7. 
     Izumi T.
     
  
  
  8. 
     Shah N. M
     
  
  
  
  
  2012
  
  
  
  Modular Genetic Control of Sexually Dimorphic Behaviors
  
  
  
  
  Cell
  
  
  148
  
  :596–607
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (44)
  
  
  
  
  1. 
     Piekarski D. J.
     
  
  
  2. 
     Boivin J. R.
     
  
  
  3. 
     Wilbrecht L
     
  
  
  
  
  2017
  
  
  
  Ovarian Hormones Organize the Maturation of Inhibitory Neurotransmission in the Frontal Cortex at Puberty Onset in Female Mice
  
  
  
  
  Curr Biol
  
  
  27
  
  :1735–1745
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (45)
  
  
  
  
  1. 
     Ishii K. K.
     
  
  
  2. 
     Hashikawa K.
     
  
  
  3. 
     Chea J.
     
  
  
  4. 
     Yin S.
     
  
  
  5. 
     Fox R. E.
     
  
  
  6. 
     Kan S.
     
  
  
  7. 
     Shah M.
     
  
  
  8. 
     Zhou C.
     
  
  
  9. 
     Navarrete J.
     
  
  
  10. 
      Murry A. D.
      
  
  
  11. 
      Szelenyi E.
      
  
  
  12. 
      Golden S.
      
  
  
  13. 
      Stuber G
      
  
  
  
  
  2024
  
  
  
  Post-Ejaculatory Inhibition of Female Sexual Drive via Heterogeneous Neuronal Ensembles in the Medial Preoptic Area
  
  
  
  
  bioRxiv
  
  :2711
  
  
  
  
  Google Scholar
  
  
  


• 
  
  (46)
  
  
  
  
  1. 
     Shimogori T.
     
  
  
  2. 
     Lee D. A.
     
  
  
  3. 
     Miranda-Angulo A.
     
  
  
  4. 
     Yang Y.
     
  
  
  5. 
     Wang H.
     
  
  
  6. 
     Jiang L.
     
  
  
  7. 
     Yoshida A. C.
     
  
  
  8. 
     Kataoka A.
     
  
  
  9. 
     Mashiko H.
     
  
  
  10. 
      Avetisyan M.
      
  
  
  11. 
      Qi L.
      
  
  
  12. 
      Qian J.
      
  
  
  13. 
      Blackshaw S
      
  
  
  
  
  2010
  
  
  
  A Genomic Atlas of Mouse Hypothalamic Development
  
  
  
  
  Nature Neuroscience
  
  
  13
  
  :767–775
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (47)
  
  
  
  
  1. 
     Romanov R. A.
     
  
  
  2. 
     Tretiakov E. O.
     
  
  
  3. 
     Kastriti M. E.
     
  
  
  4. 
     Zupancic M.
     
  
  
  5. 
     Häring M.
     
  
  
  6. 
     Korchynska S.
     
  
  
  7. 
     Popadin K.
     
  
  
  8. 
     Benevento M.
     
  
  
  9. 
     Rebernik P.
     
  
  
  10. 
      Lallemend F.
      
  
  
  11. 
      Nishimori K.
      
  
  
  12. 
      Clotman F.
      
  
  
  13. 
      Andrews W. D.
      
  
  
  14. 
      Parnavelas J. G.
      
  
  
  15. 
      Farlik M.
      
  
  
  16. 
      Bock C.
      
  
  
  17. 
      Adameyko I.
      
  
  
  18. 
      Hökfelt T.
      
  
  
  19. 
      Keimpema E.
      
  
  
  20. 
      Harkany T
      
  
  
  
  
  2020
  
  
  
  Molecular Design of Hypothalamus Development
  
  
  
  
  Nature
  
  :1–7
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (48)
  
  
  
  
  1. 
     Gegenhuber B.
     
  
  
  2. 
     Tollkuhn J
     
  
  
  
  
  2020
  
  
  
  Signatures of Sex: Sex Differences in Gene Expression in the Vertebrate Brain
  
  
  
  
  WIREs Developmental Biology
  
  
  9
  
  :e348
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (49)
  
  
  
  
  1. 
     McCarthy M. M.
     
  
  
  2. 
     Arnold A. P.
     
  
  
  3. 
     Ball G. F.
     
  
  
  4. 
     Blaustein J. D.
     
  
  
  5. 
     De Vries G. J
     
  
  
  
  
  2012
  
  
  
  Sex Differences in the Brain: The Not So Inconvenient Truth
  
  
  
  
  Journal of Neuroscience
  
  
  32
  
  :2241–2247
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (50)
  
  
  
  
  1. 
     Shansky R. M
     
  
  
  
  
  2019
  
  
  
  Are Hormones a “Female Problem” for Animal Research?
  
  
  
  
  Science
  
  
  364
  
  :825–826
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (51)
  
  
  
  
  1. 
     Chen J. Y.
     
  
  
  2. 
     Campos C. A.
     
  
  
  3. 
     Jarvie B. C.
     
  
  
  4. 
     Palmiter R. D
     
  
  
  
  
  2018
  
  
  
  Parabrachial CGRP Neurons Establish and Sustain Aversive Taste Memories
  
  
  
  
  Neuron
  
  
  100
  
  :891–899
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (52)
  
  
  
  
  1. 
     Heymann G.
     
  
  
  2. 
     Jo Y. S.
     
  
  
  3. 
     Reichard K. L.
     
  
  
  4. 
     McFarland N.
     
  
  
  5. 
     Chavkin C.
     
  
  
  6. 
     Palmiter R. D.
     
  
  
  7. 
     Soden M. E.
     
  
  
  8. 
     Zweifel L. S
     
  
  
  
  
  2020
  
  
  
  Synergy of Distinct Dopamine Projection Populations in Behavioral Reinforcement
  
  
  
  
  Neuron
  
  
  105
  
  :909–920
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (53)
  
  
  
  
  1. 
     Hashikawa Y.
     
  
  
  2. 
     Hashikawa K.
     
  
  
  3. 
     Rossi M. A.
     
  
  
  4. 
     Basiri M. L.
     
  
  
  5. 
     Liu Y.
     
  
  
  6. 
     Johnston N. L.
     
  
  
  7. 
     Ahmad O. R.
     
  
  
  8. 
     Stuber G. D
     
  
  
  
  
  2020
  
  
  
  Transcriptional and Spatial Resolution of Cell Types in the Mammalian Habenula
  
  
  
  
  Neuron
  
  
  106
  
  :743–758
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (54)
  
  
  
  
  1. 
     Mayer C.
     
  
  
  2. 
     Acosta-Martinez M.
     
  
  
  3. 
     Dubois S. L.
     
  
  
  4. 
     Wolfe A.
     
  
  
  5. 
     Radovick S.
     
  
  
  6. 
     Boehm U.
     
  
  
  7. 
     Levine J. E
     
  
  
  
  
  2010
  
  
  
  Timing and Completion of Puberty in Female Mice Depend on Estrogen Receptor α-Signaling in Kisspeptin Neurons
  
  
  
  
  Proceedings of the National Academy of Sciences
  
  
  107
  
  :22693–22698
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (55)
  
  
  
  
  1. 
     Korsunsky I.
     
  
  
  2. 
     Millard N.
     
  
  
  3. 
     Fan J.
     
  
  
  4. 
     Slowikowski K.
     
  
  
  5. 
     Zhang F.
     
  
  
  6. 
     Wei K.
     
  
  
  7. 
     Baglaenko Y.
     
  
  
  8. 
     Brenner M.
     
  
  
  9. 
     Loh P.
     
  
  
  10. 
      Raychaudhuri S.
      
  
  
  
  
  2019
  
  
  
  Fast, Sensitive and Accurate Integration of Single-Cell Data with Harmony
  
  
  
  
  Nat Methods
  
  
  16
  
  :1289–1296
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (56)
  
  
  
  
  1. 
     Huynh-Thu V. A.
     
  
  
  2. 
     Irrthum A.
     
  
  
  3. 
     Wehenkel L.
     
  
  
  4. 
     Geurts P
     
  
  
  
  
  2010
  
  
  
  Inferring Regulatory Networks from Expression Data Using Tree-Based Methods
  
  
  
  
  PLOS One
  
  
  5
  
  :e12776
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (57)
  
  
  
  
  1. 
     Morabito S.
     
  
  
  2. 
     Miyoshi E.
     
  
  
  3. 
     Michael N.
     
  
  
  4. 
     Shahin S.
     
  
  
  5. 
     Martini A. C.
     
  
  
  6. 
     Head E.
     
  
  
  7. 
     Silva J.
     
  
  
  8. 
     Leavy K.
     
  
  
  9. 
     Perez-Rosendahl M.
     
  
  
  10. 
      Swarup V
      
  
  
  
  
  2021
  
  
  
  Single-Nucleus Chromatin Accessibility and Transcriptomic Characterization of Alzheimer’s Disease
  
  
  
  
  Nat Genet
  
  :1–13
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (58)
  
  
  
  
  1. 
     Chen E. Y.
     
  
  
  2. 
     Tan C. M.
     
  
  
  3. 
     Kou Y.
     
  
  
  4. 
     Duan Q.
     
  
  
  5. 
     Wang Z.
     
  
  
  6. 
     Meirelles G. V.
     
  
  
  7. 
     Clark N. R.
     
  
  
  8. 
     Ma’ayan A
     
  
  
  
  
  2013
  
  
  
  Enrichr: Interactive and Collaborative HTML5 Gene List Enrichment Analysis Tool
  
  
  
  
  BMC Bioinformatics
  
  
  14
  
  :128
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (59)
  
  
  
  
  1. 
     Kuleshov M. V.
     
  
  
  2. 
     Jones M. R.
     
  
  
  3. 
     Rouillard A. D.
     
  
  
  4. 
     Fernandez N. F.
     
  
  
  5. 
     Duan Q.
     
  
  
  6. 
     Wang Z.
     
  
  
  7. 
     Koplev S.
     
  
  
  8. 
     Jenkins S. L.
     
  
  
  9. 
     Jagodnik K. M.
     
  
  
  10. 
      Lachmann A.
      
  
  
  11. 
      McDermott M. G.
      
  
  
  12. 
      Monteiro C. D.
      
  
  
  13. 
      Gundersen G. W.
      
  
  
  14. 
      Ma’ayan A
      
  
  
  
  
  2016
  
  
  
  Enrichr: A Comprehensive Gene Set Enrichment Analysis Web Server 2016 Update
  
  
  
  
  Nucleic Acids Res
  
  
  44
  
  :W90–W97
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (60)
  
  
  
  
  1. 
     DePasquale E. A. K.
     
  
  
  2. 
     Schnell D. J.
     
  
  
  3. 
     Van Camp P.-J.
     
  
  
  4. 
     Valiente-Alandí Í.
     
  
  
  5. 
     Blaxall B. C.
     
  
  
  6. 
     Grimes H. L.
     
  
  
  7. 
     Singh H.
     
  
  
  8. 
     Salomonis N
     
  
  
  
  
  2019
  
  
  
  DoubletDecon: Deconvoluting Doublets from Single-Cell RNA-Sequencing Data
  
  
  
  
  Cell Reports
  
  
  29
  
  :1718–1727
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (61)
  
  
  
  
  1. 
     Butler A.
     
  
  
  2. 
     Hoffman P.
     
  
  
  3. 
     Smibert P.
     
  
  
  4. 
     Papalexi E.
     
  
  
  5. 
     Satija R
     
  
  
  
  
  2018
  
  
  
  Integrating Single-Cell Transcriptomic Data across Different Conditions, Technologies, and Species
  
  
  
  
  Nat. Biotechnol
  
  
  36
  
  :411–420
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (62)
  
  
  
  
  1. 
     Haghverdi L.
     
  
  
  2. 
     Lun A. T. L.
     
  
  
  3. 
     Morgan M. D.
     
  
  
  4. 
     Marioni J. C
     
  
  
  
  
  2018
  
  
  
  Batch Effects in Single-Cell RNA-Sequencing Data Are Corrected by Matching Mutual Nearest Neighbors
  
  
  
  
  Nature Biotechnology
  
  
  36
  
  :421–427
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (63)
  
  
  
  
  1. 
     Tasic B.
     
  
  
  2. 
     Yao Z.
     
  
  
  3. 
     Graybuck L. T.
     
  
  
  4. 
     Smith K. A.
     
  
  
  5. 
     Nguyen T. N.
     
  
  
  6. 
     Bertagnolli D.
     
  
  
  7. 
     Goldy J.
     
  
  
  8. 
     Garren E.
     
  
  
  9. 
     Economo M. N.
     
  
  
  10. 
      Viswanathan S.
      
  
  
  11. 
      Penn O.
      
  
  
  12. 
      Bakken T.
      
  
  
  13. 
      Menon V.
      
  
  
  14. 
      Miller J.
      
  
  
  15. 
      Fong O.
      
  
  
  16. 
      Hirokawa K. E.
      
  
  
  17. 
      Lathia K
      
  
  
  18. 
      Tieu M.
      
  
  
  19. 
      Larsen R.
      
  
  
  20. 
      Casper T.
      
  
  
  21. 
      Barkan E.
      
  
  
  22. 
      Kroll M.
      
  
  
  23. 
      Parry S.
      
  
  
  24. 
      Shapovalova N. V.
      
  
  
  25. 
      Hirschstein D.
      
  
  
  26. 
      Pendergraft J.
      
  
  
  27. 
      Sullivan H. A.
      
  
  
  28. 
      Kim T. K.
      
  
  
  29. 
      Szafer A.
      
  
  
  30. 
      Dee N.
      
  
  
  31. 
      Groblewski P.
      
  
  
  32. 
      Wickersham I.
      
  
  
  33. 
      Cetin A.
      
  
  
  34. 
      Harris J. A.
      
  
  
  35. 
      Levi B. P.
      
  
  
  36. 
      Sunkin S. M.
      
  
  
  37. 
      Madisen L.
      
  
  
  38. 
      Daigle T. L.
      
  
  
  39. 
      Looger L.
      
  
  
  40. 
      Bernard A.
      
  
  
  41. 
      Phillips J.
      
  
  
  42. 
      Lein E.
      
  
  
  43. 
      Hawrylycz M.
      
  
  
  44. 
      Svoboda K.
      
  
  
  45. 
      Jones A. R.
      
  
  
  46. 
      Koch C.
      
  
  
  47. 
      Zeng H.
      
  
  
  
  
  2018
  
  
  
  Shared and Distinct Transcriptomic Cell Types across Neocortical Areas
  
  
  
  
  Nature
  
  
  563
  
  :72–78
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (64)
  
  
  
  
  1. 
     Zeisel A.
     
  
  
  2. 
     Hochgerner H.
     
  
  
  3. 
     Lönnerberg P.
     
  
  
  4. 
     Johnsson A.
     
  
  
  5. 
     Memic F.
     
  
  
  6. 
     Zwan J.
     
  
  
  7. 
     Häring M.
     
  
  
  8. 
     Braun E.
     
  
  
  9. 
     Borm L. E.
     
  
  
  10. 
      La Manno G.
      
  
  
  11. 
      Codeluppi S.
      
  
  
  12. 
      Furlan A.
      
  
  
  13. 
      Lee K.
      
  
  
  14. 
      Skene N.
      
  
  
  15. 
      Harris K. D.
      
  
  
  16. 
      Hjerling-Leffler J.
      
  
  
  17. 
      Arenas E.
      
  
  
  18. 
      Ernfors P.
      
  
  
  19. 
      Marklund U.
      
  
  
  20. 
      Linnarsson S
      
  
  
  
  
  2018
  
  
  
  Molecular Architecture of the Mouse Nervous System
  
  
  
  
  Cell
  
  
  174
  
  :999–1014
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (65)
  
  
  
  
  1. 
     Zappia L.
     
  
  
  2. 
     Oshlack A
     
  
  
  
  
  2018
  
  
  
  Clustering Trees: A Visualization for Evaluating Clusterings at Multiple Resolutions
  
  
  
  
  GigaScience
  
  
  7
  
  :giy083
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (66)
  
  
  
  
  1. 
     Davie K.
     
  
  
  2. 
     Janssens J.
     
  
  
  3. 
     Koldere D.
     
  
  
  4. 
     Waegeneer M. D.
     
  
  
  5. 
     Pech U.
     
  
  
  6. 
     Kreft Ł.
     
  
  
  7. 
     Aibar S.
     
  
  
  8. 
     Makhzami S.
     
  
  
  9. 
     Christiaens V.
     
  
  
  10. 
      González-Blas C. B.
      
  
  
  11. 
      Poovathingal S.
      
  
  
  12. 
      Hulselmans G.
      
  
  
  13. 
      Spanier K. I.
      
  
  
  14. 
      Moerman T.
      
  
  
  15. 
      Vanspauwen B.
      
  
  
  16. 
      Geurs S.
      
  
  
  17. 
      Voet T.
      
  
  
  18. 
      Lammertyn J.
      
  
  
  19. 
      Thienpont B.
      
  
  
  20. 
      Liu S.
      
  
  
  21. 
      Konstantinides N.
      
  
  
  22. 
      Fiers M.
      
  
  
  23. 
      Verstreken P.
      
  
  
  24. 
      Aerts S
      
  
  
  
  
  2018
  
  
  
  A Single-Cell Transcriptome Atlas of the Aging Drosophila Brain
  
  
  
  
  Cell
  
  
  174
  
  :982–998
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (67)
  
  
  
  
  1. 
     Otis J. M.
     
  
  
  2. 
     Namboodiri V. M. K.
     
  
  
  3. 
     Matan A. M.
     
  
  
  4. 
     Voets E. S.
     
  
  
  5. 
     Mohorn E. P.
     
  
  
  6. 
     Kosyk O.
     
  
  
  7. 
     McHenry J. A.
     
  
  
  8. 
     Robinson J. E.
     
  
  
  9. 
     Resendez S. L.
     
  
  
  10. 
      Rossi M. A.
      
  
  
  11. 
      Stuber G. D
      
  
  
  
  
  2017
  
  
  
  Prefrontal Cortex Output Circuits Guide Reward Seeking through Divergent Cue Encoding
  
  
  
  
  Nature
  
  
  543
  
  :103–107
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  (68)
  
  
  
  
  1. 
     Rossi M. A.
     
  
  
  2. 
     Basiri M. L.
     
  
  
  3. 
     McHenry J. A.
     
  
  
  4. 
     Kosyk O.
     
  
  
  5. 
     Otis J. M.
     
  
  
  6. 
     Munkhof H. E.
     
  
  
  7. 
     Bryois J.
     
  
  
  8. 
     Hübel C.
     
  
  
  9. 
     Breen G.
     
  
  
  10. 
      Guo W.
      
  
  
  11. 
      Bulik C. M.
      
  
  
  12. 
      Sullivan P. F.
      
  
  
  13. 
      Stuber G. D
      
  
  
  
  
  2019
  
  
  
  Obesity Remodels Activity and Transcriptional State of a Lateral Hypothalamic Brake on Feeding
  
  
  
  
  Science
  
  
  364
  
  :1271–1274
  
  
  
  
  PubMed
  
  
  Google Scholar
  
  
  


• 
  
  
  
  1. 
     Hashikawa K
     
  
  
  
  
  2025
  
  
  
  scRNAseq of MPOA during development
  
  
  
  
  NCBI Gene Expression Omnibus
  
  
  
  ID GSE172177
  
  
  
  
  


• 
  
  
  
  1. 
     Moffitt J.R
     
  
  
  2. 
     et al.
     
  
  
  
  
  2018
  
  
  
  Data from: Molecular, spatial and functional single-cell profiling of the hypothalamic preoptic region
  
  
  
  
  Dryad