Dysregulation of circRNA expression in the peripheral blood of individuals with schizophrenia and bipolar disorder

Circular RNAs (circRNAs) are head-to-tail back-spliced RNA transcripts that have been linked to several biological processes and their perturbation is evident in human disease, including neurological disorders. There is also emerging research suggesting circRNA expression may also be altered in psychiatric and behavioural syndromes. Here, we provide a comprehensive analysis of circRNA expression in peripheral blood mononuclear cells (PBMCs) from 39 patients with schizophrenia and bipolar disorder as well as 20 healthy individuals using deep RNA-seq. We observed systematic alternative splicing leading to a complex and diverse profile of RNA transcripts including 8762 high confidence circRNAs. More specific scrutiny of the circular transcriptome in schizophrenia and bipolar disorder, compared to a non-psychiatric control group, revealed significant dysregulation of 55 circRNAs with a bias towards downregulation. These molecules were predicted to interact with a large number of miRNAs that target genes enriched in psychiatric disorders. Further replication and cross-validation to determine the specificity of these circRNAs across broader diagnostic groups and subgroups in psychiatry will enable their potential utility as biomarkers to be established. • We identified 8762 high confidence circRNAs with systematic alternative splicing in human PBMCs. • CircRNAs were dysregulated in schizophrenia and bipolar disorder, compared to a non-psychiatric control group. • The DE circRNAs were predicted to interact with miRNAs with target genes enriched in psychiatric disorders. • Some circRNAs have the potential to serve as biomarkers in psychiatry.


Introduction
Schizophrenia (SZ) and bipolar disorder (BD) are complex neuropsychiatric conditions with a range of severe symptoms, usually emerging in adolescence and early adulthood, such as delusion, hallucination, cognitive deficit, depression and mania [1,2]. While they have unique aspects, SZ and BD share common clinical features, genetic determinants and treatment response. Evidence from family, twin and adoption studies has supported that both SZ and BD exhibit strong genetic components and high heritability at 60-80%, with overlapping genetic risk and gene expression profiles [3,4]. Genetic and epigenetic changes that affect RNA level in particular can play a substantial role in the etiological pathways by directly modifying the mRNA template or by altering posttranscriptional regulation of gene expression. Non-coding RNAs including small microRNA and longer lncRNA play an important part in regulating gene expression and appear to be altered by the genetic and environmental exposures associated with psychotic disorders, suggesting they may be essential for normal physiology and homeostasis of the nervous system [5][6][7].
While the biological significance of lncRNA is now well established, a novel subclass of covalently closed circular lncRNA or circRNA has also emerged [8]. This conserved family of RNA transcripts is generated by an alternative splicing mechanism, known as back-splicing, in which a donor splice site is joined to an acceptor splice site. These molecules represent several specific features that make them distinguishable in mammalian cells. CircRNAs can be spliced from a combination of host exons resulting in the production of different variants from a single gene [9]. The expression profile of the circRNA has demonstrated a cell type-, tissue-and stage-specific pattern [10][11][12]. These molecules can function as competitive endogenous RNAs (ceRNAs) by sponging miRNA molecules to control their accessibility, thus indirectly regulating their target's expression [13,14]. They can also directly regulate RNA transcription by associating with RNA polymerase II transcription complexes [15]. Their ability to sequester RNA-binding proteins (RBPs) and form RNAprotein complexes allows circRNA to be involved in intracellular localization and transport of RBPs and associated mRNAs [16,17]. More recently, some circRNAs were shown to be capable of translating into functional proteins [18]. There is growing evidence suggesting circRNA involvement in the pathogenesis of human diseases, such as autoimmunity [19], cancer [20], cardiovascular diseases [21] and central nervous system disorders including SZ, major depressive disorder (MDD), multiple sclerosis (MS) and Alzheimer's disease (AD) [22]. Recent studies have shown perturbation of these molecules in human fluids including peripheral blood and suggest that some circRNAs might be potential candidates for the prediction of disease states or treatment responses [23]. For example, 406 circRNAs were dysregulated in peripheral blood of multiple sclerosis patients with two displaying predictive validity [24]. Similarly, a global alteration of circRNAs was reported in amyotrophic lateral scleros i s ( A L S ) , a m o n g w h i c h h s a _ c i r c _ 0 0 2 3 9 1 9 , hsa_circ_0063411 and hsa_circ_0088036 were supported as potential biomarkers [25]. CircRNA_103636 expression was associated with successful antidepressant treatment in blood of individuals with MDD, and therefore maybe predictive of treatment response [26].
CircRNA has also been implicated in the regulation of neuronal activity and plasticity, synaptic transmission and neuronal depolarization, which are all critical processes in pathophysiology of psychiatric disorders [27][28][29]. While circRNA appears to be dysregulated in the brain, the genetic and environmental exposures that give rise to this may also manifest in peripheral tissues, such as blood. As these cells are highly accessible through liquid biopsy, their circRNA profile could provide useful insights into aspects of the molecular pathology of psychiatric conditions that could guide further advances related to their diagnosis and treatment. In this study, we provide a comprehensive analysis of circRNA expression patterns and characteristics in the PBMCs of SZ and BD, as well as in healthy controls (HC), and examine their dysregulation in these psychotic disorders.

