Divergent landscapes of A-to-I editing in postmortem and living human brain

Global Alu editing is elevated in postmortem prefrontal cortex

Since most RNA editing occurs in Alu elements^3–4, we first computed an Alu editing index (AEI) for each sample. The AEI is quantified by measuring the total number of edited adenosines over all adenosines with supporting RNA-sequencing read coverage in Alu elements across the entire transcriptome and is a metric of global Alu editing activity (see Materials and Methods) (Supplementary Table 1). A significant increase in the AEI was observed in postmortem relative to living DLPFC (p=4.3×10^-75; Cohen’s d =2.88) (Figure 2A). A transcriptome-wide comparative analysis compared postmortem to living DLPFC and revealed that this shift was accompanied by heightened expression of ADAR (q-value=9.3×10^-87), ADARB1 (q-value=3.5×10^-32), and ADARB2 (q-value=2.5×10^-21) in postmortem DLPFC (Figure 2B). Notably, ADAR was the 15^th most differentially expressed gene in postmortem DLPFC and was strongly correlated with the AEI (r=0.65) (Figure S1). Next, a linear mixed model quantified the fraction of global Alu editing variance explained by known biological and technical factors. Differences between living and postmortem tissues explained the largest amount of Alu editing variability (∼72%) (Figure 2C), while minimal variance was explained by other known factors, including differences by medical diagnosis (<0.5%), brain banks (<0.5%) and estimated neuronal cell type proportions (<0.5%) (Figure S2), extended postmortem interval (PMI; <0.5%) and RNA integrity (RIN; <0.5%).

• Download figure
• Open in new tab

Figure 2.

Global Alu editing across living and postmortem DLPFC. (A) Alu editing index (AEI; y-axis) computed on bulk RNA-seq from living and postmortem (PM) DLPFC. Two-sided linear regression was used to test for significance. (B) ADAR, ADARB1 and ADARB2 normalized expression profiles on bulk including RNA-seq between living and PM. All boxplots show the medians (horizontal lines), upper and lower quartiles (inner box edges), and 1.5× the interquartile range (whiskers). Reported BH adjusted p-values were derived from a moderated t-test comparing transcriptome-wide gene expression between living and postmortem tissue. (C) Linear mixed model explaining AEI variance by eleven known biological and technical factors. (D) UMAP dimension reduction analysis of snRNA-seq classified nine unique cell populations. Values in brackets indicate the number of cells per sub-population: excitatory (EXC) and inhibitory (INT) neurons, astrocytes (AST), microglia (MG), oligodendrocytes (OLI), OLI precursor cells (OPCs), endothelial cells (ENDO). (E) The mean frequencies for each cell population quantified between living and PM. Two-sided linear regression was used to test for significance. (F) Hierarchical clustering of scaled ADAR, ADARB1 and ADARB2 expression across all cell populations. (G) Cell type-specific ADAR, ADARB1 and ADARB2 expression for living and PM. (H) Alu editing index computed for each cell population for each donor and compared across living and PM samples (bottom). PM-induced effect sizes calculated by Cohen’s d for each cell population (top). (G-H) Two-sided linear regression was used to test for significance (*denotes p< 0.05). All boxplots show the medians (horizontal lines), upper and lower quartiles (inner box edges), and 1.5× the interquartile range (whiskers). Two-sided linear regression was used to test for significance. (I) Pearson’s correlation coefficient between the mean Alu editing index and mean normalized ADARB1 expression for each cell population according to living and PM samples. Standard error bars capture group-wise variance within living and postmortem tissues, respectively. RNA-seq analysis encompassed 164 and 233 biologically independent samples from living and postmortem sources, respectively. Single-nucleus RNA-seq was conducted on 31 living and 21 PM biologically independent samples. All box plots in this figure show the medians (horizontal lines), upper and lower quartiles (inner box edges), and 1.5 × the interquartile range (whiskers).

To further explore the possible influence of differences related to PMI and RNA degradation on Alu editing in living and postmortem tissues, two supporting analyses were conducted. First, the minimal effect of extended PMI on the AEI was validated by studying 2841 independent transcriptome samples across four large-scale postmortem brain consortia (GTEx project [n=1,129]; Mount Sinai Brain Bank [n=876], PsychENCODE [n=251]; BrainSpan [n=585]). These secondary postmortem analyses confirmed weak associations between PMI and the AEI (r² = 0.006, r² = 0.03, r² = -0.007, r² = -0.019 respectively) (Figure S3), indicating that elevated global Alu editing in postmortem tissue is not likely driven by extended PMI. Second, we measured A-to-I editing dynamics using an existing molecular degradation assay of the human DLPFC³¹, whereby no significant changes in ADAR expression nor the AEI were observed throughout increasing degradation stages of the DLPFC (Figure S4). These results suggest that potential differences in PMI and RNA degradation do not fully account for the observed changes between postmortem and living DLPFC.

