CircRNA are differentially expressed in psychiatric disorders
We identified a total of 8,762 high confidence circRNAs in all the 59 samples (Supplementary Table S2). To explore the implication of PBMC circRNA in the psychiatric disorders we performed differential expression in SZ and BD versus HC individuals. To control for the potentially confounding impacts of sex and age, these factors were included as covariates in our model. We also restricted analyses to the broadly expressed circRNAs, as defined by those circRNAs that were expressed in at least 10 samples for each comparison analysis (SZ vs HC and BD vs HC), to ensure the reliability of findings. A total of 22 and 33 circRNAs were found to be significantly changed in SZ and BD, respectively, as compared to the healthy control group (P <0.05, FC >2) (Fig. 1A; Supplementary Table S3). Notably, 71% of these circRNAs were downregulated, including 14 circRNAs in SZ and 26 in BD group. There were two circRNAs, circNCF1 and circLINC00969 common between the SZ and BD differential expression lists. Annotation of the DE circRNAs showed that 49 circRNAs are exonic, 2 intronic and 2 intergenic circRNAs. Also, 45 of the DE circRNAs are derived from protein coding genes, while 2 and 3 are from lincRNA and pseudogenes, respectively. In addition to the conventional case control analysis, we also combined the psychosis groups and compared to controls and observed three differentially expressed circRNAs including circUGP2, circNCF1C and circKLHDC4 (P <0.05, FC >2) (Supplementary Table S3). In addition, we compared the circRNA expression between the psychosis groups (SZ vs BD) and observed six differentially expressed circRNAs including four that were upregulated and two downregulated in SZ, one of which, circNCF1, was also differentially expressed in the SZ/BD vs HC groups analyses (Supplementary Table S3).
We next determined whether the ratio of circRNA over linear RNA is changed in psychiatric conditions. A linear model including sex and age as covariates was used for the circRNA subset used in DE analysis. The analysis showed there was 33 circRNAs with ratios significantly changed in SZ, compared to the healthy controls (P <0.05, FC >1.5) (Fig. 1B; Supplementary Table S3). However, no significant alteration of the ratio was observed in BD group.
To validate the results of circRNA sequencing, we randomly selected eight significantly dysregulated circRNAs in SZ and BD and conducted qPCR using PBMCs from 21 SZ and 21 BD patients compared to 21 sex- and age-matched healthy controls to analyse the expression levels. Our data confirmed that the tested circRNAs were differentially expressed in psychiatric patients compared to healthy individuals, indicating that the sequencing results are reliable (Fig. 2).
Functional significance of dysregulated circRNAs
CircRNA have emerged as active molecules with biological significance that contribute to the key regulatory programs through interacting with miRNA and RNA binding proteins (RBP), or by translating into peptides [16,17]. In order to gather insight into potential function of the dysregulated circRNAs, we investigated their association with miRNA and RBP, and also their potential for hosting translation elements. We found 45 circRNAs contain at least two binding sites for miRNA which was supported by two prediction tools, miRanda [35] and TargtScan [36] that suggest these molecules could function as competing endogenous RNA for the target mRNA. Collectively we identified 351 circRNA-miRNA interactions, including 148 unique miRNAs (Supplementary Table S4). Interestingly, we found targets of these miRNA were significantly enriched in neurodevelopmental and neurobehavioral disorders including SZ, BD, autism, and mental retardation. GO analysis of circRNA interacting miRNA target genes showed an enrichment of neuronal terms such as neurogenesis, synapse, and neuronal differentiation and development (Supplementary Table S4). We also detected RBP sites within all the circRNAs, with a minimum of two to hundreds of binding sites for each circRNA (Supplementary Table S4). The most frequent circRNA binding factors were ESR1, TRPS1 and BRD4. While our analysis predicted six circRNAs have the potential to translate into proteins (coding probability >0.9) with open reading frames between ~300-4,000 bp (Supplementary Table S4), we also identified 20 differentially regulated circRNA that have been observed to be m6A modified [47]. This is potentially functionally significant given that m6A modification of circRNAs has been suggested to enhance the initiation of protein translation [47] (Supplementary Table S4). These analyses suggest that the observed alteration in circRNAs could have functional significance in the PBMCs, which may contribute to the psychopathological conditions.
Alternative splicing profiles in diagnosis groups
We investigated the internal structure of all the identified circRNA by looking at the alternative patterns of splicing of circRNA exons in SZ, BD and HC groups. Our results revealed a total of 2,658 AS events occurred in all the groups (883 unique events) (Supplementary Table S5). We found 545 AS events in the disease groups that were not present in the healthy controls, of which 150 events were common between SZ and BD. There were also 120 events exclusively occurred in the HC (Fig. 3A). However, the composition of the AS events showed a similar pattern across the groups, with the A3SS showed the highest and IR the lowest fractions (Fig. 3B). Next, we analysed the relative abundance of the events across the studied groups using “percentage spliced in” (PSI) values calculated by the CIRI-AS [33], and the heatmap results showed the top 26 AS cirexons were broadly expressed across the diagnosis groups, while the remaining AS cirexons were sporadically present in certain samples (Fig. 3C). Differential analysis using PSI showed four and three AS events were significantly altered in SZ and BD, respectively, as compared to HC (P<0.05) (Supplementary Table S5). Collectively, these results suggested distinct circular alternative splicing patterns in the diagnosis groups and also their potential function in regulatory pathways in PBMCs.
Validation and characterization of the identified circRNA
We compared all the identified circRNA with ~140,000 cirRNAs in circBase [48], and also with ~100,000 circRNAs detected in our previous studies [49,29] and found an additional 754 circRNAs not previously annotated in either database. To validate the authenticity of the identified circRNAs, we randomly selected 15 circRNAs with different expression levels, mostly from the newly identified candidates (circKLRC4, circFAM107B, circSKAP1, circUTRN, circATP8B4, circPAK2, circABCA13, and circTMEM154), and performed qPCR using specific divergent primers (Supplementary Table S1). The results validated all the 15 tested circRNAs that showed 10- to 50-fold more resistance compared to the linear transcripts, following the RNase R treatment (Fig. 4A). In order to glean insight into the possible biological significance of all the detected circRNAs, we performed pathway overrepresentation analysis for their parental genes and the results indicated a significant enrichment of pathways that are relevant to immune system such as B cell activation, T cell activation, and interleukin signalling (Fig. 4B). The circRNAs were predicated to colocalized in multiple cellular locations, with more circRNAs observed in cytosolic (7,162) followed by nuclear and membrane fractions, in contrast to chromatin that contain the least number of circRNAs (1,248) (Supplementary Table S6). Analysis of the back-splicing frequency showed 167 circRNAs with a circular to linear reads (CLR) > 0.66, indicating these circRNAs had higher expression levels than their linear counterparts. (Supplementary Table S6). As shown in the heatmap plot (Fig. 4C), most of these circRNAs showed high CLR across the studied groups, suggesting their potential biological significance in gene expression regulation.