We recruited 63 patients With PD and 60 healthy individuals From the First Affiliated Hospital of Guangxi Medical University Between January 2019 and June 2021. Patients With PD Were diagnosed independently by Two neurologists based on the United Kingdom Parkinson’s Disease Brain Bank Criteria (Hughes et al., 1992). the control groups consisted of healthy people matched With PD in terms of gender and age and Had no known neurological disease, comorbidities, or PD family history. the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University approved this study. All participants provided informed consent. the clinical and demographic characteristics of all participants Are outlined in Table 1. A total of 123 peripheral blood samples Were obtained From patients With PD and healthy individuals and stored vertically at −80°C until total RNA isolation.

RNA Isolation and Real-Time Quantitative Polymerase Chain ReactionTRIzol reagent kit (Invitrogen, Takara, Japan) Was used to extract total RNA From blood samples according to the protocol provided by the manufacturer. After extraction, total RNA Was reverse transcribed Into cDNA (RR047A, Takara, Japan). the RT-qPCR Was performed on the ABI PRISM 7500 (RR820A, Takara, Japan). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) Was used as an internal reference to normalize the qPCR data. the 2^–ΔΔCT method Was used to calculate the relative expression of circRNAs. the primers Are described in Supplementary Table 1.

Sample preparation and microarray hybridization were performed according to the protocol of the Arraystar. KangChen Bio-tech Inc. (Shanghai, China) carried out the Arraystar Human circRNA Array version 2 analysis. To summarize, circRNAs from each sample were enriched by removing linear RNAs and then amplified following the manufacturer’s instructions. The labeled circRNAs were hybridized onto the Arraystar and scanned by the Agilent Scanner. CircRNAs exhibiting a fold change (FC) ≥1.5 with a p-value ≤ 0.05 were considered significantly different. Volcano plot filtering and fold change filtering were used to identify differentially expressed circRNAs. The distinct circRNA expression patterns among samples were analyzed using hierarchical clustering. The circRNA microarray raw data have been stored in the Gene Expression Omnibus (GEO) (series number: GSE198273).

The area under the curve (AUC) of the receiver operating characteristics (ROC) is used to assess the accuracy of biomarkers. The AUC is near 1, which means an ideal model (Pepe et al., 2009). The following guidelines are often used to estimate the specificity and sensitivity of a biomarker, e.g., excellent scored 0.9–1.0; good scored 0.8–0.9; acceptable scored 0.7–0.8; while poor scored 0.6–0.7 and fail scored 0.5–0.6 (Xia et al., 2013; Ren et al., 2018). The circRNA level with an AUC > 0.70 was selected to proceed forward with bioinformatics analysis to investigate potential functions.

We constructed circRNA–miRNA–mRNA networks in Cytoscape based on circRNA microarray data and the RT-qPCR results. The circRNA–miRNA interaction was predicted using Arraystar’s miRNA target prediction software. The Cytoscape plug-in ClueGO (2.5.8) and CluePedia (1.5.8) were used to predict target genes of the top 5 putative miRNAs for the candidate circRNAs, and the top 200 target genes were screened for the construction of the circRNA–miRNA–mRNA networks.

Gene Ontology (GO) enrichment analysis was used to establish functional annotations of miRNA target genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was used to explore the functions of these differentially expressed genes. Statistical analyses and visualizations were conducted in R 3.6.3 by using the clusterProfiler package and the org.Hs.eg.db package.

The protein–protein interaction (PPI) network was created by using the Retrieval of Interacting Genes (STRING) online tool to evaluate the interactions among the target genes that were identified in the circRNA-miRNA-mRNA networks. The hub genes were screened from the PPI network using the cytoHubba plugin.

To examine the upstream regulation mechanism of the candidate circRNAs, transcription factors for circRNAs were screened using the TRCirc database¹.

Statistical analyses were conducted using SPSS 20.0 software (IBM, Armonk, NY, United States). The continuous variables were presented as mean ± standard deviation (SD). The differences in age and sex between patients with PD and healthy controls were assessed using the chi-square (χ²) test or Student’s t-test. ROC curves were used to assess the potential biomarkers of the candidate circRNAs. Two-tailed p-values < 0.05 were considered significantly difference.

Human circRNA microarray was used to screen dysregulated circRNAs from three patients with PD and three healthy controls. The mean age was 62.67 ± 8.36 years for the PD group and 63.00 ± 8.67 years for controls. Both groups have two male participants and one female participant. The expression features of the dysregulated circRNAs were depicted in a two-dimensional hierarchical clustering heatmap (Figure 1A). A volcano plot (Figure 1B) and a scatter diagram (Figure 1C) were also utilized to show the differences in circRNA expression between both groups. According to microarray data analyses, a total of 10,145 circRNAs were identified, with 189 antisense (1.86%), 8,518 exonic (83.96%), 73 intergenic (0.72%), 700 intronic (6.90%), and 665 sense overlapping (6.55%) (Figure 1D). After scanning and normalization, 139 differentially expressed circRNAs were identified with a threshold of fold-change (FC) absolute value ≥1.5 and P < 0.05. Of them, 78 circRNAs were upregulated and 61 were significantly downregulated. The top 5 upregulated and downregulated circRNAs are presented in Table 2.