We next replicated features of RNA editing across 206,568 single nuclei sequenced from a non-overlapping sample of 31 living and 21 postmortem DLPFC tissues (Supplemental Data 1). Pseudo-bulk snRNA-seq pools were used to confirm heighted global Alu editing levels (p=5.5×10^-7) as well as increased mRNA expression per cell for ADAR (p=0.01), ADARB1 (p=0.0004), and ADARB2 (p=0.02) in postmortem relative to living DLPFC (Figure S5A-B) (see Materials and Methods). Unsupervised dimensionality reduction applied to all data identified nine discrete cell type clusters (Figure 2D). Fewer oligodendrocytes were detected in postmortem DLPFC (OLI, p=1.8×10^-7), while less excitatory neurons (EXC1, p=1.8×10^-7) and inhibitory neurons (INT1, p=1.5×10^-9) were detected in living DLPFC (Figure 2E). ADAR was ubiquitously expressed across all cellular populations, ADARB1 was expressed in a subset of inhibitory and excitatory neurons, and ADARB2 was uniquely expressed in oligodendrocytes and a small subset of inhibitory cells (Figure 2F). Pseudo-bulk pools were generated for each cell type per donor (Figure S5C) and used to compare the expression of ADAR enzymes and the AEI between living and postmortem DLPFC within each cellular population. Using this approach, most cell populations in postmortem DLPFC displayed significantly increased expression of ADAR, ADARB1 and ADARB2 (Figure 2G) together with an increased AEI relative to living DLPFC (Figure 2H). Notably, the largest cell type increases in global Alu editing in postmortem DLPFC occurred in microglia, endothelial cells (ENDO) and oligodendrocyte precursor cells (OPCs) (Figure 2H). These increases were concordant with increased expression of ADAR (r=0.31), ADARB1 (r=0.57) (Figure 2I) and less so ADARB2 (r=0.09) (Figure S5F).

Accurate detection of A-to-I sites in living and postmortem cortex

To study individual sites underlying these global changes, we catalogued high-confidence RNA sites using two complementary site calling techniques followed by a series of comprehensive detection-based thresholds to safe-guard against false positives (see Materials and Methods). Here, a de novo caller was used to uncover high-quality RNA editing sites not already cataloged in existing databases, together with a supervised approach applied to three large lists of known sites (Figure S6A, Supplemental Data 1). A mean of 193,195 editing sites were detected per sample in living DLPFC and 295,343 sites were detected per sample across postmortem tissues. Importantly, these sites showed hallmark characteristics of ADAR-mediated RNA editing, as the majority of sites: (1) were A-to-I sites (∼93% living, ∼95% postmortem); (2) mapped to Alu elements (∼82% living, ∼83% postmortem) (Figure S6A); (3) were predominately known sites cataloged in editing databases (∼85% living, ∼84% postmortem) (Supplemental Data 1); (4) were sites with low editing levels (20-40%; Figure S6B); and (5) commonly mapped to introns and 3′ UTRs (Figure S6C) (Supplemental Data 1). Further, while site discovery was largely correlative with sequencing depth, no major differences in library depth were observed between living and postmortem DLPFC (Figure S6D-F), suggesting it is not a driver for the observed differences. Subsequent analyses examined exclusively A-to-I sites.

Highly dynamic A-to-I editing and RNA recoding between living and postmortem DLPFC

Two approaches probed RNA editing differences between living and postmortem DLPFC. First, analysis of A-to-I sites detected at significantly different frequencies between postmortem and living DLPFC identified 28,417 sites preferentially enriched in postmortem tissues and 1,436 sites enriched in living tissues, the majority in intronic regions (Figure S7A, Supplemental Data 2). Enrichment patterns were not explained by changes in gene expression levels (Figure S7B). Second, A-to-I sites with significantly different mean editing levels were queried between living and postmortem DLPFC, focusing analyses on 54,825 A-to-I sites found across all samples in this study. Principal component analysis (PCA) distinguished living from postmortem tissues based on editing levels for these sites and PC1 was strongly correlated with differences between these two groups (r=0.92, p=9.3×10⁻¹⁶³) (Figure 3A). Statistical analysis identified differentially edited sites: 41,044 showed higher editing in postmortem tissue (“postmortem-biased”) and 1,449 showed higher editing in living tissue (“living-biased”) (Figure 3B, Supplemental Data 2). Sites annotated as postmortem-biased were edited ∼37% within the living DLPFC sample (Figure 3C). Both postmortem- and living-biased sites predominately mapped to non-coding regions, with a moderate enrichment of living-biased sites mapping to 3’UTRs and coding regions (Figure 3D). Overall, editing levels for these sites were most strongly correlated with increased expression of ADAR (mean r=0.29) in postmortem DLPFC followed by changes ADARB1 (mean r=0.18), with little effect explained by ADARB2 (mean r=0.05) (Figure 3E).

• Download figure
• Open in new tab

Figure 3.

Dynamically regulated A-to-I sites between living and postmortem DLPFC. (A) Principal component analysis of editing levels for 54,825 high-confidence sites detected across all samples in the current study. (B) Differential editing analysis compare delta editing levels (%; x-axis) and strength of significance (-log₁₀ adjusted P, y-axis) for each site between living and postmortem (PM) DLPFC. Sites are colored by novelty (i.e. detection in REDIportal) and shaped uniquely by genic region. Reported BH adjusted p-values were derived from a moderated t-test comparing transcriptome-wide A-to-I editing levels between living and postmortem tissue. (C) Frequency distribution of mean editing levels in living DLPFC (x-axis) based on PM biased sites (y-axis). (D) The fraction of genic regions for all living biased and PM biased sites. (E) Frequency distributions of Pearson’s correlation coefficients (x-axis) between the expression for ADAR, ADARB1, ADARB2 relative to editing levels for 54,825 sites. An additional analysis modeled ADAR+ADARB1-ADARB2 to capture ADAR and ADARB1 effects. The total number of sites with significant correlations are listed in the top right corner of each histogram. (F) Dynamic recoding sites: 27 living-biased recoding sites and 31 postmortem-biased recoding sites ordered by their effect size differences (y-axis, lower) and plotted alongside with the mean editing levels (y-axis middle). The strength of evolutionary conservation (phastCons) was measured for each site and the probability of being loss of function intolerant (pLI) was measured for each gene (top). (G) Frequency distributions demonstrating that living biased recoding sites are often more strongly evolutionarily conserved and map to genes with higher pLI relative to PM biased recoding sites. Mann Whitney U-test was used to test for significance. Living Brain Project data encompassed 164 and 233 biologically independent samples from living and postmortem sources, respectively.

