Untangling the brain's neuroinflammatory and neurodegenerative transcriptional responses

A common approach to understanding neurodegenerative disease is comparing gene expression in diseased versus healthy tissues. We illustrate that expression profiles derived from whole tissue RNA highly reflect the degenerating tissues' altered cellular composition, not necessarily transcriptional regulation. To accurately understand transcriptional changes that accompany neuropathology, we acutely purify neurons, astrocytes and microglia from single adult mouse brains and analyse their transcriptomes by RNA sequencing. Using peripheral endotoxemia to establish the method, we reveal highly specific transcriptional responses and altered RNA processing in each cell type, with Tnfr1 required for the astrocytic response. Extending the method to an Alzheimer's disease model, we confirm that transcriptomic changes observed in whole tissue are driven primarily by cell type composition, not transcriptional regulation, and identify hundreds of cell type-specific changes undetected in whole tissue RNA. Applying similar methods to additional models and patient tissues will transform our understanding of aberrant gene expression in neurological disease.

A common approach to understanding neurodegenerative disease is comparing gene expression in diseased versus healthy tissues. We illustrate that expression profiles derived from whole tissue RNA highly reflect the degenerating tissues' altered cellular composition, not necessarily transcriptional regulation. To accurately understand transcriptional changes that accompany neuropathology, we acutely purify neurons, astrocytes and microglia from single adult mouse brains and analyse their transcriptomes by RNA sequencing. Using peripheral endotoxemia to establish the method, we reveal highly specific transcriptional responses and altered RNA processing in each cell type, with Tnfr1 required for the astrocytic response. Extending the method to an Alzheimer's disease model, we confirm that transcriptomic changes observed in whole tissue are driven primarily by cell type composition, not transcriptional regulation, and identify hundreds of cell type-specific changes undetected in whole tissue RNA. Applying similar methods to additional models and patient tissues will transform our understanding of aberrant gene expression in neurological disease.
O ne approach to better understand the molecular mechanisms of neurodegenerative disease is to compare gene expression profiles from diseased versus control tissues and draw inferences about which biological pathways and cellular processes are altered in the disease state. However, the cellular complexity of central nervous system (CNS) tissue, in which glial cell types including microglia and astrocytes are interspersed among neurons of many subtypes, limits the utility of this approach. Expression profiles derived from whole tissue RNA represent each gene's average expression among all cells but do not reveal which cell types are responsible for a gene's normal or altered expression in healthy or diseased tissues. Lacking such information, the genes and pathways implicated by profiling whole tissues cannot be readily incorporated into cellular models of neurodegenerative disease. Moreover, changes in a gene's expression that occur in a specific cell type may be undetected in whole tissue RNA as the difference may be masked by the overall signal from all cell types.
To circumvent these shortcomings, researchers have developed methods to acutely isolate individual cell types from adult brain tissue. Most commonly, brain tissue is dissociated into single cells, from which microglial/macrophage-type cells-specifically labelled genetically (for example, Cx3cr1::GFP expression) or biochemically (for example, anti-CD11b)-are purified by fluorescence-activated cell sorting (FACS) or other antibodybased methods 1,2 . Using similar methods, researchers have also isolated astrocytes, neurons, endothelial cells and other brain cell types [3][4][5][6] , yet these significant advances have certain limitations. First, most dissociation methods involve enzymatic treatment at warm or ambient temperatures 1,7-9 , allowing stress-induced changes in RNA profiles to occur throughout the procedure. Second, genetic labelling methods require extra resources and time to obtain the desired cell type labelled at the appropriate disease stage and in the proper genetic background, and may also interfere with normal biology 10,11 . Third, researchers often focus on a cell type of particular interest rather than study multiple cell types from the same brain, so correlative cell type analyses within specimens cannot be performed. Fourth, samples are often pooled to increase RNA yield and detection, obscuring animal-to-animal variability and increasing the required number of specimens. Fifth, many gene expression studies have used microarrays or other technologies that are becoming outmoded by the advent of high-throughput RNA sequencing, which has enabled more comprehensive transcriptomic analyses.
Here we utilize an approach that avoids some of the abovementioned limitations 12 and adapt it further to isolate populations of neurons, astrocytes and microglia from single adult brain specimens and analyse their transcriptomes by RNA amplification and sequencing. To our knowledge, this is the first report of these three cell populations being purified simultaneously from the brain of an adult mouse and analysed by RNA sequencing (RNA-Seq). The method does not require incubations at warm temperatures for enzymatic dissociation, genetic labelling of any cell type or pooling of samples. Using peripheral endotoxemia as an acute neuroinflammatory model to establish the method's utility, we demonstrate the diversity and specificity of each cell type's transcriptional and RNA processing responses. We observe correlations in animal-toanimal variability between cell types and investigate the tumournecrosis factor (TNF) pathway's contribution to the brain's endotoxemia response.
We also use cell type-specific sequencing data to probe existing data sets of gene expression in neurodegenerative disease tissues from human patients and/or animal models of frontotemporal dementia (FTD), amyotrophic lateral sclerosis (ALS) and Alzheimer's disease (AD). We provide evidence that disease-related changes in expression profiles from whole tissue RNA are often not due to transcriptional regulation but rather the altered cell type composition of disease tissue samples. Finally, we apply our sorting method to a mouse model of AD to reveal many pathology-related transcriptional changes in purified cell types that were not observed in whole tissue RNA.

