Brain-wide alterations revealed by spatial transcriptomics and proteomics in COVID-19 infection

Abstract

Understanding the pathophysiology of neurological symptoms observed after severe acute respiratory syndrome coronavirus 2 (SARS-CoV2) infection is essential to optimizing outcomes and therapeutics. To date, small sample sizes and narrow molecular profiling have limited the generalizability of findings. In this study, we profiled multiple cortical and subcortical regions in postmortem brains of patients with coronavirus disease 2019 (COVID-19) and controls with matched pulmonary pathology (total n = 42) using spatial transcriptomics, bulk gene expression and proteomics. We observed a multi-regional antiviral response without direct active SARS-CoV2 infection. We identified dysregulation of mitochondrial and synaptic pathways in deep-layer excitatory neurons and upregulation of neuroinflammation in glia, consistent across both mRNA and protein. Remarkably, these alterations overlapped substantially with changes in age-related neurodegenerative diseases, including Parkinson's disease and Alzheimer's disease. Our work, combining multiple experimental and analytical methods, demonstrates the brain-wide impact of severe acute/subacute COVID-19, involving both cortical and subcortical regions, shedding light on potential therapeutic targets within pathways typically associated with pathological aging and neurodegeneration.

Extended Data Fig. 1 |. Gene expression profile from the frontal grey matter correlates with nCounter bulk tissue direct gene expression.

a) Bar plots showing the number of differentially expressed molecules in different assays (BH-adjusted p < 0.1, MSstatsTMT for proteomics; DESeq2 for nCounter direct gene expression; lme4 for spatial transcriptomics). b) Venn diagrams and Pearson’s correlation of overlapped differentially expressed genes from nCounter direct gene expression of the frontal lobe and spatial transcriptomics of frontal grey matter (p = 1.60E-11, GeneOverlap, Pearson’s correlation, two-sided, r = 0.77, p = 1.62E-33) and frontal white matter (p = 0.60, GeneOverlap, Pearson’s correlation, two-sided, r = −0.42, p = 0.72). c) 3 SARS-CoV2 transcripts in nCounter direct gene expression (basal ganglia) indicated there was no increased expression of viral mRNA (N = 6 for controls, N = 6 for COVID-19 patients, two-tailed unpaired t test). d) Bar plots showing the results of cell type enrichment analysis (EWCE) of NRGN neuron from Yang, et al, 2021 and its celltype markers’ enrichment when comparing to the major cell types (level 1) and subclass cell types (level 2) in Wamsley, et al, 2024 and Bakken, et al, 2021. The red dashed line represents a threshold of 2.5 for the sd_from_mean. Cell types with a sd_from_mean value > 2.5 and a BH-adjusted p < 0.05 are considered statistically significant.

Extended Data Fig. 2 |. Increased ACE2 and FURIN protein expression in the frontal cortex of severe acute/subacute COVID-19 patients.

Anti-ACE2 (a) and Anti-FURIN (b) immunohistochemistry highlight the increased protein expression at the vascular/perivascular structures (magenta arrowhead) in the white matter of the COVID-19 patients (N = 9) comparing to the non-COVID controls (N = 8), one-tailed unpaired t-test of area percentage for positive immunoreactivity, p = 3.95E-2 (ACE2) and p = 3.00E-3 (FURIN).

Extended Data Fig. 3 |. Increased NRP1 protein expression in the frontal cortex of severe acute/subacute COVID-19 patients.

a) Anti-NRP1 immunohistochemistry in the frontal cortex of non-COVID controls (N = 8) and COVID-19 patients (N = 9) are displayed. Within the non-COVID controls, NRP1 immunoreactivity has a predominant astrocytic expression pattern (cyan arrowhead) in five out of eight cases, while a predominant neuronal expression pattern (purple arrowhead) is observed in the remaining three. In the COVID-19 patients, NRP1 immunoreactivity is seen more frequently in both neurons and astrocytes. Moreover, increased neuronal immunoreactivity is found in the nucleus of five out of nine patients, and in the cytosol of two out of nine patients. The overall immunopositive area percentage is increased in the COVID-19 patients (p = 3.20E-3, one-tailed unpaired t-test). b) Anti-Fibrin immunohistochemistry in the frontal white matter of non-COVID controls (N = 5) and COVID-19 patients (N = 4). Increased intravascular and perivascular fibrin deposition (yellow arrowhead), indicative of thrombi formation and blood brain barrier/neurovascular unit defect, is seen in both non-COVID controls and COVID-19 patients (one-tailed unpaired t-test of area percentage for positive immunoreactivity, p = 2.68E-1).

