Whole-Genome and RNA Sequencing Reveal Variation and Transcriptomic Coordination in the Developing Human Prefrontal Cortex

1The human nervous system develops slowly over several decades, starting during embryogenesis and extending postnatally through infancy, childhood, adolescence, and young adulthood. Over this time, myriads of functionally distinct cell types, circuits, and regions are formed. To produce distinct structures and circuits, neural cells are born in an immature state and undergo a variety of molecular and morphological changes as they differentiate, migrate, and establish circuits. Consequently, the characteristics of a given cell and brain region at a given time offer only a snapshot of organogenesis and brain function, necessitating consistent profiling across development.

2The molecular and cellular processes underlying development of the nervous system rely on the diversity of transcripts and their precise spatiotemporal regulation. Functional genomic analyses of the developing human brain have revealed highly dynamic gene expression and epigenetic changes during prenatal and early postnatal development versus comparative stability over several decades of adulthood. Disruption of developmentally dynamic regulatory processes is likely to contribute to neurodevelopmental and neuropsychiatric disorders. In keeping with this expectation, spatiotemporal expression patterns have implicated mid-fetal brain development as a vulnerable process and the prefrontal cortex as a vulnerable region for autism spectrum disorder (ASD) and schizophrenia risk genes. More generally, atypical trajectories of brain maturation have been described in ASD, schizophrenia, and other neuropsychiatric traits and disorders. Given that neuropsychiatric disorders have discrete ages of onset and progression and may arise because of genetic or environmental insults at various times during the life of an individual, there is a clear need to examine gene expression and neuropsychiatric risk across the span of human brain development.

3In addition to spatiotemporal variation, genetic sequence variants also affect gene expression levels, which can contribute to differences in brain structure, function, and behavior. Several laboratories and consortia have systematically identified such expression quantitative trait loci (eQTLs) in numerous tissues, including the brain, but fewer include the developing human brain. Therefore, developmentally regulated eQTLs are sparsely represented in the current catalog of human brain eQTLs, highlighting the need for additional resources. Such eQTL catalogs offer the potential to gain insight into the functional consequences of the hundreds of coding and noncoding genetic loci that have been associated with neuropsychiatric traits and disorders, including developmental delay, ASD, educational attainment, schizophrenia, major depressive disorder, and Alzheimer’s disease.

4To help fill this gap, we generated BrainVar, a unique resource of whole-genome sequencing (WGS) paired with bulk tissue RNA sequencing (RNA-seq) of 176 samples from the human dorsolateral prefrontal cortex (DLPFC) across development, from 6 post-conception weeks to young adulthood (20 years). We focused our analyses on the DLPFC because of its importance in higher-order cognition and the observation that many risk genes for ASD and schizophrenia are co-expressed in the DLPFC during mid-fetal development. We present a systematic description of this resource, including demographics, gene expression across development, gene co-expression modules, and eQTLs. We describe interactions between these factors and comparisons with the BrainSpan dataset, cell-type-specific genes, and loci associated with neuropsychiatric traits and disorders. Our analysis replicates the late-fetal transition observation, a dramatic shift in gene expression between mid-fetal development and infancy; refines the timing of this event; and delineates the degree to which each gene is involved. We also identify 252,629 cis-eQTLs affecting 8,421 genes and classify their effects as prenatal-predominant, postnatal-predominant, or constant across brain development. Finally, we identify eQTLs that co-localize with genome-wide association study (GWAS) loci, linking specific genes to neuropsychiatric phenotypes.

5To characterize gene expression across prenatal and postnatal development of the human DLPFC and to identify genetic variants associated with expression changes, post mortem tissue was obtained from 176 de-identified, clinically unremarkable donors (genotypic sex: 104 male, 72 female) without known neuropsychiatric disorders or large-scale genomic abnormalities, ranging between 6 post-conception weeks and 20 years of age (Figure 1; Table S1). In keeping with prior analyses, we assign these samples to 12 developmental periods, which we group into four developmental epochs (Figures 1 and 2). Gene expression data were generated using RNA-seq from tissue dissected from the DLPFC (corresponding mainly to Brodmann area 46) or from the frontal cerebral wall (donors younger than 10 post-conception weeks). WGS data (31.5× median coverage) were generated simultaneously from DNA isolated from the same individuals.

6RNA-seq reads were aligned and converted to log base 2 counts per million (log2CPM) per gene (STAR Methods), with 23,782 genes meeting minimum expression criteria. We restricted further analysis to these 23,782 cortically expressed genes, of which 16,296 (68.5%) encode proteins, whereas 7,486 (31.5%) are noncoding, including long noncoding RNA (lncRNA) (12.6% of total) and antisense (9.2% of total) genes (Table S2). For the 14 samples also profiled in BrainSpan, gene expression was highly correlated per sample and per gene (Figure S1). In both datasets, the first principal component of gene expression is strongly correlated with developmental age (Figures 2A and S1). All samples were genotypically concordant between the WGS and RNA-seq data. Ancestry correlated strongly between principal-component analysis clusters and self-report (STAR Methods).

7Prior analysis of the 40 brains in the BrainSpan cohort identified developmental age as the greatest source of between-sample variance in gene expression, especially during a “late-fetal transition” between 22 post-conception weeks and 6 postnatal months. We replicate these findings in BrainVar. The first principal component explains 42% of the variance in gene expression and is highly correlated with developmental age (partial R2 = 0.88; Figures 2A and S2; similar results when excluding the 14 overlapping samples), with the greatest changes occurring in late fetal development and early infancy (Figure 2A).

8Using the increased resolution from the 176 brains in BrainVar, we show that the late-fetal transition begins around 19 post-conception weeks (start of period 6) and that the most dramatic changes are complete by 6 postnatal months (end of period 8); we label this transitional phase as epoch 2 (Figure 2A). Considering the nine samples younger than 10 post-conception weeks (periods 1–2), we also observe an “early-fetal transition,” i.e., a coordinated shift in embryonic and early fetal development, which we label epoch 0 (Figure 2A).

9To identify the specific genes that change in the late-fetal transition, we performed a trajectory analysis on the 167 samples in epochs 1–3; we excluded epoch 0 because of the sparse sampling before and during the early-fetal transition. Remarkably, over half of the genes expressed in the cortex exhibit a persistent, progressive, and statistically significant expression variance across this late-fetal transition (Figure 2B). We identified three distinct trajectories, with 6,934 “rising” genes (higher postnatal expression), 5,143 “falling” genes (higher prenatal expression), and 11,705 “non-transitional” genes (no statistically significant change). Considering more than three trajectories further split gene sets by variance rather than developmental profile (STAR Methods). Similar trajectories are observed for these three gene lists in the BrainSpan DLPFC data (Figure S2). Emphasizing the magnitude of this transition, the first principal component of the 11,705 non-transitional genes explains only 18.4% of the variance in gene expression and is weakly correlated with developmental age (partial R2 = 0.3; Figure S2).

10The magnitude of the changes in individual genes’ expression levels across late-fetal transition can be estimated by calculating the difference in log2CPM expression between epoch 3 and epoch 1, ranging from 12.1 (OPALIN, a component of myelin) to −9.2 (IGF2BP1, an IGF2 binding protein); for context, the median epoch 3-to-epoch 1 changes in log2CPM values were 2.1, 0.1, and −1.1 for rising, non-transitional, and falling, respectively. The majority of changes in gene expression reflected relative amplification or attenuation of expression levels rather than binary presence/absence of expression, with only 621 rising genes and 95 falling genes specific to epoch 3 or 1 (defined as log2CPM ≤ −5 in the other epoch; Figure 2C; Table S2).

11Compared with rising and non-transitional genes, falling genes had the highest median expression in epoch 1 (p < 2 × 10−16) and epoch 3 (p < 2 × 10−16, Wilcoxon rank-sum test [WRST]; Figure 2C) and the highest fraction of protein-coding genes (p < 2 × 10−16, Fisher’s exact test [FET]; Figure 2D) and were highly enriched for genes with high probability loss-of-function intolerant (pLI) scores (p = 5 × 10−11, WRST; Figure 2E). High pLI scores reflect detection of fewer protein-truncating variants than expected, suggesting that loss-of-function mutations in the gene are disfavored by natural selection (i.e., the gene is haploinsufficient). Rising genes had a similar proportion of protein-coding genes as falling genes (p = 0.94, FET; Figure 2D) but were depleted for genes with high pLI scores (p = 1 × 10−7, WRST; Figure 2E). If the timing of a gene’s highest expression corresponds to the timing of its most critical functions, then the pLI difference between falling and rising genes suggests that prenatal development is especially sensitive to haploinsufficiency.

