Loss of Trem2 in microglia leads to widespread disruption of cell coexpression networks in mouse brain

Rare heterozygous coding variants in the triggering receptor expressed in myeloid cells 2 (TREM2) gene, conferring increased risk of developing late-onset Alzheimer's disease, have been identified. We examined the transcriptional consequences of the loss of Trem2 in mouse brain to better understand its role in disease using differential expression and coexpression network analysis of Trem2 knockout and wild-type mice. We generated RNA-Seq data from cortex and hippocampus sampled at 4 and 8 months. Using brain cell-type markers and ontology enrichment, we found subnetworks with cell type and/or functional identity. We primarily discovered changes in an endothelial gene-enriched subnetwork at 4 months, including a shift toward a more central role for the amyloid precursor protein gene, coupled with widespread disruption of other cell-type subnetworks, including a subnetwork with neuronal identity. We reveal an unexpected potential role of Trem2 in the homeostasis of endothelial cells that goes beyond its known functions as a microglial receptor and signaling hub, suggesting an underlying link between immune response and vascular disease in dementia.

Loss of Trem2 in microglia leads to widespread disruption of cell coexpression networks in mouse brain

Introduction
Genome-wide association studies and genome sequencing have identified more than 25 Alzheimer's disease (AD) risk loci, including common, low-risk variants and rare moderate risk variants, in addition to the classical risk variants in the apolipoprotein E (APOE) gene (Guerreiro et al., 2013b;Hollingworth et al., 2011;Jonsson et al., 2013;Lambert et al., 2013;Naj et al., 2011;Pimenova et al., 2018;Steinberg et al., 2015).Although the identity of many of the associated disease genes and the mechanisms by which they increase risk remain unclear, there is evidence that they cluster around the immune system, protein and lipid metabolism, especially inflammatory response, endocytosis, and amyloid precursor protein (App) metabolism (Hardy et al., 2014;Naj and Schellenberg, 2017;Villegas-Llerena et al., 2016).Many implicated genes encode proteins that are highly expressed in microglia (ABI3, PLCG3, TREM2, SPI1, BIN1, CD33, INPP5D, MS4A6A) and/or have a role in the innate immune system in the brain (Huang et al., 2017;Pimenova et al., 2018).One of these, TREM2 (triggering receptor expressed on myeloid cells 2), has an AD-associated risk allele (R47H) with an effect size in AD similar to that of ApoE4 (odds ratio 2.90e5.05)(Guerreiro et al., 2013b;Jonsson et al., 2013).
TREM2 function can be compromised as a result of rare nonsynonymous variants that cause Nasu-Hakola disease when both alleles are affected (Numasawa et al., 2011;Ridha et al., 2004;Soragna, 2003) or significantly increase the risk of developing AD (Guerreiro et al., 2013b;Jonsson et al., 2013;Sims et al., 2017), behavioral variant frontotemporal dementia (Guerreiro et al., 2013a;Le Ber et al., 2013), semantic variant of primary progressive aphasia or occasionally Parkinson's disease (PD) (Borroni et al., 2014;Liu et al., 2016;Rayaprolu et al., 2013) when 1 allele is affected.Evidence suggests TREM2 may be important for normal brain remodeling, which peaks around adolescence (Chertoff et al., 2013) and in later life in response to age-related damage or pathologies.It is not clear whether the same signals or brain regions are affected by TREM2 activity throughout life.In people where TREM2 is compromised by complete loss of function, symptoms begin in adolescence and tend to implicate frontal lobe dysfunction, whereas in people with preservation of some normal TREM2 activity, symptoms appear much later in life and implicate the hippocampus (Giraldo et al., 2013).Although in both cases, progressive white matter changes in the brain and dementia symptoms occur.
TREM2 is a receptor highly expressed on macrophages, including microglia in the brain (Butovsky et al., 2014;Hickman et al., 2013;Paloneva et al., 2002;Schmid et al., 2002).It appears to be an important damage sensing receptor.It can respond to lipid and lipoprotein species such as phosphatidyl serine (Cannon et al., 2012;Wang et al., 2015), clusterin, and apolipoproteins, including APOE (Atagi et al., 2015;Poliani et al., 2015;Yeh et al., 2016) and nucleotides and anionic species such as heparin sulphate, proteoglycans, or other negatively charged carbohydrates (Daws et al., 2003;Kawabori et al., 2015;Kober et al., 2016), and is required for efficient bacterial clearance (N'Diaye et al., 2009).TREM2 signaling is propagated through the adaptor protein DAP12, which activates a number of pathways including Syk, P13 K, and mitogen-activated protein kinase, which culminate in increased phagocytosis and expression of an antiinflammatory phenotype in microglia (Kleinberger et al., 2014;Neumann and Takahashi, 2007;Peng et al., 2010;Poliani et al., 2015;Takahashi et al., 2007).Loss of Trem2 function culminates in a decrease in the number and activation of microglia in mouse models of AD or in mice treated with cuprizone to damage myelin (Cantoni et al., 2015;Ulrich et al., 2014;Wang et al., 2015).Trem2-deficient dendritic cells secrete more TNF-a, IL-6, and IL-12 compared with wild-type (WT) cells, particularly when activated with lipopolysaccharides suggesting there may be a shift toward cells expressing a proinflammatory and potentially more damaging phenotype in the absence of Trem2 (Hamerman et al., 2006;Turnbull et al., 2006).However, not all findings are consistent.A recent report demonstrated reduced AD pathology in an amyloid mouse crossed with a Trem2 knockout (KO) mouse (Jay et al., 2015).
The main AD-associated TREM2 variant R47H has been shown to alter glycosylation and trafficking of the TREM2 protein between the golgi and endoplasmic reticulum resulting in fewer functional TREM2 receptors in the cell membrane and thus loss of TREM2 function (Park et al., 2015).The presence of this variant also reduces the cleavage of full-length TREM2 to a soluble extracellular fragment and, in both TREM2 risk variant carriers and in people with AD, less soluble TREM2 is present in cerebrospinal fluid (Kleinberger et al., 2014) suggesting TREM2 dysfunction may be a common feature in AD and not just in those AD patients carrying a loss of function variant.Notwithstanding this, recent findings suggest that soluble form of the innate immune receptor TREM2 levels may be reduced or increased depending on the stage of AD and variant (Brendel et al., 2017).
TREM2 and other late-onset AD susceptibility genes MS4A4A/ 4 E/6A, CD33, HLA-DRB5/DRB1, and INPP5D are all part of a distinctive brain coexpression module, which also contains the signaling partner for TREM2, TYROBP, or DAP12 (Forabosco et al., 2013;Hawrylycz et al., 2012;Zhang et al., 2013).This module appears to represent a biological network active in microglial cells with an innate immune function.It is significantly perturbed in AD brain (Forabosco et al., 2013;Hawrylycz et al., 2012;Zhang et al., 2013) and remarkably contains fewer than 150 genes.This module shares identity with peripheral macrophages (Forabosco et al., 2013), and many of the genes in the module are also altered in AD blood cells (Lunnon et al., 2012).Elevated expression of many of these genes, particularly TREM2 appears to be associated with the emergence of amyloid rather than Tau pathology in AD mouse models (Jiang et al., 2014;Matarin et al., 2015).Both TREM2 and Tyrobp have also been identified as major hubs in human APOE-expressing mice following traumatic brain injury (Castranio et al., 2017).
Prominent voices in the field of AD research are proposing that to fully understand the etiology of AD, we have to go beyond reductionist approaches and the amyloid cascade linearity and that there is a need for studies that address the complex cellular context of the disease, which involves interactions between different cell types as the disease progresses across time and tissue (De Strooper and Karran, 2016).It is suggested that temporal resolution can be obtained from cohorts of mice at different stages of the disease.Previous studies using a systems level approach on TREM2 have lacked the aforementioned temporal resolution.Furthermore, the use of microarrays to measure gene expression could have meant that subtle effects are not detected.In our study, we tried to address these concerns by profiling Trem2 KO mice gene expression using RNA-Seq at 2 time points and tissues.We used brain cell markers to infer cell type specificity and detected gene coexpresion dysruption affecting a module with endothelial identity at an early stage that causes widespread disruption of other cell-type subnetworks, including a subnetwork with neuronal identity.

Design
Brain tissue samples were obtained from male Trem2 KO and WT control mice at 2 time points: 4 months and 8 months.These time points span the onset and late disease stages in well-established AD mouse models (Matarin et al., 2015).Hippocampus and cortex were selected because they represent tissues affected in AD at early and late stages, respectively (Mastrangelo and Bowers, 2008;Matarin et al., 2015).RNA-Seq was used to profile the transcriptomes for each sample.Two technical replicates were obtained for each sample.Expression data analyzed in this study are available at Gene Expression Omnibus data repository from the National Center for Biotechnology through accession number GSE104381.
Differential expression (DE) analysis allowed us to detect changes in expression between time points and tissues.Coexpression analysis was performed to detect higher-level disturbances in gene expression networks.Enrichment analysis of the results allowed us to detect functions and pathways more altered in the absence of Trem2.Finally, the integration of cell type markers enabled us to go further and not only detect time and tissue-specific changes but also uncover how the interactions between different cell types were affected and at which time point and tissue these changes were occurring.Fig. 1A displays an overview of the experimental and analytical workflow.

Generation of Trem2-/-mice
The Trem2-/-mouse model (Trem2tm1(KOMP)Vlcg) was generated by knocking a LacZ reporter cassette into the endogenous Trem2 locus in place of exons 2 and 3 and most of exon 4, resulting in a loss of Trem2 function and expression of the LacZ reporter under the control of the Trem2 promoter, as described previously (Jay et al., 2015).The mouse line was originally generated by the trans-NIH KnockOut Mouse Project (KOMP).Frozen sperms were obtained from the UC Davis KOMP repository, and a colony of mice was established at Taconic, Cambridge City, MA, USA.Mice were maintained on a B6 background and shipped to Eli Lilly and Co, Indianapolis, IN, USA, for subsequent tissue collection.

Sample description and RNA-Seq Library construction
A total of 64 samples were used in the study.The hippocampus and cortex were dissected from the left hemisphere of each mouse brain and snap frozen in liquid nitrogen.Tissues were homogenized in Qiazol reagent and total RNA isolated using the RNeasy Plus Universal Mini Kit according to the manufacturer's protocol (Qiagen).RNA integrity was measured for all samples using the Bioanalyser Agilent 2100 series using RNA nano chip.All sequencing libraries analyzed were generated from RNA samples measuring an RNA integrity score 8.9.The Illumina TruSeq Fig. 1.Experimental design, analysis workflow and KO assessment.(A) Experimental design and analysis workflow.Transcriptome profiling of 8 WT and 8 Trem2 KO mice using RNA-Seq.Read alignment to quantify transcripts was done using Kallisto.Different expression analysis was done using Sleuth.Coexpression network analysis was done using WGCNA.Enrichment analysis was done with Enrichr.To estimate cell type associations with coexpression modules, gene markers were obtained from an external source (Zhang et al., 2014).All the results were contrasted to infer time and tissue-specific changes caused by the absence of Trem2.(B) Trem2 and Treml1 isoforms expression levels.Trem2 ensembl transcripts consist of 2 coding transcripts, ENSMUST00000024791 and ENSMUST00000113237, and 2 noncoding, ENSMUST00000132340 and ENSMUST00000148545.Treml1 only has one annotated coding transcript ENSMUST00000024792.Abbreviations: KO, knockout; Trem2, triggering receptor expressed in myeloid cells 2; WT, wild type.
mRNA stranded protocol was used to obtain poly-A mRNA from all samples.Hundred nanogram of isolated mRNA was used to construct RNA-seq libraries.Libraries were quantified and normalized using the NEBNext mRNA Library prep kit for Illumina and sequenced as paired-end 76 bp reads on the Illumina NextSeq 500 platform using technical replicates.

RNA-Seq analysis
RNAseq reads were aligned to the mouse GRCm38 rel 79 reference transcriptome using Kallisto version 0.42.4 (Bray et al., 2016).The reference transcriptome was downloaded from the Kallisto project website ("http://bio.math.berkeley.edu/kallisto/transcriptomes/" n.d.).A batch of quality control methods implemented in the R sleuth package were applied, including principal component analysis, MA plots, qq-plots, and varplot (bootstrapping) (See Supplementary File 1).Quality control revealed consistent quality, and no samples were removed.Differential gene expression analysis was performed on transcripts per kilobase million-normalized counts with the R sleuth package version 0.28.1 (Pimentel et al, 2016;Pimentel et al., 2017).Multiple biological replicates were used for all comparative analysis and a technical replicate for each of the samples.Significant genes, q-value 0.01, were used for gene ontology and pathway analysis.Gene ontology and pathway analysis were performed using the Enrichr database (Chen et al., 2013;Kuleshov et al., 2016;Newhouse, 2017) and Process Network Tool in (Thomson Reuters).Plots and graphics were generated using the sleuth (Pimentel et al., 2016;Pimentel et al., 2017) and ggplot2 version 2.2.1 (Wickham, 2009) R packages.R 3.3.0version was used.
The number of reads per replicate ranged from 16 to 40 million reads (see Supplementary Table 1).The efficiency of the mapping varied from 75% to 87% with the number of reads mapped ranging from 13 million to 30 million.Each sample was profiled for 44,131 transcripts with 16,455 genes detected after quality control filtering.Principal component analysis revealed that samples clustered by age and tissue (see Supplementary File 1).Quality control revealed consistent quality, and therefore, all samples were used in subsequent analyses.

Gene network construction and module detection
For the coexpression analysis, we included all genes that were expressed in at least one time point, condition, or tissue in the final set and with at least 10 reads in more than 90% of samples (Langfelder and Horvath, 2008).Using these criteria, Trem2 is retained for analysis, even though it is not expressed in the Trem2 KO mice.This left us with 14,920 genes (see Fig. 1A).
We used WGCNA version 1.51 to identify modules of coexpressed genes within gene expression networks (Langfelder and Horvath, 2008).To construct the network, biweight midcorrelations were calculated for all possible genes pairs.Values were entered in a matrix, and the data were transformed so that the matrix followed an approximate scale-free topology (see Supplementray File 2 for detailed information).A dynamic tree cut algorithm was used to detect network modules (Langfelder et al., 2008).We used signed networks to preserve the direction of correlations.We ran singular value decomposition on each module's expression matrix and used the resulting module eigengene, which is equivalent to the first principal component (Langfelder and Horvath, 2007), to represent the overall expression profiles of the modules.
Cytoscape v3.5.1 was used to generate network maps of the most highly correlated genes and modules (Shannon et al., 2003).

Module matching between networks
Modules from the Trem2 KO network were matched to those of the WT network using a hypergeometric test to identify the WT module that has the most significant gene overlap and its color was then assigned to the KO module.

Module preservation statistics
We applied both the module preservation Z summary and the medianRank statistics to assess the module preservation from different expression datasets (Langfelder et al., 2011).Unlike the cross-tabulation test, Z summary not only considers the overlap in module membership (MM) but also the density and connectivity patterns of modules.In addition, for our study, network-based preservation statistics only require MM that is identified in the original datasets, reducing the variation coming from various parameters setting to identify new modules in validation datasets.We converted the transcript-level measurements into gene-level measurements using the collapseRows R function (Miller et al., 2011).We adapted the method used to select the most representative probe to select the most representative transcript.The transcript isoform within a gene that had the highest average expression was used to represent that gene as it has been described that it leads to best between-study consistency.Overall, 14,920 genes were retained in our preservation calculation.

Module membership statistics
MM measures correlations between each gene and each module eigengene, effectively measuring how strong is the connection of a gene with the module it has been assigned to and to the other modules too.We compared the correlation of the MM values between the WT and KO networks using 2 approaches using the methods described by the WGCNA authors (Langfelder and Horvath, 2008).The first one uses all the genes to assess overall conservation of the modules.The second one only includes the genes assigned to the modules to assess "hub" conservation (i.e., are the genes more strongly correlated to a module in one of the networks also strongly correlated to the module in the other network).

DE tests for modules
We used student t-tests to statistically compare expression differences between modules as suggested by the WGCNA authors (Langfelder and Horvath, 2008).

Integration of cell type markers
The top 500 most highly expressed genes for different brain cell types (neurons, endothelial cells, astrocytes, microglia, oligodendrocyte precursor cells, newly formed oligodendrocytes, and myelinating oligodendrocytes) were selected from an RNA-Sequencing transcriptome database (Zhang et al., 2014).Of these, 411 astrocytes, 386 endothelial cells, 316 microglia, 407 myelinating oligodendrocytes, 394 newly formed oligodendrocytes, 390 oligodendrocyte precursor cells, and 430 neurons marker genes could be mapped to our dataset and were used in subsequent analyses.The number of markers per cell type was divided by the total number of markers present for that given cell type to take into account for any bias in total number of cell markers available for different cell types.The obtained values were multiplied by 100 to avoid working with very small numbers.

Ethics statement
All animal procedures and experiments were performed in accordance with the Institutional Animal Care and Use Guidelines for Eli Lilly and Company.Protocol number 14e067.

Analysis workflow and RNA-Seq assessment
In this study, we have applied an integrative approach for a systems level characterization of how the absence of Trem2 affects transcriptional expression (see the overview of the experimental and analytical workflow in Fig. 1A).This approach aims to make our results comparable with previous models of AD progression across tissue and time (Matarin et al., 2015).For this purpose, we used hippocampus and cortical regions to model progression across tissue and 2 time points, 4 and 8 months, to model progression with age.We compared the transcriptomes of WT and KO mice identifying differentially expressed genes, and we used coexpression analysis to reveal changes in patterns of expression.To understand the functional implications of those changes, we integrated term enrichment analysis and brain cell-type markers obtained from an external source (Zhang et al., 2014).

Trem2 KO mouse assessment shows that the Trem2 gene is effectively silenced and the Treml1 gene is upregulated
To assess the efficiency of the Trem2 KO process, the levels of expression of the Ensembl (Ruffier et al., 2017) annotated Trem2 transcripts were explored (Fig. 1B).Both protein-coding transcripts and a small processed transcript were silenced in the KO.A noncoding transcript containing a retained intron had higher expression in the KO, particularly at 4m in hippocampus.Treml1, the gene immediately downstream of Trem2, was silent in the WT but upregulated in the KO.The difference in Treml1 expression between WT and KO is the largest effect seen, even more so than that of Trem2, at all time points and tissues (Fig. 1B).It was recently described that this could be an artifact caused by the floxed neomycin selection cassette (Kang et al., 2018) employed in the generation of Trem2 KO mice used for this study.It is yet to determine if this artifact is functionally meaningful but we will interpret our results with caution because it is not possible to unmistakably determine if the effects seen in the KO are due to the absence of Trem2, the upregulation of Treml1, or the combined effect of both.

The effect on gene expression is temporal with a substantially bigger number of DE genes at 4m
Initially, we explored gene expression differences across 4 and 8 months aged mice (see Fig. 1A, Supplementary Files 5 and 6 and Supplementary Tables 2e5).A full list of differentially expressed genes is available in Supplementray Tables 2e5.We found that the number of DE genes between the WT and KO mice was 16-fold higher at 4 months cortex, where a total of 1274 genes were found to be differentially expressed, compared to 8 months cortex, where the expression of 79 genes was altered.The same trend was seen in the hippocampus, although fewer genes were differentially expressed at the early time point 495 at 4 months compared with time point 126 at 8 months.
The only substantial overlap in DE genes across time and tissue was found at 4 months, where the hippocampus and cortex shared 277 DE genes (see Supplementary File 6).This represents more than half of the DE genes in the hippocampus but less than a quarter of the DE genes in the cortex.Among these are genes previously associated with AD.For example, Rbm3, which has been shown to mediate structural plasticity and have protective effects on cooling in neurodegeneration (Peretti et al., 2015), is downregulated in the KO.Nfya, which is also downregulated, has been described to have a role in angiogenesis (Xu et al., 2016).Treml2, which was upregulated in the cortex at 4 months and in both tissues at 8 months, has been reported to contain variants with a protective effect in AD (Benitez et al., 2014).
Besides the overlap at 4 months, there is a reduced number of genes that are differentially expressed in more than one time point and tissue (see Supplementary Tables 6 and 7).Only Trem2 (as expected), Treml1, Grp17, and Fam107a are differentially expressed at all time points and tissues.Interestingly, Grp17 has been previously associated with Nasu-Hakola disease (Satoh et al., 2017), a disease characterized by homozygous loss of function of Trem2 or its signaling partner Tyrobp (Paloneva et al., 2002), and has been successfully targeted with a marketed antiasthmatic drug that reduces neuroinflammation and elevates hippocampal neurogenesis resulting in improvements to learning and memory in older animals (Marschallinger et al., 2015).Fam107a, also known as Drr1, has been described as a stress-induced actinebundling factor that modulates synaptic efficacy and cognition (Schmidt et al., 2011) and has also been connected to brain cancer (Dudley et al., 2014;Le et al., 2010).
3.4.Gene ontology enrichment analysis of DE genes shows timeand tissue-specific enrichment of processes related to AD Ontology analysis showed enrichment for genes involved in response to unfolded proteins and ubiquitination-related functions at 4 months in both tissues (Table 1; also see S8 Table ).Both functions were driven mainly, but not exclusively, by the differential expression of genes of the DnaJ family (DnaJA2 and DnajB2) and heat shock proteins, for example, HspA8 and Hsph1.Enrichment for "response to heat" was significant in the hippocampus and was driven by Camk2g and Camk2b, members of the calmodulindependent protein kinase subfamily involved in calcium signaling.A similar trend was seen in the cortex.
At 8 months, both tissues showed enrichment for membrane targeting-related functions, driven mainly by ribosomal proteins such as Rpl10, but there were also differences between them.In hippocampus, there was enrichment for "regulation of macroautophagy" driven by Uchl1, which has been previously implicated in AD, and also for "organelle transport along microtubule" driven by Cdc42 and CopG1.In the cortex, enrichment for "cellular response to vascular endothelial growth factor stimulus" was driven by Flt1 (also know as VEGFR1) and Dll4, which have been both described as having a role in endothelial sprouting (Eilken et al., 2017;Pitulescu et al., 2017).

Trem2 WT and KO gene coexpression networks analysis identifies modules with cell type and/or functional identities
We generated 2 multitissue and multiage networks, one for WT and one for the KO, and matched the latter to the former (see Methods and Fig. 2).The WT network consisted of 8 modules, each containing between 127 and 6477 gene members.The KO coexpression network comprised 9 modules ranging from 204 to 5213 genes.
Gene ontology and pathway analysis were used to assign biological significance to each module (see S9 Table).In the majority of cases, WT and KO modules were enriched for similar terms.In both WT and KO, the blue module was enriched for functions related to synaptic transmission and neurotransmitter secretion.Similarly, the brown module was enriched for mRNA transport and processing functions in both networks.The tan module was enriched for functions related to translation, protein targeting to membrane, and mRNA catabolism.The turquoise module was strongly associated with ncRNA, mitochondrial-related functions, and house-keeping functions.The red module was associated with RNA splicing and cell cycle-related functions.Conversely, the midnight blue module showed different functions in the WT and KO networks.In the WT, the module was associated with deubiquitination, and in the KO, most of its associated functions were related to synaptic transmission.
We used brain cell-type markers (Zhang et al., 2014) as a proxy to test cell-type enrichment for each module (see Methods and Figs 3, 4).The WT and KO gene expression modules have similar cell type profiles (Supplementary File 3).The results of the enrichment for ontology terms support those of the cell type enrichment analysis, for example, the blue module is enriched for neuron cell markers and for neuronal functions like synaptic transmission.The modules that had a clearer cell type identity were the mentioned blue module with a neuronal identity, the light cyan module with endothelial identity, and the tan module with a microglial identity.

KO-specific coexpression modules emerge notwithstanding broad conservation between WT and Trem2 KO coexpression networks
After identifying the modules in the 2 networks, we wanted to compare how well preserved they were.Given the cell type and functional assignments done previously, we could potentially identify interesting functional or cell-type perturbations between the 2 by looking at the differences between the network modules.Fig. 2A highlights the direct comparison of the WT and KO coexpression networks.It represents the number of genes for each module in each of the coexpression networks and the number of those genes that overlap with the modules in the other network.To formally assess the preservation of the WT modules in the KO network and vice versa, we calculated their Z summary preservation scores and median rank preservation values (Fig. 2B, C) (Langfelder et al., 2011) (see Methods section).
Most WT modules have a KO counterpart, with the exception of the WT pink module.Similarly, all the KO coexpression modules had a WT counterpart, except the KO magenta and yellow modules.
The WT light cyan module was ranked as the least preserved module.The WT tan module had one of the lowest Z summary scores and rankings, as does its KO counterpart.The tan module was of particular interest because it contained Trem2 and its signaling partner Tyrobp and, consistently with Trem2 being expressed in microglia, was assigned a microglial identity using cell type markers.
The KO magenta, midnight blue, yellow, and brown had relatively lower Z summary values but we could not consider any module as not "preserved" because they all had Z summary scores higher than 10 (Horvath, 2011).
Z summary scores were influenced by module size, which can explain the apparent contradiction of the midnight blue module being at the same time the best ranked for preservation and the 1 with the lower Z summary score.This module is markedly different between the WT and KO, partially because the color assignment of the KO midnight blue module is based on 26 overlapping genes (only 5% of the 529 genes in the module) with the WT midnight blue module.
To further study the module preservation between the 2 networks, we used MM values.The modules with the lowest correlations in MM between the WT and KO were the light cyan and tan modules (see Supplementary Files 7 and 8), the same ones that had overall lowest preservation scores.MM effectively measures how strong is the connection of a gene with the module it has been assigned to, calculating the correlation between the gene and the module eigengene.The genes with the highest MM are the most highly connected genes within the module and are referred to as "hubs".This correlation is calculated not only for the module it has been assigned to but also for the other modules too so it is useful to compare different networks where the genes may have been assigned to different modules (more details in the Methods section).Despite broad conservation, there were modules that could not be matched between the 2 networks.The WT pink module had no KO counterpart, but it had a relatively high Z summary score compared with the rest of the modules(see Fig. 2B) and was one of the modules with the highest MM correlation with the KO network (cor 0.89, p-value<0.01;see Supplementary File 8).Interestingly, it was enriched for extracellular matrix organization and exosomerelated functions and both endothelial and microglia cell type markers.Furthermore, the pink module top hub gene Col8a2, encodes a protein whose absence has been described to produce progressive alterations in endothelial cell morphology and cell loss (Jun et al., 2012).Npc2, one of the genes causing frontotemporal dementia when mutated (Zech et al., 2013), is also part of this module.
In the KO network, the magenta and yellow modules were not matched to any WT modules and had low Z summary scores.Seventy-two percent of magenta module genes overlap with the light cyan endothelial module, one of the least preserved modules.Nearly 20% of the genes of the yellow module ovelap with the tan module, the other least preserved module, and 60% of them overlap the blue module with neuronal identity.We considered that this deserved further inspection and we will look at them in more detail in the following sections.Adjusted p-values obtained using a Fisher exact test with Benjamini-Hochberg method for correction for multiple hypotheses using Enrichr (Chen et al., 2013).Key: ATF6, activating transcription factor 6; BP, biological process; GO, gene ontology; SRP, signal recognition particle.3.5.2.The KO yellow module is enriched for proteasome catabolic processes, microglia cell markers, and genes previously implicated in different dementias The KO yellow module contained genes overlapping mostly with the WT blue and WT tan modules (see Fig. 2A).We looked at the different gene ontology and cell type-enrichment profiles and top hub genes of these modules to further understand the yellow module emergence in the KO (see Fig. 3A).
The WT blue module was enriched for neuron cell markers and terms associated with neuronal cell functions, such as neurotransmitter secretion and synaptic transmission.Its top hub gene, Prkar1b, has been previously linked with neurodegeneration (Wong et al., 2014) and the second one, Clstn1, has been shown to be associated with pathogenic mechanisms of AD (Vagnoni et al., 2012).
The WT tan module was enriched for microglia cell-type markers and for translation-related ontology terms and its top hub gene was the ribosomal protein, Rpl18a.
The KO yellow module that emerged as parts of the WT blue and WT tan module had a cell marker profile very similar to the WT tan Fig. 3. Summary of KO emerging modules top hubs and cell type and ontology enrichment.For each module considered, the list of the genes with the strongest correlation within its module, top hubs, is shown together with a list of ontology terms representative of the module function enrichment with its adjusted p-value (see Methods).The cell type identity of the network modules is inferred using gene markers enrichment.Each panel displays the proportion of markers for each cell type present in the module.The counts of the number of markers per cell type was divided by the number of markers present for that given cell type to take into account a possible bias (see Methods for further details).The method allows us to approximate the proportion of each cell type that conform a module.For example, the WT, the blue module is the module with the highest proportion of neuron markers, although it has also a high proportion of other cell-type markers, and the light cyan is enriched for endothelial cell markers.The cell types considered are astrocytes (AS), endothelial cells (ET), microglia (MG), myelinating oligodendrocytes (MO), newly formed oligodendrocytes (NO), oligodendrocyte precursors cells (OP), and neurons (NE).Gene colors represent the module where they have been assigned.Some genes appear as not connected to a module they are assigned to because their correlation to other members of the module is lower than 0.5.Modules are labeled according to their enrichment of cell-type markers.The light cyan module, enriched for endothelial cell markers, is closely connected with the blue module, enriched for neuron cell markers, and in turn connects it to the modules enriched in glia markers in the WT.In the KO, this connection is weakened when the majority of endothelia cell marker genes become members of a newly formed magenta module that is no longer strongly correlated with the blue module.The most correlated gene within each module, the top hub, is shown in an ellipse filled with the color of the module it belongs to.The top hub of the light cyan in the WT is Ttbk1, Tau Tubulin Kinase 1, which regulates phosphorylation of tau, and in the KO its place is taken by App, the gene that codes for the amyloid precursor, both proteins central to the Alzheimer's pathology.The blue module expression association with age is decreased in the KO network but not in the WT (2 tailed t-test with Welch modification, pval ¼ 0.0036) and the midnight blue module expression is increase with age in the KO (2 tailed t-test with Welch modification, pval ¼ 0.017) but not in the WT.(B) Differential expressed genes per module across time and tissue.In each of the 4 panels, the percentage of geneswithin each module that are DE respect the total number of genes in the module is represented for that tissue and time point.The light cyan module is the module with the higher proportion of DE genes at 4 months in both tissues in the WT.The magenta module, which can be considered a light cyan spinoff because the majority of its genes have a counterpart in that module (see Fig. 2A), is in turn the module with the higher proportion of DE genes in the KO.Differences in expression affect different modules in different tissues at 8 months, the tan module being the most affected in the hippocampus while the changes in the cortex are predominantly affecting the magenta module as it is also the case at 4 months.Abbreviations: App, amyloid precursor protein; DE, differential expression; KO, knockout; WT, wild type.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)module (microglia) but was also enriched for proteasome-, ubiquitin-, and catabolism-related processes.
Interestingly, the KO yellow module top hub gene was Atp6v0c.This gene encodes a component of vacuolar ATPase (V-ATPase), a multisubunit enzyme that is necessary for intracellular processes receptor-mediated endocytosis and synaptic vesicle proton gradient generation (Lee et al., 2010).Atp6v0c has been shown to alter autophagy-lysosome pathway function and metabolism of proteins that accumulate in neurodegenerative disease (Mangieri et al., 2014) and has been described as a promising target for therapeutic development (Higashida et al., 2017).Other genes previously associated with different forms of dementia form part of this module.Mapt, the gene that encodes the tau protein, which has been associated with numerous types of neurodegenerative disorders (Spillantini and Goedert, 2013), becomes strongly associated with the module.Pink1, a gene that causes PD when mutated, is also strongly associated, as it is Sqtsm1, a gene involved in sporadic amyotrophic lateral sclerosis pathology (Fecto, 2011).
Altogether, the characteristics of the yellow module relate to known mechanisms involved in the pathology of different types of dementias.

The KO magenta module is enriched for endothelial cell markers and functions
The KO magenta module emerged from the disruption of the WT light cyan module as can be seen from the module overlap in Fig. 2A.Over 72%, 148 of 204, of the genes in the KO magenta module overlapped with those of the WT light cyan module.As we did before for the yellow module, we looked at the different enrichment profiles and top hub genes of these modules to further understand the module emergence in the KO (see Fig. 3B).
The WT light cyan module was enriched for endothelial cell markers and terms associated with endothelial cell functions, like angiogenesis, and also for other functions related to AD such as phagocytosis and response to unfolded protein.Its top hub gene Ttbk1(tau tubulin kinase 1), whose encoded protein is a neuronspecific serine/threonine and tyrosine kinase that regulates phosphorylation of tau, has been implicated in both the pathology of amyotrophic lateral sclerosis (Liachko et al., 2014) and AD (Vázquez-Higuera et al., 2011).
When the light cyan module broke up its KO counterpart became enriched for ubiquitin-and catabolism-related processes.The WT hub gene Ttbk1, lost its hub status in the KO module, where the top hub became App, the amyloid precursor gene.It is also worth mentioning that in the KO, one of the top hub genes, the third one, (see Fig. 3B) is Atp6v0a1, which is one of the components of the vacuolar ATPase together with Atp6v0c, the top hub of the yellow module.
The module was not enriched for endothelial cell-type markers anymore but shifted toward an oligodendrocyte identity.This was caused by the vast majority of endothelial cell markers moving to the new magenta module, which was significantly enriched for genes with blood vessel development function.The magenta module top hub, Ptprb, is also reportedly involved in the establishment of endothelial cell polarity and lumen formation (Hayashi et al., 2013).Interestingly, Flt1 (also known as vascular endothelial growth factor receptor 1) was another top hub in the magenta module.Flt1 has been described to have a role in motor neuron degeneration (Poesen et al., 2008).Most importantly, the disruption of the light cyan module in the KO network suggested that endothelial cells become dyscoordinated in the absence of Trem2, the overexpression of Treml1, or when both effects are combined.

Dyscoordination of the endothelial module weakens neuron module and glia module coexpression
To help us visualize better the differences between the WT and KO networks, we plotted the coexpression networks with genes colored by their module assignment, including only correlations stronger than 0.5 (see Fig. 4A).In this way, we could see how the disruption of the endothelial module in the KO network weakens the connection of the whole network in general and in particular the connections between the neuronal subnetwork and the modules with glial identity.
In the WT, the light cyan module (enriched for endothelial cell functions and markers) was connected on one end with the blue module (enriched for neuron functions and cell markers) and on the other end, the light cyan endothelial module connected with other modules enriched for glia markers and functions.In particular, the WT endothelial subnetwork was connected with the brown module, which is mainly enriched for myelinating oligodendrocytes and astrocytes.In the KO, these connections were weakened by the endothelial module disruption.The KO light cyan module counterpart remained strongly connected to the neuron module and App became its top hub gene.The KO-specific magenta module, which has the biggest proportion of endothelial gene markers lost to the light cyan module in the KO, lost the strong connection to the neuron module but remained correlated with the brown module.
Interestingly, the new KO yellow module, enriched for proteasome functions and genes previously associated with different types of dementias, was also strongly correlated to the neuron module.
We can also see in Fig. 4A how the midnight blue module is markedly different between the WT and KO and is strongly associated with the blue neuronal module in the KO.The KO midnight blue module was enriched for neuron markers but also for the gamma-aminobutyric acid signaling pathway (see S9 Table ).While the blue module expression decreases with age in the KO (t-test with Welch correction, p ¼ 0.004), the expression of the midnight blue increases (t-test with Welch correction, p ¼ 0.02) (see Supplementary File 4).Interestingly, gamma-aminobutyric acid signaling has been described to be increased in 5xFAD Alzheimer's mouse models and its reduction has been described to rescue the impairment of long-term potentiation and memory deficit (Wu et al., 2014).
Mapping the differentially expressed genes to the modules obtained in the previous section, we were able to visualize which modules are affected at different ages in different tissues (see Fig. 4B).At 4 months, most of the changes in expression in both hippocampus and cortex affected the light cyan module in the WT or the magenta module, product of the WT light cyan disruption, if we look at the KO network.This points to the endothelial module disruption being an early event.The light cyan module is enriched for response to unfolded protein gene ontology terms in the WT and for ubiquitination in the KO.This is consistent with what we found in the enrichment analysis for genes being differentially expressed at 4 months (described in Table 1).At 4 months, the rest of the modules were broadly affected to a similar extent in both tissues.Notable exceptions are the WT-specific pink module affected specifically in the cortex at 4 months and the yellow KO module that shows larger differences at 4 months in the cortex, although not as dramatically as the pink module.
At 8 months, there was a very different picture with changes in expression affecting different modules in different tissues.The tan module is the more affected in the hippocampus.Again, consistent with the enrichment analysis performed previously for the differentially expressed genes, we find enrichment for macroautophagy in both the tan module and list of differential expressed genes at 8m in the hippocampus.The magenta endothelial module is the most affected in the cortex.This again agrees with the enrichment for "cellular response to vascular endothelial growth factor stimulus" specifically in the list of differentially expressed genes at 8 months cortex.
Overall, we found that the absence of Trem2, or potentially the overexpression of Treml1 or the combined effect of both, in microglia produces a ripple effect that provokes the disruption of a subnetwork of genes enriched for endothelial cell gene markers and functions at an early time point, aged 4 months, with the effect being stronger in the cortex at this early time point.At a later stage, aged 8 months, cortex and hippocampus are affected differently with the subnetwork of genes with microglial identity being more perturbed in the hippocampus.These are time-and tissue-specific effects that evolve with age.

Discussion
Many risk genes contributing to dementia appear to have pleiotropic effects.Rare loss of function variants in the TREM2 gene can cause Nasu-Hakola disease when both alleles are affected (Paloneva et al., 2000) or increase the risk to develop AD, PD, or frontotemporal dementia in heterozygous carriers (Guerreiro et al., 2013b;Jonsson et al., 2013).Further study of the complex molecular consequences of the absence of TREM2 could have a big impact on our understanding of AD pathogenesis from a systems level view (De Strooper and Karran, 2016).We generated RNA-Seq data to profile the hippocampus and cortex at aged 4 and 8 months to explore the temporal and spatial consequences of Trem2 absence.Differentially expressed genes were enriched for functions related to AD and other forms of dementia.Furthermore, WGCNA analysis coupled with enrichment for cell-type markers highlighted the disruption of a coexpression module with a strong endothelial identity.The absence of Trem2 in microglia appears to generate a ripple effect that causes the disruption of an alignment between endothelial cells and microglia affecting the coordination of neurons and glia.Interestingly, this suggests a new role for Trem2 beyond its known functions as a microglial receptor and signaling hub, linking the immune response with vascular consequences as potential underlying causes of AD.
Our study involved expression profiling of 8 Trem2 KO and 8 WT mice using RNA-Seq.We sampled the hippocampus and cortex at 4 months and 8 months of age.This allowed us to model the disease progression across tissue and time and enables results to be compared with previous expression profiling of transgenic mice that develop amyloid or tau pathology (Matarin et al., 2015).Surprisingly, the level of expression of Trem2 in the WT mice increases from 4 months to 8 months in both tissues, at least for the most highly expressed isoform.The differences may be attributed to the use of RNA-Seq, which let us consider each isoform independently.In fact, of the 2 protein-coding isoforms included in Ensembl, we only detect significant differential expression between 4 months and 8 months for Trem2-201.As previously noted, microglia exclusively express TREM2 in the brain (Ulrich et al., 2017), and the increased levels of expression may reflect an increase in microglia density with age (Conde and Streit, 2006).Other differential effects between 4 and 8 months were also found.Nearly 1300 genes had significant differential expression at 4 months, which was reduced to only 79 genes by 8 months.Interestingly, the number of DE genes was greater in the hippocampus than in the cortex at 8 months suggesting a shift in dominance in Trem2 activity between brain regions.These alterations follow a regional distribution that is consistent with that of Nasu-Hakola, whose pathology affects center in the cortical region earlier in life than AD, which has a much later onset and centers at least initially in the hippocampus (Ridha et al., 2004).
Surprisingly, the Treml1 gene appears to be upregulated across all times and tissues.This may be because Trem2 normally inhibits Treml1 expression or, alternatively, is an artifact of knockdown as has been previously suggested (Kang et al., 2018).Nevertheless, it is one of the few effects that are consistent across ages and tissues.Furthermore, the upregulation of Treml1 is the largest effect we see.Therefore, our results should be interpreted as either the consequences of knockdown of Trem2, increase in Treml1, or both together.This is an important finding that should be taken into account on all studies using the same Trem2 KO mice used in our study.To completely disentangle if the effects, especially in the endothelial or other nonmicroglial cell types, can be interpreted as an effect of Trem2, the distribution of expression and effects of overexpressing Treml1 should be established.There is further justified interest in the study of TREML1 expression as it has been recently described that its upregulation may have a protective effect in AD with high levels of both TREM2 and TREML1 associated with decreased disease risk (Carrasquillo et al., 2016).Interestingly, the authors of the study propose, as we do here, a possible link between neuroinflammation and vascular homeostasis as related mediators of neuronal protection and injury in AD that involves both TREML1 and TREM2.
Our WGCNA coexpression analysis coupled with cell brain celltype markers (Zhang et al., 2014) used as a proxy of cell-type enrichment helped us to identify modules with both functional and cell type identities.We found a module, blue, enriched for synaptic transmission ontology terms had the highest enrichment for neuronal markers.We also identified a module that contained Trem2 and Tyrobp and was mostly enriched for microglia and myelinating oligodendrocytes.The module containing Trem2 and Tyrobp in our study is larger and has only modest gene overlap compared to what was found in previous brain coexpression network analyses (Forabosco et al., 2013;Zhang et al., 2013).Differences are to be expected.Previous studies involved the analysis of postmortem human control and/or AD brain tissue, and in both cases, expression was profiled using microarray data.Nevertheless, some of the key genes such as Fcer1g and Ly86 identified in the AD immune module (Zhang et al., 2013) and Cxcr1 are present in the Tyrobp module in our study or, as is the case of Cxcr1, in the phagocytic yellow module that is specific to the KO network.
This emerging KO yellow module, which draws its genes from the module with neuronal identity and the module containing Trem2 and Tyrobp, was enriched in catabolic, proteasome, and degradation associated functions.Interestingly, its top hub is Atp6v0c, a gene that encodes a component of vacuolar ATPase, a multisubunit enzyme that is necessary for intracellular processes receptor-mediated endocytosis and synaptic vesicle proton gradient generation (Lee et al., 2010), and has been described to alter autophagy-lysosome pathway function and metabolism of proteins that accumulate in neurodegenerative disease (Mangieri et al., 2014) and has been proposed as a possible target gene for therapy (Higashida et al., 2017).Furthermore, a number of genes that have previously associated with dementia form part of this module.Mapt, the gene encoding for the tau protein associated with different types of neurodegenerative disorders (Spillantini and Goedert, 2013), Pink1, a gene that can cause PD when mutated, and Sqtsm1, a gene involved in sporadic amyotrophic lateral sclerosis pathology (Fecto, 2011), are all part of this module suggesting that all these neurodegenerative pathologies may share an underlying molecular mechanism in which Trem2 plays a central role.
The WT light cyan module (enriched for phosphorylation regulation, angiogenesis and protein processes, and junctionrelated cellular compartments) was disrupted in the KO, showing the lowest preservation and breaking up to form the magenta module as a consequence.The KO-specific magenta module was enriched in angiogenesis, similarly to the WT light cyan module with whom it shares most of its genes, but it was not enriched for phosphorylation regulation.The top hub gene in the WT light cyan module was Ttbk1 (tau tubulin kinase 1) but in the disrupted KO module App, the App becomes the top hub status, while the top hub in the magenta module is Ptprb, which plays an important role in blood vessel remodeling and angiogenesis.This disruption suggests that Trem2 provides a link between endothelial cells and neurons that appears essential to the survival of neurons.When the coordination between the 2 falters, App drives a subnetwork of genes involved in catabolic processes that is closely associated with the neuronal subnetwork.As pointed out before, it is also possible that Treml1, or both Trem2 and Treml1 acting together, may have a role in the effects we see.This finding also points to a role of Trem2 in the mechanisms behind the increased risk of dementia caused by vascular damage (Mantzavinos and Alexiou, 2017) and could also help explain the brain blood vessel damage described in cases of Nasu-Hakola disease (Kalimo et al., 1994).
The changes affecting the module with endothelial identity described previously occurred mainly at 4 months in both tissues, but the effect was strongest in the cortex.These findings are consistent with human studies on AD where early changes are described in the prefrontal cortex region, for example, changes in monoamine oxidase A and B (Kennedy et al., 2003) and these correlate with Braak stage (Braak and Braak, 1997).At 8 months, the changes are different from 4 months.A more phagocytosisdominated module is altered in the hippocampus perhaps coinciding with the emergence of AD pathologies such as plaque deposition.Somewhat surprisingly, a number of complement genes are differentially expressed at 4 months in this module, C1qb in the hippocampus and C1qa in the cortex, but not at 8 months, where they lose their previous association with the module having neuron identity.The complement system becomes activated in response to the presence of amyloid and is therefore a pathological hallmark of AD (Lui et al., 2016).
To conclude, we have found that the timing and magnitude of the effects described in the Trem2 KO mice is different for different brain regions.Overall, there were more changes in the cortex in younger mice capturing what is observed in Nasu-Hakola disease, which is characterized by frontotemporal dementia and blood vessel dysfunction (Kondo et al., 2002), while the hippocampus in older mice had more changes mirroring (R. Guerreiro et al., 2013b;Jonsson et al., 2013)the vulnerability of the hippocampus in later onset AD (Mastrangelo and Bowers, 2008;Matarin et al., 2015).WGCNA analysis and brain cell-type marker integration revealed links between several molecular processes previously described to play a part in AD pathology.In our Trem2 KO mice, a central role for amyloid processing emerged, as did changes in phagocytosis and altered vascularity.This may reflect a more complex interplay between an absence of TREM2 and its effects in different tissue and time and underlines the importance of considering the impact of disease and risk genes on the whole brain network when exploring cause, effect, and developing treatments.It will be important to see if these findings can be replicated in people carrying the R47H-risk variant associated with increased risk for AD (Guerreiro et al., 2013b;Jonsson et al., 2013).

Fig. 2 .
Fig.2.Module preservation between WT and Trem2 KO networks.(A) Correspondence of WT set-specific modules andTrem2 KO-matched modules.Each row of the table corresponds to one WT set-specific module (labeled by color as well as text), and each column corresponds to one KO module.Numbers in the table indicate gene counts in the intersection of the corresponding modules.Coloring of the table encodes $-ylog(p)$, with p being the hypergeometric test p-value for the overlap of the 2 modules.The stronger the red color, the more significant the overlap is.The table indicates that most WT set-specific modules have a KO counterpart.(B) Preservation of WT modules in KO network.The left panel shows the composite statistic median rank versus module size.The higher the median rank the less preserved is the module relative to other modules.Since median rank is based on the observed preservation statistics (as opposed to Z statistics or p-values), it is independent on module size.The right panel shows the composite statistic Z summary.If Z summary >10 there is strong evidence that the module is preserved(Langfelder et al., 2011).If Z summary ytextless 2, there is no evidence that the module preserved.Note that Z summary shows a strong dependence on module size.(C) Preservation of KO modules in WT network.The same as (B) but in this case, we assess the preservation of the modules obtained in the KO network to see how preserved they are compared to the ones obtained in the WT network.Abbreviations: KO, knockout; Trem2, triggering receptor expressed in myeloid cells 2; WT, wild type.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig.3.Summary of KO emerging modules top hubs and cell type and ontology enrichment.For each module considered, the list of the genes with the strongest correlation within its module, top hubs, is shown together with a list of ontology terms representative of the module function enrichment with its adjusted p-value (see Methods).The cell type identity of the network modules is inferred using gene markers enrichment.Each panel displays the proportion of markers for each cell type present in the module.The counts of the number of markers per cell type was divided by the number of markers present for that given cell type to take into account a possible bias (see Methods for further details).The method allows us to approximate the proportion of each cell type that conform a module.For example, the WT, the blue module is the module with the highest proportion of neuron markers, although it has also a high proportion of other cell-type markers, and the light cyan is enriched for endothelial cell markers.The cell types considered are astrocytes (AS), endothelial cells (ET), microglia (MG), myelinating oligodendrocytes (MO), newly formed oligodendrocytes (NO), oligodendrocyte precursors cells (OP), and neurons (NE).(A) Emergence of KO yellow module by an increase in expression correlation of parts of the WT blue and tan modules.(B) Emergence of KO magenta module by the disruption of the WT light cyan module.Abbreviations: KO, knockout; WT, wild type.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Table 1
Summary of ontology enrichment for the differentially expressed genes across time points and tissues