Extended Data Fig. 4 |. Network dendrogram, EWCE results of WGCNA and GO enrichment analysis of hdWGCNA.

a) Consensus WGCNA network dendrogram from co-expression topological overlap of genes in frontal grey matter and pontine nuclei spatial transcriptomics. b-d) Dot plots illustrating cell type enrichment analysis (EWCE) results of consensus modules, using major cell types from Wamsley, et al, 2024(b), subclass cell types from Bakken, et al, 2021(c), and Yang, et al, 2021(d). The color intensity corresponds to the sd_from_mean of each cell type. The greyscale of the dot edge corresponds to the FDR-adjusted p value for each cell type. The dot size in each row represents the relative abundance across all the modules. e) Top GO enrichment terms of all the excitatory neuron modules (EnrichR).

Extended Data Fig. 5 |. Additional GO enrichment, PPI and MCODE network analyses of differentially expressed proteins in eight regions.

a) The graph reflects the same networks depicted in Fig. 6b (downregulated GO cluster, left panel) and Fig. 6c (upregulated GO cluster, right panel), but with the regional information displayed as a pie chart for each gene. Each GO term is proportioned within the pie chart and color-coded based on the number of proteins under the term from the corresponding region (key below). Most of the terms displayed are shared by more than three regions. b) Merged protein-protein interaction and MCODE network analyses of differentially expressed proteins in eight regions. Top GO terms for each MCODE network are selected for labeling.

Extended Data Fig. 6 |. Brain-wide differential gene expression in severe acute/subacute COVID-19 patients shows minimal overlap with hypoxia-related changes.

a) Volcano plots depict the NES values from Gene Set Enrichment Analysis (GSEA^,, GeneSet Hallmark) and FDR-q values across 10 areas of the frontal lobe, pons, and midbrain from spatial transcriptomics results. Notably, hypoxia-related genes are not significantly enriched in any of the profiled areas. b) Venn diagrams illustrate the results of comparison between the top 500 differentially up-and down-regulated genes observed in primary cortical neurons under hypoxia conditions and those identified in the grey matter of COVID-19 patients (BH-adjusted p < 0.05) in this study, no significant overlap are identified (for upregulated genes, the comparisons between Zhang, et al and frontal grey matter, pontine nuclei, and substantia nigra, p = 9.04E-1, 8.13E-1, and 1.78E-1, respectively; for downregulated genes, the corresponding p = 9.35E-1, 1.00, 3.71E-1, GeneOverlap). We also compared the 500 upregulated genes of mouse brains under hypoxic conditions for 2 days (gene pattern 1) and 7 days (gene pattern 2) to the upregulated genes in frontal grey matter, finding no significant overlap (p = 4.70E-1 and 1.88E-1, GeneOverlap).

Extended Data Fig. 7 |. Increased amyloid deposition including senile plaques and cerebral amyloid angiopathy, not Tau deposition, in severe acute/subacute COVID-19 patients.

a) Anti-Aβ42 immunohistochemistry highlights the amyloid deposition including senile plaques (orange arrowhead) and cerebral amyloid angiopathy (CAA, purple arrowhead) in the COVID-19 patients (N = 9) which are not seen in the non-COVID controls (N = 8) (one-tailed unpaired t-test of area percentage for positive immunoreactivity, p = 4.61E-2). b) Anti-phospho-Tau immunohistochemistry does not detect significant positive immunoreactivity and difference in positive immunoreactivity area between non-COVID control (N = 8) and COVID-19 patient (N = 9) (one-tailed unpaired t-test, p = 4.84E-1) groups. Rare neurons containing neurofibrillary tangles and neuropil threads (cyan arrowhead) are observed. As positive controls, Aβ42 immunopositivity highlights the characteristic senile plaques and cerebral amyloid angiopathy of advanced stage Alzheimer’s disease (Braak stage V-VI). Anti-phospho-Tau antibody highlights abundant neurofibrillary tangles and neuropil threads (cyan arrowhead), and neuritic plaques (yellow arrowhead) in AD patient. Anti-Aβ42 and phospho-Tau IHCs were repeated once.

