sPD patients and control subjects (n = 530 each) from the HPMC were recruited at the First Affiliated Hospital of Nanchang University. A modified version of the United Kingdom Parkinson’s Disease Society Brain Bank Clinical Diagnostic Criteria (Gibb and Lees, 1988; Albanese, 2003) was used for sPD diagnosis. The patients exhibited at least two of following signs: resting tremor, rigidity, bradykinesia, and gait or postural disorder. All participants underwent the same clinical and laboratory examinations. We also reviewed general medical history, PD family history, related diseases in first-degree relatives, and potential exposure to toxicants related to sPD; and we performed magnetic resonance imaging of the brain and spinal cord to exclude other neurologic diseases with clinical manifestations similar to those of sPD such as tumors, demyelinating disorders, hydrocephalus, and cervical myelopathy. All procedures were approved by the institutional review board of the human ethics committee of the First Affiliated Hospital of Nanchang University. Written informed consent was obtained from all participants or their legal guardians before patient enrollment. All experimental methods were in compliance with the principles of the Declaration of Helsinki.

We downloaded the genome-wide p-value data for sPD-related GWASs from the dbGap database of National Center for Biotechnology Information (NCBI).¹ The GWAS data were from four studies of American and Caucasian (Germany and the United Kingdom) populations (Maraganore et al., 2005; Fung et al., 2006; Simón-Sánchez et al., 2009; Pankratz et al., 2011) involving 7,025 sPD cases (including 816 familial PD cases) and 10,631 controls (Figure 1). We performed pairwise comparisons of the gene sets of the four studies based on the p-values. A total of 2,779 SNPs in 52 genes were identified (Supplementary Table 1). We compared the 52 genes with our previous sALS data from a GWAS of the HPMC (Xie et al., 2014) and found 23 that were common to both sPD and sALS. We selected 27 SNPs in 10 candidate genes for allelotyping and genotyping (Supplementary Tables 2–4).

DNA was isolated from white blood cells of sPD patients and control subjects using a DNA extraction kit (Fastgene Co., Shanghai, China) according to the manufacturer’s instructions. The samples were verified by electrophoresis on a 1.5% agarose gel to ensure that there was no RNA contamination or DNA degradation. Samples of intact genomic DNA were selected for allelotyping and genotyping analyses. For the former, the false-positive report probability along with the relative allele frequency was calculated for loci surrounding candidate genes to estimate the confidence interval (CI) and p-value corresponding to the odds ratio (OR) score (Wacholder et al., 2004). Candidate polymorphic loci were genotyped using the Sequenom Mass ARRAY iPLEX platform (Agena Bioscience, San Diego, CA, United States). All experimental procedures including primer design and synthesis and genotyping were carried out by I-conoene Biotechnology Co. (Wuhan, China). The methods of genotyping were performed by Sequenom MassArray. For quality control, 5% of the DNA samples were randomly selected; and the assay was repeated, with a concordance rate of 100% (Su et al., 2015).

Demographic data (continuous variables) of sPD patients and control subjects are expressed as mean ± standard deviation (SD). Differences between sPD patients and control subjects were evaluated with the Student’s t-test, χ² test, or Fisher’s exact test using SPSS v17.0 (SPSS Inc., Chicago, IL, United States). The allele frequency of loci was estimated using the ratio of signal intensities of each allele. Differences in allele frequencies between sPD patients and controls were evaluated with the combined Z test. Statistical analyses were performed with SAS v9.1.3 (SAS Institute, Cary, NC, United States). A p-value < 0.05 was considered statistically significant.

Human long non-coding RNA (lncRNA) sequences in the GRCh38 assembly were downloaded from the Ensembl Genome Browser website² and formatted as a database. Intron sequences from 50 bp upstream to 50 bp downstream of DAB1 rs17115303 and PTPRT rs6030462 were extracted from GRCh38 and searched in the database by BLAST with an e-value cutoff of 1,000 to obtain high-scoring segment pair (HSP) hits with the complementary intron sequence.