Results
Altered cellularity in neurodegenerative expression profiles. To better understand existing data sets for gene expression in neurodegenerative diseases, we used recent RNA-Seq data from different cell types acutely isolated from postnatal mouse brains 9 to predict which cell types were most likely responsible for disease-related changes observed in whole tissue. We also used our own RNA-Seq data from neurons, astrocytes and microglia acutely isolated from adult mouse brain (see next section for details) to further inform the results.
We first looked at several hundred genes whose expression changed at least twofold (adjusted Pr0.05) in a microarray analysis of frontal cortex samples from FTD patients with mutations in the progranulin (GRN) gene, relative to non-diseased control samples 13 . When we mapped each gene to the cell type(s) most likely responsible for its expression, the overall trend was striking. Most of the transcripts with increased abundance in FTD cortex mapped to non-neuronal cell typesastrocytes and microglia, as well as endothelial cells-whereas most of the transcripts with decreased abundance mapped to neurons (Fig. 1a,b). This suggested that many of the changes in gene expression were influenced by the altered cellular composition of end stage FTD tissue, with extensive neuronal loss and prominent gliosis.
We next applied this type of analysis to our microarray data from the PS2APP transgenic AD model, in which mutant alleles of APP and PSEN2 give rise to age-dependent amyloid plaque pathology 14 . In 13-month aged mice, 53 genes (including Thy1 and Prnp transgene elements) showed Z2-fold increased expression (adjusted Pr0.05) in PS2APP versus non-transgenic cortex. Most of these differentially expressed genes were normally expressed by microglia, and some by astrocytes (Fig. 1c,d). Unlike in FTD cortex, we did not see Z2-fold downregulation of any neuronal genes (or any genes) in PS2APP cortex. These observations were consistent with histopathological features of PS2APP brains, including plaque-associated gliosis but no notable loss of neurons, suggesting that many instances of increased gene expression in PS2APP cortex were due, at least in part, to increased glial cell numbers rather than transcriptional upregulation.
Finally, we performed similar analyses on data sets from a murine ALS model and from ALS patients. RNA samples from intact spinal cords of mice expressing mutant human SOD1 (hSOD1 mut ) showed hundreds or dozens of genes whose expression was increased or decreased, respectively, at least twofold (adjusted Pr0.05) compared with samples from non-transgenic mice 15 . As with FTD samples, most of the transcripts with increased abundance in hSOD1 mut mice mapped to astrocytes and microglia, as well as endothelial cells, whereas the majority of transcripts with decreased abundance mapped to neurons (Fig. 2a,b)-likely reflecting the gliosis and neuronal loss that progressively occur with ageing in hSOD1 mut spinal cords. Consistent with this interpretation, most expression differences between hSOD1 mut and non-transgenic spinal cords magnified as the mice aged (Fig. 2a). Turning our attention to the human disease, we examined gene expression data from laser-captured motor neurons microdissected from spinal cords of ALS or control patients 16 . Surprisingly, most of the 26 genes with Z2-fold increased RNA abundance (adjusted Pr0.05) in ALS motor neuron samples were predicted to be expressed in microglia or astrocytes rather than in neurons (Fig. 2c,d). Thus, despite meticulous efforts to identify transcriptional changes in ALS motor neurons, most of the observed differences may have derived from contaminating gliotic cells, or fragments thereof, co-excised with the neurons. Moreover, no genes met the cutoff of Z2-fold decreased expression in human ALS motor neuron samples, further implying that depletion of certain neuronal RNAs in hSOD1 mut spinal cords reflects neuronal loss, whereas gene expression in extant motor neurons is relatively unchanged.
To further explore how cellular composition may influence neurodegenerative disease transcriptional profiles, we took a complementary approach by looking at whether RNAs specifically expressed by microglia or astrocytes in normal tissue showed overall enrichment in disease tissues, and whether neuron-specific RNAs showed overall depletion (Fig. 3). Microglia-specific RNAs were indeed enriched in cortex from GRN mutant FTD patients, spinal cord from hSOD1 mut mice and cortex from PS2APP mice. Astrocyte-specific RNAs were similarly enriched in human FTD cortex and mouse hSOD1 mut spinal cord, but not in mouse PS2APP cortex. In contrast, neuron-specific RNAs were clearly depleted in FTD cortex and hSOD1 mut spinal cord, whereas depletion was minimal in PS2APP cortex. These observations were perfectly consistent with known histopathology in these tissues-neuronal loss and gliosis in FTD and hSOD1 mut -driven ALS models, and gliosis, but little if any neuronal loss, in amyloid-based AD models. The astrogliosis observed in AD models and in human AD tissue involves morphological changes but not cellular proliferation, whereas microglia clearly proliferate 17,18 . Expression profiles in neurodegenerative disease tissues are thus considerably influenced by altered cellular composition, and profiles from whole tissue RNA may be misleading if changes in disease are assumed to result from transcriptional modulation.
Purifying neurons, astrocytes and microglia from adult brain. The above observations illustrate the need to isolate specific cell types in order to understand transcriptional changes in diseased or injured CNS tissue. Using the method of Guez-Barber et al. 12 for purifying neurons or endothelial cells from rat brain as a starting point, we set out to dissociate and immunopurify cortical neurons, astrocytes and microglia from a single adult mouse brain, with no genetic cell type labelling, for analysis by RNA sequencing. To avoid stress-related transcriptional responses that occur during warm enzymatic incubations, we performed all procedures at 4°C or on ice, including Accutase treatment for 30 min followed by mechanical dissociation. The suspension was filtered (70 mm) and then centrifuged through Percoll to enrich for cells and reduce debris. The pelleted fraction was resuspended, subjected to mild ethanol fixation and immunolabelled using anti-CD11b to label microglia, anti-glial fibrillary acidic protein (GFAP) to label astrocytes, and anti-NeuN to label neurons. populations were collected by FACS (Fig. 4a), followed by RNA extraction from purified cell types. The RNA levels for the immunolabelling targets CD11b (Itgam), GFAP and NeuN (Rbfox3) were appropriately enriched or depleted in each sorted cell population when analysed by reverse transcription-quantitative PCR (RT-qPCR; Supplementary Fig. 1a). To verify appropriate cell type enrichment in sorted populations, we also analysed additional RNAs encoding the cell type markers Cx3cr1, Iba1 (Aif1) and CD45 (Ptprc) for microglia; Aldh1l1, Aqp4 and GLT1 (Slc1a2) for astrocytes; and Grin1, Cx3cl1 and nNOS (Nos1) for neurons. RNA for each marker was enriched in the expected cell population relative to the other two populations (Supplementary Fig. 1b), confirming the successful purification of microglia, astrocytes and neurons from adult cortex.