A small fraction of differentially edited sites was cataloged as RNA recoding sites (∼0.13%, nL=L58 sites), which introduce nonsynonymous substitutions in protein-coding regions (Figure 3F). Ranking these recoding sites by their effect sizes revealed that sites with the largest changes in editing levels between living and postmortem DLPFC typically exhibited high editing levels (>30%). This also confirmed that 14 out of the 27 living-biased recoding sites were part of a collection of well-known functional sites on excitatory, inhibitory, and G-coupled protein receptors (e.g., CYFIP2, NOVA1, GABRA3, GRIA2)². Living-biased recoding sites were also more strongly evolutionary conserved (phastCons) and mapped to genes with higher probability of being loss-of-function intolerant (pLI) relative to postmortem-biased recoding sites (p=0.002, p=0.002, respectively) (Figure 3G), underscoring their physiological relevance.

Annotating dynamically regulated sites across cellular, developmental and disease scales

We next sought to annotate the cellular context as well as the developmental and disease relevance of the differentially regulated sites between living and postmortem DLPFC. Given the scarcity of known cell-specific A-to-I sites in the brain and the technical limitations of quantifying cell-specific sites from snRNA-seq data (Figure S8, see also¹³), we leveraged an independent resource of deeply sequenced neuronal and non-neuronal nuclei isolated from ten biological replicates across five postmortem cortical regions³². These data were subjected to RNA editing calling methods described above. In doing so, we generated an expanded catalogue of cell-specific RNA editing sites across the cortex (see Supplemental Note 1 for details, Figure S9-11). Subsequently, these findings, together with previously generated collections of editing sites catalogued as cell type-specific, temporally regulated across brain development as well as sites disrupted in neurological disorders, were used to annotate differentially edited sites between living and postmortem DLPFC (Supplemental Data 4).

Postmortem-biased sites were enriched for RNA editing sites catalogued as non-neuronal cell type-specific (p=8.9×10^-7), including oligodendrocytes (p=2.6×10^-9) (Figure 4A-B). Conversely, living-biased sites were significantly enriched for sites catalogued as neuronal cell type-specific (p=0.005), including GABAergic neurons (p=0.04) (Figure 4A-B). A more focused analysis dissected the cellular specificity of differentially edited RNA recoding sites between living and postmortem DLPFC (Figure 4C). This re-affirmed that living-biased recoding sites were neuronal-specific and featured several well-known functional and highly conserved recoding sites, whereas postmortem-biased recoding sites were non-neuronal and mechanistically less understood. Living-biased sites were also significantly enriched for several additional functional categories, including enrichment for: (i) A-to-I sites disrupted in postmortem brain tissues from individuals with autism spectrum disorder (ASD) (cortex p=1.1×10^-6; temporal cortex p=1.0×10^-4; cerebellum p=1.5×10^-4) and schizophrenia (SCZ) (ACC p=1.5×10^-4; DLPFC p=7.6×10^-4; DLPFC HBCC p=4.8×10^-4); (ii) A-to-I sites with precise spatiotemporal regulation across human prenatal and postnatal brain development (DLPFC p=0.02; cerebrum p=0.006; cerebellum p=0.02) (Figure 4B); and (iii) living-based sites preferentially mapped to genes enriched for postsynaptic organization and density, as well as presynaptic activity genes (Figure 4D). Considering the high levels of editing in living tissues, these findings imply a connection between A-to-I editing sites biased towards living tissues and their functional importance and endorse the value of postmortem case/control studies to explore RNA editing in brain health and disease.

• Download figure
• Open in new tab

Figure 4. Annotating dynamically regulated sites between living and postmortem DLPFC. (A)

CAMERA enrichment scores (y-axes) for three candidate sets of A-to-I editing sites along a ranked list of differentially edited sites (t-statistics; x-axis) between living and postmortem (PM) DLPFC, from highly living-biased (right; pink) to highly PM-biased (left; blue) (x-axes). Enrichment plots for non-neuronal sites (top, PM biased), sites disrupted in autism spectrum disorder (ASD) cortex (middle, living biased) and those disrupted in schizophrenia (SCZ) ACC (bottom, living biased). (B) Summary of all multiple test corrected p-values (-log₁₀) for all sets of RNA editing sites across cell types, neurodevelopmental disorders, and brain development. (A-B) CAMERA gene set enrichment p-values, quantifying the statistical significance of overrepresented A-to-I sites within the ranked living versus postmortem data. (C) Pearson’s correlation and scatterplot of delta editing rates for cell-specific recoding sites (y-axis) versus delta editing rates for living/PM differences (x-axis). Y-axis description: Fluorescence activated nuclei sorted (FANS) neurons and non-neuronal cell populations were collected from 10 postmortem donors across five cortical regions (see Supplemental Note 1). (D) SynGO synaptic enrichment (-log10 q-value) for genes harboring living-biased editing sites (top) and genes harboring postmortem-biased sites (bottom). RNA-seq analysis encompassed 164 and 233 biologically independent samples from living and postmortem sources, respectively.

