Identifying lncRNA–miRNA–mRNA networks to investigate Alzheimer’s disease pathogenesis and therapy strategy

Alzheimer’s disease (AD), the most common cause of dementia, leads to neuronal damage and deterioration of cognitive functions in aging brains. There is evidence suggesting the participation of noncoding RNAs in AD-associated pathophysiology. A potential linkage between AD and lncRNA-associated competing endogenous RNA (ceRNA) networks has been revealed. Nevertheless, there are still no genome-wide studies which have identified the lncRNA-associated ceRNA pairs involved in AD. For this reason, deep RNA-sequencing was performed to systematically investigate lncRNA-associated ceRNA mechanisms in AD model mice (APP/PS1) brains. Our results identified 487, 89, and 3,025 significantly dysregulated lncRNAs, miRNAs, and mRNAs, respectively, and the most comprehensive lncRNA-associated ceRNA networks to date are constructed in the APP/PS1 brain. GO analysis revealed the involvement of the identified networks in regulating AD development from distinct origins, such as synapses and dendrites. Following rigorous selection, the lncRNA-associated ceRNA networks in this AD mouse model were found to be mainly involved in synaptic plasticity as well as memory (Akap5) and regulation of amyloid-β (Aβ)-induced neuroinflammation (Klf4). This study presents the first systematic dissection of lncRNA-associated ceRNA profiles in the APP/PS1 mouse brain. The identified lncRNA-associated ceRNA networks could provide insights that facilitate AD diagnosis and future treatment strategies.

It was reported that 80% of the human genome is transcribed as noncoding RNAs (ncRNAs) [3], which can manipulate most of the potential biological functions [4]. Over the past two decades, it was demonstrated that ncRNAs, with their specific spatiotemporal expression patterns across various species, are extensively involved in numerous biological processes, including epigenetic regulation, chromatin remodeling, transcription control, and posttranscriptional processing [5]. More and more ncRNAs have been identified for their important roles in the pathogenesis of neurodegenerative disorders.
NcRNAs include microRNAs (miRNAs, ~20 nucleotides in length) and long noncoding RNAs (lncRNAs, larger than 200 nucleotides in length), which play significant roles in the regulation of various biological functions [6]. There is evidence that miRNAs and lncRNAs participate in AD pathophysiology, which includes the formation and development of β-amyloid (Aβ) plaques, neurofibrillary tangles, synaptic loss and neuronal death [7][8][9]. The ncRNAs have been detected by RNA-sequencing (RNA-seq) or microarrays in many organisms, which presents exciting implications for understanding the regulation of basic biological systems and pathophysiological conditions, as well as the development of new therapeutic treatment of many diseases [10].
Recent studies suggest that lncRNAs, circRNAs, pseudogenes and mRNAs may function as miRNA sponges [17][18][19]. They compete with each other through miRNA response elements (MREs) and modulate the progress of many diseases [20,21] including AD [9,6,13]. In addition to post-transcriptional regulation, epigenetic modifications may also play a crucial role in AD pathogenesis [7,22]. These refined and complicated regulatory networks may help to explain why the isolation of a single component (e.g., β-secretase and apolipoprotein E4) failed to fully account for the whole pathogenesis process of AD [23]. It suggests that the underlying mechanisms for the involvement of competing endogenous RNA (ceRNA) in AD remain to be determined.
The rational strategy to gain insight into neurodegenerative diseases such as AD would make the study of lncRNA-associated ceRNA networks comprehensively. The elucidation of lncRNAassociated ceRNA networks in AD might help to develop new therapeutic targets for AD. In APP/PS1 mice which express APP695swe and PS1-dE9 mutations, Aβ can be detected in 6-month old mice, and then extracellular Aβ deposits in the cortex are apparent in 9-month old mice. Moreover, synaptic transmission and long-term potentiation are clearly impaired in 9month old mice [24]. Thus, deep RNA-seq was performed in this study to find lncRNA-associated ceRNA networks in the brain of APP/PS1 mice together with the wild-type (WT) control mice at the 6-and 9month-old stages. RNA-seq is widely used to determine the differential gene-expression profiles that underlie phenotypic differences [25,26]. The data by RNA-seq identifies lncRNA-associated ceRNA networks in the APP/PS1 mouse model of AD (Figure 1), which can contribute to the development of new therapeutic targets and novel diagnostic methods for AD.

