Dissecting the role of non-coding RNAs in the accumulation of amyloid and tau neuropathologies in Alzheimer's disease.

BACKGROUND: Given multiple studies of brain microRNA (miRNA) in relation to Alzheimer's disease (AD) with few consistent results and the heterogeneity of this disease, the objective of this study was to explore their mechanism by evaluating their relation to different elements of Alzheimer's disease pathology, confounding factors and mRNA expression data from the same subjects in the same brain region.
METHODS: We report analyses of expression profiling of miRNA (n = 700 subjects) and lincRNA (n = 540 subjects) from the dorsolateral prefrontal cortex of individuals participating in two longitudinal cohort studies of aging.
RESULTS: We confirm the association of two well-established miRNA (miR-132, miR-129) with pathologic AD in our dataset and then further characterize this association in terms of its component neuritic β-amyloid plaques and neurofibrillary tangle pathologies. Additionally, we identify one new miRNA (miR-99) and four lincRNA that are associated with these traits. Many other previously reported associations of microRNA with AD are associated with the confounders quantified in our longitudinal cohort. Finally, by performing analyses integrating both miRNA and RNA sequence data from the same individuals (525 samples), we characterize the impact of AD associated miRNA on human brain expression: we show that the effects of miR-132 and miR-129-5b converge on certain genes such as EP300 and find a role for miR200 and its target genes in AD using an integrated miRNA/mRNA analysis.
CONCLUSIONS: Overall, miRNAs play a modest role in human AD, but we observe robust evidence that a small number of miRNAs are responsible for specific alterations in the cortical transcriptome that are associated with AD.

Background

Late onset Alzheimer’s disease (AD) is an age-dependent neurodegenerative disorder characterized clinically by cognitive decline and pathologically by the accumulation of neuritic β-amyloid plaques (NP) and neurofibrillary tangles (NFT) in the brain. Currently genetic [1], epigenomic [2] and transcriptomic studies [3, 4] coupled with advances in imaging techniques [5, 6] have begun to sketch the sequence of events in the causal chain linking risk factors to a syndromic diagnosis of AD dementia. One of these events may be the dysregulation of gene expression by alterations in the expression of microRNA (miRNA) and long intergenic non-coding RNA (lincRNA) molecules.

miRNA are a class of small regulatory RNA that modulate gene expression via a multiprotein complex which facilitates the interaction between an miRNA and its complementary elements in the 3’UTR of target mRNAs to initiate transcript degradation and repression of protein production [7, 8]. Aberrant expression of miRNA and/or its target mRNAs have been implicated in abnormal neuron function [9] and in several neurodegenerative disorders [10, 11]. Recently, certain miRNA, such as miR-132 have been associated with pathologic AD [12-16]. However, these studies were conducted in a modest number of subjects with limited phenotypic information, and few results are consistent across these studies [17].

lincRNA are RNA that are longer than 200 nucleotides and do not code for proteins [18]. As with most long non-coding RNA and unlike miRNA there is no clear common functional mechanism for lincRNA [18]; some may have a structural role in protein/RNA complexes. There is still debate over what percentage of lincRNA may be functional at all [19]. However, focusing on the lincRNA that lie in the same locus as protein coding genes may provide insight into their functional correlates. Given the reported association of the long non-coding RNA BACE1-AS [20] with AD and the lack of investigation of this class of non-coding RNAs in a large sample size, we investigated the potential role of long non-coding RNA in AD and its component pathologies.

We first evaluated the role of miRNAs previously associated with AD with measures of both neuritic amyloid plaque (NP) and neurofibrillary tangles (NFT) since these are key neuropathologic features of AD and allow us to explore the mechanism by which an AD-associated non-coding RNA contributes to disease. We secondarily expanded this effort to evaluate other miRNAs and lincRNAs to discover new associations. In addition, we leveraged transcriptome-wide RNA sequencing profiles, generated from the same RNA samples that were used to generate the miRNA profiles, to identify the functional consequences of altered miRNA expression on protein-coding genes in the human dorsolateral prefrontal cortex (DLPFC) and to identify additional cases where miRNA and mRNA AD associations converge in the human cortex.

Methods

Total RNA, including miRNA, was extracted from approximately 100 mg sections of frozen postmortem brain tissue from the dorsolateral prefrontal cortex (BA 9/46) of subjects from two previously described longitudinal cohorts of aging, Religious Order Study (ROS) and Rush Memory Aging Project (MAP) [2, 21–23]. Tissue was thawed partially on ice and between 50 and 100 mg of gray matter was dissected from the section then transferred immediately to 1 mL of Trizol. The tissue was then quickly homogenized using the Qiagen TissueLyser and a 5 mm stainless steel bead, for 30 s at 30 Hz. The foam was settled with a quick spin, and the sample incubated for a minimum of 5 mins at room temperature. Debris was pelleted at 4 °C at 12,000 g for 10 min and Trizol was transferred to a new 1.5 mL tube. We continued preparation of samples following the instructions of Qiagen’s miRNeasy Mini kit, with volumes adjusted for 1 mL instead of 700uL of Trizol until wash steps. RNA was eluted from the miRNeasy spin columns in 75uL of elution buffer, and quality tested by Nanodrop and Bioanalyzer RNA 6000 Nano Agilent chips. RNA yields averaged about 25μg and RIN scores ranged from 2 to 9 with an average of 6.5. RNA was normalized to 33 ng/ul and plated into 96w plates for Nanostring processing using the nCounter Human miRNA Expression assay kit version 1 with reporter library file: NS_H_miR_1.2.rlf. The data collection from 733 post mortem brain samples was done at the Broad Institutes Genomics Platform (Broad Institute of Harvard and MIT, Cambridge, USA). Subjects from different diagnostic categories were distributed across experimental batches to reduce batch effects. To minimize variability at the ligation step, processing of the annealing and ligation steps was performed on the same thermocycler. Two thermocyclers were used for the purification steps, but all samples were placed in the same thermocycler for denaturation and hybridization steps. Two nCounter Prep Stations were utilized, but all samples were then scanned on the same Digital Analyzer. The data was collected in 8 batches of 96 samples and a single sample technical replicate was introduced as control in every single 96 well plates and sometimes twice in one single plate in two different cartridges.

