Identification of evolutionarily conserved gene networks mediating neurodegenerative dementia.

Identifying the mechanisms through which genetic risk causes dementia is an imperative for new therapeutic development. Here, we apply a multistage, systems biology approach to elucidate the disease mechanisms in frontotemporal dementia. We identify two gene coexpression modules that are preserved in mice harboring mutations in MAPT, GRN and other dementia mutations on diverse genetic backgrounds. We bridge the species divide via integration with proteomic and transcriptomic data from the human brain to identify evolutionarily conserved, disease-relevant networks. We find that overexpression of miR-203, a hub of a putative regulatory microRNA (miRNA) module, recapitulates mRNA coexpression patterns associated with disease state and induces neuronal cell death, establishing this miRNA as a regulator of neurodegeneration. Using a database of drug-mediated gene expression changes, we identify small molecules that can normalize the disease-associated modules and validate this experimentally. Our results highlight the utility of an integrative, cross-species network approach to drug discovery.

Introduction

Frontotemporal dementia (FTD) is a highly heritable common cause of presenile dementia often caused by dominantly-acting mutations in the microtubule associated protein tau (MAPT; tau), and in two other genes, GRN and C9orf72[1]. Tau pathology is also a core feature of numerous other neurodegenerative syndromes, including Alzheimer’s Disease (AD) and Progressive Supranuclear Palsy (PSP), and is associated with the pattern of cognitive decline in AD (PSP; for a review see[2-4]). Thus, identifying the mechanisms by which tau and other dementia-causing mutations lead to neurodegeneration is of foremost importance in developing new therapies for dementia.

Translating mechanistic studies in mouse models of dementia to human clinical trials has proven challenging[5,6]. We reasoned that one contributing factor, aside from the evolutionary distance between mouse and humans[7], is that common laboratory mouse strains are highly inbred. Although genetic background profoundly impacts the biochemical and behavioral repertoires of mouse models of AD[8,9], virtually all studies of neurodegeneration in mice rely on a single inbred strain[10,11].

Here, we take a multi-step, systems biology approach to identify robust, human disease-relevant gene networks. To minimize the likelihood that findings were dependent on a single genetic background, we crossed transgenic mice harboring the FTD-causing P301S mutation onto three distinct backgrounds and identified early transcriptomic changes observed across mutant F1 progeny, comparing vulnerable regions, such as the cortex and brainstem, to regions less affected in the human tauopathies, such as the cerebellum[12]. Through gene network analysis, we identify two major co-expression modules that are conserved across genetic backgrounds and in humans. These modules delineate key pathways deregulated in disease and are differentially enriched for genetic drivers of human FTD, AD and PSP. Using an unbiased integrative genomics approach, we uncover genome-wide relationships between miRNAs and their biological targets, and experimentally validate these relationships by demonstrating that overexpression of a predicted driver, miR-203, recapitulates the mRNA expression patterns observed with disease and promotes neurodegenerative pathways. As a proof of principle, we also show that one can leverage these highly conserved co-expression modules to identify tool compounds that mitigate neurodegeneration induced by miR-203.

Results

We analyzed F1 offspring resulting from crosses of TPR50 mice expressing the human TauP301S transgene[13], onto FVB, DBA and parent C57BL6 lines (Fig. 1a). All three F1 crosses share key features of the disease[13,14], including decreased survival (Fig. 1b) and body weight (Supplementary Fig. 1a–b), progressive accumulation of hyperphosphorylated tau in the cortex and hippocampus (AT8 staining; Fig. 1c), and prominent astrogliosis and microgliosis (Fig. 1d–f). Consistent with previous reports that FVB mice can be more vulnerable to neurodegeneration, we observed that F1 TauP301S transgenics on the FVB/C57 (TPR50xFVB/NJ) background had decreased survival and increased pathological tau as compared to the C57/C57 (TPR50xC57BL/6J) and DBA/C57 (TPR50xDBA/2J) progeny (Fig. 1b) (Supplementary Fig. 1c–h). Importantly, no significant decrease in the weight of any of the affected brain areas was detected in transgenic animals at six months (Supplementary Fig. 1i), consistent with previous data showing no cell loss at this stage[13].

Figure 1:

Experimental Design and Characterization of TPR50 Tau Transgenic Mice in Divergent Genetic Backgrounds.

(a) Schematic of the experimental design highlighting that TPR50 mice were crossed with three genetically divergent mouse strains and that tissue from the cortex, hippocampus, brain stem and cerebellum from the resulting F1 crosses was isolated at three and six months for mRNA- and miRNA-seq and downstream weighted co-expression network analysis (WGCNA). (b) Kaplan-Meier survival curve showing significantly decreased survival of all transgenic mice compared to their wild type littermates, as well as decreased survival of the FVB/C57-Tg mice compared to DBA/C57-Tg and C57/C57-Tg mice (two-sided log rank test, p=0.0025, n=180; 15 males and 15 females/condition). (c) Representative images from three independent experiments of phospho-tau specific AT8 staining and hematoxylin counterstaining in coronal and sagittal brain sections of three- and six-month-old transgenic mice of all three strains (n=4 mice/group; sagittal scale bar = 4μm, coronal scale bar = 3μm). (d) Representative images from three independent experiments showing cortical sections of six-month-old wild type and transgenic mice of all three strains stained against Iba1 (red) and GFAP (green) and the nuclear counterstain DAPI (blue) (scale bar=50μm). (e-f) Quantification of GFAP (e) and Iba1 (f) positive cells from 6-month-old cortical slices (n=6 images/mouse, 3 mice per genotype; unpaired two-sided t-test, error bars= SEM). The center line represents the mean.

We next investigated whether we could identify a robust, disease-specific transcriptional signature by sampling brain regions varying in disease vulnerability, reasoning that key disease-related gene expression patterns would overlap across affected areas (cortex, hippocampus, brainstem), but be absent in cerebellum, which is unaffected[12]. We performed RNA sequencing (RNA-seq) of four brain regions and two-time points across the three genetic backgrounds (36 samples per region, per time point; total of 288 RNA-seq samples; Fig. 1a; Online Methods). The differentially expressed genes from transgenic and wild type mice clearly separate tau mutant from wild type genotypes within affected brain regions (Supplementary Fig. 2a–c), but not in the cerebellum (Supplementary Fig. 2a), and are consistent are across strains (Supplementary Fig. 2d, e). There was significant overlap in differential expression (Online Methods; Supplementary Fig. 3a) between the cortex and hippocampus, and significant correlation between fold changes in differentially expressed genes (Rank-rank hypergeometric test; p<0.05) from these two regions (R2=0.76, slope=0.92, Supplementary Fig. 3a). We find lesser, but significant correlation in gene expression between the cortex and brain stem (R2=0.33, slope=0.82, Supplementary Fig. 3a). Consistent with the absence of significant cell loss at six months in this model[13], positive regulators of apoptosis (GO ID:0043065) do not show significant changes across all 3 strains (Supplementary Fig. 3b–e) and negative regulators of apoptosis (GO ID :0043066) show modest, more significant changes (Supplementary Fig. 3f–i).We note that although we do not observe frank apoptosis, we do see evidence for initial activation of cellular pathways associated with inflammation in cerebral cortex by six months, including increases in GFAP and Iba-1 immunoreactivity in cortex (Fig. 1d–f).

We observed no significant overlap in gene expression between the cortex and cerebellum (Supplementary Fig. 3j–k) and no correlation between differentially expressed genes in cortex and cerebellum (R2=0.01, slope=0.14, Supplementary Fig. 3a). Even more striking is that the expression changes of the top upregulated genes in the cerebellum are reversed in the cortex (Supplementary Fig. 3l), implying that protective changes may potentially occur in the cerebellum, as has been previously suggested[12].

Identification of disease-relevant mRNA modules

To place gene expression changes in a systems level framework, we performed consensus Weighted Gene Co-expression Network Analysis (cWGCNA[15,16]) across all three heterozygous F1 mouse strains to identify relationships not observed in only a single background (Supplementary Fig. 4a, Supplementary Table 1). We identified four mRNA modules significantly correlated with the transgenic condition in cortex – which we initially labeled the salmon, turquoise magenta and pink modules per WGCNA conventions (Bonferroni-corrected p < 0.05, Fig. 2a). Three of these modules are downregulated (salmon, magenta and pink), while one module is upregulated in transgenic mice (turquoise) (Fig. 2a). Cell-type enrichment analysis demonstrated that the salmon module is enriched in neuronal markers, the magenta module for oligodendrocyte markers, and the turquoise module for microglial, astrocyte and endothelial markers (Fig. 2b) across multiple different cell type specific gene expression datasets (Supplementary Fig. 4b). Gene ontology (GO) analyses is consistent with the cell type enrichments, showing that the salmon module is enriched for synaptic pathways, magenta for DNA repair and transcription, while the turquoise module is enriched for immune and inflammatory categories (Supplementary Fig. 4c–d).

Figure 2:

mRNA Consensus Co-expression Network Analysis.

(a) Signed association (Pearson correlation) of the mRNA module eigengenes with transgenic condition. Modules with positive values indicate increased expression in transgenic mice; modules with negative values indicate decreased expression in transgenic mice. Red dotted lines indicate Bonferroni corrected p<0.05 for multiple comparisons (n=15 modules, n=36 mice/region) using p-values obtained from two-sided test for Pearson correlation. (b) Cell-type enrichment of modules (average n=200 genes) using mouse genes in mRNA modules (fisher’s two-sided exact test, ***FDR<0.005). (c) Co-expression PPI network plot of the NAS module. Top 30 hub genes are indicated by name. The edges between nodes represent both gene co-expression and PPI permitting us to focus on hub genes observed at both the RNA and protein level. (d) GO term enrichment of the NAS module using 794 NAS modules genes (permutation test Z-score). (e-i) Trajectory of the NAS module eigengene in the cortex (e) and cerebellum (i) (unpaired two-sided Wilcox rank test, n=6 mice/group). (f) Co-expression PPI network plot of the NAI module. Top 30 hub genes are indicated by gene name. (g) GO term enrichment of the NAI module using 1833 NAI module genes. (h-j) Trajectory of the NAI module in the cortex (h) and cerebellum (j) (Unpaired two-sided Wilcox rank test, n=6 mice/group). Boxplots in e, h, i and j, the upper and lower lines represent the 75th and 25th percentiles, respectively. The center line represents the median.

We next calculated the correlation between expression changes and an established early disease marker, the burden of hyperphosphorylated tau, which has been shown to correlate with disease progression in human FTD and AD[17,18]. Genes correlated highly with phospho-tau levels (R2>0.6; Online Methods) were found to be significantly enriched in the salmon (42-fold) and turquoise modules (9.5-fold), but not in the magenta and pink modules (Supplementary Fig. 4e). We therefore focused our subsequent analysis on the salmon and turquoise modules. Given the consistent cell type enrichments and strong gene ontology enrichments (Fig. 2b, d-g; and Supplementary Fig 4b–d) we relabel the salmon module, the Neurodegeneration-Associated Synaptic (NAS) module, and the turquoise module, the Neurodegeneration-Associated Inflammatory (NAI) module, so as to provide more intuitive names for these core disease-associated modules.

