Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2022 Jun;606(7912):153-159.
doi: 10.1038/s41586-022-04686-1. Epub 2022 May 4.

Gene regulation by gonadal hormone receptors underlies brain sex differences

B Gegenhuber  1   2 , M V Wu  1 , R Bronstein  1 , J Tollkuhn  3
Affiliations

Gene regulation by gonadal hormone receptors underlies brain sex differences

B Gegenhuber et al. Nature. 2022 Jun.

Abstract

Oestradiol establishes neural sex differences in many vertebrates1-3 and modulates mood, behaviour and energy balance in adulthood4-8. In the canonical pathway, oestradiol exerts its effects through the transcription factor oestrogen receptor-α (ERα)9. Although ERα has been extensively characterized in breast cancer, the neuronal targets of ERα, and their involvement in brain sex differences, remain largely unknown. Here we generate a comprehensive map of genomic ERα-binding sites in a sexually dimorphic neural circuit that mediates social behaviours. We conclude that ERα orchestrates sexual differentiation of the mouse brain through two mechanisms: establishing two male-biased neuron types and activating a sustained male-biased gene expression program. Collectively, our findings reveal that sex differences in gene expression are defined by hormonal activation of neuronal steroid receptors. The molecular targets we identify may underlie the effects of oestradiol on brain development, behaviour and disease.

PubMed Disclaimer

Conflict of interest statement

The authors declare no competing interests.

Figures