A total of 10 candidate circRNAs (five upregulated and five downregulated) were retrieved for further verification in an independent cohort of 63 patients with PD and 60 healthy controls. The clinical characteristics of this cohort are described in Table 1. The mean age was 60.16 ± 10.51 years for the PD group and 59.25 ± 10.20 years for controls. Statistical analyses showed no significant difference in sex or age between both groups (χ² = 0.004, p = 0.951; t = 0.486, p = 0.628, respectively). The average UPDRS III score was 31.30 ± 12.43 for the PD group, and the disease duration was 4.44 ± 3.61 years. Enrolled patients had an average age of onset of 60.16 ± 10.51 years, and the mean clinical Hoehn and Yahr stage was 2.37 ± 0.77.

The ten candidate circRNAs were verified using divergent primers rather than the standard convergent primers via RT-qPCR. Of the 10 circRNAs, 4 were confirmed as differentially expressed between both groups (Figures 2A–J). The circRNAs hsa_circRNA_103730, hsa_circRNA_101275, and hsa_circRNA_038416 (Figures 2A–C) were significantly upregulated, while hsa_circRNA_102850 (Figure 2D) was downregulated in patients with PD when compared to healthy controls. Considering that the incidence of male patients with PD is greater than female patients, we performed a stratification analysis based on sex. However, the results did not show a significant difference between male and female patients with PD (Supplementary Figure 1A–J).

The diagnostic value of the differentially expressed circRNAs in PD was performed using the ROC curve. The AUC for hsa_circRNA_103730, hsa_circRNA_102850, hsa_circRNA_101275, and hsa_circRNA_038416 for PD was 0.701, 0.801, 0.843, and 0.674, respectively (2K–N). A circRNA panel combining the four candidate circRNAs showed a higher diagnostic ability to distinguish patients with PD from healthy controls (AUC = 0.938) (Figure 2O).

Generally, ncRNAs act as competing endogenous RNAs (ceRNAs) that regulate target RNA transcripts by binding competitively to shared miRNAs. It has been well established that circRNAs have many miRNA-binding sites that may act as a ceRNA sponge for miRNAs. Considering the circRNA with an AUC > 0.70 was considered with sufficient specificity and sensitivity, therefore, hsa_circRNA_102850, hsa_circRNA_101275, and hsa_circRNA_103730 were selected for the circRNA-regulatory network prediction. Interactions between the three candidate circRNAs and their relationships with miRNAs and target genes were predicted by building a unified interaction-network model. The interactions between circRNAs and miRNAs were predicted by the miRNA target prediction software. The top five putative miRNAs of the three candidate circRNAs were selected to predict target genes, and the top 200 target genes were used to construct the circRNA–miRNA–mRNA interaction networks (Figure 3). The results indicated that the candidate circRNAs may act as endogenous RNAs for modulating target gene expression.

The TRCirc was used to establish a transcription factor–circRNA (TF–circRNA) network to explore the precise molecular mechanism of transcriptional regulation of the candidate circRNAs. Hsa_circRNA_103730 (circBase ID: hsa_circ_0005654) and hsa_circRNA_101275 (circBase ID: hsa_circ_0030428) were predicted to be regulated by targeting TH-relevant transcription factors such as GATA2 and GATA3 (Polanski et al., 2010).

The GO and KEGG enrichment analyses were conducted to investigate the biological function of the candidate circRNAs. The GO analysis indicated that the 200 mRNAs were mainly enriched in negative regulation of phosphorylation, cellular response to environmental stimulus, response to a steroid hormone, transcription factor complex, collagen-containing extracellular matrix, and nuclear transcription factor complex (Figure 4A). The PI3K–Akt signaling pathway and MAPK signaling pathway were found to have significant enrichment in the KEGG pathway analysis (Figure 4B).

We constructed a PPI network according to circRNA–miRNA–mRNA interaction analysis. The PPI networks identified 107 proteins and 200 edges of the candidate circRNA networks that may be involved in PD (Figure 5A). We used the cytoHubba plugin to find the critical genes in the process of PD pathogenesis and obtained a sub-network containing 10 genes, including ESR1, PTEN, SHC1, IGF1R, SMAD2, KRAS, MDM2, HIF1A, BMP4, and ACVR2B (Figure 5B). The results of PPI interaction networks showed that ESR1, PTEN, and SHC1 may be involved in PD as the most important targets.