Cellular specificity of CNS endotoxic transcriptional response.
To test the method's suitability for detecting treatment-dependent, genome-wide, transcriptional changes within each cell type, we injected mice peripherally with saline or with the endotoxin lipopolysaccharide (LPS), which induces an inflammatory response in the CNS [19][20][21] . One day post injection, we dissociated and FACS-purified cell populations from perfused cortex, followed by RNA purification, amplification and sequencing analysis. RNA-Seq data showed appropriate marker enrichment in sorted cell populations (Fig. 4b), confirming our earlier RT-qPCR data. The expression profiles we obtained for microglia, astrocytes and neurons from control adult mice were generally similar to profiles reported elsewhere for the same cell types isolated from postnatal mouse brains using different methods 9 (see Supplementary Data 1 for details). However, our data set may better reflect unperturbed microglia, as certain LPS-induced genes in our study including Cst7, Cd300lf and Ccl5 were expressed only at background levels (similar to neurons) in microglia from saline-injected mice but showed high microgliaspecific expression in the data set of Zhang et al. 9 ( Supplementary  Fig. 2a), possibly due to transcriptional induction during their enzymatic treatment of tissues at warm temperature. In addition, our neuron expression data reflect mature rather than immature stages of development. Thus, certain genes highly expressed in adult but not fetal human cortex 22 show neuron-specific expression in our data from adult mice but not in the data of Zhang et al. from postnatal mice. Conversely, certain genes highly expressed in fetal but not adult human cortex show neuronspecific expression in the data of Zhang et al., but not in our adult data ( Supplementary Fig. 2b).
Unsupervised clustering of our samples using the 2.5% most variably expressed genes segregated the samples into distinct groups according to cell type (Fig. 5a). The peripheral endotoxemia had a profound effect on gene expression in both microglia and astrocytes, as samples within these cell types clustered naturally into treatment and control groups. Neuron samples did not cluster into separate treatment and control groups at this level of analysis, indicating that the brain's most pronounced transcriptional responses to endotoxemia occurred within glial cell types. Next, we plotted the mean expression level for every gene in each cell type from LPS-injected versus saline-injected mice ( Fig. 5b; Fig. 1a,b, this heat map (colour range: À 3.34rZr3.68) and triangle plot analyse cell type expression for genes differentially expressed in bulk spinal cords from the SOD1 G93A ALS model versus wild-type mice at 18 weeks of age (GSE18597, imported). (c,d) Similar to Fig. 1a,b, this heat map (colour range: À 2.75rZr3.14) and triangle plot analyse cell type expression for genes differentially expressed in laser-captured motor neuron samples from ALS patient spinal cords versus controls (GSE18920, imported). Expression data from anterior horn material remaining after motor neuron excision are also shown. a/astro., astrocytes; Cx, cortex; e, endothelial cells; m/micro., microglia; mo, myelinating oligodendrocytes; n, neurons; nfo, newly formed oligodendrocytes; opc, oligodendrocyte progenitor cells; wt, wild type. transcriptional response, in contrast to glial cells. Microglia and astrocytes each showed more than a thousand genes with Z4-fold change in expression (adjusted Pr0.05), compared with only 75 such changes in neurons. Microglia and astrocytes are both widely known to mediate CNS inflammatory signalling, but the endotoxic transcriptional responses in these cell types were remarkably distinct, with only B15% of the genes with altered expression in either cell type being similarly regulated in the other ( Fig. 6a; see also Supplementary Data 2 to compare LPS responses of these cell types with that of neurons). Even among genes similarly affected in both cell types, some were primarily expressed by one cell type in terms of absolute expression level, further emphasizing cellular specificity in the brain's endotoxic transcriptional response (Fig. 6b). Compared with previous reports of astrocytic 3 and microglial 23 responses to endotoxemia, our data revealed hundreds of additional LPS-responsive genes in each cell type ( Supplementary Figs 3 and 4).
An attractive prospect of this method is detecting changes within individual cell types that are masked in whole tissue RNA. To demonstrate this advantage and to replicate our initial findings using a different technology, we selected several LPS-responsive genes from our sorted cell RNA-Seq data that seemed unlikely to appear as LPS-responsive in whole tissue RNA because of constitutive expression in other cell types. We injected new groups of mice, and 24 h later purified whole tissue RNA from one hemicortex and individual cell types from the other. Testing the RNA samples by RT-qPCR validated our RNA-Seq data and verified our predictions, confirming that genes including Gas7, Kcna3, Adora2a and Ppargc1b were LPS-induced in microglia; that genes including Nav1, Vipr1, St6gal1, Sema4g, Dock4 and Plxna4 were LPS-repressed in microglia; and that all of these changes were obscured in whole tissue RNA ( Supplementary Fig. 5). In other examples, Arap3 and Slc7a8 expression increased in neurons and astrocytes from LPS-injected mice but not in whole tissue, being offset by endotoxic repression in microglia ( Supplementary Fig. 5).
Another advantage of this method is the ability to correlate expression in different cell types from the same animal. For example, we found that the animal LPS6, and to a lesser extent LPS5, showed a weak LPS response in multiple cell types ( Supplementary Fig. 6). Thus, this method turns what would otherwise be isolated outliers into a consistent signal-animals in which the global LPS response, for whatever reason, was attenuated in multiple cell types.
Brain RNA processing altered in endotoxemia. RNA biology is increasingly recognized as a potential contributor to the aetiology of neurological disease 24 . To demonstrate the feasibility of  Genes specifically expressed by microglia, astrocytes or neurons were selected as those with Z20-fold enrichment compared with the other two cell types when isolated from normal adult mouse cortex and analysed by RNA-seq. For the subsets of these genes or their human homologues that were probed by microarray, distributions are shown for fold-changes in bulk cortex from FTD patients, bulk spinal cord from murine ALS model or bulk cortex from murine AD model, relative to controls, in unaffected (blue) and affected (red) regions and times. The shift of a red distribution to the right is consistent with higher content of that cell type in the affected tissue (for example, microglia in AD model and astrocytes in ALS model), whereas a shift to the left suggests a depletion of that cell type (for example, neurons in FTD cortex). NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11295 ARTICLE NATURE COMMUNICATIONS | 7:11295 | DOI: 10.1038/ncomms11295 | www.nature.com/naturecommunications profiling differential RNA processing in CNS cell types, we used our RNA-Seq data to quantify the inclusion frequencies of annotated and novel RNA processing events 25 and test for LPS-induced changes in those events 26 . To our knowledge, little is known about the effect of endotoxemia on RNA processing in the CNS.
Surprisingly, we identified significant changes for hundreds of known and novel events in all three cell types (Supplementary Data 3). In contrast with gene expression changes, neurons showed as many RNA processing changes as the glial cell types. Many of these changes were in common between the cell types, suggesting a partially shared RNA processing response to peripheral endotoxemia among CNS cell types. For example, we identified a differentially spliced cassette exon in the Sltm gene ( Fig. 7a), which encodes a protein suggested to play a role in global transcriptional regulation as well as apoptosis 27 . In all three cell types from vehicle-treated animals, nearly all of the splice junction reads observed near this exon aligned to the inclusion isoform, whereas in LPS-treated animals many reads aligned to the skipping isoform (Fig. 7a). The extent of skipping in LPS-treated animals was variable and correlated with the overall transcriptional response as the LPS6 animal, which exhibited the weakest LPS response by gene expression ( Supplementary Fig. 6), showed no skipping of this exon in either microglia or neurons (Fig. 7b).
To validate these findings, we selected 27 LPS-induced splicing events that were significant in one or more cell types and for which we could design qPCR assays (see the Methods for details), including the Sltm event, and we analysed a new cohort of vehicle-and LPS-treated animals. For each event, we designed assays targeting the flanking constitutive region (C) as well as the skipping (S) and inclusion (I) isoforms. Using this new cohort and different technology, the overall pattern of LPS-induced splicing for these 27 events was validated in all three cell types ( , and neurons and microglia were isolated as two distinct populations corresponding to NeuN þ CD11band NeuN À CD11b þ cells, respectively (third panel). The NeuN À population was further used to isolate GFAP þ astrocytes (last panel). (b) RNA-Seq data confirmed that NeuN þ sorting enriched for cells expressing neuronal markers, GFAP þ sorting enriched for cells expressing astrocytic markers, and CD11b þ sorting enriched for cells expressing microglial markers (n ¼ 5 animals; one astrocyte sample was excluded from the analysis due to evidence of neuronal contamination). Bars represent mean ± s.d. (Prism).
same trend by RT-qPCR in microglia (13 of 13 events), neurons (17 of 17 events) and astrocytes (16 of 18 events), with 85, 76 and 17% of the events reaching significance in their respective cell types by RT-qPCR (adjusted Pr0.05). We excluded the possibility that our findings were an artefact of our dissociation or fixation methods by confirming that most of the events were LPS-induced even in whole cortex RNA extracted directly from intact, opposite hemisphere of the same brains (see Supplementary Data 4 for details). Thus, peripheral endotoxemia produces a robust and specific RNA processing response in the CNS.
TNFR1 required for astrocytic response to endotoxemia. To demonstrate this method's utility for dissecting CNS inflammatory signalling mechanisms, we explored whether the TNF pathway is required for the brain's endotoxic response in each cell type. Tnf was expressed at low levels in control microglia and escalated in response to endotoxemia (Fig. 8a). As for the receptors, Tnfrsf1a (Tnfr1) was expressed most prominently in microglia and at a lower level in astrocytes, which increased during endotoxemia to a level approximating that of microglia. Tnfrsf1b (Tnfr2) expression was highly enriched only in microglia. Notably, neuronal expression of Tnfr1 and Tnfr2 was negligible, suggesting that reported neurotoxic and neuroprotective activities of Tnfr1 and Tnfr2, respectively 28,29 , likely occur through indirect signalling mechanisms. The induction of Tnf in microglia led us to test whether TNF receptors contributed to the brain's endotoxic transcriptional response. We injected LPS into wild type, Tnfrsf1a À / À , Tnfrsf1b À / À or Tnfrsf1a À / À ;Tnfrsf1b À / À mice (colony mates); FACS-purified cortical cell types 24 h post injection; and analysed expression by RT-qPCR for a set of genes with a variety of LPS responses in different cell types. The microglial response was largely Tnfr-independent, as only a few of the genes we analysed showed impaired responses in microglia lacking Tnfr1 or Tnfr2 (Fig. 8b). In contrast, astrocytes from mice lacking Tnfr1 were markedly impaired in their transcriptional response to endotoxemia, whereas this response was unimpaired in mice lacking only Tnfr2 (Fig. 8b). It is tempting to speculate that the Tnfr1-dependent astrocytic response resulted from direct TNF/Tnfr1 signalling in astrocytes, which express Tnfr1 but not Tnfr2. However, some of the dependence on Tnfr1 may be indirect since neurons, in which Tnfr1 expression is negligible, also displayed certain LPS-induced changes that were Tnfr1-dependent (data not shown). Data from the Tnfr1/Tnfr2 double knockout mice were highly similar to the Tnfr1 knockout data for all cell types. In summary, the brain's transcriptional response to endotoxemia showed a strong requirement for Tnfr1 in the case of astrocytes but not microglia and was Tnfr2-independent in all cell types tested.
Differential expression by cell type in Alzheimer's model. Finally, we explored the utility of this method in the context of neurodegenerative disease by obtaining RNA-Seq data from both bulk tissue and FACS-purified cell types from 7-and 13-monthold mouse brains in the PS2APP model of AD (Fig. 9a). In 13-month-old PS2APP whole cortex, we detected 85 genes (in addition to the transgene elements) with Z2-fold increased RNA abundance (adjusted Pr0.05)-a larger number than the 51 genes detected by microarray (Fig. 1c). Of these 85 genes, 66 showed specific or highly enriched expression in microglia compared with astrocytes and neurons sorted from either non-transgenic or PS2APP cortex of that age. However, only a fraction of these genes (13 out of 66) showed increased expression (adjusted Pr0.05) in PS2APP microglia (for example, Ccl3, Clec7a, Treml2), whereas the remainder showed no clear changes (for example, C1qc, Cybb, Gpr84, Ptprc, Trem2), with 20 out of 53 trending downward and 33 out of 53 trending upward (see Fig. 9b for examples). These data directly support our earlier conclusions that the changes in RNA profiles observed in whole tissue samples for neurodegenerative diseases are driven more by relative cell type abundance than by disease-related transcriptional modulation.
Of the remaining 19 genes with Z2-fold higher detection in PS2APP whole cortex, microglia were responsible for 8 of them, including Gpnmb and Mamdc2, which normally show astrocyteenriched expression (see Gpnmb in Fig. 9c). In all, microglia accounted for 87% (74 out of 85) of the changes observed in whole tissue, because of increased microglial cell number, increased expression levels within microglia or both. Astrocytes accounted for most (8 out of 11) of the remaining changesroughly 10% of total-with 5 genes including Gfap, Bcl3 and complement 4 reaching significance (Fig. 9c) and 3 trending strongly upward in sorted PS2APP astrocytes (not shown). Three changes in whole tissue could not be conclusively assigned to any cell type, but all three were expressed selectively by glia. Of note, neurons accounted for none of the expression changes observed in PS2APP tissue except for that of the transgene. (PS2APP neurons as well as microglia showed elevated Gfap expression, with fold changes higher than in astrocytes, but the increased Gfap signal in whole tissue was assigned primarily to astrocytes due to their higher absolute expression level.) More interestingly, in sorted PS2APP microglia we observed more than 200 genes with Z2-fold change in expression (adjusted Pr0.05) that were undetected or less apparent in whole tissue (Fig. 10a). Many of these were constitutively expressed in neurons and/or astrocytes, explaining why their upregulation in microglia was difficult to observe in whole cortex (Fig. 10b). Increased expression of Apoe, Lpl and Ldlrad3 in PS2APP microglia (example in Fig. 10c) suggested that microglia No. of events tested by qPCR   ARTICLE from plaque-ridden cortical tissue have altered lipoprotein metabolism and may contribute to pathogenic amyloid precursor protein (APP) processing and amyloid deposition [30][31][32][33] . PS2APP microglia also upregulated the genes encoding Igf1, the transcription factor Hif1a and its partner Arnt2 (ref. 34), and the transcriptional repressor Bhlhe40 and its partner Arntl (ref. 35; example in Fig. 10c). Finally, among the genes with altered expression in PS2APP microglia were 34 instances of genes with Z2-fold decreased expression. These included microglia-specific transcripts such as Clec4a3, Irf4, Mgl2, Mrc1 and Tgm2, none of which showed decreased abundance in whole tissue RNA as the decreases in expression were offset by increased microglial cell number in PS2APP cortex (example in Fig. 10c). These examples clearly illustrate the advantageous view of neurodegenerative disease that cell type-specific expression profiles can provide, which cannot be appreciated using whole tissue expression profiles. Looking at microglial gene expression more broadly, we compared the pattern of differential expression in purified PS2APP cortical microglia with that of microglia purified from spinal cords of hSOD1 mut mice 23 to look for evidence of similar gene regulation. Despite the unique aetiologies and locations of the two pathologies, and the fact that little to no neuronal loss occurs in PS2APP brains, we observed significant overlap between the 'neurodegeneration signature' described by Chiu et al. and the differential expression pattern in PS2APP microglia. Over 40% of the upregulated genes (Z2-fold, Pr0.05) in PS2APP microglia were also upregulated in hSOD1 mut microglia (Fig. 10d). Conversely, only 11% of the genes upregulated in hSOD1 mut microglia were also upregulated in PS2APP microglia, presumably reflecting the more aggressive pathological and  Lag3 Plek Lgals3bp .05) in 13-month-old PS2APP whole cortex (GSE75357, this study) to expression data for astrocytes, microglia and neurons acutely isolated from aged transgenic (Tg) and wild-type (wt) mice (GSE75431, this study). Most genes with increased abundance in PS2APP cortex were preferentially expressed by microglia, but only a fraction of these were transcriptionally upregulated in microglia. Expression was Z-score normalized for each gene, separately within each of the two data sets (colour range: À 2.03rZr4.00). One astrocyte sample from PS2APP mice was excluded from the analysis due to evidence of cell type contamination. (b,c) Expression patterns of indicated genes analysed by RNA-Seq for whole cortex and for sorted cells, from PS2APP versus non-transgenic cortices, with FC and adjusted P-values (Wald test, DESeq2; n ¼ 5 animals per genotype) shown and bars representing mean±s.d. (Prism). Note: Absolute nRPKM values between whole cortex and sorted cell studies should not be directly compared, due to different library preparation methods.
neurodegenerative phenotypes in the ALS model compared with the AD model.