12Compared with RNA-seq data from 53 adult tissues, falling genes were only enriched in non-cortical tissues (driven by genes related to RNA transcription and cell division; Table S2), whereas rising genes were enriched for many brain regions, including the adult cortex and excluding cerebellum (Figure 2F), highlighting the distinctions between the fetal and adult cortex. Non-transitional genes had the lowest proportion of protein-coding genes and were expressed ubiquitously across adult tissues (Figure 2).

13To capture the contribution of changing cell type proportions to gene expression profiles, we assessed expression trajectories of genes specific to each of ten cortical cell types from prenatal and postnatal human brain. The estimated profiles of all ten cell types vary dramatically across epoch 2, with radial glia/neural progenitor cells and fetal neurons decreasing as mature neurons and other glial cells increase (Figures 2G and 2H); this pattern is replicated in BrainSpan DLPFC samples (Figure S2). These analyses support the hypothesis that varying cell type proportions are major contributors to the late-fetal transition in the DLPFC, but distinguishing cellular composition effects from differential expression within a cell type will require single-cell data from across this age range.

14Tofurthercharacterizetherelationshipsbetweenthe23,782cortically expressed genes, we applied a weighted gene co-expression network analysis (WGCNA) to define 19 consensus modules that included 10,459 genes (Figures 3A and S3; Table S3). As expected, genes within each module shared functional roles (Figure 3B), temporal trajectories of geneexpression(Figures3C and3D),regulatory transcriptionfactors (Figure S3), and cell type enrichment (Figures 3E and 3F). Module preservation analysis using BrainSpan data identified similar co-expression patterns across brain regions, especially independent DLPFC samples (Figure 3G).

15Similar to transitional genes (Figure 2), multidimensional scaling of the module eigengenes demonstrated that developmental age accounted for 44.7% and 36.1% of the variance in the first two dimensions. Considering the position of the 19 modules along these two dimensions and the developmental trajectories of the genes in each module, we identified five groups of related modules (Figures 3A and 3B). Group 1 modules (M1 black, M2 royal blue, M3 green-yellow, and M4 yellow) are enriched for falling genes, whereas group 5 modules (M16 blue, M17 silver, M18 light cyan, and M19 turquoise) are enriched for rising genes (Figure 3D). The remaining three groups (2, 3, and 4) are enriched for non-transitional genes (Figure 3D).

16Five modules are of particular note. The M2 royal blue module (group 1) captures cell cycle ontology and is enriched in neuroprogenitor cells, radial glia, and intermediate progenitor cells. The M4 yellow module (group 1) is enriched for numerous ontology terms related to neuronal development, contains genes specific to neuronal stem cells (e.g., NCAM1/ncam and PROM1/cd133), and is highly enriched for genes related to maturing excitatory and inhibitory neurons. The M8 red module (group 3) is enriched for ontology terms relating to cell fate and morphogenesis, is highly enriched for noncoding genes, and has an expression peak in early fetal development, capturing many genes that are involved in early-fetal transition (Figure 2A). Several genes associated with regional patterning in non-cortical tissues, including the hindbrain (e.g., UNCX and CCDC140) and hypothalamus (e.g., DMBX1 and SOX14), are expressed at high levels in this module. The M18 light cyan and M19 turquoise modules (group 5) are strongly enriched in glial and other nonneuronal cell clusters; accordingly, both modules are enriched for ontology terms related to immune responses. The M19 turquoise module is also enriched in excitatory neurons in the postnatal cortex and ontology terms relating to synaptic signaling and neurotransmitter transport.

17We next considered the intersection between genes associated with developmental trajectories, modules, or cell types and genes associated with ten human traits and disorders. For ASD and developmental delay with and without seizures, we used gene lists derived from exome association studies of rare and de novo variants. For educational attainment, attention deficit hyperactivity disorder (ADHD), schizophrenia, major depressive disorder, multiple sclerosis, Parkinson’s disease, and Alzheimer’s disease, we used genes within 10 kb of the lead SNP detected in GWASs. For our analyses, we excluded genes within the major histocompatibility complex on chromosome 6 because of the complicated nature of this region.

18Developmental delay, ASD, and educational attainment genes were enriched for falling genes (p = 8.5 × 10−6, p = 5.1 × 10−3, and p = 2.4 × 10−4 respectively; FET adjusted for 30 comparisons), consistent with a prenatal origin for aspects of their neurobiology. A non-significant trend toward enrichment for rising genes was observed for Parkinson’s disease and Alzheimer’s disease (Figure 4A; Table S4). The M4 yellow module was enriched for ASD and educational attainment genes, including NRXN1, TCF4, and BCL11A (Figure 4B; Table S4), and the M9 brown module (enriched for chromatin organization Gene Ontology terms and non-transitional genes; Figure 3) was enriched for genes associated with developmental delay and educational attainment, including CDK13, PACS1, and EP300 (Figure 4B; Table S4).

19Across the ten CNS traits and disorders, five cell type clusters (C) showed significant enrichment in fetal brain (Figure 4C; Table S4) and none in adult brain (Figure S4; Table S4). ASD genes were enriched for C18 excitatory newborn neurons and C1 striatal interneurons (Figure 4C), in keeping with a role of excitatory and inhibitory lineages. Both lineages were also enriched in educational attainment, specifically C3 early excitatory neurons and C6 medial ganglionic eminence (MGE)-derived interneurons, whereas genes associated with developmental delay with seizures were enriched in C15 caudal ganglionic eminence (CGE)-derived interneurons (Figure 4C). We observed a nominally significant trend toward enrichment of C19 microglia genes in multiple sclerosis and Alzheimer’s disease.

20We identified 6,573,196 high-quality SNPs and insertions or deletions (indels) from the WGS data using methods described previously, with an allele frequency of at least 5% in our prenatal (periods 1–6, n = 112) and postnatal (periods 8–12, n = 60) samples (Figure 1). To identify eQTLs within 1 Mb of a gene (eGene), we used linear regression for adjusted expression level (STAR Methods), with developmental period, sex, and the first five principal components for ancestry as covariates. Results were corrected for multiple comparisons using Benjamini-Hochberg (false discovery rate [FDR] ≤ 0.05). To distinguish temporal-predominant eQTLs, we performed three cis-eQTL analyses: all 176 samples (complete sample, 216,026 eQTLs of 5,728 eGenes), 112 prenatal samples (periods 1–6, 154,440 eQTLs of 4,378 eGenes), and 60 postnatal samples (periods 8–12, 51,528 eQTLs of 2,199 eGenes). These discovery rates are in line with similarly sized cohorts (Figure S5). The union of these three analyses identified 252,629 eQTLs of 8,421 eGenes (Table S5). As expected, the eQTLs are enriched for markers of active transcription derived from the human brain (Figure S5;;; ). We find that eQTL effect size and direction are correlated with prenatal whole brain (Pearson’s r = 0.73, p ≤ 1 × 10−16; Figure 5A) and postnatal (adult) frontal cortex (Pearson’s r = 0.73, p ≤ 1 × 10−16; Figure 5A) from independent datasets.

21We leveraged the consistently processed prenatal and postnatal data in BrainVar to identify eQTLs with differing effect sizes across development (Figure 5B; STAR Methods). The majority of eQTLs were constant, reaching nominal significance in all three analyses with the same direction of effect (161,923 eQTLs, 64.1% of the total). Many eQTLs were prenatal-predominant, with significantly greater prenatal than postnatal effect sizes (24,760 eQTLs, 9.8%). Fewer eQTLs were postnatal-predominant, with significantly greater postnatal than prenatal effect sizes (9,352 eQTLs, 3.7%). The remaining 56,593 eQTLs (22.4%) showed a trend toward stronger prenatal effects (19.8%) or postnatal effects (2.6%). With larger sample sizes, we would expect a greater fraction of constant eQTLs to show some degree of temporal specificity, especially postnatal. Although the magnitude of effect varied across development for many eQTLs, we did not observe a single eQTL with opposing prenatal and postnatal directions of effect.