Sample collection
Participants provided signed informed consent. This study was approved by the Human Research Ethics committees of the University of New South Wales (HC12384), St Vincent's Hospital (HREC/10/SVH/9) and the South East Sydney and Illawarra Area Health Service (HREC 09/081). PBMCs were collected from 20 participants with schizophrenia, 19 with bipolar disorder and 20 non-psychiatric controls following the application of exclusion criteria and quality control processes. (For additional details, see Supplementary Information.)

RNA extraction, library preparation and RNA-seq
The PBMC samples were RNA extracted and quality controlled, and 300 ng RNA was subjected to rRNA depletion followed by fragmentation and library preparation before 100 cycles of paired-end sequencing (2×100 bp). (For additional details, see Supplementary Information.)

CircRNA detection
Sequencing reads were quality checked and subjected to CIRIquant pipeline [30] to predict circRNA species. A total of 108,865 circRNAs were initially identified in all samples; however, in order to obtain high confidence circRNAs, we used a filtering cut-off minimum of two junction reads in at least five samples, which allowed a minimum of 10 back-splice junction reads (BSJ) per circRNA. This criterion resulted in 8762 unique circRNAs in all samples, and we used these high confidence circRNAs for all the analyses performed in this study.

Differential expression analysis
Differential expression (DE) of circRNA was performed using DE pipeline in CIRIquant that implements a modified version of edgeR [31] fed with BSJ reads and total mapped reads. We used sex and age as covariates. To detect significant differences in the ratio of circRNA to linear RNA between the conditions, we used a linear model adjusted for sex and age using the Voom package [32]. (For additional details, see Supplementary Information.) Detection of circular to linear RNA ratio (CLR) We used CIRIquant [30] to obtain the circular to linear transcribes ratio using the forward-splice junction reads (FSJs) and BSJs,and the following formula: 2*bsj/(2*bsj+fsj). (For additional details, see Supplementary Information.)

cDNA synthesis and quantitative PCR
The RNase R-treated total RNA was reverse transcribed with random hexamers. Real-time PCR was performed using divergent primers. (For additional details, see Supplementary Information.) We used ΔCt method to normalize the data using GAPDH and HMBS as internal references. Student's t-test (one-tailed) was used to compare expression levels. Primer sequences are provided in Supplementary Table S1.

Detection of alternative splicing
We applied CIRI-AS pipeline [33] to identify alternative splicing (AS) events using paired-end reads. This approach employs a de novo assembly algorithm which is able to specifically detect circRNA exons and AS events. (For additional details, see Supplementary Information.)

Prediction of circRNA interactions
Sequences of the circRNAs were obtained from the circAtlas 2.0 database [34]. Then, miRNA-binding sites were predicted using miRanda [35] and TargtScan [36], and circRNAs predicted by both tools were used for circRNA-miRNA interaction analysis. The potential associations of the circRNA and RBP were determined by the "network" tool from circAtlas [34]. miRNA target gene prediction and enrichment analysis are detailed in Supplementary Information.

Protein-coding potential, modification and disease involvement
To analyse the protein-coding features of the circRNAs, we used the Circbank database [37]. This was also used to identify the m6A modification within the circRNA transcripts. We also used circRNA disease databases circ2disease [38], circad [39], circAtlas [34] and circRNADisease [40] to intersect our deferentially expressed circRNAs to those reported in previous studies.

Cellular localization prediction and functional enrichment
To predict cellular localization of all the identified circRNAs, we used the CSCD database [41]. Overrepresentation analysis was conducted using Gene Ontology (GO) annotation [42,43]. Pathway overrepresentation analysis was additionally run using PANTHER pathway database version 15.0 [44]. (For additional details, see Supplementary Information.)