Fig. 1 |. Upregulation of antiviral defense genes and dysregulation of cell type markers in the frontal lobe and pons of patients with severe acute/subacute COVID-19.

a, Schematic workflow and experimental design of this study. b,c, Multiplex immunofluorescence images of frontal lobe (b) and pons (c) in spatial transcriptomics (magenta: anti-NeuN; greenyellow: anti-GFAP; blue: nuclei; cyan: anti-IBA1) (biologically independent samples n = 11 for frontal lobe, n = 10 for pons). White squares highlight the ROI selected in the frontal gray and white matter (n = 121) and in the superior pons (pontine nuclei, pontocerebellar tracts, corticospinal tract and tegmentum, n = 119). Scale bar in whole tissue section, 1 mm. Scale bar in ROI, 50 μm. d, Venn diagram showing significant overlapping of DE genes (BH-adjusted P < 0.1, LMM, lme4 (ref. 33)) in the frontal gray matter and pontine nuclei (P = 2.07 × 10⁻⁷¹, GeneOverlap; Pearson’s correlation: r = 0.81, P = 0.00 × 10⁰, two-sided). e, Dot plot illustrating differential expression of SARS-CoV2 viral probes and viral-related genes in spatial transcriptomics in frontal lobe and pons. Yellow represents upregulation; blue represents downregulation. The color intensity represents the effect size (log₂FC). The grayscale of the dot edge corresponds to the adjusted P value. The dot size represents the relative abundance across six areas. f, Dot plot depicting the EWCE results of DE genes in the frontal lobe and pons, using the cortical snRNA-seq data from Wamsley et al.. Deep, deep layer; ExNeu or EXT, excitatory neuron; ODC, oligodendrocyte; BBB, blood brain barrier. g, Dot plots depicting the EWCE results of DE genes in the frontal gray and white matter using the cortical snRNA-seq data of Yang et al. and Bakken et al.. CT, corticothalamic cell; Endo, endothelial cell; ET, extratelencephalic cell; IT, intratelencephalic cell; NP, near-projecting; Oligo, oligodendrocyte; Sst, somatostatin; Vip, vasointestinal peptide; VLMC, vascular and leptomeningeal cells. In f and g, red and blue represent enrichment of the upregulated and downregulated genes. The color intensity represents the sd_from_mean of each cell type. The grayscale of the dot edge corresponds to the FDR-adjusted P value. The dot size represents the relative abundance across the assessed areas. FC, fold change; GM, gray matter; WM, white matter.

Fig. 2 |. Consensus network analysis revealed shared mechanisms overlapping with neurodegenerative disorders.

Four disease-related consensus modules—MEturquoise, MEgreenyellow, MEgreen and MEblue—exhibiting distinct cell-type-specific signatures (Extended Data Fig. 4b) were identified using WGCNA. a,c,e,g, Volcano plots illustrating DE genes in the frontal gray matter and pontine nuclei with genes from each module highlighted. Selected top GO enrichment terms of each module are listed. MEturquoise (a) is enriched in deep-layer excitatory neurons (FDR-adjusted P < 1.00 × 10⁻⁵, EWCE) and PD pathways (BH-adjusted P = 1.08 × 10⁻⁴⁵, Metascape). MEgreenyellow (c) is enriched in astrocytes (FDR-adjusted P < 1.00 × 10⁻⁵, EWCE), and metallothioneins bind metal-associated pathways (BH-adjusted P = 4.31 × 10⁻¹⁰, Metascape). MEgreen (e) is enriched in the BBB (FDR-adjusted P = 1.69 × 10⁻², EWCE), oxidative stress-induced senescence and NABA matrisome-associated pathways (BH-adjusted P = 7.59 × 10⁻⁴ and P = 9.74 × 10⁻³, Metascape). MEblue (g) is enriched in microglia/BBB (FDR-adjusted P = 2.4 × 10⁻⁴ and 6.60 × 10⁻³, EWCE), and pathways associated with lipid transport, detection of chemical stimulus, S. aureus infection, GPCR binding and response to bacterium (BH-adjusted P = 6.04 × 10⁻³, 1.84 × 10⁻², 4.92 × 10⁻², 4.95 × 10⁻² and 5.57 × 10⁻², respectively, Metascape). b,d,f,h, PPI networks of hub genes and their associated genes within each module are displayed (STRING PPI network enrichment P = 1.00 × 10⁻¹⁶ in b; 1.34 × 10⁻¹⁰ in d; 2.7 × 10⁻⁶ in f; 1.00 × 10⁻¹⁶ in h). Hub genes from each module are circled in black and highlighted in yellow.

