Single-cell analyses reveal increased gene expression variability in human neurodevelopmental conditions

Abstract

Interindividual variation in phenotypic penetrance and severity is found in many neurodevelopmental conditions, although the underlying mechanisms remain largely unresolved. Within individuals, homogeneous cell types (i.e., genetically identical and in similar environments) can differ in molecule abundance. Here, we investigate the hypothesis that neurodevelopmental conditions can drive increased variability in gene expression, not just differential gene expression. Leveraging independent single-cell and single-nucleus RNA sequencing datasets derived from human brain-relevant cell and tissue types, we identify a significant increase in gene expression variability driven by the autosomal aneuploidy trisomy 21 (T21) as well as autism-associated chromodomain helicase DNA binding protein 8 (CHD8) haploinsufficiency. Our analyses are consistent with a global and, in part, stochastic increase in variability, which is uncoupled from changes in transcript abundance. Highly variable genes tend to be cell-type specific with modest enrichment for repressive H3K27me3, while least variable genes are more likely to be constrained and associated with active histone marks. Our results indicate that human neurodevelopmental conditions can drive increased gene expression variability in brain cell types, with the potential to contribute to diverse phenotypic outcomes. These findings also provide a scaffold for understanding variability in disease, essential for deeper insights into genotype-phenotype relationships.

Keywords: CHD8; NPCs; autism; autosomal aneuploidy; gene expression; neural development; transcription; trisomy 21; variability.

Graphical abstract

Figure 1

T21 drives increased gene expression variability in human NPCs (A) Schematic of isogenic euploid and T21 NPCs used for scRNA-seq analysis. (B) UMAP clustering of euploid and T21 NPCs (n = 18,262 cells) identifies six clusters: progenitor 1, progenitor 2, cycling progenitor, neuronal progenitor, gliogenic progenitor, and unknown. (C) Cell type variability analysis between T21 and euploid NPCs for each of the six clusters. -log(p_adj) is shown on the x axis, with the gray shaded area indicating a lack of significance. More variable in T21 is indicated in red. The percentage of cells in each cluster for each genotype is shown on the right. Significance was determined by the two-sided Wilcoxon rank-sum test followed by Benjamini-Hochberg correction for multiple comparisons. (D) GO analysis of HVGs and LVGs using the total list of genes present in three uncommitted progenitor clusters from euploid (left) and T21 (right) datasets (progenitor 1 [T21: 3,558 HVGs and 10,092 LVGs; euploid: 4,641 HVGs and 7,661 LVGs], cycling progenitor [T21: 2,519 HVGs and 12,713 LVGs; euploid: 2,432 HVGs and 13,928 LVGs], and progenitor 2 [T21: 3,366 HVGs and 9,256 LVGs; euploid: 2,252 HVGs and 13,496 LVGs] clusters). The color scale indicates p_adj values, and the circles are scaled to the gene ratio. (E) Constraint analysis performed using gnomAD (v.4.1.0) to identify the LOEUF (loss-of-function observed/expected upper bound fraction) for the top 150 LVGs and HVGs from progenitor 2 cluster with a 0.7 mean-normalized expression cutoff (p = 0.0016). Of 300 genes in total, 264 had constraint data (HVG n = 121 and LVG n = 143). Significance was determined using the two-sided Wilcoxon test. (F) Correlation between differential gene expression (log₂FC) and differential variability (residual T21 − residual euploid) for the progenitor 2 cluster. Significant genes are shown in red. The Pearson correlation was 0.09. All statistical significance was determined using an alpha threshold of 0.05. The schematic was created in BioRender.

Figure 2