Visualization
All the plots and graphs were generated in R. Venn diagrams were produced with the VennDiagram package [45]. Heatmaps were constructed using the heatmap.2 function from the gplots package (http://CRAN.R-project.org/ package=gplots). The remaining plots and graphs were generated with ggplot2 package [46].

CircRNA are differentially expressed in psychiatric disorders
We identified a total of 8762 high confidence circRNAs in all the 59 samples (Supplementary Table S2). To explore the alteration of PBMC circRNA in psychotic 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 in SZ and 26 in the BD group. There were two circRNAs, circNCF1 and circLINC00969, in 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 them to controls, and observed three differentially expressed circRNAs including ci rcUGP2, 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 upregulated and two downregulated in SZ, one of which, circNCF1, was also differentially expressed in the SZ/BD vs HC analyses (Supplementary Table S3).
We next determined whether the ratio of circRNA over linear RNA (CLR) is changed in psychiatric conditions. A linear model including sex and age as covariates was used for the circRNA subset used in the DE analysis. The analysis showed there were 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 has emerged as active molecules with biological significance that contribute to the key regulatory programmes 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 that were supported by two prediction tools, miRanda [35] and TargtScan [36], suggesting these molecules could function as competing endogenous RNA for 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  Table S4), we also identified 20 differentially regulated circRNAs that have been observed to be m 6 A modified [47]. This is potentially functionally significant given that m 6 A modification of circRNA 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 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 2658 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 that exclusively occurred in the HC (Fig. 3a). However, the composition of the AS events displayed a similar pattern across the groups, with the 3' AS site (A3SS) showing the highest and intron retention (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). In differential splicing analysis using PSI, we observed four significantly altered AS events in SZ and a further three in BD, compared to the control group (P<0.05) (Supplementary Table S5). Collectively, these results suggested there were distinct circular alternative splicing patterns in the diagnosis groups and that they may have potential function in regulatory pathways in PBMCs.

Validation and characterization of the identified circRNA
We compared all the identified circRNAs with~140,000 cirRNAs in circBase [48] and~100,000 circRNAs detected in our previous studies [29,49] 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, CircRNA expression was normalized using GAPDH and HMBS as control genes. Statistical significance was performed by Student's t-test (one-tailed). *P < 0.05; **P < 0.01; ***P < 0.001 circFAM107B, circSKAP1, circUTRN, circATP8B4, circPAK2, circABCA13 and circTMEM154), and performed qPCR using specific divergent primers (Supplementary Table S1). The results validated all 15 tested circRNAs, with 10-to 50-fold greater 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 systems such as B cell activation, T cell activation and interleukin signalling (Fig. 4b). The circRNAs were   Table S6). Analysis of the back-splicing frequency showed 167 circRNAs with a 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.