Quality control and dataset pre-processing

All data is available at https://www.synapse.org/#!Synapse:syn3219045. The miRNA from the Nanostring RCC files were re-annotated to match the definitions from the miRBase v19. The raw data from the Nanostring RCC files were accumulated and the probe-specific backgrounds were adjusted according to the Nanostring guidelines with the corrections provided with the probe sets. After correcting for the probe-specific backgrounds, a three-step filtering of miRNA and sample expressions was performed. First, miRNA that had less than 95% of samples with a missing expression level were removed. This is followed by removing samples that had less than 95% of miRNAs with missing expressions. Thus, the call-rates for the samples and the miRNA are set at 95%. Finally, all miRNA whose absolute value is less than 15 in at least 50% of the samples were removed to eliminate miRNA that had negligible expression in brain samples. After the miRNA and sample filtering, the dataset consisted of 309 miRNAs and 700 samples. A combination of quantile normalization and Combat [24], specifying the cartridges as batches for the miRNA data, was used to normalize the data sets. The strong association observed between miRNA expression and RNA-integrity was validated via qRT-PCR (Additional file 1: Figure S3) and verified not to be specific to the Nanostring platform.

The mRNA sequencing data have been described elsewhere [25]. For the lincRNA analysis, a non-gapped aligner Bowtie [26] was used to align reads to a hg19 lincRNA reference [27] and then RSEM [28] applied to estimate expression levels for all lincRNA. A combination of quantile normalization and Combat [24] to account for batch was used for normalization. After filtering out lowly expressed lincRNA with an expected count less than 5, a final dataset of 454 lincRNA measured in 540 samples was used.

Identifying differentially expressed miRNA or lincRNA

Simple linear regression analysis was used to associate the expression levels of miRNA and lincRNA to several variables that measured the neuropathology in the AD brain. These included either the numbers of neuritic plaques, neurofibrillary tangles or a binary variable representing the pathologic diagnosis of AD on autopsy according to the NIA Reagan criteria [29]. All the associations were adjusted for age, sex, study (ROS or MAP), the proportion of neurons in the tissue [2], RNA Integrity number (RIN) and post-mortem interval (PMI). The proportion of neurons in the tissue is estimated from DNA methylation data available from the same brain region of each individual, as described in our recent study [2]. A Bonferroni corrected significance threshold of 0.05 was used to account for multiple comparisons.

Constructing micro-RNA and linc-RNA networks

Linear models were used to identify miRNA or lincRNA that were associated with either NP, NFT or AD. To do this the models included age, neuronal composition (NNLS), sex, study of origin (Study, ROS or MAP), post-mortem interval (PMI), RNA integrity number (RIN), NP, NFT and AD., NP, NFT and AD as covariates. Using these models a miRNA or lincRNA were included in the networks if there was evidence that any of effects sizes of NP, NFT or AD were non-zero (nominal p-value from an F-test from less than 0.05). For each of the included miRNA or lincRNA, forward stepwise variable selection was used with a Bayesian information criteria (BIC) to select which edges between miRNA or lincRNA and explanatory variables should be included in the network. As RIN is associated with all the miRNA and lincRNA its edges are excluded from the networks.

Pathway analysis

Our approach for pathway analysis of the miRNA data using pMim [30] involved analyses of the miRNA and mRNA data on the 525 samples which had both miRNA and gene expression data. The mRNA analyses are designed to focus on sets of genes that are co-expressed and are predicted to be a target of one of the tested miRNAs. Targetscan v6.2 [31] was used for prediction of miRNA targets and the GO biological processes [32] were used for pathway annotation. Figure 3a outlines how pMim is used to construct and test “miR-pathways”, where a miR-pathway consists of genes that are targeted by a miRNA and lie in a common pathway. Specifically, pMim identifies sets of genes that (1) are predicted to be targeted by an miRNA associated with AD, (2) lie in a common biological pathway and (3) are also associated with AD diagnosis themselves in terms of mRNA expression. All 309 miRNA and their corresponding miR-pathways (sets of genes that are targeted by the same miRNA and share a common biological process) were tested; a joint statistic summarizes and ranks the evidence for both the miRNA and mRNA analyses testing if both a miRNA and one of its corresponding miR-pathways are associated with AD. This joint statistic is calculated with two significance combination methods. Within a miR-pathway, Stouffer’s method was used to combine significance of genes associated with AD. A one-sided Pearson’s method was used to combine the significance between the miR-pathway gene summary statistic and the association of its miRNA with AD. Under very strict assumptions these joint statistics could be considered as p-values, however due to the large amount of correlation within an annotated pathway they are highly inflated. Hence we consider these joint statistics for ranking purposes only.

Fig. 3