To investigate further if the expression changes were dependent on changes in cell type composition, or reflected cell intrinsic signaling changes, we applied a multivariate linear regression model using the top one hundred cell-type specific marker genes for five major cell types[19,20] (Online Methods). Both modules retained their significant association with the transgenic condition (salmon, r= −0.73, p= 8.6e-07; turquoise, r= 0.71, p=3e-06), indicating that they are independent of changes in cell-type composition caused by neuronal loss or gliosis. Moreover, both the NAI and NAS module eigengenes show changes across affected brain regions (cortex, Fig. 2e, h; hippocampus, Supplementary Fig. 4f–g; brain stem, Supplementary Fig. 4h–i) but not cerebellum (Fig. 2i, j), prior to decreases in brain weight that would be indicative of cell loss (Supplementary Fig. 1i). This is further supported by analysis of expression profiles from sorted neurons and glia from the Tg4510 and PS2APP AD model which show cell intrinsic changes in the NAI and NAS module trajectories (Supplementary Fig. 5a–d). The NAI module is highly preserved in both microglia and astrocytes and significantly upregulated in these populations in transgenics (Supplementary Fig. 5a, c, d–l), while the neuronal NAS module is highly preserved in sorted neurons and downregulated in transgenics (Supplementary Fig. 5b, e–l).

We next used experimentally derived databases of human protein-protein interactions (PPI) from Inweb[21] and Biogrid[22] to create an integrated co-expression-PPI network[23] to functionally annotate network edges and provide independent validation[24] (Fig. 2c, f).

One of the hubs of the NAS module, Stx6, a SNARE protein[25], is also one of the few known genetic risk factors for tauopathy, having been previously implicated by GWAS and eQTL analyses in PSP. Other hubs, such as Syt4, Nlgn1, Nrgn also play important roles in synaptic maintenance and function[26,27]. Functional annotation of the NAS module with GO and KEGG terms, confirms its broader relationship to synaptic function (Fig. 2d and Supplementary Fig. 4c). In contrast, functional annotation of the NAI module highlights pathways associated with astrogliosis and inflammatory changes[28], including T cell activation and Toll-like receptor signaling pathways (Fig. 2g and Supplementary Fig. 4d). Among the most central hub genes of the NAI module are Tnip1 and Malt1, important inflammatory signaling genes[29]. The integrated co-expression-PPI network derived from the NAI module also contains the gene Fus, mutations in which account for approximately 4% of familial ALS, and which forms aggregates in Amyotrophic Lateral Sclerosis (ALS) and FTD[30].

Reproducibility of disease-specific modules across mouse models and human disease

Further analysis of 4 independent datasets from mouse models harboring pathological tau mutations including a replication dataset that we generated using the Tg4510 tau model (Methods, Supplementary Fig 5f) confirm the preservation of the NAI and NAS modules (Supplementary Fig 5e), and the expected changes in gene expression early in the disease course (Supplementary Fig. 5f-h, l). Comparison of the NAS and NAI modules in mice harboring 4 different AD and FTD risk mutations (PS2APP, CRND8, APP/PS1 and GRN - see Online Methods) was also performed to assess to what extend these modules represent convergent pathways across models harboring different pathological proteins. The NAS and NAI models are highly preserved and their disease associations are maintained and both the NAI and NAS module expression changes occur with similar temporal profiles across each of these models (Supplementary Fig. 5e, i–l).

Considering this cross-mouse model conservation, we next assessed module preservation in human disease models and post mortem tissue samples to validate their human relevance (See Online Methods, Supplementary Table 2a for list and characteristics of all human datasets, Supplementary Table 2b–c). The NAS and NAI modules are dysregulated in patient cortex, but not in cerebellum, both in the previously published microarray-based data and in our newly generated RNA-seq dataset (Fig. 3a–d and Supplementary Fig. 5e, m–n) consisting of both tau-positive and tau-negative FTD (Fig. 3c). Moreover, the NAS module and its disease-trait relationship are preserved in iPSC-derived neurons from human FTD patients carrying GRN mutations (Supplementary Fig. 5e, o).

Figure 3:

Transcriptomic and Proteomic Analyses in Human FTD samples

(a) Scatterplot showing Pearson correlation of sub-sampled discovery (Control n=7, Tau-positive FTD n=5) and replication FTD (Control n=7, Tau-positive FTD n=5) dataset. P-values obtained from two-sided test for Pearson correlation are shown. (b) Module preservation in human FTD (cortex) using module definitions from strain independent transgenic mouse network analysis. (c-d) NAS and NAI module eigengene expression in human FTD and control samples in the cortex (c, Control n=14, Tau-positive FTD n=10, Tau-negative FTD n=6) and cerebellum (d, Control n=10, Tau-positive FTD n=7). Tau-positive FTD (FTD-tau Pos.) and tau-negative FTD (FTD-Tau Neg.) are shown (unpaired two-sided Wilcox rank test). (e) Log2 fold change of the top 20 NAS and NAI module genes at the mRNA and protein level. (f) NAS and NAI module eigengene (ME) in human FTD and control protein samples from the cortex. Progranulin-positive FTD (FTD GRN Pos.) and progranulin-negative FTD (FTD GRN Neg.) are shown (unpaired two-sided Wilcox rank test). (g-o) NAS and NAI module eigengene (ME) in various neurological diseases – (g) Human AD Temporal Cortex (Control n=52, AD n=52; Allen et al., 2016), (h) Human AD frontal cortex (Control n=308, AD n=157, Zhang et al., 2013), (i) Human AD frontal cortex proteomics (Control n=15, AD n=20; Seyfried et al., 2017), (j) Human Pathological Aging temporal cortex (Control n=70, Pathological Aging n=30; Allen et al., 2016), (k) Human ALS frontal cortex (Control n=9, C9orf ALS n=8, Sporadic ALS n=10; Prudencio et al., 2015), (l) Human PSP temporal cortex (Control n=73, PSP n=83; Allen et al., 2016), (m) Human Major depressive disorder (MDD) (Control n=67, MDD n=66; Chang et al., 2014), (n) Human Schizophrenia (Control n=167, Schizophrenia n=131; Fromer et al., 2016) and (o) Bipolar disorder (Control n=65, Bipolar disorder n=40; Fromer et al., 2016); unpaired two-sided Wilcox rank test. In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively. The center line represents the median. (p) Mean scaled enrichment of GWAS hits (MAGMA calculated p-value < 0.05) from FTD GWAS (Ferrari et al., 2014), PSP GWAS (Hoglinger et al., 2011) and AD GWAS (Lambert et al., 2013) in various TPR50 modules (n=15 modules). NAI module enrichment for AD risk genes was still significant after omitting APOE from the analysis (Supplementary Table 4d).

We also performed mass-spectrometry-based, label-free, quantitative proteomics in an independent set of GRN-positive and GRN-negative FTD frontal cortex (Supplementary Fig. 5p, Supplementary Table 2c–d, Online Methods). We found similar NAI and NAS expression patterns at the mRNA and protein level (Fig. 3e) both in GRN-positive and GRN-negative FTD samples, as compared with age-matched controls (Fig. 3f). Taken together, our analyses show that the NAS and NAI co-expression modules are conserved across multiple model systems, and generalize from mouse to humans, reflecting convergent RNA co-expression and protein-level changes in FTD.

Examination of RNA sequencing and proteomic post mortem brain datasets from AD patients showed preservation of both the NAI and the NAS modules at the mRNA (Supplementary Fig. 5e, Fig. 3g–h, Supplementary Table 4a), and at the protein level (Fig. 3i). Notably, both modules are not dysregulated in cases of pathologic aging without dementia (p=0.77) (Fig. 3j), indicating that they are not related to the presence of Aβ amyloid plaques alone[31]. The NAS and NAI modules are also dysregulated in C9orf72-related and sporadic ALS (Fig. 3k) and PSP (Fig. 3l), but not in major depression (Fig. 3m), schizophrenia (Fig. 3n) or bipolar disorder (Fig. 3o). These results support the relevance and specificity of these changes to allied human neurodegenerative syndromes, but not other non-degenerative forms of neuropsychiatric disorders.

Assessment of Genetic Risk within Modules

Transcriptomic changes may be the cause or consequence of disease, so integrating these data with genetic risk data provides a means to further understand their relationship to potential causal mechanisms (e.g.[32,33]). Several AD candidate genes mapped to the NAI module including APOE, CLU, PICALM, C1q and TREM2 (Online Methods; Fig. 3p, Supplementary Table 1i); In contrast to AD, FTD and PSP GWAS data both showed significant enrichment in the NAS module, which contains risk genes such as SLC32A1, NSF and ELAVL2 (Fig. 3p). These differential module enrichments indicate that genetic risk for AD, FTD and PSP likely operates via distinct causal pathways that may converge on common downstream neuronal and glial-immune processes.

Identification of potential miRNA drivers

Based upon the strong co-expression relationships observed in the NAS and NAI modules, we hypothesized potential co-regulation by miRNA. We identified two miRNA modules associated with the transgenic condition (Bonferroni corrected P < 0.05, Fig. 4a, Supplementary Table 3): the miM12 module that is anti-correlated, and the miM16 module (Fig. 4b), which is positively correlated with the transgenic condition in affected regions, but not cerebellum (Fig. 4c and Supplementary Fig. 6a–e). The miM16 module eigengene is strongly anti-correlated with the NAS module eigengene (Fig. 4d and Supplementary Fig. 6f), and predicted targets of miR-203, an miM16 hub, are highly enriched (see Online Methods) in the NAS module (Fig. 4e and Supplementary Fig. 6g, h). Consistent with its potential regulatory role, miR-203 is robustly upregulated in transgenic mice at six months in disease-affected areas (Fig. 4f and Supplementary Fig. 6i–k), and in the frontal cortex in both tau-positive and tau-negative FTD (Fig. 4g, h). Protein levels of its targets are also significantly decreased in both granulin positive and negative FTD (p<0.005, Fig. 4i). These data across mouse and humans suggested that miR-203 may be acting as a driver of the NAS neurodegeneration-associated transcriptional program.

Figure 4:

miRNA Co-expression Network Analysis

(a) Signed association (Pearson correlation) of the miRNA module eigengenes with transgenic condition. Modules with positive values indicate increased expression in transgenic mice; modules with negative values indicate decreased expression in transgenic mice. Red dotted lines indicate Bonferroni corrected p<0.05 for multiple comparisons (n=16 modules) using p-values obtained from two-sided test for Pearson correlation. (b) miRNA co-expression network plot of the miM16 module showing hub miRNAs in the center. Large nodes indicate top five hub miRNAs. (c) Trajectory of the miM16 module in the cortex (unpaired two-sided Wilcox rank test, n=6 mice/group). (d) Multidimensional scaling plot illustrating correlations between module eigengenes of the mRNA and miRNA modules. Colors indicate bi-weighted mid-correlation (R) values. (e) Enrichment of selected miM16 module miRNA predicted targets in mRNA modules. All enrichment values (odds ratio) with FDR<0.05 and OR>2 are shown (fisher’s two-sided exact test, ***FDR<0.005). For a full list of enrichments refer to Supplementary Fig. 6g. Targetscan database was used for miRNA target prediction. (f) Trajectory of log2 expression of miR-203 in the cortex (unpaired two-sided Wilcox rank test, n=6 mice/group). (g,h) miR-203 expression in the cortex (g, Control n=14, Tau-positive FTD n=10, Tau-negative FTD n=6) and cerebellum (h, Control n=10, Tau-positive FTD n=7) of human FTD and control samples. Tau-positive FTD (FTD-Tau Pos.) and tau-negative FTD (FTD-tau Neg.) are shown are shown for cortex. (Unpaired two-sided Wilcox rank test). (i) Module eigengene of predicted targets of miR-203 expression in human FTD and control protein samples in the cortex. Control=8, FTD GRN-Pos=6, FTD GRN-Neg=10, Progranulin-positive FTD (FTD GRN Pos.) and progranulin-negative FTD (FTD GRN Neg.) are shown (unpaired two-sided Wilcox rank test). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median.