A nervous system regulatory network was constructed using transcription factors that bind to the promoter region of DAB1 [V-rel avian reticuloendotheliosis viral oncogene homolog A (RELA), myoblast determination protein 1 (MYOD1), and MDS1 and EVI1 complex locus (MECOM)] and PTPRT [Kruppel-like factor 4 (KLF4) and Ras-responsive element binding protein 1 (RREB1)] genes along with risk genes for ALS [fusion, derived from T(12;16) malignant liposarcoma (FUS), uncoordinated 13 homolog A (UNC13A), activator protein complex subunit beta (APB), peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PPARGC1A), hemoprotein 2B (HMP2B), peripherin (PRPH), matrin 3 (MATR3), triggering receptor expressed on myeloid cells 2 (TREM2), coiled-coil-helix-coiled-coil-helix domain-containing 10 (CHCHD10), dynactin subunit 1 (DCTN1), optineurin (OPTN), TANK-binding kinase 1 (TBK1), ubiquilin 2 (UBQLN2), D-amino acid oxidase (DAO), neurofilament heavy (NEFH), FIG4, retinoblastoma-binding 4 (RBB4), SOD1, HNRNPA1, valosin-containing protein (VCP), chromosome 9 open reading frame 72 (C9ORF72), TAR DNA-binding protein (TARDBP), sequestosome 1 (SQSTM1), angiogenin (ANG), paraoxonase 1 (PON1), PON2, PON3, profiling 1 (PFN1), and ataxin 2 (ATXN2)] obtained from the OMIM database (OMIM 105400³).

All 10 identified risk genes common to sPD and sALS in our study are listed in the Gene Ontology (GO) database.⁴ Functional prediction analyses were performed by Gene For Health Biotechnology Co. (Shanghai, China). We selected “hsapiens” as the organism and “hsapiens_gene_symbol” as the gene ID type when uploading the genes of interest and reference gene set. We selected “hsapiens_genome” as the reference set, p < 0.01 as the significance level, and hypergeometric as the statistical method.

Studied populations were composed of a total of 530 sPD and 530 controls. The average age at the onset (sPD) was 58.72 ± 0.91. The average age at the enrollment (Controls) was 59.56 ± 0.92. The percent of males was 51.70%; that of control was 60%. The median (range) age of sPD was 60 (32–83) years, and that of control was 60 (31–83) years (Supplementary Table 5).

The published genetic data of GWAS dataset from four origins of non-HPMC composed of the United States, Germany, and the United Kingdom were from performed enrichment analysis, which showed that pha000004 was drastically enriched at the top of list from pha002840 [enrichment score (ES) = 0.409, normalized ES (NES) = 1.455, nominal (NOM) p-value < 10^–3]. The significant enrichment of pha002840 in the pha000004 data further indicated the greater consistency across two genetic datasets (ES = 0.398, NES = 1.602, NOM p-value < 10^–3) (Figure 2A); pha000004 was drastically enriched at the top of list from pha002865 (ES = 0.393, NES = 1.555, NOM p-value < 10^–3). The significant enrichment of pha002865 in the pha000004 data further indicated the greater consistency across two genetic datasets (ES = 0.397, NES = 1.610, NOM p-value < 10^–3) (Figure 2B); pha000004 was drastically enriched at the top of list from pha003128 (ES = 0.335, NES = 1.739, NOM p-value < 10^–3). The significant enrichment of pha003128 in pha000004 data further indicated the greater consistency across two genetic datasets (ES = 0.372, NES = 1.520, NOM p-value < 10^–3) (Figure 2C); pha002840 was drastically enriched at the top of list from pha002865 (ES = 0.390, NES = 1.521 NOM p-value < 10^–3). The significant enrichment of pha002865 in the pha002840 data further indicated the greater consistency across two genetic datasets (ES = 0.398, NES = 1.410, NOM p-value < 10^–3) (Figure 2D); pha002840 is drastically enriched at the top of list from pha003128 (ES = 0.357, NES = 1.825, NOM p-value < 10^–3). The significant enrichment of pha003128 in pha2840 data further indicated the greater consistency across two genetic datasets (ES = 0.397, NES = 1.415, NOM p-value < 10^–3) (Figure 2E). pha002865 was drastically enriched at the top of list from pha003128 (ES = 0.361, NES = 1.862, NOM p-value < 10^–3). The significant enrichment of pha003128 in the pha002865 data further indicated the greater consistency across two genetic datasets (ES = 0.405, NES = 1.608, NOM p-value < 10^–3) (Figure 2F). We obtained 2,779 SNPs of 52 genes from the enrichment results (Supplementary Table 1; Figure 3). Then we conducted an integrative analysis with our previously reported candidate genes of sALS from our GWAS data (Xie et al., 2014). Of these, our results showed that 23 genes were overlapped genes in both sPD and sALS; they might be the potential common pathogenic genes of both sPD and sALS.