An integrated analysis of miRNA and their target mRNA. a This diagram presents an overview of (1) how the miR-132:Protein acetylation miR-pathway was constructed from the KEGG pathways database and Targetscan miRNA binding target prediction database and (2) how a score was derived to assess the association of this miR-pathway combination with pathological AD diagnosis. The approach is described in greater detail elsewhere [30]. b This heatmap illustrates the relation of miR-132 with the outcome measures (AD, NP, and NFT) and the mRNA expression levels of five of its target genes that are member of a protein acetylation miR-pathway. Each column presents data from one subject. Scaled expression values for a given variable are colored from red (to denote high expression) to blue (low expression). c The relationship between miR-132 and one of its targets, EP300, is demonstrated in a scatterplot. Each dot represents one subject. d The standardized effect sizes () capturing the strength of the relationship between the five protein acetylation genes and AD are plotted in this panel. We highlight the relative amount of the effect size (β) which can and cannot be explained by miR-132 expression, as miR-132 is predicted to be regulating each of these genes

Results

Demographic features and nature of our non-coding RNA data

The NanoString nCounter assay was used to measure miRNA expression from the DLPFC of each subject. At the conclusion of a rigorous quality control and preprocessing pipeline, expression measures from 309 miRNAs in 700 subjects were retained for downstream analysis. The subjects profiled in this study are participants in one of two prospective cohort studies of aging, the Religious Order Study (ROS) and Memory and Aging Project (MAP) which are designed to be merged for joint analyses [33, 34]. These studies enroll non-demented individuals and include detailed, annual ante-mortem characterization of each subject’s cognitive status as well as prospective brain collection and a structured neuropathologic examination at the time of death. The study design of ROS and MAP yields an autopsy sample that includes a range of syndromic diagnoses and neuropathologic findings that are common in the older population. While this is not a true population-based sample, it captures the diversity of the older population at the time of death and is different than the collection of subjects typically used for age-matched case-control studies [12, 13, 15, 16]. The subjects in the study have an average age of 88 at death, 61% meet criteria for pathologic AD by NIA Reagan criteria [29] and 64% are female (Additional file 1: Table S1). To evaluate the nature of our data in relation to published results, we first assessed whether the expression of miRNAs are associated with a pathological diagnosis of AD according to the NIA Reagan criteria [29], focusing on an evaluation of previously reported miRNA associations [12, 13, 15, 16].

miRNA associated with a pathologic diagnosis of AD

We used a linear regression model controlling for age, sex, study, proportion of neuronal cells, post-mortem interval and a measure of RNA quality (RNA integrity (RIN) score) to evaluate associations of miRNA with pathologic AD. In this analysis of the human DLPFC, two previously reported associated miRNAs are significantly associated with pathologic AD, exceeding a threshold defined using a Bonferroni correction for the number of miRNA tested, p < 1.6 × 10−4 (Table 1, Fig. 1). Specifically, miR-132 (P = 2 × 10−8) and miR-129-5p (P = 3.9 × 10−6) were found to be diminished in expression in subjects with AD, consistent with prior studies [12, 13]. Thus, our data confirm well-validated miRNA associations with AD. To provide perspective on effect size, miR-132 explains 6.7% of the variance in AD, which is similar to the 6.1% variance explained by APOE ε4 in our data and 6% in previous literature [35]. In contrast, the miRNA effect is more than the <1% of variance explained by genetic susceptibility variants in ABCA7, CD33, and CR1 that were originally discovered in genome-wide association studies [35].

Table 1

miRNA and lincRNA associated with AD, NP or NFT

	coefficient: AD	p-value: AD	coefficient: NP	p-value: NP	coefficient: NFT	p-value: NFT
miR-132	−0.26	2.00E-08	−0.29	1.40E-11	−0.31	4.10E-08
miR-129-5p	−0.11	3.90E-06	−0.15	6.90E-11	−0.17	2.20E-08
miR-129-3p	−0.095	0.0013	−0.13	1.80E-06	−0.12	0.00085
miR-99b	0.081	0.0034	0.075	0.0037	0.14	6.00E-05
linc-CTSD-3	−19	0.0012	−22	6.90E-05	−33	8.60E-06
linc-BRD9–1	32	0.00025	35	2.40E-05	23	0.036
linc-ADC	3.8	0.18	3.8	0.15	14	6.70E-05
linc-RNFT2–1	3.3	0.36	3.4	0.31	18	7.10E-05

In bold are the tests which meet the Bonferroni correction threshold of p < 0.00016

Fig. 1

Validation of miRNAs associated with AD in other studies. The relationship between average miRNA expression (x-axis) and association with pathological AD in ROSMAP (y-axis) is plotted for each miRNA. The association between a miRNA and AD is calculated as a z-score after correcting for other factors; a negative z-score denotes decreased miRNA expression in the context of AD. Previously reported miRNAs [12] are highlighted by colored triangles, with red denoting higher expression in AD and blue reduced expression in AD in the original publication. The dotted horizontal line marks the z-score threshold of significance that is equivalent to a Bonferroni p-value cut-off of 0.05; none of the miRNAs with positively correlated expression are significantly associated with AD. The gray circles represent all other miRNAs tested in our study

We next systematically reviewed our results for those miRNA that have previously been associated with AD in various brain regions [13, 15, 16]: they are catalogued in Table 2. Figure 1 illustrates the strength of the association of miRNAs with pathologic AD in our cohort, highlighting 11 miRNAs that were deemed to be significantly associated with AD in a prior study of moderate sample size (n = 49 subjects) [12]. Of the 11 miRNA previously reported to be associated with AD in the prefrontal cortex, only miR-132 and miR-129-5p replicated. In addition, miR-129-3p (P = 0.0013), miR-200a (P = 0.018), miR-30e (P = 0.0048), miR-100 (P = 0.011), let-7i (P = 0.045) and miR-185 (P = 0.013) display suggestive association and may warrant further investigation. The Bonferroni threshold (p < 0.00016) that we employed may be overly conservative, but it provides a useful way to distinguish the most robustly associated miRNA. Given our large sample size (n = 700), we have good power to confirm results with small effect sizes: we have 80% power to detect an effect small enough to only explain 2.2% of the variance in relation to AD (which is much smaller than the effect of mi132 described above). Thus, if other miRNAs that we tested have a role in AD, they are more in line with the magnitude of effect found with AD SNPs.

