Whole transcriptome profiling of Late-Onset Alzheimer’s Disease patients provides insights into the molecular changes involved in the disease

Alzheimer’s Disease (AD) is the most common cause of dementia affecting the elderly population worldwide. We have performed a comprehensive transcriptome profiling of Late-Onset AD (LOAD) patients using second generation sequencing technologies, identifying 2,064 genes, 47 lncRNAs and 4 miRNAs whose expression is specifically deregulated in the hippocampal region of LOAD patients. Moreover, analyzing the hippocampal, temporal and frontal regions from the same LOAD patients, we identify specific sets of deregulated miRNAs for each region, and we confirm that the miR-132/212 cluster is deregulated in each of these regions in LOAD patients, consistent with these miRNAs playing a role in AD pathogenesis. Notably, a luciferase assay indicates that miR-184 is able to target the 3’UTR NR4A2 - which is known to be involved in cognitive functions and long-term memory and whose expression levels are inversely correlated with those of miR-184 in the hippocampus. Finally, RNA editing analysis reveals a general RNA editing decrease in LOAD hippocampus, with 14 recoding sites significantly and differentially edited in 11 genes. Our data underline specific transcriptional changes in LOAD brain and provide an important source of information for understanding the molecular changes characterizing LOAD progression.


Results
Transcriptomic Profiling of Hippocampus of LOAD patients. We used RNA-seq to examine the transcriptomic profile of the hippocampal CA1 region of a cohort of six patients affected by LOAD (samples AD1-6), six patients affected by PD (samples PD1-6 samples) and six cognitively normal controls (samples Ctrl1-6) ( Table 1). All sampled individuals were Caucasian males with a comparable age of death, while LOAD patients had comparable dementia status (Braak V or VI). All frozen tissues preserved total RNA integrity as indicated by RIN values (5 to 8). PD patients were included as "disease-control patients", to highlight transcriptomic changes that are LOAD specific, rather than common to general neurodegenerative processes. We chose the hippocampal region CA1 as the area of investigation for its relevance to memory processes and as this area is, together with the entorhinal cortex, the first affected by the pathogenic mechanisms associated with LOAD.
Directional RNA sequencing of rRNA-depleted total RNA generated an average of ∼140 million reads per sample, of which 88 to 92% could be aligned to the reference genome, suggesting excellent coverage and sequencing depth with low ribosomal RNA contamination (Supplementary Table S1). RNA-seq data analysis identified 25,272 annotated genes as expressed in, at least, one of the sequenced samples. Principal Component Analysis (PCA) indicated anomalous behavior for Ctrl4 and AD2 samples, which were thus discarded from all further analyses (Supplementary Figure S1), while no outlier was observed among PD samples.   CuffDiff2 16 identified a total of 2,122 genes as differentially expressed (p-value ≤ 0.05 after false discovery rate correction) between LOAD samples and controls, of which 2,075 were protein coding genes and 47 were lncR-NAs (of which 23 lincRNAs). Of the 2,075 differentially expressed protein coding genes, 789 were upregulated and 1,286 were downregulated, and of the 47 differentially expressed lncRNAs, 19 were upregulated (of which 11 lincRNAs) and 28 were downregulated (of which 12 lincRNAs) in LOAD samples. The complete list of deregulated genes is reported in Supplementary Table S2. Nineteen protein-coding genes were found to be differentially expressed between hippocampal PD samples and controls (Supplementary Table S3). Eleven of these genes (ANKRD22, HAMP, HLA-DRA, HSPA6, CD14, FCGBP, HSPB1, HSPA7, BAG3, SERPINH1 and TNFRSF1B) were also present in the list of the 2,075 protein coding genes deregulated in LOAD patients. We analyzed these shared deregulated genes using the STRING tool 17,18 , and noted that they were involved in neurodegenerative and inflammatory pathways such as defense response, response to unfolded proteins, positive regulation of immune response and toll-like receptor (TLR) signaling. Accordingly, the common LOAD and PD differentially expressed protein coding genes were removed from the initial list of 2,075 LOAD deregulated genes and, in order to determine which pathways are affected in LOAD patients, the resulting 2,064 genes that were differentially expressed in a LOAD-specific manner, were analyzed using the IPA tool 19 . As shown in Table 2, we found that activities related to the regulation of important neurological functions, including synaptic long-term potentiation, neuronal signaling, axonal guidance signaling and mitochondrial dysfunction, were significantly enriched among this differentially expressed gene set. Interestingly, we noted that the products of several genes among this set are involved in pathways regulated by the miR-132/212 cluster, known to be down-regulated in AD [20][21][22][23][24][25][26] . In particular, we found the up-regulation of some direct targets of miR-132, such as ITPKB (log2FC: 1.33, padj: 0.0024), involved in Tau phosphorylation 27 , TLR6 (log2FC: 1.36, padj: 0.0024), IL6R (log2FC: 1.37, padj: 0.0014) and IRAK3 (log2FC: 1.33, padj: 0.0048), involved in inflammatory signaling 28 . Moreover, we found BDNF, a positive upstream regulator of miR-132/212 transcription 29 , to be down-regulated (log2FC: −2.6, padj: 0.0041) as was RASAL1, which regulates angiogenesis through mir-132 30   We validated the hippocampal RNA-seq differential expression data by RT-qPCR using an enlarged cohort of nine controls (five original controls plus four new controls) and nine LOAD patients (five original LOAD patients plus four new LOAD patients), chosen to fit criteria of selection of the original cohort used in the RNA-seq analysis. Given that the exclusion of Ctrl4 from RNA-seq analyses caused a reduction of the average age of control samples with respect to LOAD samples, we selected control samples with a higher average age than the original control groups (see Supplementary Table 4). From the list of deregulated genes, we selected 21 protein  Table S5). RT-qPCR analysis confirmed the differential expression levels for all 8 up-regulated genes (CPLX3, NR4A2, GRIK3, TESPA, SLCO4A1, SERPINA5, ADAM33, SERPINA1), both in the original and extended cohorts (p-value ≤ 0.05) ( Fig. 1a and Supplementary Figure S2). Regarding the down-regulated genes, RT-qPCR analysis confirmed the expression trend resulting from RNA-seq data for all 14 genes (BHLHE22, PRSS12, NEUROD6, PCDH8, NRN1, DUSP4, CAMK1D, NEUROD1, GRIA1, SYTL5, PRKCG, ARC, SCN11A and LOC400891), in the original and/or in the extended cohort (p-value ≤ 0.05), although two genes (NEUROD1 and SCN11A) did not show a statistically significant reduction of expression in the RT-qPCR analyses ( Fig. 1b and Supplementary Figure S3). In addition, for all 22 deregulated genes, a positive and highly significant correlation was also found between the estimates of fold change in expression level from RNA-seq data and RT-qPCR results (R 2 : 0.9187) (Fig. 1c).