Fig. 3 |. hdWGCNA indicated shared pathways between severe acute/subacute COVID-19 and neurodegenerative disorders.

a, Projection of consensus turquoise and greenyellow modules identified in spatial transcriptomics of frontal gray matter and pontine nuclei to the frontal cortical snRNA-seq data from Yang et al.. UMAP plots illustrating the major cell types in Yang et al.. Subclass cell types are highlighted in module projected UMAP plots. b, Dendrogram of co-expression topological overlap of genes in excitatory neurons from hdWGCNA analysis. c, Unsupervised module feature UMAP plots of 10 modules identified. d, Volcano plot illustrating DE modules between the non-flu/non-COVID controls and patients with COVID-19 based on module eigengene values. Highlighted in red is the top downregulated module, ExNeuN-M3. e, Violin and box plots of module eigengene (ME) value (38,217 cells) of ExNeuN-M3 module genes show their preferential expression in excitatory neurons and NRGN neurons. The box plots depict the interquartile range (IQR) of ME, with the bottom and top edges representing the 25th and 75th percentiles. Whiskers extend from minimum to maximum within 1.5 times the IQR from the quartiles. f, Top GO enrichment terms of ExNeuN-M3 module are related to neurodegenerative diseases, oxidative phosphorylation and synaptic vesicle pathways (Metascape). g, Hub gene network display of top hub genes in the ExNeuN-M3 modules. The red squares highlight the shared hub genes identified in the consensus turquoise module in Fig. 2b.

Fig. 4 |. Pervasive downregulation of neurodegeneration-associated mitochondrial and synaptic vesicle pathways in the brain of patients with severe acute/subacute COVID-19.

a,d, Circos plot of shared genes of differentially downregulated (a) and upregulated (d) genes (BH-adjusted P < 0.1, lme4 (ref. 33)) in the frontal gray matter (GM), white matter (WM), pontine nuclei (PN), pontocerebellar tract (PCT), corticospinal tract (CST) and tegmentum (Teg). Blue and orange lines link genes shared by multiple areas. b,e, Top 20 GO enrichment terms of downregulated (b) and upregulated (e) genes are displayed (Metascape). c, Merged PPI network and MCODE analyses of differentially downregulated genes of six areas in the frontal lobe and pons (physical network, P < 0.001, STRING). Each node is color proportionate to the hits of genes from each area and shows that most of the genes are shared by multiple areas. Top five MCODE clusters with their GO enrichment terms are displayed, including synaptic vesicle cycle, oxidative phosphorylation, prion disease, endocytosis, PD and L1CAM interactions (BH-adjusted P = 1.00 × 10⁻¹⁴, 2.00 × 10⁻⁹, 1.00 × 10⁻⁶, 3.16 × 10⁻⁵, 3.16 × 10⁻⁶ and 2.00 × 10⁻⁶, respectively, Metascape).

Fig. 5 |. Bulk tissue RNA and proteomic profiling of the frontal lobe and pons confirm the biological processes implicated by spatial transcriptomics.

