We searched the GEO database and the Affymetrix Human Genome U133 Plus 2.0 Array chip platform to identify gene expression profile data using the keyword “Alzheimer's disease.” Samples sizes with ≥30 datasets and 7 datasets were obtained. In addition, the dataset GSE16759, which simultaneously detected the expression profiles of miRNA and mRNA, included eight sets of 629 samples. Forty samples of undefined AD were excluded, and finally 589 samples were included as shown in Table 1. For the expression profile data, the raw expression level data (.CEL files), for each sample was downloaded and normalized using the robust multi-array (RMA) of the R package affy (Gautier et al., 2004). Next, the batch effect was removed using the combat of the R package SVA (Leek et al., 2012). The microarray data were derived from the GEO (accession GSE133349). We matched the probes with their corresponding genes. In cases where multiple probes corresponded to one gene, the median was expressed. Where one probe corresponded to multiple genes, such a probe was excluded.

To obtain the latest lncRNA expression profile, lncRNA re-annotation was performed using the HG-U133 2.0 array. First, the GENCODE Release 21 version of lncRNA transcript sequences was downloaded from GENCODE. Then, probe clusters sequences were aligned to lncRNA sequences using seqmap (Hawkins et al., 2011), and mismatch was set to 0. At least 11 clusters were selected and compared to probes on the same lncRNA to confirm the success of re-annotation. Finally, 3,943 Affymetrix probe sets (3,158 lncRNA) were included in subsequent analyses.

StarBase V3.0 database (Li et al., 2014) is designed to integrate large-scale CLIP-Seq (HITS-CLIP, PAR-CLIP, iCLIP, and CLASH) data to decode interactive networks. Five predictive algorithms were used to predict miRNA target genes, including TargetScan, miRanda, Pictar, PITA, and RNA22. Data of 1,055,319 miRNA-mRNA interactions for 484 miRNAs and 15,064 mRNAs and of 63,698 miRNA-lncRNA interactions for 642 miRNAs and 3,789 lncRNAs were downloaded from the starBase V3.0 database.

As described by Fan et al. (2018), the hypergeometric model was used to construct a ceRNA network. In summary, the miRNAs shared by mRNA and lncRNA were used to calculate the interactive probability of mRNA-lncRNA networks. The threshold was FDR < 0.05. When more than three miRNAs were found to be shared between mRNA and lncRNA, they were selected. The value of p was calculated using the following formula:

where “m” represents the total number of miRNAs in starBase database, “t” represents the number of miRNAs that interacted with the mRNA, “n” represents the number of miRNAs that interacted with the lncRNA, and “r” represents the number of miRNAs shared between mRNA and lncRNA.

Correlations among mRNA-lncRNA pairs were calculated in the expression profiles of normal and AD samples, and the correlation coefficient > 0.25 and FDR < 0.05 were selected to identify the ceRNA global network.

We compared the global networks of ceRNAs in normal and AD samples and selected the intersection of two networks. Furthermore, the difference between the Pearson's correlation coefficients of mRNA-lncRNA in two ceRNA networks was calculated. With 1,000 random readings as the statistical background, the probability (p ≤ 0.05) of Pearson's correlation coefficient difference was applied to identify dysregulated mRNA-lncRNA.

Each sample was subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of ssGSEA using the R package GSVA38 to identify pathways associated with lncRNA expression, applying a correlation > 0.5.

To facilitate the reproduction of our results, the eight datasets were uploaded to the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi), with the ID: GSE133349.