MicroRNA Profiling in LOAD Brain regions.
To gain insights into the complexity and specificity of transcriptomic changes in LOAD, we next extended our analysis to miRNA profiling. In order to identify miRNAs whose deregulation is related to the spatio-temporal progression of the disease, we analyzed three different vulnerable brain areas, the hippocampal region CA1, the middle temporal gyrus Brodmann's Area 21 and the middle frontal gyrus Brodmann's Area 46, in the same LOAD patients and controls used in the RNA-seq experiment. This choice was determined by existing evidence that these brain regions are involved in distinct stages of the neurodegenerative process in terms of development of the pathology and functional outcome, with the hippocampus affected first, followed by the temporal lobe and, finally, the frontal lobe 4,5 . In the hippocampal analysis, we also included PD patients, in order to identify shared deregulated miRNAs between LOAD and PD hippocampal samples.
Small RNA-seq generated an average of ∼4.5 millions of reads per sample, with 90% of bases showing Q scores above 30 (Supplementary Table S6). Mapping of reads from hippocampal samples identified 1,018 annotated miRNAs expressed in, at least, one of the sequenced samples. Considering a coverage of at least 10 reads/miRNA, this number was reduced to 407. Differential expression analysis identified 4 miRNAs (miR-184, miR-34c-3p, miR-375 and miR-132-5p), that were differentially expressed (all down-regulated) in LOAD patients respect to control samples (padj ≤ 0.05) ( Table 3). In addition, 40 miRNAs were differentially expressed in hippocampal PD samples with respect to controls (padj ≤ 0.05) (11 up-regulated and 29 down-regulated) (Supplementary Table S7), but none of these were present in the list of deregulated LOAD hippocampal miRNAs. PD deregulated miRNA targets are enriched in some neuronal and inflammatory pathways also identified for LOAD miR-NAs, consistent with the involvement of such pathways in general responses to brain injury (p-value ≤ 0.05) (Supplementary Table S8).  Mapping of reads from middle temporal gyrus samples identified 1,205 annotated miRNAs as expressed in, at least, one of the sequenced samples. Considering a coverage of at least 10 reads/miRNA, this number was reduced to 442. Differential expression analysis identified 10 miRNAs deregulated (padj ≤ 0.005) in LOAD patients. Five of these (miR-501-3p, miR-10a-5p, miR-320a, miR-28-3p and miR-30a-3p) were upregulated and 5 (miR-539-5p, miR-132-5p, miR-212-3p, miR-132-3p and miR-212-5p) were downregulated (Table 3).
KEGG pathway annotations using DIANA miRPath tool 31 allowed us to identify 6 pathways that are enriched in genes with predicted targets for the miRNAs deregulated in every brain region studied (p-value ≤ 0.05). These pathways include MAPK and neurotrophin signaling, axon guidance, long-term potentiation, glutamatergic and cholinergic synapses (Table 4).
miRNA profiling confirmed the down-regulation of the miR-132/212 family -previously identified as down-regulated in LOAD 20-26 -in each of the three different regions of studied LOAD brains, consistent with a role in LOAD pathogenesis for these miRNAs and supporting the consistency of our analyses. No other miRNA showed differential expression in all three regions of LOAD brains, although the identified deregulated miRNAs of each region affected common pathways, involved in the regulation of different aspects of neuronal cellular function.
RT-qPCR assays were used to validate the deregulated miRNAs of each brain region. For the hippocampal region, we used both the original and the enlarged cohorts, as employed for the RNA-seq data validation (Table 1 and Supplementary Table S4). We confirmed the downregulation of all miRNAs identified in miRNA-seq  Table 4. Pathways identified by DIANA miRPath analysis, affected by the deregulated miRNAs from hippocampus, temporal and frontal gyrus of LOAD patients. The Fisher-exact test P-value is reported in each column.  Figure S4). For the middle temporal gyrus and middle frontal gyrus, we used the original cohort of 5 controls and 5 LOAD samples. In both brain regions, for all deregulated miR-NAs, RT-qPCR assay confirmed the expression trend identified in miRNA-seq (in the middle temporal gyrus: p-value ≤ 0.05 for miR-10a-5p, miR-28-3p miR-539-5p and miR-132/212 cluster members; in the middle frontal gyrus: p-value ≤ 0.001 for miR-132/212 cluster members) ( Fig. 2b and Supplementary Figure S4).