Discussion
CircRNA are highly conserved molecules with distinct and specific expression patterns across tissues and species [50,51]. Their covalently closed ends also endow these molecules with high resistance to degradation such that their half-life, over 48 h, is much longer than comparative linear RNA [52]. These special characteristics make circRNA favourable candidates as biomarkers for human disease, particularly in relation to disorders where access to tissues is limited and liquid biopsy from circulating tissue or fluids such as blood, exosomes and extracellular vesicles [53][54][55][56][57] are preferable. In support of this approach, we observed 22 and 33 circRNAs that were significantly altered in circulating PBMCs from individuals with SZ and BD, respectively, compared to non-psychiatric controls, three of which were also observed to be altered in the combined psychosis group when compared with controls. Most of the differentially expressed circRNAs were downregulated including two (circNCF1 and circLINC00969) that were shared between both disorders. One of these molecules, circDPYD, was recently shown to be dysregulated in plasma samples from coronary artery disease patients and was suggested to regulate TRPM3 through suppression of miR-130a-3p [58]. Interestingly, this circRNA arises from DPYD, a gene located at the genomewide significant risk loci for schizophrenia [59]. Another circRNA, circARHGEF12, was also reported to be significantly downregulated in gastric cancer tissues compared to the paired nontumorous tissue [60]. It is worth noting that we crossreferenced the altered circRNAs to those expressed in the human brain and found the majority of altered circRNAs (48 out of 53) are expressed in the brain. In previous work in psychiatry, Yao et al. observed 9 differentially expressed circRNAs in PBMCs from a small sample of individuals with SZ (n=9) and nonpsychiatric controls (n=9), using microarray [61]. While this finding was supported in a larger validation sample by qPCR, we could not replicate these results in our cohort, which may reflect heterogeneity, which is particularly noticeable in the different sample sizes. Also, different expression profiling platforms may contribute to the observed inconsistency.
The circRNA interaction analysis of differentially expressed circRNAs suggested they modulate a diverse profile of target miRNA in PBMCs. Six of the interacting miRNA including miR-564, miR-572, miR-107, miR-1275, miR-576-5p, miR-342-5p and miR-330-3p were previously reported to be dysregulated in the peripheral blood of patients with SZ and/or BD [62][63][64][65]. Our analysis of gene targets of the interacting miRNA revealed a substantial enrichment of the targets in brain disorders such as SZ and BD, also in GO neuronal clusters including synapse, axon guidance and neurogenesis. These data suggest that the dysregulated circRNAs might play an intermediate regulatory role in pathogenesis of these disorders by modulating the miRNA levels. Nevertheless, due to limitations with in silico analyses, such as prediction errors, more experimental studies are warranted to validate these bioinformatic-based findings in biological systems. Association of circRNA with miRNA was shown to be critical for normal physiology, and the perturbation of this interaction may lead to a pathological condition. For example, dysregulated CDR1as altered miR-7 accessibility and its target genes, resulted in deficits in synaptic neurotransmission and abnormal brain function associated with psychiatric disorders [28]. circHomer1a, a neuronal-enriched circRNA associated with neural plasticity [27], is decreased in the prefrontal cortex (PFC) and stem cell-derived neuronal cultures from patients with SZ and BD [66]. Furthermore, circHomer1a expression levels in dorsolateral prefrontal cortex (DLPFC) and orbitofrontal cortex (OFC) were positively associated with the age of onset of schizophrenia [66]. Functional role of circMTO1 was initially reported in hepatocellular carcinoma (HCC), where suppression of HCC progression was regulated by circMTO1 acting as the sponge of oncogenic miR-9 to promote p21 expression, suggesting its potential as a HCC treatment target [67]. Following studies indicated circMTO1 roles in various cancer pathways through inhibiting its target miRNAs including miR-9, miR-223, miR-630, miR-92 [68,69]. This evidence confirms the functional association of the circRNA transcripts in the pathogenesis of human diseases, mostly through regulating miRNA accessibility.
We also detected RBP sites within the changed circRNAs, with some circRNAs containing more than a hundred binding sites for the interacting proteins. Interaction with RBP may contribute to the transportation, preservation or transcription ability of the target proteins [16,17]. It has been recently shown that circRNAs are capable of producing functional proteins [18]. As such, we assessed the translation features in our altered circRNAs and found six of which have an open reading frame greater than 300 bp with the potential to encode a protein. Interestingly, more than a third of our list of dysregulated molecules were among circRNAs previously shown to be the subject of m 6 A modification, which is thought to be associated with circRNA translation [47].
Our analyses, along with current evidence, suggest that circRNAs could have functional significance in the PBMCs, through which they could associate with neuropathological features. This is due to the immunological link to the pathophysiology of SZ and BD supported by several genome-wide association studies and gene expression reports that showed pathways involved in immunity, apoptosis and inflammation are significantly enriched in these disorders [65]. In addition, maternal infection and birth in the peak infection seasons were suggested to play a critical role in SZ and associated cognitive deficits, by impairing the immune functions [65]. These lines of evidence support immunological manifestation of the disorders, suggesting that PBMCs reflect genomic alterations observed in patients with psychiatric disorders; and thus they can be insightful in investigations of the etiology of these diseases.
Furthermore, we observed some circRNAs have a higher relative back-splicing frequency compared with their linear counterparts, suggesting their potential biological function in regulation of gene expression. Our investigation of potential exon usage within PBMC circRNA structures revealed a diverse array of alternative splicing involved in circRNA biogenesis, with a bias towards the A3SS event. While some AS events were shared between the diagnosis groups, some were exclusively expressed in one group. We also found several AS events were significantly altered in SZ and BD, as compared to HC. These data suggest that alternative splicing in circRNA might be essential for PBMCs physiology, and the splicing pattern of some circRNAs could have important pathological significance.
In summary, we characterized circRNA transcript expression in human PBMCs from SZ and BD patients and observed changes compared to non-psychiatric controls. Further bioinformatic analyses of psychosis-associated circRNA showed they may function as ceRNA for mRNA and form a miRNA interacting platform with the potential to modulate genes enriched in psychiatric conditions. While these preliminary findings are tantalizing, large-scale replication and validation are also warranted to confirm their biological role and potential application as biomarkers with relevance to psychiatric disease and treatment.
Author contribution EM was involved in the research design, methodology and data analysis, and writing the first draft. MG contributed to the study design, data interpretation and revising the manuscript. MC was involved in the supervision of the project, design of the study, interpretation of data and revising the manuscript. All authors read and approved its final version. Availability of data and material As the sequencing data belong to patients, we cannot make them available at this stage.
Code availability Not applicable.

Declarations
Ethics approval This study was approved by the Human Research Ethics committees of the University of New South Wales (HC12384), St Vincent's Hospital (HREC/10/SVH/9) and the South East Sydney and Illawarra Area Health Service (HREC 09/081).

Consent to participate Participants provided signed informed consent.
Consent for publication Not applicable.

Conflict of interest
The authors declare no competing interests.