TRIzol reagent (ThermoFisher Scientific, 15596026) was used to extract total RNA from serum samples. The synthesis of total RNA into cDNA was performed according to the instructions of the Reverse Transcription Kit (TshermoFisher Scientific, BTK1622). The ABI 7500 real-time PCR instrument was used to perform the real-time PCR reaction using the quant one-step real-time quantitative reverse transcription PCR (qRT-PCR) Kit [FP303, Tiangen Biochemical Technology (Beijing) Co., Ltd.]. The primers used were as follows: SNHG14 forward, ATGAGCTGACAACCTACTCC and reverse, AAGTCATCTTCTGCAAGGGT; PART1 forward, CCCTTTCACTATGAAGGACC and reverse, ATTTACCCGTCCAGTTCTG; NNT-AS1 forward, AGAAACAGGTCTAAAGACCCT and reverse, TTCTTGGCATCTCTGAGCA; AC093010.3 forward, TGATGTGTTTGCTATCTGCT and reverse, TGTTAACAGCTAGCCATTCAG; ARMCX5-GPRASP2 forward, AAGAGAAGGGATAGAGTGGTG and reverse, CTTCTGTCATAGAAATTTCCCTCTC; GAPDH forward, TCAAGATCATCAGCAATGCC and reverse, CGATACCAAAGTTGTCATGGA; GAPDH were used as internal controls. 2^−ΔΔCT was used to calculate the relative expression.

