Identifying the potential genes in alpha synuclein driving ferroptosis of Parkinson’s disease

Parkinson’s disease (PD) is a common neurodegenerative disease with aggregation of α-synuclein (α-syn) in substantia nigra (SN). The association between the α-syn and ferroptosis in PD remains unclear. GSE49036 was obtained from the Gene Expression Omnibus (GEO) database and intersected with ferroptosis genes. Bioinformatics analysis was used to identify the potential differentially expressed genes (DEGs) included the development of Gene set enrichment analysis (GSEA), Kyoto Encyclopedia of Genes and Genomes (KEGG) and protein–protein interaction (PPI) network. We screened 8 key genes were modulated and crosslinked by 238 miRNAs. Additionally, 5 hub genes were predicted and 38 lncRNAs targeting 3 key miRNAs were revealed. Finally, 3 hub genes (PIK3CA, BRD4, ATM) and the key lncRNA (NEAT1) were verified in neurotoxic PD models. The in vitro experiments showed that PIK3CA and ATM were significantly upregulated or the BRD4 was downregulated in the rotenone treatment and they could be rescued by the specific ferroptosis inhibitor, liproxstatin-1. The expression of the key lncRNA NEAT1 were consistent with the hub genes in same models. This study identified the proposed NEAT1-PIK3CA/ATM ceRNA network may be a specific biomarker in α-syn driving ferroptosis as well as to predict clinical outcomes and therapeutic targets in PD patients.


PPI network construction
The DEGs were processed using the STRING database, an online tool that retrieved the interaction between proteins.We used to predict the correlation among the DEGs (confidence score > 0.4) and the interaction relationship was visualized by Cytoscape v3.9.1 software.The nodes represented proteins and the edges represented the links between them.Moreover, the molecular complex detection (MCODE) and the CytoHubba, two plugins of Cytoscape were carried out to screen hub genes.The MCODE was used for clustering analysis of key gene sub-networks and the CytoHubba was employed to select the top 10 hub genes via 11 maximal clique centrality (MCC) methods for the further study.

Gene-miRNA network construction
The hub 9 genes, selected according to the PPI network, were analysed to predict targeted pivotal miRNAs and build the gene-miRNA interaction relationship networks.We applied miRWalk 2.0 to integrate miRNA databases and intersected the predicted results of the miRDB database to ensure the accuracy of the results.We visualized the networks by employing Cytoscape v3.9.1 and analysed their biological pathways by Funrich, an functional enrichment analysis tool special for miRNAs.

miRNA-lncRNA prediction and ceRNA network construction
Based on StarBase v2.0, the upstream target lncRNAs of the predicted miRNAs were identified and then intersected with the differential lncRNAs.We selected the cross-linked hub genes showing the interaction of the miRNAs relevant lncRNAs.The highest reliability (very high stringency, > 5) was selected as the criterion.The ceRNA network of the hub genes-miRNAs-lncRNAs was constructed by the Cytoscape v3.9.1 software.