Replicating postmortem-induced effects on A-to-I editing in independent transcriptomic resources

Altogether, 72,356 A-to-I sites were catalogued as either detected at different frequencies or altered in mean editing levels between postmortem and living DLPFC, here defined as ‘LIV-PM sites’. To validate these findings, we first replicated the strong postmortem-bias in editing levels for these sites leveraging pseudo-bulk snRNA-seq data from a non-overlapping sample of 31 living and 21 postmortem DLPFC samples (r=0.56) (Figure S12). Next, we asked what fraction of the LIV-PM sites identified above are commonly detected across four independent large-scale postmortem brain transcriptome consortia, which contain a diverse collection of anatomical regions, neurological disorders, and age ranges (GTEx project [n=1,129]; Mount Sinai Brain Bank [n=876], PsychENCODE [n=251]; BrainSpan [n=586]). Collectively, these data have been extensively studied for their RNA editing properties and serve as cornerstone resources for RNA editing sites in the human brain^{13,16,19,20,33} (Supplemental Data 3). Overall, we found a significant over-representation of LIV-PM sites routinely detected across these independent postmortem transcriptomic resources, whereby LIV-PM sites represented ∼15-31% of all commonly detected A-to-I sites (Figure S13). LIV-PM sites also exhibited stably high editing levels (∼40-46%), which were significantly higher (∼7-20%) than all other detected A-to-I sites (Figure S13), implicating a systematic postmortem-induced effect for these sites across diverse cohorts, anatomical regions and ages.

Profound increases of Alu editing in postmortem tissues explained by interferon activation and hypoxia

We next explored the biological processes that may best explain the profound postmortem-biases in RNA editing. Gene set variation analysis computed single-sample scores for 10,493 Gene Ontology Biological Processes for each bulk tissue RNA-seq sample, which were regressed onto the AEI to identify biological processes predictive of alterations in global Alu editing (Figure S14A). A total of 1,688 biological processes were positive predictors of global Alu editing (FDR <5%) and were broadly enriched for categories of innate immune and inflammatory responses, hypoxia, intracellular signaling, apoptosis and cellular metabolism (Figure 5A-B). Notably, biological processes that predicted the AEI were also strong predictors of living versus postmortem DLPFC (r=0.86) (Figure S14B-C). For example, the expression of genes that subserve the following biological processes stratified living from postmortem samples and were strongly predictive of changes in the AEI: ‘inositol trisphosphate metabolism’ (GO:0032957; t-statistic=14.6, q-value=1.2×10⁻³⁶), ‘desensitization of G-protein coupled receptor (GPCR) signaling pathway’ (GO:0002029; t-statistic=12.1, q-value=8.4×10⁻²⁷), ‘regulation of cell cycle’ (GO:0051726; t-statistic=11.8, q-value=1.5×10⁻²⁵), ‘response to hypoxia’ (GO:0001666; t-statistic=9.3, q-value=4.8×10⁻¹⁷) and ‘positive regulation of interferon-gamma (IFN-γ) signaling pathway’ (GO:0060355; t-statistic=8.7, q-value=7.6×10⁻¹⁶), among other innate immune responses (Figure 5D). Furthermore, the gene expression profiles underlying these biological processes were also significantly elevated in postmortem relative to living DLPFC (Figure 5C).

• Download figure
• Open in new tab

Figure 5.

Biological processes that explain differences in Alu editing. (A) Single-sample scores were generated for 10,493 Gene Ontology Biological Processes. These pathway scores were regressed onto the AEI while covarying for known biological and technical factors. Reported BH adjusted p-values and t-statistics were derived from a moderated t-test comparing pathway scores between living and postmortem tissue. The t-statistics (x-axis) for each biological process relative to strength of association with the AEI (y-axis; -log₁₀ adjusted p-value). Pathways with an absolute t-statistic > 7 and FDR adjusted p-value < 0.05 were deemed significant (blue, negative association; red, positive association). (B) A density distribution of t-statistics illustrates most pathways are positive predictors of the AEI (top) and REVIGO semantic similarity representation of the top positive 1688 pathways (bottom,). Multiple broad groupings emerge that map to intra-cellular signaling, apoptosis, hypoxia, cellular metabolism and innate immune/inflammatory responses. Colors indicate the Adjusted P-value of the enriched GO terms. The size of each bubble shows the GO terms with more significant P-values. (C) Single sample scores that represent top pathways from each cluster also predict differences between living and postmortem (PM) tissues. (D) Pearson’s correlation coefficient illustrates regressions of single-sample pathway score onto the AEI for the top six candidate pathways. All Living Brain Project data encompassed 164 and 233 biologically independent samples from living and postmortem sources, respectively. (E) Validating the interaction between interferon-γ and the AEI. Two-dimensional (2D) hiPSC-derived neural progenitor cells (NPCs; day 18) and mature neurons (day 30) treated with interferon-gamma (IFN-γ) (PMID: 32875100). Data was generated by bulk RNA-seq from n=3 biological replicates. (F) Validating the interaction between hypoxia and the AEI. Three-dimensional (3D) model of human cortical spheroids (hCS) exposed to differing degrees of hypoxia (PMID: 31061540). Data was generated by bulk RNA-seq from n=8 biological replicates. (E-F) Two-sided linear regression was used to test for significance. All boxplots in this figure show the medians (horizontal lines), upper and lower quartiles (inner box edges), and 1.5× the interquartile range (whiskers).