Overview of lncRNA and mRNA-seq data
A total of 1,145,853,722 raw reads were generated. 317,603,558 raw reads and 276,260,804 raw reads were generated for 6-month-old WT and APP/PS1 mice respectively, while 286,203,120 raw reads and 265,786,240 raw reads were generated for 9-month-old WT and APP/PS1 mice respectively. After discarding the reads with adapters, poly-N > 10%, or any other potential contaminants, 1,119,187,680 clean reads were obtained. 311,741,312 clean reads and 266,033,410 clean reads were obtained for 6-month-old WT and APP/PS1 mice respectively, while 278,938,832 clean reads and 262,474,126 clean reads were obtained for 9month-old WT and APP/PS1 mice respectively. Both reference genome and gene model annotation files were downloaded directly from the genome website. An index of the reference genome was constructed with bowtie2 v2.2.8, and paired-end clean reads were aligned to the reference genome with HISAT2 v2.0.4 [27]. HISAT2 was run with "--rna-strandness RF," while other parameters were set as default. Especially the mapping rates was 91.53% and 91.74% for APP/PS1 and WT mice, respectively. The transcripts were filtered out with coding potential prediction with the function of CNCI (Coding-Non-Coding-Index) (v2) [28], CPC (Coding Potential Calculator) (0.9-r2) [29], PfamScan (v1.3) AGING [30], and PhyloCSF (phylogenetic codon substitution frequency) (v20121028) [31]. 9,299 lncRNAs (including 2,463 annotated lncRNA and 6,836 novel lncRNA) and 48,718 protein-coding transcripts (mRNA) were identified and used for subsequent analyses.

Overview of miRNA-seq data
A total of 99,137,837 raw reads were generated. 28,672,205 raw reads and 20,760,197 were generated for 6-month-old WT and APP/PS1 mice respectively, while 23,016,993 raw reads and 26,688,442 raw reads were generated for 9-month-old WT and APP/PS1 mice respectively. After the removal of low quality and adapter sequences, 96,028,757 clean reads were obtained. 27,710,151 clean reads and 20,281,596 clean reads were obtained for 6-month-old WT and APP/PS1 mice respectively, while 22,393,001 clean reads and 25,644,009 clean reads were obtained for 9-month-old WT and APP/PS1 mice respectively. These clean reads were filtered by length (18-35 nucleotides), and a majority of the selected reads for both groups were 22 nucleotides in length. The selected reads were then mapped to the mouse reference sequence with Bowtie [32], and the mapping rate was 94.61% and 94.78% for APP/PS1 and WT mice respectively. Then, the mapped tags were then annotated and classified by alignment with noncoding small RNAs (including rRNA, tRNA, small nuclear RNA, and small nucleolar RNA), repeatassociated RNA, exon-and intron-associated RNAs in GenBank before sequenced in the miRBase v.20.0. In addition, both miREvo [33] and miRDeep2 [34] software was used to predict previously unidentified miRNAs. Ultimately, 1,411 mature miRNAs (1,312 known and 99 previously unknown) were detected and used for subsequent analyses.

AGING
(TPM) values, 32 significantly dysregulated miRNAs were identified between the 6-month-old groups, which included 9 upregulated miRNAs and 23 downregulated miRNAs in APP/PS1 mice (Figure 2D, Supplementary  Table 3). 42 miRNAs were significantly dysregulated between the 9-month-old groups, which included 20 upregulated miRNAs and 22 downregulated miRNAs in APP/PSI mice ( Figure 2E, Supplementary Table 4). The cluster analysis of miRNAs expression was performed, and then a heatmap was generated ( Figure 2F). Finally, the FPKM values (fragments per kilobase of exons per million fragments mapped) were used to estimate the expression levels of mRNA transcripts. 310 mRNA transcripts were significantly dysregulated, including 132 upregulated transcripts and 178 downregulated transcripts in APP/PS1 mice at 6 months ( Figure 2G, Supplementary Table 5), while 226 mRNAs were significantly dysregulated, with 108 and 118 upregulated and downregulated in APP/PS1 mice at 9 months ( Figure 2H, Supplementary Table 6). Once, cluster analysis for the expression of mRNAs was performed and a heatmap was generated ( Figure 2I).

qPCR validation
The differential expression identified by RNA-seq experiments were confirmed with qPCR. 24 differentially expressed transcripts were randomly selected, including 6 lncRNAs, 9 miRNAs and 9 mRNAs. All of the selected transcripts were detected in the brain of 2-9-month-old APP/PS1 and WT mice. Besides, there was statistical difference between the two groups ( Figures 3-5). Overall, the qPCR results were highly consistent with the RNA-seq data.

