A Parkinson's disease CircRNAs Resource reveals a link between circSLC8A1 and oxidative stress

Abstract

Circular RNAs (circRNAs) are brain-abundant RNAs of mostly unknown functions. To seek their roles in Parkinson's disease (PD), we generated an RNA sequencing resource of several brain region tissues from dozens of PD and control donors. In the healthy substantia nigra (SN), circRNAs accumulate in an age-dependent manner, but in the PD SN this correlation is lost and the total number of circRNAs reduced. In contrast, the levels of circRNAs are increased in the other studied brain regions of PD patients. We also found circSLC8A1 to increase in the SN of PD individuals. CircSLC8A1 carries 7 binding sites for miR-128 and is strongly bound to the microRNA effector protein Ago2. Indeed, RNA targets of miR-128 are also increased in PD individuals, suggesting that circSLC8A1 regulates miR-128 function and/or activity. CircSLC8A1 levels also increased in cultured cells exposed to the oxidative stress-inducing agent paraquat but were decreased in cells treated with the neuroprotective antioxidant regulator drug Simvastatin. Together, our work links circSLC8A1 to oxidative stress-related Parkinsonism and suggests further exploration of its molecular function in PD.

Keywords: AGO2; CircRNA; CircSLC8A1; Parkinson's disease; RNA sequencing.

Figure EV1. Tissues retrieved and experimental in‐house RNA sequencing methodology

1. A

   Brain region origin of tissues retrieved from the NBB detailing each brain area.

2. B

   Average RIN plot of RNA prepared from the frozen tissues.

3. C

   Library preparation key steps (see details in Materials and Methods).

4. D

   Library plot in a TapeStation run.

5. E–G

   Number of circRNAs detected in each sample of all 3 regions. Samples with < 5,000 circRNAs were removed from the analysis.

6. H

   Number of samples from PD and control from each tissue that were included in the analysis.

7. I

   Number of total reads and circRNAs from each tissue.

Figure 1. Brain region‐specific mRNA profiles indicate multicellular PD‐related processes

1. A

   Brain region origin of tissues: SN (blue), MTG (green), and AMG (red).

2. B

   Non‐supervised clustering heatmap plot indicates substantial clustering of transcripts from each tissue but not of transcripts from patients compared to controls. PD samples indicated in red.

3. C

   PCA clustering of all RNA molecules in the SN control (green) and PD samples (blue).

4. D

   qPCR validation of TH levels in SN samples of PD and controls, normalized to beta‐actin mRNA. t‐test *P = 0.025. Data presented as mean ± SD. n = 18 for CT and 24 for PD.

5. E

   Top enriched GO pathways of DE genes in the SN of PD vs controls and the corresponding logP values, Wald test (DEseq2 analysis).

6. F, G

   Venn diagrams demonstrating cell type‐specific DE genes in the SN of PD vs control tissues, up/down arrows indicate gene groups that were up/down regulated.

7. H

   Module–trait relationship of WGCNA analysis showing correlation and corrected P values (in brackets) for each gene module (indicated by colors) as related to the different external traits selected (disease condition, brain region, age, and sex). Correlation test p value calculated by WGCNA package.

8. I

   Module membership vs. gene significance in selected modules as reflecting external sample traits: disease condition (reflected in the brown and red modules), brain region (the blue module), and sex (the cyan module), regression‐based p value calculated by WGCNA analysis. n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control, and 13 for MTG PD, 10 for SN control and 15 for SN PD.

9. J

   Significant GO terms that emerged as associated with the blue, red, and brown modules and the corresponding log P values, Fisher's exact P value.

Figure EV2. PCA analysis of all tissues and of PD vs control samples, WGCNA analysis QC, clustering of the samples and modules

1. A

   PCA of all mRNA samples from all tissues.

2. B

   Substantia nigra sample heatmap, relative expression of each gene is represented by colors, marked by the color key: Reduction in expression is represented by light green and elevation is represented by blue.

3. C, D

   PCA plot of control vs PD samples in the amygdala and MTG.

4. E

   Sample dendrogram showing sample clustering and the different classification patterns of external traits for each sample.

5. F

   Cluster dendrogram of all genes.

6. G, H

   log₂ fold change and −log P value of known DE splicing factors in PD SN vs control SN. n = 10 for SN control and 15 for SN PD. Data presented as mean ± SD, Walt test (DEseq2 analysis).

