Characterization of human mosaic Rett syndrome brain tissue by single-nucleus RNA sequencing

Abstract

In females with X-linked genetic disorders, wild-type and mutant cells coexist within brain tissue because of X-chromosome inactivation, posing challenges for interpreting the effects of X-linked mutant alleles on gene expression. We present a single-nucleus RNA sequencing approach that resolves mosaicism by using single-nucleotide polymorphisms in genes expressed in cis with the X-linked mutation to determine which nuclei express the mutant allele even when the mutant gene is not detected. This approach enables gene expression comparisons between mutant and wild-type cells within the same individual, eliminating variability introduced by comparisons to controls with different genetic backgrounds. We apply this approach to mosaic female mouse models and humans with Rett syndrome, an X-linked neurodevelopmental disorder caused by mutations in the gene encoding the methyl-DNA-binding protein MECP2, and observe that cell-type-specific DNA methylation predicts the degree of gene upregulation in MECP2-mutant neurons. This approach can be broadly applied to study gene expression in mosaic X-linked disorders.

Figure 1.

Single-cell SNP sequencing in a female mouse model of Rett syndrome. A) Flow chart of single-cell SNP sequencing pipeline. Single-cell RNA sequencing was performed on visual cortex from five female Mecp2^+/− mice followed by graph clustering to identify the group of excitatory neurons (Slc17a7 +). Allele-specific SNPs in genes expressed in cis with the Mecp2 mutation were identified by variant calling and then used to assign the corresponding transcriptotype to the individually sequenced cells. B) Heatmap of reads per analyzed cell (rows of the heatmap) that map to wild-type (WT)- or knockout (KO)-specific SNPs (columns of the heatmap). C) Violin plots of Mecp2 mRNA counts in cells that were grouped based on their SNP-identified transcriptotype (WT, Mecp2^+/− wild-type excitatory neurons, KO, Mecp2^+/− mutant excitatory neuron, tails represent min and max of data) or by randomly assigned transcriptotypes (Random 1, Random 2). Mecp2 expression was significantly higher in the WT cells (sampled n = 593) compared to KO cells (n = 593) (Kruskal-Wallis test, H = 210, ****P < 0.0001, + indicates mean) and the populations with randomly assigned transcriptotypes (Random 1, n = 593, Random 2, n = 593; ****P < 0.0001). The groups with randomly assigned transcriptotypes had similar levels of Mecp2 expression (P > 0.9999). For the transcriptotyped excitatory neurons, we obtained an average of 7,634 transcripts per cell representing 3,879 distinct genes. D) The number of significantly misregulated genes (FDR < 0.1, monocle2) when comparing gene expression differences between groups of mutant and wild-type excitatory neurons (KO v WT, 734 genes) or two groups of randomly assigned transcriptotypes (Random, 4 genes). E) The mean fold-changes of the misregulated genes described in D (KO v WT, Random) are displayed as a function of excitatory neuron gene body DNA methylation (mCA/CA) (KO v WT, Pearson’s r = 0.38, Random, Pearson’s r = 0.04). The correlation between MeCP2-dependent gene expression and mCA/CA was significantly greater in KO v WT than Random (permutation test, P < 0.001). F) The fold-change of genes in D (KO v WT, Random) binned by gene body MeCP2 ChIP enrichment over input. The correlations between MeCP2-dependent gene expression and two MeCP2 ChIP replicates from purified cortical excitatory neurons (ChIP1, Pearson’s r = 0.41, ChIP2, Pearson’s r = 0.31) are significantly greater than the correlations observed in the Random controls (Random ChIP 1, Pearson’s r = 0.06, Random ChIP 2, Pearson’s r = 0.04) (permutation test, P < 0.001). G) Mean fold-change in gene expression of mutant excitatory neurons (KO) compared to wild-type excitatory neurons (WT) from Mecp2+/− mice, with genes separated into groups of highly methylated genes (normalized expression > 0.1, high mCA, top 25%) or lowly methylated genes (normalized expression > 0.1, low mCA, bottom 66%) and binned by their gene length. MeCP2-dependent gene expression and gene length were significantly more correlated in KO v WT than Random for high mCA genes (KO v WT, Pearson’s r = 0.10, Random, Pearson’s r =0.00, permutation test P < 0.001). The correlations between MeCP2-dependent gene expression and gene length were not statistically different between KO v WT and Random for low mCA genes (KO v WT, Pearson’s r = 0.04, Random, Pearson’s r = 0.02, permutation test P = 0.23). In E-G, the lines represent mean fold-change in expression for genes binned according to gene length (250 gene bins, 25 gene step), methylation (100 gene bins, 10 gene step), or MeCP2 enrichment (100 gene bins, 10 gene step); the ribbon displays s.e.m. of each bin.

Figure 2.