From 23 genes of overlapped genes in both sPD and sALS, we chose 27 SNPs in 10 candidate genes for further allelotyping and genotyping analyses (Supplementary Tables 2–4). The result determined that two novel SNPs, DAB1 rs17115303 and PTPRT rs6030462, shared in both sPD and sALS. They both were SNPs of introns. Gene DAB1 base sequence is located in the chr1:57194802 position; DAB1 rs17115303 changed from C to A at the reverse strand (Figure 4A). Gene PTPRT base sequence is located in the chr20:42771589 position; PTPRT rs6030462 changed from G to A at the reverse strand (Figure 4B).

Allelic and genotypic frequencies are summarized in Tables 1, 2. Two novel SNPs were identified, as follows: rs17115303 in DAB1 gene and rs6030462 in PTPRT gene. The genotype frequency of rs17115303 differed significantly (p = 0.047) between patients and controls. The minor allele frequency (MAF) of rs17115303 was higher in the cases (8.8%) than in the controls (6.8%), and the difference was significant (p = 0.039). The minor A allele of rs17115303 (OR 1.42, 95% CI 1.024–1.970) exhibited an increased risk for the disease. C allele of rs17115303 (OR 0.704, 95% CI 0.508–0.976) exhibited a decreased risk for the disease. p-Value of rs6030462 genotype frequency was significant (p = 0.019) between patients and controls. The MAF of rs6030462 was higher in cases (44.2%) than in controls (41.3%), and the difference was not significant (p = 0.306). The minor A allele of rs6030462 (OR 1.125, 95% CI 0.899–1.407) exhibited an increased risk for the disease. G allele of rs6030462 (OR 0.889, 95% CI 0.711–1.112) exhibited a decreased risk for the disease. Minor alleles of rs17115303 in DAB1 gene and rs6030462 in PTPRT gene were correlated factors with the disease, and they significantly increased the risk of disease.

The sequence from upperstream 100 bp to downstream 100 bp of DAB1 rs17115303 mutation position was analyzed using NHRscan software. The position from 98 to 117 of sequence in the fragment of red frame in Figure 5A, ER8: GGAGATGCGCTTTGAGGACT, was predicted to be the binding site fragment, and the mutation site was at the fourth base position of binding fragment (Figures 5A,B). The sequence in the upstream 1,500 bp of DAB1 rs17115303 was applied to predicted transcription factor binding sites using the Mscan software. The binding site was predicted at distances from 1,300 to 1,420 bp and 900 to 1,100 bp upstream of the binding site; the predicted results based on the Mscan software showed that the DAB1 rs17115303 mutation position might be in the binding block (Figures 5E,F).

The secondary structure of sequence from the upperstream 100 bp to the downstream 100 bp of DAB1 rs17115303 mutation position was predicted using the Mfold software and visualized potential binding sites. The intron sequence from the upperstream 100 bp to the downstream 100 bp of DAB1 rs17115303 mutation position was extracted from the GRCh38 of human genome version. The result showed that the binding region formed a hairpin-like structure in the secondary structure. By overlaying observed SNPs and predicted binding sites on the secondary structure of intron, it was possible that the intron used a hairpin structure for the protein binding. The observed intronic SNP was located in the cap region, possibly harmed the stability of hairpin, and in turn influenced the splicing event (Figure 6A).

The intron sequence from the upperstream 50 bp to the downstream 50 bp of DAB1 rs17115303 mutation position was extracted from the version GRCh38 of human genome. The each sequence of lncRNA was applied to blast against the database. Blastn with e-value cutoff 1,000 was performed to obtain short HSP hits with the complemented intron sequence. The result obtained three lncRNA hits. Those hits with the opposite orientation compared with intron sequences (Plus/Minus) were filtered as the lncRNA that could bind the transcriptional intron, but not directed on the DNA strand (Figure 7A).