To functionally validate the observed associations, we quantified the AEI across two existing cellular and mechanistically related in vitro models^34,35. First, we disentangled the influence of IFN-γ signaling on global Alu editing in hiPSC-derived neural progenitor cells (NPCs) and mature neurons treated with IFN-γ³⁴. Here, IFN-γ significantly induced global Alu editing across acutely treated NPCs and chronically treated mature neurons relative to untreated conditions (p=0.0003, p=0.02, respectively) (Figure 5E). Next, we quantified the effect of hypoxia on global Alu editing in human cortical spheroids (hCS) exposed to hypoxic conditions (<1% O₂) and after 72 hours of reoxygenation versus unexposed (21% O₂)³⁵. Hypoxic conditions significantly increased global Alu editing from 24 to 48 hours (p=0.0003) and returned to baseline levels after 72 hours of reoxygenation (Figure 5F). Notably, both models revealed significant increases in ADAR1 expression profiles (Figure S15). These in vitro results support a mechanistic model, whereby IFN-γ and hypoxia induce global Alu editing in human neuronal models, perhaps explaining some of the observed increases in the AEI in postmortem DLPFC.

Mapping context-dependent edQTLs between living and postmortem DLPFC

We next sought to elucidate RNA editing quantitative trait loci (edQTLs). Paired WGS data were used to detect SNPs that could influence A-to-I editing levels for 155 living and 195 postmortem DLPFC samples. Two cis-edQTL analyses were performed: (1) A primary cis-edQTL analysis fit A-to-I editing levels to SNPs while covarying for differences between living and postmortem tissues, sex, estimated neuronal content, RNA-seq QC metrics as well as eleven surrogate variables (PEER factors); and (2) An interaction cis-edQTL analysis tested for context-dependent effects between living and postmortem tissues (see Materials and Methods). A 1Mb window (±) was defined to search for SNP-editing pairs of an editing site and identified a total of 4,858 and 2,362 unique editing sites with cis-edQTLs (eSites) from the primary and interaction analyses at FDR < 5%, respectively (Supplemental Data 5). eSites were largely non-overlapping between the two analyses and together comprised a total of 6,895 unique editing sites. Each lead SNP was located close to their associated editing site (±200 kb) (Figure 6A). eSites from the primary and interaction analyses mapped to introns and 3’UTRs with few in protein-coding sequences (Figure 6B). Furthermore, editing levels for eSites, especially those from the interaction analysis, displayed higher variability in postmortem relative to living DLPFC (Figure 6C). Most cis-edQTLs from the primary analysis explained a large percentage of editing level variability per site, explaining up to 90% editing level differences for a given site (Figure 6D). cis-edQTLs from the interaction analysis often had smaller effects, with ∼83% illustrating postmortem-specific effects (Figure 6E). Notably, ∼2% eSites identified from the interaction analyses were commonly reported across previous cis-edQTL investigations in the postmortem brain^13,33,36 (Figure S15), indicating that postmortem-related effects do not widely contribute to our current view of edQTLs in the brain.

• Download figure
• Open in new tab

Figure 6. Context-dependent cis-edQTLs between living and postmortem DLPFC.

(A) Distribution of the cis-edQTL associations evaluating the distance between eSites and common variants. The gray box indicates ±200Lkb relative to the editing site. Inset Venn diagram depicts the overlap of eSites between the primary and interaction analyses. (B) eSite discovery (y-axis) according to genic regions (x-axis) for primary (top) and interaction (bottom) cis-edQTL analyses. (C) Editing level variance within living and postmortem DLPFC parsed by sites with primary and interaction cis-edQTLs. (D) The fraction of eSites from the primary (top) and interaction (bottom) analyses that are either neuronal or non-neuronal cell type specific. (E) Two examples of primary cis-edQTLs, in which editing levels (y-axes) are associated with common genotypes (x-axes). (F) Two examples of interaction cis-edQTLs, in which different common genotypes (x-axis) are associated with differing editing levels between living and postmortem DLPFC (y-axes). (C-E) Two-sided linear regression was used to test for significance. Analyses encompassed 164 and 233 biologically independent samples from living and postmortem sources, respectively. All boxplots in this figure show the medians (horizontal lines), upper and lower quartiles (inner box edges), and 1.5× the interquartile range (whiskers).

Ethics declaration

All human subjects research was carried out under STUDY-13-00415 of the Human Research Protection Program at the Icahn School of Medicine at Mount Sinai and was approved by the Icahn School of Medicine at Mount Sinai’s Institutional Review Board. The study design and conduct complied with all relevant regulations. Research Participants in the living cohort provided informed consent for sample collection, genomic profiling, clinical data extraction from medical records and public sharing of de-identified data.

Experimental model and subject details

Living Brain Project Sample Batch Assignment, RNA processing and RNA-sequencing