Fig. 1
Fig. 1. Genomic targets of ERα in sexually dimorphic neuronal populations.
a, Coronal sections containing sexually dimorphic brain areas used for ERα CUT&RUN. MPOA, medial pre-optic area; BNSTp, posterior bed nucleus of the stria terminalis; MeAp, posterior medial amygdala. b, Line plots (top) and heatmaps (bottom) of mean IgG and ERα CUT&RUN (C&R) CPM ±1 kb around E2-induced ERα CUT&RUN peaks (DiffBind edgeR, Padj < 0.1). The heatmaps are sorted by E2 ERα CUT&RUN signal. Colour scale is counts per million (CPM). Veh, vehicle. c, Cross-tissue ERα comparison, showing the proportion of ERα peaks detected specifically in brain. d, Top Gene Ontology biological process terms associated with genes nearest to brain-specific or shared (≥4 other tissues) ERα CUT&RUN peaks (clusterProfiler, Padj < 0.1). e, Combined sex E2 versus vehicle RNA-seq in BNSTp Esr1+ cells; light grey and red dots (DESeq2, Padj < 0.1), dark grey and red dots (DESeq2, P < 0.01), purple dots (validated by in situ hybridization (ISH)). FC, fold change. Positive FC is E2-upregulated, negative FC is E2-downregulated. f, Images (left panels) and quantitative analysis (right panels) of ISH for select genes induced by E2 in both sexes. Boxplot centre, median; box boundaries, first and third quartiles; whiskers, 1.5 × IQR from boundaries. Two-way analysis of variance: Brinp2 P = 0.0373, Rcn1 P = 0.0307, Enah P = 0.0003, Tle3 P = 0.0001; n = 4 per condition; scale bar, 200 µm. g, MA plot of E2-regulated ATAC–seq peaks in BNSTp Esr1+ cells; red dots are E2-open peaks (DiffBind edgeR, log2[FC] > 1, Padj < 0.05), grey dots are E2-close peaks (DiffBind edgeR, log2[FC] < −1, Padj < 0.05). h, Example ERα peaks at E2-induced genes. Top left number is the y-axis range in CPM. Shaded band indicates peak region. Source Data
Fig. 2
Fig. 2. Sex differences in cell type abundance and gene regulation in BNSTp Esr1+ cells.
a, Uniform manifold approximation and projection (UMAP) visualization of BNSTp Esr1+ snRNA-seq inhibitory neuron clusters, coloured by identity (left), sex (middle) and Esr1 expression (right). b, Proportion of BNSTp Esr1+ nuclei in each BNSTp Esr1+ inhibitory neuron cluster per sex. Higher proportions of i1:Nfix (Padj = 0.002) and i3:Esr2 (Padj = 0.002) neurons are in males than females. Boxplot centre, median; box boundaries, first and third quartiles; whiskers, 1.5 × IQR from boundaries, n = 7, **Padj < 0.01, one-sided, Wilcoxon rank-sum test, adjusted with the Benjamini–Hochberg procedure. c, BNSTp immunofluorescence (IF) staining for GFP (left micrographs) and Nfix (middle micrographs) in P14 female and male Esr1Cre/+;Sun1–GFPlx/+ animals (scale bar, 100 µm), with combined images (right micrographs) and their quantification (boxplots; right). Boxplot centre, median; box boundaries, first and third quartiles; whiskers, 1.5 × IQR from boundaries, n = 6, P = 0.0422, *P < 0.05, two-sided, unpaired t-test. d, Heatmap of median MetaNeighbor area under the receiver operating characteristic curve (AUROC) values for BNSTp Esr1+ clusters and cortical/hippocampal GABAergic neuron subclasses. The colour bar indicates the developmental origin of GABAergic subclasses. CGE, caudal ganglionic eminence; MGE, medial ganglionic eminence. e, Top: heatmap of MetaNeighbor AUROC values for BNSTp and MPOA Esr1+ clusters. Bottom: average expression of i1:Nfix marker genes across BNSTp and MPOA Esr1+ clusters. Dotted box indicates shared identity of i1:Nfix and i20:Gal.Moxd1 cells. n = 297 i20:Gal.Moxd1 cells, 2,459 i1:Nfix cells. Boxplot centre, median; box boundaries, first and third quartiles; whiskers, 1.5 × IQR from boundaries. f, Number of differentially expressed genes (DEGs) between females and males (DESeq2, Padj < 0.1) per BNST neuron snRNA-seq cluster. g, R2 between percentage of TF gene expression and number of sex DEGs per cluster across snRNA-seq clusters. The inset shows correlation for the top-ranked TF gene, Esr1. The error band represents the 95% confidence interval. h, Differential ATAC sites between gonadectomized (GDX), vehicle-treated females and males (top) and gonadally intact females and males (middle). Blue dots (edgeR, log2[FC] > 1, Padj < 0.05), red dots (edgeR, log2[FC] < −1, Padj < 0.05). Bottom: enrichment analysis of sex-biased ATAC peaks at sex DEGs. i, Top: k-means clustering (c1–c4) of differentially accessible ATAC peaks across four conditions(edgeR, Padj < 0.01). Bottom: dotplot showing the percentage of sites per cluster overlapping E2-open ATAC loci and motif enrichment analysis of peaks in each cluster (AME algorithm). ARE, androgen response element. j, Example ATAC peaks in k-means clusters 1 and 2. Top left number is the y-axis range in CPM. Shaded band indicates peak region. Source Data
Fig. 3
Fig. 3. Neonatal ERα genomic binding drives a sustained male-biased gene expression program.
a, Heatmap of P4 BNST Esr1+ ATAC, P0 IgG CUT&RUN and P0 ERα CUT&RUN CPM ±1 kb around 1,605 NE-open and 403 NE-close ATAC peaks (edgeR, Padj < 0.1). ERα+, Sun1–GFP+ nuclei; ERα−, Sun1–GFP− nuclei. b, UMAPs of adult (left) and neonatal (middle left) BNST Esr1+ snRNA-seq clusters; neonatal snRNA-seq clusters coloured by sex (middle right) and time point (right). c, Left: UMAPs of Nfix expression (top left), gene activity score (top right), motif chromVAR deviation score (bottom left) and CUT&RUN chromVAR deviation score (bottom right). Right: neonatal single-nucleus ATAC (snATAC) and adult BNSTp Nfix CUT&RUN tracks at the Nfix locus. Top left number is the y-axis range in CPM. Shaded band indicates peak region. Peak–RNA correlation indicates correlation coefficient for snATAC peaks correlated with Nfix expression. d, Heatmap of differential snATAC CPM between males (M) and females (F) at 1,605 NE-open sites, scaled across snRNA-seq clusters and grouped using k-means clustering. The barplot indicates the percentage of overlap for each k-means cluster with total and E2-induced BNSTp Nfix CUT&RUN peaks. e, Top: number of sex DEGs (MAST, Padj < 0.05) in P4 multiome clusters. Bottom: heatmaps indicating RNA log2[FC] of P4 sex DEGs (left) and Pearson’s correlation coefficient of NE-open (red) and NE-close (blue) ATAC peaks (right) linked to sex DEGs in each cluster. Genes without significant differential expression or correlation coefficients (not significant (NS)) are shown in white. f, Cyp19a1/aromatase expression on P4. g, Left: NE-open ATAC peaks correlating with Lrp1b expression in Cyp19a1− clusters, i2:Tac2 and i12:Esr1. Top left number is the y-axis range in CPM. Shaded band indicates peak region. Right, sex difference in Lrp1b expression in i2:Tac2 (n = 260 female, 153 male, Padj = 2.13 × 10−8), i4:Bnc2 (n = 437 female, 373 male, Padj = 5.62 × 10−37), i12:Esr1 (n = 803 female, 507 male, Padj = 1.09 × 10−12) cells. ***Padj < 0.001, MAST. h, Proportion of P4 sex DEGs detected as sex biased on P14. i, Top: i1:Nfix-specific, NE-open ATAC peaks at Fat1 and Scg2 loci on P4 and P14. Top left number is the y-axis range in CPM. Shaded band indicates peak region. Bottom: Sex difference in i1:Nfix Fat1 and Scg2 expression on P4 (Fat1, Padj = 1.28 × 10−37; Scg2, Padj = 1.54 × 10−46; n = 887 female, 676 male) and P14 (Fat1, Padj = 1.13 × 10−11; Scg2, Padj = 1.52 × 10−5; n = 554 female, 829 male). ***Padj < 0.001, MAST.
Fig. 4
Fig. 4. ERα is required for sustained sex differences in gene expression.
a, UMAPs of adult (top left) and P14 (top right) BNST Vgat+ snRNA-seq clusters; P14 Vgat+ snRNA clusters coloured by group (bottom left) and Esr1+ status (bottom right). b, Top: number of female versus male sex DEGs (MAST, Padj < 0.05) in P14 snRNA clusters (black bar). Number of female versus male sex DEGs detected in female versus male KO comparison (grey bar). Bottom: heatmap of mean expression of i1:Nfix sex DEGs, scaled across control males, control females and conditional ERα-KO males. c, Neonatal ERα activation drives a sustained male-typical gene expression program.
Extended Data Fig. 1
Extended Data Fig. 1. Validation of ERα CUT&RUN in MCF-7 cells.
a, Heatmap of mean MCF-7 ERα CUT&RUN CPM ±1Kb around 12,995 17β-estradiol (E2)-induced MCF-7 ERα ChIP-seq peaks (DiffBind, DESeq2, padj < 0.01) for individual replicates (n = 2 per condition and antibody). b, pA-MNase-cut footprint (CUT&RUNTools) around ESR1 motif sites (FIMO) detected in ERα ChIP-seq peaks. c, MA plots of differential ERα CUT&RUN peaks (DiffBind, DESeq2, padj < 0.1) for ERα antibody #1 (Santa Cruz sc-8002) and ERα antibody #2 (EMD Millipore Sigma 06-935). d, Pearson correlation coefficient of CPM-normalized CUT&RUN signal within the consensus peak matrix across ERα CUT&RUN samples. Red text indicates E2 treatment group. e, ERα CUT&RUN (both antibodies) and ChIP-seq tracks at canonical MCF-7 ERα target genes (TFF1, GREB1). f, (left) Top enriched motifs (AME) and (right) top de novo motifs (DREME) within ERα ChIP-seq peaks, Ab #1 ERα CUT&RUN peaks, and Ab #2 ERα CUT&RUN peaks. % TP = % of peaks called as positive for the indicated motif. De novo motifs were classified into motif families using TomTom.
Extended Data Fig. 2
Extended Data Fig. 2. Additional analysis of adult brain ERα CUT&RUN dataset.
a, MA plot of differential ERα CUT&RUN peaks (DiffBind, edgeR, padj < 0.1) in adult mouse brain. red dots=E2-induced peaks, grey dots=E2-down peaks. b, Heatmap of mean brain ERα CUT&RUN CPM ±1Kb around 1930 E2-induced ERα CUT&RUN peaks (see also Fig. 1b) for individual replicates (n = 2 per condition). c, ESR1 motif footprint in ERα peaks (CUT&RUNTools). d, Top enriched motifs (AME) in (left) E2-induced ERα peaks and (right) E2-down ERα peaks. e, Heatmap of mean brain IgG and ERα CUT&RUN CPM ±1Kb around 185 E2-down ERα peaks. f, Number of overlaps between E2-induced ERα peaks and 7 external ERα ChIP-seq peaksets: intersected peaks of uterus 1 and uterus 2, intersected peaks of liver 1 and liver 2, aorta, efferent ductules, and mammary gland. Red indicates brain-specific ERα peaks. g, Log-odds motif scores (FIMO) for the ESR1 motif (MA0112.3, left) and ESR2 motif (MA0258.2, right) in brain-specific (red) and shared (pink) ERα peaks. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. n = 1304 brain-specific ESR1, 139 shared ESR1, 1276 brain-specific ESR2, 157 shared ESR2. p-values from two-sided, Wilcoxon rank-sum test. h, Top Hugo Gene Nomenclature Committee (HGNC) gene families (clusterProfiler, padj < 0.1) enriched within brain-specific ERα peak-associated genes. i, Top Disease Ontology terms associated with genes nearest to brain-specific ERα peaks (DOSE, padj < 0.1). j, Brain-specific ERα peak-associated genes within each enriched Disease Ontology (DO) term (clusterProfiler, padj < 0.1), colored by term k, Example brain-specific (Cntnap2, Ntrk2), shared (Rybp, Myrip), and disease-associated (Drd3, Htr1a, Grin2b) ERα peaks.
Extended Data Fig. 3
Extended Data Fig. 3. Additional BNSTp ISH validation.
In situ hybridization (ISH) validation of additional BNSTp E2-regulated genes identified by RNA-seq (see also Fig. 1e, f). Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. p-value from 2-way ANOVA test, n = 4, scale = 200 um.
Extended Data Fig. 4
Extended Data Fig. 4. Additional analysis of adult BNSTp ATAC-seq.
a, GFP immunofluorescence staining in an adult male Esr1Cre/+; Sun-GFPlx/+ mouse, scale = 1mm. Inset shows Sun1-GFP signal at nuclear membrane, scale = 10um. b, Fluorescence-activated cell sorting (FACS) gating strategy for isolating BNSTp GFP+ nuclei for ATAC-seq. c, Proportion of E2-open ATAC peaks (red), total E2 ATAC peaks (black), and total Veh ATAC peaks (black) annotated to promoters (±1Kb around TSS), exons, introns, and intergenic regions. E2-open ATAC peaks have a significantly lower proportion of`peaks annotated to gene promoters than total vehicle (11% vs 1%, Fisher’s Exact Test, p = 4.6 x 10−260) and total E2 (11% vs 1%, Fisher’s Exact Test, p = 4.3 x 10−267) peaks. ***p < 0.001. d, Heatmap of mean BNSTp Esr1+ ATAC CPM ±1Kb around 7293 E2-open ATAC peaks for individual female and male replicates (n = 3 per condition) (see also Fig. 1g). e–g, Top Gene Ontology (GO) Biological Process terms (e), HGNC gene families (f), and DO terms (g) enriched within E2-open ATAC peak-associated genes (clusterProfiler, padj < 0.1). h, Top motifs enriched in E2-open ATAC-seq peaks (AME). % of peaks containing motifs determined with FIMO. i, Overlap in E2-open ATAC peaks containing the ESR1 motif (blue) and the ESR2 motif (red), identified using FIMO. The majority of peaks (4434/6479) containing either motif are the same. j, (left) Overlap between brain ERα CUT&RUN peaks and E2-open ATAC peaks (log2FC > 1). (right) Overlap between remaining 777 brain ERα CUT&RUN peaks and log2FC > 0 E2-open ATAC peaks.
Extended Data Fig. 5
Extended Data Fig. 5. Integration of adult RNA-seq, ATAC-seq, and CUT&RUN datasets.
a–b, BETA enrichment of E2-open ATAC peaks (a) and brain ERα CUT&RUN peaks (b) at E2-induced and E2-down genes identified by RNA-seq (DESeq2, p < 0.01) relative to a background of non-differential, expressed genes. c, Top enriched motifs (AME) in E2-open ATAC peaks ±350Kb around E2-down genes (identified with BETA). d, Example E2-open ATAC peaks/ERα peaks at E2-repressed genes, Nr2f1 and Astn2. e, Normalized counts for example genes (Ccdc134, Zfp804b) with a sex-dependent response to E2 treatment. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. n = 4. f, (Left) Volcano plot of sex-dependent, E2-responsive genes; light blue and red dots (DESeq2, padj < 0.1), dark blue and red dots (DESeq2, padj < 0.1). (Right) Mean, normalized expression of sex-dependent, E2-responsive genes (DESeq2, padj < 0.1), grouped by k-means clustering. g, Lack of significant enrichment of E2-open ATAC peaks (top) and ERα peaks (bottom) at sex-dependent, E2-reponsive genes relative to a background of non-differential, expressed genes (BETA). h, Volcano plots of sex-dependent, E2-responsive ERα CUT&RUN peaks (edgeR, padj < 0.1) (left) and ATAC peaks (edgeR, padj < 0.1). i, MA plot of male E2 vs. female E2 ERα CUT&RUN peaks (DiffBind, edgeR, padj < 0.1); red dots=male E2-biased peaks, grey dots=female E2-biased peaks. j, Heatmap of mean ATAC CPM, split by sex and treatment, ±1Kb around E2-induced ERα peaks.
Extended Data Fig. 6
Extended Data Fig. 6. Characterization of a shared BNSTp/MPOA transcriptomic cluster.
a, Dotplot of top marker genes for each adult BNSTp Esr1+ GABAergic cluster (Wilcoxon rank-sum test, padj < 0.05). b, Differentially-expressed genes between the i1:Nfix cluster and the other six BNSTp Esr1+ inhibitory neuron clusters (DESeq2, log2FC > 2, padj < 0.01). c, ISH of adult gonadectomized, Veh-treated male (Esr1) and adult male (Nfix) mouse. Arrows denote Nfix ISH staining in BNSTp (dorsal) and POA (ventral). Scale = 1 mm. d, Co-expression of top i1:Nfix marker genes (St18, Moxd1, Nfix, Cplx3) in individual BICCN cortical and hippocampal scRNA-seq GABAergic clusters, colored by subclass. Co-expression defined as % of cells per cluster with non-zero counts for all 4 marker genes. e, Mean expression of Lamp5+ subclass marker genes (Wilcoxon rank-sum test, avg_log_FC > 0.75, <40% expression in non-Lamp5+ neurons, padj < 0.05) in BNSTp (left) and MPOA (right) Esr1+ clusters, scaled across clusters within each brain region. f, Normalized expression of top marker genes (Moxd1, St18, Nfix, Cplx3, Gpd1, Prox1) shared between i1:Nfix and i20:Gal.Moxd1 (labeled in red). Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. n = 297 i20:Gal.Moxd1 cells, 2459 i1:Nfix cells. g, UMAP visualization of integrated BNSTp and MPOA Esr1+ clusters, demonstrating shared Nfix expression across datasets (see also Fig. 2e). h, GFP (left) and Nfix (middle) immunofluorescence staining in an adult male Esr1Cre/+;Sun-GFPlx/+ mouse. Solid white circle indicates BNSTp; dotted white circle indicates SDN-POA. Scale = 100 um.
Extended Data Fig. 7
Extended Data Fig. 7. Additional analysis of BNSTp sex DEGs and gonadally intact ATAC-seq.
a, Normalized, pseudo-bulk expression of sex DEGs identified within each BNSTp Esr1+ cluster (DESeq2, padj < 0.1). Each heatmap column corresponds to a pseudo-bulk sample (gene counts aggregated across cells in sample). b, Hierarchical clustering of log2FC values for sex DEGs called as significant in at least one BNSTp Esr1+ cluster. Sex DEGs with non-significant differential expression colored in white. c, Example sex DEGs with significant differential expression in a single Esr1+ cluster (Dlg2, Kctd16) and in multiple Esr1+ clusters (Tiparp, Socs2). Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. n = 4-7 female pseudo-bulk replicates, 6-8 male pseudo-bulk replicates. d, Top HGNC gene families (clusterProfiler, padj < 0.1) enriched within female-biased DEGs (left) and male-biased DEGs (right) relative to non-differential, expressed genes. e, Number of sex DEGs per cluster in Esr1+ clusters annotated to the BNST posterior (BNSTp, n = 7 clusters) and anterior (BNSTa, n = 7 clusters) subregions. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. p-value from two-sided, Wilcoxon rank-sum test. f, Dotplot of sex hormone receptor (HR) expression across BNSTp Esr1+ clusters. g, Barplots of (top) male-biased DEGs ranked by male-biased ATAC peak regulatory potential score and (bottom) female-biased DEGs ranked by female-biased ATAC peak regulatory potential score. Higher score indicates higher density of sex-biased ATAC peaks around the TSS of sex DEGs. h, Example sex DEGs (Fkbp5, Epha6) with high density of sex-biased ATAC peaks. *sex-biased ATAC peak. i, Principal component analysis (PCA) of gonadally intact and gonadectomized (GDX), Veh-treated ATAC CPM values within the consensus peak matrix. j, Heatmap of ATAC CPM for gonadally intact (n = 2 per condition) and GDX, Veh-treated ATAC samples (n = 3 per condition) at differential peaks (edgeR, glmQLFTest, padj < 0.01), grouped by k-means clustering (see also Fig. 2i). k, Example differential ATAC peaks in k-means clusters c3 (left) and c4 (right).
Extended Data Fig. 8
Extended Data Fig. 8. P4 ATAC-seq and P0 ERα CUT&RUN analysis.
a, Heatmap of mean ATAC CPM for P4 NV male, NV female, and NE female individual replicates (n = 3 per condition) at differential peaks (edgeR, glmQLFTest, padj < 0.1), grouped by hierarchical clustering (cutree, k = 6). Clusters c3 and c5 correspond to NE-open and NE-close sites, respectively, shown in Fig. 4a. b, Top enriched motifs (AME) in NE-open ATAC peaks. c, (left) Overlap between P4 NE-open ATAC peaks and adult E2-open ATAC peaks (log2FC > 1). (right) Overlap between remaining 509 P4 NE-open ATAC peaks and log2FC > 0 E2-open ATAC peaks. d, Example P4 NE-open ATAC peaks not detected as E2-induced in adult E2-open ATAC peakset. e, MA plot of P0 female E2 vs. female Veh ERα CUT&RUN peaks (DiffBind, DESeq2, padj < 0.01); red dots=E2-induced peaks, grey dots=E2-down peaks. f, Top enriched motifs (AME) in P0 E2-induced ERα peaks. g, Heatmap of mean P0 ERα CUT&RUN CPM ±1Kb around 8102 E2-induced ERα peaks for individual replicates (n = 2 per condition). h, (left) Overlap between P4 NE-open ATAC peaks and P0 E2-induced ERα peaks. (right) Overlap between P4 NE-close ATAC peaks and P0 E2-induced ERα peaks. i, (top) Top GO Biological Process terms (clusterProfiler, padj < 0.1), (middle) DO terms (clusterProfiler, padj < 0.1), and (bottom) HGNC gene families (clusterProfiler, padj < 0.1) enriched within P4 NE-open peak-associated genes. j, Example P0 ERα peaks overlapping P4 NE-open peaks at high-confidence ASD candidate genes, Scn2a1 and Slc6a1.
Extended Data Fig. 9
Extended Data Fig. 9. Comparison of P4 and adult Esr1+ ATAC-seq.
a, (left) Overlap between P4 NE-open ATAC peaks and gonadally intact adult male-biased ATAC peaks. (right) Overlap between P4 NE-close ATAC peaks and gonadally intact adult female-biased ATAC peaks. b, Dotplot of BETA enrichment p-values for P4 NE-open ATAC peaks (top) and NE-close ATAC peaks (bottom) at adult BNST snRNA-seq sex DEGs relative to a background of non-differential, expressed genes (see also Fig. 2h). c, Histogram of mean distance between P4 NE-open peaks and nearest gonadally intact adult male-biased ATAC peak (red line) vs. nearest chromosome-matched, non-differential adult ATAC peak (n = 1000 permutations) (blue histogram). Mean distance between P4 NE-open peaks and adult male-biased peaks is significantly smaller than the expected distribution (Permutation test, p = 0.007). d, Example adult male-biased DEGs (Prlr, Cckar, Pdzrn4, Tiparp) with neighboring P4 NE-open (highlighted in yellow) and adult male-biased ATAC peaks (highlighted in purple).
Extended Data Fig. 10
Extended Data Fig. 10. Additional analysis of neonatal BNST Esr1+ single-nucleus multiome dataset.
a–b, RNA (a) and ATAC (b) quality control (QC) metrics for neonatal (P4, P14) single-nucleus multiome experiments, split by timepoint and sex. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=minimum and maximum values. n = 4265 P14_male, 3148 P14_female, 3128 P4_male, 4295 P4_female. c, UMAPs of de novo clustering of neonatal multiome snRNA data, colored by cluster identity (top left), sex (top right), timepoint (bottom left), and Esr1 expression (bottom right). d, (left) Prediction scores of adult-to-neonatal label transfer for each adult BNST Esr1+ reference cluster, split by timepoint. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=minimum and maximum values. n = 14836 cells. (right) % of nuclei in each neonatal de novo cluster that mapped to each adult BNST Esr1+ cluster. e, UMAPs of neonatal marker gene module expression in neonatal dataset (left) and adult dataset (right). f, (left) Heatmap of pseudo-bulk ATAC CPM at 18783 marker peaks for neonatal multiome clusters. (right) Top three motifs enriched in marker peaks for each multiome cluster. g, (left) Correlation analysis of TF expression and motif accessibility across cells. Putative identity regulator TFs colored in pink. (right) TF RNA expression, activity score, and motif deviation UMAPs of example putative BNST Esr1+ neuron identity regulators, Nr4a2 and Zfhx3 (see also Fig. 3c, d). h, Heatmap of mean cortical IgG and BNSTp Nfix CUT&RUN CPM ±1Kb around 32,578 consensus Nfix peaks. i, Top motifs enriched (AME) in (top) 30,825 mHypoA cell Nfix CUT&RUN peaks (MACS2, q < 0.01) and in (bottom) 32,578 consensus BNSTp Nfix CUT&RUN peaks (MACS2, q < 0.01; peaks intersected across treatment and sex). %TP=% of peaks called as positive for the indicated motif.
Extended Data Fig. 11
Extended Data Fig. 11. Sex differences in single-nucleus multiome dataset.
a, Hierarchical clustering of log2FC values for P4 sex DEGs detected in Esr1+ inhibitory neuron clusters (see also Fig. 3e). Sex DEGs with non-significant differential expression colored in white. b, Neonatal E2 (NE) vs. neonatal vehicle (NV) female nuclear RNA-seq on P4 BNST Esr1+ cells; grey, red dots (DESeq2, padj < 0.1). c, (top) Overlap between NE-induced genes and P4 multiome male-biased genes. (bottom) Overlap between NE-downregulated genes and female-biased genes. d, (left) Pearson’s correlation coefficient values for non-differential (grey) and NE-regulated (gold) ATAC peaks that correlate with P4 sex DEG expression. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=1.5*IQR from boundaries. n = 5169 non-differential, 244 NE-regulated. p-value from two-sided, Wilcoxon rank-sum test. (right) Distance between non-differential (grey) and NE-regulated (gold) ATAC peaks to P4 sex DEG transcription start sites (TSS). p-value from Kolmogorov-Smirnov test. e, Example P4 sex-biased genes that are also NE-regulated, Htr4 (top) and Csgalnact1 (bottom). (left) n = 3, (right) n = 887 i1:Nfix female cells, 676 i1:Nfix male cells, 404 i3:Esr2 female cells, 550 i3:Esr2 male cells. f, Tracks for NE-open ATAC peaks that correlate with NE-regulated, sex-biased targets, Htr4 and Csgalnact1. g, Different NE-open ATAC peaks across i1:Nfix, i3:Esr2, and i4:Bnc2 neurons correlated with a common male-biased target, Arid1b. h, Heatmaps indicating (left) RNA log2FC of P14 sex DEGs and (right) Pearson’s correlation coefficient of NE-open (red) and -close (blue) ATAC peaks linked to sex DEGs within each cluster. Non-significant genes and correlation values colored in white.
Extended Data Fig. 12
Extended Data Fig. 12. Additional analysis of P14 BNST Vgat+ snRNA-seq dataset.
a, RNA QC metrics for P14 BNST Vgat+ snRNA-seq experiment, split by sample. a and b refer to technical replicates. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=minimum and maximum values. n = 6355 ko_male_a, 5614 ko_male_b, 7184 wt_female_a, 6367 wt_female_b, 6881 wt_male_a, 6561 wt_male_b. b, Prediction scores of adult-to-P14 label transfer for each adult BNST Vgat+ reference cluster, split by group. Boxplot center=median, box boundaries=1st and 3rd quartile, whiskers=minimum and maximum values. n = 38962 cells. c, Heatmap of P14 marker gene mean module score in P14 BNST Vgat+ clusters (left) and adult BNST Vgat+ clusters (right). d, Proportion of total P14 Vgat+ nuclei in each Vgat+ cluster, separated by group. Adult male-biased Esr1+ clusters i1:Nfix and i3:Esr2 are indicated in grey. e, Heatmap of mean expression of sex DEGs in (left) i3:Esr2 and (right) i4:Bnc2 clusters, scaled across experimental groups (see also Fig. 4b). f, % of sex DEGs in P14 Vgat+/Esr1+ clusters with that are also detected as sex-biased on P4 in corresponding multiome clusters.
See this image and copyright information in PMC

Comment in

Similar articles

See all similar articles

Cited by

See all "Cited by" articles

References

    1. MacLusky NJ, Naftolin F. Sexual differentiation of the central nervous system. Science. 1981;211:1294–1303. doi: 10.1126/science.6163211. - DOI - PubMed
    1. McEwen BS. Neural gonadal steroid actions. Science. 1981;211:1303–1311. doi: 10.1126/science.6259728. - DOI - PubMed
    1. McCarthy MM. Estradiol and the developing brain. Physiol. Rev. 2008;88:91–124. doi: 10.1152/physrev.00010.2007. - DOI - PMC - PubMed
    1. Wharton W, Gleason CE, Olson SRMS, Carlsson CM, Asthana S. Neurobiological underpinnings of the estrogen - mood relationship. Curr. Psychiatry Rev. 2012;8:247–256. doi: 10.2174/157340012800792957. - DOI - PMC - PubMed
    1. Mauvais-Jarvis F, Clegg DJ, Hevener AL. The role of estrogens in control of energy balance and glucose homeostasis. Endocr. Rev. 2013;34:309–338. doi: 10.1210/er.2012-1055. - DOI - PMC - PubMed

LinkOut - more resources