The nervous system regulation network building from transcription factors that bind to the promotor region of DAB1 gene showed that DAB1 gene mainly participated in the regulation of protein phosphorylation, the glial cell differentiation, and the regulation of neuron death, which are majorly involved in SOD1, PPARGC1A, TBK1, VAPB, ANG, and SQSTM1 genes; their related genes including MATR3, ATXN2, UBQLN2, TARDBP, FUS, OPTN, FIG4, HNRNPA1, DCTN1, CHMP2B, and VCP genes; and KLF4, RELA, and MYOD1 transcription factors (Figure 8A).

The sequence from the upperstream 100 bp to the downstream 100 bp of PTPRT rs6030462 mutation position was analyzed using the NHRscan software. The position from 83 to 102 of sequence in ER8: GGAAGTGGTGCAGGCTGTCG was predicted to be the binding site fragment, and the mutation site was at the 19th base position of binding fragment (Figures 5C,D). The sequence in the upstream 1,500 bp of PTPRT rs6030462 was applied to transcription factor binding sites prediction using the Mscan software. The binding site was predicted at distances from 725 to 531 bp upstream of binding site; predicted results based on the Mscan software showed that the PTPRT rs6030462 mutation position might be in the binding block (Figure 5G). The secondary structure of sequence from upperstream 100 bp to downstream 100 bp of PTPRT rs6030462 mutation position was predicted using the Mfold software and visualized potential binding sites. The intron sequence from the upperstream 100 bp to the downstream 100 bp of PTPRT rs6030462 mutation position was extracted from the version GRCh38 of human genome. The result showed that the binding region formed a hairpin-like structure in the secondary structure. By overlaying observed SNPs and predicted binding sites on the secondary structure of the intron, it was possible that the intron used a hairpin structure for the protein binding. Observed intronic SNPs were located in the cap region, which possibly harmed the stability of hairpin and in turn influenced the splicing event (Figure 6B).

Human lncRNA sequences were downloaded from the version GRCh38 of Ensemble website and formatted as the database. The intron sequence from the upperstream 50 bp to the downstream 50 bp of PTPRT rs6030462 mutation position was extracted from the version GRCh38 of human genome. Each sequence of lncRNA was applied to blast against the database. Blastn with e-value cutoff 1,000 was performed to obtain short HSP hits with the complemented intron sequence. The result found 23 lncRNA hits, and those hits with the opposite orientation compared with the intron sequence (Plus/Minus) were filtered as the lncRNA, which could bind the transcriptional intron, but not directed on the DNA strand (Figure 7B).

The nervous system regulation network building from transcription factors that bind to the promotor region of PTPRT gene showed that PTPRT gene mainly participates in the positive regulation of programmed cell death and dephosphorylation, which are majorly involved in SOD1, SQSTM1, VCP, and FIG4 genes; their related genes including MATR3, HNRNPA1, TARDBP, ATXN2, UBQLN2, DCTN1, FUS, OPTN, VAPB, NEFH, PRPH, CHMP2B, ANG, TBK1, and PPARGC1A genes; and KLF4, RELA, and MYOD1 transcription factors (Figure 8B).

Ten candidate genes were involved in 40 of the biological processes. Of them, DAB1-related biological processes included 36 of the biological processes, which mainly consisted of the negative regulation of cell adhesion, the phosphate-containing compound metabolic process (dephosphorylation), the cell development (neuron, cerebral cortex, forebrain, neuron projection, pallium, dendrite, and brain), the regulation of cell morphogenesis involved in differentiation, the regulation of neuron differentiation, the regulation of response to external stimulus, the regulation of neurogenesis, the locomotory and single-organism behavior, the differentiation, migration, generation recognition of neurons, neurogenesis, and the regulation of axonogenesis. PTPRT-related genes included six of the biological processes, which were the biological adhesion, the cell adhesion, the cell–cell adhesion, the phosphate-containing compound metabolic process, the phosphorus metabolic process, and dephosphorylation, which mainly participated in the adhesion and phosphorus metabolic processes (Supplementary Table 6 and Supplementary Figure 1). The commonly shared biological processes were the biological adhesion, the cell adhesion, the cell–cell adhesion, the phosphate-containing compound metabolic process, and the phosphorus metabolic process. It indicated that DAB1 and PTPRT genes existed in very similar biological processes.