A map of transcriptional heterogeneity and regulatory
variation in human microglia
Adam MH Young
1-3*
, Natsuhiko Kumasaka
2*
, Fiona Calvert
2
, Timothy R. Hammond
4,5
, Andrew
Knights
2
, Nikolaos Panousis
2
, Jeremy Schwartzentruber
6
, Jimmy Liu
7
, Kousik Kundu
2
, Michael
Segel
1
, Natalia Murphy
1
, Christopher E McMurran
1
, Harry Bulstrode
3
, Jason Correia
3
, Karol P
Budohoski
3
, Alexis Joannides
3
, Mathew R Guilfoyle
3
, Rikin Trivedi
3
, Ramez Kirollos
3
, Robert
Morris
3
, Matthew R Garnett
3
, Helen Fernandes
3
, Ivan Timofeev
3
, Ibrahim Jalloh
3
, Katherine
Holland
3
, Richard Mannion
3
, Richard Mair
3
, Colin Watts
3,8
, Stephen J Price
3
, Peter J
Kirkpatrick
3
, Thomas Santarius
3
, Nicole Soranzo
2
, Beth Stevens
4,5
, Peter J Hutchinson
3
, Robin
JM Franklin
1$
& Daniel J Gaffney
2$
.
1. Wellcome Trust MRC Stem Cell Institute, University of Cambridge, Cambridge, UK, CB2
0QQ. 2. Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridgeshire,
UK, CB10 1SA. 3. Division of Neurosurgery, Department of Clinical Neurosciences,
Cambridge University Hospitals, Cambridge, UK, CB2 0QQ. 4. FM Kirby Neurobiology Center,
Boston Children's Hospital, Harvard University, Boston, USA. 5. Howard Hughes Medical
Institute, Broad Institute of Harvard and MIT, Boston, USA. 6. EMBL-EBI, Wellcome Genome
Campus, Hinxton, Cambridgeshire, CB10 1SD. 7. Biogen, Cambridge, MA, 02142, USA. 8.
Institute of Cancer and Genomic Sciences, College of Medical and Dental Sciences,
Birmingham UK, B15 2TT $ Corresponding author
Abstract
Microglia, the tissue resident macrophages of the CNS, are implicated in a broad range of
neurological pathologies, from acute brain injury to dementia. Here, we profiled gene
expression variation in primary human microglia isolated from 141 patients undergoing
neurosurgery. Using single cell and bulk RNA sequencing, we defined distinct cellular
populations of acutely in vivo-activated microglia, and characterised a dramatic switch in
microglial population composition in patients suffering from acute brain injury. We mapped
expression quantitative trait loci (eQTLs) in human microglia and show that many disease-
associated eQTLs in microglia replicate well in a human induced pluripotent stem cell (hIPSC)
derived macrophage model system. Using ATAC-seq from 95 individuals in this hIPSC model
we fine-map candidate causal variants at risk loci for Alzheimer’s disease, the most prevalent
neurodegenerative condition in acute brain injury patients. Our study provides the first
population-scale transcriptional map of a critically important cell for neurodegenerative
disorders.
preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for thisthis version posted December 20, 2019. ; https://doi.org/10.1101/2019.12.20.874099doi: bioRxiv preprint
Introduction
Microglia are tissue resident macrophages of the central nervous system and play critical roles
in neurological immune defence, development and homeostasis (Schafer and Stevens 2015;
Q. Li and Barres 2018; Salter and Stevens 2017). These highly dynamic cells are challenging
to study in the laboratory and are strongly influenced by different experimental environments
(Gosselin et al. 2017). Genetic studies also strongly implicate microglial dysfunction in
neurodegeneration (Guerreiro et al. 2013; Jonsson et al. 2013; Tansey, Cameron, and Hill
2018; Gjoneska et al. 2015) particularly in the context of the injured brain (Johnson and
Stewart 2015). Single cell transcriptomics has suggested that microglial function may vary
across age, sex and brain region (Olah et al. 2018; Keren-Shaul et al. 2017; Hammond et al.
2019; Masuda et al. 2019; Mrdjen et al. 2018; Mathys et al. 2017). Previous studies have used
frozen post-mortem tissue from existing brain banks or fresh surgical samples typically from
restricted patient groups, typically temporal lobe resections for epilepsy or peritumoral
sampling. However, variability in the post-mortem index produces substantial variation in
cellular expression (Welch et al. 2019). Because of this, studies of microglial activation in
humans have relied on ex vivo stimulation with no available data from acutely injured human
brains. The challenge of sampling also means that large scale genetic studies of microglia
have not been attempted to date. Population studies have demonstrated that individuals
subject to mild brain trauma are 5-fold more likely to develop Alzheimer’s Disease (Mackay et
al. 2019). Consequently, it is of particular importance to understand the activation of human
microglia in the context of acute brain injury together with the underpinning genetic contribution
to neurodegeneration.
Characterisation microglial cell populations
Here, we describe the analysis of human microglia isolated from 141 patients undergoing a
range of neurosurgical procedures (Figure 1a). We recruited patients from a range of
pathologies, including 51 individuals with acute brain injury (haemorrhage and trauma), who
sustained substantial parenchymal injury, enabling us to observe in vivo microglial activation.
For each individual, we isolated CD11b-positive cells and performed both single cell
(SmartSeq2) (Picelli et al. 2014) and bulk RNA-seq on each individual. After QC, we retained
112 bulk RNA-seq samples, and 9,538 single cells from 129 patients (Figure 1b). All but three
of our bulk RNA-seq samples formed a single cluster with microglia from two previous studies
(Y. Zhang et al. 2016; Gosselin et al. 2017), and were distinct from both GTEx brain and
BLUEPRINT monocytes (Figure 1c). We then compared our single cell data to public datasets
of 68K PBMCs isolated from a healthy donor (Zheng et al. 2017) and 15K brain cells from 5
GTEx donors (Habib et al. 2017). A total of 8,662 cells formed a cluster with the microglia
population found in GTEx samples and distinct from PBMCs (Figure 1d) and expressed a
range of known microglial marker genes, including P2RY12, CX3CR1 and TMEM119, to a
high level (Extended Data Figure 1a). We defined this population of cells as microglia for the
remainder of our analysis. We found three less common populations of cells that closely
resembled other blood cell types, including NKT cells, monocytes or B-cells that comprised
8.4%, 0.5% and 0.3% of our single cell dataset, respectively. These cell types may reflect
either infiltration of immune cells as a result of blood-brain barrier breakdown or intravascular
contamination within the tissue. In support of the former hypothesis, we also found that the
abundance of infiltrating cells strongly correlated with patient pathology, with trauma patients
in particular enriched (OR=7.6, Fisher exact test P=1.2x10
-155
) (Figure 1d). We also found a
preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for thisthis version posted December 20, 2019. ; https://doi.org/10.1101/2019.12.20.874099doi: bioRxiv preprint
significant effect of age on the abundance of infiltrating cells (3.4% increase per year, Wald
test P=0.014) after adjusting for all known confounding factors, which could reflect blood brain
barrier degeneration over the lifespan (Extended Data Figure 1b).
Within microglia, we found four subpopulations of cells (Figure 2a). Two populations (C and
D) were common in patients with acute brain injury (25-76% of cells) but rare in other
pathologies (<5% of cells). Population B was enriched in tumor patients (OR=4.9, P=7.6x10
-
169
) while population A was most common in control and hydrocephalus patients (Figure 2b,
c). In populations A and B, we observed higher expression of microglial markers including
P2RY12 and CX3CR1. Cells from B, C and D also demonstrated an upregulation of general
immune response and cell activation (IL1B, CD83 & CCL3) (Figure 2d, e; Extended Data
Figure 2; Supplementary Table 1). Cells from C and D exhibited additional upregulation of
acute immune response pathways, including NF-kappa B, STAT3, RUNX1 as well as MHC-I
expression. Population C also showed differential expression of genes associated with stress
induced senescence and DNA damage (HIST1H2BG), populations D expressed genes
associated with cell proliferation (FLT1) and chemotaxis (CCL4, CXCL8, CXCL16), the latter
of which is shared with population B. Population B additionally showed strong upregulation of
catabolic process and metabolism (GPX1) and phagocytosis (TREM2). Our cells partially
overlapped with the transcriptional signatures of disease-associated microglia established in
previous literature (Keren-Shaul et al. 2017; Xue et al. 2014) (Figure 2f, g). Taken together,
these results suggest that our data contain a mix of naive microglia (population A), with three
distinct states of activation that, in part, are driven by patient pathology.
Biological drivers of microglial expression
Our sampling design enabled us to explore the relative importance of a wide range of biological
factors in driving microglial gene expression while controlling for important technical
confounders, using variance components analysis. Of the biological factors we examined,
clinical pathology explained more variation than all other factors combined, although all factors
except sex, including age, brain region, dominant hemisphere and ethnicity, explained a
fraction of variation that was significantly different from zero (LR test FWER 0.05) (Figure 3a).
Patient explained the most variability of any single factor in the model. Although this factor
captured the contribution of genetic background, it is also likely to reflect unmeasured
technical effects, such as variability in cell dissociation and surgical sampling, which are
confounded with patient in the model. The cellular pathways that differed between patient
pathologies closely resembled the differences we observed between different subpopulations
of microglia (Extended Data Figure 3a, b). We also detected 260 genes that varied
significantly by patient age, showing upregulation of inflammation (CLEC7A, CIITA and TLR2)
and downregulation of cell identity (P2RY12, CX3CR1), motility and proliferation (CSF1R) with
increasing age (Figure 3b-e; Extended Data Figure 3c; Supplementary Table 2). Although
sex explained little variation globally, we found 97 genes that were differentially expressed
between males and females (Figure 3f). These included multiple genes in the complement
pathway and synaptic pruning mechanisms (C1QA, C1QC and C3) that were more highly
expressed in females than males (Figure 3g; Extended Data Figure 3d; Supplementary
Table 3). Anatomical region of sampling also had a subtle effect on transcriptional variation,
with cerebellar microglia, which are known to exhibit a distinct, less ramified morphology
upregulating multiple recruitment chemokines (CCL4, CCL3, CCL4L2, CCL3L3) (Figure 3h;
Extended Data Figure 3e; Supplementary Table 4).
preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for thisthis version posted December 20, 2019. ; https://doi.org/10.1101/2019.12.20.874099doi: bioRxiv preprint
eQTL mapping in human microglia and neurodegeneration
We constructed a map of expression quantitative trait loci (eQTLs) in primary human microglia.
After excluding samples with low genotyping quality or substantial non-European ancestry, we
mapped eQTLs using our bulk RNA-seq data from 93 individuals, and detected 401 eQTLs,
summing over hierarchical model posteriors (585 eQTLs at FDR 5% using linear model). The
low number of eQTLs reflected the high between-sample heterogeneity in microglia, compared
with other cell types (Extended Data Figure 4a). We tested for colocalization of risk loci from
18 genome wide association studies (GWAS) with microglia eQTLs (Figure 4a), including five
previous studies of Alzheimer’s disease (AD), and our own meta-analysis of these five studies
for comparison (Online Methods). Across all AD GWAS, we found up to 11 risk loci with a
posterior probability of colocalisation (PP4) greater than 0.5 (Table 1). These included well-
known AD loci, such as BIN1, and less well-studied AD associations, for example EPHA1-
AS1. We repeated the analysis using microglia eQTLs mapped by RASQUAL to support the
colocalisation result using the allele specific expression signature (Supplementary Table 5).
This analysis retrieved colocalisations at other AD GWAS loci, such as CD33 and CASS4.
However, the test statistics may be inflated due to the additional overdispersion in our data
set (Extended Data Figure 4a).
Next, we compared AD risk loci from our meta-analysis with eQTLs from the GTEx project
(v7), in circulating blood monocytes and in a novel dataset of IPSC-derived macrophages
(IPSDMs) from 133 healthy individuals (Online Methods) (Figure 4b). We found more
colocalised AD eQTLs in microglia than in any GTEx brain region. We also observed many
AD risk loci that colocalised with eQTLs in blood monocytes and IPSDMs. To explore the level
of cell-type specificity, we mapped eQTLs jointly analysing data from microglia, monocytes
and IPSDMs using a three-way Bayesian hierarchical model (Extended Data Figure 4a, b;
Online Methods). Using this approach, we discovered 855 eQTLs, of which 108 were
microglia-specific, 449 were found in all three cell types, and 192 were shared with IPSDMs
but not monocytes. We also used the three-way model to evaluate the extent of sharing
between the microglia eQTLs, IPSDM or monocyte eQTLs, and AD risk loci. Many colocalised
AD loci, including BIN1, are found in both microglia and IPSDMs, but absent in monocytes
(Figure 4c, d). There were also multiple AD loci where an eQTL was only detectable in
circulating monocytes (e.g., CASS4 locus), although this is likely to reflect primarily the
differences in power between the monocyte (n=193) and microglia data sets.
IPS models of AD risk loci are an invaluable resource for the development of future
therapeutics. We next identified three AD association signals (BIN1, the EPHA1 locus and
PTK2B) that colocalised with both microglia and IPSDM eQTLs (Figure 4d). The association
for EPHA1-AS1 was shared across many cell types (Extended Data Figure 5a, b), while the
direction of effect at the PTK2B locus was inconsistent (Extended Data Figure 5c). To fine-
map causal variants we generated ATAC-seq data from 5 primary microglia and 89 IPSDMs.
Colocalisation analysis revealed that the AD association signal at BIN1 was highly cell type
specific (Figure 5a). The lead SNP of this association signal, rs6733839:C>T, was located in
a region of open chromatin in both microglia and IPSDMs in which the AD risk allele.
rs6733839:C>T was also associated with a significant change in chromatin accessibility
(Figure 5b, P<6.1x10
-10
), and the association signal for chromatin also colocalised
(PP4=0.996) with the AD association signal (Figure 5c-f). The AD risk allele at
rs6733839:C>T created a predicted high-affinity binding site for the MEF2A transcription factor
preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for thisthis version posted December 20, 2019. ; https://doi.org/10.1101/2019.12.20.874099doi: bioRxiv preprint
(Extended Data Figure 5d). We found that, although BIN1 and MEF2A are broadly expressed
in many tissues, co-expression of both genes was found only in primary microglia and IPSDMs
(Extended Data Figure 5e).
Discussion
Here we present a population-level study of human primary microglia. By sampling cells from
living donors, we defined transcriptional signatures of in vivo microglial activation avoiding
artefacts from post mortem index and cell culture. We identified multiple microglial
subpopulations and showed how these populations are shaped by insult, injury and other life
history factors. We also created the first map of eQTLs in microglia, identified high confidence
causal genes and variants underlying risk loci for Alzheimer’s disease, and identified a subset
that replicated in a scalable IPS model system.
Our results underscore the variability between microglia cells from different individuals. One
implication of the variation we observed between different patient pathologies is that the full
spectrum of microglial function is not well cannot be captured by small studies of a single
patient population. The most obvious example of this are the populations of activated microglia
we identified that account for less than 5% of cells in non-trauma patients. Our results also
provide a picture of the function of microglia following severe trauma, producing cell
populations that exhibit a mixture of a proinflammatory and chemotactic phenotypes. Notably,
although animal models of acute brain injury suggest rapid expansion of microglia following
trauma (Vela et al. 2002), we only observed one population we identified had a proliferative
phenotype, and both showed downregulation of CSF1R. Also in contrast to previous reports
(Olah et al. 2018), we found relatively subtle effects of age on microglial transcription. The
modest changes we did detect were consistent with increased inflammatory senescence in
microglia over the lifespan. Likewise, differences in microglia expression between males and
females were relatively small, although we did observe increased complement activity in
females, perhaps suggesting a role for complement pathways in the higher incidence of AD in
women.
Our eQTL analysis revealed a number of candidate risk genes for AD that function in microglia.
This included well-known genes, such as BIN1, and a number of less well understood loci.
One example we discovered was the EPHA1-AS1 locus, where AD risk appeared to be driven
by a change in the expression of a long noncoding RNA, rather than the neighbouring protein-
coding gene EPHA1 (Extended Data Figure 5a, b). We did not detect some well-known AD
risk loci, such as CD33, with suspected function in myeloid cells. In the case of CD33, analysis
of splicing patterns did indeed reveal a splice QTL at exon 2 (Extended Data Figure 6a-c),
consistent with previous studies (Raj et al. 2014). In other cases, we found strong
colocalisation between AD risk loci and monocyte, but not microglia, eQTLs. While it is
tempting to conclude that this reflects a monocyte-specific function, we believe it is more
plausible that this reflects lower power in our microglia dataset and that, with an increased our
sample size, many of these eQTLs would be found to be shared between the two cell types.
One example of this is the CASS4 locus, where the minor allele frequency of the GWAS lead
variant (rs6014724:A>G) was >10% in the monocyte data, but <5% in our microglia data set
(Extended Data Figure 6d). Other examples of apparently spurious monocyte-specific
preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for thisthis version posted December 20, 2019. ; https://doi.org/10.1101/2019.12.20.874099doi: bioRxiv preprint