To minimize batch effects, all living and postmortem samples were processed together for RNA sequencing at the Icahn School of Medicine in New York City. First, a randomization algorithm was employed to create batches for RNA extraction, cDNA library preparation, and RNA sequencing. The algorithm aimed to minimize correlations between batch assignments and living and postmortem status, ensuring an even distribution of samples across processing steps. Next, approximately 5-10 milligrams of each LIV and PM sample were used for RNA extraction. A reference piece of postmortem brain tissue was used to standardize the aliquoting process, eliminating the need to weigh individual samples and prevent thawing. To maintain a uniform temperature, a cryostat set at -20°C was used, and all equipment was treated with RNase Zap to prevent RNA degradation. After aliquoting, samples were homogenized in TRIzol, and RNA extraction was primarily carried out with the RNeasy Kit. RNA integrity was assessed, and only specimens with a RIN greater than 4.0 were sent for sequencing. Following, CDA libraries and RNA sequencing were performed at Sema4 in Stamford, CT. Libraries were prepared using the TruSeq Stranded Total RNA with Ribo-Zero Globin Kit, and sequencing was conducted on the NovaSeq 6000 System. Two S4 flow cells were used for each sequencing batch to achieve the targeted depth of 100 million paired-end reads per sample library. Sequencing was carried out in two waves, with the first wave reaching a depth of 50 million paired-end reads and the second wave, over a year later, sequencing the remaining samples to a depth of 100 million, with an additional 50 million reads generated for samples from the first wave. Further detailed materials and methods regarding sample collection and data generation for LBP samples are outlined in Liharska et al.³⁰.

RNA-sequencing data and quality control