22Most eGenes have more than one eQTL (5,538 of the 8,421, 65.8%). Defining the top eQTL per eGene as that with the lowest FDR-significant p value in any of the three sample sets (Table S5), we identified 2,977 (35.4% of total) constant eGenes, 1,691 (20.1%) prenatal-predominant eGenes, and 1,145 (13.6%) postnatal-predominant eGenes (Figure 5B). The remaining 2,608 eGenes (31.0%) trend toward prenatal (25.1%) or postnatal (5.9%) effects. Because of linkage disequilibrium (LD), cis-eQTLs for an eGene are likely to have a similar direction and magnitude of effect; accordingly, the temporal category of the top eQTL matched the majority of eQTLs for 88.2% of all eGenes (7,425 eGenes; Figure S5).

23To validate the prenatal- and postnatal-predominant eQTLs, we evaluated their performance in independent datasets. In prenatal whole brain, we observed stronger correlation for the effects of prenatal-predominant (r = 0.54, p = 3.9 × 10−22) than postnatal-predominant eGenes (r = 0.36, p = 1.1 × 10−3; Figure 5C). In contrast, in postnatal frontal cortex, stronger correlations were observed for the effects of postnatal-predominant (r = 0.60, p = 4.0 × 10−24) than prenatal-predominant eGenes (r = 0.37, p = 1.1 × 10−6;Figure 5C).

24The top eQTL for constant eGenes is closer to the transcription start site than that for temporal-predominant eGenes (median: 92,223 bp versus 403,702 bp, p ≤ 2 × 10−16, two-sided WRST; Figure 5D). However, −log10(P) values increased with proximity to the transcription start site for constant and temporal-predominant eQTLs (for both sets, p ≤ 2 × 10−16, linear regression), as did eQTL effect size to a small degree (constant p = 1.0 × 10−12, temporal-predominant p = 0.03, linear regression). Compared with constant eGenes, temporal-predominant eGenes also included a higher proportion of protein-coding genes (odds ratio [OR] = 1.56, 95% confidence interval [CI]: 1.40–1.74), p = 3.7 × 10−16, two-sided FET; Figure 5E), genes with high pLI scores (pLI ≥ 0.995; OR = 1.45 [95% CI: 1.12– 1.90], p = 0.004, two-sided FET; Figure 5F), and greater connectivity in protein-protein interaction (PPI) networks (median Z score 1.45 versus 0.91, p ≤ 2 × 10−16, two-sided WRST;Figure 5H).

25Given the dynamic expression profiles over development (Figure 2), we expected prenatal-predominant eGenes to be enriched for falling genes and postnatal-predominant eGenes to be enriched for rising genes. We did not observe these effects (Figure 5I). Instead, the prenatal-predominant eGenes are enriched for rising genes (OR = 1.5 [95% CI: 1.3–1.6], p = 1.1 × 10−10, two-sided FET; for example, Figure 5J), and pre- and postnatal-predominant eGenes are depleted for falling genes (prenatal OR = 0.79 [95% CI: 0.69–0.9], p = 4.2 × 10−4; postnatal OR = 0.85 [0.73–0.997], p = 0.04; two-sided FET). Instead, we observed coordination between the timing of eQTLs’ strongest effects and the timing of eGenes’ greatest expression variation between samples (STAR Methods). Prenatal-predominant eGenes are strongly enriched for genes with greater prenatal variance (OR = 2.3, p = 4.0 × 10−52, two-sided FET) and depleted for genes with greater postnatal variance (OR = 0.36, p = 1.6 × 10−15, two-sided FET). Postnatal-predominant eGenes show a complementary but weaker pattern of enrichment for genes with greater postnatal variance (OR = 1.1, p = 0.25, two-sided FET) and depletion of genes with greater prenatal variance (OR = 0.76, p = 7.0 × 10−5, two-sided FET). Considering the role of selective pressure, we observed that genes with higher pLI scores also had lower eQTL effect sizes (p = 4.0 × 10−36; two-sided WRST; Figure 5J) as well as lower expression variance between samples (Figure 5K).

26The differences in constant and temporal-predominant eGenes led us to consider how genes associated with neuropsychiatric traits and disorders (Figure 4) fit into this classification. Gene sets associated with traits by GWAS loci or exome sequence association followed the patterns of temporal-predominant eGenes but to a greater extent, with a higher proportion of protein-coding genes (Figure 5E), higher pLI scores (Figure 5F), and stronger clustering within PPI networks (Figure 5H).

27At the variant level, we expect GWAS loci to be enriched for eQTLs in relevant tissues. Using a permutation-based method accounting for LD structure, minor allele frequency (MAF), and gene density (STAR Methods), we tested four of the larger GWASs and observed eQTL enrichment for educational attainment, schizophrenia, and multiple sclerosis but not Alzheimer’s disease (Figure S6; Table S6). We did not see evidence of the reverse hypothesis that eGenes are enriched for GWAS signals (Figure S6). Using a colocalization analysis, we looked for overlap between specific eQTL loci with educational attainment and schizophrenia GWAS loci using a posterior probability of colocalization threshold of 0.8. In the schizophrenia GWAS, 13 of 108 loci (12.0%) showed evidence of colocalization, including two prenatal-predominant and two postnatal-predominant eQTLs (Table S6). A lower proportion of educational attainment loci showed evidence of colocalization (4.1%, 52 of 1,271), including 14 prenatal-predominant and two postnatal-predominant eQTLs (Table S6). Focusing on multigenic loci with the strongest evidence of colocalization, we implicate specific genes and expression changes as the likely mechanism underlying the GWAS loci. SNPs associated with educational attainment at a chromosome 14 locus (Figure 6A) colocalized only with eQTLs for the lncRNA LOC101926933 (also called RP11–298I3.1, AL132780.1, or ENSG00000257285; Figure 6B). Across the locus, the p values for the GWAS SNPs and LOC101926933 eQTLs are highly correlated, resulting in a posterior probability for colocalization of 0.92 (Figure 6C); for this locus, eQTL effect size is similar across development (Figure 6D). We also observe colocalization of prenatal-predominant eQTLs with educational attainment or schizophrenia GWAS loci. For example, SNPs contributing to an educational attainment GWAS locus on chromosome 12 (Figure 6E) overlap specifically with eQTLs for the protein-coding gene RHEBL1, which encodes a brain-enriched G-protein activator of the mTOR pathway (Figure 6F). GWAS and RHEBL1 eQTL prenatal p values are highly correlated and result in a posterior probability for colocalization of 0.97 (Figure 6G). We see a significantly greater eQTL effect size in prenatal compared with postnatal samples (p = 0.003; Figure 6H), with higher RHEBL1 expression associated with increased educational attainment.

28In this manuscript, we describe BrainVar, a unique resource of paired genome (WGS) and transcriptome (bulk tissue RNAseq) data derived from 176 human DLPFC samples across prenatal and postnatal development (Figure 1). We identify 23,782 genes expressed during human cortical development, gene lists relating to developmental trajectories and co-expression, and common variants that alter gene expression (eQTLs). Our analyses show how these datasets relate to each other and to gene expression in cell types derived from single-cell RNA-seq data and to CNS traits and disorders derived from genomic analyses (exome sequencing and GWAS). In addition to developing a resource with utility for future studies of human development, neurobiology, and neuropsychiatric disorders, we also describe key biological insights, including the nature of the late-fetal transition in gene expression (Figures 2 and 3), an early-fetal transition (Figures 2 and 3), developmental processes and cell types implicated in CNS traits and disorders (Figure 4), eQTLs split by effect size across development (Figure 5), differing characteristics of genes with constant versus temporal-predominant eQTLs (Figure 5), and the application of this dataset to implicate specific genes at GWAS loci (Figure 6).

29Principal component analysis identifies developmental age as the most important factor underlying the variance in gene expression in this dataset. The majority of this temporal variance occurs in two transitional phases (Figure 2), the early-fetal and late-fetal transitions. The early-fetal transition is a coordinated decrease in expression of multiple genes in early development (epoch 0; periods 1–2; 6–10 post-conception weeks) that coincides with the establishment of regional identity across the brain. Concordant with a possible role of the early-fetal transition in this process, the expression of several genes associated with non-cortical tissues (e.g., UNCX and DMBX1) is decreased during this period. In addition, we found that the early-fetal transition is captured in the M8 red module (Figure 3), which is enriched for lncRNA transcripts and Gene Ontology terms related to morphogenesis and cell fate.