Single-nucleus SNP sequencing of human Rett brain tissue. A) Single-nucleus RNA sequencing of occipital cortex from three females with Rett syndrome. Graph clustering nuclei from the three individuals together according to their respective brain cell types. B) Flowchart for the identification and assignment of allele-specific SNPs for each Rett donor. Single nuclei suspensions from each Rett donor were sorted based on their level of immunoreactivity to a C-terminal MeCP2 antibody (MECP2^high and MECP2^low). The weak staining observed in MECP2^low nuclei represents background immunofluorescence. cDNA from the MECP2high and MECP2^low nuclei was Sanger sequenced to confirm that the sorted populations expressed the expected MECP2 allele. Deep high-throughput RNA sequencing of these populations followed by variant calling identified the allele-specific SNPs that were used to assign transcriptotypes to each nucleus from the single-nucleus RNA sequencing dataset shown in A. C) Heatmap of reads per cell (rows of the heatmap) that map to WT- or MECP2 mutant (MT)-specific SNPs (columns of the heatmap) for each of the three donors. D) The number of total nuclei, excitatory neuronal nuclei, and VIP interneuronal nuclei that could be transcriptotyped from the single-nucleus RNA sequencing dataset of Rett donors. E) The number of significantly misregulated genes (FDR < 0.01, monocle2, R255X v WT, 3158 genes in excitatory neurons, 237 genes in VIP interneurons) identified when comparing gene expression differences between groups of mutant and wild-type neurons, or two groups of neurons with randomly assigned transcriptotypes (Random, 2 genes in excitatory neurons, 10 genes in VIP interneurons). The difference in number of misregulated genes between excitatory and inhibitory neurons is largely explained by the number of cells analyzed (Supplementary Fig. 8B). The number of excitatory neuronal nuclei and VIP interneuronal nuclei used for differential expression analysis is shown in D.

Figure 3.

Cell-type-specific DNA methylation patterns predict gene misregulation in Rett syndrome. For each graph in A-B and D-E, mean fold-change in gene expression of R255X MECP2 nuclei compared to WT nuclei (R255X v WT) or of two groups of the respective cell type that were randomly assigned transcriptotypes (Random) is binned according to the fraction of gene body DNA methylation (mCH/CH). Gene expression changes (FDR < 0.01, monocle2) from R255X v WT or Random excitatory neurons are compared to patterns of DNA methylation from (A) excitatory neurons (Pearson’s r = 0.22) or from (B) VIP interneurons (Pearson’s r = −0.01) (250 gene bins, 25 gene step). R255X v WT is significantly more correlated with excitatory neuron mCH/CH than Random (A) (permutation test, P < 0.001) and significantly more correlated with excitatory mCH/CH patterns than mCH/CH patterns from VIP interneurons (B, R255X v WT (A) correlation compared to R255X v WT (B), permutation test, P < 0.001). Gene expression changes (FDR < 0.25) from R255X v WT or Random VIP interneurons are compared to DNA methylation patterns from (D) excitatory neurons (Pearson’s r = 0.05, R255X v WT; Pearson’s r = 0.03 Random) or from (E) VIP interneurons (Pearson’s r = 0.18, R255X v WT; Pearson’s r = −0.05 Random) (50 gene bins, 5 gene step). R255X v WT is significantly more correlated with mCH/CH than Random in (E) (permutation test, P < 0.001) but not in (D) (permutation test, P = 0.71). In VIP interneurons, the correlation of R255X v WT with mCH/CH is significantly greater for mCH/CH patterns from VIP interneurons (E) than mCH/CH patterns from excitatory neurons (D) (permutation test, P = 0.008). C,F) Mean fold-change in gene expression of R255X v WT excitatory neuronal nuclei (C) or VIP interneuronal nuclei (F) for expressed genes (> 0.1 normalized counts) with high mCH (top 25% mCH/CH) or low mCH (bottom 66% mCH/CH) binned according to gene length (250 gene bins, 25 gene step). MECP2-dependent gene expression and gene length were significantly more correlated in R255X v WT than Random for high mCH/CH genes (R255X v WT: C, Pearson’s r = 0.07; F, Pearson’s r = 0.08; Random: C, Pearson’s r = 0.02; F, Pearson’s r = 0.00, C, permutation test, P = 0.007, F, permutation test, P < 0.001) and significantly more anti-correlated for low mCH/CH genes (R255X v WT: C, Pearson’s r = −0.07; F, Pearson’s r = −0.09; Random: C, Pearson’s r = 0.00; F, Pearson’s r = 0.00, C, F, permutation test, P < 0.001). The lines represent mean fold-change in expression for genes binned as described; the ribbon is s.e.m. of each bin.

Figure 4.

Characterization of MECP2-regulated genes in human and mouse A) Venn diagram of the number of overlapping significantly up-regulated (left) or down-regulated (right) genes (FDR < 50.1, monocle2) between R255X MECP2 mutant and wild-type nuclei in excitatory neurons of each donor. P-values describing the significance of overlap between pairs of up- or down-regulated gene lists were calculated by hypergeometric testing. B) Boxplot of the gene body DNA methylation level (mCH/CH) of the 537 overlapping up-regulated genes and 395 overlapping down-regulated genes in the 3 donors, as well as all other expressed genes (****P < 0.0001 (Dunn’s), Kruskal-Wallis test H(2) = 146.6) C-D) Lists of the most highly significant gene ontology terms (Fisher’s Exact test with FDR) enriched in the 537 overlapping genes that are up-regulated (C) or the 395 overlapping genes that are down-regulated (D) between R255X MECP2 mutant and WT excitatory neurons. E) Venn diagram of the genes that are commonly up-regulated (top, P = 2.1 × 10⁻¹², hypergeometric test) or down-regulated (bottom, P = 1.9 × 10⁻³⁹, hypergeometric test) in mutant MECP2 compared to wild-type excitatory neurons in human and female heterozygous mice. F) Boxplot of the fraction of gene body DNA methylation (mCH/CH) of the 58 overlapping up-regulated genes and 84 overlapping down-regulated genes between human and mouse (****P < 0.0001 (Dunn’s), Kruskal-Wallis test H(2) = 52.35). Boxplots show the median (line), 25^th to 75^th percentiles (box), and 1.5X the interquartile range (whiskers).