Cell cultures and drug treatment
Human neuroblastoma SH-SY5Y cells (#CL-0208) were purchased from Procell (Wuhan, China) and cultured in Minimum Eagle Medium (MEM, Procell, Wuhan, China) supplemented with 10% fetal bovine serum, 1% penicillin and streptomycin.Cells were incubated at 37 °C and 95% humidity with 5% CO 2 .To contribute PD cell culture model in vitro, SH-SY5Y cells were seeded in 1 ml media at a density of 1 × 10 5 per well in a 12-well plate for approximately 24 h., then were treated with indicated drugs respectively.In control group, cells were exposed to medium (containing 0.1% dimethyl sulfoxide).In model group, cells were treated with 1 µmol/L (µM) rotenone (#S2348, Selleck, Houston, USA) with or without 5 µmol/L (µM) liproxstatin-1 (#S7699, Selleck, Houston, USA).After being treated with 24 h, the cells were digested and washed with phosphate buffer saline and subjected to further treatment and analysis.

RNA extraction and quantitative real time-PCR (qRT-PCR) analysis
The total SH-SY5Y cell RNA was extracted using the TRIzol reagent, the quality was tested using Quantus Fluorometer/spectrometer and the assessment of purity relied on the A260/A280 ratio.The cDNA synthesis was performed using reverse transcription reaction kit (Vazyme) following the manufacturer's instructions.Quantitative real-time PCR was carried out using the Hieff Universal qPCR SYBR Green Master Mix (UNICON, Yeasen Biotechnology, Shanghai, China).Primer sequences used for qRT-PCR assays were listed in Table 1.

Lipid peroxidation, GSH and iron content assay
SH-SY5Y cells were seeded in 12 well plates at a density of 1 × 10 5 cells/well and treated with indicated drugs: normal group and rotenone group with or without liproxstatin-1.The cells from each sample were washed by PBS buffer, then homogenized and harvested by centrifugation at 1000 rpm at 4 °C for 5 min.The quantifications of lipid peroxidation (MDA) (#A003-1, Jiancheng Bioengineering Insitute, Nanjing, China), GSH (#A006-1-1, Jiancheng Bioengineering Insitute, Nanjing, China) and iron (#ab83366, Abcam, Cambridge, USA)in the collected supernatant were performed respectively with commercial assay kits according to the manufacturer's instructions.

Statistical analysis
All data are presented as mean ± standard deviation (SD).Two-tailed Student's t test or one-way ANOVA followed by Tukey's post hoc analysis was used to compare the differences between groups.Each experiment was replicated at least three times.Statistical analysis and visualization were performed using GraphPad Prism 8.The significance of differences was set at the level of p < 0.05.

Ethics statement
All the information of patients were recruited from a public database and the study involving none of ethics issue.

Identification of DEGs
In total, the mRNA expression data of 16 substantia nigra tissue generated from GSE49036 with 8 PD patients (Braak α-synuclein Stage 5-6) and 8 controls (Braak α-synuclein Stage 0) were obtained.Compared with normal samples, a total of 3342 DEGs in the dataset were found (Fig. 1A,B).We also gained the dataset including 259 genes from FerrDb and intersected them with GSE49036 to identify 50 ferroptosis DEGs (Fig. 1C).Among the DEGs, 27 up-regulated and 23 down-regulated DEGs in the substantia nigra tissue were identified and they were further classified as ferroptosis driver, ferroptosis suppressor, ferroptosis marker or unclassified via the FerrDb online tool (Tables 2, 3).

GSEA analysis of DEGs
We used GSEA enrichment analysis in the KEGG pathway by WebGestalt software to uncover the underlying mechanisms of the DEGs.The results indicated that DEGs were mainly enriched in the various cancer pathways included microRNAs in cancer pathway.In addition, they were also significantly enriched in splicesome, regulation of actin cytoskeleton, synaptic vesicle cycle, calcium signaling pathway and neuroactive ligand-receptor interaction (Fig. 2).

Enrichment analysis of DEGs
We investigated the biological mechanisms of 50 ferroptosis-related genes mentioned above using the GO and KEGG pathway enrichment analyses (Fig. 3A,B).Analysis of Metascape demonstrated that these genes was significantly activated in the cellular response to chemical stress, autophagy, NRF2 pathway and Kaposi sarcomaassociated herpesvirus infection (Fig. 3C).It was shown that the biological processes were remarkably enriched in glutathione biosynthetic process, nonribosomal peptide biosynthetic process, response to iron ion, cellular www.nature.com/scientificreports/modified amino acid biosynthetic process and regulation of transcription from RNA polymeraseIIpromoter in response to stress, the cell components were enriched in autolysosome and secondary lysosome, the molecular function were enriched in NAD + or NADP + activity and estradiol 17-β-dehydrogenase activity (Fig. 3D).www.nature.com/scientificreports/

PPI network construction and analysis of DEGs
To further explore the interactions of DEGs, we imported the 50 genes into the STRING database and constructed a PPI network by the Cytoscape visualization software.We isolated 42 nodes and 122 edges and showed that the larger the degree of nodes, the redder the colour and the larger the points (Fig. 4A).Moreover, we used MCODE to obtain a closest cluster of PPIs (15 nodes, 49 edges, MCODE score ≥ 7) consisting of ATG3, ULK1, MAP1LC3A, LAMP2, MTOR, NRAS, HRAS, STAT3, NFE2L2, HMGB1, BRD4, RB1, PIK3CA, ATM, NF2 (Fig. 4B).Meanwhile, using the CytoHubba to select module and the top 10 genes were taken as the hub genes consisting of PIK3CA, ATM, RB1, BRD4, HRAS, NRAS, STAT3, MTOR, HMOX1, NFE2L2 (Fig. 4C).The PIK3CA, ATM, RB1, BRD4, HRAS, NRAS, STAT3, MTOR, NFE2L2 mined by the two modules were identical.Furthermore, the functional enrichment analysis for the cluster using Metascape and HiPlot showed that the hub genes were mainly involved in autophagy, phosphatidylinositol kinase activity, Kaposi sarcoma-associated herpesvirus infection and Thyroid stimulating hormone (TSH) signaling pathway (Fig. 4D,E).

Gene-miRNA network construction and analysis
We screened 9 hub genes and performed gene-miRNA network by using miRWalk 2.0 software.8 genes, except MTOR, were modulated and crosslinked by 238 miRNAs (Fig. 5A).The miRNAs with higher amounts of crosslinked target genes (≥ 2) are shown (Table 4).Functional enrichment analysis of 238 miRNAs by Funrich indicated that the molecular function was significantly enriched in transcription factor or regulator activity, protein serine/threonine kinase activity and a series of receptors or transporters activity (Fig. 5B).The biological pathways were enriched in the Integrin family cell surface interaction, TRAIL signaling pathway and proteoglycan syndecan-mediated signaling events (Fig. 5C).

ceRNA network construction and analysis
The corresponding lncRNAs of miRNAs were predicted with StarBase.Even though 15 miRNAs had higher amounts of cross-linked genes (≥ 2), only hsa-miR-1287-5p, hsa-miR-6893-3p and hsa-let-7i-5p had the upstream lncRNAs with higher reliability.The miRNAs cross-linked target lncRNAs and their type are shown (Tables 5,  6, 7).After cross-linking, 5 hub genes were predicted and 38 lncRNAs targeting 3 key miRNAs were revealed, 6 lncRNAs were predicted targeting two or three differential miRNAs included AL022311.1,AC016876.2,MALAT1, NEAT1, XIST and KCNQ1OT1.Based on the above, only NEAT1 served as a key lncRNA that predicted interacting with 3 miRNAs and 5 hub genes, the ceRNA network was constructed (Fig. 5D).www.nature.com/scientificreports/visibly upregulated and the BRD4 were slightly downregulated to the control.After the ferroptosis special inhibitor liproxstatin-1 added, the higher levels of PIK3CA and ATM with rotenone could be reversed as well as the lower levels of BRD4 (Fig. 6A-C).Interestingly, compared with the hub genes, the expression levels of NEAT1 were significantly upregulated in the neurotoxic groups and could not be rescued by liproxstatin-1 (Fig. 6D).Moreover, the NEAT1_2 levels in rotenone group also showed a higher trend and could be rescued by the liproxstatin-1 (Fig. 6E).We next determined the effects of stimulating the α-syn activation in rotenone models.The results showed that in comparison with control group, the expression of α-syn were increased in rotenone group.Liproxstatin-1 treatment also attenuated rotenone induced the aggregation of α-syn (Fig. 6F,G).To further confirm the ferroptosis in rotenone-treated SHSY5Y cells, GPX4 and FTH1 expressions, iron, MDA and     www.nature.com/scientificreports/ the expressions of GPX4 and FTH1 were significantly consistent decreased.Taken together, these results further confirmed that the NEAT1-PIK3CA/ATM ceRNA network may be activated and directly contacted with the ferroptosis induced by α-syn in PD models.

Discussion
Ferroptosis as a new cell death pathway has been reported the potential links to loss of Dopaminergic neurons and a number of ferroptosis features are commonly described in living or postmortem PD patients 13 .The aggregated of α-syn in PD pathology has long been investigated, but there is a strong rationale for involvement of the funtion of α-syn in ferroptosis recently.While iron can bind to α-syn oligomers in neurons induce irondependent oxidation 14 , the α-syn may also mediates iron metabolism by interacting with the related proteins.The mechanism in which α-syn is proposed to modulate iron import by facilitating the uptake of transferrin1-bound iron and upregulating the iron transport protein DMT1 15,16 .Besides, the aggregated of α-syn may mobilize the cellular iron stores through impairing ferritinophagy 17 .Furthermore, multiple studies report the lipid peroxidation related composition, such as PUFAs,4-HNE and AA, also link with the modification or formation of toxic oligomers 6 , suggesting a possible role of α-syn in lipid metabolism.However, the precise molecular mechanism of α-syn leading to dopaminergic neuronal death through ferroptosis remains unclear.
The present study used microarray data and bioinformatic techniques further explored the mechanisms of α-syn promoting ferroptosis in PD patients.We first obtained the 50 DEGs from the intersection of the datasets GSE49036 and FerrDb, including 27 upregulated genes and 23 downregulated genes.GO and KEGG enrichment analysis, eliminated the interferences, revealed that these DEGs were significantly involved in up or downstream whereby GSH and NRF2 pathway.The phosphatidylinositol kinase activity and the pathways related with lipid metabolism may play a pivotal role 18 .Moreover, the regulatory ceRNA network were successfully constructed, including 5 hub genes,3 miRNAs and 38 lncRNAs.The changes in those genes expression after rotenone and liproxstatin-1 by rescue were identified by qRT-PCR, suggested that the α-syn may influence ferroptosis in the ceRNA network and validated the reliability of the integrated data analysis.We also provided several genes that had not yet been mentioned in the area of PD and ferroptosis, which might be helpful to further illustrate the mechanism of synucleinopathy from the perspective of bioinformatics.
The results showed that the hub genes were significantly enriched in NRF2 signalling pathway which is considered to be an important regulatory factor for ferroptosis.Increased oxidative stress can promote the nuclear translocation of NRF2 and interacts with antioxidant response elements (AREs) resulting in transcriptional activation of ferroptosis moleculars, including NAD(P)H, HO-1, GPX4, SLC7A11/xCT and several iron transporters 19,20 .Our study also provided the evidence that the GSH synthesis and iron or lipid metabolsm related to NRF2 were enriched consistent with the above.Previous study indicated that Nrf2 KO with α-Syn stereotaxic delivery mice showed the dystrophic dendrites and increased neuroinflammation which is largely linked to the PD pathological feature 21 .Meanwhile, the NRF2 loss may enhances the α-Syn phosphorylation, which causes oligomerization and eventual convert the α-Syn into insoluble aggregates 22 .The discovery that there is a vicious cycle of α-syn overexpression, NRF2 suppression, and enhanced ferroptotic neuronal death that could be a key driver of parkinsonian phenotypes 23 .This may explain why the α-syn can accentuate ferroptosis of PD pathology via suppression of NRF2 and if there is a physical interaction between NRF2 and α-syn will be interesting to pursue.www.nature.com/scientificreports/ The recently affirmed role of selective autophagy such as ferritinophagy and lipophagy in driving cells toward ferroptotic demonstrated the accuracy of our results.The mechanism of the selective autophagy is the specific autophagy receptors function by interacting simultaneously with pidated LC3 or GABARAP subfamily proteins 24 .Excessive autophagy may impair lysosomal activity and promote ferroptosis, For example, the NCOA4-mediated ferritinophagy can induce the iron accumulation, the RAB7A-mediated lipophagy, SQSTM1-mediated clockophagy and HSP90-mediated CMA (D) may promote the lipid peroxidation in ferroptosis 25 .The interaction of them is also modulated by the upstream regulators of PIK3C3 complex contains the catalytic subunit of phosphatidylinositol 3-kinase 26 .Meanwhile, we found in the present study that the key modules identified by MCODE of PPI were significantly enriched in macroautophagy, autophagosome membrane or the regulation of autophagy and involved in the activity of phosphatidylinositol 3-kinase.As well known, the PI3K/Akt/mTOR signaling pathway was closely related with oxidative stress and the Akt can mediate the PI3K to play a key role in the pathogenesis of PD 27 .Abnormal α-syn aggregation participates in PD is associated with the involvement of TLR4/ PI3K/Akt/NF-κB signaling in neuroinflammation 28 .Therefore, the PI3K/Akt pathway may play an key role in ferroptosis following the α-syn stimulation.
The genes of interest were identified by the key module analyses with the highest MCODE scores and the top10 of CytoHubba, nine overlapped genes, included PIK3CA, ATM, BRD4, RB1, STAT3, NRAS, NFE2L2, HRAS and MTOR were selected as the most significant hub genes.PIK3CA is a major regulators of numerous protein kinases which plays key roles in glutamate oxidative stress-induced cell death in primary cerebrocortical neurons.It has been demonstrated that the protective roles of PIK3CA inhibitors against ferroptosis triggered by glutamate oxytosis and elucidated the potential therapeutic value in PD patients 29 .The DNA damage response serine/threonine kinase ATM was also found to be essential for iron metabolism.ATM induced ferroptosis by decreasing the expression of iron regulators involved in iron storage (ferritin heavy and light chain, FTH1 and FTL) and export (ferroportin, FPN1), moreover, an unexpected ATM-Ferritin/FPN1 regulatory axis as original determinants of ferroptosis through regulating labile iron levels were identified 30 .BRD4 belongs to the bromodomain and extraterminal domain (BET) protein family which can recognize acetylation sites and recruit transcription factors.A previous study suggested that targeting BRD4 may confer a clinical benefit in proliferation of cancer cells via the increase in ferritinophagy and enhanced the ferroptosis 31 , suppression of BRD4 may be associated with the pathogenesis of PD.RB1 and STAT3 are all signal key factors of transcription, in which RB1 is a corepressor of G1/S checkpoint while the STAT3 serves as a activator to enhance the inflammatory response.The negative status of two proteins promotes the occurrence of ferroptosis, especially STAT3, could improve the processes associated with ferroptosis by regulating SLC7A11 32,33 .NRAS and HRAS are the key components of RAS signaling pathways originating from the surface receptors.Emerging studies have provided evidence that the expression of RAS family related genes were significantly more susceptible to loss of neurons viability upon treatment with ferroptosis inducers 34 .Our results showed that the expression of PIK3CA, ATM, RB1, STAT3 and NFE2L2 was upregulated, while BRD4, NRAS, HRAS and MTOR expression was downregulated.The PIK3CA or ATM belongs to the drivers of ferroptosis in the FerrDB database and the BRD4 belongs to the suppressor.These data suggest that the three hub genes may play a significant role in ferroptosis upstream that caused by the α-syn.
The miRNAs are short chain non-coding RNA that can inhibit or degrade the translation of the downstream target genes.In our study, we found eight hub genes could successfully forecast 238 miRNAs and also identified 15 miRNAs targeting at least two genes.Among these miRNAs, only hsa-miR-1287-5p, hsa-miR-6893-3p and hsa-let-7i-5p could predict the relevant lncRNAs.Previous studies on these miRNAs and PD were limited, miR-1287-5p plays critical roles in multiple cancers such as hepatocellular carcinoma, osteosarcoma and breast cancer that involved in pro-inflammatory and antioxidant pathway like GPX4 35 .The miR-6893-3p may regulate IL-6, IGF-1 and VEGF in multiple myeloma or monoclonal gammopathy 36 .Several reults demonstrated the hsa-let-7i-5p could target the fatty acid and lipid metabolism pathway 37 which consistent with our findings on interaction with PIK3CA.Funrich enrichment analysis indicated that the molecular function was significantly enriched in transcription factor or regulator activity, protein serine/threonine kinase activity and a series of receptors or transporters activity.The biological pathways enriched included the Integrin family cell surface interaction, TRAIL signaling pathway and proteoglycan syndecan-mediated signaling events.Further researches should be carried out to pinpoint the unique molecular mechanism.
lncRNAs refers to a type of long length noncoding RNA that affect the downstream miRNAs but not translate into proteins.In this study, five hub genes (PIK3CA, ATM, BRD4, STAT3, NRAS) and three miRNAs (hsa-miR-1287-5p, hsa-miR-6893-3p, hsa-let-7i-5p) were successfully constructed in the ceRNA network.They were highly correlated with 38 lncRNA and six of them targeting at least two miRNAs, included AL022311.1,AC016876.2,MALAT1, NEAT1, XIST and KCNQ1OT1.The nuclear paraspeckle assembly transcript 1 (NEAT1) has been showed to be involved in the regulation of cell cycle, proliferation and migration of tumor, apoptosis and ferroptosis.NEAT1 has two transcripts contains NEAT1_1 (3.7 kb) and NEAT1_2 (23 kb) which act as sponges to absorb miRNA thereby promoting the target genes express 38 .A previous research found that NEAT1 could competitively bind miR362-3p and thus leads to inhibition of Myo-inositol oxygenase (MIOX), the axis was regulated by the p53 and thereby increasing the sensitivity to ferroptosis 39 .Our findings identified that only NEAT1 predicted three target miRNAs and served as a maincenter linking five key genes together.It may be potentially used as a novel biomarker of α-syn promoting ferroptosis in PD patients.
In order to further affirm our results, we performed experiments in vitro by generating the PD models with neurotoxins.We used SH-SY5Y cells and developed by introducing rotenone to validate the integrated data.The qRT-PCR analysis in our own samples showed that PIK3CA and ATM were significantly upregulated in rotenone group, whereas BRD4 were downregulated.Research studies have manifested that in the rotenone model, α-Syn were observed in the surviving dopaminergic neuron, however, the MPP + and other neurotoxin models lack the formation of Lewy Bodies 40 .Our western blot analysis consistent showed that the rotenone could stimulate the aggregation of α-Syn thereby causing the ferroptosis.Furthermore, liproxstatin-1, a specific ferroptosis inhibitor, can ameliorate the cells injury and reduce the expression of PIK3CA and ATM.It may be that ferroptosis inhibitor could suppress the toxic α-syn oligomers via the mechanism of aggregate-membrane interaction 8 .Therefore, PIK3CA and ATM may serve as new hub genes for α-Syn promoting ferroptosis.We also found the key lncRNAs NEAT1, whether NEAT1_1 or NEAT1_2, was upregulated in the rotenone groups.The expression of NEAT1, especially NEAT1_2, in rotenone group were significant higher compared with the control.Together, we proposed α-Syn as a key contributor to ferroptosis progression which may involve in NEAT1-PIK3CA/ATM ceRNA network.These demostrated the hub genes and the ceRNA may act independently on ferroptosis but not other cell death mechanism and broaden new insights that targeting the ceRNA network can be developing strategies for PD therapy (Supplementary Figure ).
Despite the significant results obtained in our study, there were several limitations.First, our data are not from high-throughput sequencing but from single microarray databases, the number of samples was not large enough which may contribute to a high false-positive rate and one-sided results.Moreover, although we analyzed the expression of related genes and lncRNAs in PD cell models, the relationship between them and other unknown molecules are still undetected, single-cell RNA-Seq or other methods are necessary for verifying our results and illustrating the roles of PIK3CA, ATM and NEAT1 in PD.In addition, functional experiments based on transgenic PD models in vivo and induced pluripotent stem cells from PD patients are need to further verify the expression of current genes protein level in related mechanisms.
In summary, our current study provided a reliable comprehensive analysis of combining the α-syn and ferroptosis based on the gene expression profiling.The identified DEGs revealed core networks and molecules significantly associated with the α-syn aggregation driving ferroptosis in SN of PD patients.The above view were all strongly correlated with NRF2, PI3K and autophagy signaling pathway.Our results will help enhance the cognition of the internal connection between the two independent death patterns, ferroptosis and autophagy.We believe that the genes verified by our findings may become specific diagnostic markers as well as provide promising targets for genetic therapy for PD patients in the future.

Figure 1 .
Figure 1.Identification of DEGs.(A)Volcano plot of DEGs between PD patients and controls (PD vs Control).(B) A heatmap of top 100 DEGs.(C) Venn diagram of ferroptosis DEGs by the GSE49036 intersected with FerrDB.

Figure 2 .
Figure 2. GSEA analysis of DEGs.(A-G) The result of the enrichment gene dataset analysis with the WebGestalt software.(H) The bar chart of the top biological pathways.

Figure 3 .
Figure 3. Enrichment analysis of DEGs.(A,B) Network of enriched terms with Metascape.software.(C) The bar chart of the top biological pathways by the Metascape software.(D) Go enrichment analysis of biological process, cell components and molecular function.

Figure 4 .Figure 5 .
Figure 4. PPI network construction and analysis of DEGs.(A) Cytoscape network visualization of the 50 genes according to the STRING online database.(B) Cluster of the PPIs identified by the Mcode.(C) Select module and the top 10 genes by the CytoHubba.(D) The bar chart of the biological pathways in hub modules.(E) Go-KEGG enrichment analysis of the hub genes.

Table 1 .
Specific primers used for quantitative qRT-PCR.

Table 4 .
The miRNAs with higher amounts of cross-linked target genes.