Construction of lncRNA-associated ceRNA networks
According to the ceRNA hypothesis, the ceRNAs can compete for the same MREs in regulatory networks. In this study, RNA-seq data were used to map ceRNA networks in the APP/PS1 brain for the first time. The differentially expressed transcripts (lncRNAs, miRNAs, and mRNAs) were split into three groups depending on the expression patterns. The 6yes9no group included transcripts differential expressed at 6 months but not differential expressed at 9 months of age, play a role in AD pathogenesis; The 6no9yes group included transcripts not differential expressed at 6 months but differential expressed at 9 months, participate in the development of AD; The 6yes9yes group included transcripts differential expressed at both 6 and 9 months, which contributed to in all stages of AD ( Figure 6A). 6yes9no group included a total of 148 lncRNAs and 376 mRNAs that were differentially expressed and shared common MRE binding sites from 33 significantly dysregulated miRNAs (Supplementary Tables 7, 8). 6no9yes significantly dysregulated group included a total of 135 lncRNAs, 526 mRNAs, and 50 miRNAs (Supplementary Tables 9, 10). Besides, 7 lncRNAs, 31 mRNAs and 2 miRNA were included in 6yes9yes group (Supplementary Table 11). The ceRNA networks included both positive and negative regulation ( Figures  6, 7). Figure 6 shows the increased lncRNAs, decreased miRNAs and increased mRNAs in APP/PS1 mice, while Figure 7 shows the decreased lncRNAs, increased miRNAs, and decreased mRNAs in APP/PS1 mice. It indicated potential critical RNA interactions involved in AD pathogenesis.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses
A lncRNA-associated ceRNA network can alter the regulation of related mRNA-encoding genes. GO analyses were performed on the genes included in the networks identified here and several GO terms were found to be significantly enriched (Supplementary  Tables 12-14). The GO terms included biological process (BP), cellular component (CC), and molecular function (MF), as shown in Figure 8. The top highly enriched terms were cytoskeleton (GO:0005856), postsynaptic density (GO:0014069), cell-cell adherens junction (GO:0005913) and dendrite (GO:0030425). A number of cognition-associated terms were also observed, such as axon (GO:0030424), synapse (GO:0045202), postsynaptic density (GO:0014069), intracellular signal transduction (GO:0035556), and neuron projection (GO:0043005). Noteworthy, the enriched GO terms for 6yes9no group were different from those for 6no9yes group, which suggested the expression changes in the functional genes during the progression of disease. For example, the Go pathways "transferase activity" was found in the "Molecular function" section of both B and C, but not in A of the Figure 8. It was explained that AD is characterized by the accumulation of intracellular and extracellular proteins, including the microtubule-associated protein Tau and the decomposition product of amyloid precursor protein β-amyloid Aβ. Tissue-type transglutaminase (tTG) is a calcium-dependent enzyme that catalyzes the cross-linking of proteins to generate isomeric peptides of the γ-glutamyl-ε-lysine structure (representing transferase activity). This covalent attachment results in protein aggregation and deposition due to its strong resistance to proteolysis. In previous studies it was demonstrated that the protein levels of tTG and isopeptides were increased in the brain of patients with advanced AD, and the activity increases noticeably with age. These findings suggest that tTG may be an important cause of abnormal protein AGING accumulation in advanced AD pathology, which may not have been changed in the early stages [35,36]. Consequently, 6yes9no group is not enriched in this pathway. In summary, the lncRNA-associated ceRNA networks might participate in the pathological prog-ression of AD at distinct stages through different mechanisms. The establishment of these networks helps to investigate the functions of the key genes in AD and to guide determination of the regulatory mechanisms between the components of the ceRNA network.

Figure 2. Expression profiles of distinct RNAs. (A-C) Expression profiles of lncRNAs. (A, B)
In the volcano plots, green, red, and black points represent lncRNAs that were downregulated, upregulated, and not significantly different in APP/PS1 mice relative to wild-type (WT) control mice at 6 and 9 months, respectively. x-axis: log2 ratio of lncRNA expression levels between AD and WT. y-axis: false-discovery rate values (-log10 transformed) of lncRNAs, P<0.05 (C) Cluster analysis of expression of lncRNAs. Red and blue: increased and decreased expression at 6 and 9 months, respectively. Expression profiles are similarly shown for (D-F) miRNAs, p<0.04 and (G-I) mRNAs, q<0.05.

