For bioinformatics analysis, raw sequencing data of RNA-seq and relative metadata were obtained from Gene Expression Omnibu (GEO) database (GSE104704, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104704), and RNA expression levels of genes in the lateral temporal lobe of the patients were sequenced by Illumina NextSeq 500 platform. In addition, the diagnosis information of participants from two groups, elderly control (Control) and AD, were obtained for further research. For clinical experiments, patients were diagnosed with AD in our hospital were enrolled [26]. Controls refer to participants who were examined due to transient neurologic symptoms. Cerebrospinal fluid was drawn through lumbar puncture from the space between L3-L4 or L4-L5 and stored at –80 °C. The study protocol was approved by the Medical Ethics Committee of Red Cross Hospital (Approval NO: 202280).

RNA-sequencing data were downloaded from the GEO database and converted to FASTQ files by using SRAtoolkit software (version 3.0.1, National Institutes of Health, Bethesda, MD, USA). Gene expression data were directly downloaded from supplemental data in GSE104704. CircRNA expression levels were predicted with CIRCexplorer software (version 2.0, CIRCexplorer, Shanghai, China). Briefly, FASTQ files were mapped to human genome HG19 with STAR aligner. The alignment process produced BAM files and Chimeric junction files. Chimeric junction files were parsed by using the CIRCexplorer2 parse module. Then, the Browser Extensible Data (BED) file containing the subsequent connection information is created. BED files were annotated with reference genome HG19, and results containing circRNA expression and annotation files were saved.

Differential expression analysis and normalization of circRNA expression data were processed using R package edgeR. Results were visualized with a volcano plot. CircRNAs with low expression were removed. Expression levels of filtered circRNAs were then visualized using a boxplot.

The ROC curve was used to evaluate the diagnostic value of differential expressed circRNAs. Briefly, 12 ROC curves were plotted using GraphPad Prism (version 8.0, GraphPad Software, San Diego, CA, USA), circRNA with the highest AUC was thought to have crucial diagnostic value for AD.