Discussion
Studying whole tissue transcriptomes from neurodegenerative disease specimens has frequently led to the conclusion that inflammation-related pathways are active in the disease. However, whether these changes are causal or consequential in the disease process has not been clear. Our examination of existing data sets from neurodegenerative disease tissues revealed a common trend-that genes with decreased mRNA abundance were typically expressed by neurons, and genes with increased mRNA abundance were typically expressed by microglia or astrocytes-suggesting that many changes resulted not from transcriptional modulation but from neuronal loss and the attendant increase in glial cell fraction. Supporting this interpretation, we observed bulk enrichment of microglia-and astrocyte-specific transcripts and depletion of neuron-specific transcripts in these tissue profiles. We also revealed potential limitations of laser-capture microdissection, as most of the transcripts with increased abundance in ALS motor neuron samples likely originated from gliotic cell fragments that were unintentionally co-excised with the neurons.
To illuminate the brain's cellular mechanisms of inflammatory signalling, we improved on a method to dissociate and immunopurify CNS cell types from adult mouse brain, followed by RNA purification, amplification and sequencing. The method avoids unwanted changes to the transcriptome that occur during enzymatic dissociation at warm temperatures, and it does not require genetic labelling of any cell type. Fixation after dissociation enables the labelling of intracellular targets and prevents additional changes to RNA that may occur during labelling and sorting. No pooling of samples is required, so treatment responses can be correlated between cell types within individual animals and compared between animals. Using peripheral endotoxemia as a model to prove this method, we demonstrated that the brain's responses in gene expression and RNA processing are exquisitely cell type specific, and we observed myriad changes that would go undetected without cell sorting.
A caveat of the method is whether the population of cells recovered after dissociation and sorting adequately represents the starting population. The number of cells that survive the harsh dissociation procedure is a fraction of the cells originally present; it seems likely that subpopulations within each cell type will have different liabilities. Moreover, only cell body-contained transcripts will be recovered, as cellular extensions such as axons and dendrites are sheared during mechanical dissociation. Another consideration is immunolabelling specificity-NeuN is expressed in most, but not all, neurons; CD11b is also expressed by peripheral myeloid cells that may infiltrate the brain; some astrocytes express little GFAP under normal conditions. We recovered similar numbers of GFAP þ cells from saline and LPS-injected animals, arguing against the notion that LPS-induced changes in astrocyte gene expression resulted from additional subtypes being collected after GFAP upregulation (data not shown). Numbers of CD11b þ cells recovered after LPS injection were also similar, although a small percentage was CD45 high , suggesting a slight contribution from peripheral myeloid cells (data not shown).
The mechanisms by which peripheral endotoxin induces a CNS transcriptional response are imprecisely understood. A pervasive misconception is that astrocytes can directly respond to LPS exposure [36][37][38] , but recent data have clarified that microglia are critical intermediaries in the astrocytic response 39 ; our expression data confirm that astrocytic expression of the LPS co-receptors Tlr4 and CD14 is extremely low (not shown). The intermediary signals between stimulation of immune cells by LPS and the endotoxic astrocyte response are unknown, but using our method we established a requirement for Tnfr1 in the indirect process of LPS-induced astrocyte activation. The microglial response was independent of TNF receptors, consistent with the possibility of direct exposure to LPS following LPS-triggered permeabilization of the blood-brain barrier 40 .
Extending the method to a murine AD model, we verified that the gene expression changes observed in whole cortex RNA resulted mostly from increased numbers of microglia in PS2APP cortex, sometimes in combination with increased expression within microglia. This was yet another indication that many of the 'hits' identified in whole tissue expression data are trivial, unavoidable consequences of altered cell type composition. Importantly, the vast majority of transcripts with altered expression in sorted PS2APP microglia were not significantly altered in whole tissue RNA, confirming the importance of obtaining cell type-specific expression data for revealing diseaserelated transcriptional changes. Comparing our data with published data from an ALS mouse model revealed that many of the microglial expression changes were conserved between models of ALS and AD, suggesting a common microglial response to distinct disease pathologies. The approach we have described will allow researchers to identify bona fide instances of disease-related gene regulation, pinpoint the cell types in which those change occur and formulate better hypotheses for the cellular mechanisms of neurodegenerative disease. Applying this method to additional models and patient tissues will potentially transform our understanding of the role of aberrant gene expression in neurological disease or injury.