Figures 3. Validation of expression of lncRNAs by using qPCR.
The identified differentially expressed transcripts (lncRNAs, miRNAs, and mRNAs) were divided into three groups. (A) 6yes9no group represents transcripts differential expressed at 6 months but not at 9 months; (B) 6no9yes group represents transcripts not differential expressed at 6 months but differential expressed at 9 months; (C) 6yes9yes group represents transcripts differential expressed at both 6 and 9 months. The expression of lncRNAs was quantified relative to Gapdh expression level by using the comparative cycle threshold (ΔCT) method. Data are presented as means ± SD (n = 3, *p < 0.05, **p < 0.01). AGING KEGG pathway analysis was conducted to determine the signaling cascades related to the identified genes. By using p < 0.05 as the threshold value, a number of significantly enriched pathways were identified ( Figure  9, Supplementary Tables 12-14), including neuroactive ligand-receptor interaction, AMP-activated protein kinase (AMPK) signaling, long-term potentiation, Hippo signaling, glutamatergic synapse, PI3K-Akt signaling, insulin secretion, focal adhesion and axon guidance.

Association study
We selectively analyzed the data for lncRNAs and miRNAs of which the target genes significant differential expressed between APP/PS1 and WT mice (corrected p < 0.05). Additionally, we also selected lncRNAs and miRNAs of which the target genes showed enrichment in the mouse brain and were associated with AD. Analyses were performed to investigate the relationships between lncRNA-associated ceRNA networks and AD. For example, LNC_000854, LNC_001450, LNC_001451, LNC_001887, LNC_002205, LNC_002746, LNC_ 003197, LNC_003206, LNC_003458, LNC_004148, LNC_004514, LNC_004707 and LNC_006482 were identified as ceRNAs of mmu-miR-122-5p, which targets Klf4. The expression of Klf4 was higher in AD mice than in WT mice. Klf4 regulates amyloid-β (Aβ)-induced neuroinflammation and plays a potential role in not only oligomeric Aβ42-induced neurotoxicity but also the pathogenesis of Alzheimer's disease [37]. LNC_000033 was found to be a ceRNA of mmu-miR-128-2-5p, AGING mmu-miR-135b-5p, mmu-miR-3097-3p, mmu-miR-31-5p, and mmu-miR-449a-5p, which target Synpo. Synpo was upregulated in the sequencing data. Synpo have been identified as crucial components in dendritic spine plasticity of the developing hippocampus. The disorder of Synpo expression affected the development of dendritic spines [38][39][40][41][42], which might also affect the early stages of AD. In addition, our analysis also revealed that miRNAs may act directly on their target genes. For example, Grsf1, targeted by mmu-miR-187-3p, mmu-miR-363-3p and mmu-miR-7004-5p, could be essential for the development of embryonic brain [43]. The additional results are listed in Supplementary  Tables 7-11. Overall, it suggested that the identified lncRNA-associated ceRNA networks may be involved in the regulation of AD. AGING DISCUSSION AD is the most common neurodegenerative disease, of which the exact pathologic mechanism remains unknown. Although the traditional research tools and expertise detect some early brain changes of AD, additional research technologies remains necessary to fine-tune the accuracy of these tools. Recently, several studies on AD have focused on the epigenetic regulation of AD pathogenesis and identified the potential targets for therapy. lncRNAs expressed in brains have been reported to be contributory to the pathophysiology of AD [44,45]. However, the roles of the lncRNAs in AD have remained mostly unknown.
MiRNAs are 22-nucleartide-long ncRNAs that can induce target gene silencing through complementary base-pairing with MREs on their 3ʹ UTRs and recruitment of RNA induced silencing complex (RISC) [46]. A total of 519 canonical miRNA genes have been identified in the human genome [47,48]. There are roughly 70% of the identified miRNAs expressed in the provisional brain and in neurons [49], which suggests that miRNAs perform critical regulatory functions in the development of central nervous system (CNS), the formation of dendritic spine, neurite outgrowth, as well as neuronal differentiation and maintenance. The deregulation of miRNA is involved in neurodegenerative disorders such as AD and Parkinson's disease (PD) [6], as well as in psychiatric disorders such as schizophrenia [50]. For example, let-7, miR-15a and miR-101 target APP, while miR-15a, miR-9 and miR-107 regulate BACE1 [8].
Over the past several years, the ceRNA hypothesis has been validated by numerous experiments. So far, ceRNA mechanisms and network construction have been mainly studied in the field of cancer research [51][52][53][54]. Only a small number of ceRNA interactions have been reported to be associated with neurodegenerative disorders. Recently, researchers have started to explore the ceRNA regulatory mechanism for specific neurodegenerative disorders in a systematic manner. Nevertheless, exciting advances have been made in our understanding of ceRNA interactions in neurodegenerative disorders. In our previous study, the lncRNA-ceRNA network was constructed in the brain AGING of 12-month-old APP/PS1 mice [55]. As the sequencing depth and coverage limit the amount of generated sequences, the compositions of networks are slightly different. However, the key pathways in AD have been discovered. In order to not only explore changes of the ceRNA networks in the early stages of AD but also identify the early biomarkers, lncRNA-ceRNAs networks were constructed in the brain of 6-and 9month-old APP/PS1 mice.
So far, it is the first comprehensive high-throughput sequencing analysis of the expression profiles of lncRNA, miRNA and mRNA in the APP/PS1 mouse model of AD. The dysregulated lncRNAs, miRNAs and mRNAs exhibited significant differential expression between AD and WT control groups, which suggested that these transcripts are associated with the pathogenesis of AD. For example, Fndc3b is essential for proliferation, adhesion, spreading and migration of nerve cells [56]. Trappc9, another differential expressed gene, plays a critical role in the development of human brain, possibly through its effect on NF-kappaB activation and protein trafficking in the postmitotic neurons of the cerebral cortex [57]. Acsl6 mRNA is highly enriched in the brain and acsl6-/-mice demonstrate motor impairments, altered glutamate metabolism, increased astrogliosis and microglia activation [58]. Mmu-miR-376a-5p is a pancreatic isletspecific miRNA that regulates insulin secretion as an important pathway for the development of AD [59], while mmu-miR-134 regulates the development of cortical neurons [60]. The level of mmu-miR-29 is also substantially reduced in AD patients and it acts to regulate BACE1 expression [61,62]. The qRT-PCR experiment confirmed the profiles from the highthroughput sequencing data, which indicated the reliability of the sequencing data.
In our study, RNA-seq was used to systematically analyze lncRNA, miRNA and mRNA profiles in the brain of 6and 9-month-old APP/PS1 mice. The transcripts of 6yes9no group might participate in AD pathogenesis, of which the stability and specific expression could make them suitable as optimal biomarkers for AD. The transcripts in 6no9yes group might function in the development of AD. It is worth noting that the transcripts in 6yes9yes group may be involved in the disease at all stages, which suggests that focus on these transcripts could facilitate the development of lncRNA-based diagnostic tools and therapeutic strategies for AD. Overall, lncRNA and miRNA molecules have a potential to act as key regulators AD in different aspects. Both lncRNAs and protein-coding mRNAs function as ceRNAs and supersponges to regulate the expression of miRNA. Therefore, we predicted the miRNA-mRNA and miRNA-lncRNA interaction with miRanda and constructed DElncRNA-DEmiRNA-DEmRNA triple networks for APP/PS1 and WT mouse brain. The selected lncRNA-associated ceRNA networks may facilitate new insights into AD and contribute novel treatments for the disease.
We performed GO enrichment and KEGG analysis of the genes in the ceRNA networks and identified not only a number of enriched terms relevant to the pathological process of AD, including cytoskeleton (GO:0005856), cell adhesion (GO:0005913), dendrite (GO:0030425), postsynaptic density (GO:0014069), axon (GO:0030424), synapse (GO:0045202), and neuron projection (GO:0043005), but also various pathways, including MAPK signaling, insulin secretion, Type II diabetes mellitus, cAMP signaling, Hippo signaling, focal adhesion, dopaminergic synapse, and PI3K-Akt signaling pathways. Analysis of the data revealed several lncRNA-associated ceRNA networks that participate in AD. Akap5 contributes to synaptic plasticity mediated by NMDARs and AMPA-type glutamate receptors (AMPARs) and plays a critical role in the progression of AD [21,63,64]. One of these networks involves the gene Akap5 and the ceRNAs, including LNC_000217, LNC_000233, LNC_000622, LNC_001498, LNC_001502, LNC_001818, LNC_ 002144, LNC_002373, LNC_002451, LNC_002620, LNC_003852, LNC_004317, LNC_005072, and LNC_005613. These ceRNAs have the potential to bind mmu-miR-679-5p, which targets Akap5. The APP/PS1 mice used in this study cannot represent the whole disease, which is primarily related to β-amyloid toxicity. Therefore, further research is necessary to better understand the regulation of these networks in AD.
Our research is just the beginning so that there remain many challenges and problems to be solved in the future. The cumulative evidence has helped to refine the dynamic ceRNA regulation [20], which reveals that various factors could contribute to the creation of miRNA and ceRNA hierarchies. These factors can be summarized as follows: miRNA target-site efficacy; shared MRE abundance; miRNA/ceRNA expression level and subcellular localization [21]; miRNA: target ratio; competition between rate-limiting molecules, such as Ago [63]; and advanced ceRNA hierarchy strategies. Moreover, attempt was made to decode the complexity of ceRNA regulatory networks. After the initial establishment of the ceRNA network, the biological mechanisms and functions mediated by the ceRNA mechanism should be predicted and verified by means of in vivo experiments in the future.
For a long time, there has been widespread recognition that ncRNAs are incapable to encode proteins [64]. Nevertheless, with the advancement of deep ribosome profiling sequencing (Ribo-Seq) technology, mass AGING spectrometry and algorithms, a subset of ncRNA have been identified to encode peptides (<100 amino acids) or proteins, such as muscle-specific lncRNAs [65][66][67] and cancer-related lncRNA HOXB-AS3 [68]. At present, only few studies focused on Alzheimer's disease. Peptides/ proteins encoded by ncRNAs might represent the drug targets or biomarkers for the prognosis of AD patients. Therefore, it's significant to summarize the characteristics of peptides/proteins encoded by ncRNAs and their outlook for small molecule peptide drugs, drug targets and biomarkers.
In conclusion, the brain lncRNA-associated ceRNA profiles of APP/PS1 and WT mice were clarified. Our findings improve the current understanding of ceRNA biology and the regulatory roles of these RNAs in the pathogenesis of AD. These new networks reveal the potential biomarkers and may offer a promising target for the development of drugs to treat AD.

