Sex significantly impacts the function of major depression-linked variants in vivo

Genome-wide association studies have discovered blocks of common variants—likely transcriptional-regulatory—associated with major depressive disorder (MDD), though the functional subset and their biological impacts remain unknown. Likewise, why depression occurs in females more frequently than males is unclear. We therefore tested the hypothesis that risk-associated functional variants interact with sex and produce greater impact in female brains. We developed methods to directly measure regulatory variant activity and sex interactions using massively parallel reporter assays (MPRAs) in the mouse brain in vivo, in a cell type-specific manner. We measured activity of >1,000 variants from >30 MDD loci, identifying extensive sex-by-allele effects in mature hippocampal neurons and suggesting sex-differentiated impacts of genetic risk may underlie sex bias in disease. Unbiased informatics approaches indicated that functional MDD variants recurrently disrupt sex hormone receptor binding sequences. We confirmed this with MPRAs in neonatal brains, comparing brains undergoing the masculinizing hormone surge to hormonally-quiescent juveniles. Our study provides novel insights into the influence of age, biological sex, and cell type on regulatory-variant function, and provides a framework for in vivo parallel assays to functionally define interactions between organismal variables like sex and regulatory variation. One-Sentence Summary Massively parallel assays in vivo identified extensive functional and sex-interacting common variants in depression risk loci.