miR-184 targets the 3'UTR of NR4A2 transcript.
We used the microRNA Target Filter, a microRNA target prioritization tool available within IPA 19 , the DIANA miRPath 31 algorithm and Miranda 32 to compare the deregulated miRNAs and transcripts identified in LOAD hippocampus. All tools predicted with high confidence, using TargetScan Human 33 as source, the interaction between miR-184 and the NR4A2 and NRN1 transcripts. Interestingly, we found an inverse expression correlation for NR4A2 and miR-184, as NR4A2 was over-expressed ( Fig. 1a and Supplementary Figure S2) and miR-184 was down-regulated ( Fig. 2a and Supplementary Figure S4) in our RNA-seq data. We verified the predicted interaction miR-184/NR4A2 using a luciferase reporter assay. We cloned the 3' untranslated region (UTR) of the NR4A2 transcript downstream of the luc2 firefly luciferase open reading frame in the pMIR-Reporter Luciferase miRNA Expression Vector (Fig. 3a). The recombinant vector was transfected into H1299 cells either with a negative control miRNA, with the miR-184 mimic alone or together with a miR-184 inhibitor (anti-miR-184). We found that miR-184 caused a significant reduction, of about 25%, of the normalized luciferase expression with respect to the control miRNA, while the co-transfection of equimolar quantity of the miR-184 mimic and the antimiR-184 didn't produce any significant difference with respect to the control miRNA (Fig. 3b). These results indicate that NR4A2 transcript can be directly targeted by miR-184, consistent with the expression levels of NR4A2 that we found in LOAD hippocampus. On the contrary, NRN1 was down-regulated ( Fig. 1a and Supplementary Figure S2) in our RNA-seq data, as miR-184. We also tested the predicted interaction miR-184/NRN1 by luciferase reporter assay. We found that miR-184 did not produce any significant reduction in luciferase expression with respect to the control miRNA, thus excluding a NRN1 regulation from miR-184 (Supplementary Figure S5b).