Tissue preparation
WT and APP/PS1 mice [originally from The Jackson Laboratory; strain B6.Cg-Tg(APPswe, PSEN1dE9) 85Dbo/Mmjax [43]] were purchased from the Model Animal Research Center of Nanjing University. The mice were housed one per cage under standard conditions (25°C, 50% humidity, a 12-hour light/dark cycle, and specific-pathogen-free environment). The mice were provided with free access to the standard diet until they met the age requirements (6 and 9 months). 3 AD and 3 wild-type male mice born in the same litter were used as experimental and control groups, of which cerebral cortex samples were collected for RNA-seq. All animal experiments were performed in accordance with animal use protocols approved by the Committee for the Ethics of Animal Experiments, Shenzhen Peking University, the Hong Kong University of Science and Technology Medical Center (SPHMC) (protocol number 2011-004).

RNA extraction and qualification
Total RNA from each sample was isolated with TRIzol reagent (Invitrogen) according to the manufacturer's instructions and separated on 1% agarose gels to assess RNA degradation and contamination. RNA purity was measured with a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured with a Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was evaluated with the RNA Nano 6000 Assay Kit of a Bioanalyzer 2100 System (Agilent Technologies, CA, USA).

RNA-seq
Details of the mRNA-seq, miRNA-seq, and lncRNAseq methods are described in Supplementary Materials.