Decreased correlation strength in the T21 highly variable gene network points to intrinsic variability in human NPCs (A) Matrices showing changes in gene-gene correlation strength for HVGs in the uncommitted progenitor clusters. Gene-gene correlations for each genotype were performed for HVGs in common between T21 and euploid conditions. Left, change in correlation strength between T21 and euploid conditions for progenitor 1 (top; mean change = 0.0017; n = 2,731 gene pairs), cycling progenitor (middle; mean change −0.0272; n = 1,592 gene pairs), and progenitor 2 (bottom; mean change = −0.0279; n = 1,442 gene pairs). Right, plotting the gene-gene pairs against the change in correlation shows the percentage of genes with increased or decreased correlation strength in T21 compared to euploid. Blue shading indicates loss in correlation strength in T21, with red shading illustrating increased correlation strength in T21. (B) Violin plots showing expression levels of example HVGs from euploid (blue) and T21 (red) datasets including RMST, LDHA, EFNA5, and TRPM3. Log₂FC and p values of differential expression analysis are displayed below each violin plot. (C) Violin plots showing expression levels of example LVGs from euploid (blue) and T21 (red) datasets including HDAC3, ATP5MC2, COX7C, and TBCA. Log₂FC and p values of differential expression analysis are displayed below each violin plot. All statistical significance was determined using an alpha threshold of 0.05.

Figure 3

Validation of T21-driven gene expression variability in human NPCs and postmortem brain tissue (A) Cell type variability analysis between T21 and euploid NPCs for each of seven clusters identified in Qiu et al.. -log(p_adj) is shown on the x axis, with the gray shaded area indicating a lack of significance. More variable in T21 is indicated in red, and more variable in euploid is indicated in blue. The percentage of cells in each cluster for each genotype is shown on the right. Significance was determined by the two-sided Wilcoxon rank-sum test followed by Benjamini-Hochberg correction for multiple comparisons. (B) GO analysis of HVGs and LVGs from euploid (left) and T21 (right) datasets from Qiu et al. The color scale indicates p_adj values, and the circles are scaled to the gene ratio. (C) Constraint analysis performed using gnomAD (v.4.1.0) to identify the LOEUF for the top 150 LVGs and HVGs from the Qiu et al. dataset (p < 0.0001). Of the total 300 genes, 274 genes had constraint data (HVG n = 140, LVG n = 134). Significance was determined using the two-sided Wilcoxon test. (D) Overlap analysis of LVGs (left) and HVGs (right) identified in progenitor clusters between our NPC dataset (Figure 1) and the Qiu et al. NPC dataset. (E) Cell type variability analysis between T21 and euploid human brain tissue for each of seventeen clusters identified in Palmer et al. -log(p_adj) is shown on the x axis, with the gray shaded area indicating a lack of significance. More variable clusters in T21 are indicated in red, and more variable clusters in euploid are indicated in blue. The percentage of cells in each cluster for each genotype is shown on the right. Significance was determined by the two-sided Wilcoxon rank-sum test followed by Benjamini-Hochberg correction for multiple comparisons. (F) GO analysis of HVGs and LVGs from euploid (left) and T21 (right) datasets from Palmer et al. The color scale indicates p_adj values, and the circles are scaled to the gene ratio. (G) Constraint analysis performed using gnomAD (v.4.1.0) to identify the LOEUF for the top 150 LVGs and HVGs from the Palmer et al. dataset (p < 0.0001). Of the total 300 genes, 260 genes had constraint data (HVG n = 113, LVG n = 147). Significance was determined using the two-sided Wilcoxon test. (H) Violin plots showing expression levels of an example HVG (CLU; left) and LVG (MACF1; right) from euploid and T21 astrocyte 1 cluster from Palmer et al. Log₂FC and p values of differential expression analysis are displayed below each violin plot. All statistical significance was determined using an alpha threshold of 0.05. The NPC and brain schematics were created in BioRender.

Figure 4