Table 2

P-values from tests of association with AD and other covariates for miRNA implicated in other studies

	ROSMAP	Lau	Cogswell	Hebert	Nunez	Study	Age	Sex	NNLS	PMI	RIN	AD	AD (no covariates)
miR-132	Down	Down	Down			0.54	0.58	0.32	0.46	0.0052	7.00E-14	2.00E-08	4.00E-12
miR-129-5p	Down	Down				0.31	0.012	0.72	0.23	6.00E-04	1.70E-11	3.90E-06	6.50E-07
miR-129-3p		Down				0.59	0.089	0.063	0.73	0.0057	3.00E-08	0.0013	0.00098
miR-136		Down				0.59	0.3	0.71	0.16	0.046	0.1	0.67	0.78
miR-370		Down				0.29	0.084	0.025	0.051	0.17	0.0036	0.67	0.57
miR-409-5p		Down				0.52	0.0064	0.22	0.26	0.52	0.071	0.33	0.53
miR-487a		Down				0.89	0.1	0.78	0.79	0.53	5.90E-09	0.8	0.99
miR-92a		Up	Up			0.36	0.81	0.039	0.00056	0.43	7.50E-10	0.95	0.1
miR-27a		Up	Up			0.39	0.69	0.46	0.47	0.53	0.56	0.56	0.6
miR-92b		Up	Up			0.14	0.44	2.70E-05	0.96	0.74	2.30E-18	0.056	0.74
miR-200a		Up				0.95	0.013	0.58	0.65	0.095	2.00E-53	0.018	0.00013
miR-9			Down	Down		0.47	0.15	0.21	0.005	0.65	1.20E-06	0.73	0.39
miR-146b-5p			Down			0.45	0.12	0.92	0.5	0.36	0.0064	0.93	0.84
miR-425			Down			0.35	0.19	0.88	0.34	0.41	0.35	0.29	0.11
miR-30e			Up		Down	0.91	0.043	0.14	0.71	0.33	0.019	0.0048	0.0043
miR-423-5p			Up			0.5	0.023	0.78	0.42	0.076	1.80E-07	0.38	0.25
miR-27b			Up			0.19	0.57	0.84	0.66	0.16	0.08	0.48	0.38
miR-100			Up			0.13	0.99	0.0091	0.11	0.16	8.70E-05	0.011	0.0083
miR-34a			Up			0.4	0.00015	0.056	0.013	0.39	0.59	0.18	0.028
miR-145			Up			0.015	0.63	0.96	0.89	0.55	5.50E-06	0.26	0.084
miR-381			Up			0.55	0.73	0.12	0.74	0.6	0.47	0.12	0.022
miR-125b			Up			0.28	0.42	0.92	0.66	0.6	3.20E-10	0.63	0.37
miR-148a			Up			0.19	0.74	0.22	0.23	0.97	0.016	0.95	0.69
miR-15a				Down	Down	0.47	0.24	0.49	0.016	0.0053	0.031	0.51	0.48
miR-29b				Down	Down	0.98	0.27	0.15	0.99	0.057	0.00017	0.34	0.019
miR-101				Down	Down	0.77	0.24	0.013	0.21	0.099	0.96	0.3	0.15
miR-181c				Down	Down	0.13	0.49	0.24	0.038	0.23	0.42	0.34	0.3
miR-363				Down		0.2	0.027	0.13	0.2	0.072	0.096	0.92	0.54
miR-19b				Down		0.72	0.05	0.069	0.34	0.18	0.65	0.064	0.032
miR-106b				Down		0.25	0.92	0.95	0.68	0.3	0.82	0.15	0.13
miR-22				Down		0.84	0.35	0.75	0.16	0.64	0.15	0.81	0.77
miR-93				Down		0.11	0.7	0.0064	0.81	0.74	0.03	0.34	0.2
miR-26b				Down		0.34	0.36	0.25	0.94	0.76	0.67	0.13	0.21
let-7i				Down		0.46	0.79	0.79	0.58	0.93	0.39	0.045	0.036
miR-320a				Up	Up	0.44	0.24	0.26	0.78	0.38	4.60E-14	0.16	0.057
miR-197				Up		0.23	0.43	0.49	0.36	0.31	0.39	0.77	0.43
miR-29c					Down	0.87	0.36	0.43	0.089	0.0082	2.00E-06	0.35	0.31
miR-494					Down	0.61	0.4	0.94	0.5	0.044	2.50E-11	0.38	0.16
miR-598					Down	0.26	1	0.25	9.70E-05	0.1	0.0034	0.2	0.11
miR-374a					Down	0.84	0.48	0.74	0.081	0.35	0.039	0.63	0.72
miR-376a					Down	0.93	0.58	0.027	0.16	0.58	0.71	0.11	0.012
miR-148b					Down	0.89	0.15	0.69	0.12	0.6	0.92	0.72	0.75
miR-95					Down	0.51	0.76	0.33	0.0017	0.71	0.056	0.86	0.74
miR-582-5p					Down	0.6	0.86	0.94	0.012	0.89	0.00018	0.35	0.18
miR-432					Up	0.59	0.32	0.99	0.16	0.084	0.53	0.12	0.1
miR-188-5p					Up	0.8	0.82	0.18	0.4	0.12	0.011	0.99	0.76
miR-382					Up	0.83	0.51	0.29	0.15	0.52	0.029	0.33	0.28
miR-185					Up	0.92	0.19	0.77	0.055	0.81	0.41	0.013	0.0044

