Single-nucleus transcriptome analysis of human brain immune response in patients with severe COVID-19
Genome Medicine volume 13, Article number: 118 (2021)
Coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, has been associated with neurological and neuropsychiatric illness in many individuals. We sought to further our understanding of the relationship between brain tropism, neuro-inflammation, and host immune response in acute COVID-19 cases.
Three brain regions (dorsolateral prefrontal cortex, medulla oblongata, and choroid plexus) from 5 patients with severe COVID-19 and 4 controls were examined. The presence of the virus was assessed by western blot against viral spike protein, as well as viral transcriptome analysis covering > 99% of SARS-CoV-2 genome and all potential serotypes. Droplet-based single-nucleus RNA sequencing (snRNA-seq) was performed in the same samples to examine the impact of COVID-19 on transcription in individual cells of the brain.
Quantification of viral spike S1 protein and viral transcripts did not detect SARS-CoV-2 in the postmortem brain tissue. However, analysis of 68,557 single-nucleus transcriptomes from three distinct regions of the brain identified an increased proportion of stromal cells, monocytes, and macrophages in the choroid plexus of COVID-19 patients. Furthermore, differential gene expression, pseudo-temporal trajectory, and gene regulatory network analyses revealed transcriptional changes in the cortical microglia associated with a range of biological processes, including cellular activation, mobility, and phagocytosis.
Despite the absence of detectable SARS-CoV-2 in the brain at the time of death, the findings suggest significant and persistent neuroinflammation in patients with acute COVID-19.
Coronavirus disease 2019 (COVID-19), caused by the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is currently the most urgent health care issue in the world. The central nervous system (CNS) is not the primary organ affected by SARS-CoV-2; however, there is increasing evidence that the brain is affected by SARS-CoV-2 by multiple mechanisms, each with potential short- and long-term impacts on neurological and cognitive function in affected individuals. As a result, SARS-CoV-2 has been associated with an unusually wide range of neurological and neuropsychiatric manifestations, ranging from asymptomatic cases to prolonged states of disability following initial recovery from COVID-19-associated delirium [1,2,3,4].
Systematically studying neurological disease in COVID-19 patients presents several challenges, including only a subset of the patient population displaying neurological symptoms, an inability to sample CNS tissues directly in living individuals, and difficulties in distinguishing direct neuro-invasion from systemic viremia within the brain. Brain autopsies and neuroimaging studies have demonstrated acute hypoxic and vascular injury [5, 6], as well as plausible SARS-CoV-2 CNS tropism [7,8,9,10]. Studies employing neuronal organoids have provided conflicting results regarding the effect of SARS-CoV-2 on neurons [9, 11, 12] or on epithelial cells in the choroid plexus [13, 14]. Overall, existing studies provide multiple mechanisms through which SARS-CoV-2 affects the human brain, including direct invasion and infection of specific types of neurons and glia, in addition to endothelial injury and vascular coagulopathy and diffuse neuroinflammatory processes or systemic inflammatory and hypercoagulable states.
To fully understand the impact of COVID-19 on the CNS, it is critical to perform a direct, unbiased, and comprehensive high-resolution exploration of the transcriptomic landscape in human brain tissue from patients with COVID-19. Here, we performed droplet-based single-nucleus transcriptome profiling in the brain regions that are potentially implicated in clinical manifestations of COVID-19: the prefrontal cortex, which is involved in higher-order cognitive function; medulla oblongata, which includes the respiratory center; and choroid plexus, which provides a toxin barrier to the brain and forms an interface between the blood and cerebrospinal fluid. The analysis described here examined the relationship between SARS-CoV-2 infection and neuroinflammation, brain tropism, and host immune response.
Samples used in this study
The study cohort consisted of tissue from 5 COVID-19 patients and 4 controls. For each donor, we studied 3 brain regions (PFC, medulla, and ChP). All specimens were obtained and de-identified by the biorepository at the Icahn School of Medicine at Mount Sinai in accordance with the policies and regulations at the Icahn School of Medicine at Mount Sinai (ISMMS) and its institutional review board. All samples were subjected to western blot, targeted RNA-seq, and RNA-FISH for viral spike protein, as well as transcriptomic analysis via snRNA-seq. We note that all COVID-19 cases were classified as severe; the clinical characteristics of donors are detailed in Additional file 1: Table S1.
SARS-CoV-2 viral load quantification in frozen specimens of postmortem human brain
Brain tissue was harvested based on a modification of published protocols  for rapid dissection and systematic neuroanatomical sampling. For each donor, we targeted 3 brain regions, which included the dorsolateral prefrontal cortex (PFC), medulla oblongata (medulla), and choroid plexus (ChP). For the immunoblotting and SARS-CoV-2-targeted RNA-seq assay, we included a single dissection with the exception of the PFC, from which we included separate dissections of cortical gray matter and white matter. Western blot analysis was performed using 100 μg of total protein and anti-SARS-CoV-2 spike glycoprotein (Abcam, Cambridge, UK. Cat no: ab272504) or anti-β-actin (Abcam, Ab4970S) primary antibodies. For the SARS-CoV-2 targeted assay, we used the AmpliSeq Library Plus and cDNA Synthesis kits from Illumina (Illumina, San Diego, CA. Cat nos: 20019103 and 20022654), and, following quantification and pooling, libraries were run on a NovaSeq 6000 S4 (Illumina) in a 2 × 150 run format. For RNA-FISH, 10-μm sections were incubated with anti-sense DNA probes against the SARS-CoV-2 spike protein RNA sequence. Slides were incubated in TrueBlack Lipofuscin Autofluorescence Quencher (Biotium), and images were acquired on a Zeiss LSM780 confocal microscope. For further details, see Additional file 2: Supplementary methods.
For immunohistochemical staining, fixed 40-μm sections of the choroid plexus from cases and controls were stained with anti-CD68 antibody (Abcam, ab955), and images were obtained on a Zeiss LSM780 confocal microscope and processed with ImageJ. For further details, see Additional file 2: Supplementary methods.
Generation of single-nucleus RNA-seq from postmortem human brain
To better control for donor batch effects in the snRNA-seq analyses, we performed two dissections for each brain region and individual. The nuclei were extracted from 50 mg of frozen tissue. Each batch included 6 tissue samples, and the nuclei were incubated with individual nuclear hashing antibodies (BioLegend, San Diego, CA. TotalSeqA MAb414 anti-Nuclear pore complex protein) to facilitate multiplexed loading of the single-cell microfluidics device (10x Genomics, Pleasanton, CA). A total of 46,560 nuclei (7760 each) were loaded, in duplicate, on 10x Genomics B chips using 3′ capture chemistry. snRNA-seq and hash-tag oligo (HTO) libraries were generated separately and were sequenced by Nova-seq (Illumina) obtaining 2 × 100 paired-end reads.
Single-nucleus data processing and demultiplexing
10x Genomics paired-end sequencing reads were processed and aligned on a pre-mRNA reference genome using cell ranger v3.1.0. HTO sequencing reads were preprocessed using kallisto indexing and tag extraction . We customized an algorithm to identify singlets and doublets using HTO barcodes (see Additional file 2: Supplementary methods for more details). Demultiplexed unique molecular identifier matrices were combined and labeled by their associated experiment batch and brain region. The nuclei with fewer than 200 expressing genes or higher than 5% mitochondrial reads were omitted. We applied DoubletFinder  to remove doublets that were missed by HTO demultiplexing. We used Seurat  to perform data normalization, UMAP dimension reduction, and graph-based clustering. Harmony  was used to adjust batch effects using the first 50 principal components of the transformed count matrix. UMAP projection and Louvain clustering were then performed based on the shared nearest neighbor graph of these 50 harmonized principal components. We annotated clusters by inspecting a selected set of canonical markers.
Further clustering was performed on the pooled population of immune cells that includes microglia, monocytes, and macrophages. The clustering and UMAP visualization were based on 2000 variable genes of the pooled population. We applied Harmony  to calibrate batch effects. After manual inspection, a few clusters that were seemingly doublets were excluded from the subclustering analysis. We identified 9 clusters which are annotated as 7 subclusters of microglia (Mic-1 to Mic-7) and 2 subclusters of monocytes/macrophage (Mo-1 and Mo-2) based on the relative abundance of Mic and Mo/MP in each cluster. Marker genes of these clusters, comparing one cluster to the other 8 clusters, were calculated by the Wisconsin rank sum test as implemented in Seurat .
Statistical analysis on cell composition
To identify cell composition changes associated with COVID-19, we computed the fraction of cells annotated with a cell type in each brain region and then fitted a linear mixed model:
Associated cell types were then identified by contrasting cases and controls for each tissue. We repeated statistical tests on the normalized cell fraction using centered log ratio transformation and obtained similar results.
Identification of differentially expressed genes
Differentially expressed genes (DEGs) were identified using a linear mixed model. We first split the combined count matrix of 10,000 variable genes into sub-matrices containing the nuclei from one annotated cell population. K-nearest neighbor smoothing  was applied to each sub-matrix to impute dropped-out genes. Genes expressed in less than 1% of cells in a cell population were further filtered. A pseudo-bulk count matrix was formed by aggregating single-cell reads per sample. Using the dream pipeline , we fitted a mixed linear model to the pseudo-bulk count matrix:
DEGs were then identified by contrasting cases and controls for each tissue. We omitted the test results on genes that were sparsely expressed (< 1% cells) in each tissue.
Cell type enrichment
We performed Fisher’s exact test to compare gene signatures of annotated cell types with those reported in other studies. To generate the set of marker genes for an annotated cell type, we used the Wilcoxon rank sum test, implemented by the FindMarker function, to compare cells of one cell type against the rest. Up to 50 marker genes, with adjusted P-value less than 0.05 and ranked by fold changes, were included to form the set of marker genes. We curated gene sets of known cell types from 4 publicly available data sets [8, 22, 23], which are further described in Additional file 2: Supplementary methods. The background gene sets were set to be the intersection of 10,000 variable genes and all possible genes from the other data source.
Construction and analysis of the transcription factor-gene network
We used pySENIC [24, 25] to identify transcription factor regulons. The count matrix of 10,000 variable genes was used. Genes expressed in less than 1% of cells were further filtered as recommended by the pySCENIC protocol. The gene co-expression network was inferred using the gradient boosting machine implemented by arboreto. Enriched motifs for a gene co-expression module were predicted using pre-computed databases from cisTargetDB and the ctx function in pySCENIC. Lastly, activity scores of inferred regulons were quantified at the single-cell level using AUCell.
To identify regulons associated with COVID-19, we first computed the averaged AUCell score per sample and fitted a linear mixed model:
Associated regulons were then identified by contrasting cases and controls for each tissue. We restricted our analysis to the 5 regulons that are most specific to a cell type ranked by the regulon-specific score .
Pseudo-temporal trajectory score
To construct a pseudo-time trajectory that captures disease progression in the microglia, we computed the principal components using the gene count matrix of the top DEGs across three brain regions. We constructed a pseudo-time trajectory along the first two principal components using Slingshot . The order of the pseudo time trajectory is chosen so that the fraction of cells from COVID patient is increasing with the pseudo-time. Gene expression profiles over the pseudo-time trajectory were analyzed using tradeSeq . Applying resampling-based sequential ensemble clustering to genes with non-monotonic profiles, we identified 22 clusters, but a strong similarity between clusters was observed. Using a step-wise approach, we merged meta-clusters into 4 categories (see also Additional file 2: Supplementary methods).
Transcriptome-wide association study and gene set enrichment
We used the “B2_ALL_eur_leave_23andme” summary statistics from the Release 4 (October 2020) of COVID-19 Host Genetics Initiative, which corresponds to the phenotype “Hospitalized covid vs. population, leave out 23andMe.” We employed blood and brain EpiXcan  tissue imputation models, trained in the STARNET [29, 30] and PsychENCODE [31, 32] cohorts, respectively. To derive gene-tissue-trait associations, transcriptomic imputation models were applied to the summary statistics following the S-PrediXcan approach . For each tissue, we kept genes with r2 of the correlation between cross-validated prediction and observed expression greater than, or equal to, 0.01. For the gene set enrichment analysis, only protein-coding genes were considered, and the enrichment was tested with Fisher’s exact test.
SARS-CoV-2 viral load quantification across multiple brain regions
To characterize the CNS effect of SARS-CoV-2, we performed viral load quantification in human brain tissue from 5 COVID-19 patients and 4 controls (Fig. 1A). We note that all COVID-19 cases were classified as severe; the clinical characteristics of donors are detailed in Additional file 1: Table S1. For each donor, we targeted 3 brain regions (PFC, medulla, and ChP). Immunoblotting was negative for the presence of viral spike S1 protein in all tissues examined (Additional file 2: Fig S1). We then performed transcriptome analysis covering > 99% of SARS-CoV-2 genome and all potential serotypes. For each brain region and donor, we included a single dissection with the exception of the PFC, from which we included separate dissections of cortical gray matter and white matter (for an illustrative example, see Additional file 2: Fig S2), and generated, on average, 1.2 million reads per library. Across all samples, none of the sequencing reads mapped to the SARS-CoV-2 genome (Additional file 3: Table S2). In addition, examination of additional brain regions (red nucleus and substantia nigra), using fluorescence in situ hybridization (FISH) for SARS-COV-2 spike protein, also failed to detect the virus (Additional file 2: Fig S3). Overall, by employing three different experimental approaches, and exploring multiple brain regions, we did not detect SARS-CoV-2 in the postmortem brain tissue.
snRNA-seq profiling across multiple brain regions in COVID-19 patients and controls
We characterized the molecular and cellular perturbations in the CNS of COVID-19 patients, independent of SARS-CoV-2 direct invasion, by performing droplet-based single-nucleus RNA sequencing (snRNA-seq) in the PFC, medulla, and ChP in the same set of 5 patients and 4 controls (Fig. 1A). To better control for donor batch effects, we performed two dissections for each brain region and individual (Additional file 2: Fig S2), and all samples per donor (n = 6) were pooled together using nuclear hashing. After preprocessing of snRNA-seq data, demultiplexing using hashtag-oligo intensity and quality control, 68,557 high-quality single-nucleus barcodes, with a median of 2817 detected genes per cell, were available for downstream analysis (Additional file 2: Fig S4). Variance in the gene expression was mostly driven by biological factors (cell type, brain regions, and donor) (Additional file 2: Fig S5).
We performed de novo taxonomy based on the graph-based clustering and uniform manifold approximation and projection (UMAP) across all brain regions and samples, and identified 15 major cell clusters (Fig. 1B). Clustering was independent of donor effect and technical variables, while differences between the brain regions were preserved (Additional file 2: Fig S6). Annotation of cell clusters based on the expression of canonical gene markers identified the following populations: excitatory neurons (Ex) that express SYT1 and SLC17A7, inhibitory neurons (In) that express SYT1 and GAD1, astrocytes (Ast1 and Ast2) that express AQP4, ependymal cells (Ep) that express CFAP299, oligodendrocyte progenitor cell (OPC) that express VCAN, oligodendrocytes (Oli) that express MOBP, epithelial cells (Epi) that express HTR2C, endothelial cells (End) that express FLT1, mesenchymal cells (Mes) that express COL1A1, pericytes (Per1 and Per2) that express PDGFRB, microglial cells (Mic) that express APBB1IP, and lymphocyte (LM) that express CD96 and monocytes/macrophage (Mo/MP), expressing CD163 (Fig. 1 and Additional file 2: Fig S7). Gene set enrichment analysis showed overlap for expected molecular pathways and functions, such as myelination for oligodendrocytes, chemical synaptic transmission for excitatory, and inhibitory neurons, and T cell activation for lymphocytes (Additional file 2: Fig S8). The expression profiles of each cell type show high concordance with previous snRNA-seq in human brain tissue [8, 22], peripheral blood cells, and brain organoids  (Additional file 2: Fig S9), indicating the robust definition of cell subpopulations in the current study.
Compositional analysis identifies changes in the proportion of immune-related populations in COVID-19 patients
We assessed the relative proportions of the 15 major cell types in COVID-19 cases compared to controls across the three brain regions. For each cell cluster, we applied a linear mixed model to detect the interaction between COVID-19 cases and brain regions, while controlling for donor effects. Among 45 combinations of cell types and brain regions, we identified 2 cell populations from choroid plexus, including Mes cells and Mo/MP, showing a significant increase in their relative proportions in COVID-19 cases (Fig. 1D). We did not detect any significant COVID-19-associated changes in the cell type composition in either the PFC or medulla (Additional file 2: Fig S10). We next sought to confirm the presence of increased numbers of Mo/MP in the ChP by staining with an antibody against CD68, which is a marker for Mo/MP cells. We observe an abundance of CD68-positive cells in ChP of COVID-19 cases, and these cells do not appear to be associated with a blood vessel (Additional file 2: Fig S11). Overall, these results suggest that, in COVID-19, Mo/MP extravasate from the blood vessels into the stroma of the ChP, which is composed of Mes cells.
To further explore the cell states of immune-related cells associated with COVID-19 disease status, we then performed subclustering on the Mic and Mo/MP populations. We identified 9 clusters that corresponded to 7 subclusters of Mic and 2 subclusters of Mo/MP (Additional file 2: Fig S12a). The relative abundance of Mic-1 and Mic-2 subclusters in PFC showed significant changes in COVID-19 patients in the opposite direction compared to controls (Additional file 2: Fig S12b), suggesting a transition from Mic-1 to Mic-2 in PFC with COVID-19. Interestingly, expression of marker genes show that Mic-2 expresses, at a higher level, the Complement C3 gene (Additional file 2: Fig S12c) which has previously been shown to drive microglia activation . The relative abundance of Mo/MP-1 in ChP increases significantly while Mo/MP-2 remains stable. Compared to Mo/MP-2, the Mo/MP-1 subcluster expresses LYVE1, which is a marker that is highly expressed in a subset of macrophages in the meninges .
Differential gene expression analysis uncovered activation of innate immune cells in COVID-19 brain parenchyma
Molecular processes and biological pathways affected by disease states can be detected by studying the differences in gene expression among cases and controls within each cell subpopulation. For each cell type and brain region, we applied linear mixed models to identify differentially expressed genes (DEGs) among COVID-19 patients and controls, while controlling for donor effects. Among the 15 cell types and 3 brain regions, the microglia in the PFC showed the highest number of perturbations, including 178 DEGs (Fig. 2A), followed by monocytes in PFC and ChP. Interestingly, the majority of DEGs were upregulated in the immune cells (Mic and Mo/MP), indicating increased immune activity in COVID-19 patients. We provide a summary of differentially expressed genes across all cell types in Additional file 4: Table S3.
We examined whether gene perturbations affected specific biological processes. To preserve sufficient statistical power, we focused on the DEG signatures detected in microglia from the PFC. Gene set enrichment analysis identified biological processes such as “macrophage activation” (8 genes, P = 9.3 × 10−6) and “phagocytosis” (14 genes, P = 1.9 × 10−6) as being enriched with the 178 DEGs in PFC microglia (Additional file 5: Table S4). To further investigate transcriptomic changes in canonical pathways, we calculated the activity scores across 186 KEGG molecular pathways (Additional file 2: Fig S13). We then applied linear mixed models and identified differences in average activity levels across COVID-19 patients and controls. We identified 16 pathways showing significant differences in activity levels in the microglia of the PFC (Additional file 6: Table S5). As an illustrative example, we show the expression levels of the four most significant upregulated pathways (Fig. 2B) and associated genes (Fig. 2C). Interestingly, the majority of significant pathways were immune-related, and all of them, except one, which was related to “steroid biosynthesis,” were upregulated in COVID-19 patients. Depletion of steroid synthesis in COVID-19 patients is consistent with the clinical evidence that steroid administration is beneficial to modulate inflammatory response . Overall, these results suggest the strong activation of innate immune cells in COVID-19 brain parenchyma.
Pseudo-temporal analysis detected transition of inflammatory response in PFC across a range of biological processes
Statistical modeling can be applied to snRNA-seq to extract temporal information and study dynamic biological processes from cross-sectional data sets. We calculated a pseudo-temporal trajectory score (PTS) in the microglia, based on the progression of the transcriptional dysregulation in COVID-19 patients compared to controls (Fig. 3A). PTS has a wide distribution, potentially indicating microglia in different stages of activation and, as expected, COVID-19 patients demonstrated higher PTS than controls (Fig. 3B). To rule out patient-specific batch effects, we performed a similar analysis stratified by donor and did not observe any effect beyond case/control status (Additional file 2: Fig S14).
We categorized 646 commonly expressed microglia genes into four groups (increasing, decreasing, early transient, and late transient) based on the expression patterns capturing progressive changes related to PTS (Fig. 3C, Additional file 7: Table S6). The majority of genes were clustered as “increasing” (579 genes), followed by late transient (36 genes), early transient (16 genes), and decreasing (15 genes). Genes within the “increasing” cluster were more perturbed in COVID-19 patients (estimated based on pi1 = 0.683), compared to the other 3 clusters (range of pi1 = 0.051 to 0.077) and were enriched for 452 biological pathways, including “regulation of immune system process” (136 genes, P = 1.75 × 10−16) and “apoptotic process” (104 genes, P = 3.19 × 10−5). In the “early transient” group, 8 out of the 16 genes, including CD83 and 3 heat shock proteins (HSP90AA1, HSPB1, HSPH1), belong to the “cell activation” pathway (P = 6.40 × 10−6), while the “late transient” group was enriched for the “cell mobility” pathway (P = 4.44 × 10−4). Overall, the pseudo-temporal analysis supports a model where myeloid-driven inflammatory response in PFC involves transition across a range of biological processes, including cellular activation, mobility, and phagocytosis.
Gene regulatory network (GRN) analysis identified activated microglia response in patients with COVID-19
Gene regulatory networks (GRNs) define the co-expression patterns of transcripts, considering the regulatory relationships between transcription factors (TFs) and their target genes. This analysis can determine cellular functions and model different systemic behaviors to uncover gene-level relationships in cells associated with disease states. To further understand the impact of SARS-CoV-2 on the TF-gene relationships, we explored the differences in GRNs among COVID-19 cases and controls. Across all cell subpopulations, we identified 131 TF modules that regulated, on average, 272 genes per module  (Additional file 8: Table S7). UMAP projection based on activity scores of GRNs reaffirmed the robustness of annotated cell types (Additional file 2: Fig S15).
We defined cell population specificity by ranking TF modules according to regulon score  and uncovered well-known cell type-specific TFs, such as PAX6 for astrocytes and IRX8 for microglial cells (Additional file 2: Fig S16). We tested whether changes in the activity level of the top 5 TFs for each cell population were associated with COVID-19. Among the cell types and brain regions, TF modules in PFC microglia were the most affected (Additional file 2: Fig S17). We prioritized 4 out of 5 microglia PFC TFs (IRF8, ATF5, SPI1, TAL1) based on the upregulation in their activity in patients with COVID-19 (Z > 2.5, P < 0.05 and FDR within cell type < 0.05) (Fig. 4A; Additional file 9: Table S8). Projecting microglia DEGs onto the GRNs of these 4 TFs showed the co-regulatory TF-gene patterns affected in COVID-19 (Fig. 4B). Collectively, these results suggest GRNs corresponding to activated microglia response in patients with COVID-19 and nominate TFs that partially regulate those transcriptome changes.
GRNs corresponding to microglia activation in acute COVID-19 patients demonstrate enrichment for genes that are predicted to be downregulated in a transcriptome-wide association study (TWAS) of severe COVID-19
To assess if acute microglia activation has a beneficial or deleterious effect, we examined whether the direction of effect of transcriptome signatures is concordant with genetically regulated changes in the gene expression that are associated with severe COVID-19. We leveraged genetic variation from the COVID-19 Host Genetics Initiative and the gene expression models from the brain [31, 32] and blood  tissues to impute the genetically regulated transcriptomic changes  associated with severe COVID19 outcomes (Additional file 10: Table S9). Since microglia gene expression models are not currently available, we included proxy tissue (homogenate brain and whole blood) models to capture the transcriptome profiling of the microglial cell lineage. Given that the microglia only represent 0.5 to 16.6% of all cells in the brain parenchyma , we also included blood tissue models corresponding to immune cells which are close to the microglial cell lineage [38, 39]. We identified 12 significant genes (at FDR 5%) in the brain and blood (AP000295.1, CCR3, CR936218.2, CRHR1, FYCO1, IFNAR2, IL10RB, IL10RB-AS1, LRRC37A4P, MAPT-AS1, OAS1, OAS3) that were associated with hospitalized COVID-19 patients, with respect to the general population. Nominally significant gene-trait associations (P < 0.05) from the imputed blood transcriptome were enriched (OR = 1.64, P = 0.030, Fisher’s exact test) with GRNs that are associated with the top 4 microglia TFs (IRF8, ATF5, SPI1, TAL1) (Fig. 4C). In a secondary analysis, we examined whether the genetic liability was associated with a predicted downregulation or upregulation of the genes participating in these TF modules. Enrichment was observed only for genes predicted to be downregulated (OR = 1.89, P = 0.035, Fisher’s exact test), compared to genes predicted to be upregulated in susceptible individuals (OR = 1.35, P = 0.25, Fisher’s exact test). In addition, OAS1, which was predicted to be downregulated in susceptible individuals (FDR-adjusted P = 0.003), was involved in 3 microglia GRNs (SPI1, IRF5, TAL1, Fig. 4B). Overall, one possible interpretation of these observations is that these TF modules, which are activated in the microglia of acute COVID-19 patients and are predicted to be hypo-active at baseline in susceptible individuals, may have a beneficial role in the host response.
Recent advances in single-cell approaches provide an opportunity to study highly complex tissues like the human brain at unprecedented resolution. In order to better understand the impact of acute COVID-19 on the CNS, we studied the effects of COVID-19 in individual cells across 3 functionally distinct regions of the human brain (prefrontal cortex, choroid plexus, and medulla oblongata). Although no virus was detected, single-nucleus gene expression analysis revealed extensive differences in the brains of COVID-19 patients when compared to controls; specifically in the ChP and PFC. We observed a relative increase in the proportions of infiltrating immune cells in the ChP, suggesting potential migration of monocytes/macrophages across the blood-brain barrier. Microglia residing in the PFC of COVID-19 patients displayed dysregulated gene expression. The majority of the microglial DEGs were upregulated, mediating a myeloid-driven inflammatory response that involved a range of biological processes, including cellular activation, mobility, and phagocytosis. This is consistent with previous studies that have also described increased inflammatory response of microglia in COVID-19 cases [8, 10]. Finally, by leveraging genetic variation to infer differences in COVID-19 susceptible individuals, we provided preliminary evidence for a potential beneficial role for microglia activation during the acute COVID-19 phase.
Although there is evidence that SARS-CoV-2 spike protein can be detected in the brain, including the cortex , choroid plexus , and medulla oblongata [10, 40], immunoblotting, immunohistochemistry, and viral genome RNA-seq indicate that the virus was not present at the time of death in the specimens included in this study. The ability to detect SARS-CoV-2 in the CNS is affected by the duration of COVID-19 infection, with a marked decrease in detectable virus by day 20 . In our study, all individuals had been infected for ≥ 14 days while the duration of infection was ≥ 30 days for 3 out of 5 COVID-19 cases. In addition, only a subset of COVID-19 patients indicates non-zero SARS-CoV-2 RNA copies in the CNS, which are more difficult to detect in the brain parenchyma compared to the olfactory mucosa . Although our focus was on immune cells, there is evidence, in addition to microglia activation, for COVID-19-related transcriptional changes in a range of brain cell types including astrocytes, oligodendrocytes, and excitatory neurons . The observed differences in the number of DEGs, and the cell types affected, might be explained by the experimental design: two versus a single dissection per brain region and individual. If we only consider a single dissection per brain region and individual in our analysis, the number of DEGs increases (data not shown) and includes perturbations among every major CNS cell type.
Taken together, these findings indicate persistent activation of the innate immune response in the brains of patients with COVID-19. Based on our results, it is possible that the inflammatory response of microglia is induced as a result of peripheral immune cells infiltrating CNS through the blood-brain barrier. Another point of entry of SARS-CoV-2 to the CNS is by crossing the neural-mucosal interface in olfactory mucosa . These two mechanisms are non-mutually exclusive and might be associated with different stages of disease progression and presentation of clinical symptoms. While shedding light on the impact of SARS-CoV-2 on the CNS, our study was insufficiently powered to fully elucidate all of the relevant cellular states associated with COVID-19, and future research efforts are required to confirm, and expand on, our findings. In conclusion, our study suggests extensive neuro-inflammation and brain immune response in acute COVID-19 patients, even in the absence of direct evidence of SARS-CoV-2 neuro-invasion.
Availability of data and materials
Processed and raw data can be downloaded from NCBI GEO (GSE164485) : https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE164485.
Coronavirus disease 2019
Severe acute respiratory syndrome coronavirus 2
Single-nucleus RNA sequencing
Fluorescence in situ hybridization
Central nervous system
Differentially expressed genes
Oligodendrocyte progenitor cell
Uniform manifold approximation and projection
pseudo-temporal trajectory score
- HSP90AA1, HSPB1, HSPH1:
Heat shock proteins
Gene regulatory networks
Genome-wide association studies
Mao L, Jin H, Wang M, Hu Y, Chen S, He Q, et al. Neurologic manifestations of hospitalized patients with coronavirus disease 2019 in Wuhan, China. JAMA Neurol. 2020;77(6):683–90. https://doi.org/10.1001/jamaneurol.2020.1127.
Frontera JA, Sabadia S, Lalchan R, Fang T, Flusty B, Millar-Vernetti P, et al. A prospective study of neurologic disorders in hospitalized COVID-19 patients in New York City. Neurology. 2021;96(4):e575–e586. https://doi.org/10.1212/wnl.0000000000010979.
McLoughlin BC, Miles A, Webb TE, Knopp P, Eyres C, Fabbri A, et al. Functional and cognitive outcomes after COVID-19 delirium. Eur Geriatr Med. 2020;11(5):857–62. https://doi.org/10.1007/s41999-020-00353-8.
Helms J, Kremer S, Merdji H, Clere-Jehl R, Schenck M, Kummerlen C, et al. Neurologic features in severe SARS-CoV-2 infection. N Engl J Med. 2020;382(23):2268–70. https://doi.org/10.1056/NEJMc2008597.
Solomon IH, Normandin E, Bhattacharyya S, Mukerji SS, Keller K, Ali AS, et al. Neuropathological features of COVID-19. N Engl J Med. 2020;383(10):989–92. https://doi.org/10.1056/NEJMc2019373.
Lee MH, Perl DP, Nair G, Li W, Maric D, Murray H, et al. Microvascular injury in the brains of patients with COVID-19. N Engl J Med. 2021;384(5):481–3. https://doi.org/10.1056/NEJMc2033369.
Puelles VG, Lütgehetmann M, Lindenmeyer MT, Sperhake JP, Wong MN, Allweiss L, et al. Multiorgan and renal tropism of SARS-CoV-2. N Engl J Med. 2020;383(6):590–2. https://doi.org/10.1056/NEJMc2011400.
Yang AC, Kern F, Losada PM, Agam MR, Maat CA, Schmartz GP, et al. Dysregulation of brain and choroid plexus cell types in severe COVID-19. Nature. 2021. https://doi.org/10.1038/s41586-021-03710-0.
Song E, Zhang C, Israelow B, Lu-Culligan A, Prado AV, Skriabine S, et al. Neuroinvasion of SARS-CoV-2 in human and mouse brain. J Exp Med. 2021;218(3):e20202135.
Meinhardt J, Radke J, Dittmayer C, Franz J, Thomas C, Mothes R, et al. Olfactory transmucosal SARS-CoV-2 invasion as a port of central nervous system entry in individuals with COVID-19. Nat Neurosci. 2021;24(2):168–75. https://doi.org/10.1038/s41593-020-00758-5.
Ramani A, Müller L, Ostermann PN, Gabriel E, Abida-Islam P, Müller-Schiffmann A, et al. SARS-CoV-2 targets neurons of 3D human brain organoids. EMBO J. 2020;39(20):e106230. https://doi.org/10.15252/embj.2020106230.
Dobrindt K, Hoagland DA, Seah C, Kassim B, O’Shea CP, Iskhakova M, et al. Common genetic variation in humans impacts in vitro susceptibility to SARS-CoV-2 infection. Stem Cell Reports. 2021;16(3):505–18. https://doi.org/10.1016/j.stemcr.2021.02.010.
Pellegrini L, Albecka A, Mallery DL, Kellner MJ, Paul D, Carter AP, et al. SARS-CoV-2 infects the brain choroid plexus and disrupts the blood-CSF barrier in human brain organoids. Cell Stem Cell. 2020;27(6):951–961.e5. https://doi.org/10.1016/j.stem.2020.10.001.
Jacob F, Pather SR, Huang W-K, Zhang F, Wong SZH, Zhou H, et al. Human pluripotent stem cell-derived neural cells and brain organoids reveal SARS-CoV-2 neurotropism predominates in choroid plexus epithelium. Cell Stem Cell. 2020;27(6):937–950.e9. https://doi.org/10.1016/j.stem.2020.09.016.
Vonsattel JP, Amaya Mdel P, Cortes EP, Mancevska K, Keller CE. Twenty-first century brain banking: practical prerequisites and lessons from the past: the experience of New York Brain Bank, Taub Institute, Columbia University. Cell and tissue banking. 2008;9(3):247–58. https://doi.org/10.1007/s10561-008-9079-y.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7. https://doi.org/10.1038/nbt.3519.
McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8(4):329–37.e4.
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell. 2019;177(7):1888–902.e21.
Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–96. https://doi.org/10.1038/s41592-019-0619-0.
Wagner F, Yan Y, Yanai I. K-nearest neighbor smoothing for high-throughput single-cell RNA-Seq data. bioRxiv. 2018:217737.
Hoffman GE, Roussos P. Dream: powerful differential expression analysis for repeated measures designs. Bioinformatics. 2021;37(2):192–201. https://doi.org/10.1093/bioinformatics/btaa687.
Lake BB, Chen S, Sos BC, Fan J, Kaeser GE, Yung YC, et al. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nat Biotech. 2018;36(1):70–80. https://doi.org/10.1038/nbt.4038.
Pellegrini L, Bonfio C, Chadwick J, Begum F, Skehel M, Lancaster MA. Human CNS barrier-forming organoids with cerebrospinal fluid production. Science. 2020;369(6500):eaaz5626.
Van de Sande B, Flerin C, Davie K, De Waegeneer M, Hulselmans G, Aibar S, et al. A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nat Protoc. 2020;15(7):2247–76. https://doi.org/10.1038/s41596-020-0336-2.
Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14(11):1083–6. https://doi.org/10.1038/nmeth.4463.
Suo S, Zhu Q, Saadatpour A, Fei L, Guo G, Yuan G-C. Revealing the critical regulators of cell identity in the mouse cell atlas. Cell Rep. 2018;25(6):1436–45.e3.
Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics. 2018;19(1):477. https://doi.org/10.1186/s12864-018-4772-0.
Van den Berge K, Roux de Bézieux H, Street K, Saelens W, Cannoodt R, Saeys Y, et al. Trajectory-based differential expression analysis for single-cell sequencing data. Nat Commun. 2020;11(1):1201.
Zhang W, Voloudakis G, Rajagopal VM, Readhead B, Dudley JT, Schadt EE, et al. Integrative transcriptome imputation reveals tissue-specific and shared biological mechanisms mediating susceptibility to complex traits. Nat Commun. 2019;10(1):3834. https://doi.org/10.1038/s41467-019-11874-7.
Franzén O, Ermel R, Cohain A, Akers NK, Di Narzo A, Talukdar HA, et al. Cardiometabolic risk loci share downstream cis- and trans-gene regulation across tissues and diseases. Science. 2016;353(6301):827–30. https://doi.org/10.1126/science.aad6970.
Gandal MJ, Zhang P, Hadjimichael E, Walker RL, Chen C, Liu S, et al. Transcriptome-wide isoform-level dysregulation in ASD, schizophrenia, and bipolar disorder. Science. 2018;362(6420):eaat8127.
Wang D, Liu S, Warrell J, Won H, Shi X, Navarro FCP, et al. Comprehensive functional genomic resource and integrative model for the human brain. Science. 2018;362(6420):eaat8464.
Barbeira AN, Dickinson SP, Bonazzola R, Zheng J, Wheeler HE, Torres JM, et al. Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nat Commun. 2018;9(1):1825. https://doi.org/10.1038/s41467-018-03621-1.
Schafer DP, Lehrman EK, Kautzman AG, Koyama RK, Mardinly AR, Yamasaki R, et al. Microglia sculpt postnatal neural circuits in an activity and complement-dependent manner. Neuron. 2012;74(4):691–705. https://doi.org/10.1016/j.neuron.2012.03.026.
Brezovakova V, Jadhav S. Identification of Lyve-1 positive macrophages as resident cells in meninges of rats. J Comp Neurol. 2020;528(12):2021–32. https://doi.org/10.1002/cne.24870.
RECOVERY Collaborative Group; Horby P, Lim WS, Emberson JR, Mafham M, Bell JL, et al. Dexamethasone in hospitalized patients with covid-19 - randomized controlled trial. N Engl J Med. 2021;384(8):693–704. https://doi.org/10.1056/nejmoa2021436.
Mittelbronn M, Dietz K, Schluesener HJ, Meyermann R. Local distribution of microglia in the normal adult human central nervous system differs by up to one order of magnitude. Acta Neuropathol. 2001;101(3):249–55. https://doi.org/10.1007/s004010000284.
Ginhoux F, Greter M, Leboeuf M, Nandi S, See P, Gokhan S, et al. Fate mapping analysis reveals that adult microglia derive from primitive macrophages. Science. 2010;330(6005):841–5. https://doi.org/10.1126/science.1194637.
Hoeffel G, Chen J, Lavin Y, Low D, Almeida FF, See P, et al. C-Myb(+) erythro-myeloid progenitor-derived fetal monocytes give rise to adult tissue-resident macrophages. Immunity. 2015;42(4):665–78. https://doi.org/10.1016/j.immuni.2015.03.011.
Matschke J, Lütgehetmann M, Hagel C, Sperhake JP, Schröder AS, Edler C, et al. Neuropathology of patients with COVID-19 in Germany: a post-mortem case series. Lancet Neurol. 2020;19(11):919-29.
Fullard JF, Lee H-C, Voloudakis G, Suo S, Javidfar B, Zhiping Shao Z, Peter C, Zhang W, Jiang S, Corvelo A, Wargnier H, Woodoff-Leith E, Purohit DP, Ahuja S, Tsankova NM, Jette N, Hoffman GE, Akbarian S, Fowkes M, Crary JF, Yuan G-C, Roussos P: The landscape of human brain immune response in patients with severe COVID-19. Gene Expression Omnibus. 2021. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE164485 .
Fullard JF, Lee H-C, Voloudakis G, Suo S, Javidfar B, Zhiping Shao Z, Peter C, Zhang W, Jiang S, Corvelo A, Wargnier H, Woodoff-Leith E, Purohit DP, Ahuja S, Tsankova NM, Jette N, Hoffman GE, Akbarian S, Fowkes M, Crary JF, Yuan G-C, Roussos P: The landscape of human brain immune response in patients with severe COVID-19. GitHub. 2021. https://github.com/howchihlee/covid_brain_sc .
We thank the patients and families who donated materials for these studies. We thank the members of the Roussos Laboratory for the thoughtful advice and critique and the computational resources and expertise provided by the Scientific Computing at the Icahn School of Medicine at Mount Sinai. This paper is dedicated to the memory of Mary Fowkes.
Supported by the National Institute on Aging, NIH grants R01-AG067025 (to P.R.) and R01-AG065582 (to P.R.) and Mount Sinai COVID-19 seed fund 0285VV12 (to S.A.).
Ethics approval and consent to participate
All autopsies were performed with written consent from the legal next-of-kin. All tissue samples were obtained de-identified under approved Institutional Review Board (IRB) protocols under consultation with the Ethics Committee of the Medical Board at the Mount Sinai Hospital. The study was performed under IRB-approved guidance and regulations to keep all patient information strictly de-identified. All research conformed to the principles of the Helsinki Declaration.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Clinical characteristics of donors.
Supplementary methods. Supplementary figures: Fig S1. Western Blot Analysis Investigating the presence of COVID-19 in Human Brain. Fig S2. Brain specimen dissections. Fig S3. SARS-COV2 RNA not detected in the postmortem midbrain of COVID-19 patients. Fig S4. Quality control of single nucleus data. Fig S5. Variance partition plot. Fig S6. Cell type composition across 15 major cell clusters. Fig S7. UMAP visualization of the distribution of canonical gene markers on annotated cell populations. Fig S8. Gene set enrichment analysis for marker genes of annotated populations. Fig S9. Comparison of marker genes of annotated populations from the current study and other published studies. Fig S10. Cell type composition among cases and controls. Fig S11. Immunohistochemistry of CD68 positive cells in the Choroid Plexus. Fig S12. Sub-clustering of microglial and monocyte/macrophage populations. Fig S13. KEGG pathway enrichment network. Fig S14. Pseudo-temporal trajectory score (PTS) analysis in microglia identifies gene expression signatures across different donors. Fig S15. UMAP projection based on the activity scores of 131 TF regulons. Fig S16. Regulon specificity scores (RSS) of each annotated cell population. Fig S17. Transcription factor activity scores across different cell types and brain regions.
Sequencing statistics for SARS-CoV-2 targeted RNA-seq.
Differentially expressed genes.
Gene Ontology term enrichment in differentially expressed genes of PFC microglia.
Differences in gene set activity scores among COVID-19 patients and controls.
Genes clustered by pseudo time analysis.
Transcription regulatory modules estimated by SCENIC.
Differences in TF regulon activity scores comparing COVID-19 patients and controls.
Imputed gene associations from COVID-19 GWAS.
About this article
Cite this article
Fullard, J.F., Lee, HC., Voloudakis, G. et al. Single-nucleus transcriptome analysis of human brain immune response in patients with severe COVID-19. Genome Med 13, 118 (2021). https://doi.org/10.1186/s13073-021-00933-8
- Gene expression
- Prefrontal cortex
- Choroid plexus