Autism-associated CHD8 haploinsufficiency drives increased gene expression variability in human brain organoids (A) Bar graph showing the number of variable clusters from 3-month organoids generated in independent genetic backgrounds (HUES66, H1, GM08330, and Mito210) and with different genotypes (CHD8, KMT5B, and ARID1B) from Paulsen et al. (B) Cell type variability analysis from an example CHD8 mutant organoid (CHD8_GM_3m_r1) from Paulsen et al. for each of 20 clusters. -log(p_adj) is shown on the x axis, with the gray shaded area indicating a lack of significance. More variable clusters in the CHD8 mutant condition are indicated in red, and more variable clusters in the isogenic WT control are indicated in blue. The percentage of cells in each cluster for each genotype is shown on the right. Significance was determined by the two-sided Wilcoxon rank-sum test followed by Benjamini-Hochberg correction for multiple comparisons. (C) GO analysis of CHD8 mutant (left) and WT control (right) from an example CHD8 mutant organoid (CHD8_GM_3m_r1) from Paulsen et al. The color scale indicates p_adj values, and the circles are scaled to the gene ratio. (D) Violin plots showing expression levels of an example HVG (FABP7; top) and LVG (BASP1; bottom) from CHD8 mutant and control CPNs from CHD8 mutant organoid (CHD8_GM_3m_r1) from Paulsen et al. Log₂FC and p values of differential expression analysis are displayed below each violin plot. (E) Bar graph showing the number of variable clusters from CHD8 mutant organoids at 20, 60, and 120 days from Villa et al. (F) Cell type variability analysis of the 120-day time point CHD8 mutant organoid from Villa et al. for each of 15 clusters. -log(p_adj) is shown on the x axis, with the gray shaded area indicating a lack of significance. More variable clusters in the CHD8 mutant condition are indicated in red, and more variable clusters in WT control are indicated in blue. The percentage of cells in each cluster for each genotype is shown on the right. Significance was determined by the two-sided Wilcoxon rank-sum test followed by Benjamini-Hochberg correction for multiple comparisons. (G) GO analysis of CHD8 mutant (left) and WT control (right) from the 120-day time point CHD8 mutant organoid from Villa et al. The color scale indicates p_adj values, and the circles are scaled to the gene ratio. (H) Overlap analysis between LVGs (top) and HVGs (bottom) from Paulsen et al. CPNs from CHD8_GM_3m_r1 and Villa et al. ENEs from 120-day organoids. All statistical significance was determined using an alpha threshold of 0.05. The organoid schematics were created in BioRender.

Figure 5

Molecular basis of gene expression variability (A) Heatmap of enriched histone post-translational modifications (PTMs) for HVGs and LVGs from the most variable cluster from each dataset (Upadhya et al. [this study], Qiu et al., Palmer et al., Paulsen et al.CHD8 mutant organoid [CHD8_GM_3m_r1], and Villa et al. 120-day CHD8 mutant organoid). The left heatmap was made with the top 150 genes based on gene residuals that met the normalized mean expression cutoff of 0.70, and the right heatmap was made with the top 150 genes based on gene residuals that met the normalized mean expression cutoff of 0.25. Data from the Epigenomics Roadmap ChIP-seq data through Enrichr were used to evaluate enrichment. The color scale indicates the -log(p_adj) level of gene set enrichment. (B) Heatmap of enriched transcription factor motifs with HVGs and LVGs from the most variable cluster from each dataset (Upadhya et al. [this study], Qiu et al., Palmer et al., Paulsen et al.CHD8 mutant organoid [CHD8_GM_3m_r1], and Villa et al. 120-day CHD8 mutant organoid). Data from the ENCODE TF ChIP-seq dataset through Enrichr were used to evaluate enrichment of the top 150 HVGs and LVGs that met the mean-normalized mean expression cutoff of 0.70. The color scale indicates the -log(p_adj) level of gene set enrichment. (C) Bar graphs showing examples of individual transcriptional regulators significantly enriched in LVGs and/or HVGs. Motifs with a minimum adjusted p value for enrichment <1e−10 are displayed. LVGs are shown in blue and HVGs in red, with -log(p_adj) on the y axis. All statistical significance was determined using an alpha threshold of 0.05.