In bold are the tests which meet the Bonferroni correction threshold of p < 0.00016

The role of confounding variables in miRNA expression

There is a substantial lack of concordance in behavior between the previously reported AD miRNA signatures and their level of association in our cohort (Table 2, Additional file 1: Figure S1). While this could be explained, in part, by differences in the brain regions that have been assayed and the technologies used in the different studies, many of these miRNAs are highly associated with age, sex, the proportion of neurons in the profiled tissue sample, post-mortem interval and/or RIN score (Table 2, Additional file 2: Table S2). Thus, explicitly accounting for all of these factors is essential in brain miRNA studies, and the RIN score appears to play an especially significant role.

Identifying miRNA associated with amyloid and tau neuropathologies

Having described the nature of our data in relation to prior reports, we proceeded to our main goal of investigating the relationship of miRNA expression with specific features of AD neuropathology: namely, we investigated whether miRNA are associated with the accumulation of neuritic amyloid plaques (NP) and/or neurofibrillary tangles (NFT) that are defined by the presence of Tau accumulation. The pathologic AD-associated miRNA, miR-132 and miR-129-5p, were associated with both NP and NFT (Table 1). These two neuropathologic features of AD are strongly correlated with one another and it can be challenging to dissect the effect of one from the other unless large sample sizes are available. In addition, we found that, among the tested miRNAs, miR-129-3p (P = 1.8 × 10−6) was associated with NP, and miR-99b (P = 6 × 10−5) was associated with NFT, exceeding our significance threshold (p < 0.00016).

Dissecting the network of miRNA associations in AD

The complexity of the relationships between miRNA expression, the burden of NP and NFT neuropathologies, a pathologic diagnosis of AD, technical variables and demographic variables are demonstrated in Fig. 2a. Here, we consider all variables in a single model to resolve the most likely primary association between each miRNA and the neuropathologic traits. These traits are all correlated with one another (Fig. 2b, Additional file 1: Table S3), making it difficult to understand where an miRNA’s primary effect may be exerted simply by looking at the univariate results: for example, while miR-132 is significantly associated with both NP and NFT in separate analyses (Table 1), it appears, in the joint network model, that its effect may be driven predominantly through NP. This is visualized by the edge linking miR-132 and NP in the network model. The effect on NFT may be secondary and result from the accumulation of NP. To elaborate a more comprehensive network model, we included all of the miRNAs with either a significant or suggestive (p < 0.05) association to one of the AD pathologies. Similar to miR-132, we see that many miRNAs are only associated with either NP, NFT or AD, suggesting that they may have a primary effect through one or the other pathology or some other aspect of AD not captured by these pathologic measures. On the other hand, miR-433, miR-487b and miR-485-3p are negatively associated with NP and also positively associated with AD in our joint model. This result captures the complexity and likely non-linear relationships of miRNA associations with pathology as miR-433, miR-487b and miR-485-3p are all negatively associated with AD in univariate analysis. This result highlights the need to evaluate the possibility of more complex relationships between miRNA and disease in future studies.

Fig. 2

Network of miRNAs associated with phenotypes. a The relationships between miRNA expression and various demographic and technical features as well as neuropathologic outcomes are displayed in a network diagram. All miRNA with a significant (p < 0.00016) or suggestive (p < 0.05) association are included in the network. These relationships were extracted via feature selection from linear models which explain an miRNA’s expression level by using age at death (Age), neuronal composition (NNLS), sex, study of origin (Study, ROS or MAP), post-mortem interval (PMI), RNA integrity number (RIN), neuritic amyloid plaques (NP), neurofibrillary tangles (NFT) and pathological AD diagnosis (AD). Each node in the diagram represents an miRNA (small circles) or a demographic or a neuropathologic variable (large circles). The significantly associated miRNA (Table 1) are represented by rectangles. A red edge represents a positive association between the miRNA and trait; a blue edge denotes an inverse association. As most of the miRNAs are associated with RIN, the associations with RIN were removed from this figure for clarity. b The correlations between the expression levels of eight miRNA and lincRNA with significant associations in our study are shown, along with correlations with the outcome variables (NP, NFT and AD). Each of the correlations are displayed in each cell of the correlation map and is colored by the strength of correlation, red for positive associations and blue for negative associations

Evaluating the functional consequences of AD-associated miRNA in human brain

We next attempted to delineate the downstream effects of the expressed miRNAs by integrating the miRNA data with gene expression data from the same brain region and the same individuals. In fact, both the miRNA and RNA sequence data were generated from the same RNA sample. Our approach involved separate analyses of miRNA and mRNA data, with the mRNA analyses being focused on sets of genes that are co-expressed and are predicted to be a target of one of the tested miRNAs. Specifically, pMim [30] was used to identify sets of genes that (1) are predicted to be targeted by an miRNA associated with pathologic AD, (2) lie in a common biological pathway and (3) are also associated with pathologic AD themselves in terms of mRNA expression. Figure 3a outlines how pMim is used to construct and test “miR-pathways” where Targetscan v6.2 [31] was used for prediction of miRNA targets and the GO biological processes [32] were used for pathway annotation. All 309 miRNA and their corresponding miR-pathways (sets of genes that are targeted by the same miRNA and share a common biological process) were tested; a joint statistic summarizes and ranks the evidence for both the miRNA and mRNA analyses. The top eight miR-pathways in which both the miRNA and the mRNA are significantly associated with AD are presented in Table 3.