a, Dot plots illustrating the cell type enrichment analysis (EWCE) of DE genes (BH-adjusted P < 0.1, DESeq2 (ref. 32)) in frontal lobe nCounter direct gene expression assay (n = 12 for non-COVID controls; n = 17 for patients with COVID-19). Subclass cell subtypes from Wamsley et al., Yang et al. and Bakken et al. are used as references. Red and blue indicate the enrichment for upregulated and downregulated genes, respectively. The color intensity corresponds to the sd_from_mean of each cell type. The grayscale of the dot edge corresponds to the FDR-adjusted P value of each cell type. Endo, endothelial cell; EXT, excitatory neuron; IT, intratelencephalic cell; MG or Micro, microglia; ODC or Oligo, oligodendrocyte; PV or Pvalb, parvalbumin; PVM, perivascular macrophage; Sst, somatostatin; Vip, vasointestinal peptide. b, Heatmaps illustrating the NES values from GSEA^, of frontal lobe and pons spatial transcriptomics and proteomics with BH-adjusted P value cutoff of 0.025 and 0.1, respectively. Shared gene sets/biological pathways between proteomics and transcriptomics are highlighted by the dotted boxes. Orange highlights mitochondria-related gene sets; blue highlights disease-related gene sets. c, Bar plots showing the number of overlapped DE proteins (P < 0.05, MSstatsTMT) in frontal lobe and pons with genes in the four consensus modules identified in the spatial transcriptomics in Fig. 2 with P values from GeneOverlap test indicated. d, Genes with top kME value in the consensus turquoise module (Fig. 2a,b) are selected. Their corresponding expression (log₂FC) in the nCounter direct gene expression and proteomics of frontal lobe and pons are displayed in the heatmap (** indicates BH-adjusted P < 0.1; * indicates P < 0.05).

Fig. 6 |. Proteomic profiling of eight brain regions of patients with severe acute/subacute COVID-19 identifies shared changes.

a, Left, diagram of eight regions profiled in the proteomics analysis (126 biological independent samples; n = 54 for non-COVID controls and n = 72 for patients with COVID-19). Right, dot plot of EWCE results of DE proteins (P < 0.05, MSstatsTMT) in eight regions. Eight major cell types from Wamsley et al. were used. Red and blue indicate the enrichment for upregulated and downregulated genes, respectively. The color intensity corresponds to the sd_from_mean of each cell type. The grayscale of the dot edge corresponds to the FDR-adjusted P value. The dot size in each row represents the relative abundance across eight regions. b,c, Shared GO enrichment term cluster network (Metascape) of the downregulated (b) and upregulated (c) proteins. The GO terms with a kappa similarity score greater than 0.3 were connected by edges to form clusters and color coded by cluster entity. Representative term(s) from each cluster were selected for labeling. The magenta cluster in b with leading GO terms of PD and TCA cycle/respiratory electron transport is further plotted proportioned to the number of gene hits in each GO term across all the regions (b inset, Metascape). Detailed functional terms are listed in Supplementary Tables 3 and 4.

Fig. 7 |. Neuronal and mitochondrial dysfunction in the substantia nigra in the midbrain of patients with severe acute/subacute COVID-19.

a,c, Midbrain ROI selection diagram in biological independent non-COVID controls (n = 4) and patients with COVID-19 (n = 4). A total of 79 ROIs (white square) were selected. Multiplex immunofluorescence highlights cerebral peduncle (CP, a), red nucleus (RN, a), periaqueduct gray matter (PeriGM, a) and substantia nigra (SN, c). Magenta, anti-NeuN; blue, nuclei; greenyellow, anti-GFAP; cyan, anti-IBA1. Scale bar in whole tissue section, 1 mm. Scale bar in ROI, 50 μm. b, Dot plot of differential expression of SARS-CoV2 viral probes and viral-related genes in four areas in the midbrain. Yellow represents upregulation; blue represents downregulation. The color intensity corresponds to the effect size (log₂FC). The grayscale of the dot edge corresponds to the BH-adjusted P value. The dot size represents the relative abundance across the assessed areas. c, Representative histological sections with well-preserved pigmentation in the DA neurons at substantia nigra pars compacta (SNpc) in both cohorts (12 controls and 17 patients examined). Scale bar, 100 μm. Multiplex immunofluorescence highlights DA neurons in the substantia nigra ROIs in both groups (25 ROIs selected). Scale bar, 50 μm. d, Volcano plot of differential expression profile (LMM, Dream) in the substantia nigra. Genes associated with neuronal markers, oxidative phosphorylation and PD are decreased and highlighted. Metallothionein (MT)-related genes are increased and highlighted. e, Volcano plots showing the differential expression of 90 oxidative phosphorylation genes in the substantia nigra and red nucleus (all BH-adjusted P < 0.05). FC, fold change.