Oxygen-induced stress reveals context-specific gene regulatory effects in human brain organoids
Abstract
The interaction between genetic variants and environmental stressors is key to understanding the mechanisms underlying neurological diseases. In this study, we use human brain organoids to explore how varying oxygen levels expose context-dependent gene regulatory effects. By subjecting a genetically diverse panel of 21 brain organoids to hypoxic and hyperoxic conditions, we identify hundreds of gene regulatory changes that are undetectable under baseline conditions, with 148 trait-associated genes showing regulatory effects only in response to oxygen stress. To capture more nuanced transcriptional patterns, we employ topic modeling, which reveals context-specific gene regulation linked to dynamic cellular processes and environmental responses, offering a deeper understanding of how gene regulation is modulated in the brain. These findings underscore the importance of genotype-environment interactions in genetic studies of neurological disorders and provide new insights into the hidden regulatory mechanisms influenced by environmental factors in the brain.
Understanding how gene regulatory variants function across different cellular and environmental contexts is essential for interpreting genetic associations with disease. Gene-by-environment (GxE) interactions occur when genetic variants influence how individuals respond to specific environmental exposures, leading to interindividual differences in phenotypes, including variability in disease susceptibility. This concept is particularly important in complex diseases, where individuals with different genetic backgrounds may exhibit varying risk profiles for conditions such as neuropsychiatric disorders (Hollander et al. 2020; Zhang et al. 2022). For instance, environmental factors such as stress (Moyerbrailean et al. 2016; Nguyen et al. 2022; Seah et al. 2023), oxygen deprivation (Brooks et al. 2009; Ward et al. 2021), or infection (Barreiro et al. 2012; Fairfax et al. 2014; Lee et al. 2014; Kim-Hellmuth et al. 2017; Alasoo et al. 2018; Randolph et al. 2021) can trigger disease-relevant gene regulatory effects that remain hidden in static, steady-state conditions.
Gene regulatory catalogs like the Genotype-Tissue Expression (GTEx) project (The GTEx Consortium 2020) provide valuable insights into how genetic variants affect gene expression across various tissues in steady-state conditions. However, the majority of disease-associated loci remain unexplained, likely due to their regulatory effects being specific to certain cell types or environmental contexts that have not been fully explored (Findley et al. 2021; Umans et al. 2021). This gap is particularly pronounced in the brain, where the complex interplay between different cell types and environmental stressors can contribute to the onset and progression of neurological and psychiatric diseases (Walker et al. 2019; Aygün et al. 2021; Jerber et al. 2021; Schumann et al. 2023; Emani et al. 2024; Wen et al. 2024).
Brain cells, especially neurons, are highly sensitive to environmental perturbations like hypoxia (oxygen deprivation), given the brain's high metabolic demand and susceptibility to oxidative damage (Cobley et al. 2018). Hypoxia is a well-known neurological risk factor throughout life, arising from conditions such as sleep apnea, high altitude, respiratory infections, and premature birth (Wehby 2013; Salmaso et al. 2014; Lavie 2015; Prabhakar et al. 2020; Erickson et al. 2021). Hypoxic exposure has profound effects on cognitive function and white matter integrity and increases the risk for neurodegeneration and psychiatric disorders (Cannon et al. 2002; Thomas et al. 2005; Halbower et al. 2006; Lin and Beal 2006; Jackson et al. 2011; Varçin et al. 2012; Dias et al. 2013; Kim et al. 2013; Blesa et al. 2015; Daulatzai 2015; Kwak et al. 2015; Giannopoulou et al. 2018; Luk 2019; Musgrove et al. 2019; Owen et al. 2019; Garcia et al. 2020; Wan et al. 2020; Ionescu-Tucker and Cotman 2021; Zacharias et al. 2021). Despite the importance of oxygen homeostasis and hypoxia to brain function, we lack comprehensive insight into how different brain cell types respond to these environmental stressors at the gene regulatory level, which limits our ability to interpret genetic associations with neurological traits (Romanko et al. 2004; Komitova et al. 2013; Kwak et al. 2015; Ren et al. 2024).
In this study, we used human brain organoids to investigate the transcriptional responses of diverse brain cell types to oxygen perturbation across 21 individuals. By applying single-cell RNA sequencing, we captured gene expression data both under baseline conditions and after exposure to varying oxygen levels. Through whole-genome sequencing of each donor, we identified genetic contributions to these responses, revealing dynamic gene regulatory effects with significant relevance to neurological and psychiatric disease susceptibility. This approach goes beyond characterizing gene regulatory variation in static, postmortem tissue and opens new avenues for studying GxE interactions in a controlled, in vitro setting.
Results
We differentiated brain organoids from the induced pluripotent stem cells (iPSCs) of 21 unrelated Yoruba individuals from Ibadan, Nigeria, which have been extensively characterized and have proven an efficient resource for mapping gene regulatory variation (Banovich et al. 2018) (see Methods). We performed oxygen manipulation experiments in two batches of 7–16 individuals, with two individuals replicated across batches to allow us to control and account for batch effects (Fig. 1A). Specifically, following 7 weeks of growth at atmospheric oxygen levels (20% O2) (see Methods), organoids were adapted to 10% O2 to mimic the physiological environment experienced by brain cells in vivo. After 1 week of culture at 10% O2, organoids were either maintained at physiological oxygen (baseline/normoxia), transferred to low oxygen (1% O2; hypoxia), or transferred to high oxygen (20% O2; hyperoxia) for 24 h. Following the treatments, we dissociated organoids in the presence of transcriptional inhibitors and multiplexed equal proportions of each sample in preparation for single-cell RNA sequencing, targeting 3000 cells per individual and oxygen condition, and a depth of 20,000 reads per cell. After demultiplexing and quality control, we retained data from 170,841 cells (normoxia: 52,671, hypoxia: 57,788, hyperoxia: 60,382; median 5666 UMI counts per cell).
A panel of brain organoids yields diverse cell types across individuals. (A) Workflow of data collection. Brain organoids were first differentiated (see Methods) and adapted to physiological oxygen before 24-h exposure to high, low, or control oxygen levels. Data shown in this figure are taken from the 10% “normoxia” condition. (B) UMAP representation of organoid single-cell transcriptomes highlighting principal cell types obtained across all individuals in our iPSC panel in the normoxia condition. (C) Proportion of cells from each parental cell line by annotation, with colors corresponding to the UMAP shown in B.
Brain organoids comprise diverse cortical and noncortical cell types
We first sought to characterize the cell type composition of brain organoids maintained at baseline oxygen levels. We annotated cell clusters using fetal brain reference data (Polioudakis et al. 2019; Trevino et al. 2021) and known marker genes, finding a variety of cortical cell types, including radial glial progenitors, intermediate progenitors, excitatory neurons, and inhibitory neurons (Fig. 1B; Supplemental Fig. S1; Methods). We also identified a substantial cluster of neurons with noncortical identities, including thalamic and midbrain inhibitory neurons and GNRH+ cells (Fig. 1B; Supplemental Fig. S1; Methods). Whereas the cell type proportions varied across individuals, out of 20 high-confidence cell types, 10 were present in over half of the individuals, with a median of 11 cell types detected per individual (Fig. 1C; Supplemental Fig. S1E).
We used propeller (Phipson et al. 2022) to assess differences in cell type proportions across treatment conditions and found that organoid composition was largely unaltered by oxygen manipulation (Supplemental Fig. S2A; Supplemental Table S1). Moreover, the cell type composition of an additional set of control organoids which we maintained at atmospheric oxygen levels for the duration of the experiment did not differ substantially from what we observed in the treatment conditions (Supplemental Fig. S1B). We observed differences in cell type composition across individuals, with the rarer cell types (such as midbrain dopaminergic-like cells, mature oligodendrocytes, vascular leptomeningeal cells, and Cajal-Retzius cells) present in only a minority of samples (Fig. 1; Supplemental Fig. S1E). However, across treatments, the cell type composition of organoids from the same individuals remained similar (with the exception of NA19144) (Supplemental Fig. S2A), suggesting that the treatments did not have a marked effect on cellular identity.
Although the treatments applied in the 24 h prior to sample collection did not result in noticeable differences in cell composition, we noted a considerable range of major cell type abundances across individuals and investigated the possible drivers of these differences using the propeller framework. We found no effect of sex or iPSC passage number on any cell type's relative abundance but observed that certain rare cell types—including oligodendrocytes, inhibitory neuron subtypes (midbrain, thalamic, and SST+), and midbrain dopaminergic-like cells—differed in proportion across individuals included in different collection batches (Supplemental Table S1). To assess the extent of confounding by batch, we merged biologically similar clusters to generate a set of 10 coarse annotations representing the principal cell types in our data (Supplemental Fig. S1D; Methods). When we examined cell type proportions among coarsely defined clusters, the batch effect disappeared, suggesting that principal cortical cell types were largely stable across experiments. Still, we performed all subsequent analyses using both annotations to assess potential bias caused by overclustering (Methods) and to ensure that the two approaches produce similar outcomes. In what follows, we report findings from the more interpretable fine-grained annotation set (results from the coarse clustering approach are provided in Supplemental Figs. S1D,E, S2D,E, S4A,B, S5B–F).
Shared transcriptional response patterns reveal cell type–specific vulnerabilities to hypoxia
To identify differentially expressed (DE) genes between baseline and either hypoxic or hyperoxic conditions, we applied a linear mixed model to pseudobulk expression data for each cell type and treatment condition (Hoffman and Roussos 2021). We identified a total of 10,230 DE genes in response to hypoxia, ranging from 91 to 5590 genes per cell type (FDR < 0.05) (Supplemental Table S2). Similarly, we identified 10,425 hyperoxia-responsive genes, ranging from 17 to 6102 genes per cell type (FDR < 0.05) (Supplemental Table S2). In at least one (any) cell type, 2703 hypoxia-responsive genes (and correspondingly, 2855 hyperoxia-responsive genes) exhibited a greater than 1.5-fold change in expression compared to baseline. As expected, we detected far more DE genes in abundant cell types (Supplemental Fig. S2C), within which 76%–92% of hypoxia-evoked transcriptional effects and 84%–93% of hyperoxia-evoked transcriptional effects were modest (<1.5-fold).
We were intrigued by the large differences in the numbers of DE genes across cell types, an observation that cannot be fully explained by cell abundance (Supplemental Fig. S3). To explore this further, we analyzed the data using multivariate adaptive shrinkage (mash) (Urbut et al. 2019), a method designed to assess patterns of sharing across studies under incomplete power, which improves detection of small but consistent signals across cell types and conditions. By combining power across cell types, we were able to detect weak DE effects that emerge in multiple cell types and accurately identify condition- and cell type–specific DE genes.
As expected, we observed similar oxygen response patterns among related cell types (Fig. 2A). However, even using mash, we found that more than half of the oxygen-responsive genes we identified (68% of hypoxia-responsive genes and 63% of hyperoxia-responsive genes, FDR < 0.05, effect size >1.5-fold) were DE in fewer than three cell types, consistent with the idea that oxygen homeostasis mechanisms are tuned to the needs of distinct brain cell types (Chi et al. 2006). We also observed a tendency for DE genes with large effects to be more cell type–specific than DE genes with smaller effects (Supplemental Fig. S2B), an observation that is counter-intuitive with respect to power considerations.
Transcriptional responses of brain organoids to oxygen changes. (A) Fraction of shared differential expression effects between cell types and conditions. Sharing was assessed from mash posterior estimates of significance and effect sizes (see Methods; Supplemental Fig. S2). Stress treatments are indicated by colored bars. Cell type classification colors are detailed in panel E and Figure 1. (B) Enrichment of gene annotations among differentially expressed genes across cell types and conditions. The 39 MSigDB Hallmark gene sets with at least one significant enrichment are shown. (C) Fractional change in hypoxia-stressed fraction of each cell type after hypoxia exposure (coarse cell classification). Cells were classified as hypoxia-stressed or hyperoxia-stressed by Gruffi using a gene set derived from treatment responses shared across all cell types. Points represent individual-level fractional changes, with cell type medians shown as crosses. See also Supplemental Figure S4. (D) Distribution of cell distances to organoid perimeters measured after immunofluorescent labeling. “Unlabeled” measurements are derived from DAPI-labeled nuclei with no immunofluorescent label. Black dots represent sample means. (E) Average topic loading for each cell type and treatment condition. Topic 7 tracks hypoxia exposure, whereas other topics reflect processes found in one or several cell types.
Most cell types exhibited a modest positive correlation between the transcriptional responses to hypoxia and hyperoxia (Fig. 2A; Supplemental Fig. S2B,E), and many gene families related to general response to stress were enriched among DE genes in both treatment conditions. For example, all oxygen-treated cell types are enriched for DE genes with roles in cellular metabolism and inflammatory responses, and most are enriched for DE genes involved in cell proliferation and DNA damage repair (Fig. 2B). Consistent with expectation, hypoxic treatment, specifically, also induced a well-recognized core regulatory program across all cell types, including the upregulation of hypoxia-inducible factor (HIF) target genes (Cummins and Taylor 2005; Semenza 2012; Chang et al. 2015). Consistent with this, the hypoxia response in our study showed greater similarity to previously reported hypoxia-induced gene regulation in brain organoids than the hyperoxia response in any cell type (Supplemental Fig. S2G).
Because our observations pointed to both a general stress response and modality-specific effects following the treatments, we sought to characterize the proportion of single cells expressing regulatory signatures of stress in baseline and treatment conditions, so that we could specifically focus on the response to the treatments. To parse cell type heterogeneity in the transcriptional response to oxygen perturbation and to differentiate between stressed and unstressed cells, we identified gene sets that capture robust and general responses to each treatment condition (Supplemental Table S3; Methods). We leveraged these gene sets, which were not specified a priori on the basis of biological annotations, in a granular filtering approach (Vértesy et al. 2022) to classify cells as either stressed or unstressed and then repeated our differential expression analysis after censoring the stressed cells. We found that, although stressed cells contribute disproportionately to the treatment expression response, they do not account for it entirely: in cell types with a higher proportion (12%–36%) of stress-censored cells, the number of DE genes (FDR < 0.05, fold-change > 1.5) decreased by as much as 64% relative to randomly censored data (Supplemental Fig. S4D). Nevertheless, stressed cells retained expression of key cell type markers, suggesting a change in state but overall preservation of cell identity (Supplemental Fig. S4C).
Cell types that include an increased fraction of stressed cells following exposure to hypoxia or hyperoxia may be especially sensitive to oxygen perturbation. We calculated the fractional change in the proportion of stressed cells between baseline and treatment conditions to assess cell type–specific sensitivity within brain organoids. We found that intermediate progenitors, immature neurons, and radial glia—three early cell types that fall in a developmental sequence—are particularly responsive to both hypoxia and hyperoxia, suggesting that early developmental transition points are especially labile to oxygen stress (Fig. 2C; Supplemental Fig. S4).
The observation of large differences in cell type sensitivity to changes in oxygen may be significant and help us better understand disease mechanisms. It is possible, however, that cell type–specific sensitivity to the treatments could be explained by variation in organoid spatial structure, as we observed a wide range of sensitivities across organoid lines. To examine this, we used antibody markers to map eight major cell types in cryosections obtained from the same organoids we used for sequencing (Methods). We then quantified the accessibility of each cell type to exogenous oxygen by measuring the distance from immunostained cells to the organoid periphery and compared the distance distribution across cell types. Whereas different cell types were distributed at various depths within each organoid, variation between organoids was comparable to the differences observed between cell types (ANOVA F-test of organoid-level medians, P = 0.342) (Fig. 2D). On the whole, the most oxygen-sensitive cell types were localized neither more superficially nor more deeply than other cell types. Thus, cell type–specific sensitivity to oxygen perturbation appears to be driven primarily by cellular identity rather than cell position within organoids.
Context-specific responses to oxygen perturbation
We observed that oxygen perturbation induced widespread transcriptional effects, many of which were shared among subsets of cells within and between cell types. Perturbed cells continued to express cell type–specific markers, simultaneously maintaining their respective identities while adopting a more general signature of oxidative stress (Supplemental Fig. S4C). This suggests that discrete cell type labels may fail to capture continuous contextual responses, including signatures shared by developmentally related or physically proximal cells. In an effort to capture the subtler transcriptional patterns in our data, and as an alternative to the clustering approach used above, we used topic modeling to decompose cellular transcriptomes into 15 groups (topics) that capture distinct sources of transcriptional variation (see Methods).
Several topics largely recapitulate discrete cell type classifications, consistent with the notion that cell identity is the primary source of transcriptional heterogeneity in our data (Fig. 2E; Supplemental S2F). Other topics, such as topic 7, recapitulate properties that were already known to us, such as the collection of DE genes that we previously used to identify hypoxia-stressed cells (Supplemental Fig. S3D). We also identified topics that captured dynamic cellular processes and developmental states. For example, topic 4 is associated with dividing radial glia and dividing intermediate progenitors and is distinct from separate topics that tightly correlate with each of these cell identities in isolation. Indeed, topic 4 is defined by elevated expression of genes involved in DNA replication and cell division, including MKI67, TOP2A, and centrosomal proteins (Supplemental Fig. S3D), capturing shared aspects of the dividing cell environment across different cell types. In turn, topic 3 shows modest loading in cortical hem progenitors and higher loading in choroid plexus cells, reflecting their shared developmental origins: signals from the cortical hem influence the differentiation and patterning of cells at the boundary of the cerebral cortex and hippocampus, including those forming the choroid plexus (Moore and Iulianella 2021). Altogether, topic modeling allowed us to recapture functional relationships and shared states that were concealed in our analysis of discrete, mutually exclusive cell type clusters, providing us with an enriched view of brain organoid dynamics.
Transcriptional responses to oxygen perturbation are genetically regulated
Having detected thousands of genes that are differentially expressed in response to oxygen perturbation, we sought to uncover potential genetic sources for inter-individual variation in treatment response. We aggregated single-cell gene expression data into pseudobulk groups defined by their unique combination of donor, cell type, and treatment condition. After excluding groups comprising fewer than 20 cells, we removed cell types that had fewer than eight individuals in each treatment condition (Methods). For each of the 13 remaining cell types, we separately mapped cis expression quantitative trait loci (eQTLs) under hypoxic, hyperoxic, and baseline conditions, including gene expression principal components as model covariates to account for the effects of sex, batch, and latent confounding factors (Methods). We expected most eQTLs to be shared across treatment conditions—indeed, our DE analysis indicates that most genes are robust to oxygen perturbation. We therefore used mash to integrate evidence for SNP-gene associations across all three treatment conditions, considering each cell type in turn. As we and others have found, this approach improves power to detect individually weak signals that emerge consistently in multiple experimental contexts and to conservatively estimate the degree of sharing among those conditions (Panousis et al. 2023; Natri et al. 2024; Popp et al. 2024).
Across 13 cell types, we tested a total of 9469 genes and found 2438 cis eQTLs in 1778 genes (local false sign rate < 0.1). Among these, we identified 1736 standard eQTLs (in 1266 eGenes), in which the eQTL effect size is of similar size and direction across all treatment conditions (we used a conservative twofold cutoff to define similar effect sizes across treatment conditions, as we consider shared effects to be the null) (Supplemental Table S4). We also identified 702 oxygen-response eQTLs in 658 genes, namely eQTLs that have significantly different effect sizes between treatment conditions. Note that the sum of eGenes associated with at least one standard eQTL and those associated with at least one oxygen-responsive eQTL is larger than 1778 because eGenes are often associated with more than a single eQTL across cell types and often with both standard and oxygen-responsive eQTLs in different cell types.
Oxygen-responsive eQTLs included 111 loci associated with distinct effects in the hypoxia condition (in at least one cell type), 144 loci associated with distinct effects in hyperoxia, and 72 loci associated with a difference between the baseline normoxia and both treatments (Fig. 3A). Of particular note, across cell types, 497 oxygen response eQTLs are not associated with a statistically significant eQTL in the baseline (normoxia) condition. The genetic effect of these loci on gene regulation in the different cell types can only be detected under the stress conditions imposed by the change in oxygen levels (Fig. 3B). Consistent with this, oxygen-responsive eQTL genes—particularly those not found under normoxia—are associated with the expression of genes that are less likely to have eQTL effects in cerebral cortex tissues assayed in the GTEx project (one-sided paired Wilcoxon test, P = 0.01636) (Fig. 3C; Supplemental Fig. S5A). We next asked whether eQTLs identified in organoids are comparable to findings from other primary tissue sources (Supplemental Fig. S5G–J). We found that a well-powered eQTL meta-analysis of cerebral cortex tissue (Sieberts et al. 2020) recovered 91% of the eGenes identified in our study. Moreover, organoid eQTL SNP-gene pairs that were previously identified in GTEx showed correlated effect sizes (Pearson's correlation 0.53 across all tissues) (Supplemental Fig. S5I,J). Taken together, by incorporating environmentally perturbed conditions, our in vitro system identifies eQTL effects in many genes that require much larger primary tissue samples.
Discovery of treatment context–specific eQTL effects across cell types. (A) Number of eGenes identified in each cell type, classified as “standard” (blue) or oxygen-responsive (red) according to contexts in which effect sizes differ by less than twofold. Conditions with shared effect sizes are connected by gray lines, leaving black “condition-specific” effects. Total eQTLs of each category are indicated by asterisks (right axis). (B) Number of eGenes identified in each cell type, classified by sharing across oxygen treatment conditions and detection in normoxia condition. Light and dark shades of red and blue categories, respectively, sum to the blue and red categories shown in A. (C) Fraction of eGenes identified in this study that are classified as eGenes in GTEx cerebral cortex tissues (see Methods). Treatment context-specific eGenes that were not detected in control conditions were less likely to be found in GTEx (P = 0.01636, one-sided paired Wilcoxon test of oxygen-insensitive category against oxygen-responsive/undetected in normoxia category).
Context-dependent genetic regulation
In the analysis of gene expression levels, topic modeling allowed us to place cells along continuous axes of variation that are not predicated upon marker gene or reference annotations. We reasoned that topics could also be deployed to identify genetic regulatory effects that emerge in contexts defined more precisely than is possible using the discrete categories of cell type and treatment. To explore the effects of cis eQTLs in an expanded set of precisely defined cellular contexts, we tested for interactions between eQTLs and topics. Rather than individually testing each eQTL-topic interaction, we used CellRegMap (Cuomo et al. 2022) to jointly test all linear combinations of topics, improving our power to detect a wide range of genotype-context interactions.
We identified 289 genes with a topic-interacting eQTL. To infer the relevant cellular context for each eQTL, we assessed the correlation between its estimated effect and the loading for each topic (Supplemental Table S5). When possible, we checked to ensure our interpretation was corroborated by results from our analysis of discrete cell types and treatment conditions. For example, topic 15 describes the specialized radial glia of the cortical hem and glial progenitor cells. Out of the top 12 eGenes associated with topic 15 (Pearson's correlation > 0.6), five were identified as significant eGenes using our pseudobulk approach, with effects exclusively in cortical hem or radial glial cells. For example, among the top 12 eGenes is the cholesterol transporter ABCA1 (Pearson's correlation 0.77), which showed modest eQTL effects in our pseudobulk analyses of cortical hem cells (Fig. 4A). Concordantly, our topic-based approach revealed a modest correlation between ABCA1 QTL estimates and topic 6 (Pearson's correlation 0.26), which is defined primarily by radial glia.
Organoid eQTLs can help interpret human disease genetics. (A) Example topic-interacting eQTLs. ABCA1 expression is correlated with the inferred cell environment, as defined by linear combinations of topics (Bonferroni-corrected P-value 0.00277), and this effect is largely driven by cortical hem and glial progenitor topic 15 in a genotype-dependent manner. The WDR45B gene-by-context eQTL effect (Bonferroni-corrected P-value 0.00184) is largely explained by hypoxia-associated topic 7. Each point corresponds to a single pseudocell used for CellRegMap topic interaction QTL mapping. (B) Fraction of eGenes in each cell type which are the nearest gene to a genome-wide significant GWAS finding among 402 brain-related traits. (C) Example of a cell type– and treatment-specific regulatory association matching a significant GWAS variant (lfsr 0.54, 0.04, and 0.66, with posterior estimated effect sizes of −0.11, −0.55, and −0.02, for control, hypoxia, and hyperoxia conditions, respectively). (D) Number of eGenes in each cell type and discovery condition for which rare loss-of-function alleles have been associated with disease (see Methods).
Using the topic eQTL analysis, we found 58 eGenes whose regulatory effects were significantly correlated with topic 7, which is associated with hypoxic stress (P < 0.05, Bonferroni correction, and Pearson's correlation > 0.1). Of these, 38 did not have any eQTLs in our standard cell type–specific analysis, including WDR45B and CRABP2, which were more strongly correlated with topic 7 than any other topic. WDR45B is a member of the WIPI protein family of autophagosomal proteins, and WDR45B mutations have been linked to numerous severe neurodevelopmental disorders (Cong et al. 2021). The retinoic acid binding protein CRABP2 regulates expression and localization of HIF1A, thereby affecting cellular metabolism and hypoxia response (Fu et al. 2024). The remaining 20 topic 7-correlated eGenes were previously identified in our pseudobulk analysis, but only five showed treatment-specific effects, suggesting that a higher resolution description of cell states can assign additional context that might be missed in a categorical analysis. We also found eQTLs correlated with more than one topic, such as CD44, a known regulatory target of HIF1A that interacts with HIF2A to modify tissue responses to hypoxia (Johansson et al. 2017; Liang et al. 2017), which was correlated with both topic 7 (hypoxia) and topic 15. Taken together, these results highlight the utility of decomposing complex cellular states into constituent programs for the analysis of gene-by-environment interactions.
Context-dependent eQTLs help to interpret the effects of disease-associated loci
The use of brain organoids allowed us to identify context-specific gene regulatory effects that are underrepresented in standard eQTL studies of postmortem tissues (Fig. 3C; Supplemental Fig. S5). We reasoned that oxygen-responsive eQTLs could help to explain uncharacterized genetic associations with disease—particularly if the eQTL effects were undetectable at baseline. To explore this possibility, we examined the overlap of eGenes with GWAS loci assembled from 402 brain-relevant traits in each cell type (Methods; Supplemental Table S7). Across cell types, we identified 571 eGenes implicated in GWAS of brain-related traits. Furthermore, eGenes with response eQTLs that were latent at baseline included a similar proportion of disease-associated genes (median 0.308 per cell type, 148 total) as standard eGenes detected at baseline (median 0.334 per cell type, 291 total; two-sided paired Wilcoxon test P = 0.2163) (Fig. 4B). Focusing on the 1029 novel eGenes that were not represented in GTEx cerebral cortex tissues, we found a median of five disease-associated eGenes per cell type (total 96) to overlap oxygen-response eQTLs that were latent at baseline. Finally, of the 218 eGenes that interact with hypoxia (i.e., topic 7), 55 correspond to a GWAS gene, including 31 genes that are not eGenes in GTEx cortical tissue. Thus, mapping eQTLs in oxygen-treated brain organoids allowed us to uncover novel, disease-relevant effects that could not be detected in primary cortical tissues.
Next, we asked whether context-specific eQTLs could revise previous interpretation of disease-associated SNPs. Focusing on oxygen-responsive lead eQTL SNPs that were associated with brain traits in GWAS, we identified 19 associations (corresponding to seven genes) in which eQTL mapping implicated a different target gene than the original GWAS report or a simple nearest-gene heuristic. For example, we identified an association between rs10141157 and the expression of COA8 in a subset of inhibitory neurons (Fig. 4C). The SNP rs10141157 is associated with neuropsychiatric disorder risk (Yao et al. 2021) and has been described as an eQTL for nine different genes, principally in blood, with brain data pointing to effects on expression of KIF26A and TRMT61A (Sieberts et al. 2020; Ghoussaini et al. 2021; Mountjoy et al. 2021). COA8 (previously known as APOPT1) is involved in mitochondrial cytochrome c release and apoptosis and with a protective role under oxidative stress (Brischigliaro et al. 2019). Whereas pathogenic mutations in COA8 have been found in patients with mitochondrial complex IV deficiency (MC4DN17 [Melchionda et al. 2014]), which include systemic and neurological features, our data suggest that subtle, context-specific regulation of COA8 expression during development may also contribute to psychiatric disease risk.
Most genes have been associated with at least one regulatory eQTL. However, protein-coding mutations are relatively rare. We reasoned that genes harboring rare deleterious coding mutations might also be regulated by common variants—albeit with subtler phenotypic effects. To determine whether context-specific eQTLs can be used to connect common GWAS variants with rare disease-causing mutations, we assembled results from five large exome studies of neurological or psychiatric traits, finding 1672 genes with rare variants that are associated with at least one disease or developmental condition (Methods). Out of these 1672 genes, 225 are eGenes in at least one cell type or condition in our data, and 123 have a significant GWAS association (Fig. 4D). At a relaxed eQTL threshold (lfsr < 0.2), we identified 16 cases (corresponding to six genes) in which the lead SNP was significantly associated with a brain-related GWAS trait (Supplemental Table S6). For example, missense mutations in POU3F3 are associated with developmental delay, autism, epilepsy, and psychiatric morbidities (Rossi et al. 2023). In dividing radial glia, we identified a novel eQTL for POU3F3 (rs1451533), which is associated with brain structural features assessed by imaging. Similarly, we identified a novel neuronal eQTL for the sodium channel gene SCN1A—rare variants in which are associated with epilepsy, developmental delay, and autism—which is a risk SNP for epilepsy, suggesting that subtle changes in expression level result in similar phenotypes as loss-of-function mutations. Overall, half of the eQTLs that surfaced in this comparison were novel and two of the six eGenes were not detected under baseline oxygen conditions, highlighting the importance of examining diverse cell types and perturbed states to identify trait-relevant regulatory effects.
Discussion
In this study, we measured cell type–specific responses to oxygen stress in a genetically diverse panel of brain organoids. Oxygen perturbation induced a robust transcriptional response across all assayed cell types, which included a common oxygen stress response signature as well as cell type–specific changes. Leveraging the genetic and cellular heterogeneity present in our organoid panel, we identified thousands of dynamic, oxygen-responsive eQTLs, many of which have effects that are undetectable at baseline and therefore are absent in data collected from postmortem cortical tissues. Moreover, the use of topic modeling allowed us to identify genetic effects that transcend categorical notions of cellular identity to regulate cell division, differentiation, and other continuous processes. By collecting functional data from biological contexts that are difficult to access in vivo, we were able to characterize the putative regulatory role of hundreds of GWAS loci, many of which had never been associated with an eQTL. Our results show that repurposing organoids for gene-environment interaction studies is a valuable and increasingly tractable approach that complements population genomic studies of primary tissues and may be particularly valuable for studies of neurological and psychiatric diseases.
Our organoid differentiation approach included a small molecule cocktail designed to promote dorsal telencephalic patterning, which we selected to mitigate the between-organoid variability previously observed in unpatterned cortical organoids (Velasco et al. 2019). We nonetheless found substantial differences among organoids derived from different donor cell lines, with some cell types present in only a minority of samples. Whereas the causes of this heterogeneity could not be assigned with our experimental design, previous studies have correlated differentiation efficiency with gene expression states of starting iPSC cultures (Jerber et al. 2021). Whether iPSC states conducive to balanced organoid formation are intrinsic to individual lines or can be fully modified by in vitro conditions is an interesting question, and future efforts to develop organoid technology would benefit from considering interindividual reproducibility in addition to intraindividual consistency. In our study, this cell type heterogeneity resulted in incomplete power to map eQTLs across cell types, reflecting the trade-off between cell type resolution and abundance inherent to single-cell data analysis. More generally, the limited sample size of this study precluded the use of colocalization methods for comparing eQTLs to GWAS results, limiting the extent of our causal conclusions about context-dependent eQTLs in disease. At the same time, the heterogeneity of our in vitro system allowed us to measure gene regulation in biological contexts that have never been examined at the population level.
Although all cell types in our study responded to hypoxia, intermediate progenitor cells were among the most sensitive to hypoxic challenge. At the same time, immature neurons showed increased expression of the intermediate progenitor cell marker TBR2/EOMES following 24-h hypoxia exposure, possibly indicating a rapid transition from intermediate progenitor to neuronal identity. Paşca et al. (2019) demonstrated that a 48-h hypoxic challenge led to an unfolded protein response (UPR)-dependent depletion of TBR2+ intermediate progenitors, which underwent an apparently precocious developmental transition to CTIP2+ neurons. Using a related model, Daviaud et al. (2019) observed a greater vulnerability among outer radial glia, including a change in cell division and delayed defects in differentiation and survival. Although we observed minimal induction of transcriptional markers of UPR and found only modest changes in the overall abundance of intermediate progenitor cells, our results are consistent with these findings in showing an increased sensitivity among both intermediate progenitor cells and radial glia, suggesting that the effects of hypoxia on intermediate progenitor cell development are visible even after a short exposure period. It is also possible that a sustained 48-h shift from 20% to <1% oxygen, used by Paşca and colleagues, induces a stronger effect with greater dependence on the UPR than our paradigm, which includes a period of adaptation to physiological oxygen and a shorter hypoxic treatment. In either case, our data suggest that the transitional states from radial glia to newborn neurons are especially labile in the face of fluctuating environmental oxygen and may play an outsized role in linking transient oxygen stress to brain phenotypes (Romanko et al. 2004; Daviaud et al. 2019; Paşca et al. 2019). Whereas our analysis of the spatial organization within organoids excludes the simplest explanation for this oxygen sensitivity—namely, that cells’ proximity to the organoid surface dictates their stress sensitivity—the local tissue microenvironment likely conditions responses to oxygen- and other exogenous stress. Increasing throughput and affordability of spatial transcriptomic platforms offers exciting opportunities to take this analysis further and estimate the contributions of spatial location and local cell type composition to stress responses at the gene level.
We characterized gene regulatory variation in the context of normoxia, hypoxia, and hyperoxia. Although acute hypoxia-induced brain injury is relatively rare, transient fluctuations in brain oxygen availability are common throughout life, suggesting that the gene regulatory effects in this study may be pervasive in the population. For example, sleep apnea, in which repeated episodes of oxygen desaturation and restoration conspire to produce persistent oxidative stress, is estimated to affect nearly a billion people worldwide (Benjafield et al. 2019). Conversely, preterm birth prematurely shifts the newborn from lower in utero to higher atmospheric oxygen levels, an ambient increase that may be countered by underdeveloped lung and breathing control systems (Salmaso et al. 2014). Outside of these pathological but common scenarios, recent studies using in vivo oxygen biosensors have also observed regions of local tissue hypoxia in rodent brains, both at rest and during demanding tasks (Butt et al. 2021; Beinlich et al. 2024). Although the significance of these minute-to-minute environmental fluctuations are not yet clear, their existence raises the possibility that hypoxia-evoked transcriptional changes may be ongoing features of the brain under normal conditions. Our results imply that these conditions elicit a host of gene expression changes with varying consequences across the population and that those differences may in turn affect complex brain-related traits. Although the experimental treatments applied in our study may provide the essential context for linking regulatory variation to disease, as both hypoxia and hyperoxia can occur in vivo, such a direct connection seems unlikely to be the general rule. Rather, the cellular states induced by our in vitro treatments may share features with in vivo pathological conditions, thereby invoking similar regulatory regimes. Thus, our ability to subject inaccessible, transient cell types to disease-relevant conditions provides a useful complement to in vivo investigations linking genetic susceptibility to neurodevelopmental phenotypes.
Our study focused on just three experimental conditions but revealed that even simple in vitro manipulations produced a range of stress responses across cells in an organoid culture. By employing a topic model, which can describe cell states in greater detail than discrete categorical labels, we identified additional eQTLs and, furthermore, connected eQTLs found in a coarse pseudobulk analysis to specific cell state topics. By extending the topic modeling framework to a wider range of treatments, future studies could determine how many of the dynamic eQTLs we discovered can also be identified in the presence of other in vitro stressors that might induce similar gene expression programs. Moreover, approaches that analytically identify the major axes of gene expression variation, like the topic modeling we employed, also offer one solution to the problem of resolving cell states that lie along a continuum, which poses a challenge for standard pseudobulk methods for analyzing single-cell data. Recent technical advances, including supplementing organoids with nonneuronal cell types (Schafer et al. 2023) and in vivo implantation (Mansour et al. 2018; Revah et al. 2022; Paşca 2024), promise to dramatically expand the scope of disease-relevant interactions that can be captured in brain organoids and further extend the approach employed here. In this vein, a recent study from Antón-Bolaños et al. (2024) used brain organoids to identify cell line–specific differences in response to neurotoxic ethanol and valproic acid exposures (Antón-Bolaños et al. 2024). The response-eQTL framework employed here can advance our understanding of such individual-specific treatment responses by connecting them to germline genetic features. We expect that future studies of regulatory variation in treatment- and stress-related contexts will help prioritize targets for in vivo experimental manipulation.
Methods
Stem cell culture and organoid formation
We generated brain organoids using 21 iPSC lines (12 male, 9 female) drawn from an extensively characterized panel of iPSCs derived from Yoruba individuals from Ibadan, Nigeria (YRI) (Banovich et al. 2018). Cells were passaged at least twice before organoid formation. Organoids were formed using a protocol modified from published methods, with detailed methods provided in the Supplemental Methods (Lancaster and Knoblich 2014; Cederquist et al. 2019; Velasco et al. 2019; Yoon et al. 2019). Briefly, cells were dissociated and passed through a 40-µm filter, then aggregated by centrifugation in 96-well ultralow attachment round-bottom plates, with 10,000 cells per well in StemFlex medium with penicillin/streptomycin, 5 µM XAV939, and CEPT (Chen et al. 2021). After overnight incubation, aggregates were treated with a dual SMAD inhibition cocktail over the course of 7 days, followed by 4 days of N2-supplemented medium. Aggregates were embedded in Matrigel droplets and transferred to ultralow attachment six-well plates in 1:1 DMEM/F12:Neurobasal medium supplemented with chemically defined lipid concentrate, N2, nonessential amino acids, Glutamax, beta-mercaptoethanol, N21, insulin, and penicillin/streptomycin. After 3 weeks, organoids were gradually transitioned to BrainPhys-based medium, in which DMEM/F12/Neurobasal base medium was replaced with BrainPhys medium, in 25% increments over the course of four feedings. Organoids were maintained for four additional weeks in BrainPhys-based medium before sample collection for a total of 8 weeks of maturation. Prior to oxygen manipulation (see Supplemental Methods), cultures were maintained at 20% oxygen (room air with 5% CO2).
Single-cell RNA sequencing sample preparation
Organoids were prepared for single-cell capture using a combination of enzymatic and mechanical dissociation. Organoid medium was replaced with 1 mL papain solution (20 U/mL in EBSS, Worthington LK003150) with DNase I (100 U/mL, Worthington) supplemented with actinomycin D (5 µg/mL, Sigma-Aldrich A9415) and TTX (1 µM, Tocris 1069), and organoids were rapidly sheared with a pair of needles. Enzymatic digestion proceeded in the incubator, with continuous shaking, for 30 min. Organoids were pipetted twice with a 700- to 800-µm fire-polished Pasteur pipet, then returned to the incubator for an additional 10 min of enzymatic digestion. Samples were gently triturated four times each with fire-polished Pasteur pipets of decreasing widths (800, 600, and 300 µm) and heavy debris was allowed to settle. Samples were transferred to tubes with 2 mL inhibitor solution (3.75 mg/mL ovomucoid, 3.75 mg/mL albumin, 100 U/mL DNase I, Worthington LK003150) and spun for 5 min at 200g. Pellets were resuspended in cold Neurobasal medium with 0.5% BSA and actinomycin D and counted using a Countess II automated cell counter (Thermo Fisher Scientific) with Trypan blue. Cells from different individuals were pooled to equal concentrations, yielding three combined samples (control, low-oxygen, high-oxygen), spun down, resuspended in cold Neurobasal medium with actinomycin D, filtered with a 40-µm filter (Flowmi), and counted using a hemacytometer. Samples were loaded onto a 10x HT chip (10x Genomics) for single-cell encapsulation according to the manufacturer's instructions, targeting approximately 3000 cells per individual per treatment condition. Sequencing libraries were prepared using the 10x Genomics 3′ HT kit v3.1, according to the manufacturer's instructions, in a single batch, and libraries were sequenced according to 10x Genomics specifications, targeting a minimum of 20,000 reads per cell, on an Illumina NovaSeq 6000 instrument at the University of Chicago Genomics Core Facility (RRID:SCR_019196).
Differential expression analysis
To identify genes differentially expressed between treatment conditions, we relied on well-established methods for analyzing bulk RNA sequencing data. Single-cell transcriptomes were summed to pseudobulk samples, each of which corresponded to one combination of individual, treatment condition, cell type, and collection batch. Oligodendrocytes and midbrain dopaminergic neurons were excluded from analysis for lack of sufficient sample sizes, and pseudobulk samples derived from fewer than 20 cells were excluded. Pseudobulk data were TMM-normalized, and genes were filtered using the filterByExpr function in the edgeR package (Robinson et al. 2010), using treatment condition as the primary comparison group. A separate linear mixed model was estimated for each cell type using dream (Hoffman and Roussos 2021), with treatment condition modeled as a fixed effect and batch and parental iPSC line as random effects. Further details about assessing patterns of differential sharing across cell types are provided in the Supplemental Methods.
Topic modeling
Topic modeling offers an alternative to fixed-category classification of cell states, allowing for cells to be described quantitatively by multiple gene expression programs. To alleviate the computational burden of topic model estimation and downstream analysis, single-cell data were first aggregated into 10,707 “pseudocells” by first clustering at high resolution (resolution = 20) and then splitting each cluster of cells by parental cell line and treatment condition. The fastTopics package (Dey et al. 2017; Carbonetto et al. 2022) was used to fit models with a range of topics (k = 10–40), and the most parsimonious model that retained a clear hypoxia-associated topic was selected. Models were estimated by Poisson nonnegative matrix factorization (fit_poisson_nmf) of the pseudocell count data using 400 expectation-maximization steps, followed by 200 stochastic coordinate descent steps. A multinomial topic model was obtained using the poisson2multinom command. Individual topics were analyzed by grade-of-membership differential expression analysis (Carbonetto et al. 2023), comparing topic-specific DE results to cell type– or treatment-specific marker genes.
Cis eQTL analysis
Cis eQTLs were identified separately for each combination of cell type and oxygen treatment condition using methods originally developed for bulk RNA-seq analysis. Briefly, pseudobulk expression measurements of cells grouped by parental iPSC line, treatment, and cell type were obtained, excluding pseudobulk samples derived from fewer than 20 cells. Cell types present in fewer than eight individuals in any treatment condition after pseudobulk filtering were excluded from subsequent analysis. QTL testing was performed using MatrixEQTL (Shabalin 2012), using filtered variants within 50 kb of a gene's transcription start site, a range limited by the increasing burden of null tests using larger windows, even in larger samples collected in a similar manner (Popp et al. 2024). Gene expression principal components, obtained using the prcomp function in R (R Core Team), were used as covariates after confirming that additional sample-level factors did not increase eQTL discovery. The number of gene expression principal component covariates was chosen for each cell type and treatment so as to explain more variance in our data than in a random permutation of the data (Kumasaka et al. 2016).
Most genetic effects on gene expression are expected to be shared across conditions. To increase power to detect subtle eQTL effects, mash was used to compare the strongest variant-gene associations across treatment conditions independently within each cell type. Note that this approach does not allow us to make rigorous statements about sharing across different cell types but rather across treatment conditions within a single cell type. Posterior eQTL effects were considered shared across two (or more) conditions if the variant-gene pair was significant (i.e., local false sign rate < 0.1) in at least one of the conditions and the posterior effect size estimates differed by a factor of less than 2. Conversely, eQTL effects were considered oxygen-responsive if they were not shared in at least one oxygen condition. See Supplemental Methods for additional details about data processing for eQTL identification.
Data access
All raw and processed sequencing data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE273907. Code used to generate the results and figures in this publication is available as Supplemental Code and at GitHub (https://github.com/bumans/organoid_oxygen_eqtl).
Competing interest statement
The authors declare no competing interests.
Acknowledgments
We thank Natalia Gonzalez, Katherine Rhodes, and Bradley Wierbowski for manuscript comments and Josh Popp, Wenhe Lin, Peter Carbonetto, Benjamin Fair, and Yunqi Yang for technical advice. Jonathan Burnett and Olivia Allen assisted with cell culture and cryosectioning. Imaging was performed at the University of Chicago Integrated Light Microscopy Core (RRID: SCR_019197), and sequencing was performed by the University of Chicago Genomics Facility (RRID: SCR_019196). This work was supported by a grant from the National Institutes of Health (R35GM131726 to Y.G.).
Author contributions: B.D.U. and Y.G. conceived the study. B.D.U. performed experiments and analysis. Y.G. supervised the study. B.D.U. and Y.G. wrote the manuscript.
Footnotes
-
[Supplemental material is available for this article.]
-
Article published online before print. Article, supplemental material, and publication date are at https://www.genome.org/cgi/doi/10.1101/gr.280219.124.
-
Freely available online through the Genome Research Open Access option.
- Received November 12, 2024.
- Accepted May 30, 2025.
This article, published in Genome Research, is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.