Table 3

Integrated analysis of miRNA and their target mRNA identifies miR:target combinations that are associated with AD

miRNA	Pathways	Direction	Score	Genes in pathway targeted by miRNA
miR-132	protein deacetylation; protein deacylation; macromolecule deacylation	Down	1.17E-09	EP300, NACC2, SALL1, SREBF1, SIRT1
miR-200a	anion transport	Up	1.07E-07	SLC20A1, SLC25A3, RIMS1, SLC35A1, SLC35D1, GLS, ATP8A2, ATP8A1, IRS2, SLC38A9, CXCL12, SLC16A7, SLC18A2
miR-129-5p	regulation of transcription, DNA-templated	Down	1.19E-07	EP300, NCOA1, ZBTB20, PBX3, LMNA, ZBTB5, ETS1, CHD7, CBX4, WWC2, YWHAB, CBX7, PKNOX1, ZBTB10, CBX6, PHIP, ZFY
miR-129-5p	cell migration; cell motility; localization of cell	Down	1.19E-07	SUN2, DLC1, LMNA, TMEM201, ETS1, RPS6KB1
miR-129-5p	positive regulation of transcription, DNA-templated;positive regulation of transcription from RNA polymerase II promoter;positive regulation of RNA metabolic process;positive regulation of RNA biosynthetic process;positive regulation of nucleic acid-templated transcription	Down	1.19E-07	EP300, NCOA1, LMNA, ETS1, PKNOX1, PHIP
miR-129-5p	regulation of RNA metabolic process;regulation of nucleic acid-templated transcription;regulation of RNA biosynthetic process	Down	1.19E-07	EP300, NCOA1, ZBTB20, PBX3, LMNA, ZBTB5, ETS1, CHD7, CBX4, WWC2, YWHAB, ACTN1, CBX7, PKNOX1, ZBTB10, CBX6, PHIP, ZFY
miR-129-5p	cardiovascular system development; circulatory system development	Down	1.19E-07	EP300, DLC1, LMNA, ETS1, CHD7, PKNOX1, TIPARP
miR-129-5p	single-organism organelle organization	Down	1.19E-07	SUN2, EP300, DLC1, NCOA1, LMNA, VAMP1, WHAB, SNX9, ACTN1, PAPD5, PHIP

Looking more closely at the top miR-pathways, we see that, in the DLPFC, miR-132 appears to be targeting a group of five genes (EP300, NACC2, SALL1, SREBF1 and SIRT1) that are related to protein acetylation. All five of these protein acetylation genes lie within the top 20 predicted miR-132 target genes whose expression is inversely associated with miR-132 abundance in our cortical samples (p-values <0.00018, Table 4). The miR-132 association with EP300 expression confirms earlier reports and in vitro studies [12, 14, 36, 37], illustrating the robustness of our results and the relevance of earlier studies. As shown in Fig. 3b-d, miR-132 explains approximately a third of the association of each gene to AD, and therefore it does not appear to be the sole driver of the role of these protein acetylation genes in AD. While miR-132 may thus be influencing AD, in part, through histone acetylation and chromatin remodeling, miR-129-5p appears to be regulating genes related to the regulation of transcription. Interestingly, both miR-132 and miR-129-5p target EP300, which encodes the histone acetyltransferase protein E1A-associated cellular p300 transcriptional co-activator, and each miRNA explains some of the variation in EP300 expression (adjusted R-squared 0.039 and 0.033 respectively). The role of these two miRNAs appears to be non-redundant: when they are considered together, they explain more of the variance in EP300 than either miRNA explains separately (adjusted R-squared 0.045), further illustrating the complex interactions involved in mRNA regulation by miRNAs.

Table 4

Predicted target genes that are negatively associated with miR-132 expression in the human cortex

Gene	P-value: association with miR-132	P-value: association with AD	In protein acetylation pathway
CLMN	1.00E-14	0.0016
TJAP1	1.60E-12	0.00059
ANKRD40	9.30E-11	0.023
CXorf36	1.00E-10	0.0082
CDK19	2.60E-10	0.019
SSH2	6.10E-09	0.39
SEMA6A	4.90E-08	0.013
HIP1R	5.50E-08	0.0017
EP300	9.60E-08	0.0024	✓
DYNC1LI2	1.30E-07	0.023
NACC2	1.60E-07	0.021	✓
SOX5	1.60E-07	0.22
SIRT1	1.70E-07	0.069	✓
SALL1	5.10E-07	0.031	✓
MAPK3	4.10E-06	0.15
BCAN	2.00E-05	0.034
EPB41	0.00017	0.8
CTGF	0.00033	0.68
MEX3C	0.00038	0.53
SREBF1	0.00051	0.058	✓

In addition to these miR-pathways related to the validated miR-132 and miR-129-5p miRNAs, it is intriguing that miR-200a, which is suggestively associated to AD in our miRNA only analyses (Table 2), is the second strongest result in this pathway analysis because of the strong association of its downstream genes: this set of thirteen predicted target genes of miR-200a are all involved in anion transport. For these thirteen genes, an average of 18% of the mRNA association with AD is accounted by miR-200a. These integrated miRNA:mRNA analyses therefore report highly significant results drawn from autopsy tissue and set the stage for exploring specific downstream hypotheses; they also suggest that miRNAs are just one element contributing to the changes in gene expression found in AD since most of the variance of the affected genes – such EP300 and the miR200a targeted genes – remains unexplained at this time.