Figure 2. CircRNA abundance is brain region‐dependent and PD‐modified

1. A

   mRNA processing either involves unidirectional splicing, creating linear mRNA transcripts, or back‐splicing generating circRNAs.

2. B, C

   Venn diagram (not to scale) showing brain region‐ and disease‐characteristic circRNA abundance patterns in analyzed brain regions and disease state.

3. D

   Heatmap of circRNA abundance, indicating significant tissue‐based SN clustering, PD samples indicated in red.

4. E, F

   Pie chart demonstrating brain‐specific expression of mRNAs (E) in comparison to circRNAs (F).

5. G

   Venn diagram demonstrating predicted cell‐type characteristics of circRNA host genes in the SN.

Figure EV3. PCA analysis of circRNA profile

1. A–D

   PCA analysis of circRNA profiles in samples according to specific tissues, age, condition, and gender.

2. E

   CircRNA abundance does not change in a gender‐specific manner in the tested tissues. n = 23 for females and 46 for males. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

Figure 3. CircRNA abundance is reduced in the PD SN and negatively correlates to RNA editing levels

1. A

   Total circRNAs detected in the SN, MTG, and AMG of both PD and control donors, *P = 0.02 *P = 0.026 and *P = 0.022 for comparison between PD and CT in SN, MTG,and AMG. t‐test for total SN vs MTG *P = 2.27E‐8 and SN vs AMG *P = 5.86E‐10. n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values. P value was calculated by t‐test.

2. B

   Global editing levels in the 3 brain regions are based on the Alu editing index, representing the weighted average editing level across all expressed Alu sequences, Wilcoxon, corrected *P = 0.038, 0.055 and 0.42 for amygdala, MTG and SN. Wilcoxon test between the Alu editing index of control and PD samples.

3. C

   Global RNA editing levels based on the Alu editing index, studied only in Alu elements within circRNA exons and their flanking introns, in the 3 brain regions, for circRNAs corrected *P = 0.025 and 0.040 for the AMG and MTG and P = 0.3 for SN. For flanking introns, *P = 0.02 and 0.040 for AMG and MTG and P = 0.34 for SN, Wilcoxon test. n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

4. D

   Negative correlation between total RNA editing and circRNA abundance, correlation = −0.37, correlation test p value = 0.0015.

5. E

   Adar1 expression levels in all tissues. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

6. F, G

   Correlations between circRNAs detected in each sample and the age of the control (green, correlation = 0.68, P = 0.032) and PD donors (blue, correlation = −0.14, P = 0.62).

7. H

   Volcano plot of DE circRNAs in PD vs control tissues, red dots indicate statistically significant DE circRNAs according to FDR correction of Wald test (DEseq2 analysis). CircSLC8A1 is marked is red.

Figure EV4. RNA Editing in PD sample and controls

1. A

   Hyper‐editing levels in all 3 tissues, PD, and controls, n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

2. B

   Global editing levels measured by Alu editing index comparing the three different tissues in PD or control, n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

3. C

   Alu editing index showing that A‐to-G mismatch is the most dominant, as other mismatches are a result of noise or other biological mechanisms and should be much lower than the A‐to-I (viewed as A‐to-G) RNA editing signal, n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

4. D

   Alu editing index in circRNAs and flanking introns, in 3 different intron windows: 1,000, 5,000, and 10,000 nucleotides, n = 8 for amygdala control, 15 for amygdala PD, 8 for MTG control and 13 for MTG PD, 10 for SN control and 15 for SN PD. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

Figure 4. CircSLC8A1 increases in the PD brain and upon PQ exposure are accompanied by suppression of the SLC8A1 protein

1. A

   Graphic representation of circSLC8A1 and its post‐transcriptional regulation in PD and under PQ exposure.

2. B, C

   CircSLC8A1 and linear SLC8A1 levels in the SN of control and PD brains, t‐test P = 0.025, P = 0.118, respectively.

3. D

   qPCR validations of circSLC8A1 and SLC8A1 mRNA changes in the SN, t‐test *P = 0.025. n = 18 for CT and 24 for PD, x defines mean values.

4. E

   No difference in circSLC8A1 levels in females and males, t‐test P = 0.583 for control and P = 0.262 for PD. n = 15 females and 16 males.. The box is drawn from Q1 to Q3 with a horizontal line drawn in the middle to denote the median and x marks the average. Whiskers mark minimum or maximum values.