Gene Set Enrichment Analysis (GSEA) and Metascape pathway enrichment analysis were processed separately. The pathways affected by circRNA signature were searched using GSEA (version 4.3.2, Broad institute, Cambridge, MA, USA). Patients with low and high circ_0001535 expression were divided into two groups by median. Then, GSEA was processed to screen for pathway activity difference. Pathways with p value 0.05 were considered significantly enriched. After screening for key pathways affected by circ_0001535, Metascape was used to analyze the relationship between mRNA and circRNA. Briefly, Pearson’s correlation score between each mRNA and circ_0001535 was calculated. Gene which significantly correlated to circ_0001535 expression was obtained and listed on the Metascape website (https://metascape.org/gp/index.html#/main/step1). The results were downloaded from the webpage of the Metascape website.

The binding motif sequence between transcription factors and targeted genes was predicted using the JASPAR webpage (http://jaspar.genereg.net/). The gene symbol of the transcription factor was first searched in the JASPAR webpage and selected for use. The transcription start site (TSS) region of targeted mRNA was downloaded from National Center for Biotechnology Information (NCBI) database and uploaded to the JASPAR database. Binding sites between transcription factor and mRNA transcription start site (TSS) region were predicted automatically.

Human neuroblastoma cells SH-SY5Y were purchased (BeNa Culture Collection, Beijing, China). SH-SY5Y cell line was authenticated at Meisen Biotech Co., Ltd. by short tandem repeat (STR) profiling. No contamination of mycoplasma has been identified by the company. SH-SY5Y cells were cultured in the Dulbecco’s modified Eagle’s medium (DMEM) (Invitrogen, Carlsbad, CA, USA) including 10% heat-inactivated fetal calf serum (FBS) and 1% Penicillin/streptomycin (P/S). Cells were cultured at 37 °C in humidified 5% CO${}_2$. In order to construct an AD cell model, cells were treated with 20 µmol/L A$\beta$25-35 (HY-P0128, MedChemExpress, Monmouth Junction, NJ, USA). Overexpression (oe)-NC and oe-circ_0001535 were all purchased from RiboBio in Guangzhou. Small interfering (si)-E2F1, si-DHFR and si-NC were all purchased from RiboBio in Guangzhou. Cell transfection was done using Lipofectamine${}^{\mathrm{\circledR }}$ 2000 (11668019, Thermo Fisher Scientific, Waltham, MA, USA).

ChIP assays were performed using a ChIP kit (ab500, Abcam, Burlingame, CA, USA), following according to the manufacturer’s instructions. Briefly, SH-SY5Y cells were harvested and sonicated to chromatin. The anti-E2F1 antibody (ab245308, Abcam, Burlingame, CA, USA) was added to form the antibody-target protein-DNA complex, and protein A-Sepharose beads were used to immunoprecipitate the complex. Reverse transcription‑quantitative PCR (RT-qPCR) was performed to detect the binding site. Primers used in ChIP-qPCR were as follows, Site 1 Primer: forward: 5’-GGTGAGTTGTGGGGGATTCT-3’ and reverse: 5’- CCATCACCTATAGGGGGCCA-3’; Site 2 Primer: forward: 5’-TACGTCAGGCCTTCTCAGAGT-3’ and reverse: 5’-ATCCCCCACAACTCACCAGA-3’. Magna RIP kit (Millipore, Billerica, MA, USA) was used to process the RIP assay. Briefly, the magnetic beads were incubated with anti-E2F1 antibodies (ab245308, Abcam, Burlingame, CA, USA) or IgG-negative control antibodies (ab172730, Millipore, Billerica, MA, USA). RT-qPCR was used to detect the expression of circ_0001535 and DHFR.

Total RNA was extracted from Cerebrospinal fluid samples or cells with QIAamp RNeasy Micro Kit (74034, Qiagen, Dusseldorf, German). Reverse transcription was performed using 1 µg total RNA as the template and Advantage RT-for-PCR Kit (639505, Takara, shiga, Japan). Relative expression levels of genes were detected by RT2 SYBR${}^{\mathrm{\circledR }}$ Green qPCR Mastermixes (330509, Qiagen, Dusseldorf, German). Primer sets for qPCR are listed in Supplementary Table 1. The relative expression of cDNA was normalized to GAPDH and quantified using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ cycle threshold method.

CCK-8 assay was used to detect cell proliferation. AD cells were seeded in 96 well plates and assayed at 0-, 24-, 48- and 72-hours using CCK-8 Kit (ab228554, Abcam, Burlingame, CA, USA), following the manufacturer’s instructions. Cell viability was detected by using a microplate reader at 450 nm.

Flow cytometry was performed by using Alexa Fluor 488 Annexin V/Dead Cell Apoptosis Kit (V13241, Thermo Fisher Scientific, Waltham, MA, USA) following manufacturer’s instructions. For cell apoptosis analysis, cells were harvested and fixed in pre-cold 70% ethanol at 4 °C overnight. The cells were stained with Annexin V-FITC and propidium iodide (PI) and subsequently the ratio of apoptotic cells was tested by flow cytometry (BD FACS Calibur, Becton Dickinson, Franklin Lake, NJ, USA).

Differences between two or multiple groups were conducted using the Student’s t-test or one-way analysis of variance (ANOVA) followed by a post hoc Tukey’s test. Bioinformatics analysis was undertaken using R software (version 3.6.3, Statistics Department of the University of Auckland, Auckland, CA, USA) and GraphPad Prism (version 8.0, GraphPad Software, San Diego, CA, USA). p 0.05 were considered statically significant.

To select circRNA associated with AD pathology, we searched the literature for studies reporting RNA-seq analysis of AD patients and matched controls, and we focused on studies of human brain tissue of the lateral temporal lobe. GSE104704 was selected for brain tissue RNA-seq analysis. In total, 22 brain samples were examined, including 10 cognitively normal-aged brains, and 12 AD brains. 42,504 circRNA were selected (Supplementary Table 2). Among them, 4553 known circRNAs (Supplementary Table 3) were identified using CiRCexvoer2 software (version 2.0, CIRCexplorer, Shanghai, China).

Using edgeR to analyze the difference in expression of the above-screened circRNAs, we set at least 2-fold change thresholds in both expression aspects to draw volcano maps. In total, 135 circRNAs were found to be higher than the control and 102 lower than the control (Fig. 1A). CircRNAs with too low expression in AD patient group or control group were subsequently removed using the automatic filtering function of edgeR, which finally left 12 circRNAs with a significantly different expression between control and AD patient group. As shown in the Table 1, the chromosomal location of the parental gene where the circRNA was located. CircbaseID and p value were recorded. Compared to the expression levels in control group, the expression levels of circ_0000497, circ_0001519, circ_0001535, circ_0002454, circ_0001380 and circ_0087960 were up-regulated in AD group (p 0.05). The expression levels of circ_0006916, circ_0131934, circ_0099634, circ_0005114, circ_0017248 and circ_0006837 were down-regulated (Fig. 1B, p 0.05).

Fig. 1.

Differential expression analysis of circRNA between Alzheimer’s disease (AD) and Health doners. (A) Volcano plot shows differential expressed circRNAs. (B) Differential expression analysis of 12 circRNAs between AD and Health groups. *p 0.05 compared with control (CTL) group. CPM, count-per-millon.

To assess the predictive ability of 12 circRNAs for AD, we used the ROC curve of gene expression abundance to calculate the AUC and significance (Fig. 2). The AUC of all 12 circRNAs was greater than 0.7. Except for circ_0005114, the p value of the other 11 circRNAs was less than 0.05. Among these circRNAs, circ_0001535 had the largest AUC value of 0.900 and a p value of 0.002 (Fig. 2). Overall, all 12 circRNAs performed well in determining AD, with circ_0001535 performing best in determining whether the sample had AD. In addition, the detection of circRNA expression levels in clinically collected samples from AD patients and healthy samples revealed that the level of circ_0001535 in the AD group was about three times higher than that in the healthy group, most significant in 12 circRNAs (Fig. 3, p 0.01).

Fig. 2.

ROC curve of 12 representative circRNAs. ROC, receiver operating characteristic curve; AUC, area under curve.

Fig. 3.

The expression levels of 12 circRNAs in cerebrospinal fluid of AD and control group were detected by RT‑qPCR. ns, no significant difference compared with Control group; *p 0.05 compared with Control group. RT‑qPCR, reverse transcription‑quantitative PCR.

In order to investigate the possible pathways affecting the differential expression of circ_0001535 in AD patients, AD patients were first divided into high and low-expression groups according to the median expression of circ_0001535 in AD patients. Then the GSEA method was used to screen out the pathways with significant enrichment differences between the high and low-expression groups. As shown in the Fig. 4, circ_0001535 was associated with six major pathways.

Fig. 4.

GSEA applied to validate the hub gene signatures which were related to circ_0001535. GSEA, gene set enrichment analysis.

In addition, we further screened the genes with a high correlation with the expression of circ_0001535. A total of 503 genes were screened which were significantly associated with circ_0001535 expression using Pearson correlation analysis at p 0.05 (Supplementary Table 4). This gene list was analyzed using Metascape Gene List Analysis for pathway and process enrichment using the following databases: KEGG Pathway, GO Biological Process, Reactome Gene Set, Canonical Pathways, CORUM, TRRUST, DisGeNET, PaGenBase, Transcription. As shown in Fig. 5A, the top 20 clusters are shown here in descending order of −Log10(p) values for their representative enriched terms (Fig. 5A). To further analyze the relationships between the clustered terms, a subset of each cluster was selected and visualized by Cytoscape to network its presentation as a network graph. Each node represents a rich term and is first colored by its cluster ID in these networks (Fig. 5B).

Fig. 5.

Metascape gene list analysis results shows pathways related to circ_0001535. (A) Network of enriched terms colored by cluster ID, where nodes that share the same cluster ID are typically close to each other. (B) Protein-protein interaction network and MCODE components identified in the gene lists. (C) Bar graph of enriched terms across input gene lists, colored by p-values. MCODE, molecular complex detection.

Futhermore, protein-protein interaction enrichment analysis was performed using the following databases: STRING, BioGrid, OmniPath and InWeb_IM9. Only physical interactions in STRING (physical score $>$0.132) and BioGrid were used. The network contains a subset of proteins that form physical interactions with at least one other member of the list (Fig. 5C). The transcription factors in the network could also be considered potential and studied further.

In order to prove the potential transcription factors regulated by circ_0001535, an analysis of transcription factors that may affect the expression of these genes using TRRUST. The results revealed that E2F1 is the key transcription factor. It may affect 7 genes: DHFR, ERCC1, FSHR, HIC1, RRM2, SIP2 and TOP2A (Table 2, Fig. 6A). DHFR was the primary downstream gene that might be affected (p = 0.038). Subsequently, the binding relationship between E2F1 and DHFR was examined using JASPAR, and 2 binding sites were found to exist on chromosome 5 (Table 3, Fig. 6B).

Fig. 6.

TRRUST and JASPAR were used to find downstream genes effected by circ_0001535. (A) Summary of enrichment analysis in TRRUST. (B) JASPAR analyze results shows two binding sites of E2F1 on DHFR.

SH-SY5Y cells were used as model cells to test the interaction relationship between circ_0001535 and its downstream regulators. A$\beta$ protein were added into SH-SY5Y cells to construct the AD cell model. Firstly, circ_0001535, E2F1 and DHFR expressions were evaluated using the RT-qPCR method. Results showed that the expressions of circ_0001535, E2F1 and DHFR were significantly higher in the AD cell model (Fig. 7A, p 0.05). Next, RIP and ChIP assays were processed to identify the binding relationship among circRNA, transcription factor and downstream effector gene in SH-SY5Y cells. RIP result showed circ_0001535 has a binding ship to E2F1 (Fig. 7B, p 0.05). ChIP result showed E2F1 could bind to DHFR promoter (Fig. 7C, p 0.05).

Fig. 7.

Identification of regulation ship between circ_0001535 and its downstream regulators. (A) The differential expression of circ_0001535, E2F1 and DHFR between SH-SY5Y cells and AD model cells group was estimated by RT-qPCR. (B) RIP assay showed the regulation ship between circ_0001535 and E2F1. (C) ChIP assay showed the bind ship between E2F1 and DHFR. *p 0.05 compared with SH-SY5Y cell group or IgG group. DHFR, dihydrofolate reductase.

Overexpression plasmids of circ_0001535 and si-RNA of E2F1 and DHFR were created to identify biological functions of the circ_0001535/E2F1/DHFR axis. First, the results of RT-qPCR showed that circ_0001535 overexpression could promote DHFR expression, while E2F1 knock down could inhibit DHFR expression (Fig. 8A, p 0.05). Then, CCK-8 results showed circ_0001535, E2F1 and DHFR could down-regulate cell viability of AD cell model (Fig. 8B, p 0.05). In apoptosis assay, similar results could be observed, circ_0001535, E2F1 and DHFR all could promote cell apoptosis of the AD cell model (Fig. 8C, p 0.05).

Fig. 8.

Identification of biological functions of circ_0001535/E2F1/DHFR axis. (A) The differential expression of circ_0001535, E2F1 and DHFR among different groups was estimated by RT-qPCR. (B) The proliferation ability of AD model cells in different groups was estimated by CCK-8 assay. (C) Flow cytometry were used to identify regulation ship of circ_0001535/E2F1/DHFR axis on AD cell apoptosis. *p 0.05. CCK8, Cell Counting Kit-8; NC, Negative control.