Deregulation of lincRNA expression in relation to AD neuropathology

Our RNA sequencing data was also used to measure lincRNA expression from the DLPFC in the ROSMAP cohorts. At the conclusion of a rigorous quality control and preprocessing pipeline, expression measures from 454 lincRNAs in 540 subjects were retained for downstream analysis. The subjects with lincRNA expression are primarily a subset of the subjects with miRNA expression; 525 subjects have both miRNA and lincRNA expression data.

To investigate the relationships between lincRNA expression and the available neuropathologic measures, we used the strategy deployed in studying miRNA: we fit three independent linear models using age, sex, study, an estimate of the proportion of neuronal cells, post-mortem interval and RNA integrity number and either NP, NFT or pathological AD as covariates to explain lincRNA expression. The lincRNA that were significantly associated with NP, NFT or AD are shown in Table 1 and Additional file 3: Table S4. No lincRNA is significantly associated with AD after Bonferroni correction. However, two are associated with NP (linc-CTSD-3 and linc-BRD9–1) and three with NFT (linc-CTSD-3, linc-ADC and linc-RNFT2–1); none of these lincRNAs have sequences that overlap the exons of protein coding genes.

We then repeated our network-building strategy to illustrate the primary associations among lincRNAs and the various demographic and pathological variables (Fig. 4a). We use linc-CTSD-3 to illustrate the fact that the pathologic AD-associated non-coding RNAs can also be independently influenced by confounding variables, highlighting the need to include them in the modeling (Fig. 4a, b, Additional file 1: Figure S2). As in the miRNA network, no lincRNA maintains a relationship between both NP and NFT after feature selection.

Fig. 4

Network of lincRNAs associated with phenotypes. a The relationships between lincRNA expression and various demographic and technical features as well as neuropathologic outcomes are displayed in a network diagram. All lincRNA with a significant (p < 0.00016) or suggestive (p < 0.05) association are included in the network. These relationships were extracted via feature selection from linear models which explain an lincRNA’s expression level by using age at death (Age), neuronal composition (NNLS), sex, study of origin (Study, ROS or MAP), post-mortem interval (PMI), RNA integrity number (RIN), neuritic amyloid plaques (NP), neurofibrillary tangles (NFT) and pathological AD diagnosis (AD). Each node in the diagram represents a lincRNA (small circles) or a demographic or a neuropathologic variable (large circles). The significantly associated lincRNA (Table 1) are represented by rectangles. A red edge represents a positive association between the lincRNA and trait; a blue edge denotes an inverse association. As most of the lincRNAs are associated with RIN, the associations with RIN were removed from this figure for clarity. b The association of linc-CTSD-3 expression with each variable is plotted after fitting a model including all variables. The standardized effect size of each of these explanatory variables and 95% confidence intervals are shown. This illustrates the fact that non-coding RNAs significantly associated with AD pathology can be strongly influenced by confounding variables, but that these effects are independent

We have therefore investigated two different types of non-coding RNAs, miRNA and lincRNA, which may play a role in regulating gene expression in AD. Interestingly, while the expression level of different pathologic AD-associated miRNAs are fairly strongly correlated to one another, they are not strongly associated with lincRNAs, and the lincRNAs are only modestly correlated with one another (Fig. 2d). This suggests that, while related, the molecular mechanisms to which these non-coding RNAs contribute in the context of AD may be largely different, with greater coordination among miRNA effects.

Discussion

The number of individuals profiled in our cohorts and the prospective nature of the brain collection in these longitudinal studies of aging make our dataset a valuable resource for exploring, in greater detail, the role of miRNAs that have previously been associated with pathologic AD: our analysis meets a need for studies in larger datasets [17, 38] and offers a potentially less biased perspective of the disease than the comparison of AD cases and controls that are pulled from a brain bank to fit certain diagnostic criteria and were not collected prospectively. We see that the role of miRNAs in AD begins to be resolved in terms of AD’s component pathologies: linking a specific miRNA with either NP or NFT. These associations with a pathologic feature relating to either amyloid or tau pathology are important to generate hypotheses that can be tested in future studies, particularly in mouse models that may recapitulate only one of these pathologies. Aside from the two validated miRNAs, we found suggestive evidence of association (p-values <0.05) for 8 of the 48 previously proposed pathologic AD-associated miRNAs, demonstrating that additional miRNAs may have smaller effects on AD. While all of the other studies age-matched their subjects, none explicitly modeled sex or other technical covariates in their analysis. In particular, we observe strong associations with RIN score, a technical measure capturing the quality of the RNA sample that led to spurious associations when not accounted for in the analysis. One study [12] reports the overall RIN scores for their hippocampal and prefrontal cortex samples; however, they do not provide a description of how this measure correlates with AD, and the incomplete RIN score data reported by another study [13] shows some potential differences between AD and controls. By accounting for the covariates that we have measured, we not only enhance confidence in the miRNA reported as being associated with AD pathology but also provide an opportunity to speculate on the reason why some miRNA may not have been validated due to the technical, clinical and demographic differences among subjects selected for different studies.

Our study has certain limitations, including the advanced average age at death (88 years) in these cohorts and the fact that they are representative of the older population but are not truly population-based, which limit the generalizability of our results. However, the high rate of autopsy (>90%) among study subjects ensures that our results are representative for the entire study population, which consists of subjects who are non-demented at the beginning of the study. Further, the data that we analyzed was generated from the cortex (gray matter). This is a practical compromise to attain large sample sizes, but it presents a challenge for future work as it is not clear which cell type may be driving a particular association. While accounting for the proportion of neurons in the tissue addresses some of the concerns that relate to the role of changes in the relative frequency of cell populations, future work in purified cell populations will be needed to resolve these questions more fully. The Nanostring technology used to measure miRNA expression could also introduce a technical bias in the replication of previously reported miRNA, and the use of probes to measure miRNA with Nanostring and of a predefined reference for aligning lincRNA data will limit the discovery of unannotated non-coding RNAs. Finally, we cannot comment on causality since we have a performed a cross-sectional analysis of brain tissue.