RNA editing analysis.
RNA editing is an important post-transcriptional process that alters the genetic blueprint of an organism by specific modifications in primary RNAs. In human, it mainly involves the deamination of adenosines to inosines by the family of adenosine deaminase acting on RNA (ADAR) enzymes acting on double RNA strands and its deregulation has been linked to a variety of neurological and neurodegenerative disorders 34,35 . We investigated A-to-I RNA editing alterations in LOAD hippocampal tissues using RNA-seq data.
First, for each sample, we calculated the global editing activity through the Alu editing index (AEI) as it represents the weighted average editing level across all expressed Alu sequences 36 . We found a slight, but non-significant, increase of editing activity in LOAD samples compared to controls (Fig. 4). Interestingly, the RNA editing activity at recoding sites, that have a functional role in brain, appeared decreased in LOAD hippocampal samples (p-value < 0.05), as attested by the Recoding Editing Index (REI) calculated as the weighted average of editing levels over all known recoding sites from the REDIportal database 37 (Fig. 4). Focusing on recoding editing sites, we found that 14 out of 1,585 sites were significantly and differentially edited (t-test followed by 10% FDR   multiple-testing correction) in LOAD samples compared to controls (Table 5). In particular, we confirmed the deficient RNA editing at the glutamate receptor subunit B, GRIA2, Q/R site, previously described in LOAD hippocampus 38 and at other five sites (COPA I/V site, GRIK1 Q/R site, GRIK2 I/V site, GRIA3 R/G site, GRIA4 R/G site), identified in a previous study that analyzed the editing levels at 180 recoding sites in different LOAD brain regions 39 . Finally, we evaluated the expression of ADAR genes in RNA-seq data. We found that ADAR1 and ADAR2 were not significantly differentially expressed in LOAD samples respect to controls although the ADAR2 expression level appeared slightly decreased in LOAD (Fig. 5). Interestingly, we found that the ADAR3 locus, encoding for an inactive adenosine deaminase protein, primarily expressed in the brain, was significantly over-expressed in LOAD samples (Fig. 5).