Introduction
Major depressive disorder (MDD) is a profoundly disruptive and sometimes lethal disorder, affecting women 2-3 times more frequently than men across countries and cultures (1).
Sex differences are present across multiple levels of the disease, from symptom profiles (2) and effective drug classes (3) to brain-wide gene expression (4, 5). Genome-wide association studies(GWASes) have identified dozens of linkage regions each containing numerous singlenucleotide polymorphisms (SNPs) associated with MDD, demonstrating its heritability (6-8).
More recently, sex-by-genotype (SxG) analyses of large GWAS cohorts have revealed that MDD 5 risk loci are for men and women, yet these loci explain up to 4-fold greater MDD heritability in females (9,10). These findings suggest that sex interacts with a common pool of SNPs to attenuate or amplify the MDD risk they confer. However, disease-associated SNPs are seldom found in protein-coding space, complicating prediction of their molecular consequences. Instead, these SNPs are found in probable regulatory elements (REs), including transcriptional-regulatory 10 sequences predicted from measures such as chromatin marks, accessibility, and conformation. Specific brain regions and cell types are enriched for such measures at-and in putative regulatory target genes of-MDD-associated loci, including the hippocampus (7,11) and excitatory neurons (12)(13)(14)(15)(16), suggesting sites of action for these REs. In particular, there has been long-standing interest in the hippocampus regarding both MDD pathology and sex differences in 15 the brain. Hippocampal volume reductions in MDD patients have been widely reported (17).
Moreover, the hippocampus is subject to influences of sex from perinatal (18,19) to adult life, presenting in MDD as sex differences in hippocampal volume loss (20) and gene expression (4).
Determining the identity of functional SNPs from MDD-associated regions is the first key step toward understanding the biological perturbations resulting from risk genotypes, which 20 can in turn enable inference of dysregulated target genes and shared regulatory programs involved across loci. However, studies connecting MDD-associated SNPs to gene expression in brain tissue, even those that considered sex effects (21), have been limited to indirect and imperfect indicators of function (e.g., chromatin state), or confounded by linkage disequilibrium (e.g., for expression quantitative trait loci (eQTLs)). In contrast, direct measurement of regulatory output of common variants associated with disease has largely been restricted to the in vitro setting. The large-scale in vitro identification of functional regulatory variants has been made possible by massively parallel reporter assays (MPRAs), a method for functionally detecting activity from thousands of REs (and their variants) simultaneously. In brief, MPRAs 5 adapt a traditional reporter assay paradigm-placing REs upstream of an optically measured reporter (e.g., luciferase)-but adds a unique, RE-identifying "barcode" sequence to the reporter's 3' untranslated region, enabling quantification of activity for thousands of REs simultaneously by RNA barcode sequencing. MPRAs have enabled identification of trait-and disease-associated SNPs affecting REs in culturable, disease-relevant cell types in vitro (22)(23)(24)(25). 10 However, the complexities of cell types interacting in the brain and of the sex hormonal milieu cannot readily be emulated ex vivo.
We overcome these prior limitations in functional regulatory SNP (rSNP) identification to interrogate the biological contexts under which MDD rSNPs act by delivering an MPRA library of MDD-associated variants (26) into the adult mouse hippocampus in vivo. Building on 15 prior brain MPRAs of enhancers (27) and an RE variant (28), our approach greatly extends in vivo MPRA methods to identify rSNPs and their sex interactions, including those which are cell type-specific. First, we combined MPRA with translating ribosome affinity purification (TRAP) to simultaneously identify MDD rSNPs in both excitatory neurons and the broader hippocampus; these experiments utilized mice of both sexes, enabling us to additionally test the hypothesis that 20 rSNPs are subject to sex-by-genotype (SxG) interactions. Finally, to further characterize the potential role of circulating hormones in sex-differentiated rSNP activity, and to functionally replicate predicted fetal brain RE enrichments suggesting a role for MDD SNPs during circuit organization (15,(29)(30)(31)(32), we likewise delivered the library to the mouse brain in utero. This allowed us to identify rSNPs neonatally, coinciding with a testosterone surge and critical period for establishing sex-specific brain circuitry (33), and test for loss of SxG effects in juveniles, when hormonal influences are quiescent. In sum, we illustrate that MPRAs can be leveraged in vivo to directly identify not only functional variants, but their context-dependence on age, sex, and cell type, while demonstrating that all three of these factors have substantial impacts on 5 MDD-associated regulatory variation.

Combining Translating Ribosome Affinity Purification (TRAP) with MPRAs
Adeno-associated viruses (AAVs) have long been used to evaluate the activity of single 10 REs in the brain, most recently in MPRA-like designs to screen multiple REs in parallel with RNA-seq-based quantification (27,34,35). This suggests it may also be possible to adapt AAV-MPRAs to study functional consequences of RE variants associated with disease. As MDD genetic risk is enriched in neuronal REs, we tested the feasibility of combining a cell typespecific profiling method, TRAP, with MPRA to attain measurements of RE activity specifically 15 in neurons. We generated 4 small AAV9 libraries ( Fig. 1A-C) expressing dsRed under the control of the hsp68 minimal promoter, with full-length human promoters (denoted pGene) with documented expression in neurons (human pCAMK2A), astrocytes (pGFAP), or all cells (pPGK2), each carrying unique 3' untranslated region (UTR) barcodes for quantification by RNA-seq. We first individually confirmed cell-type specificity by immunofluorescence (IF) 20 ( Fig. 1D-F), and then delivered (36) a mixed pool of the four barcoded AAV9 libraries into the brain of postnatal day 2 (P2) neuron-specific TRAP mice (Snap25-Rpl10a-eGFP) (37). TRAP was then used to compare total brain and neuronal activity levels of the three promoters and the hsp68 promoter alone. We prepared MPRA sequencing libraries from 1) the delivered AAV pool (DNA), 2) total brain (input) RNA, and 3) TRAP (neuronal) RNA and assessed expression of each barcode by calculating a ratio of RNA counts to DNA counts. Results were highly replicable at the level of barcodes (Fig. 1G-H) in both RNA fractions. Moreover, neuronal TRAP fractions demonstrated increased expression of barcodes driven by pCAMK2A and lower 5 expression of pGFAP-paired barcodes (Fig. 1I). These results demonstrated that cell typespecific effects, even of relatively small magnitudes, can be detected using a combined MPRA-TRAP approach. We then turned to applying this method to test the effects of human variants associated with psychiatric disease. Amplicons were cloned into an AAV plasmid. C) Further restriction cloning added a reporter cassette containing a minimal hsp68 promoter, dsRed, and an RNA-stabilizing 3' UTR hepatitis 5 E "woodchuck" (WPRE) element. Barcoded pools with each promoter were packaged into AAV9 separately. (D-F) IF of P27 mouse brain after P2 injection with a single AAV9 barcode We delivered the AAV9 library bilaterally into the hippocampus of Vglut1-TRAP mice ages P60-P80 (n=6 per sex), followed by hippocampal dissection and TRAP ( Fig. 2A) to identify rSNPs and their shared regulatory features (Fig. 2B). IF of hippocampi from additional mice confirmed robust hippocampal expression of the dsRed reporter 28 days after injection (Fig. 2C).
We first confirmed by qPCR that TRAP RNA was depleted for glial marker genes several-fold 5 as expected (Fig. 2D), then conducted MPRA sequencing ( Table S1) and analysis of both input RNA (hippocampus) and TRAP RNA (Vglut1 + ) for each sample. After quality control, n=5 male and n=5 female hippocampal samples (both the input and 10 TRAP RNAs) were retained for MPRA analysis covering ~1,000 SNPs (see Online Methods). These showed replicability at the level of barcode counts (pairwise Pearson's r 0.7-0.89) (Fig.   2E), RE-level expression (r 0.85-0.96) (Fig. 2F), and allelic fold changes within each condition (r 0.90-0.94) (Fig. 2G), confirming we were able to reliably measure SNP-mediated regulatory effects and differences from a defined cell type in vivo. 15

rSNPs in the adult hippocampus and hippocampal Vglut1 + neurons
We first assessed ~1,000 SNPs for allelic effects in each individual sex and RNA fraction, using linear mixed models (LMMs) fitting barcode expression as a function of allele with random barcode effects (see Online Methods and Fig. S2-S4). We calculated empirical p (pemp) values using 50,000 simulated 'allelic' comparisons between subsets of random barcodes 20 coupled to the minimal promoter alone (26, 51) to account for technical and barcode-mediated noise. Significance was called at FDR-corrected pemp<0.2, a stringency comparable to a recent study of sex-interacting eQTLs (21). Supplemental Data S1 provides allelic beta (log2FC) values and significance status for variants at 10 different significance thresholds.
In total hippocampus, we identified 36 (male) and 31 (female) rSNPs, 34 and 31 of which 25 were from MDD loci, respectively. While male and female total hippocampi had similar numbers of rSNPs, we observed a striking sex difference in the number of rSNPs in Vglut1 + cells specifically-only 7 (male) compared to 58 (female), indicating that within excitatory neurons, a higher proportion of MDD SNPs have discernible allelic effects in females. . Moreover, all 7 male rSNPs were also functional in females. Notable rSNPs from ≥1 condition included rs2563323 and rs250427, putative brain and hippocampal (52) eQTLs for SRA1, a noncoding RNA which activates nuclear receptors even in the absence of their ligand (53). Also notable were rs301806-an MDD GWAS index SNP (40)-and rs301807 ( Fig. 3A-B), both of which  10 target genes in human tissues ascertained by Hi-C. A plot of rSNP effects, colored by most significant condition, is embedded with its x-axis in human genomic (hg19) coordinates; chromatin contacts between the SNPs and distal gene promoters are illustrated below (56)(57)(58)(59)(60)(61). Curve heights correspond to -log10pemp for the plotted rSNP. iPSC: Induced pluripotent stem cell; DG: dentate gyrus; Ctx: cortex. *: pemp-derived FDR < 0.25; **: <0.2; ‡ < 0.15; ‡ ‡ < 0.1; & < 0.05. 5 We therefore performed combined-sex LMM analysis of the Vglut1 + and hippocampus results, identifying 41 sex-allele interaction rSNPs in each. Notably, while only 1 SxG rSNP was shared between total hippocampus and Vglut1 + , sex-interacting SNPs originated from the same GWAS loci across both analyses. The tag locus rs1193510 (39), for example, contained 11 sex- 10 interacting rSNPs, including several unique to hippocampus ( Fig. 3C) or Vglut1 + (Fig. 3D). This region is rich in human neural chromatin contacts (56-61), implicating several target genes of the identified rSNPs. Interestingly, SxG effects in hippocampus and Vglut1 + segregated into distinct portions of this LD region (Fig. 3E). We additionally examined three SNPs in our assay that were recently reported significant in sex-genotype interaction GWASes of over 500 traits 15 (62). One of these reported variants, rs2400075, was found to have several significant sexinteracting associations to body traits; this variant is located near LIN28B-which shows sexdifferential expression in mouse norepinephrine neurons (63)-and was a near-significant SxG rSNP in Vglut1 + (0.20<FDR<0.25). 20 We next asked whether there were any shared transcriptional-regulatory mechanisms underlying MDD rSNP effects in the hippocampus. We tested whether rSNPs perturbed specific transcription factor (TF) binding motifs more frequently than expected by chance (defined by their rates in rSNPs vs. non-effect SNPs in the assay) (26). We assessed motif disruptions using motifbreakR (64) and RSAT var-tools (65) (see Online Methods) defining rSNPs at a nominal 25 LMM pemp of 0.05. This resulted in sets of 80-110 MPRA-identified rSNPs per condition (pemp 13 FDR levels of 0.22-0.28; Data S2A), ensuring adequate depth for enrichment analysis-akin to other methods which loosen significance thresholds to analyze higher-order relationships of granular molecular measures, e.g., polygenic risk scoring or gene ontology. To refine these results, we filtered the enriched TFs (FDR<0.05) to those with altered putative binding sites ≥4 rSNPs and expressed in Genotype-Tissue Expression atlas (v8) hippocampus in the 5 corresponding sex.

Transcriptional-regulatory systems shared across hippocampal rSNPs
Altogether, we identified 34 enriched TFs in male total hippocampus rSNPs and 19 in female. The hippocampal TFs identified were largely distinctive between sexes; for example, KLF family TFs were unique to male hippocampal rSNPs, while nuclear receptor (NR) TFs were mostly unique to female rSNPs (Fig. 4A). Among Vglut1 + rSNPs, we identified 8 TFs in female 10 and 16 in male, many of which were shared (e.g., DLX1, POU3F1/2/3). POU3F2 has been previously shown to be a highly centralized, cross-disorder hub gene in postmortem brain coexpression analysis by PsychENCODE (66).
To understand integrative biological functions of these TF sets, we utilized the tool Enrichr (67) for each TF set, including both tissues per sex, to identify putative upstream 15 regulators, co-interacting TFs (protein-protein interactions (PPIs)), enriched disease gene sets, ontologies, and brain regions expressing the TFs (Data S2B-C). The most striking enrichments

Transcriptional-regulatory systems implicated in SxG interactions at MDD rSNPs in hippocampus
Using a similar approach to that above, we investigated TFs enriched at sex-allele interaction rSNPs (pemp<0.05, corresponding to FDRs of 0.34-0.37) relative to sex-agnostic 10 rSNPs (combined-sexes LMM allele main effect pemp<0.05), and separately enriched to nonfunctional SNPs, then combined the two enrichment outputs to form TF sets (see Online Methods) (Fig. 4A). We identified only 8 TFs enriched among total hippocampal SxG rSNPs, including ZNF410 and five TFs with a shared binding motif for FOX(C1/J2/J3/O4/P1).

Identification of rSNPs in developing whole mouse brain
As sex differences in brain structure and transcriptional regulation are established in part by the effects of sex hormones, including the perinatal testosterone surge (33), we sought to investigate whether MDD risk variants were subject to sex-differential regulation during early 15 development. To be able to assess the brain during this period, we delivered the AAV library intracerebroventricularly to embryonic day 15 (E15) mice , followed by whole brain collection at postnatal day 0 (P0) and P10 (see Online Methods). P0-while not amenable to regionallytargeted viral assays-is in the midst of the testosterone surge, during which both acute activational effects and transcriptional-regulatory organizational effects occur; by contrast, P10 20 is a timepoint where sex hormones are effectively absent in normal development. We first verified by IF that dsRed expression was detectable at P0 and P10 following in utero delivery. Clear, widespread reporter expression was apparent at both timepoints despite the relatively short incubation time (Fig. 5A-D; Fig. S6). IF at P10 demonstrated prominent 15 expression of the reporter in the hippocampus (Fig. 5D)-a structure neither present at E15 nor well-developed at P0-consistent with prior observations of AAV9 expression ultimately occurring in hippocampus when delivered to the perinatal brain (73). We subsequently collected the whole brain (except cerebellum) for RNA isolation and MPRA sequencing. Additionally, we isolated DNA from n=4 brain samples (one per age and sex) to profile the transduced library 20 contents, verifying that the distribution of delivered MPRA barcodes was similar both between replicates (r 0.86-0.94) and to input virus (r 0.88-0.91) (Fig. E). Ultimately, we analyzed 15 samples for P0 (6 female, 9 male) and 13 for P10 (6 female, 7 male; see Supplemental Text).
Replicates from each condition had well-correlated barcode counts (PCC 0.86-0.98) and RE expression values (PCC 0.58-0.96) (Fig. 5F-H). 25 Within single sexes at P0, we identified 5 rSNPs in females and 12 rSNPs in males, respectively (pemp FDR < 0.2), 4 of which were shared between conditions with consistent effect direction. By contrast, we identified 105 female and 72 male rSNPs at P10, with 42 rSNPs identified in both conditions with consistent direction of effect. 4 of these shared rSNPs were the same rSNPs shared between P0 sexes, with consistent effect direction in all four conditions (two of which are illustrated, Fig. 5I, J). Three of these four rSNPs also have rich chromatin contact evidence supporting gene regulatory roles in fetal, adult, and cultured human neural cell types (for the illustrated rSNPs, Fig. 5K and 5L, respectively), consistent with their detection as rSNPs in whole brain tissue, highlighting in utero MPRA delivery as a robust method for detecting functional variation in the developing brain. 5

Sex-allele interactions are widespread neonatally but absent during hormonal quiescence
To identify SxG interactions occurring during neurodevelopment, we tested for SxG interactions as before, now within age groups. Again, the minimal promoter-only control was not sex-differentially expressed (t-test of barcode expression, p>0.1) (Fig. S7). At pemp FDR <0.2, we identified 31 rSNPs with sex interactions in the P0 brain (e.g., Fig. 5M-N). By contrast, we 10 identified no SxG interactions at pemp FDR<0.25 among the 930 analyzed SNPs at P10. We confirmed this result-despite similar n and inter-sample correlations to P0-by repeating the SxG analysis with the least variable 5 samples per sex (removing two males and one female) (Data S3). 15

development
We examined single-sex, single-age rSNP sets for enrichment of TF motif perturbations using motifbreakR and RSAT approaches as before, and again using nominally significant (pemp <0.05) rSNPs as the set tested for enrichment. We did not filter enriched TFs for expression in analogous human tissue, as fetal whole-brain gene expression profiles are unavailable; we 20 instead required a more stringent minimum of 5 rSNPs to be found at significantly enriched (FDR<0.05) TF motifs.
Male-specific transcription factors again included ZBTB7A and NR4A3 (Fig. 6A). Femalespecific TFs included a variety of zinc finger TFs, Krüpell-like factors (KLFs), endothelialdevelopmental regulator SOX17, and the neurodevelopmental TF SMAD3 (Fig. 6A). 5 Despite the absence of sex interactions in the P10 brain, we also found largely distinct sets of TFs in each sex at this age (only 4 shared) (Fig. 6A) We then looked at these TF sets as before to identify convergent regulators and functions among them ( Fig. 6B; Data S2D-E). Of the few enriched annotations for male P0 TFs, we notably found that 4 corresponded to genes downregulated by Rara knockout, and 6 were putative regulatory targets of PBX3, which shows widespread subcortical and midbrain 15 expression in E18.5 and P4 mouse brain (74). Female P0 TFs were enriched in Allen Atlas expression signatures for cortical layers 1, 3, and 5, and, as found in hippocampus, in PPI targets of HDAC2 and ESR1. P10 male TFs showed the greatest extent of overlap (12/42) with generegulatory targets of ZBTB7A and were enriched for PPI targets of HDAC2, RXRA and RARA.
The only enrichment found for P10 female TFs was via a modest (3/15) set of FOXA1 regulatory 20 targets.

Transcriptional-regulatory systems implicated in sex-differential rSNPs postnatally
Given the absence of FDR-corrected sex-genotype interactions at P10, we only analyzed P0 interaction rSNPs for TF motif perturbation enrichment. TFs enriched at interaction SNPs were comprised largely of EGR, KLF, and SP family members, as well as PAX5 and neurodevelopmental factor SMAD3 (Fig. 6A). Annotation of SxG TFs revealed a broader extent of hormonal roles in functional variation than observed in either sex alone at P0: we found SxG 5 TFs were again enriched for PPI targets of ESR1 (as had been P0 female TFs), but additionally enriched in PPI targets of AR and ESR2; the latter was otherwise absent from gene set analyses of TFs (Fig. 6B, Data S2D-E).

neurons
Finally, we wanted to test the hypothesis that the number of SxG interactions from the assayed SNPs exceeded what was expected by chance. We randomly scrambled sex labels and repeated our analyses to define a distribution for the null expectation regarding number of SxG rSNPs at a pemp FDR of 0.2. While our P0 SxG findings did not exceed chance expectations ( Fig.   15 6C), our P10 results were consistent with a true absence of sex interactions (Fig. 6D). Most strikingly, the number of adult total hippocampal (Fig. 6E) and Vglut1 + (Fig. 6F) SxG interactions we identified were both significantly greater than would be anticipated. Color key shared with panel B. *: locus from all-female MDD GWAS (42); **: near-genome-wide significant locus in males with MDD developing after age 50 (41); ‡: collapsed results from two loci near CELF4 with tags less than 150kb apart; ‡ ‡: collapsed results 5 from LD partners of five tag SNPs, comprising two GWAS significant tags and several weaker association peaks covering a span of ~5Mb around the gene TCF4 (6).

Risk loci are characterized by multiple functional variants across age and cell type
Overall, we found that it was the norm for loci to contain multiple functional, context-10 dependent SNPs, in contrast to the concept of a singular "causal variant" driving a given GWAS association. Altogether, our analyses identified 280 rSNPs from 31 LD regions (28 depressionassociated), with up to 13 rSNPs in a single locus found in a single condition (P10 female, tag SNP rs11707582 (6)) (Fig. 6G). In terms of age and cell type, we identified dozens of rSNPs with allele or SxG effects specific to one timepoint or tissue region: 26 P0-specific, 101 P10- 15 specific, 34 total hippocampus-specific, and 55 specific to Vglut1 + cells of the hippocampus.

Discussion and Conclusion
We have directly measured sex-genotype interactions across MDD in the adult hippocampus and sexually differentiating brain, thus demonstrating the existence of a genetic component of sex differences in MDD and the regulatory architecture underlying these 25 differences across space, time, and genome. We have uncovered functional differences between sexes at particular MDD-associated SNPs in the hippocampus, its excitatory neurons, and the brain during its sexual differentiation, expanding on observed genetic and clinical sex differences in MDD from general heritability to direct identification of sex-interacting variants.
Our data are supported by a wide variety of orthogonal datasets, including reporter assays and human brain epigenomic datasets beyond those used for variant prioritization. MDDassociated SNP rs1467013 was previously demonstrated to be functional in a classical luciferase 5 reporter assay in three different cell lines (75). A prior in vitro MPRA identified rs301807, but not rs301806, as an rSNP (26), while both were identified as functional here. This highlights the importance of in vivo context for obtaining relevant insights about functional variation within GWAS loci-in this case, revealing two rSNPs in close (~2kb) proximity that likely influence expression of the same target (RERE). 10 Human datasets further support the translatability of our mouse approach in identifying regulatory SNPs and sex interactions. Four rSNPs identified here, rs7244124 (Vglut1 + SxG), rs76931017 (P10 female), rs827187 (P10 female), and rs4482931 (male hippocampus) were recently identified as chromatin accessibility QTLs in human midfetal neural progenitors and/or neurons (76), consistent with the allele-differential regulatory activity we observed. Intriguingly, 15 an early attempt to identify sex-interacting GWAS loci for MDD (77) found suggestive significance for rs1345818, near TMEM161B and MEF2C, a locus which has since been sexagnostically associated to MDD (6). Our assay did not include the reported sex-interacting GWAS tag, but we identified three variants in LD with rs1345818 that showed sex interactions, confirming that this risk locus indeed has sex-dependent regulatory activity: rs1814149 20 (hippocampal SxG, FDR<0.2, R 2 with rs1345818=0.67), rs5869417 (hippocampal SxG, FDR<0.1, R 2 =0.21), and rs6452770 (Vglut1 + SxG, FDR<0.2, R 2 =0.69).
Downstream analyses aimed at identifying regulatory programs involved across rSNPs provided both support for the rSNP findings themselves, while also identifying novel candidate TFs underlying sex interactions at MDD loci. Our SxG-enriched TF sets were especially rich in sex hormone receptors and interactors, consistent with expectations for an in vivo assay detecting sex interactions, and indicating a role of sex hormone receptors in co-regulation of MDD risk 5 variants. Our hippocampal analyses revealed male-specific roles for AR: male rSNPs were enriched for binding sequences of ZMIZ1, an AR co-activator, while an AR co-repressor, ZBTB7A (78), was identified as a shared upstream regulator of these TFs. Notably, ZBTB7A also regulates human non-coding RNA LINC00473 (79), which was recently been demonstrated to have sex-differentiated effects on depressive mouse behaviors when overexpressed in cortex 10 (80). Our TF analyses of P0 SxG rSNPs also identified regulatory programs consistent with the critical period for sexualization of the brain. P0 SxG rSNPs were enriched for TFs interacting with ESR1, ESR2, and AR, consistent with the regulatory landscape necessary for accommodation of sex hormonal signals during the perinatal testosterone surge. Additionally, PAX5 motif disruptions were unique to P0 SxG variants; interestingly the PAX5 motif was 15 recently shown to be enriched in promoters of sex-differentially expressed genes in adult brain (21). Our neurodevelopmental rSNP-enriched TFs likewise recapitulated aspects of recent preclinical studies of sex and depression: or example, P0 female rSNPs were enriched in motifs of endothelial marker SOX17, consistent with demonstration of sex-differential changes in mouse brain vascular permeability after stress (81). 20 Our assay has several limitations. Notably, our P0/P10 assays could not be compared directly to our adult findings to look for sex-by-age effects, as adult experiments were limited to the hippocampus while developmental assays were brain-wide. Unfortunately, it is impractical to region-specifically deliver AAV9 in utero, and hippocampal morphogenesis is postnatal.
Likewise, episomal AAV delivery may not capture all regulatory information of a genomic delivery, though benchmark studies indicate a strong correlation (82). Finally, results might be influenced by use of either a different minimal promoter or longer fragments if nearby elements interact with the rSNPs. However, none of these limitations would be expected to create spurious sex effects, indicating the surprising and widespread SxG interactions are likely to be robust. 5 The presented in vivo MPRA approach indicates that critical biological and environmental factors involved in brain gene regulation and regulatory variation can be studied using a high-fidelity model of development, cell types, and biological signals. Our approach provides a framework for direct, functional study of psychiatric risk genetics and their interactions with biological and environmental factors that are imperfectly modeled in vitro, 10 including cell type, sex, and brain development. This same approach could be readily used in the future to directly identify variants subject to genetic-environmental interactions with other key psychiatric risk factors, such as early life adversity and chronic stress.

References and Notes
Acknowledgments: Scott Lee for assistance with immunofluorescence of pilot hippocampal  Competing interests: Authors declare that they have no competing interests. 10 Data and materials availability: The annotation of library SNPs performed at the time of design are available at https://bitbucket.org/jdlabteam/n2a_atra_mdd_mpra_paper/src/master/Library%20Design %20Epigenomic%20Annotations/. Code is available at https://bitbucket.org/jdlabteam/paper-resources-mdd-in-vivo-mpras/src. Raw sequencing 15 files and tabulated barcode counts per sample will be made available through GEO accession GSE186348 on publication. Outside datasets used for annotation are detailed in Supplementary Materials. 20