All RNA-sequencing data were processed and mapped to the GRCh38 primary assembly with GENCODE gene annotations v30 using STAR (v.2.7.2a)⁴⁹. Picard v2.22.3 tools marked duplicate reads and gathered RNA-sequencing short read metrics and distributions (https://github.com/broadinstitute/picard/). STAR produced a coordinated-sorted mapped BAM file for each sample, which was used to study features of RNA editing. For gene expression analyses, featureCounts (v1.6.3)⁵⁰ quantified gene expression and unspecific filtering removed lowly expressed genes with less than 1 count per million in at least 10% of samples. The resulting raw counts were normalized with voomWithDreamWeights() from the variancePartition R package (v1.20.0)⁵¹ and used for downstream analyses. Linear mixed models from the R package varianceParition were also used to characterize and identify biological and/or technical drivers that may affect the observed Alu editing index (AEI). This approach quantifies the main sources of AEI variation attributable to differences in biological factors (e.g. clinical diagnoses, age, sex ethnicity) and technical factors (e.g. RIN, batch). For further details on bulk tissue RNA-sequencing data generation and pre-processing, we refer the reader to Liharska et al.³⁰.

snRNA-seq data generation and quality control

For isolating single cells from fresh brain tissue samples (n=31), we utilized the Adult Brain Dissociation Kit (Miltenyi Biotech, #130107677), which is particularly suited for processing live tissue. This kit enables the enzymatic breakdown of the tissue while maintaining the integrity and viability of the cells, which is crucial for subsequent cell-specific analyses. The process included washing the tissue, enzymatic digestion, and mechanical dissociation to ensure a high yield of viable cells. After dissociation, the cells were further purified using filtration and centrifugation steps provided in the kit’s protocol, including the removal of myelin and debris, to ensure a clean sample for downstream applications. For extracting nuclei from frozen brain samples (n=21), we chose the Minute Single Nucleus Isolation Kit (Invent Biotechnologies, #BN-020), which is designed to handle the challenges posed by frozen tissue. This kit’s protocol is optimized to protect and isolate nuclei, which are less susceptible to damage from the freeze-thaw cycle than whole cells. The kit facilitates the gentle extraction of nuclei while preserving RNA integrity, which is essential for accurate transcriptomic profiling from postmortem tissue. Subsequent gene expression profiling for both isolated cells and nuclei was carried out using the Chromium Single Cell 3’ Gene Expression platform with Next GEM reagents (10X Genomics, #CG000204 Rev D), which allows for high-throughput analysis and precise barcode assignment to individual cells or nuclei. The integrity of the synthesized cDNA was verified through meticulous quality control checks. The libraries produced were then sequenced on the NovaSeq system using the NovaSeq 6000 S2 Reagent Kit (Illumina, #20028315), adhering to a sequencing protocol that ensures comprehensive coverage. For more details on snRNA-seq data generation, we refer the reader to Vornholt et al⁵².

For all snRNA-seq data, CellRanger software (v7.0) performed genome alignment using 3’ gene expression chemistry against the GRCh38 primary assembly, generated barcode/UMI counts, and cell filtering to create mapped BAM files and feature-barcode matrices. Reads mapped to introns were incorporated into final count matrices to include both pre-mRNA and mRNA, which is representative of nuclear RNA populations. SoupX⁵³ was applied to remove contaminating ambient RNA. To optimize cell classification and reduce unwanted variance, we quality filtered, normalized, and scaled data according to Seurat’s guidelines⁵⁴. For these data, we used a set of previously implemented methods consisting of the following steps. First, a cell was excluded if the number of expressed genes was less than 200, with the number of UMI less than 200, or the percentage of mitochondria reads more than 1%. The normalization method was LogNormalize with a scale factor of 10000. The linear regression was performed using the percentage of mitochondrial reads as a variable. Second, the FindVariableFeatures(selection.method=vst, nfeatures=2000) function from Seurat R package (v4.2.0)⁵⁴ was used to identify highly variable gene features. Third, principal component analysis and uniform manifold approximation and projection (UMAP) performed unsupervised clustering of each cell type. Hierarchical clustering was manually checked along with the top-ranked genes in each cell cluster to determine cellular specificity based on well-known gene markers to verify the assignment of cell types and subtypes. Finally, cells with inconsistent assignments were pooled into their corresponding cell type cluster based on shared transcriptome-wide profiles using CellSelector() from Seurat R package (v4.2.0). Ultimately, we report on nine major cell type clusters, defined by canonical cell-type markers.

Generating an Alu editing index

For RNA-seq data, the AEI method v1.0 computed the Alu editing index (AEI)²⁹ using a STAR mapped BAM file as input. The AEI is computed as the ratio of edited reads (A-to-G mismatches) over the total coverage of adenosines in Alu elements and is a robust measure that retains the full Alu editing signal, including editing events residing in low-coverage regions with a low false discovery rate. The resulting metric is multiplied by 100, so the index describes the percentage level of editing. The predetermined genomic regions were set to all SINE/Alu repeats using the Alu bed table of the UCSC genome browser, where most A-to-I editing occurs in mouse. Common genetic variation was also discarded using coordinates from UCSC genome browser (hg38 CommonGenomicSNPs150). Notably, we have applied these metrics to hundreds of independent samples and this method has proven to be highly robust and scalable across postmortem brain RNA-seq samples from unique studies and library preparation protocols^29,16.

For snRNA-sequencing data, the AEI method v1.0 was applied two different ways. First, the AEI was quantified ignoring all cell barcoded information, thereby analyzing each mapped snRNA-seq BAM file as a pseudo-bulk tissue. Second, using the nine cell type annotations identified via unsupervised dimensionality reduction, cell-specific barcodes were used to parse each pseudo-bulk snRNA-seq BAM file into cell type specific BAM file (i.e., pseudo-bulk pooling of cell types per donor). The AEI was then computed for each cell type-specific mapped BAM file for each donor.

Cellular deconvolution of bulk tissue

Cell type deconvolution of the bulk tissue RNA-sequencing data was performed on raw gene count data using dtangle⁵⁵ and a scRNA-seq cell type reference panel from the DLPFC⁵⁶, including GABAergic and glutamatergic neurons, oligodendrocytes, astrocytes, and microglia. The sum of GABAergic and glutamatergic neuronal cell type proportions was used as proxy of total neuronal fraction.

RNA editing site detection and annotation

A two-step approach was used to quantify high-quality A-to-I sites from sorted mapped bam files:

• 1. We first quantified A-to-I sites de novo to facilitate the discovery of A-to-I sites not yet catalogued in existing RNA editing databases. Here, JACUSA2⁵⁷ was applied with the following parameters: -p 10 –a D, M, Y, E, -m 20. All analyses considered read strandedness when appropriate.

• 2. Next, we applied a supervised approach to query nucleotide coordinates for A-to-I sites already catalogued through REDIportal⁵⁸, A-to-I sites cataloged across human brain cell types¹³ and an extensive list of A-to-I recoding sites⁵⁹. Here, the samtools mpileup function was used to query editing levels of known sites, as shown prior^13,16,19. This secondary supervised approach was applied to ensure identification and inclusion of well-known sites into downstream analyses.

Subsequently, filtering steps were applied to retain only high-quality, high-confident bona fide A-to-I sites¹⁶. Briefly, the following sites were removed: i) multi-allelic events; ii) sites mapping to homopolymeric regions or black-listed genomic regions in the genome⁶⁰; iii) sites mapping to common genomic variation in dbSNP(v150) and those in gnomAD with minor allele frequency greater than 0.05; iv) sites mapping to high confidence heterozygous or homozygous genomic calls using paired WGS data; v) de novo called sites adjacent to read ends and splice sites; vi) de novo called sites with coverage was below ten reads, edited read coverage was below three reads and an editing ratio below 1%; vii) supervised sites with coverage below five reads and the number of edited reads below three. Following, remaining sites were annotated using ANNOVAR⁶¹ to gene symbols using RefGene, repeat regions using RepeatMasker v4.1.1, known RNA editing sites using the most recent version of REDIportal and conservation metrics were gathered using phastCons from the PHAST package⁶².

Detecting dynamically regulated A-to-I sites

Two different approaches were used to identify dynamically regulated A-to-I sites. First, sites observed at significantly different population frequencies were computed using a two-proportions z-test via the prop.test() function in R. This result returns the value of the Pearson’s chi-squared test statistic, a p-value, 95% confidence intervals and an estimated probability of success. We required that all significant sites must have more than 20% difference in detection rates between living and postmortem DLPFC and no more than 10% detection levels in the comparison group for which the A-to-I site is depleted. Second, to identify sites with significantly different mean editing levels between living and postmortem DLPFC, linear modeling via the limma R package⁶³ was implemented and adjusted for the possible influence of the following covariates: RNA editing levels ∼Neurons + Sex + Age + RIN + 3’ read bias + percent mRNA bases + median insert sizes + strand balance + batch. The duplicateCorrelation() function was used to model donor (i.e. technical replicates) as a repeated measure. Additional models were fit covarying for the influence of ADAR, ADARB1 and/or ADARB2. All significance values were adjusted for multiple testing using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). Sites passing a multiple test corrected p-value < 0.05 were labeled significant. This approach was applied to all bulk RNA-seq and pseudo-bulk snRNA-seq data from the Living Brain Project and FANS-derived neuronal and non-neuronal nuclei.