Methods
Mice. All protocols involving animals were approved by Genentech's Institutional Animal Care and Use Committee, in accordance with guidelines that adhere to and exceed state and national ethical regulations for animal care and use in research. Eight-to ten-week-old C57BL/6 males were injected intraperitoneally with LPS (0111:B4 L3024, Sigma) at 10 mg kg À 1 . Control animals were injected with an equal volume of PBS. Animals were monitored for 24 h and perfused with PBS before dissection. Animals were perfused and processed in pairs of one LPS-and one PBS-injected mouse at a time, and sorted in two BD FACSAria sorters simultaneously. TNF receptor knockout mice were obtained from Jackson (stock no. 003243) and maintained as a colony to ensure production of both knockout and wild-type mice. For cell sorting studies using PS2APP 14 model, female cohorts (except one male in 13-mo non-transgenic group and two males in 13-mo transgenic group) were perfused and processed at 7 and 13 months, usually in pairs of same-age transgenic and non-transgenic mice. For microarray studies using PS2APP model, male cohorts were perfused and processed for whole hemicortex RNA purification at 3, 7 and 13 months. For bulk tissue RNA-Seq studies using PS2APP, female cohorts were similarly processed at 7 and 13 months. In general, using 4-6 mice per treatment group or genotype was enough animals to detect robust changes in gene expression as statistically different.
Preparing dissociated cell suspensions from adult mouse brain. Immediately following perfusion with PBS, cortical caps (including hippocampi) were rapidly dissected on ice, minced into smaller pieces using a fresh chilled razor blade, and transferred to a round-bottom 2 ml Eppendorf tube containing 1.6 ml of Accutase (SCR005, Millipore). The tissue was rotated at 4°C for 20 min and spun at 2,000 relative centrifugal force (r.c.f) for 1 min in a refrigerated centrifuge. The supernatant was discarded and 1.5 ml Hibernate A Low Fluorescence medium (Brain-Bits) was added to the tissue. The tissue was manually and carefully triturated using first a blunted 1 ml tip five times followed by an uncut 1 ml tip seven or eight times. The cell suspension was kept on ice for 10-15 s to allow the larger tissue to settle. The top cloudy 1 ml cell suspension was removed and passed through a 70-mm prewetted filter. A fresh 1 ml of Hibernate A was then added to the cell suspension and the trituration was repeated until a total of 4 ml of cell suspension was obtained. This cell suspension was pipetted onto a discontinuous Percoll gradient (described in Guez-Barber, et al. 12 ) in a 15-ml Falcon tube, high density on bottom, and centrifuged at 430 r.c.f. for 4 min in a tabletop centrifuge, followed by removal and discarding of the top 2 ml of cloudy liquid. The remaining 5 ml was centrifuged at 550 r.c.f. for 4 min. The resulting supernatant was discarded and the cell pellet was resuspended gently in 1 ml of ice cold Hibernate A. One millilitre of ice-cold 100% ethanol was added drop by drop while gently swirling the falcon tube containing cells on ice to uniformly and gently fix the cells. After the cells were mixed to a final concentration of 50% ethanol, the falcon tube was kept on ice for 15 min. The fixed cells were centrifuged at 550 r.c.f. for 5 min, and the resulting supernatant was discarded. The cell pellet was resuspended gently in 10 ml of ice-cold Hibernate A to wash residual ethanol, centrifugation was repeated and cells were finally resuspended in 1 ml of ice-cold Hibernate A and aliquoted for immunostaining (50 ml cells þ 950 ml Hibernate A for isotype and unstained control cells; remaining 900 ml of cells þ 100 ml of Hibernate A for antibody staining).
Dissociated cell immunostaining, FACS and RNA isolation. The following antibodies were used for immmunostaining neurons, astroctyes and microglia, respectively: AlexaFlour 488-conjugated anti-NeuN (MAB377X, Millipore) at 1:1,000, PE-conjugated anti-GFAP (561483, BD Biosciences) at 1:50 and APC-conjugated anti-CD11b (561690, BD Biosciences) at 1:250. 1 ml of DAPI (1 mg ml À 1 stock) was added to all the tubes. The cells were rotated at 4°C for 20 min. Cells were centrifuged at 2,000 r.c.f. for 2 min in a refrigerated centrifuge followed by two brief washes in 1 ml of Hibernate A. Cells were resuspended in 3 ml of Hibernate A and passed through a 40-mm filter before sorting. Cells were sorted using DAPI þ signal as a gate to select for singlet cell bodies and ignore doublets and cell debris. Using this group as the parent gate, further FACS gates were determined relative to isotype controls. CD11b þ cells were distinct from other cell types. CD11bcells separated into distinct NeuN þ and NeuNclusters, and neurons were gated from the NeuN þ cluster. Astrocytes were gated from the GFAP high cells that were also NeuN-negative. All cell types were finally gated against other fluorescent gates to get the final population. In 13 months PS2APP versus non-transgenic mice, the CD11b þ population included cells with increased NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11295 ARTICLE NATURE COMMUNICATIONS | 7:11295 | DOI: 10.1038/ncomms11295 | www.nature.com/naturecommunications GFAP þ signal (although still lower than astrocytic GFAP þ signal) but these were gated out to ensure that collected microglia did not include contaminating astrocytic material. Cells were sorted into protein LoBind tubes (022431081, Eppendorf) containing 50 ml of Hibernate A medium. Cells were centrifuged at 2,600 r.c.f. for 7 min (neurons) or 6,600 r.c.f. for 7 min (other cell types) and the supernatant was discarded. RNA was extracted from the cell pellet using Qiagen RNeasy Plus Micro kit in a final volume of 25 ml. We typically obtained 200-300 ng RNA from neurons, 10-20 ng from microglia and 2-5 ng from astrocytes. RNA quality was quantified using the Bioanalyzer.
Reverse transcription and pre-amplification. RNA (typically, 7.5 ml in a 20 ml total reaction volume) from sorted cells was reverse transcribed using the High-Capacity cDNA Reverse Transcription kit (4368814, Applied Biosystems). Standard Taqman FAM-MGB probes were ordered from Applied Biosystems for all qPCR assays, except our assays to analyse LPS-induced alternative splicing were designed and synthesized at Genentech using FAM-BHQ1 for probe detection. For pre-amplification, up to 100 qPCR assays (primer/probe sets in 20 Â stock concentration) were pooled and diluted to a 0.2 Â concentration. For splicing analyses, separate assay pools were constructed for constitutive, inclusion and skipping isoforms so that different isoforms from the same gene were never in direct competition for primer/probe interaction during pre-amplification. Five microlitres of pooled assays were combined with 2.5 ml of sample cDNA, 10 ml of 2 Â TaqMan PreAmp Master Mix (1410056, Applied Biosystems) and 2.5 ml of TE buffer, and pre-amplified for 14 cycles using the cycling conditions recommended by Applied Biosystems, in a 96-well plate (N8010560, Applied Biosystems) sealed tightly with MicroAmp Clear Adhesive films (4306311, Applied Biosystems). No cDNA pre-amplification was performed for qPCR reactions derived from whole tissue RNA samples.
Fluidigm qPCR analyses. Fluidigm reactions were performed using the 96 Â 96 or 48 Â 48 chips and included 2-3 technical replicates for each combination of sample and assay, except splicing analyses had 6 technical replicates. For sample mixtures, 2.5 ml pre-amplification product from sorted cell samples or 1.5 ml undiluted reverse transcription product from whole tissue samples, was combined with 20 Â GE Sample Loading Reagent (85000735, Fluidigm), 2 Â PCR master mix (58003365-01, Applied Biosystems) and TE buffer in a 10-ml volume, of which 5 ml was loaded into sample wells. For assay mixtures, equal volumes of TaqMan assay (or custom assay for splicing analyses) and 2 Â Assay Loading Reagent (85000736, Fluidigm) were combined, and 5 ml of the resulting mixture was loaded into multiple assay wells for technical replicates. Data from Fluidigm runs were manually checked for reaction quality before analysis, and C t values for each gene target were normalized to C t values for housekeeping genes. Data analysis for gene-level expression was performed using Microsoft Excel software and plotted using GraphPad Prism 6 software. Data analysis for splicing detection assays is described below.
RNA library prep and sequencing. Total RNA extracted from sorted astrocytes, microglial and neuronal cells was subjected to RNA-Seq. The concentration of RNA samples was determined by NanoDrop 8,000 (Thermo Scientific). The integrity of RNA samples was tested using Bioanalyzer RNA 6,000 Pico Kit (Agilent). RNA Integrity Number (RIN) varied between the samples from 3 to 8.5 where RNA from neurons showed lower quality (average RIN B4) compared with RNA from astrocytes (average RIN B5.5) and microglia (average RIN B6.5). cDNA was generated from up to 25 ng of total RNA using Nugen's RNA-Seq method for low-input RNA samples, Ovation RNA-Seq System V2 (NuGEN). (Per manufacturer's instructions, total RNA was neither depleted of rRNA nor polyA-selected.) Generated cDNA was sheared to 150-200 bp size using LE220 focused ultrasonicator (Covaris). Following cDNA shearing, the size of samples was determined using Bioanalyzer DNA 1,000 Kit (Agilent). In addition, sheared cDNA samples were quantified by Qubit dsDNA BR Assay (Life Technologies). One microgram of sheared cDNA was taken into further processing, starting at end repair step, using Illumina's TruSeq RNA Sample Preparation Kit v2 (Illumina). Generated libraries were amplified using six cycles of PCR. Size of the libraries was confirmed using Tapestation (Agilent).
For sequencing of bulk tissue RNA, total RNA was extracted from intact cerebrocortical tissues of PS2APP and non-transgenic mice at 7 and 13 months age. Concentration of RNA samples was determined using NanoDrop 8,000 (Thermo Scientific) and RNA integrity was determined by Fragment Analyzer (Advanced Analytical Technologies). 0.5 mg of total RNA was used as an input material for library preparation using TruSeq RNA Sample Preparation Kit v2 (Illumina). Size of the libraries was confirmed using Fragment Analyzer (Advanced Analytical Technologies).
Library concentrations were determined by qPCR-based method using Library quantification kit (KAPA). The libraries were multiplexed and then sequenced on Illumina HiSeq2500 (Illumina) to generate 30 M of single end 50 base pair reads per library.
RNA-Seq alignment and feature counting. HTSeqGenie 41 was used to perform filtering, alignment and feature counting. HTSeqGenie uses GSNAP 42 to align reads to the genome. We used version mm9 of the mouse genome, and gene models from our internal database, mostly based on RefSeq. The coordinates of these models are available in Supplementary Data 5 (mouse). Only reads with unique genomic alignments were analysed. Such reads whose alignments overlapped any exon of the gene model (even by a single base) were counted towards that gene. (For comparison, we also ran our pipeline with mm10, and the results were found to be similar.) In the LPS study (GSE75246), two astrocyte samples-one from saline and one from LPS-injected-were excluded from the analysis because of evidence of cell type contamination. In the PS2APP study (GSE75431), one astrocyte sample from the PS2APP group was excluded from the analysis for the same reason.
RNA-Seq normalization. nRPKM (normalized Reads Per Kilobase gene model per Million total reads) values were used as a normalized measure of gene expression (except for heat maps, see below). This statistic is an attempt to combine the best of DESeq sizeFactor 43  Note that this is a single normalization based on size factors, with an overall adjustment to approximate total reads, but is not a double normalization. For projects where the size factor is highly correlated with total reads, nRPKM will be very close to the RPKM. For projects where there is a discrepancy, typically due to fluctuations in a few very highly expressed genes, the two statistics will be different. In either case, the nRPKM values for one gene are exactly proportional to the 'normalized count' statistic used by some other authors, which is just number of reads overlapping a gene divided by the size factor. Furthermore, we remind the reader that care should be taken in comparing nRPKM values between projects (that is, between groups of samples which were not normalized together). In particular, our amplification method results in large numbers of intronic reads, so the M value is relatively higher (and nRPKM values relatively lower) than a typical polyA-selected RNA-Seq library.
For cufflinks/cuffdiff users, our 'traditional' RPKM statistic is analogous to the classic-fpkm statistic generated by cufflinks, whereas the nRPKM statistic is analogous to the FPKM statistic provided by the cuffdiff package with geometric normalization (see http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/#librarynormalization-methods for details). For an in-depth comparison of these methods, see Supplementary Data 6. For DESeq2 users, our nRPKM gives a matrix which is the same as the one returned by the fpkm function, multiplied by a constant which is the ratio of total uniquely alignable reads (which we use to calibrate the M) to the sum of all gene counts (which DESeq2 uses to calibrate the M).
RNA-Seq differential gene expression. Differential gene expression was performed with DESeq2 (ref. 43). A pre-filter was applied: only genes with at least ten counts in at least three samples (of either condition) were analysed. P-values for other genes were simply set to 1 and log-fold-changes to 0 for visualization purposes, but such genes were not included in the multiple testing correction.
RNA-Seq splice isoform prediction. Splicing events were predicted and quantified by SGSeq 25 . SGSeq yields 'splicing events' (or just 'events'), which are composed of multiple (but usually just 2) variants. For example, a classical cassette alternative exon would be considered a single event, and the inclusion and skipping isoforms considered the 'splice variants' (or just 'variants'). The variants in turn are composed of 'splice graph features' (or just 'features'). For visualization purposes (for example, see Fig. 7), variants were quantified with the variantFreq statistic, which ranges from 0 (fully skipped) to 1 (fully included). variantFreq was summarized with R's box-plot function: boxes show median and interquartile range, and whiskers show full data. RNA-Seq differential splicing analysis. To identify differential splicing events, first a single count for each variant in each sample was determined by adding together the 5 0 and 3 0 counts for each variant, or, if they represented the same features, by simply taking the single unique value. These counts were then analysed using a variation on the DEXSeq package 26 . Normally (for example, looking at the Bioconductor vignette or manuscript), DEXSeq analyses differential usage of exons across a single gene. In our case, we analysed differential usage of variants across a single event. In fact, within the DEXSeq vignette and manuscript one can find suggestions of using the tool in such a manner, so although not typical, this is not a new idea. As with the differential gene expression, a pre-filter was applied: only variants with at least five counts in at least three samples (of any condition) were analysed. Then, any events for which only a single variant remained were considered to be effectively constitutive and that remaining variant was discarded. For most comparisons these two filters greatly reduced the total number of variants tested, to around 6,000. Then, DEXSeq was run on the remaining data using the SGSeq event as the DEXSeq 'gene' and the SGSeq variant as the DEXSeq 'exon'. The model used was sample þ exon þ condition : exon. After running DEXSeq, to limit the number of tests, the first variant of each event (generally a 'skipping' variant) was discarded. As our data had high levels of intronic reads, possibly due to the Nugen amplification, retained intron events were also discarded.
Fluidigm splicing assay primer design. For each target event, primers were designed using the Primer3 programme 44 (Supplementary Table 1). 'Constitutive' assays targeted constitutively spliced regions within 300 bases of, but not overlapping, the alternative event. For 'Inclusion' assays, one primer was designed to anneal to the inclusion region, the other to the flanking constitutive region, and the detection oligo anywhere between (not necessary overlapping a splice junction). For 'Skipping assays', one primer was designed to anneal to the splice junction, with the default required overlaps (that is, at least 4 bases on the 3 0 end and at least 7 on the 5 0 end). In addition, the following non-default physical parameters were provided to primer 3: qPCR splicing data analysis. Amplification curves were visually inspected and pass/fail calls were adjusted as appropriate. Technical replicates were averaged together, omitting fail calls, to get basic C t values for each assay/sample pair. Assays that failed on all technical replicates for a sample were imputed at C t max þ 1, where C t max was the highest C t (passing) value achieved for that assay in any sample. These imputed values are indicated with different symbols in the plots. À DC t values were generated by subtracting the average C t values for the Actb loading control and taking the inverse, and DDC t values were calculated as ( À DC t inclusion ) À ( À DC t skip ). Imputed À DC t and DDC t values are indicated with different symbols in the plots. Student's t-test was performed for each cell type comparing the DDC t values from vehicle and LPS-treated samples (including imputed values), and P-values were adjusted for multiple testing using the Benjamini-Hochberg method. In Fig. 7, these statistics were summarized with R's boxplot function: boxes show median and interquartile range, and whiskers show full data range (generated using R's boxplot function).
Whole tissue RNA purifications. For comparison of LPS response in whole cortex versus sorted cell populations, fragments of cortical tissue were weighed, minced and added to appropriate volume of lysis buffer RLT, and 350 ml of lysate was used for purification according to the RNeasy Micro Plus kit and protocol (Qiagen). For comparison of PS2APP transgenic versus non-transgenic RNA from intact cortical tissue, hemicortices were rapidly frozen on dry ice and stored until a full collection was obtained. RNA was extracted using the RNeasy Lipid Tissue extraction kit with on-column DNase digestion (Qiagen). Quality and quantity of total RNA samples were determined using ND-1000 spectrophotometer (Thermo Scientific) and Bioanalyzer 2100 (Agilent), respectively.
Microarray RNA processing, hybridizations and scanning. The method for preparation of Cy-dye-labelled cRNA and array hybridization was provided by Agilent. Briefly, 1 mg total RNA was converted to double-stranded cDNA and then to Cy5-labelled cRNA using Quick Amp Labeling Kit (Agilent). The labelled cRNA was purified using RNeasy mini kit (Qiagen). cRNA yield and Cy5 incorporation were determined using ND-1000 spectrophotometer. 750 ng of labelled cRNA was fragmented and hybridized to Agilent's Whole Mouse Genome 4 Â 44Kv2 arrays as described in the manufacturer's hybridization kit, against an equal amount of Cy3-labelled universal mouse reference (Stratagene). Following hybridization, microarrays were washed, dried and scanned on Agilent's G2505C scanner. Agilent's Feature Extraction 11.5 software was used to analyse acquired array images.
Microarray analysis. For Affymetrix data, CEL files were downloaded from Gene Expression Omnibus (GEO) and normalized using the rma algorithm 45 . Probesetgene mappings were derived from bioconductor. For Agilent data, TXT were processed with a custom R script: ControlType weights were set to 0, spots with background channel more than 50 above the test channel were set to the median background intensity, background correction was performed with limma::backgroundCorrect() using the 'normexp' method and an offset of 50, limma::normalizeWithinArrays() with the 'loess' method, and finally limma::normalizeBetweenArrays() with the 'Aquantile' method (all from the limma Bioconductor package). Finally, gene expression values were quantified as 'expression ratio', the ratio of the normalized signals in the test and reference channels. For both platforms, the probeset (Affymetrix) or probe (Agilent) with highest interquartile range (IQR) across all samples (without regard to sample meta-data) was selected for further analysis, and used as the 'official' level of the gene in that data set 46 . For all microarray data, differential expression analysis was performed with the limma package 47 .
Heat maps. For Fig. 5, the count data were first transformed using the 'regularized log transformation' (rlogTransformation function from DESeq2). After this, the 2.5% most variable genes were selected (by standard deviation, without regard to sample meta-data) and used for the heat-map. Hierarchical clustering was performed using the hclust function, with the 'ward.D' clustering method and Euclidean distance function. For the heat maps in Figs 1, 2 and 9, nRPKM values were simply log2-transformed, 'floored' at À 4 (that is, any log2(nRPKM) value below À 4 was set to À 4), then Z-score transformed within each data set. In Fig. 9, all Z-scores greater than 4 were coloured at 4.00 for visualization purposes. For microarrays, log2-scale probeset or probe intensities were Z-score transformed. The selected genes in each plot were those differentially expressed at adjusted P-valuer0.05 and log2-fold change Z2, unless otherwise noted. As necessary, orthologues 48 were displayed in subsequent data sets.
Triangle plots. The triangle plots in Figs 1, 2 and 9 were created starting with the three-dimensional vector of the mean nRPKM of each gene in each of the three cell types. Then these means were normalized to 1, so that they all lay within the plane x þ y þ z ¼ 1. The representation shows their position on the triangle, which is the intersection of this plane with the first octant. This is achieved by multiplication with the matrix [ [ À 1 1 0] [ 0 0 sqrt (3)] ].