5. F, G

   CircSLC8A1 and linear SLC8A1 mRNA levels are correlated in control, R ² = 0.47 t‐test P = 0.0017 but not in PD samples, R ² = 0.016, P > 0.05. n = 18 for CT and 23 for PD. Data presented as mean ± SD. Note that y‐axis scales are different to allow better visualization.

6. H

   Increased CircSLC8A1 RNA levels in PQ‐exposed SH‐SY cells (normalized to Tubb3 and RPL19) while other circRNAs remain unchanged, t‐test *P = 0.028 and *P = 0.0215 for 25 μM and 50 μM, respectively. n = 4 biological replicas for each condition. Data presented as mean ± SD.

7. I, J

   Immunoblot decline (Western blot and quantification) of the SLC8A1 protein in PQ‐exposed SH‐SY cells, with Tubulin as a loading control, t‐test *P = 0.015. n = 3 biological replicas for each condition. Data presented as mean ± SD.

Figure 5. CircSLC8A1 is reduced under neuroprotection

1. A

   qPCR measurements of circSLC8A1 and SLC8A1 mRNA after SH‐SY treatment of statins or the PF LRRK2 inhibitor, t‐test **P = 0.007 and **P = 0.01 for PF‐06447475 and *P = 0.047 and *P = 0.0136 for statins for cirSLC8A1 and SLC8A1 mRNA, respectively, n = 3 biological replicas for each condition. Data presented as mean ± SD.

2. B, C

   Protein gel and quantification of SLC8A1 and Tubulin as a loading control in statin and PF‐treated neuronal cell cultures, t‐test *P = 0.027 for statins and *P = 0.011 for PF‐06447475, n = 3 biological replicas for each condition. Data presented as mean ± SD.

3. D

   SLC8A1 immunostaining of statin‐treated SH‐SY cells.

Figure 6. Cytoplasmic circSLC8A1 binds Ago2 and regulates miR‐128 targets

1. A

   Gene structure and conservation of the human and murine SLC8A1 gene. The box marks the second exon that circularizes to create circSLC8A1.

2. B

   qPCR quantification from nuclear/cytoplasmic fractionations of SH‐SY cells, t‐test P = 0.005 for circSLC8A1, P = 0.009 for NEAT1 and P = 0.009 for Rpl19. n = 3 biological replicas for each condition. Data presented as mean ± SD.

3. C

   qPCR quantification of Ago2‐bound RNA, immunoprecipitated (IP)/input values normalized to negative control for a known Ago2 target, circCDR1as, and circSLC8A1, t‐test P = 0.029 for circSLC8A1 and P = 0.001 for circCDR1as, n = 3 biological replicas for each condition. Data presented as mean ± SD.

4. D

   Schematic representation of circSLC8A1 with Ago2‐binding sites from CLIP‐seq and predicted miRNA binding sites.

5. E

   miRNAs predicted to bind circSLC8A1 with 3 sites and more.

6. F

   Schematic representation of circSLC8A1‐Ago2 complex binding and its effect on miR‐128 targets.

7. G

   Relative expression of known miR‐128 targets from SN RNA‐seq, **P = 0.0004 for BMI1, *P = 0.011 for SIRT1 and *P = 0.005 for AXIN1.

8. H

   Relative expression of known miR‐128 targets in SH‐SY cells treated with Paraquat, t‐test *P = 0.03, *P = 0.02 and *P = 0.002 for AXIN1, BM1 and SIRT1, respectively. Data presented as mean ± SD with 4–6 biological replicas for each condition.

9. I

   miR enrichment analysis of DE genes in SN from validated miR‐targets.

Figure EV5. Gene expression of circSLC8A1 knockdown

1. A

   circSLC8A1 expression after siRNA treatment (t‐test *P = 0.0075), n = 3 biological replicas for each condition. Data presented as mean ± SD.

2. B

   siRNA targeting circSLC8A1 does not affect SLC8A1 mRNA expression after siRNA treatment (t‐test P > 0.05). n = 6 for each condition. Data presented as mean ± SD.

3. C

   Volcano plot of gene expression profile of siCircSLC8A1 compared to siControl, blue dots indicate corrected P < 0.05, FDR correction of Wald test (DEseq2 analysis).

4. D

   Enriched Go terms of DE genes in siCircSLC8A1 compared to siControl analysis, FDR correction Wald test (DEseq2 analysis).