30The late-fetal transition between mid-fetal development and infancy involves over 12,000 genes with similar numbers rising and falling (Figure 2). Prior reports of humans and primates associated this transition with a reduction in intra- and inter-regional variation evident at the levels of bulk tissue and individual cell types. Our data similarly suggest that this transition represents a combination of changes in the relative proportions of various cell types and biological processes within these cells (Figures 2 and 3). Critically, the larger BrainVar sample set allowed us to define 19 post-conception weeks as the inflection point at which the late-fetal transition begins (Figure 2), further distinguishing the late-fetal transition from previously reported organotypic changes.

31Although previous analyses have identified eQTLs in human brain tissue postnatally and prenatally, no prior study has assessed the effect of genomic variation on gene expression across the whole of brain development, from embryogenesis through fetal development, infancy, childhood, and adolescence and into young adulthood. Consequently, we were able to identify temporal-predominant eQTLs that have a greater effect on expression prenatally or postnatally (Figure 5B). The eQTLs identified here were highly correlated with prior eQTL catalogs despite differing cohorts, methods, and analysis (Figures 5A and 5B). Furthermore, comparison with these independent catalogs support our temporal categorization of eQTLs, with prenatal-predominant eQTL effects more correlated in prenatal whole-brain and postnatal-predominant eQTL effects more correlated in the postnatal frontal cortex (Figure 5C).

32Across multiple metrics, we observe dramatic differences between eGenes with constant and temporal-predominant eQTLs (Figures 5D–5G). Compared with other genes expressed in the cortex, genes affected by constant eQTLs are more likely to be noncoding and have low pLI scores and few protein-protein interactions. In contrast, genes regulated by eQTLs with a degree of temporal specificity are similar to genes for which we did not detect eQTLs. Critically, we find that pLI score, a measure of sensitivity to variation in genetic sequence, and eQTL effect size are inversely related (Figure 5J). Furthermore, prenatal-predominant eGenes are more common among rising genes, which have their highest expression during postnatal time. These observations suggest that developmental and evolutionary constraints limit the population frequency or effect of eQTLs on key developmental processes, a hypothesis that might be testable in future studies as additional information concerning spatiotemporal and cell type specificity of enhancers and eQTLs becomes available for a variety of tissues. Under this model, constant eQTLs with high effect sizes tend to influence genes that tolerate variation in expression (e.g., non-rate-limiting metabolic steps) or are non-critical to brain function, whereas temporal-predominant eQTLs tend to influence genes with critical roles that are sensitive to variation in genetic sequence but only to a small degree or at a stage in development when variation in expression of the gene is tolerated.

33The eQTLs identified here also provide insights into CNS traits and disorders, with co-localization in 13 of 108 GWAS loci for schizophrenia and 52 of 1,271 GWAS loci for educational attainment, including the lncRNA LOC101926933 and the protein-coding gene RHEBL1 (Figure 6). LOC101926933 remains largely uncharacterized, whereas RHEBL1 (Ras homolog enriched in brain-like 1) is a highly conserved G-protein that activates mTOR, a pathway that has been implicated previously in neurodevelopmental and neurodegenerative disorders. Our results suggest that higher expression of RHEBL1, which may lead to greater mTOR activation, is associated with increased educational attainment. Of note, RHEBL1 has a pLI score of 0, suggesting that loss of one allele does not lead to a selective disadvantage. Higher-resolution datasets across development, including single cells, additional brain regions, and larger sample sizes, along with complementary analyses of brains of individuals with neuropsychiatric disorders and rare genetic disorders, are likely to provide additional insights. The combination of genomic and transcriptomic data across development allows us to interrogate human cortical development from a molecular perspective at a higher resolution than before. Understanding patterns of temporal and cell type specificity, along with eQTL colocalization to resolve GWAS loci, has already provided insights into the pathology underlying neuropsychiatric disorders. Further delineation of these patterns is likely to be critical for a detailed understanding of etiology as a foundation for therapeutic development.

34Further information and requests for resources and reagents should be directed to and will be fulfilled by Stephan Sanders (stephan.sanders@ucsf.edu).

35This study did not generate new unique reagents.