Discussion
LOAD is an age-related neurodegenerative disorder, for which whole transcriptome profiling represents an informative approach to characterize end-stage alterations resulting from interactions among multiple genetic, epigenetic and environmental factors.  49 ), different disease stages and brain areas, variable post-mortem interval (PMI) affecting RNA quality, and finally different bioinformatics protocols to analyze data. An accurate set of genes and pathways deregulated in LOAD is, thus, far from established.
In an attempt to overcome these issues, we have performed, using NGS technology, a comprehensive transcriptomic profiling, including coding and non-coding RNAs (lncRNAs and miRNAs) of post-mortem hippocampal tissues from three cohorts, one consisting of LOAD patients, one consisting of cognitively normal controls and the last consisting of PD patients. To restrict the variability, samples were selected to minimize gender, ethnicity, age at death and disease stage differences between patients and controls and to have a PMI not longer than 23 hours. Indeed, it has been demonstrated that all categories of biomolecules remain stable up to 48 hours post-mortem 50 , consistent with RIN values of RNA obtained from all samples. We used the 'whole' frozen brain tissue as laser-micro-dissection, targeting exclusively the neuronal soma, misses a variable proportion of transcripts that are transported for pre-or post-synaptic translation 51 . We used both grey and white matter since several studies, which examined grey and white matter separately, observed similar deregulated level of expression of genes involved in LOAD 40 . We included PD patients as "neuro-inflammatory disease" controls in order to identify LOAD specific transcriptomic changes, not shared with general neuro-inflammation processes. In fact, considering that affected regions also contain inflammatory cells, linked to neurodegeneration and not normally associated with healthy tissue (such as macrophages and T lymphocytes), the comparison involving the control group may involve tissues with a quite distinct cellular composition. Considering that our LOAD patients had high Braak score (advanced disease), to exclude that the differential expressed genes might be affected by an imbalance between the original cellular populations between LOAD and control samples, we analyzed also the expression level of some cell-type specific genes in our RNA-seq data. We did not find any statistically significant changes in the expression of neuronal markers: DCX (log2FC: −0.43, padj: 0.50), RBFOX3 (log2FC: −0.37, padj: 0.46) and TUB (log2FC: −0.43, padj: 0.43); astroglial markers: GFAP (log2FC: 0.23, padj: 0.94), AQP4 (log2FC: 0.74, padj: 0.52) and SLC1A3 (log2FC: 0.59, padj: 0.37); microglial markers: AIF1 (log2FC: 1.21, padj: 1) and CX3CR1 (log2FC: −0.07, padj: 0.94). This observation indicates that LOAD tissues analyzed in our study preserved the original cellular population, and that the results of the differential analysis were not affected by differential expression of cell-type specific genes.
RNA-seq identified a set of 2,064 specific differentially expressed transcripts in the hippocampus of LOAD patients, which includes coding and long non-coding RNAs. The deregulated genes are enriched in roles related to the regulation of important neural processes, such as neurogenesis, synaptic vesicle trafficking, long-term potentiation, neurite outgrowth and hyper-phosphorylation of Tau protein, consistent with their involvement in LOAD pathogenesis. For the data validation by RT-qPCR, we enlarged both the control and the LOAD sample groups, including control samples with a higher average age with respect to the original control groups (Supplementary  Table S4), as we were aware that, with the exclusion of Ctrl4 in the RNA-seq analysis, the ages of control samples became slightly smaller than those of LOAD samples. The deregulated expression of the selected group of transcripts was confirmed both in the original and in the enlarged cohorts, thus supporting the consistence of our bioinformatics analysis and suggesting that the difference of age between the two groups did not significantly affect our results.
In order to identify miRNAs differently expressed during LOAD neurodegenerative progression, with the aim of uncovering new diagnostic and prognostic markers for the disease as well as new potential drug targets, we extended the miRNA profiling from the hippocampus to the middle temporal gyrus (Brodmann's Area 21) and middle frontal gyrus (Brodmann's Area 46) of the same LOAD patients. These specific brain regions were chosen according to evidence that they are involved in distinct stages of the neurodegenerative process, with the hippocampus affected first, followed by the temporal lobe and at last the frontal lobe 4,5 . Our data further sustain the role of the miR-132/212 cluster in AD pathogenesis 52 , as its expression was deregulated in all three brain regions of the LOAD patients. Additionally, sequencing and RT-qPCR analyses demonstrate the differential expression of other specific miRNAs in each of the three brain regions considered, three miRNAs in the hippocampus, six miRNAs in middle temporal gyrus and three miRNAs in middle frontal gyrus of the same LOAD patients. These results suggest that miRNA expression changes in LOAD might be brain area specific and could differently contribute to disease progression.
Interplay between LOAD hippocampal deregulated genes and miRNAs was supported by in silico predictive tools that revealed a possible interaction between the NR4A2 and NRN1 transcripts and miR-184. Both NRN1 and NR4A2 predicted miRNA target sites can be considered canonical (7-8 nucleotides that match perfectly with miR-184 seed region) ( Fig. 3a and Supplementary Figure S5a). Nevertheless, it is well known that many factors beyond pairing capacity determine functional interactions in vivo and, consequently, the presence of a putative miRNA binding site in a transcript is not sufficient to assert its functional relevance in the absence of experimental validation 53,54 . Accordingly, we performed luciferase assays that showed that miR-184 is able to target the 3'UTR of the NR4A2 transcript but not of the NRN1 transcript indicating that only NR4A2 expression is regulated by miR-184 in vivo (Fig. 3b and Supplementary Figure S5b). In order to understand the different target specificity in the two transcripts, we considered the mammalian orthologs of NR4A2 and NRN1 genes. We observed higher conservation of the seed-matching region in NR4A2 than in NRN1 binding site, consistent with the observation that the conservation level of the seed matched site in 3′UTRs of orthologous genes is a relevant feature in predicting miRNA interaction 55 . Notably, our sequencing data recover an inverse correlation of expression between NR4A2, that was over-expressed, and miR-184, that was down-regulated in LOAD patients, while NRN1, like miR-184, was strongly down-regulated in LOAD patients.
MiR-184 was previously shown to trigger neural stem cell (NSC) proliferation and inhibit differentiation by repressing the NSC fate-regulator Numblike 56 . Its deregulation could negatively affect ongoing neurogenesis in the hippocampal neurogenic niches of adult brains, and which are known to be compromised in AD 57 . Several gene targets for miR-184 have been described, including mediators of neurological development and apoptosis, such as E2F1 and DP 58 .
NR4A2, also known as Nurr1, encodes an orphan nuclear receptor, belonging to the nuclear receptor subfamily 4A (NR4A), which also includes NR4A1 and NR4A3 59 . It is expressed in several brain areas and plays important roles in various neural processes including cognitive functions and/or long-term memory in forebrain areas and hippocampal formation. Immune-fluorescence staining with a specific NR4A2 antibody in 5XFAD mice showed that the NR4A2 protein was prominently localized in brain areas with Aβ plaque accumulation, highlighting a possible involvement of this protein in AD pathology 60 . The potential regulation of NR4A2 by miR-184 requires further functional testing but leads us to suggest that the loss of expression of miR-184 in LOAD hippocampus may contribute to NR4A2 accumulation in Aβ plaques.
While the down-regulation of miR-132/212 cluster in AD is well established by the current study (Table 3 and Fig. 2) and by other published studies 21,22 , we also identified deregulation of several genes whose products are involved in miR-132/212 related pathways. BDNF (brain derived neurotrophic factor), which positively regulates miRNA-132/212 cluster transcription via the ERK1/2-dependent (but CREB-independent) pathway 29 was down-regulated in LOAD brains. MiR-132 regulates Tau phosphorylation by targeting ITPKB (inositol-trisphosphate3-kinase B), which is recovered as up-regulated in LOAD brains in the current study. ITPKB is known to enhance both ERK1/2 kinases and BACE1 activity, and to be elevated in LOAD brains 27 . Finally, miR-132 is a negative regulator of inflammation 28 , and a positive regulator of vasculogenesis 30 . In our SCIeNtIfIC REPORTS | (2018) 8:4282 | DOI:10.1038/s41598-018-22701-2 study, miR-132 down-regulation is correlated with the up-regulation of inflammatory signals (such as TLR6, IL6R and IRAK3) and the down-regulation of RASAL1, a positive effector of angiogenesis.
Finally, we investigated RNA editing levels in the LOAD hippocampus since alterations in RNA editing have been implicated in AD 39 and in other neurodegenerative diseases 34,35 . We found a general decrease of RNA editing activity at recoding sites (Fig. 4), in agreement with a previous study that reported decreased editing levels in different regions of AD patients' brains, analyzing 180 recoding sites by a targeted resequencing approach supplemented by a microfluidic-based high-throughput PCR coupled with next-generation sequencing 31 . We identified 14 recoding sites with significant and differential editing levels in LOAD hippocampal tissues compared to controls (Table 5). Among these sites, we found down-regulation of RNA editing levels in several glutamate receptors: GRIA2, GRIA3, GRIA4, GRIK1, GRIK2. Regarding the Q/R site of GRIA2, it has been demonstrated that editing of this site is essential in neurons, as incorporation of unedited GRIA2 subunits results in an increase in glutamate receptors that are permeable to calcium, leading to excitotoxicity and cell death of neurons 61 . Our finding of low RNA editing levels at the Q/R site of GRIA2 is in line with previous findings that reported this deficiency in the hippocampus of AD patients, suggesting that it may be an early event that precedes neuronal demise 38 . Three ADAR enzymes are present in mammals: ADAR (ADAR1), ADARB1 (ADAR2) and ADARB2 (ADAR3). ADAR1 and ADAR2 are the main catalytic enzymes, accountable for all A-to-I editing events, while ADAR3 has not been shown to have deaminase activity in vitro and has no known in vivo target 62 . Interestingly, while changes in ADAR1 and ADAR2 expression levels were not observed in our RNA-seq data, we recovered a significant increase of ADAR3 transcript levels in LOAD hippocampal tissues (Fig. 5). Recently, it has been reported that ADAR3 is able to bind directly to GRIA2 pre-mRNA transcripts, thus competing with ADAR2 activity and inhibiting the RNA editing at the Q/R site of GRIA2. This mechanism could explain the reduced editing levels at the Q/R site of GRIA2 observed in human glioblastoma, where ADAR2 is normally expressed while ADAR3 expression is increased 63 . Analogously, we can suggest that the inhibition of ADAR2 activity by ADAR3 may be responsible for decreased RNA editing levels at the Q/R site of GRIA2 as well at other recoding sites in the LOAD brain.
In conclusion, post-mortem samples do not permit the study of transcriptomic changes during early stages of LOAD and can only evaluate the transcriptome alterations that characterize the final stages of neurodegeneration. Furthermore, the use of bulk RNA does not allow the identification of cell specific transcriptional deregulation, but currently, single-cell transcriptome analysis is not feasible on frozen post-mortem tissues. Despite these limitations, we believe that our study, providing a list of transcriptomics alterations in LOAD hippocampus, may represent an important resource for further molecular investigations aimed at better understanding the molecular mechanisms involved in LOAD and for the potential identification of novel therapeutic targets in LOAD.