miR-203 regulates the neuronal mRNA module and causes cell death

To test the predicted causal relationship between miR-203 expression and NAS module down-regulation, we acutely over-expressed miR-203 (Fig. 5a) in primary cortical mouse neuronal cultures (Online Methods, Supplementary Fig. 7a), observing downregulation of its predicted targets (Fig. 5b) and the NAS module more broadly (Fig. 5c, Supplementary Fig. 7b, Supplementary Table 4a). We also observed increased apoptosis (Fig. 5d), peaking at DIV8 (1.85-fold; Fig. 5e), concurrent with peak NAS module downregulation. We further experimentally validated several predicted direct targets of miR-203 (Bcl2l2, Dgkb, Mapk10, Vsnl1) by luciferase reporter assays (Fig. 5f) and western blot (BCL2L2 and VSNL1l; Fig. 5g, see Supplementary Fig 9 for uncropped blots).

Figure 5:

Overexpression of miR-203 in vitro and in vivo

(a) Schematic of lentiviral vector used for in vitro studies. (b-c) Trajectory of the miR-203 target genes (b) and NAS module eigengene (c) in uninfected primary cortical cultures or cultures infected with either miR-203- or sc-miRNA-lentiviral construct (n=4/group, unpaired two-sided Wilcox rank test). (d) Representative images of TUNEL staining from three independent experiments in mouse primary cortical neurons overexpressing miR-203 or sc-miRNA (control) at DIV6, DIV8 and DIV10 days. Green: GFP (infection), Red: TUNEL and Blue: Hoechst. Scale bar=25μm. (e) Quantification of TUNEL staining. The percent of TUNEL-positive cells in miR-203 overexpressing cultures were normalized at each time point to the average percent of TUNEL-positive cells in sc-miRNA control (error bars = SEM, unpaired two-tailed t-test, n=60 cells for DIV6, n=50 for DIV8 and n=60 cells for DIV10). The center line represents the mean. (f) Normalized luminescence of luciferase reporter assay. Luciferase vectors containing 950 bps of the 3’UTR sequence of Bcl2l2, Dgkb, Mapk10, Vsnl1 and Npepps genes were co-transfected with 20nM miR-203 or control mimics in HEK293T cells and assayed after 24–48hrs (error bars = SEM, unpaired two-tailed t-test, n=12/group). The center line represents the mean. (g) Representative immunoblots from three independent experiments and quantification of BCL2L2 and VSNL1 protein levels from total cell lysates isolated from primary mouse cortical neurons overexpressing miR-203 or sc-miRNA (error bars = SEM, paired two-tailed t-test). The center line represents the mean. Uncropped blots are shown in Supplementary Fig. 9. (h) Schematic representation of experimental design and timeline to overexpress or inhibit miR-203 in C57BL/6 wild type or Tg4510 tau transgenic mouse frontal cortex using AAV2/9 system. (i) NAS module eigengene expression in GFP-positive cells overexpressing miR-203 or sc-miRNA at 3 or 6 weeks after AAV injection in C57BL/6 wild type mice (Unpaired two-tailed Wilcox rank test, n=6/group). (j) Expression of genes involved in positive regulation of apoptosis (GO term ID: 0043065) and negative regulation of apoptosis (GO term ID: 0043066) in GFP-positive cells overexpressing miR-203 or sc-miRNA at 3 and 6 weeks after AAV injection in C57BL/6 wild type mice (unpaired two-tailed Wilcox rank test, n=6/group). (k) Caspase-8 (CASP8) intensity of GFP-positive cells overexpressing miR-203 were normalized at each time point to the average Caspase-8 intensity of GFP-positive cells overexpressing sc-miRNA control (unpaired two-tailed Mann-Whitney test, n=464 cells for 3-weeks control, n=937 cells for 3-weeks miR-203, n=629 cells for 6-weeks control, n=1441 cells for 6-weeks miR-203; three independent biological replicates/condition). The center line represents the mean and error bars showing SEM (l-m) Module eigengene of miR-203 targets (l) or NAS module (m) in GFP-positive cells overexpressing miR-203 6 weeks after AAV injection (unpaired two-tailed Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). (n) Caspase-8 (Casp8) intensity of GFP-positive cells overexpressing miR-203 in Tg4510 Tau transgenic mice were normalized at each time point to the average Caspase-8 intensity of GFP-positive cells overexpressing sc-miRNA control (unpaired two-tailed Mann-Whitney test, n=507 cells for sc-miRNA and 2793 cells for miR-203; three independent biological replicates/condition). The center line represents the mean and error bars showing SEM. (o-p) Module eigengene of miR-203 targets (o) or NAS module (p) in GFP-positive cells overexpressing sc-TuD control or miR-203 TuD 6 weeks after AAV injection in Tg4510 Tau transgenic mice (Unpaired two-tailed Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median. (q) C-fos intensity of GFP-positive cells expressing TuD-miR203 in Tg4510 Tau transgenic mice were normalized at each time point to the average c-fos intensity of GFP-positive cells expressing sc-TuD control (Unpaired two-tailed Student’s T-test, n=50 images (control) and n=60 images (miR203); n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). The center line represents the mean and error bars showing SEM.

To control for potential confounding by supra-physiological over-expression in vitro, we moderately overexpressed miR-203 (approximately 2-fold; Online Methods) or a scrambled miRNA control (sc-miRNA) (Fig. 5h) targeting neurons in wild-type mice in vivo (Supplementary Fig. 7c–f). The NAS module is highly preserved (Supplementary Fig. 7g, Supplementary Table 4a) and down-regulated in neurons overexpressing miR-203 after six weeks (Fig. 5i). Strikingly, the level of miR-203 expression correlates remarkably well with downregulation of both miR-203 predicted targets (R2=0.81, slope=−2.3; Supplementary Fig. 7h) and NAS module genes (R2=0.81, slope=−2.2; Supplementary Fig. 7i). Furthermore, apoptotic pathways are altered in miR-203-overexpressing neurons (, Figure 5j, Supplementary Table 4b), confirmed by Casp8 protein expression (Figure 5k, Supplementary Fig. 7j). Additionally, miR-203 over-expression in one-month-old Tg4510 Tau transgenic cortex (Supplementary Fig. 8a–b) causes downregulation of predicted miR-203 targets (Fig. 5l), the NAS module (Fig. 5m), a significant increase in CASP8 protein (Fig. 5n) and apoptotic pathways (Supplementary Fig. 8c–d) and genes involved in calcium signaling and neuroactive ligand receptors (Supplementary Fig. 8e–f).

To further investigate the causal role of miR-203, we inhibited miR-203 expression in neurons using a Tough Decoy (TuDs) AAV system[34] in vivo in Tg4510 mice (Fig. 5h). Both short term (6-weeks, Supplementary Fig. 8g–h) and longer term (4-months, Supplementary Fig. 8i–j) inhibition of miR-203 oppose the pathologic transcriptional changes in predicted miR-203 targets (Supplementary Fig. 8k, Fig. 5o) and the NAS module broadly (Supplementary Fig. 8l, Fig. 5p) including genes involved in calcium signaling and neuroactive ligand receptor pathways (Supplementary Fig. 8m–n), and apoptosis (Fig. 5q, Supplementary Fig. 8o–p). Reducing miR-203 did not affect phospho-tau immunostaining (Supplementary Fig. 8q–r), consistent with its acting downstream of tau pathology. Previous work has shown that five-month-old Tg4510 mice show reduced neocortical network activity based on c-fos immunoreactivity[35]. C-fos intensity was also significantly increased in mutant neurons expressing the miR-203 TuD construct, as compared to those expressing the control construct (Fig. 5q, Supplementary Fig. 8s, p-value <0.0001), suggesting functional rescue. These findings demonstrate that inhibition of miR-203 can reverse expression patterns of both the NAS module and apoptotic pathway genes associated with disease state in vivo, while overexpression of miR-203 promotes downregulation of the NAS module and a neurodegenerative phenotype.

Pharmacological manipulation of NAS module genes

We reasoned that if the patterns of gene expression robustly associated with disease state across mice and humans represented causal phenomena as suggested by the GWAS enrichment, then their reversal should ameliorate the miR-203 mediated cell death phenotype. To test this, we screened the Connectivity Map (CMAP), a public compendium of cell line gene expression responses to drugs, to identify small molecules predicted to reverse NAS or NAI module changes observed across pathological conditions, albeit in non-neuronal cells. Four of the top 10 hits were histone deacetylase inhibitors (HDACs, Supplementary Table 4c), which was significant (permutated p< 10−5).

We chose two of the predicted molecules, scriptaid, which was the top hit in CMAP (enrichment score = −0.969, Fig. 6a) and SAHA (vorinostat; enrichment score = −0.846), which is in human clinical trials[36]. As predicted, addition of scriptaid to miR-203-overexpressing neurons decreases neuronal death caused by miR-203 overexpression to levels equivalent to controls (unpaired t-test, Fig. 6b–c). We also found that 0.5μM SAHA decreases miR-203 induced cell death (Fig. 6d–e). SAHA upregulated NAS module genes (Fig 6f); however, SAHA showed toxicity at 1μM even in control neurons (Fig. 6f). Next, we treated human iPSC-derived neurons from control and FTD patients with SAHA, observing that SAHA upregulates NAS module genes in a dose-dependent manner (Fig. 6g). We therefore conclude that miR-203 directly regulates many NAS module hub genes and that the HDAC inhibitors scriptaid and SAHA, by targeting the NAS and NAI module genes and normalizing their expression, can ameliorate miR-203-induced pro-neurodegenerative changes.

Figure 6:

Small Molecule Inhibition of miR-203 Induced Cell Death in vitro.

(a) Co-expression based network plot of scriptaid targets from the Connectivity Map database with genes in the NAS and NAI modules. Genes which are upregulated with scriptaid treatment are connected by red edges and those which are downregulated are connected by green edges. Nodes size represents centrality within the network. (b, d) Representative images from three independent experiments of TUNEL staining in mouse primary cortical neurons overexpressing miR-203 or sc-miRNA (control) at DIV7 treated with DMSO control and 1uM or 2.5uM scriptaid (b) or DMSO control and 0.5uM or 1uM SAHA (d) for 24hrs prior (green: GFP, red: TUNEL, blue: Hoechst, scale bar=25μm). (c, e) Quantification of TUNEL positive cells from three independent experiments treated with scriptaid (c) and SAHA (e) (unpaired two-tailed t-test, n=15/group for SAHA experiments and n=20/group for scriptaid experiments). Scriptaid treatment does not alter miR-203 overexpression levels (Supplementary Fig. 9, indicating that the decrease in cell death is not a result of changes in viral infection, miRNA processing or regulation of miR-203 by scriptaid. The center line represents the mean and error bars showing SEM. (f) Trajectory of the NAS module eigengene in cultures infected with either miR-203- or sc-miRNA-lentiviral construct treated with either DMSO, 0.5 or 1uM SAHA (unpaired two-tailed Wilcox rank test, n=6/group). (g) Trajectory of the NAS module eigengene in human iPSC derived neurons from control and A152T Tau patients treated either DMSO, 0.5, 1 or 2.5uM SAHA (unpaired two-tailed Wilcox rank test, n=6/group). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median.

Discussion

The mechanisms of neurodegenerative diseases are still poorly understood, a factor that impedes the design of therapeutic interventions to treat these disorders. Rather than focusing on a priori candidate genes, we pursued a systematic, genome-wide, functional genomics approach. Such studies have not been realized in the study of most neurodegenerative diseases, with a few notable exceptions[37-39]. By taking genetic background into consideration, representing diverse causal mutations in transcriptomic network analyses, and validation in a dozen independent data sets, we identify disease-relevant gene networks representing specific molecular pathways. It is noteworthy that modules identified in the initial individual F1 crosses that are not conserved and dysregulated across all three strains, are also not conserved in other mouse models, or human brain (Supplementary Fig 10a–c), supporting the approach that we have taken.

In functional genomic studies, each experimental system on its own, whether post mortem human tissue, in vitro or in vivo mouse models, has significant limitations[40]. Here we show how a multi-stage systems biology analysis coupled with substantial in vitro experimental validation, provides a framework for discovery of new disease mechanisms and therapeutic targets. We leverage the benefits of animal models of human disease, which permit temporal assessment of changes prior to substantial atrophy and cell loss, while mitigating the potential bias introduced by a single genetic background. We also consider the principles of regional vulnerability and disease trajectories[12,41] and find that gene networks associated with transgenic status are altered in brain areas that show signs of neurodegeneration. Conversely, in the cerebellum, a brain region that is spared in AD and FTD, the expression of the NAS and NAI modules, as well as miR-203, do not change significantly.

A key issue in any analysis of gene expression in disease is that changes in gene expression may be either a cause or consequence of the disorder. To advance our work towards identifying causal drivers, we experimentally validate a putative miRNA driver of the altered transcriptional networks, miR-203, which has previously not been implicated in neurodegeneration, and we show reversal by gene network-predicted pharmacological regulators. We also find significant enrichment of common genetic risk for FTD and PSP in the neuronal NAS module, and enrichment for AD GWAS signals in the NAI module, which is enriched in glial and immune genes. These data further support the potential causal role of these modules, and suggest different causal pathways in AD versus FTD and PSP, the latter of which are considered primary tauopathies[3,4]

It is important to note that both expression changes in the disease associated gene modules NAS and NAI as well as the upregulation of miR-203, though clearly a point of convergence, are most likely downstream of initial pathological insults, whether dysfunction of tau or other major risk genes. From this perspective, these modules represent robust, highly conserved biological processes that provide a crucial link in the chain between an initial causal genetic insult and neuronal death and inflammation. As we show, long-term inhibition of miR-203 in the Tg4510 model, while upregulating anti-apoptotic pathways, downregulating pro-apoptotic pathways and increasing neuronal activity as measured by c-fos activation, does not reduce the phospho-tau burden (Supplementary Fig. 8q–r). Though pharmacologically targeting these gene networks, or miR-203 itself, may not reverse tau aggregation or Aβ plaque formation, it may provide a means to uncouple dysfunction of aggregation prone proteins from inflammation, neurodegeneration and their behavioral outcomes – dementia. Supporting this possibility is the fact that neither the NAS nor the NAI modules are dysregulated in patients that show pathological aging, in which significant amyloid plaque density, but no dementia, is observed. Therefore, developing therapies that reverse the changes in the disease associated gene co-expression modules may be an opportunity to impede neurodegeneration across a broad spectrum of dementias without having to target unique upstream pathological events. We hypothesize that regulators of these highly conserved transcriptional programs, such as miR-203, may therefore be new therapeutic targets. Efforts to target miRNAs to treat diseases are already underway[42] and traditional anti-sense oligo based targets have reached clinical trials[43], while small molecule modulators represent a new frontier in miRNA-based therapeutics[44]. In this regard, identifying regulators of the NAI module is also likely to be a fertile avenue for future studies aimed at therapeutic development.

Methods

Animals

TPR50 mice were generated in the C57BL/6J background and maintained at the Shonan research center, Fujisawa as described previously[13]. TPR50 mice overexpress the TauP301S (4R2N isoform) under the control of mouse prion promoter. Mice were outbred to DBA/2J and FVB/NJ strains. Heterozygous F1 male mice (at three- and six-months-of-age) were used to harvest brain tissues (cortex, hippocampus, brain stem and cerebellum) for immunohistochemical staining and RNA isolation. Wild type littermates of same genetic background and age were used as controls. All procedures involving animals were performed in accordance with the University of California, Los Angeles (UCLA) animal research committee, the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals and by the Experimental Animal Care and Use Committee of Takeda Pharmaceutical Company Limited, Japan. Wild type C57BL/6J mice were directly purchased from Jackson Labs (Bar Harbor, USA). Methods describing total RNA isolation from brain tissue, tissue fixation and staining are described in the Online Methods.

Human postmortem brain samples

Post-mortem brain tissue from human FTD and control individuals used for mRNA- and miRNA-seq were acquired from the Emory ADRC brain bank at the Emory University School of Medicine. Human postmortem tissues were acquired under Emory University Institutional Review Board (IRB) with consent from family. Sample acquisition protocols approved by UCLA IRB were followed and samples were de-identified prior to acquisition. Brain samples and individual level metadata is available in Supplementary Table 2b. Proteomics samples were obtained from Penn ADRC as part of the Accelerating Medicine Partnership for Alzheimer’s Disease (AMP-AD) and sample level metadata is available in Supplementary Table 2c.

Immunohistochemistry/Immunofluorescence

Paraffin-embedded mouse brain samples were briefly deparaffinized and blocked in 5% serum. Sections were stained overnight with Phospho-tau (Ser202, Thr205) AT8 clone (Thermo Fisher, MN1040) at 1:250 dilution followed by incubation in biotinylated secondary antibody and imaged under microscope. For immunofluorescence, sections were stained with an antibody raised against Glial Fibrilary Acid protein (GFAP, Dako, Z033429, 1:500) or Iba-1 (Wako, 019–19741, 1:500) followed by incubation in Alexa Fluor 488 or 594 secondary antibody (1:500, Invitrogen) for 1h at room temperature and counterstained with DAPI. Images were obtained either using Zeiss Axio-Imager or on a digital slide scanner NanoZoomer 2.0-HT, Hamamatsu Photonics K.K.

mRNA and miRNA sequencing

1μg of total RNA was used for mRNA- and miRNA-seq. For mRNA-seq polyA-selected mRNAs were processed for unstranded libraries using the TruSeq v2 kit (Illumina) that underwent 50bp paired-end sequencing on an Illumina HiSeq 2500 machine. Paired-end reads were mapped to the reference mouse NCBIM37 genome using Tophat2 and Ensembl release 67 annotations (May 2012 data freeze). Count level data was quantified using union gene models with HTseq-counts. For miRNA sequencing a second sample from the same total RNA fraction was used and libraries were prepared using Truseq small-RNA library prep kit (Illumina). The 50bp single-end reads were processed for miRNA expression using the miRDeep2 analysis pipeline. For additional information on sequencing, read alignment parameters and data-preprocessing, please see the Online Methods.

Co-expression network analysis

Co-expression network analysis was performed in R using WGCNA package[16]. To identify mRNA modules that are conserved across various transgenic mouse strains, we used a consensus network analysis (cWGCNA) approach. To identify common expression patterns across genetic backgrounds we created a signed consensus network for each brain region individually by calculating component-wise minimum values for topological overlap across different mouse strains. For miRNA data we performed WGCNA analysis using signed networks. Please see Online Methods for more details.

Mouse embryonic cortical cell culture

Cortical tissue from E15 C57BL/6J embryos was harvested, dissected and washed in ice-cold HBSS (Invitrogen; 14170–112). Tissue was incubated in 0.25% trypsin (Invitrogen; 15090–046) in the presence of DNase I (Roche; 10104159001) at 37°C for 10min. Tissue was washed with cold HBSS and titurated in plating media (Neurobasal Media (Invitrogen; 21103–049), 20% Horse Serum (Invitrogen; 26050–088), 25mM Sucrose and 0.25% GlutaMax (Invitrogen; 35050–061)) in the presence of DNase I. Dissociated cells were centrifuged at 125g for 5min at 4°C, resuspended in plating media, counted and plated in poly-L-lysine (Sigma P1274) coated plates at a density of 300,000 cells/mL. If cells were infected with miRNA-203 or scrambled control lentiviral constructs (Systems Bioscience; MMIR203-PA-1), 10 IFU/cell was added to the suspended cells prior to plating. Plating media was replaced 24hrs after plating with Neurobasal Medium supplemented with 1% B27 (17504–044) and 0.25% GlutaMax. For virus production and titration methods please see Online Methods.

Stereotactic Injection of AAV2/9 into mouse brain

Mice were anesthetized using isoflurane and head-fixed over a heating pad set to 37°C. Craniotomy was performed according to approved procedures (animal protocol number#2000–159): Eyes were covered with artificial tear drops, the scalp was bisected and the exposed area sterilized on a stereotactic frame (Model 900, David Kopf Instruments). The skull was thinned using a dremel drill and 1μl of AAV 2/9 (1 × 1013 vector genomes (Vg)/ml of hSyn-miR203-hSyn-eGFP; 1 × 1013 Vg/ml of hSyn-scrambled:miRNA-hSyn-eGFP (Signagen Labs, Rockville, USA)) was injected bilaterally into the frontal cortex (anterior/posterior: +2.1mm; mediolateral: ±1.5mm; dorsal/ventral: −2.5mm) at a rate of 100nl/min rate using Hamilton 5ul syringe (87930 Hamilton , Reno, USA. Before each injection the pipette was lowered 0.1mm beyond intended target depth and held in place for 3 minutes to create space for injected solution, while after each injection the pipette was held in place for 7 min before retraction to prevent leakage. The incision was sterilized, glued and post-operative antibiotics (Amoxil, 50mg/ml) were administered for 7 days following surgery.

Data and Code Availability

miRNA- and mRNA-seq data from TPR50 Tau mice, microarray data on PS19 hippocampus, microarray data on overexpression of miR-203 in vitro, RNA-seq on sorted mouse neurons, and RNA-seq data with SAHA available at the NCBI Gene Expression Omnibus database under Gene Expression Omnibus accession number #GSE90696. Human FTD miRNA- and mRNA-seq data is available at https://www.synapse.org/#!Synapse:syn7818788. Human UPenn FTD Proteomics data is available at https://www.synapse.org/#!Synapse:syn9884357. Custom code used for the analysis can be accessed using this link in github - https://github.com/dhglab/Identification-of-evolutionarily-conserved-gene-networks-mediating-neurodegenerative-dementia

Online Methods

Sample Collection and RNA Isolation

TPR50 mouse at 3- and 6-months-of-age were sacrificed using cervical dislocation. Cortex, hippocampus, brain stem and cerebellum were dissected on ice-cold PBS and flash frozen to minimize RNA degradation. RNA was isolated from different brain regions of interest from each sample using the miRNeasy kit with on-column DNase digestion (Qiagen) using a QIAcube automated system. For each RNA sample, RNA concentration was determined using Ribogreen assay (Invitrogen) and RNA integrity was quantified using the RNA Integrity Number (RIN)[1] on an Agilent 2200 TapeStation.

Frozen human brain samples obtained from Emory brain bank were dissected on dry ice in a dehydrated dissection chamber to reduce degradation effects from sample thawing or humidity. 50–100mg tissue samples were used to isolate RNA using miRNeasy kit with on-column DNase digestion and samples with RIN > 2 were further used.

Mouse primary cortical neurons were washed with ice-cold RNase-free PBS and Qiazol was immediately added to the 6-well plates. RNA was isolated using miRNeasy kit following the manufacture’s protocol.

Several datasets have been used throughout the study, details of which can be found in Supplemental Table 2a. Mouse datasets consist of – Tg4510[2], PS2APP[3], CRND8[4], APP/PS1[5], GRN[6]. Human postmortem data consist of AD Temporal cortex[7], AD Frontal cortex[8], AD proteomics[9], Pathological aging without dementia[7], ALS[10], PSP[7], major depression[11], schizophrenia and bipolar disorder[12]

ELISA

Total tau and pT231 tau content were measured by commercial tau ELISA kits according to the manufacturer’s instructions (total tau - KHB0041; pT231 tau - KHB8051, Invitrogen). Briefly, standards, RIPA-soluble or sarkosyl insoluble samples were applied to the ELISA plate. After washing, a biotin-conjugated detection antibody was applied. The positive reaction was enhanced with streptavidin-HRP and colored by TMB. The absorbance at 450 nm was then measured and the concentration of tau protein was calculated from the standard curve.

RNA-seq Library Preparation

For TPR50 mouse data, 1μg of total RNA was used to obtain poly-A selected mRNAs, which were processed for unstranded libraries using TruSeq RNA Sample preparation (v2) kit. Briefly, poly-A tail containing mRNA was isolated using oligo-dT attached magnetic beads. mRNA is fragmented, and first strand cDNA is generated using SuperScript II reverse transcriptase followed by second strand cDNA generation. Barcodes and adapters were added and subsequent steps followed the TruSeq protocol to generate fragment sizes (150–500bp, mean 250bp). The libraries were quantified with the Quant-iT PicoGreen assay (Life Technologies) and validated on an Agilent 2200 TapeStation system. Libraries were multiplexed 24 samples per lane and each lane was sequenced several times to get an average read-depth of 40–50 million reads per sample on a HiSeq2500 instrument using high output mode with standard chemistry and protocols for 50bp paired end reads.

For postmortem human data, 1μg of total RNA was used for rRNA depletion with the RiboZero Gold kit (Illumina). Remaining RNA was size selected using AMPure XP beads (Beckman Coulter) and standard libraries were prepared following Illumina’s TruSeq protocols for 50bp paired end reads. Libraries were sequenced at an average read-depth of 60–70 million reads per sample on a HiSeq2500 instrument using rapid mode.

We used SMART-Seq v4 ultra low input RNA kit (Clontech) for the library preparation of FACS sorted adult mouse neurons. RNA from FAC sorted cells (~1000 cells) were isolated using Nucleospin RNA XS kit (Clontech) with on-column DNase digestion following the manufacture’s protocol and the RNA was eluted in 10μl of RNase-free water. rRNA was depleted from the eluted total RNA and cDNA was synthesized using Primescript Reverse transcriptase. Sequencing adapters were added and library was amplified using SeqAmp DNA polymerase. Strand-specific library was size selected using AMPure XP beads and validated for fragment sizes (200–1000bp, peaking at 300bp). Libraries were multiplexed 12 samples per lane and each library was sequenced several times to get an average read-depth of 30–40 million reads per sample on a HiSeq2500 instrument using rapid mode.

We used TruSeq Stranded RNA-seq (Illumina) for the library preparation of vorinostat-treated neurons. Briefly cell pellet was collected from 6-well plates and RNA isolated using miRNeasy kit with on-column DNase digestion (Qiagen). 100ng of total RNA was used for downstream library preparation using manufacture’s recommended protocol. Libraries were multiplexed 24 samples per lane, and each library was sequenced several times to get an average read-depth of 30–40 million reads per sample on a HiSeq4000 instrument.

mRNA-seq Read Alignment and Processing

To analyze TPR50 mouse mRNA-seq data, the paired-end raw reads were mapped to the reference mouse NCBIM37 genome using Ensembl release 67 (May 2012 data freeze) annotations with Tophat2 (v2.0.5)[13] with novel junction option turned off. Aligned reads were sorted and alignments mapped to different chromosomes were removed from the BAM file using samtools[14]. Gene expression levels were quantified for all the samples using union gene models with HTSeq-Counts (v0.5.4)[15] which uses uniquely aligned reads. Genes were included in the analysis if they were expressed in 80% of the samples with >10 reads. It is a common practice in RNA-seq experiments to filter low expressed gene[16,17]. There are a few reasons for this: 1) the low expressed genes are not reliably quantified and thus vary a lot between biological replicates. 2) None of the models such as the negative binomial distribution from DESeq2 (or edgeR) or log2 transformation followed by linear regression (which we used) can effectively model such highly variable genes. Due to these reasons, we remove these genes mainly to reduce noise in gene-expression from RNA-seq experiments.

The resulting read counts were normalized for library size using cpm function of egdeR package in R[18] to get fragments per kilobase million mapped reads (FPKM) values. Log2 transformed FPKM values were quantile normalized using betweenLaneNormalization function of EDAseq package in R[19].

For postmortem human data, the paired-end strand-specific reads were mapped to the reference human GRCh37.73 (Gencode 19) genome using rna-STAR aligner (v2.4.2a)[20]. Count level data were quantified using union gene and union exon models with HTSeq-Counts (v0.6.1) and genes were kept if they were expressed in 80% of the samples with >10 reads with HTSeq union exon and union gene quantifications. Read counts were log2 transformed and GC content, gene length, and library size normalized FPKM values using the cqn package in R[21].

FACS sorted neuronal RNA-seq data were processed similar to the human RNA-seq data except the paired-end reads were mapped to the reference mouse NCBIM37 genome using Ensembl release 67 annotations using rna-STAR aligner (v2.4.2a).

miRNA-seq Read Alignment and Processing

Mouse and human miRNA-seq read alignments and quantifications were performed using miRDeep2 package[22]. We used mapper.pl script to map the 50bp single end reads to the reference genome – NCBIM37 for mouse and GRCh37 for human. We then used quantifier.pl script of the package to map the reads to annotated miRNAs, which were downloaded from miRbase (v20 for mouse data and miRbase v21 for human data) website. The script gives a count-level quantification of each known miRNA and we used the quantile-normalized values as provided by the script for downstream analyses.

Microarray - mouse embryonic cortical cells

Total RNA isolated from E15 cortical cultures using QIAGEN miRNeasy mini kit (QIAGEN; 217004) were processed on an Illumina mouse Ref8 v2 beadchip microarrays following the manufacture’s protocol. Microarray data analysis was performed using R and Bioconductor packages. Raw expression data were log2 transformed and normalized by quantile normalization. Probes were considered robustly expressed if the detection P value was 0.05 for at least half of the samples in the data set. Probes were re-mapped to mouse Ensembl gene IDs (v67; May 2012 data freeze) for comparison with RNA-seq data.

Differential Gene Expression

We used principal component (PC) analysis of the normalized gene-expression data to understand the biological and technical covariates effecting the data. For TPR50 dataset, we found that none of the first five gene-expression PCs correlated significantly with batch or sequencing PCs, and thus we only performed differential expression within a brain region and separately for each age-group using the normalized FPKM values in the linear regression model in R as follows:

For all other mouse datasets, RIN, first two sequencing PCs, batch, experimental batch were used as covariates in the linear regression model. For human data, additional traits like age, sex and PMI were also used for the analysis. We also performed DGE for TPR50 dataset using DESeq2 R package[23] using raw counts data from HTSeq-Counts following the default parameters and taking Transgenic Condition, Mouse Strain, Sequencing Lane Batch and RIN as covariates for DESeq2 regression.

mRNA Weighted Co-expression Network Analysis

Co-expression network analysis was performed using a user-friendly R WGCNA library[24]. We wanted to investigate mouse brain co-expression networks that are disease-specific but are independent of genetic backgrounds (C57BL6/J, DBA and FVB). We thus used consensus network analysis approach, which provided us with groups of co-expressed genes (or modules) that are not affected by mouse genetic background. TPR50 mouse were generated using a BAC clone as transgene, which overexpressed P301S mutation in human MAPT gene (4R2N isoform). However, the BAC clone had 3 additional confounding genes – Prnd (ENSMUSG00000027338), Erv3 (ENSMUSG00000037482) and Rassf2 (ENSMUSG00000027339) in the construct, which were also overexpressed. To remove the effect of the confounding genes, we took the 1st principal component of the expression of the confounding genes and removed the genes, which are significantly correlated at FDR p<0.05. This resulted in removal of about 1000 genes, which are highly correlated with the confounding genes. We only used the 6-month time-point to construct the initial network. The resulting expression data for each brain region, except cerebellum, at the 6-month time-point was used to create consensus network separately and matchLabels function of WGCNA package was used to have consistent network labeling between brain regions.

In short, biweighted mid-correlations were calculated for all pairs of genes, and then a signed similarity matrix was created. In the signed network, the similarity between genes reflects the sign of the correlation of their expression profiles. The signed similarity matrix was then raised to power β to emphasize strong correlations and reduce the emphasis of weak correlations on an exponential scale. The resulting adjacency matrix was then transformed into a topological overlap matrix as described here[25].

Since we are primarily interested in finding co-expression patterns conserved across different mouse genetic backgrounds, we created a consensus network to identify common co-expression patterns across strains, following published methods[26]. After scaling for each individual network (consensus scaling quantile=0.2) a thresholding power of 12 was chosen (as it was the smallest threshold that resulted in a scale-free R2 fit of 0.8) and the consensus network was created by calculating the component-wise minimum values for topologic overlap (TO). Using 1 − TO (dissTOM) as the distance measure, genes were hierarchically clustered. Initial module assignments were determined by using a dynamic tree-cutting algorithm (cutreeHybrid, using default parameters except deepSplit = 4, cutHeight = 0.999, minModulesize = 100, dthresh=0.1 and pamStage = FALSE).

The resulting modules or groups of co-expressed genes were used to calculate module eigengenes (MEs; or the 1st principal component of the module). Modules annotated using the GOElite package[27]. We performed module preservation analysis using mRNA module definitions. MEs were correlated with different biological and technical traits like transgenic condition, strain, RIN, etc. to find disease-specific modules. Module hubs were defined by calculating module membership (kME) values which are the Pearson correlations between each gene and each ME. Genes with kME < 0.7 were removed from the module. Network visualization was done using iGraph package in R[28]. Module definitions from the network analysis were used to create synthetic eigengene for the 3-month time point and were used to understand the trajectory of various modules across time-points. For the cerebellum a data synthetic module eigengene was calculated using cortex co-expression definitions and plotted.

Cell-type Specific Regression Analysis

To ensure that the changes in module expression identified were not solely a product of neurodegeneration and gliosis, we employed a multivariate linear regression model to regress gene expression levels against transgenic condition and the first principal component (PC1) of the top one hundred cell-type specific marker genes[29] for five major cell types: neurons, astrocytes, oligodendrocytes, microglia and endothelial cells. This method removes the confounding effect of changes in cell composition (e.g. see[30,31]).

Module Preservation Analysis

We used module preservation analysis to validate co-expression in independent mouse and human datasets. Module definitions from TPR50 consensus network analysis were used as reference and the analysis was used to calculate the Zsummary statistic for each module. This measure combines module density and intramodular connectivity metrics to give a composite statistic where Z > 2 suggests moderate preservation and Z > 10 suggests high preservation[32].

Enrichment Analyses for Gene Sets

Gene set enrichment analysis was performed using a two-sided Fisher exact test with 95% confidence intervals calculated according to the R function fisher.test. We used p values from this two-sided approach for the one-sided test (which is equivalent to the hypergeometric p-value) as we do not a priori assume enrichment[33]. To reduce false positives, we used FDR adjusted p-values[34] for multiple hypergeometric test comparisons. For cell-type enrichment analysis we used already published mouse brain dataset[29]. The background for over-representation analysis was chosen as brain region expressed data from our RNA-seq data.

For miRNA binding site/target enrichment analysis, we downloaded predicted miRNA targets from mouse TargetScan database (v6.2) using only conserved predicted targets[35]. For background we used all genes expressed in the mouse genome.

Genes in network modules were characterized using GO Elite (version 1.2.5) using the brain region expressed genes as background[27]. GO Elite uses a Z-score approximation of the hypergeometric distribution to assess term enrichment, and removes redundant GO or KEGG terms to give a concise output. We used 10,000 permutations and required at least 10 genes to be enriched in a given pathway at a Z score of at least 2. We report only biological process and molecular function category output.

Protein-Protein Interaction Analysis

We used two protein-protein interactions (PPIs) resources - InWeb[36] and BioGRID[37]. Basic analysis was performed similar to methods published elsewhere[38]. Using the union of the two resources, the subset of compiled PPIs between genes in a module was extracted and all edges were counted. The PPI dataset was matrix multiplied with co-expression data from RNA-seq (as above) and edges, which were present both in PPI and co-expression datasets were eventually kept. This approach allowed us to infer tissue and species specificity in the PPI network. PPI network was then visualized using iGraph package in R. We also assessed whether the modules are interconnected by PPIs above chance using DAPPLE[36], which uses a within-degree within-node permutation method that allows us to rank PPI hubs by p value.

miRNA Weighted Co-expression Network Analysis

miRNA co-expression networks for each brain region, except the cerebellum, at the 6 month time point was constructed separately using the R package WGCNA as described previously[39]. A thresholding power of 12 was chosen (as it was the smallest threshold that resulted in a scale-free R2 fit of 0.8). The network was created by calculating the component-wise minimum values for topologic overlap (TO) and miRNAs were hierarchically clustered. Initial module assignments were determined by using a dynamic tree-cutting algorithm (cutreeHybrid, using default parameters except deepSplit = 2, cutHeight = 0.999, minModulesize = 40, dthresh=0.1 and pamStage = FALSE). Other steps were similar to mRNA co-expression analysis described above. Module definitions from the network analysis were used to create synthetic eigengenes for the three-month time point and were used to understand the trajectory of various modules across time-points.

Label-free Quantitative Proteomics

Label-free Quantitative proteomics were performed at the Emory Proteomics Core, Emory University, USA. Detailed methods were published elsewhere[9]. Briefly, postmortem frozen human brain samples were homogenized in urea lysis buffer (8M urea, 100 mM NaHPO4 buffer system, pH 8.5), including 5 μL (100x stock) HALT protease and phosphatase inhibitor cocktail (ThermoFisher, Cat# 78440) and further diluted with 50 mM NH4HCO3 to a final concentration of less than 2M urea and then treated with 1 mM dithiothreitol (DTT) at 25°C for 30 minutes, followed by 5 mM iodoacetamide (IAA) at 25°C for 30 minutes in the dark. Protein was digested with 1:100 (w/w) lysyl endopeptidase (Wako) at 25°C for 2 hours and further digested overnight with 1:50 (w/w) trypsin (Promega) at 25°C. Resulting peptides were desalted with a Sep-Pak C18 column (Waters) and dried under vacuum.

Brain derived tryptic peptides (2 mg) were resuspended in peptide loading buffer (0.1% formic acid, 0.03% trifluoroacetic acid, 1% acetonitrile) containing 0.2 pmol of isotopically labeled peptide calibrants (ThermoFisher, #88321). Peptide mixtures were separated on a self-packed C18 (1.9 um Dr. Maisch, Germany) fused silica column (25 cm × 75 mM internal diameter; New Objective, Woburn, MA) by a NanoAcquity UHPLC (Waters, Milford, FA) and monitored on a Q-Exactive Plus mass spectrometer (ThermoFisher Scientific, San Jose, CA). Elution was performed over a 120 min gradient at a rate of 400 nL/min with buffer B ranging from 3% to 80% (buffer A: 0.1% formic acid and 5% DMSO in water, buffer B: 0.1% formic and 5% DMSO in aceto- nitrile). The mass spectrometer cycle was programmed to collect one full MS scan followed by 10 data dependent MS/MS scans. The MS scans (300–1800 m/z range, 1,000,000 AGC, 150 ms maximum ion time) were collected at a resolution of 70,000 at m/z 200 in profile mode and the MS/MS spectra (2 m/z isolation width, 25% collision energy, 100,000 AGC target, 50 ms maximum ion time) were acquired at a resolution of 17,500 at m/z 200. Dynamic exclusion was set to exclude previous sequenced precursor ions for 30 s within a 10ppm window. Precursor ions with +1, and +6 or higher charge states were excluded from sequencing. The label free quantitation (LFQ) algorithm in MaxQuant[40] was used for protein quantitation. The quantitation method only considered razor and unique peptides for protein level quantitation. The LFQ intensities were log2 transformed for downstream analyses.

Lentivirus Production

10cm dishes were coated in poly-L-ornithine (Sigma) diluted in PBS and incubated at 37°C for 1hr-overnight. Plates were washed twice with PBS. HEK293T cells were seeded in DMEM supplemented with 10% Fetal Bovine Serum (FBS; Invitrogen). 24hrs after plating, media was replaced with prewarmed DMEM and HEK293T cells were transfected with 9ug purified lentiviral construct DNA, 9ug purified PAX2 DNA and 0.9ug purified VSV-G DNA diluted in 1mL OptiMEM (Invitrogen) and 60uL lipofectamine 2000 (Invitrogen) diluted in 1mL OptiMEM. 4–6hrs post-transfection, media was removed and replaced with 10mL prewarmed DMEM + 30% FBS. 16hrs after media exchange, media was replaced with prewarmed DMEM + 10% FBS. 48hrs following second media exchange, media was collected, centrifuged for three minutes at 3000g at room temperature and filtered through a .45um filter. Virus was concentrated using Lenti-X concentration (Clontech) and resuspended in Neurobasal A Medium (Invitrogen) supplemented with 5ug/mL polybrene. Virus was titered using Lenti-X p24 Rapid Titer Kit (Clontech), aliquoted and stored at −80°C until used.

miRNA expression

Expression levels of individual miRNAs of interest were quantified using the QIAGEN miscript kits. For low input miRNA quantification, miRCURY LNA Universal RT microRNA PCR (Exiqon) were used using primers specifically designed for miR-203.

TUNEL staining – mouse embryonic cortical cells

Embryonic cortical cultures were fixed with 4% Paraformaldehyde (PFA) (EMS; 15710), washed with Phosphate Buffered Saline (PBS) and permeabilized in 0.25% triton X-100 in PBS for 20min at room temperature. For TUNEL staining, cover slips were washed with deionized water and double-stranded DNA breaks were labeled using the Click-iT TUNEL AlexaFluor Imaging Assay (Invitrogen, C10246). Nuclei were labeled using the Hoechst 33342 dye provided in Click-iT TUNEL kit. Cover slips were mounted in Fluoromount-G (SouthernBiotech; 0100–01) and imaged using Zeiss Axio-Imager. The percentage of TUNEL+ cells in each condition was determined as the number of TUNEL+ cells over total nuclei detected.

Isolation of infected adult mouse neurons using flow cytometry.

Frontal cortex of adult male mice was dissected in ice-cold Hibernate A-lacking calcium (BrainBits; HA-Ca). Tissue was roughly chopped in HA-Ca, moved to 15mL Falcon tube and incubated in 5ml of activated papain (Worthington; LK003178; resuspended in HA-Ca) in the presence of DNase I for 30min at 37°C with continuous agitation. Papain digested tissue was triturated and cell debris was removed using ovomucoid (Worthington; LK003182; resuspended in HA-Ca) discontinuous density gradient centrifugation. The cell pellet was resuspended in 1.8mL Hibernate A-low fluorescence (BrainBits; HA-LF) to create a mononuclear cell suspension. To further reduce the amount of debris, the Miltyne myelin removal kit was used. Briefly, 200uL myelin removal beads (Miltyne; 130–096-733) were added to the cell suspension, incubated at 4°C for 15 minutes, 2mL of HA-LF was added and the cell suspension was centrifuged at 300g for 10 minutes at 4°C. The supernatant was aspirated and the pellet was resuspended in 1mL of HA-LF and applied to prepared LS columns (Miltyne; 130–042-401). Flow-through, as well as two 1mL washes with Ha-LF, were collected, centrifuged at 600g for 5min at 4°C and resuspended in 750uL HA-LF. Myelin depleted samples were labeled with live cell marker DRAQ5 (Fisher; 6225; 1uL/1mL of cell suspension) and dead-cell marker NucBlue® (Invitrogen; R37606; 2 drops/mL of cell suspension). The cells were FAC sorted on BD FACS Aria cell sorter gating for DAPI-/DRAQ5+/GFP+ cells. For each sample a maximum of 1000 cells were isolated and directly collected in 100ul of RA1 buffer with 2ul tris(2-carboxyethyl)phosphine (TCEP) and low mass RNA was isolated with on-column DNase digestion using Nucleospin RNA XS kit (Clontech) following the manufacture’s protocol.

Immunohistochemistry

Three or six weeks after AAV-mediated infection, mice were anesthetized with Pentasol (sodium pentobarbital, 40 mg/kg body weight) and perfused intracardially, first with PBS, followed by 4% PFA in PBS. The brains were removed, post-fixed overnight, and cryoprotected by immersion in 30% sucrose in PBS for 48hrs. Brains were embedded in optimum cutting temperature (OCT) compound, sectioned at a thickness of 20 μm on a cryostat, mounted and used for IHC using standard protocols. Incubations with primary antibodies were performed overnight at 4°C and with secondary antibodies for 60 min at room temperature. For primary antibodies, we used rabbit anti-GFP (1:500, Invitrogen) and rabbit anti-Casp8 (1:800, Cell Signaling). Images were acquired with a Zeiss LSM 780 laser scanning confocal microscope (fluorescence) or Zeiss Axio-Imager.

Western blotting

E15 cortical cultures were washed in PBS, dissociated in 0.25% prewarmed trypsin, washed in PBS and pelleted by centrifugation at 1,500g for 1min at room temperature. Cell pellets were resuspended in 50–70uL TNT buffer (150mM NaCl, 1mM EDTA, 50mM Tris pH 7.4, 1% Triton X-100) supplemented with protease and phosphatase inhibitors (Roche, 4693159001 & 4906845001). Lysates were incubated at 4°C for 30 minutes with constant agitation and centrifuged at 13,000g for 5 minutes at 4°C. Supernatant was transferred to new tube, protein concentration was measured using Bradford Reagent (BioRad, 500–0006) and equal volumes of Laemmle sample buffer (BioRad; 161–0737) supplemented with β-mercaptoethanol were added to the lysate. Lysates were briefly heated to 95°C, loaded onto 4–15% mini-PROTEAN TGX Precast gels (BioRad, 456–1086) and run at 100V for approx. 1.5hrs. Samples were transferred to PLVD membranes (BioRad, 162–0177), blocked with 5% milk or 3% BSA in PBST for 1hr at room temperature, probed with primary antibody overnight at 4°C, washed with PBST and probed with secondary antibody for 1hr at room temperature. Primary antibodies used: VSNL1 (1:1000; Millipore MABN762), DGKB (1:1000; ThermoScientific PA5–15416), β-actin (1:2000, Sigma A1978), BCL2L2 (1:1000, Origene TA302168). Mouse brain samples were dissected, washed in PBS and snap frozen in liquid nitrogen. Samples were homogenized in 0.5% sarkosyl buffer and sarkosyl soluble and insoluble fractions were separated. Western blots were probed with either total tau (Ab-3 antibody, RB-1429, Invitrogen) or phospho-tau (AT8 antibody, pS202/pT205, MN1020, Invitrogen).

Luciferase Assay

HEK293T cells were plated on poly-L-ornithine (Sigma, P3655) coated 96-well plates. When 70–80% confluent, cells were transfected with pmirGLO Dual-Luciferase miR Target Expression vector (Promega, E1330) containing 950pb target 3’UTR (centered around proposed miR-203 binding site) and 20nM mimic (Dharmacon; C-310523–05, C-310459–07) using lipofectamine 2000 (ThermoFischer Scientific; 11668027). 24–48hr after transfection, Dual-Glo Luciferase Assay (Promega, E2920) was used to measure both firefly and renilla luciferase-mediated luminesence using BioTek Synergy 2 plate reader (1 second integration time, 100 sensitivity). Firefly luciferase luminescence was normalized to renilla luciferase luminescence to control for transfection efficiency. Samples treated with miR-203 mimic were normalized to samples treated with scrambled mimic (control).

Connectivity Map (CMAP) Analysis

Top 300 salmon and turquoise hub genes (according to descending kME values) were used as input for CMAP database (https://portals.broadinstitute.org/cmap/). This signature was used to query drugs or small molecules and the similarity between this signature and more than 7,000 expression profiles for 1,309 compounds (reference signatures) in the cmap database were evaluated[41]. Enrichment of both salmon and turquoise module hub genes in the profiles of each treatment instance were estimated with a metric based on the Kolmogorov-Smirnov statistic, (a nonparametric, rank-based pattern-matching strategy) as described[41] and combined to produce a “connectivity score.” Mean connectivity scores across several cell lines were then ranked by increasing order of connectivity. Gene-expression data of the top scoring drug, scriptaid was extracted from CMAP database and used for finding genes differentially affected with the drug treatment.

Generation of the Human iNgn2 Neurons and Compound Treatment

The control 8330–8-RC1 and the Tau-A152T 19–5-RC6 neural progenitor cell (NPC) lines (derived from 8330–8 and 19–5 iPSC (generation of these iPSC lines are described in[42-44], respectively) were stably transduced with an inducible Neurogenin 2 (iNgn2) construct containing puromycin and blasticidin selection markers and a transactivator construct (derivation method described in Cheng et al. (2017)) resulting in the iNgn2-NPC stable cell lines. On day 0, each of the iNgn2-NPC cell line was plated at 9×105 cells/well into 6-well plates precoated with 20 μg/ml poly-ornithine (Sigma cat# P3655–50MG), and 5 μg/ml laminin (Sigma cat# L-2020–1MG) in neural medium (48% Neurobasal medium (Gibco cat# 21103–049), 48% DMEM/F12 (Gibco cat# 11330), 1% B27 (Gibco cat# 17504–044), 0.5% N2 supplement (Gibco cat# 17502–048), 0.75% GlutaMax (Gibco cat# 35050–061), 1% penicillin/streptomycin (Gibco cat# 15140–122), 0.5% MEM NEAA (Gibco cat# 11140), 50 μM 2-mercaptoethanol (BioRad cat# 161–0710), 0.2% bovine serum albumin (Sigma, cat# A7906–100 G), 2 μg/ml doxycycline (Clontech cat# 631311), 10 ng/ml BDNF (Peprotech cat# 450–02), 10 ng/ml NT3 (Peptrotech cat# 450–03) to start the neural induction. The cells were re-fed every other day with half replacement of the neural medium, with 1 μg/ml puromycin (Sigma cat# P8833–25MG) added on day 2 and 4, and cytosine arabinoside (AraC; Sigma cat# C6645) on day 6. Between days 13–14, the cells were treated for 24hrs with either 0.05% DMSO (vehicle control), 500nM, 1uM or 2.5uM SAHA in neural medium minus the growth factors and doxycycline. After the 24hr treatment, the 14-day iNgn2 neurons were harvested in TRIzol (Ambion cat# 15596018) with each replicate collected from two wells and stored at −80°C.

FACS sorting of adult mouse neurons.

Cortical cell suspension from adult mouse neurons were labeled with live cell marker DRAQ5 (Fisher; 6225; 1uL/1mL of cell suspension) and dead-cell marker NucBlue® (Invitrogen; R37606; 2 drops/mL of cell suspension). The cells were FACS sorted on BD FACS Aria cell sorter gating for DAPI-/DRAQ5+/GFP+ cells. For each sample a maximum of 1000 cells were isolated and directly collected in 100ul of RA1 buffer with 2ul tris(2-carboxyethyl)phosphine (TCEP) and low mass RNA was isolated with on-column DNase digestion using Nucleospin RNA XS kit (Clontech) following the manufacture’s protocol.

Supplementary Figure 1: Characterization of TPR50 tau Transgenic Mouse (a) Body weight (g=grams) of male wild type and transgenic mice at three and six months of age (error bar = SEM, unpaired two-tailed t-test, n=10/group). (b) Body weight (g=grams) of female wild type and transgenic mice at three and six months of age (error bar = SEM, unpaired two-tailed t-test, n=10/group). (c) Quantification of Western blots of soluble pS202/pT205 tau as detected by AT8 antibody (error bar = SEM, n=6/group). (d) Quantification of Western blots of soluble total tau as detected by Aβ−3 tau antibody (error bar = SEM, n=6/group). (e) Quantification of Western blots of Sarkosyl-insoluble pS202/pT205 tau as detected by AT8 antibody (error bar = SEM, unpaired two-tailed t-test, n=6/group). (f) Quantification of Western blots of Sarkosyl-insoluble total tau as detected by Aβ−3 tau antibody (error bar = SEM, unpaired two-tailed t-test, n=6/group). (g) Phospho-tau levels as measured by ELISA using pT231 tau antibody at three months (error bar = SEM, n=6/group). (h) Phospho-tau levels as measured by ELISA using pT231 tau antibody at six months (error bar = SEM, unpaired two-tailed t-test, n=6/group) (i) Brain weight (mg) of male wild type and transgenic mice at six months of age. (error bar = SEM, unpaired two-tailed t-test, n=6/group). In all the plots. the center line represents the mean and error bars showing SEM.

Supplementary Figure 2: Transcriptome-wide Differential Gene Expression in Different Brain Regions and Mouse Strains (a-c) 2-D scatterplots of the first three principal components of differentially expressed (DE) genes with FDR <0.1 (p values obtained using linear regression) between wild type and transgenic in the different brain regions indicated. (d) Scatterplot showing the fold change using qPCR and RNA-seq for C57 TPR50 transgenic mice compared to wild-type controls (n=6/group, Pearson correlation). (e) Scatterplot between p-values obtained using DESeq2 and linear regression model between wild-type and transgenic mice (n=6/group, Pearson correlation).

Supplementary Figure 3: Differential Expression and Regulators of Apoptosis in Different Brain Regions and Mouse Strains. (a) Rank-rank hypergeometric overlap (RRHO) plots comparing wild type and transgenic gene expression in the cortex to the hippocampus, brain stem and cerebellum. Each pixel represents the overlap between brain regions, color-coded according to the −log10 p-value of a hypergeometric test. Scatterplot of correlations between DE genes with p<0.005 in the cortex (n=36 mice) as compared to the hippocampus (left panel, n=36 mice), brain stem (left panel, n=36 mice) and cerebellum (right panel, n=36 mice). (b-i) Expression of genes involved in positive regulation of apoptosis (GO: 0043065) in transgenic and wild type mice in different brain regions, at different time points across strains (n=6/group, unpaired two-tailed Wilcox rank test). (i-l) Expression of genes involved in negative regulation of apoptosis (GO: 0043066) in transgenic and wild type mice in different brain regions, at different time points across strains (n=6/group, unpaired two-tailed Wilcox rank test). (j) Venn diagram showing DE gene overlap (genes with p<0.005) between brain regions (n=6/group). (k) Barplot of the number of DE genes with p<0.005 (n=6/group). DE analysis was conducted independently for each brain region, mouse strain and time point. In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median. (l) Log2 fold change of the top twenty most upregulated genes in the cerebellum of TPR50 transgenic mice (all strains) and their respective expression in the cortex (n=6/group). All show opposite patterns in cerebellum compared with cortex, making them potential protective changes.

Supplementary Figure 4: mRNA Consensus Co-expression Network Analysis. (a) cWGCNA dendrogram showing the mRNA modules in the cortical samples of transgenic mice at 6 months. Color bars below give information on correlation of transgenic condition of RNA integrity number (RIN) with expression of a particular gene, red is positive correlation and blue is anti-correlation (n=6/group). (b) Cell-type enrichment using various cell-type specific mouse and human datasets in mRNA modules (n=15 modules, fisher’s exact two-tailed test, **FDR<0.01, ***FDR<0.005). (c-d) KEGG term enrichment in the NAS (c) and NAI (d) modules (n=15 modules). (e) Enrichment of genes correlated with phospho-tau levels (pT231) in transgenic samples of all strains at 6 months (n=18/strain) (fisher’s two-sided exact test, ***FDR<0.005). (f-g) Trajectory of the NAS (f) and NAI (g) module eigengenes in the hippocampus of TPR50 mice (Unpaired two-sided Wilcox rank test, n=6 mice/group). (h-i) Trajectory of the NAS (h) and NAI (i) module eigengenes in the brain-stem of TPR50 mice (Unpaired two-sided Wilcox rank test, n=6 mice/group). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median.

Supplementary Figure 5: NAS and NAI Module Expression in published gene-expression datasets. (a) NAI module eigengene (ME) expressions from sorted microglia from Tg4510 mouse (unpaired two-sided Wilcox rank test, n=4 mice/group). (b-d) NAS ME expression from sorted neurons (b) and NAI ME from sorted microglia (c) and astrocytes (d) from PS2APP mice (Unpaired two-sided Wilcox rank test, n=6 mice/group). (e) Summary heatmap of module preservation of all datasets described using module definitions from transgenic mouse network analysis. (f-l) NAS and NAI module eigengene (ME) expressions in Tg4510 replication RNA-seq set (f, n=5/group), (g) Tg4510 (Wes et al., 2014, n=4/group; tTA-control is Camk2a-tTA control for Tg4510 bigenic mice), (h) Tau-P301L (Matarin et al., 2015, n=6/group), (I) CRND8 APP (Matarin et al., 2015 n=5/group), (j) APP/PS1 (Matarin et al., 2015, n=6/group), (k) Progranulin knockout (Lui et al., 2016, n=6/group); Unpaired two-sided Wilcox rank test. (l) NAS ME for PS19 mouse hippocampus.(unpaired two-sided Wilcox rank test, n=6/group). (m-n) NAS and NAI module eigengene (ME) expression in human cortical (m) and cerebellar (n) control samples (n=8); Progranulin-positive FTD cases (FTD GRN Pos., n=6) and progranulin-negative FTD cases (FTD GRN Neg., n=10) are shown (unpaired two-sided Wilcox rank test). (o) NAS module eigengene expression in human iPSC-derived neurons from FTD patients and controls; Progranulin-positive FTD (FTD GRN Pos.) and progranulin-negative FTD (FTD GRN Neg.) are shown, unpaired two-sided Wilcox rank test, n=5/group. In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median. (p) Description of workflow for RNA and protein sample processing in human FTD samples.

Supplementary Figure 6: miRNA Co-expression Network Analysis. (a-b) Trajectory of the miM12 module eigengene in the cortex (a) and cerebellum (b) (unpaired two-sided Wilcox rank test, Bonferroni corrected p-values, n=6 mice/group). (c-e) Trajectory of the miM16 module eigengene in the hippocampus (c), brain stem (d), and cerebellum (e) (unpaired two-sided Wilcox rank test, Bonferroni corrected p-values, n=6 mice/group). (f) Correlation of mRNA and miRNA module eigengenes, Pearson correlation test p-value (g) Enrichment of all miM16 module miRNA predicted targets in mRNA modules (fisher’s two-sided exact test, ***FDR<0.005). Targetscan database was used for miRNA target prediction. (h) Conservation of miR-203 across species. (i-k) Trajectory miR-203 expression in the hippocampus (i), brain stem (j) and cerebellum (k) (unpaired two-sided Wilcox rank test, Bonferroni corrected p-values, n=6 mice/group). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median.

Supplementary Figure 7: miR-203 overexpression (a) miR-203 expression levels in infected primary cortical neurons (error bar = SEM, n=4 DIV6, n=8 DIV8 and DIV10 samples/condition). The center line represents the mean. (b) Trajectory of module preservation in infected primary cortical neurons using module definitions from transgenic mouse network analysis (n=4 DIV6, n=8 DIV8 and DIV10 samples/condition). (c) Representative images from three independent experiments of GFP fluorescence in coronal slices of stereotaxically injected C57BL/6 wild type 3 and 6-weeks after AAV2/9 injection (n=6/group), scale bar = 500μm. (d-e) Representative from four independent experiments flow cytometry dot plots of sc-miRNA (d) and miR-203 (e) infected neurons showing the various gatings used. DRAQ5 is used as a live cell marker and DAPI as a nuclear dead cell marker (n=6/group). (f) miR-203 expression levels in GFP-positive FACS sorted neurons from C57BL/6 wild type 3 and 6-weeks after AAV2/9 injection (unpaired two-sided Wilcox rank test, n=6/group). (g) Module preservation in GFP-positive neurons overexpressing sc-miRNA (left) or miR-203 (right) using module definitions from TPR50 mouse network analysis (n=6/group). (h-i) Correlation between the module eigengene of predicted miR-203 targets and miR-203 expression (h) and between the NAS module eigengene and miR-203 expression (i) in GFP-positive neurons from C57BL/6 wild type mice overexpressing miR-203 (n=6/group, Pearson correlation). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median. (j) Representative images from three independent experiments of Caspase-8 staining in C57BL/6 wild type mouse frontal cortex infected to overexpress miR-203 or sc-miRNA (control) for three or six weeks (n=6/group). Green: GFP (infection), Red: Casp8 and Blue: Hoechst. Scale bar=20μm.

Supplementary Figure 8: Modulation of miR-203 expression using AAV system in Tg4510 Tau transgenic mice. (a) Representative images from three independent experiments of GFP fluorescence in coronal slices of stereotaxically injected Tg4510 tau transgenic mice 6-weeks after AAV2/9 injection (n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice), scale bar = 500μm. (b) miR-203 expression levels in GFP-positive FACS sorted neurons from one-month Tg4510 mice 6-weeks after AAV2/9 injection (unpaired two-sided Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). (c-d) Expression of genes involved in positive regulation of apoptosis (GO: 0043065; c) or negative regulation of apoptosis (GO: 0043066; d) in GFP-positive cells overexpressing miR-203 6 weeks after AAV injection (Unpaired two-sided Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). (e-f) Expression of genes involved in calcium signaling (KEGG ID: mmu04020; e) and neuroactive receptor ligand mediated pathway (KEGG ID: mmu04080; f) in GFP-positive cells overexpressing miR-203 6 weeks after AAV injection in Tg4510 tau transgenic mice (unpaired two-sided Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). (g,i) Representative images from three independent experiments of GFP fluorescence in coronal slices of stereotaxically injected Tg4510 tau transgenic mice 6-weeks (g) or 4-months (i) after AAV2/9 injection with miR-203 TuD, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice; scale bar = 500μm. (h,j) miR-203 expression levels in GFP-positive FACS sorted neurons from five-month old Tg4510 mice 6-weeks (h) or 4-months (j) after AAV2/9 injection (Unpaired two-sided Wilcox rank test, n=5 scrambled (sc)-miRNA infected mice and n=6 miR203 infected mice). (k-p) Expression of targets of miR-203 (k), NAS ME (l), genes involved in calcium signaling (KEGG ID: mmu04020; m) and neuroactive receptor ligand mediated pathway (KEGG ID: mmu04080; m), genes involved in positive regulation of apoptosis (GO: 0043065; o) or negative regulation of apoptosis (GO: 0043066; p) in GFP-positive cells expressing sc-TuD control or miR-203 TuD in five-month old Tg4510 mice for 6-weeks (unpaired two-sided Wilcox rank test, n=6/group). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median. (q) Representative images from three independent experiments of hyperphosphorylated tau AT-8 staining in Tg4510 tau transgenic mice frontal cortex infected to express TuD-miR203 for 4 months. Green: GFP (infection), Red: AT8 and Blue: Hoechst. Scale bar=20μm, (n=6/group). (r) AT8 Tau intensity of GFP-positive cells expressing TuD-miR203 in Tg4510 Tau transgenic mice were normalized at each time point to the average AT8 intensity of GFP-positive cells expressing sc-TuD control (n=55 images from 5 mice for sc-TuD control and n=106 images from 6 mice for miR-203 TuD). The center line represents the mean and error bars showing SEM. (s) Representative images from three independent experiments of c-fos staining in Tg4510 tau transgenic mice frontal cortex infected to inhibit miR-203 for 4 months (n=6/group). Green: GFP (infection), Red: c-fos and Blue: Hoechst. Scale bar=20μm.

Supplementary Figure 9: Uncropped Western Blots. Uncropped images of the four Western blots shown in Figure 5g. The upper panels are an overlay of the color images of the blots, showing the ladder indicating molecular weight markers, and the luminescent images, showing the protein band detected. The lower panels are the luminescent images showing the protein bands. All images show the entire section of the blot. The red squares indicate the cropped sections depicted in Fig5g. Full ladder shown separately. Representative images from three independent experiments.

Supplementary Figure 10: Strain-aware Network Analysis. (a) Summary heatmap of module preservation of all datasets described in Supplementary Fig. 5e using module definitions from strain-aware network analysis. (b-c) Trajectory of strain-aware network – M1 module (b) and M5 module (c) at six-months in cortex (unpaired two-tailed Wilcox rank test, Bonferroni corrected p-values, n=6/group). In all the boxplots, the upper and lower lines represent the 75th and 25th percentiles, respectively, while the center line represents the median.

Supplementary Table 1: TPR50 mRNA-Seq Analysis (a-d) Differential gene expression and co-expression network module assignments from mRNA-seq for cortex (a), hippocampus (b), brain stem (c) and cerebellum (d), n=6/group. Linear regression was used to calculate p-values and FDR calculated p-values are also reported. (e-h) Module eigengenes for all samples and sample characteristics for TPR50 data in cortex (e), hippocampus (f), brain stem (g) and cerebellum (h). (i) Significant genes (with p<0.05) risk loci in GWAS data as obtained using MAGMA for AD, PSP and FTD GWAS datasets. (j) Metadata of TPR50 dataset along with sequencing statistics obtained from Picard Tools.

Supplementary Table 2: Human transcriptomics and proteomics data. (a) Sample size and accession numbers of different human datasets used in the manuscript. (b) Metadata for postmortem human brain samples used for RNA-seq from cortex and cerebellum along with biological and technical information. (c) Metadata for postmortem human brain samples used for proteomics. (d) Protein groups quantified by Maxquant and used for downstream analyses. (e) Overlap of proteins identified in human FTD proteomics dataset and Human FTD transcriptomics dataset.

Supplementary Table 3: TPR50 miRNA-seq Analysis. (a-d) Differential gene expression and co-expression network module assignments from miRNA-seq for cortex (a), hippocampus (b), brain stem (c) and cerebellum (d). (e-h) Module eigengenes for all samples and sample characteristics for TPR50 data in cortex (e), hippocampus (f), brain stem (g) and cerebellum (h).

Supplementary Table 4: Module preservation and connectivity Map. (a) Module preservation across all modules are shown. For human FTD mRNA-seq (Control n=14, Tau-positive FTD n=10, Tau-negative FTD n=6), primary cortical neurons (n=6 for DIV6, DIV8 and DIV10 each), FACS AAV GFP and miR203 (n=6/group). P-values are calculated by transforming the Z summary statistics. (b) List of genes that are positive and negative regulators of apoptosis that are also in NAS and NAI module. (c) Connectivity Map Output. (d) IGAP AD GWAS enrichment analysis without APOE gene across all modules , Two-sided p-values are calculated using a permutation test.