36This study was conducted using postmortem human brain specimens from tissue collections at the Department of Neuroscience at Yale University School of Medicine. Additional specimens were procured from the Birth Defects Research Laboratory at the University of Washington, Advanced Bioscience Resources Inc., Human Brain Collection Core (HBCC), the Brain and Tissue Bank at the University of Maryland, the MRC-Wellcome Trust Human Developmental Biology Resource at the Institute of Human Genetics, University of Newcastle, UK, and the Human Fetal Tissue Repository at the Albert Einstein College of Medicine (AECOM). Tissue was collected after obtaining parental or next of kin consent and with approval by the institutional review boards at the Yale University School of Medicine, the National Institutes of Health, and at each institution from which tissue specimens were obtained. Tissue was handled in accordance with ethical guidelines and regulations for the research use of human brain tissue set forth by the NIH (https://oir.nih.gov/sites/default/files/uploads/sourcebook/documents/ethical_conduct/guidelines-biospecimen.pdf and the WMA Declaration of Helsinki (https://www.wma.net/policies-post/wma-declaration-of-helsinki-ethical-principles-for-medical-researchinvolving-human-subjects/).

37All available non-identifying information was recorded for each specimen in Table S1. In total, 176 postmortem brain specimens (104 male, 72 female; postmortem interval of 21.7 ± 15.9 (mean ± SD) hours and pH, 6.41 ± 0.35) ranging in age from 6 post-conception weeks to 20 postnatal years (Figure 1; Table S1) were included in this study. Fetal age was extrapolated based on the date of the mother’s last menstruation, characteristics of the fetus noted upon ultrasonography scanning, foot length of the fetus, and visual inspection. The postmortem interval (PMI) was defined as hours between time of death and time when tissue samples were frozen.

38Tissue was dissected as described previously. Samples collected from 6 – 9 post-conception weeks specimens contained the entire thickness of the cerebral wall. Samples collected from 12 to 22 post-conception weeks specimens contained the cortical plate. Samples from 35 post-conception weeks to 20 postnatal years were dissected such that the entire gray matter (layer 1–6) and part of the underlying subplate (4 – 12 postnatal months) or white matter (1 – 20 postnatal years) were collected.

39Total RNA was extracted using mirVana kit (Ambion) with some modificationsto the manufacturer’s protocol, as described below. Each tissue sample was pulverized with liquid nitrogen in a prechilled mortar and pestle and transferred to a chilled safe-lock microcentrifuge tube (Eppendorf). Per tissue mass, equal mass of chilled stainless-steel beads (Next Advance, catalog # SSB14B) along with one volume of lysis/binding buffer were added. Tissue was homogenized for 1 min in Bullet Blender (Next Advance) and incubated at 37°C for 1 min. Another nine volumes of the lysis/binding buffer were added, homogenized for 1 min, and incubated at 37°C for 2 min. One-tenth volume of miRNA Homogenate Additive was added, and extraction was carried out according to the manufacturer’s protocol. RNA was treated with DNase using TURBO DNA-free Kit (Ambion/Life Technologies) and RNA integrity was measured using Agilent 2200 TapeStation System.

40Barcoded libraries for RNA-seq were prepared with 5ng of RNA using TruSeq Stranded Total RNA HT Sample Prep Kit with Ribo-Zero Gold kit (Illumina) per manufacturer’s protocol. Paired-end sequencing (100 bp × 2) was performed on HiSeq 4000 sequencers (Illumina) at Yale Center for Genome Analysis.

41Genomic DNA was isolated using the QIAamp DNA Mini Kit (QIAGEN). In detail, approximately 25 mg of brain tissue was transferred to a chilled safe-lock microcentrifuge tube (Eppendorf) and equal mass of chilled stainless-steel beads (Next Advance, catalog # SSB14B) along with 90 μl of buffer ATL were added. Tissue was homogenized for 1 min in Bullet Blender (Next Advance) and incubated at 37°C for 1 min. Another 90 ul of buffer ATL was added and blended for an additional minute. After incubation on ice for 5 min, tubes were gently centrifuged to collect beads at the bottom. Supernatant was transferred to a new tube and 20 μl of Proteinase K was added. Sample was incubated at 56°C for 3 hours in a shaking heat block. After incubation, genomic DNA was further purified following the manufacturer’s protocol. DNA was eluted in nuclease free water and concentration was estimated by nanodrop.

42DNA library preparation and sequencing were carried out at GENEWIZ (New Jersey). Before library preparation, the concentration of the DNA was measured using a fluorescent assay and DNA quality was assessed by visualization on agarose gels. PCR-free DNA library preparation was performed and resulting libraries were sequenced at 2×150 bp to achieve mean coverage of 30x (Table S1).

43Using the pipeline from the Centers for the Common Disease Genomics project, FASTQ reads were aligned to the GRCh38 reference from the 1000 Genomes Project using BWA-MEM version 0.7.15. Reads were sorted and duplicates were removed with Picard, version 2.17.5 (https://github.com/broadinstitute/picard/); base quality score recalibration was then performed with the Genome Analysis Toolkit (GATK), v3.8–0-ge9d806836. Variant calling and joint genotyping were done with the Haplotyper and Genotyper tools from Sentieon v201711.01, a toolkit containing modules that are mathematically equivalent to their counterparts in the GATK. SNP and indel recalibration were performed on the joint genotyped VCF file. Variant Quality Score Recalibration (VQSR) metrics were created from a training set of highly validated variant resources: dbSNP build 138, HapMap 3.3, 1000 Genomes OMNI 2.5, and 1000 Genomes Phase 1. For the following analyses, we excluded: variant calls with any VQSR tranches (keeping “PASS” only), variants located in low-complexity regions, variants located on non-canonical chromosomes (decoy chromosomes or contigs), indels > 50 bp, single nucleotide variants (SNVs) with allele balance > 0.78 or < 0.22 (indels > 0.8 or < 0.2), variants with < 90% call rate, and variants and genotypes that did not meet high quality thresholds as identified in an ROC-based optimization procedure using family-based WGS data.

44Variants with a minor allele frequency of ≥ 5% in both the prenatal (periods 1–6; N = 112) and postnatal (periods 8–12; N = 60) samples and Hardy Weinberg equilibrium p value ≥ 1×10−12, were included in downstream expression quantitative trait locus (eQTL) analysis (N = 6,573,196 variants). For annotation and subsequent analyses, we converted the final VCF into Variant Dataset Format using Hail version 0.1. SNPs and insertions/deletions (up to 50bp) annotation based on the GENCODE comprehensive version 21 using Ensemble VEP version 90.

45RNA-seq reads were aligned to the human genome (hg38/GRCh38) using STAR aligner and gene-level read counts were calculated using HTSeq based on GENCODE v21 annotation.

46Read counts per gene were then converted to counts per million (CPM), which were logarithmically scaled to base 2 (log2CPM). Of the 60,155 genes assessed, 23,782 were defined as being cortically expressed, based on CPM ≥ 1 in at least 50% of samples of either sex in at least one of the 12 developmental periods (Table S2). For these 23,782 genes, the median log2CPM ranged from −5.9 to 12.3, with a median of 2.3.

47The read count data matrix of 176 samples by 60,155 genes (Table S2) was normalized as follows: Step 1: Read count data matrix converted to counts per million (CPM).Step 2: Genes with CPM ≥ 1 CPM in at least 50% of the samples in any one sex in any one period were included; 23,782 genes passed these criteria.Step 3: CPM values were transformed to log2(CPM) using the voom function in the limma R package.Step 4A: In order to correct for the technical artifacts, we performed hidden covariate analysis on the residuals of the expression matrix after developmental period and sex were subtracted from the log2(CPM) data matrix using the hidden covariate analysis method (HCP).Step 4B: In parallel, we performed surrogate variable analysis (SVA) on residuals of the expression matrix after developmental period and sex were subtracted from the log2(CPM) data matrix using the SVA R package.Step 5: In each of Step 4A and 4B, we subtracted contributions from 20 hidden covariates (HCP) and 2 surrogate variables (SVA) from the log2(CPM) data matrix.

48Step 1: Read count data matrix converted to counts per million (CPM).

49Step 2: Genes with CPM ≥ 1 CPM in at least 50% of the samples in any one sex in any one period were included; 23,782 genes passed these criteria.

50Step 3: CPM values were transformed to log2(CPM) using the voom function in the limma R package.

51Step 4A: In order to correct for the technical artifacts, we performed hidden covariate analysis on the residuals of the expression matrix after developmental period and sex were subtracted from the log2(CPM) data matrix using the hidden covariate analysis method (HCP).

52Step 4B: In parallel, we performed surrogate variable analysis (SVA) on residuals of the expression matrix after developmental period and sex were subtracted from the log2(CPM) data matrix using the SVA R package.

53Step 5: In each of Step 4A and 4B, we subtracted contributions from 20 hidden covariates (HCP) and 2 surrogate variables (SVA) from the log2(CPM) data matrix.

54Unadjusted log2CPM gene expression data for 176 samples by 23,782 genes was used for expression trajectory (Figure 2) and WGCNA analyses (Figure 3), while the adjusted (HCP/SVA) values were used for eQTL analysis.

55To confirm that the WGS and RNA-seq data from each sample were of sufficient quality for downstream analysis and corresponded to the same individual, a series of quality metrics and checks were performed (Table S1).

56For the WGS data, coverage metrics were assessed using PicardTools (v2.18.1). Mean coverage per sample ranged from 22.7–65.4x, with a cohort median of 31.5x. Across all samples, 92.1%–93.6% of the mapped genome was covered at 10x or greater (mean of 92.8%). The FREEMIX metric from VerifyBamId (v1.1.3; ) was used to identify samples with potential contamination, with a maximum observed FREEMIX score of 0.064, suggesting no contamination.

57Sample identity was verified by comparing sex and genotype between the WGS and RNA-seq data. In the WGS data, sex was determined from chromosome X heterozygosity using Peddy (v0.3.2;), with the Peddy hg19.sites converted to GRCh38 using the UCSC Genome Browser LiftOver utility. High-quality variants with an allele frequency ≥ 1% were exported from the VCF using Hail for input into Peddy. In the RNA-seq data, sex was determined from the expression levels of XIST and the 18 most highly expressed genes on chromosome Y: KDM5D, DDX3Y, ZFY, TBL1Y, PCDH11Y, PRKY, USP9Y, RPS4Y1, TXLNGY, NLGN4Y, TTTY14, UTY, EIF1AY, GYG2P1, TTTY10, TTTY15, KALP. Based on gene-specific expression thresholds determined by visual inspection of bi-modal expression histograms, each sample’s sex was predicted according to the expression level of all 19 genes. Sex was consistent in the WGS and RNA-seq data for all 176 samples and matched the recorded sex in 132 out of 134 samples with such data (55/56 females, 77/78 males).

58To confirm identity by genotype, we compared the genotypes from 289 common, coding SNPs with high fixation index (FST), called from both the WGS and RNA-seq data. Genotypes were callable in both data types for 118–206 SNPs per sample (40.8%–71.3% of 289 SNPs; median = 177, 61.2%). SNP variant genotypes were highly concordant between the WGS and RNA-seq data (median 87.4% concordance between WGS and RNA-seq for the corresponding sample; lowest concordance 73%), with corresponding samples showing higher concordance than comparisons between all discordant samples. There was no evidence of duplicate or closely related samples (SNP-based relatedness coefficients from Peddy: −0.000332 to 0.1481).

59To confirm the approximate accuracy of samples’ reported age, the expression level of the doublecortin gene (DCX) was examined. DCX is involved in neuron migration, and is expressed most strongly during prenatal development, with distinctly decreased postnatal expression. All 176 samples showed the expected DCX expression levels given samples’ reported age. Similar results across all expressed genes were obtained using principal component analysis (below).

60Ancestry was estimated using principal component analysis of common SNPs and indels in the WGS data, run alongside 3,804 additional individuals of known ancestry with WGS data (parents from the Simons Simplex Collection; ). From 10,688,106 SNPs and indels with allele frequency ≥ 5% in either this dataset, one of the three batches of Simons Simplex Collection data, or GnomAD genomes, variants were pruned for independence with linkage disequilibrium r2 < 0.1 and then randomly downsampled to 118,849 variants. Principal component analysis was run using Hail 0.1. The first two principal components were used to classify samples by ancestry, and the first five principal components were used as covariates in the identification of eQTL loci (Figure S1).

61Using the first two principal components of SNP-based ancestry estimated above, we identified clusters corresponding to African American ancestry (42 samples, 24% of cohort), European ancestry (82 sample, 47%), and Asian ancestry (4 samples, 2%), while 48 samples (27%) were outside of these clusters and enriched for individuals who identified as Hispanic, Alaskan native, or mixed ancestry (Figure S1; Table S1). Analyzing the first principal component of ancestry shows ancestry groups were differentially represented across developmental periods (F = 3.4, df = 5, p = 0.006, ANOVA), with post hoc analysis showing a greater representation of African American samples in later developmental periods compared with samples of Hispanic, Alaskan native, or mixed ancestry (p.adj = 0.02, TukeyHSD).

62PCA was performed on the covariance matrix of 23,782 cortically expressed genes in 176 samples. A secondary PCA was performed on the 11,705 Non-transitional genes in 167 samples (excluding Period 1 and 2 samples) to assess the extent to which removing the late-fetal transition accounted for variance in gene expression. For each PCA, the variance explained by each principal component was assessed (Figure S2).

63To quantify the relative contributions of biological and technical covariates, we calculated the partial R2 of each covariate with each principal component using the rsq R package, in a generalized linear model where loadings of principal components are considered as a response and biological and technical covariates (such as developmental period, sex, sequencing batch, sequencing depth, RNA integrity number (RIN), mitochondrial RNA proportion, ribosomal RNA proportion, intronic reads proportion, intergenic reads proportion) are considered as predictor variables (Figure S2).

64While the BrainVar bulk tissue RNA-seq dataset catalogs gene expression in the DLPFC of 176 donors across development (6 post-conception weeks to 20 years of age), BrainSpan catalogs gene expression across 16 brain regions, including DLPFC, in 40 brains ranging from 8 post-conception weeks to 40 postnatal years. Of note, data were generated for 14 brains in both BrainVar and BrainSpan. To compare the data from the 14 samples profiled in both BrainSpan and BrainVar, we reprocessed the entire BrainSpan dataset using the BrainVar RNA-seq analysis pipeline, transforming raw read counts to log2CPM values. We then filtered to genes that passed BrainVar minimum expression criteria and calculated sample-to-sample and geneto-gene Pearson correlation coefficients. We performed principal component analysis on a covariance matrix of the 14 commonly sequenced samples and the 23,782 BrainVar-expressed genes for both datasets. Despite differences in library preparation, with BrainSpan using poly-A priming compared to TruSeq random priming in BrainVar, gene expression was highly correlated per-sample and per-gene.

65BrainSpan also generated ChIP-seq datasets of histone 3 lysine 27 acetylation (H3K27ac), a marker of active genes, in DLPFC samples. These H3K27ac data mirror the expression profiles in BrainVar, with Falling genes enriched for fetal-biased H3K27ac peaks and depleted for adult-biased H3K27ac peaks, while Rising genes show the opposite pattern of enrichment (Figure S2).

66Statistically, the temporal dynamics of the expression of genes can be modeled as a mixture of K distinct trajectories, each with Gaussian noise. To delineate the temporal dynamics for K groups of genes, we used the Flexmix R package, which provides the expected trajectory for each group and the soft group assignments of individual genes to groups.

67The expression of 23,782 genes was transformed as log2(CPM) and normalized by the interquartile range. The samples in epoch 0 were excluded to avoid biased estimation due to very few samples. To identify the overall trend of expression over age, first we fitted the model on all the 23,782 genes assuming there are three groups and that the expected trajectories for each group can be represented with degree-4 polynomials on age. Three typical trajectories were identified, including a group of 6,941 genes with Rising expression levels, a group of 5,173 with Falling expression levels, and a group of 11,705 genes with roughly flat (fitted) expression over time, which we called Non-transitional (Figure 2B).

68We considered whether fewer or more trajectory groups described the data better. If we input more than four groups, adjacent classes are automatically combined due to the estimated priors falling below the minimum threshold leaving a maximum of four groups. Using the Akaike information criterion (AIC) and Bayesian information criterion (BIC) we observe the lowest values, indicating a better fit to the data, for four groups, followed closely by three groups, but not for two groups. However, the additional group does not reveal a new profile (e.g., genes rising to period 8/9 then falling for 10–12), instead it distinguishes two Rising groups based on gene variance. Since our objective was to identify temporal gene trajectories, not differences in variance, we selected three groups for the final model to avoid overfitting and ensure interpretability.

69For functional enrichment, we characterized genes sets for each trajectory using the R package, gProfiler. The pathway enrichment test was performed using Gene Ontology Biological Process terms, which contain between 10 and 2,000 genes, and all 23,782 cortically expressed genes were used as background. Enrichment tests were subject to the “moderate” hierarchical filtering parameter, and the FDR multiple correction in the gProfiler (Figure 2F).

70To assess enrichment across tissues, 27,546 transcripts with an RPKM ≥ 0.5 in 80% of samples one or more tissues in GTEx (https://gtexportal.org/home/) were defined and log-transformed (log2[RPKM+1]). For each gene, expression between each tissue and all other tissues was assessed using a moderated t-test (R package limma), with models adjusted for age, RIN, gender, and surrogate variables. The Benjamini and Hochberg method was used to estimate false discovery rate (FDR) and tissue enriched genes were defined as: log fold-change > 0.5 and FDR < 0.05. The enrichment of the Falling, Non-transitional, and Rising genes was assessed with using Fisher’s exact test with 23,782 cortically expressed genes as a background (Figure 2).

71To identify cell type-enriched genes we calculated a tau metric from the per gene log2TPM of log2UMI values for genes within cell type clusters in prenatal forebrain and postnatal cortex. Each gene was also ranked across all clusters within the dataset on the basis of TPM/UMI, so that the cluster with the highest TPM/UMI was ranked as “1,” while the cluster with the second highest TPM/UMI for that gene was ranked as “2,” etc. The genes were sorted by TPM/UMI rank and visually inspected for cell type specificity or substantial enrichment in single cell RNA-seq data from the prenatal and postnatal human cortex (https://cells.ucsc.edu/?ds=cortex-dev#). For each cell type, the top ten genes that showed clear specificity/enrichment were selected.

72H3K27ac peaks present in more than two samples of fetal or adult dorsal frontal cortex in BrainSpan were tested for fetal versus adult temporal bias using DESeq2. Temporally biased genes were defined as adjusted p < 0.01 and fold change ≥ 2. A category of “non-temporal” H3K27ac peaks was generated with peaks showing p > 0.05. All peaks were annotated using the gene with the closest transcription start site in Gencode v21. Genes were classified as only-fetal or only-adult if they were associated with fetal or adult-only H3K27ac peaks, respectively. Enrichments of each category of H3K27ac-genes in each category of eGenes were tested by means of a Fisher Exact’s test and p-values were adjusted using Benjamini-Hochberg, using genes associated to non-variant H3K27ac peaks as a reference background.

73To assess the functional topology in cortical samples, we applied Weighted Gene Co-Expression Network Analysis (WGCNA) to 23,782 cortically expressed transcripts. Network analysis was performed with WGCNA (version 1.63) using a signed network, choosing a soft-threshold power, the mean connectivity less than 50, and scale-free topology greater than 0.8. To reduce the bias driven by a few sample outliers, we applied the blockwiseConsensusModules function, which detects consensus modules across 100 subsampled networks. We used the average linkage hierarchical clustering of the topological overlap dissimilarity matrix (1-TOM) to generate the network dendrogram. Modules were defined as branches of the dendrogram using the hybrid adaptive tree cut with the following parameters: minimum module size = 200, negative pamStage, height cut = 0.999, and deep split = 2. Modules were summarized by their first principal component (ME, module eigengene), followed by merging modules with high correlations (eigengene value ≥ 0.9).

74For functional enrichment, we characterized WGCNA module genes using the gProfiler R package, as described above in the analysis of temporal trajectory genes. To identify WGCNA module genes that are regulatory targets, we searched for transcription factor binding targets using the ChEA, and TRANSFAC, and microRNA using the mirTarbase database.

75To assess whether 19 co-expression modules in our samples were preserved in other, independent DLPFC or frontal cortex expression datasets, we compared our dataset with non-overlapping samples from the BrainSpan dataset and applied the module preservation function from the WGCNA R package. From the BrainSpan dataset, we selected DLPFC (n = 30), 10 non-DLPFC neocortical regions (n = 317), and subcortical regions excluding the cerebellum (n = 140). Given the original co-expression network constructed above, we used modulePreservation to calculate module preservation statistics from 100 permutations (Table S3).

76To examine functional association of a group of genes/proteins, we performed a clustering analysis of a protein-protein interaction network implemented in SANTA R package. Ripley’s K function provides a measure of whether points are clustered together or randomly dispersed (homogeneous) in a network and the SANTA R package reformulated Ripley’s K function for a protein-protein interaction network. Clustering can be indicative of functional association between the genes/proteins under consideration. To test departure from homogeneity of a given gene set, we drew an empirical null distribution of clustering from 1,000 random samples of matching sized gene sets from the BioGRID protein-protein interaction network data (v3.4.132). We reported departure from null distribution as a Z-score (i.e., Z > 0: level of clustering of a gene set greater than null expectation and Z < 0: level of clustering of a gene set less than null expectation).

77Cis-eQTLs were identified for all high quality, common variants (N = 6,573,196) within 1 Mb of a gene boundary using the linreg function in Hail 0.1, with period, sex, and the first five principal components of common variant ancestry as covariates. This analysis was run on three cuts of the BrainVar dataset: complete sample (N = 176, periods 1–12), prenatal-only (N = 112, periods 1–6), and postnatal-only (N = 60, periods 8–12). Separately for the results of each analysis, false discovery rate (FDR) was calculated for all gene-variant pairs using the Benjamini-Hochberg procedure.

78We then classified all gene-variant pairs with FDR ≤ 0.05 from at least one analysis into groups defined by the temporal specificity of their eQTL effects. To do this, we first identified one variant per gene with the smallest, FDR-significant p value, from any of the three analyses. We then used a Z-test to compare the regression coefficients for these variant-gene pairs from the prenatal and postnatal analyses: Using the results from the eQTL analyses and from this prenatal-postnatal comparison, we then classified each of these top (one per gene) gene-variant pairs and their corresponding target gene (eGene) into one of five groups: Constant eQTLs/eGenes, characterized by consistent effects across this developmental dataset: FDR ≤ 0.05 in the complete sample analysis, same direction of effect and unadjusted p ≤ 0.05 in both the prenatal and postnatal analysesPrenatal-predominant eQTLs/eGenes, with strongest effects during prenatal development: FDR ≤ 0.05 in the prenatal analysis, unadjusted p > 0.05 in the postnatal analysis, pre-post comparison Z-test FDR-adjusted p ≤ 0.05Postnatal-predominant eQTLs/eGenes, with strongest effects during postnatal development: FDR ≤ 0.05 in the postnatal analysis, unadjusted p > 0.05 in the prenatal analysis, pre-post comparison Z-test FDR-adjusted p ≤ 0.05Prenatal-trending eQTLs/eGenes, which did not fit into earlier categories, but had higher prenatal effects (BPre > BPost)Postnatal-trending eQTLs/eGenes, which did not fit into earlier categories, but had higher prenatal effects (BPost > BPre)

79Constant eQTLs/eGenes, characterized by consistent effects across this developmental dataset: FDR ≤ 0.05 in the complete sample analysis, same direction of effect and unadjusted p ≤ 0.05 in both the prenatal and postnatal analyses

80Prenatal-predominant eQTLs/eGenes, with strongest effects during prenatal development: FDR ≤ 0.05 in the prenatal analysis, unadjusted p > 0.05 in the postnatal analysis, pre-post comparison Z-test FDR-adjusted p ≤ 0.05

81Postnatal-predominant eQTLs/eGenes, with strongest effects during postnatal development: FDR ≤ 0.05 in the postnatal analysis, unadjusted p > 0.05 in the prenatal analysis, pre-post comparison Z-test FDR-adjusted p ≤ 0.05

82Prenatal-trending eQTLs/eGenes, which did not fit into earlier categories, but had higher prenatal effects (BPre > BPost)

83Postnatal-trending eQTLs/eGenes, which did not fit into earlier categories, but had higher prenatal effects (BPost > BPre)

84All FDR-significant variants associated with the expression of a single gene were classified into one of these five groups according to the classification of the top variant for the same gene (Figure 5).

85Many eGenes are associated with multiple eQTLs, each of which could individually meet criteria for any one of the five temporal categories. As described above, we categorized eGenes into temporal categories based on the performance of their top eQTL (smallest p-value), but we assessed the performance against three alternative approaches: (1) eGene assigned to the same category as a majority of their eQTLs (“majority eQTL” approach), with ties assigned in the order Constant, Prenatal-predominant, Postnatal-predominant, Prenatal-trending, Postnatal-trending, (2) for eGenes with ≥ 1 eQTL, category assignment based on the performance of the second most significant variant (“second eQTL” approach), and (3) each eQTL individually assigned to categories (“individual eQTL” approach). For each of these alternative approaches, we calculated the percent of eGenes or eQTLs from each top variant-based category that were assigned to each category using the majority eQTL, the second eQTL, or the individual eQTL approach (Figure S5).

86To validate that our eQTL discovery analyses were adequately adjusted for differences in sample ancestry, we calculated genomic control, or lambda, an estimate of inflation of genetic association signal, from the results of each of the three analyses (complete sample, prenatal-only, postnatal-only). To do this, we randomly sampled 500 gene-variant pairs that occur on the same chromosome but are located 10–100 Mb apart, under the assumption that these distant gene-variant pairs would be unlikely to be enriched for true gene expression association. Using the eQTL signals for these 500 gene-variant pairs from the three analyses, we calculated lambda for analysis. We repeated this procedure for 100 random selections of 500 distant gene-variant pairs. We find that the median lambda values across 100 permutations for each are near 1, in keeping with a test has been properly corrected for population structure (specifically, 1.01 for prenatal, 0.99 for postnatal, and 1.00 for the complete sample).

87To assess the sensitivity of our cis-eQTL discovery analysis relative to previous work, we evaluated the relationship between sample size and eGene discovery for: the BrainVar prenatal, postnatal, and complete sample analyses as run using the HCP- and SVA-adjusted expression data, GTEx v7 analyses by tissue (gtexportal.org), postnatal human frontal cortex by the CommonMind Consortium, and prenatal human whole brain. Using the sample size reported in each analysis, or for each tissue (GTEx), and the number of genes with at least one eQTL reaching significance of FDR ≤ 0.05, we plotted the relationship between eGene discovery and sample size (Figure S5). We observe a strongly positive relationship across BrainVar and the published analyses, in keeping with prior reports that eQTL and eGene discovery is positively associated with sample size.

88We evaluated the performance of the eQTLs that we identified in our analyses with published sets of eQTLs identified in the human postnatal frontal cortex and in the human prenatal whole brain. For the postnatal frontal cortex data, we downloaded significant variant-gene pairs from the GTEx v7 data release from gtexportal.org and used R to write out the variant locations to a bed file format. We then used the command line LiftOver utility from the UCSC Genome Browser to convert the hg19 variant positions to GRCh38. For the prenatal brain data, we downloaded the eQTL summary statistics, results for the top eQTLs per gene, and SNP positions bed file from the study data repository on Figshare (https://figshare.com/articles/Summary_statistics_for_expression_quantitative_trait_loci_in_the_developing_human_brain_and_their_enrichment_in_neuropsychiatric_disorders/6881825). We then used each gene’s nominal significance threshold from the top eQTLs file to identify the full set of variant-gene pairs meeting significance.

89Using variant (GRCh38 position, reference, and alternate alleles) and gene (Ensembl gene IDs) identifiers, we matched significant variant-gene pairs separately from GTEx frontal cortex and prenatal brain to the variant-gene pairs meeting FDR significance (%0.05) in the BrainVar analyses. For all significant eQTLs in BrainVar that overlapped with the reference datasets, we then compared the effect of the variant on the expression of its associated gene to determine the percentage of overlapping eQTLs with concordant direction of effect, as well as the Pearson correlation between the eQTL effects (beta from BrainVar, slope from GTEx or prenatal whole brain). The significance of this correlation was evaluated using the cor.test function in R.

90The distance between each significant eQTL and the transcription start site (TSS) of its associated eGene was calculated by comparing the variant position to the TSS position and strand of the gene according to Gencode v21. Positive distances indicate variants downstream of the TSS, negative distances indicate upstream variants. Comparison between groups of eQTLs was run using only the top eQTL per gene and the absolute value of distance from the associated gene’s TSS.

91We tested the global overlap between eQTLs and H3K27ac from human fetal, infant and adult dorsal frontal cortex and cerebellum and embryonic cortex from BrainSpan and. Intersection between sets of coordinates were performed using Bedtools. We tested three sets of variants: (1) best eQTL per eGENE, (2) all significant eQTLs with FDR < 0.05 in the corresponding eQTL category of Prenatal-predominant, Postnatal-predominant or Constant, and (3) a background group composed by all variants tested for eQTL, excluding those with a p < 0.05 to any gene at any time period tested. 95% confidence intervals were obtained by bootstrapping variants 100 times (Figure S5).

92We tested the enrichment of different categories of eQTLs in sets of genomics elements using GREGOR. We tested all best eQTL per gene in: (1) dorsal frontal cortex H3K27ac peaks from fetal and adult brain samples, and (2) 18 chromatin states whole-genome segmentation of sample E073-Medial Frontal Cortex Lobe from The Roadmap Epigenomics Project. We reported observed/expected number of overlapping eQTLs and BH adjusted p-values (Figure S5).

93We used an F test to compare prenatal variance (periods 1–6) and postnatal variance (periods 8–12) in gene expression level, correcting for age (period), sex, and five ancestry principal components within each stage. At a Benjamini-Hochberg-adjusted p value ≤ 0.05, we identified 8,094 genes with greater prenatal variance, and 1,752 genes with greater postnatal variance.

94To compare with 23,782 cortically expressed genes, we created the list of gene sets for previous trait and disorder association and functional properties. Gene identifiers were converted between studies based on the complete HUGO Gene Nomenclature Committee dataset. Autism spectrum disorder (ASD) risk genes were obtained from, an exome sequencing based gene discovery refining to high-confidence genes (n = 99) at a false discovery rate (FDR) ≤ 0.1. Genes associated with developmental delay (n = 93) were selected from the exome sequencing analysis of the Deciphering Developmental Disorders project. From, we chose 33 genes as high-confidence epilepsy candidates, where multiple de novo variants were seen.

95We used significant loci from the genome-wide association studies (GWAS) of attention deficit hyperactivity disorder (ADHD), Alzheimer’s disease, educational attainment, schizophrenia, major depressive disorder, multiple sclerosis, Parkinson’s disease. We selected loci from a summary statistics file if publicly available, otherwise we used the table of genome-wide significant loci from each study. We retrieved genes harboring loci or within 10kb from loci. We excluded the extended major histocompatibility complex region on the chromosome 6, known to have a substantial number of genes due to high linkage disequilibrium for downstream enrichment analyses.

96Constrained genes were defined as probability of loss-of-function intolerant (pLI) score ≥ 0.995 in the ExAC database. Genes specific to cell types in the mid-fetal cortical development were selected from and BrainSpan. For all gene lists, see Table S2 and Table S4.

97We tested for enrichment of DLPFC eQTLs among GWAS significant SNPs using a permutation-based procedure. GWAS SNPs were taken from published GWAS summary statistics at a significance threshold of p < 5 × 10−8. The following procedure was repeated separately for summary statistics from GWAS of four phenotypes: schizophrenia, educational attainment, multiple sclerosis, and Alzheimer’s disease. First, SNPs in the summary statistics and the list of SNPs tested for eQTL discovery were filtered to those in 1000 Genomes Project data. SNPs in the GWAS summary statistics were then filtered to variants also tested for eQTL discovery, using direct (same SNP) or proxy (r2 > 0.8 in CEU 1000 Genomes samples) as defined by the PLINK 1.9 – r2 command. To estimate the proportion of eQTL SNPs among phenotype associated SNPs versus among all SNPs (or a sample of null SNPs), we recognized three important factors that could differ between null SNP sets and the phenotype-associated SNPs: LD structure, MAF distribution, and gene density. First, to account for LD structure, we used PriorityPruner version 0.1.4 (http://prioritypruner.sourceforge.net) to LD prune SNPs supervised by GWAS p-value in order to preferentially retain as many phenotype-associated SNPs as possible while adequately removing SNPs in high LD (r2 > 0.7 within a sliding 500kb window) using in 1000 Genomes CEU sample data. Second, the remaining SNPs were grouped according to their MAF decile. Third, each remaining SNP was grouped into decile of gene density to allow for differential opportunity to be identified as an eQTL. Gene density was determined by the number of genes within the 1MB eQTL detection window as defined by the annotation package “TxDb.Hsapiens.UCSC.hg38.knownGene” from R Bioconductor. SNPs in the ~3.2 MB HLA region (hg38 coordinates: chr6:29,751,784–32,915,731) as defined by the “GWASTools” R Bioconductor package and UCSC genome browser were excluded from enrichment testing. Next, 1 million null SNP sets were drawn by matching each phenotype-associated SNP (GWAS p < 5 × 10−8) to a random SNP matched on both MAF and gene density. Enrichment fold statistics were computed as the proportion of eQTLs in the phenotype-associated set divided by the mean proportion of eQTLs across null sets. P values were calculated as the proportion of null set fold-enrichment statistics as or more extreme than the observed phenotype-associated fold enrichment statistic. This permutation procedure was repeated for each of six eQTL SNP lists: all, Constant, Prenatal-predominant, Prenatal-trending, Postnatal-predominant, and Postnatal-trending.

98To assess whether eQTL targets (eGenes) are enriched for GWAS signal, we performed competitive gene set enrichment analysis for each group of eGenes using the MAGMA software. We input the eGene lists from Prenatal-predominant, Postnatal-predominant, Constant, Prenatal-trending, and Postnatal-trending (Table S5) and summary statistics from published GWAS of schizophrenia, autism spectrum disorder, educational attainment, multiple sclerosis, Alzheimers disease, triglycerides, and height. First, SNPs from the GWAS summary statistics files were annotated to NCBI protein-coding genes that passed RNA-seq QC in our DLPFC expression data with a 10kb window on either side of the gene boundaries. Next, a gene-level analysis was done to determine the strength of association for each gene with phenotype of interest. To assess whether genes in the eGene gene-sets are more strongly associated with the phenotype of interest than other genes, gene-based z-scores are regressed on a gene-set indicator variable and MAGMA default covariates (gene size, gene density, sample size, 1/MAC, and the log of each of these). The beta coefficient for the gene-set indicator variable is tested for significance HA: β1 > 0. Results from this analysis are reported in Table S6. We did not observe significant enrichment of GWAS signal from any of the six phenotypes tested in any of the five temporally assigned eGene gene sets. However, we repeated a similar test, annotating GWAS summary statistics SNPs to NCBI protein-coding genes with 10kb flanking region that passed RNA-seq QC in our DLPFC expression data and had an assigned pLI score, and found that a gene-set defined as pLI score > = 0.995 showed significant enrichment for stronger GWAS association in all tested phenotypes except for multiple sclerosis, compared to genes with pLI score < 0.995.

99Coloc was used to formally test for co-localization of GWAS signal from schizophrenia and educational attainment summary statistics with our five eQTL categories. Coloc tests five hypotheses (H0: no association, H1: GWAS association only, H2: eQTL association only, H3: both but not co-localized, H4: both and co-localized) and returns a posterior probability for each hypothesis in each region. Posterior probably of H4 > = 0.8 is strong Bayesian evidence of co-localization.

100Open source scripts used in this study are referenced throughout. The pipeline for whole-genome analysis is available online at: https://github.com/sanderslab/psychcore-compute-platform. The accession number for the raw RNA-seq and WGS data reported in this paper, along with processed files, is PsychENCODE Knowledge Portal: syn21557948 on Synapse.org (https://www.synapse.org/#!Synapse:syn4921369).