Methods
Brain Samples. Frozen post-mortem brain tissue samples from cognitively normal controls, LOAD patients and PD patients were obtained from several non-profit brain banks: the NICHD Brain and Tissue Bank for Developmental Disorders (University of Maryland School of Medicine, Baltimore, USA), the London Neurodegenerative Diseases Brain Bank (King's College Hospital, London, UK) and the Parkinson's UK Brain Bank (Imperial College London, UK) ( Table 1 and Supplementary Table S4). All patients provided their informed consent to the brain bank. All methods were carried out in accordance with relevant guidelines and regulations and the study was approved by the Institutional Review Board of the Institute of Biomembranes, Bioenergetics and Molecular Biotechnologies, National Research Council. Three brain regions, hippocampus (CA1 region), middle temporal gyrus (Brodmann's Area 21) and middle frontal gyrus (Brodmann's Area 46) were used in this study. All tissues contained both grey and white matter. Control subjects did not have a history of neurological disease or indications of brain abnormalities at tissue level, as determined at autopsy. LOAD patients were positive for Aβ plaques and NFTs and were selected according to dementia status (Braak V or VI). PD patients were positive for Lewy bodies in the substantia nigra.

RNA Extraction.
Total RNA and small RNA fractions were selectively extracted using the mirVana ™ miRNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions. RNA was quantitatively and qualitatively evaluated using NanoDrop 2000c (Thermo Fisher Scientific) and Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA), respectively. For total RNA, a RNA Integrity Number (RIN) ranging from 5.2 to 7.9 was obtained, which was considered acceptable as deriving from post-mortem tissues. Preprocessing and analysis of RNA-seq data. Reads with a Phred quality score (Q) > 20, a length higher of 50 nt and homopolymeric tract lower of the 50% of the total read length, were selected. In detail, reads in FASTQ format were inspected using FASTQC program (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/) whereas adaptors and low quality regions were trimmed using Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Cleaned reads were aligned onto the complete latest human genome (assembly hg19) by means of GSNAP 64  in the binary BAM format by SAMtools 65 and basic statistics were calculated using Picard tools (CollectRnaSeq Metrics.jar) (http://picard.sourceforge.net/). Transcriptome quantification and differential expression of coding and non-coding RNA was performed using CuffDiff2 (http://cufflinks.cbcb.umd.edu/) 16 software version 2.1.1 (using as main parameters:-library-type fr-secondstrand-labels Ctl,Als -u -b hg19.fa -M rRNA.gtf refgenes.gtf). The Principal Component Analysis (PCA) was conducted by ClustVis (http://biit.cs.ut.ee/clustvis/) 66 on gene expression values obtained by CuffDiff2 on RNA-seq aligned reads. After the exclusion of outlier samples (Ctrl4 and AD2), CuffDiff2 analysis was repeated on the remaining samples. Reference human transcriptome was obtained from iGenomes repository (http://cole-trapnell-lab.github.io/cufflinks/igenometable/) and annotations for rRNA genes were downloaded from UCSC genome browser selecting the RepeatMask table. Differentially expressed genes were selected if |log2(FC)| > 1 and corrected P value < 0.05.

RNA
Differentially expressed coding RNAs were used as input to perform pathway enrichment analysis by IPA system (Ingenuity ® Systems, www.ingenuity.com), STRING tool (http://string-db.org) 17  Preprocessing and analysis of miRNA-seq data. Adapter sequences were trimmed from raw reads using a custom Python script. Trimmed reads were mapped to human mature miRNA sequences, downloaded from miRBase (v21), using another ad hoc Python script, which required a best alignment with overlap of at least 17 bp, without mismatches, with an annotated mature miRNA sequence to assign miRNA identity. Only miRNAs with a read count higher than 10 were considered in the differential expression analysis. To assess the coherence of miRNA expression profiles, miRNA hierarchical cluster analysis was carried out using the R package pvclust (https://CRAN.R-project.org/package=pvclust) 68 and computing pairwise-distances among samples. Differentially expressed miRNAs were identified using both DESeq2 69 and in-house script methodologies, which generated substantially congruent conclusions. The in-house script method, although based on the log 2 of fold change of expression levels, differs from DESeq2 method since: it does not use mathematical models to estimate the variance associated with the level of expression, it uses the upper quartile normalization system, it compares the log 2 of fold change of expression levels intra-and inter-conditions, and it applies the T test (instead of the Wald test). Potential miRNA-targeted genes and their impact on the biological pathways were assessed using DIANA-mirPath v2.0 (http://diana.imis.athena-innovation.gr) 70 .

Reverse Transcription Quantitative PCR Analysis. Coding and non-coding RNA Expression
Analysis. 1 µg and 4 µg of total hippocampal RNA were used in the reverse transcription reaction (RT) for coding and non-coding RNA expression levels analysis, respectively, using the iScript ™ Advanced cDNA Synthesis Kit (Bio-Rad Laboratories Ltd, Berkeley, California, USA), according to the manufacturer's instructions. Quantitative PCR (qPCR) reactions were performed using 1 µl of the diluted cDNA (1:4) as input, specific TaqMan Gene Expression Assays and following the TaqMan Fast Advanced Master Mix Protocol (Thermo Fisher Scientific). IDs of the TaqMan Gene Expression Assays used are listed in Supplementary Table S5. RT-qPCR experiments were performed in triplicate for each sample on StepOnePlus Real-Time PCR System (Thermo Fisher Scientific). The expression of five housekeeping genes, CYC1, EIF4A2, GAPDH, HPRT1 and β-actin was also analyzed using two software applet for Microsoft Excel, named GeNorm 71 and NormFinder 72 in order to define the most stable housekeeping genes among samples, that resulted in selecting CYC1, EIF4A2 GAPDH genes. To normalize the data, the arithmetic Ct mean of housekeeping genes was evaluated for each sample; then, the transcript expression levels were calculated according to the ΔΔCt method, setting the arithmetic ΔCt mean of control group, as calibrator. The resulting data were represented as log2(ΔΔCt) and were expressed as the mean ± SD. Two-tailed Student's T tests were performed to assess the statistical significance of gene expression levels differences observed between normal and LOAD samples. A p-value < 0.05 was considered statistically significant.
miRNA Expression Analysis. 10 ng of small RNA fraction isolated from middle temporal gyrus, middle frontal gyrus and hippocampal samples were used to perform RT reactions using the TaqMan ™ Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific), according to the manufacturer's instructions. The resulting cDNA, diluted 1:10, was used as input for qPCR using specific TaqMan Advanced miRNA Assays and the TaqMan Fast Advanced Master Mix Protocol (Thermo Fisher Scientific). IDs of the TaqMan Advanced miRNA Assays used are listed in Supplementary Table S9. GeNorm and NormFinder tools were used to define the most stable miRNAs to be used as endogenous controls, that resulted to be: miR-423-3p (478327_mir), miR-181b-5p (478583_mir) and miR-191-5p (477952_mir). To normalize the data, the arithmetic Ct mean of the three stable miRNAs was evaluated for each sample; then, the miRNA expression levels were calculated according to the ΔΔCt method, setting the arithmetic ΔCt mean of control group as calibrator. The resulting data were represented as log2(ΔΔCt) and were expressed as the mean ± SD. Two-tailed Student's T tests were performed to assess the statistical significance of miRNA expression levels differences observed between normal and LOAD samples. A p-value < 0.05 was considered statistically significant. RNA editing analysis. RNA editing candidates were detected using REDItools 74 . The Alu editing index, the weighted average editing level across all expressed Alu sequences, was calculated using custom scripts and according to the methodology described in Bazak et al. 75 . RNA editing levels in recoding sites were assessed using REDItools and providing a list of 1585 known positions from REDIportal database 37 in which RNA editing causes amino acid change. Only positions supported by RNAseq reads in at least three samples per group and showing a median editing level higher than 0.1 were used to calculate dysregulated RNA editing. A t-test was used to assess editing dysregulation for each RNA site. P-values were corrected using Benjamini-Hochberg with FDR = 0.1. ADAR expressions in LOAD and control samples were calculated using CuffDiff2 16 . Statistics. Differentially expressed genes in RNA-seq analysis were defined by using CuffDiff2. miRNA-seq data were statistically analyzed by applying DESeq2 both for data normalization and comparison (Wald test) and in-house procedure relying in the upper quartile normalization and T test for comparisons. For RT-qPCR and luciferase analyses, the statistical significance was assessed by Two-tailed Student's T tests and results were expressed as the means ± SD. For RNA editing analyses, the statistical significance was assessed by using T test followed by Benjamini-Hochberg procedure for multiple test p-values correction. P-values less than 0.05 were considered to indicate statistical significance.