Quantifying A-to-I sites across large-scale postmortem brain consortia

To query A-to-I sites from existing large-scale postmortem transcriptomic resources outlined above (GTEx, MSBB, PsychENCODE, BrainSpan), raw FASTQ files were processed and mapped to the GRCh38 primary assembly with GENCODE gene annotations v30 using STAR (v.2.7.2a)⁴⁹. Picard v2.22.3 tools marked duplicate reads and gathered RNA-sequencing short read metrics and distributions. STAR produced a coordinated-sorted mapped BAM file for each sample, which was used to study features of RNA editing. Next, because these cohorts have been extensively studied for the RNA editing properties, we applied a supervised approach (samtools mpileup) to quantify millions of known sites from three main resources: REDIportal⁵⁸, A-to-I sites cataloged across human brain cell types¹³ and an extensive list of A-to-I recoding sites⁵⁹. The supervised also required that all A-to-I sites must have a coverage of at least 5 reads and 2 edited reads to be considered for downstream analyses.

Quantifying A-to-I sites in postmortem neuronal and non-neuronal nuclei

It is acknowledged that individual A-to-I site identification using snRNA-seq presents difficulties due to limitations such as low capture efficiency and sequencing depth, resulting in partial coverage of the genome⁶⁴. Our prior research has discussed these technical challenges in depth, particularly in the context of human cortex snRNA-seq data¹³. Furthermore, snRNA-seq data are inclined to exhibit an overrepresentation of intronic editing due to the predominance of nuclear RNA, which contrasts with bulk tissue sequencing that includes both nuclear and cytoplasmic fractions^65–67. To this end, we incorporated deep RNA-seq from florescence activated nuclei sorted (FANS)-derived neuronal and non-neuronal cortical nuclei from a prior study³², which were mapped to the GRCh38 primary assembly with GENCODE gene annotations v30 using STAR (v.2.7.2a)⁴⁹. This generated a coordinated-sorted mapped BAM file for downstream analyses. Next, we applied both a de novo caller (JACUSA2) together with a supervised approach (samtools mpileup), as described above for all the Living Brain Project bulk RNA-sequencing data. Because of the deep level of sequencing performed on these cell populations, we applied the de novo caller to ensure capture of cell-type specific A-to-I sites that are not already catalouged in existing databases. All subsequent filters and thresholds as described for the Living Brain Project bulk RNA-sequencing data were also applied here.

Enrichment for cellular, developmental and disease-related sites

Correlation adjusted mean rank (CAMERA) gene set enrichment⁶⁸ was performed using the resulting sets of differential editing summary statistics between living and postmortem tissues. Here we used CAMERA to perform a competitive editing set (i.e. a set of curated A-to-I sites) rank test to assess whether the sites in each editing set were highly ranked in terms of differential editing relative to sites that are not in the editing set. For example, CAMERA first ranks editing level differences in living cortical tissues relative to postmortem cortical tissues. Next, CAMERA tests whether the user-defined editing sets are over-represented toward the extreme ends of this ranked list. After adjusting the variance of the resulting editing set test statistic by a variance inflation factor that depends on the site-wise correlation (which we set to default parameters, 0.01) and the size of the set, a p-value is returned and adjusted for multiple testing. We used this function to test for enrichment of editing sets derived from: three major brain cell types¹³, human brain development¹⁶, schizophrenia¹⁹, Fragile X Syndrome and ASD²⁰.

Single sample pathway activation scores

The gene set variation analysis (GSVA) R package⁶⁹ was applied to VOOM normalized gene expression data to generate gene set-centric activation scores for each transcriptome sample, converting a matrix of genes to gene sets. Gene set activation scores were generated across a well-curated list of 10,493 Gene Ontology (GO) Biological Processes. Subsequently, linear modeling via the limma R package⁶³ tested single-sample pathway activation scores for associations with the AEI, as well as for differences in these scores between living and postmortem samples while adjusting both analyses for the covariates described above. Gene sets passing a multiple test-corrected p-value < 0.05 with an absolute t-statistic > 7 were labeled significant. Gene sets annotated as positive and negative predictors of global Alu editing were subjected to REVIGO semantic similarity⁷⁰ to reveal consensus groups of gene sets with similar gene content.

Identification of RNA editing quantitative trait loci

Cis-edQTLs were identified for all high-quality common variants within 1 Mb (±) of an editing site using the fastQTL permutation-based analysis⁷¹. Two different analyses were run: (1) A primary edQTL analysis testing relationships between editing levels and SNPs while covarying for differences between living and postmortem tissues, sex, estimated neuronal content, the top eleven PEER factors and a series of RNA-sequencing metrics (median 3’ bias, percent mRNA bases, batch, median insert size, strand balance). This analysis was run using a total of 10,000 permutations: (2) An interaction analysis testing context-dependent effects between living and postmortem tissues while adjusting for the covariates above. This analysis was run using a total of 1,000 permutations. For the results of each analysis, all SNP-variant pairs with p-value < 0.05 obtained via beta approximation were deemed significant and used for downstream analyses.