Expression analysis
We calculated the FPKM values of transcripts by using Cuffdiff (v.2.1.1) to evaluate the expression levels of protein-coding genes and lncRNA in each sample [69]. The expression levels of miRNAs were estimated as TPM values as described [69]. Transcripts with p values less than 0.05 were regarded as being differentially expressed between APP/PS1 and WT mice. Normalized expression = (mapped reads)/(total reads) * 1,000,000.

CeRNA network analysis
The expression levels of lncRNAs, miRNAs, and mRNAs differed significantly between APP/PS1 and WT mice. We searched the sequences of the lncRNAs and mRNAs to identify potential MREs. We used miRanda (http://www.microrna.org/microrna/) to predict miRNA-binding seed-sequence sites, and the presence of the same miRNA-binding sites in both lncRNAs and mRNAs indicated potential lncRNA-miRNA-mRNA interaction.

GO annotations and KEGG pathway analyses
The DAVID (https://david.ncifcrf.gov/summary.jsp) database was used to analyze lncRNA-miRNA-enriched genes. GO and KEGG terms with p values less than 0.05 were considered significantly enriched.

Construction of lncRNA-associated ceRNA networks
The lncRNA-associated ceRNA networks were constructed and visually displayed by using Cytoscape software V3.5.0 (San Diego, CA, USA) based on the analysis of high-throughput sequencing data, as described above. In the figures, distinct shapes and colors are used to represent different RNA types, and regulatory relationships.

Real-time qPCR validation
Total RNA was extracted by using TRIzol reagent (Sigma) according to the manufacturer's protocol. RNA quantity was measured by using a NanoDrop 2000 (Thermo Fisher Scientific). Quantitative RT-PCR was performed by using the GoScript TM Reverse Transcription System (Promega), in a C1000 Thermal Cycler (Bio-Rad). The glyceraldehyde-3-phosphate dehydrogenase gene (Gapdh) and U6 were used as an AGING internal control. Relative gene-expression levels were calculated using the 2 −ΔCt method (n=3).

Statistical analysis
Two normally distributed groups were compared by using t tests. Parameters for the high-throughput sequencing-related data were calculated, and statistical computing was performed by using R software. All data are expressed as means ± SD; a value of p < 0.05 was considered statistically significant.

Data access
All raw and processed sequencing data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm. nih.gov/geo/) under accession number GSE132177.

AUTHOR CONTRIBUTIONS
JW, WZ, and NM designed the study. NM, CT and BY participated in the animal experiments including tissue collection and RNA/protein extraction. NM performed the experiments and analyzed the data. NM, JW, and WZ analyzed the results and wrote the manuscript. All authors read and approved the final manuscript. AGING

RNA isolation, quantification and qualification
RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the Nano Photometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for small RNA sequencing
A total amount of 3 μg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, NEB 3' SR Adaptor was directly, and specifically ligated to 3' end of miRNA, siRNA and piRNA. After the 3' ligation reaction, the SR RT Primer hybridized to the excess of 3' SR Adaptor (that remained free after the 3' ligation reaction) and transformed the single-stranded DNA adaptor into a double-stranded DNA molecule. This step is important to prevent adaptor-dimer formation, besides, dsDNAs are not substrates for ligation mediated by T4 RNA Ligase 1 and therefore do not ligate to the 5´ SR Adaptor in the subsequent ligation step. 5´ends adapter was ligated to 5´ends of miRNAs, siRNA and piRNA. Then first strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H-). PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for illumina and index (X) primer. PCR products were purified on a 8% polyacrylamide gel (100V, 80 min). DNA fragments corresponding to 140~160 bp (the length of small noncoding RNA plus the 3' and 5' adaptors) were recovered and dissolved in 8 μL elution buffer. At last, library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips.

Clustering and sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50bp.

Quality control
Raw data (raw reads) of fastq format were firstly processed through custom perl and python scripts. In this step, clean datas(clean reads) were obtained by removing reads containing ploy-N, with 5′ adapter contaminants, without 3′ adapter or the insert tag, containing ploy A or T or G or C and low quality reads from raw data. At the same time, Q20, Q30, and GCcontent of the raw datas were calculated. Then, chose a certain range of length from clean reads to do all the downstream analyses.

Reads mapping to the reference sequence
The small RNA tags were mapped to reference sequence by Bowtie (Langmead et al. 2009) without mismatch to analyze their expression and distribution on the reference.

Known miRNA alignment
Mapped small RNA tags were used to looking for known miRNA. miRBase20.0 was used as reference, modified software mirdeep2 (Friedlander et al. 2011) and srna-tools-cli were used to obtain the potential miRNA and draw the secondary structures. Custom scripts were used to obtain the miRNA counts as well as base bias on the first position of identified miRNA with certain length and on each position of all identified miRNA respectively.

Remove tags from these sources
To remove tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker, Rfam database or those types of datas from the specified species itself.

Novel miRNA prediction
The characteristics of hairpin structure of miRNA precursor can be used to predict novel miRNA. The available software miREvo ) and mirdeep2 (Friedlander et al. 2011) were integrated to predict novel miRNA through exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps. At the same time, custom scripts were used to obtain the identified miRNA counts as well as base bias on the first position with certain length and on each position of all identified miRNA respectively.

Small RNA annotation summary
Summarizing all alignments and annotations obtained before. In the alignment and annotation before, some small RNA tags may be mapped to more than one category. To make every unique small RNA mapped to only one annotation, we follow the following priority rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > NAT-siRNA > gene > novel miRNA > ta-siRNA. The total rRNA proportion was used a marker as sample quality indicator. Usually it should be less than 60% in plant samples and 40% in animal samples as high quality.

miRNA editing analysis
Position 2~8 of a mature miRNA were called seed region which were highly conserved. The target of a miRNA might be different with the changing of nucleotides in this region. In our analysis pipeline, miRNA which might have base edit could be detected by aligning all the sRNA tags to mature miRNA, allowing one mismatch.

miRNA family analysis
Exploring the occurrence of miRNA families identified from the samples in other species. In our analysis pipeline, known miRNA used miFam.dat (http://www.mirbase.org/ftp.shtml) to look for families; novel miRNA precursor was submitted to Rfam (http://rfam.sanger.ac.uk/search/) to look for Rfam families.

Target gene prediction
Predicting the target gene of miRNA was performed by miRanda (Enright et al, 2003) for animals.

Quantification of miRNA
miRNA expression levels were estimated by TPM (transcript per million) through the following criteria (Zhou et al. 2010): Normalization formula: Normalized expression = mapped readcount/Total reads*1000000

Differential expression of miRNA
For the samples with biological replicates: Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.8.3). The Pvalues was adjusted using the Benjamini and Hochberg method. Corrected P-value of 0.05 was set as the threshold for significantly differential expression by default.

RNA isolation, quantification and qualification
RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the Nano Photometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for lncRNA sequencing
A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Firstly, ribosomal RNA was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA), and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations. Briefly, fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase(RNaseH-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with sizeselected, adaptor-ligated cDNA at 37° C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using AGING TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated.

Quality control
Raw data(raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data(clean reads) were obtained by removing reads containing adapter, reads on containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content of the clean data were calculated. All the down stream analyses were based on the clean data with high quality.

Mapping to the reference genome
Reference genome and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using bowtie2 v2.2.8 and paired-end clean reads were aligned to the reference genome using HISAT2 (Langmead et al.) v2.0.4. HISAT2 was run with '--rna-strandness RF', other parameters were set as default.

Transcriptome assembly
The mapped reads of each sample were assembled by StringTie (v1.3.1) (Mihaela Pertea et al. 2016) in a reference-based approach. StringTie uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus.

Coding potential analysis
CNCI CNCI (Coding-Non-Coding-Index) (v2) profiles adjoining nucleotide triplets to effectively distinguish protein-coding and non-coding sequences independent of known annotations (Sun et al. 2013). We use CNCI with default parameters.

CPC
CPC (Coding Potential Calculator) (0.9-r2) mainly through assess the extent and quality of the ORF in a transcript and search the sequences with known protein sequence database to clarify the coding and non-coding transcripts (Kong et al. 2007). We used the NCBI eukaryotes' protein database and set the e-value '1e-10' in our analysis.

PhyloCSF
PhyloCSF (phylogenetic codon substitution frequency) (v20121028) examines evolutionary signatures characteristic to alignments of conserved coding regions, such as the high frequencies of synonymous codon substitutions and conservative amino acid substitutions, and the low frequencies of other missense and nonsense substitutions to distinguish protein-coding and non-coding transcripts (Lin et al. 2011). We build multispecies genome sequence alignments and run phyloCSF with default parameters. Transcripts predicted with coding potential by either/all of the four tools above were filtered out, and those without coding potential were our candidate set of lncRNAs.

Conservative analysis
Phast (v1.3) is a software package contains much of statistical programs, most used in phylogenetic analysis (Siepel, et al. 2005), and phastCons is a conservation scoring and identificating program of conserved elements. We used phyloFit to compute phylogenetic models for conserved and non-conserved regions among species and then gave the model and HMM transition parameters to phyloP to compute a set of conservation scores of lncRNA and coding genes.

Quantification of gene expression level
Cuffdiff (v2.1.1) was used to calculate FPKMs of both lncRNAs and coding genes in each sample . Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. FPKM means fragments per kilo-base of exon per million fragments mapped, calculated based on the length of the fragments and reads count mapped to this fragment.

Differential expression analysis
The Ballgown suite includes functions for interactive exploration of the transcriptome assembly, visualization of transcript structures and feature-specific abundances for each locus, and post-hoc annotation of assembled features to annotated features (Alyssa C. Frazee et al. 2014). Transcripts with an P-adjust <0.05 were assigned as differentially expressed. Cuffdiff provides statistical AGING routines for determining differential expression in digital transcript or gene expression data using a model based on the negative binomial distribution . Transcripts with an P-adjust <0.05 were assigned as differentially expressed. AGING

Supplementary Tables
Please browse Full Text version to see the data of Supplementary Tables 1, 2 Supplementary