The availability of transcriptome-wide mRNA data from the same RNA samples in a large (n = 540) subset of the ROSMAP subjects profiled for miRNA also provides a rare opportunity to directly explore the relation of miRNA with their putative target mRNAs and of miRNAs with a different class of non-coding RNAs, lincRNAs, in human tissue. Some of the lincRNAs are associated with AD pathology but their expression appears to be largely independent of the miRNAs. On the other hand, our large autopsy-derived mRNA sequence data has identified several different molecular pathways whose component genes have mRNA levels that are associated with AD and are targeted by AD-associated miRNA. The robustness of these results is nicely demonstrated by our lead miRNA, miR-132, which has been validated to be associated with AD in prior targeted studies [12, 13, 39] and for which selected putative target genes have been evaluated in brain samples, including EP300 and SIRT1 [12, 14]. Here, we not only refine earlier observations by showing that the effect of miR-132 is mediated by the accumulation of amyloid pathology but also expand prior targeted studies of downstream effects by organizing the target genes in pathways to highlight cellular functions, such as protein acetylation that appear to be targeted by alterations in miR-132.

The functional consequences of lincRNAs remain poorly understood, and we therefore could not repeat our pathway analysis with this subset of non-coding RNAs. Since lincRNAs may influence the coding RNAs found in the same locus, we did evaluate the association of the neighboring protein-coding genes with AD and the neuropathologic outcomes, but none of the coding transcripts were significantly associated. We also note that the previously reported BACE1-AS lincRNA [20] was not in the reference used in this study, and therefore could not be evaluated in our study. The positive association of BACE1 mRNA expression with AD was not replicated in our cohort (P = 0.63).

With our analyses, we have therefore begun to use autopsy data from a large set of well-characterized human subjects that capture the heterogeneity of older human brains to resolve which aspect of AD-related pathology is influenced by each miRNA of interest. The two well-validated miRNAs in AD illustrate this well: while miR-132 and mir-129-5p are strongly associated with the correlated amyloid and tau pathology measures, both miRNAs are more strongly associated with amyloid pathology than with the accumulation of Tau pathology when the pathologies are included in the same model. This leads to very different experimental paths to further dissect the mechanism of these miRNAs in AD. In addition, using these quantitative pathologic traits that are more precise than a categorical diagnosis of AD, we find that some new miRNAs, such as miR-99b, that may have a stronger effect on a specific pathology, such as Tau/NFT. lincRNAs also appear to be involved, but the downstream consequences of these non-coding RNAs remain unclear. Nonetheless, the lincRNA associations bring another dimension to the broad narrative that emerges from our report: that molecular changes associated with AD include an important alteration in the regulation of cortical transcription, which is consistent with prior reports of epigenomic changes in certain model systems such as Drosophila melanogaster DNA methylation profiles and the involvement of REST in AD [40, 41]. This narrative is also illustrated by miR-200 where the simple analysis of the miRNA alone is suggestive but not convincing of association with AD; however, an integrated analysis that also considers alterations of miR-200 target genes prioritizes this miRNA and downstream transcriptional changes in anion transporters for further evaluation. Such integrated analyses of complementary data may be helpful to resolve the broader perspective on alterations in cellular function in AD. With this manuscript, we therefore provide a robust foundation of detailed neuropathologic associations that set the stage for a new generation of integrative analyses that consider different molecular measures generated from the same subjects and allows for the direct modeling of the complex phenotypic and molecular heterogeneity of the aging population at risk of AD and other dementia.

Conclusions

By studying cortical levels of non-coding RNA in a large prospectively recruited cohort of well-characterized human subjects we had a unique opportunity to explore their associations with measures of both neuritic plaques and neurofibrillary tangles as well as the target mRNAs that they may be regulating. This work provides a robust foundation for future hypothesis-driven work to further dissect the mechanism of the reported associations to specific neuropathologic features of AD.

Demographic information. Table S3. A table of the correlations between each of the pertinent covariates, Study, Age, Sex, NNLS, PMI, RIN, NP, NFT and AD. Figure S1. Behavior of miRNAs implicated in other experiments. Figure S2. Differentially expressed miRNAs and lincRNAs. Figure S3. Validation of RIN association. (DOCX 483 kb)

A table to assess whether the various demographic and neuro-pathological variables may be confounding the relationship between miRNA expression and pathological AD. A linear model is used to model miRNA expression by Study, Age at death, Sex, neuronal proportion (NNLS), post-mortem interval (PMI), RNA Integrity number (RIN) and pathological AD. The effect of pathological AD on the miRNA is reported (coef.AD) as well as the probability of observing this effect or something more extreme assuming that it should be zero (pval.AD). Then for each of the covariates a model is fit without including this covariate in the model and the effect of AD on miRNA expression is calculated again; e.g. pval.AD (−Age) is from the model not including Age at death. (CSV 28 kb)

Reported are the p-value and beta coefficient for the association of each lincRNA with AD, NP and NFT from a series of linear models which include Study, Age at death, Sex, neuronal proportion (NNLS), post-mortem interval (PMI) and RNA Integrity number (RIN) as covariates. (CSV 69 kb)