Statistical analyses were performed using R 3.4.3, and all analyses (non-specific description) were performed with default parameters. Group comparisons were performed using the Student's t-test, correlations were determined using Pearson's correlation coefficient, g:profiler (Reimand et al., 2019) was used for gene enrichment analysis, and Cytoscape (Shannon et al., 2003) (http://www.cytoscape.org/) was used for network visualization.

To investigate the potential role of miRNA-mediated ceRNA networks in AD, we analyzed the ceRNA landscape in normal samples and AD samples based on the expression profiles (Figure 1A). Using the hypergeometry test, we identified 164,494 miRNA-mediated mRNA-lncRNA pairs with 7,484 genes. To construct a sample-specific ceRNA network, the correlation of these mRNA-lncRNA pairs in the expression profiles of normal and AD samples was determined. The competitive endogenous RNAs (ceRNAs) network was established using the correlation coefficient > 0.25 and FDR < 0.05. A total of 2,571 mRNA-lncRNA pairs and of 255 lncRNAs for 2,551 genes were identified in normal samples (Figure 1B), while 3,353 mRNA-lncRNA pairs and of 227 lncRNAs for 1,619 genes were identified in AD samples (Figure 1C). The two ceRNA network distributions exhibited a consistent law distribution (Figure 1D) and were consistent with biological network characteristics, indicating that ceRNA plays an important role in AD that is similar to other biological regulatory networks.

To explore the co-expression relationship of the mRNA-lncRNAs in the ceRNA network, we analyzed the relationship between network degree and mRNA-lncRNA co-expression. The network degree was significantly positively correlated with mRNA-lncRNA co-expression (Figures 2A,B). Genes with a high network degree showed higher co-expression relationships, and many miRNAs were found to be shared by mRNA and lncRNA. A high correlation coefficient was found between the corresponding mRNA and lncRNA (Figures 2C,D). Studies have shown that the expression levels of miRNA determine the activity of ceRNA networks (Ala et al., 2013). We, therefore, studied the effect of miRNA expression on ceRNA networks in AD. Dicer is a key enzyme that regulates miRNA processing and maturation. Herein, samples were divided into two groups based on the expression level of Dicer: Dicer-low and Dicer-high. In normal samples, the co-expression pattern in the Dicer-high group was higher than that in the Dicer-low group, while opposite results were obtained in AD samples (Figures 2E,F), suggesting that moderate expression of miRNAs in AD is necessary to maintain the ceRNA network, but high expression of miRNAs may lead to degradation of many genes.

Comparison of the ceRNA networks in normal and AD samples showed that the number of ceRNAs in AD is only 53.5% to that of normal samples, suggesting a large number of ceRNA regulation abnormalities in AD samples, compared to normal samples. Notably, 2,575 (41.1%) ceRNAs were unchanged while 3,696 (58.9%) ceRNAs disappeared from the AD samples (Figure 2G). Analysis of genes and lncRNAs showed that 1,489 (58.4%) genes were shared between the two networks while 1,062 (41.6%) genes disappeared from the ceRNA network of AD (Figure 2H). These findings imply a high dysregulation of ceRNA in AD. Furthermore, we randomly divided the samples into two queues, the training queue and the verification queue, to determine the dysregulation of ceRNA networks. In the training queue, Pearson's correlation coefficients of mRNAs in two ceRNA networks were calculated, and those with p ≤ 0.05 were selected as the background of the random network. A total of 294 dysregulated mRNAs-lncRNAs were identified in training queue, and a total of 333 dysregulated mRNA-lncRNAs were identified in the verification cohort by the same method. In addition, 258 dysregulated mRNA-lncRNAs were selected from all samples by the same method. By comparing the dysregulated mRNA-lncRNAs in the three cohorts, a high degree of overlap was found between them (Supplementary Figure 1A). This showed the stability of dysregulated mRNA-lncRNAs in different samples. Finally, in all samples, we identified 255 dysregulated mRNA-lncRNAs. Among them, 103 pairs were downregulated in AD, including 96 genes and 30 lncRNAs, while 152 pairs were upregulated in AD, including 146 genes and 53 lncRNAs (Supplementary Figure 1B). Notably, 71 lncRNAs in the statistical network were co-expressed with other mRNAs, with an average proportion of the same type accounting for 95.1%. The co-expression of lncRNA and other mRNAs tended to be the same type of aggregation, suggesting that lncRNA expression levels play a key role in the dysregulation of ceRNA networks in AD.

To elucidate on the functional implications of the ceRNA network, functional enrichment analysis was performed on all 231 genes in the ceRNA network. Results showed that these genes were enriched in two gene ontology (GO) terms and were not associated with any pathway (Figure 3A), implying that the dysfunctional ceRNA network lacked specific functions to be considered a large module. These genes were enriched in 24 transcription factors and 17 miRNAs, indicating that they were highly shared between transcription factors and miRNAs. The most significantly enriched miRNA, hsa-miR-106b-5p, was found to be associated with AD susceptibility and is a potential peripheral blood biomarker for AD (Yilmaz et al., 2016; Figure 3B). The most significantly enriched transcription factor, M038071 (SP2), is an important marker in cancer chemotherapy (Vizcaino et al., 2015; Figure 3C). Thus, ceRNA is used as a small module in AD instead of a single or large module.

We further used the paired samples of miRNA, mRNA, and lncRNA expression profiles from the GSE16759 dataset to identify dysregulated mRNA-miRNA and miRNA-lncRNA pairs in AD samples. Using p < 0.05 as the cut-off value for selection, 86 pairs of miRNA-lncRNA and 111 pairs of mRNA-miRNA, including 19 miRNAs, 131 genes, and 37 lncRNAs, were identified. Interestingly, unlike mRNA-lncRNA, these dysregulated miRNA-target pairs were downregulated in AD samples (Figure 4A), and this phenomenon is consistent with negative regulation for miRNA and mRNA. The number of miRNAs shared by mRNA-lncRNA in the ternary network revealed that 71.4% of mRNA-lncRNA pairs contain only one dysregulated miRNA (Figure 4B). We constructed a ternary interaction network in which more than one mRNA-lncRNA pair of dysregulated miRNAs was shared for 40 pairs of mRNA-lncRNA, which contained 39 genes, 18 lncRNAs, and 17 miRNAs (Figure 4C). Notably, three of the four miRNAs in ENSG00000162437 (RAVER2) and ENSG00000212978 (LOC339803) pairs that shared the highest number of dysregulated miRNAs were significantly downregulated in AD (p < 0.05; Supplementary Figure 2), indicating that the shared dysregulated miRNAs were more likely to be related to ceRNAs of mRNA-lncRNA.

To identify AD-related lncRNAs, we analyzed the differential expression of 18 lncRNAs in AD from more than one mRNA-lncRNA pair sharing dysregulated miRNAs. We found nine (52.9%) lncRNAs to be significantly differentially expressed (Supplementary Figure 3). Similarly, we examined the differential expression of 39 genes in AD and found that 32 (82.1%) genes were significantly differentially expressed (Supplementary Figure 4), suggesting that the constructed dysregulated ceRNA network could effectively screen for disease markers. There were 24 differential pairs with simultaneous differential differences of mRNA and lncRNA. These 24 pairs had significant correlations. Among them, the average correlation coefficients of ENSG00000117242 and ENSG00000224078 and their corresponding five genes ranged between 0.48 and 0.27, while the average correlation coefficient of ENSG00000248092 and ENSG00000259976 and their respective four genes was between 0.3 and 0.28 (Figure 5A). In summary, the lncRNA-mRNA pairs that were found to be dysregulated in AD were significantly co-expressed, while 9 lncRNAs and 24 mRNAs may be potential biomarkers in AD. Most of these 24 genes have been reported to be associated with AD, such as low-density lipoprotein receptor (LDLR) overexpression which increases β-amyloid clearance and decreases amyloid deposition (Krishnan et al., 2020), HSPA8 is downregulated in multiple brain regions in AD (Silva et al., 2014), and SLC9A6 mutation causes intellectual disability (Garbern et al., 2010). These nine lncRNAs have not been reported to be associated with AD, and based on this, we further analyzed them. Moreover, KEGG pathway analysis was performed on the samples using ssGSEA to determine the pathways associated with the nine lncRNAs. KEGG pathways with a correlation exceeding 0.5 were selected and enriched in 43 pathways. These pathways were divided into two categories, one of which is positively correlated with the highly expressed lncRNAs in AD, the other category is positively related to lncRNAs, which are lowly expressed in AD. Pathways associated with downregulated lncRNAs were mainly those related to senile disease-related pathways and tumor, as well as immune-related pathways. These pathways include AD, Parkinson's disease, and type II diabetes mellitus (Figure 5B). These pathways play important roles in the development of AD, and the nine lncRNAs regulate the pathogenesis of AD.

To investigate the diagnostic value of lncRNA in AD, a receiver operating characteristic (ROC) curve classification was performed for the nine lncRNAs in AD samples and normal samples. Results showed that each of the nine lncRNAs could diagnose AD with an average area under the curve (AUC) of 0.6 (Figure 6A). However, the diagnostic value of single molecules may be limited. Thus, we randomly combined nine genes to obtain 502 combinations and used these combined expression profiles to establish support vector machine (SVM) models for the diagnosis of AD. We found that the accuracy of the lncRNA combinations gradually increased as the number of lncRNAs increased (Figure 6B). The combination with the highest accuracy contained five lncRNAs: ENSG00000224078 (SNHG14), ENSG00000152931 (PART1), ENSG00000248092 (NNT-AS1), ENSG00000259976 (AC093010.3), and ENSG00000271147 (ARMCX5-GPRASP2). The SVM was used to establish a diagnostic model, which was tested using a 10-fold cross-validation method. This test showed that the classification accuracy of the diagnostic model was 69%; 409 of the 589 samples were correctly classified. The sensitivity of the diagnostic model for AD was 71.3% while its specificity was 68%. We further used the GSE5281 dataset to verify the accuracy of the diagnostic model, with results showing a classification accuracy rate of 78.3%; 126 of the 161 samples were correctly classified. The sensitivity of the model for AD was 77% while the specificity was 79.7% (Figure 6C). These results indicate that the diagnostic model could effectively distinguish between patients with AD and normal control populations. Five lncRNAs showed a good potential to be reliable biomarkers for the diagnosis of AD. To confirm the role of five lncRNAs, qRT-PCR was used to detect the expression levels of lncRNA in the serum samples of patients with AD. PART1 and ARMCX5-GPRASP2 expressions were found to have similar trend expressions with bioinformatics analysis (Supplementary Figure 3). However, only PART1 was downregulated while SNHG14 was upregulated in the serum samples of patients with AD when compared to normal samples (Figure 7).