Bioinformatics Analysis of circRNA Expression and Construction of “circRNA-miRNA-mRNA” Competing Endogenous RNAs Networks in Bipolar Disorder Patients

Bipolar disorder (BD) is a severe mood disorder disease in China, and its underlying pathogenesis remains unknown. Circular RNAs (circRNAs) have been reported to play a key role in mental disorders and can be used as competitive endogenous RNAs (ceRNAs). However, little is known about the correlation of circRNAs with BD. In this study, Deep RNA sequencing was used to identify differentially expressed circRNAs (DE-circRNAs) and differentially expressed mRNAs (DE-mRNAs) between BD patients and a control group. Real-time quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was used to validate the differentially expressed RNAs (DE-RNAs). In all 9,593 circRNAs and 20,030 mRNAs were found in the two groups of specimens, among which 50 DE-circRNAs and 244 DE-mRNAs were significantly upregulated, and 44 DE-circRNAs and 294 DE-mRNAs were significantly downregulated. Based on the regulatory mechanism of ceRNAs, circRNAs can directly bind microRNAs (miRNAs) to affect mRNA expression, and the expression trends of circRNAs and mRNAs are consistent. According to this mechanism, we constructed two ceRNA networks by using the RNA sequencing data. The function of these DE-circRNAs was further elucidated by enrichment analysis. In summary, the present study showed that the circRNA expression profile of BD patients is altered, and a ceRNA regulatory network was constructed, which provided a hypothesis about the pathogenesis of BD.


INTRODUCTION
Bipolar disorder is a mood disorder disease of unknown etiology with a series of serious symptoms. the average age of onset is 30 years old (McIntyre et al., 2020). Evidence from genetic studies indicates that BD shows high heritability, and the associated genetic risk and gene expression profiles overlap (Shao and Vawter, 2008;Cardno and Owen, 2014). In particular, genetic and epigenetic changes can affect gene expression by directly modifying mRNA templates or by regulating post transcriptional translation, which plays an important role in the pathogenesis of BD. Non-coding RNAs (ncRNA), including microRNAs, tRNA, and circRNAs, involve in regulating gene expression and appear to be altered by genetic and environmental exposure in ways that are associated with mental illness, suggesting that they play a critical role in maintaining normal physiological function and homeostasis of the nervous system (Geaghan and Cairns, 2015;Mahmoudi and Cairns, 2017;Rusconi et al., 2020).
Circular RNA is a newly hot spot of endogenous ncRNA that is evolutionarily conserved in eukaryotic species and highly expressed in mammals (Jeck et al., 2013). CircRNAs are formed from special reverse complementary precursor mRNAs, which are cut at variable positions and cyclized head to tail. This special structure makes them more stable and not degraded by RNase R. CircRNAs are increasingly recognized as major epigenetic regulators. Their function is similar to that of a sponge of endogenous RNA or miRNA and can competitively inhibit the transcriptional regulation of miRNA. CircRNAs can also regulate variable splicing or transcription by inhibiting transcription initiation sites (Li et al., 2015), isolate RNA-binding proteins (RBPs) and form RNA-protein complexes affects the localization and transportation of RBPs and their related mRNAs (Fischer and Leung, 2017;Haque and Harries, 2017). Recently, some circRNAs were proven to be successfully encode functional proteins (Legnini et al., 2017). Increasing evidence shows that circRNAs are involved in the occurrence and progression of many human diseases, such as cardiovascular diseases, cancer, and central nervous system disorders, including schizophrenia (SZ), major depressive disorder (MDD), amyotrophic lateral sclerosis (ALS), and Alzheimer's disease (AD) (Su et al., 2019;Jaé and Dimmeler, 2020;Mehta et al., 2020). The study also revealed the perturbation of ncRNAs in tissues and body fluids and suggested that they might serve as diagnostic biomarkers or therapeutic targets (Mehta et al., 2020). For example, 22 circRNAs were dysregulated in the peripheral blood samples of SZ patients. One of those circRNAs, circDPYD, was dysregulated in blood samples from coronary artery disease patients and was proven to increase TRPM3 expression via inhibiting miR-130a-3p (Pan et al., 2017). Interestingly, the parental gene of this circRNA is DPYD, which is located at a genome-wide significant risk locus for SZ. Similarly, a major alteration of circRNAs has been found in ALS patients, among which circ_0023919, circ_0063411, and circ_0088036 are considered candidate diagnostic biomarkers (Dolinar et al., 2019). The Cdr1 as contains more than 70 conserved miRNA target sites, which can sponge miR-7 and thus inhibit the expression of the ubiquitin carboxyl terminal hydrolase L1 (UCHL1) gene; this leads to decreased β-amyloid precursor protein (APP) and β-site APP cleaving enzyme 1 (BACE1) protein levels, indicating that ciRS-7 may be a potential therapeutic target for AD (Zhao et al., 2016;Shi et al., 2017). However, few reports on the circRNA-related ceRNA network in BD patients. Thus, the identification of BD-related circRNA profiles and the characterization of these circRNAs and their functions are urgent goals that will provide very useful information.
This research aimed to investigate differentially expressed circRNAs and mRNAs (DE-circRNAs and DE-mRNAs) in BD patients vs. controls and to comprehensively analyze the "DEcircRNA-miRNA-DEmRNA" ceRNA network involved in BD. The findings of this work may put forward some candidate targets for the development of novel treatment policies for BD.

Samples and RNA Isolation
From January 2019 to December 2019, according to the structured clinical interview of DSM-IV, two psychiatrists diagnosed 20 patients aged 19-27, and collected their blood samples and 20 normal people aged 19-25 as controls (basic clinical data are shown in Table 1). These collected blood samples were applied in two parts. Four blood samples from BD patients and four normal blood samples were used for high-throughput sequencing, and the remaining blood samples were used for subsequent qPCR to verify the sequencing results. Total RNA was isolated from peripheral blood samples using a miRVana TM RNA isolation Kit according to the reagent instructions. RNA concentration and OD260/OD280 were tested with a NanoDrop 2000 instrument (Thermo Fisher Scientific, United States). RNA integrity was tested by agarose gel electrophoresis. This protocol was approved by the Ethics Committee of Jiangxi Mental Hospital, Nanchang, China. All participants gave informed consent for the collection and use of their samples for this study.

Illumina High-Throughput Sequencing
Ribosomal RNA was extracted using the TruSeq Stranded Total RNA with Ribo-Zero Gold following the reagent instructions. RNA libraries were constructed by using rRNA-depleted RNAs with the Illumina Library Prep Kit according to the reagent instructions. The constructed RNA library was qualified by Agilent 2100 Bioanalyzer and sequenced by Illumina sequencer (HiSeq X Ten) in PE150 mode following the manufacturer's instructions. All sequencing raw reads subjected to quality control (Q40). Next, the reads of sequencing difference was Frontiers in Genetics | www.frontiersin.org cleared using Cutadapt software (V1.9.3). The clean reads were used for the analysis of circRNAs and mRNAs by using Hisat2 software (Kim et al., 2015) and the circBase database (Glažar et al., 2014). All data analyses were operated by Shanghai OE Biotech Co., Ltd. (China).

Identification of DE-circRNAs and DE-mRNAs
The number of circRNA counts in each sample was standardized by using DEseq (Anders and Huber, 2012) software, the fold change (log 2 FC) was calculated, and the different significance testing of reads number was carried out by a negative binomial distribution test. HT-seq count software (Anders et al., 2015) was performed to obtain the number of mRNA reads in all samples, and Cufflinks software (Roberts et al., 2011) was used to calculate the fragments per kilobase per million (FPKM) values of the mRNAs. circRNA with p-values lower than 0.05 and log 2 fold change (FC) value higher than 1 were considered differentially expressed, while mRNA with p-values lower than 0.05 and log 2 fold change (FC) value higher than 0.58 were considered differentially expressed.

qRT-PCR Validation of DE-circRNAs and DEmRNAs
Quantitative reverse transcription-polymerase chain reaction was performed using a PCR thermocycler (ABI 7900TH, United States) to validate DE-circRNAs and DE-mRNAs, which included three upregulated and three downregulated circRNAs and mRNAs. The circRNA-specific primer sequences were designed to span the reverse-spliced sequences of circular RNAs but not mRNA of the same sequence (The primer sequences are presented in Table 2). The expression levels of the DE-circRNAs and DE-mRNAs were standardization to (input the reference gene, e.g., β-actin) and were calculated via the 2 − Ct method (1) Ct con = mean CT value of target gene -mean CT value of internal reference gene Ct BD = mean CT value of target gene -mean CT value of internal reference gene, (2) Ct = Ct BD -Ct con , (3) Change Fold = 2 − Ct (con, normal control group; BD, BD group).

Gene Ontology and Pathway Enrichment Analyses
The functional annotation of the target mRNAs of DE-circRNAs was performed by Gene Ontology (GO) 1 term and Kyoto Encyclopedia of Genes and Genomes (KEGG) 2 pathway analyses. Corrected P-values < 0.05 were considered to indicate significant enrichment.

Statistical Analysis
The obtained statistical data were presented as the mean ± SEM by using SPSS 23.0 and were compared via Student's t-test. A twotailed P < 0.05 was considered statistical significance.

Distribution Profiles of circRNAs
A total of 9,593 circRNAs were detected between the normal control and BD groups by RNA sequencing. According to the position of circRNA in their parental gene, it can be divided into different types, as shown in Figure 1A. A total of 4,665 circular RNAs (48.6%) were identified in the circBase database from previous studies, while 4,928 (51.4%) circRNAs were not identified in the database ( Figure 1B). Most of 9,593 circular RNAs in the two groups ranged in length from 100 to 3,500 BP ( Figure 1C). The histogram represents the total cirRNAs distribution of the test sample, and the ordinate represents the number of transcriptional parent genes of the number of cirRNAs on each chromosome ( Figure 1D).

Identification of DE-circRNAs in BD
We used DEG-seq software to analyze the DE-circRNAs between the BD and control groups. A total of 94 DE-circRNAs were obtained according to the significance thresholds of a | log2FC| ≥ 1 and a P-value ≤ 0.05, including 50 upregulated and 44 downregulated circRNAs (Figure 2A). A Manhattan plot shows the distribution of DE-circRNAs on chromosomes, and the threshold P-value ≤ 0.05 is converted it into log value for display, but not FC value. We found that DE-circRNAs were distributed in all chromosomes except chromosome 16 and 21, but not in mitochondria (MT) ( Figure 2B). Furthermore, volcano plot and heat map were presented to identify DE-circRNAs between the two groups ( Figures 2C,D).

Validation of circRNA Expression by Quantitative Real-Time PCR (qRT-PCR)
To validate the expression profiles of those DE-circRNAs, we selected three upregulated and three downregulated circRNAs from the top 10 upregulated DE-circRNAs and the top 10 downregulated DE-circRNAs ( Figure 3A). Figures 3B-G

Identification and Validation of DE-mRNAs in BD
All transcriptional sequencing identified 20,030 mRNAs in the control and BD groups. A total of 538 significantly DE-mRNAs were obtained according to the thresholds of a | log 2 FC| ≥ 0.58 and a P ≤ 0.05, among those 244 were upregulated DE-mRNAs and 294 were downregulated DE-mRNAs ( Figure 5A). We  generated a volcano plot of the DEGs between the two groups ( Figure 5B). In the above results, we predicted 362 miRNAs of upregulated TOP 10 cirRNAs, based on which we further predicted 4,334 mRNAs. We also predicted 772 miRNAs of downregulated TOP 10 cirRNAs, based on which predicted 3317 mRNAs. We intersected the mRNA predicted by the upregulated and downregulated top 10 cirRNAs with the differential genes in our own mRNA sequencing data. Finally, we obtained obtain 52 upregulated and 44 downregulated mRNAs ( Figure 5C). Among these mRNAs, we selected three upregulated and three downregulated mRNAs that we were interested to valid their expression profile. As shown in Figures 5D-I, IL1B, MAFB, and GLUL were the three selected upregulated mRNAs. While ABO, SP6, and TIMP3 were the three selected downregulated mRNAs.

Functional Enrichment Analysis
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway functional enrichment analyses were executed derived from the 52 overlapping upregulated and 44 overlapping downregulated mRNAs in Figure 5C. These overlapping upregulated mRNAs were mainly associated with the following terms: plasma membrane, filopodium, and rough endoplasmic reticulum, in the cellular component (CC) category ( Figure 6A); regulation of cell growth, inflammatory response, and apoptotic mitochondrial changes in the biological process (BP) category ( Figure 6B); and protein heterodimerization activity, transcription cofactor activity, and cytokine activity in the molecular function (MF) category ( Figure 6C). These overlapping downregulated mRNAs were mainly enriched in the following terms: semaphoring receptor complex in the CC category ( Figure 7A); negative regulation of interleukin-8 production in the BP category ( Figure 7B); and superoxidegenerating NADPH oxidase activator activity in the MF category ( Figure 7C). KEGG pathway enrichment analysis showed that these upregulated genes were mainly associated with the NF-kappa B signaling pathway, the Toll-like receptor signaling pathway and antifolate resistance (Figure 6D), while the downregulated genes were mainly associated with glycosphingolipid biosynthesis-lacto, cholinergic synapse and Alzheimer's disease ( Figure 7D).

Visual Graph of ceRNA Networks
Based on the regulatory mechanism of ceRNAs, circRNAs, and mRNAs can competitively bind to miRNA binding sites and act as ceRNAs. Moreover, the expression trends of circRNAs and mRNAs are consistent. Therefore, we established two ceRNA regulatory networks of DE-circRNAs and DE-mRNAs. As shown in Figure 8A, the ceRNA network included the top 10 downregulated DE-circRNAs, 44 downregulated DE-mRNAs ( Figure 5C) and 56 miRNAs predicted by both. The ceRNA regulatory network presented in Figure 8B comprises the top 10 upregulated DE-circRNAs, 52 upregulated DE-mRNAs ( Figure 5C) and 68 miRNAs predicted by both.

DISCUSSION
The incidence of BD, which seriously affects human health and quality of life, is gradually increasing (McIntyre et al., 2020). Although some research progress has been acquired recently, the etiopathogenesis of BD have not yet been complete elucidated. The development of bioinformatics has accelerate the progress of disease mechanisms and treatment strategy . In this work, we comprehensively analyzed circRNA and mRNA expression pattern of BD patients and constructed circRNAassociated ceRNA network.
Circular RNAs are involved in the adjustment of neuronal activity, plasticity, depolarization, and synaptic transmission, which are key steps involved in the pathophysiology of mental diseases (Piwecka et al., 2017;Mahmoudi et al., 2020). Recently, Ebrahim et al., identified 22 and 33 circRNAs showing significant changes in peripheral blood mononuclear cells from SZ and BD patients, respectively, compared to normal controls (Mahmoudi et al., 2021). Although these findings were validated by qPCR in a large sample, we were unable to replicate these results in our cohort, likely because of the different sample sizes. In addition, different sequencing platforms may lead to inconsistent observations. Moreover, these previous authors did not establish the ceRNA regulation mechanism of the DE-circRNAs. Therefore, the expression and function of circRNAs in BD patients are in need of further study. In this research, we implemented an analysis of DE-circRNAs and DE-mRNAs between BD patients and healthy controls and constructed circRNA-associated ceRNA network.
We gained 538 DE-mRNAs and 94 DE-circRNAs between the two groups. Based on the ceRNA theory, circRNAs and mRNAs serve as ceRNAs, and they share the same miRNA binding sites. We constructed two ceRNA regulatory networks. The network of the top 10 upregulated ceRNAs was composed of 52 DE-mRNAs, 10 DE-circRNAs, and 68 target miRNAs and included 134 circRNA-miRNA-mRNA ceRNA circular pathways. The network of the top 10 downregulated ceRNAs was composed of 44 DE-mRNAs, 10 DE-circRNAs, and 56 target miRNAs and included 226 circRNA-miRNA-mRNA ceRNA circular pathways.
To assess the functions of the differentially expressed genes (DEG), GO, and KEGG signaling pathway analyses were executed derived from the 52 overlapping upregulated mRNAs and 44 overlapping downregulated mRNAs. We found that these DE-mRNAs were mainly involved in the regulation of cell growth, immune imbalance, inflammatory response and mitochondrial apoptosis. It has long been reported that BD is a chronic disabling medical disease associated with immune imbalance (Barbosa et al., 2014). During manic or depressive episodes in BD patients, activated T cell counts increase, regulatory T cell counts decrease, and plasma proinflammatory cytokine levels increase. The increased levels of inflammatory cytokines may be involved in the occurrence and development of diseases by regulating neuronal metabolism, synaptic plasticity, and activation of hypothalamic pituitary adrenal (HPA) axis (Raison et al., 2006;Miller et al., 2009). Although the mechanism responsible for the immune disorder is still unclear to a large extent, recent studies have shown that the signaling cascade of innate and adaptive cells is altered (Qian and Cao, 2013). In BD patients, the activity of some signaling pathways involved in immune response increased, such as NF-κB and MAPK signaling pathway (Barbosa et al., 2013). In addition, increases in the NLRP3 inflammasome and caspase-1 levels are associated with increases in serum IL-1B and IL-18 levels in patients with major depression. Toll-like receptor (TLR) family is significant mediators of the inflammatory response that connect innate immunity and acquired immunity. TLR involves a complex intracellular cascade including MyD88 (the common adaptor protein of all TLRs), IRAK, and TRAF family members (Majewska and Szczepanik, 2003;Akira and Takeda, 2004;Ishii et al., 2005). Therefore, immune imbalances and the inflammatory response may be important pathogenic mechanisms of BD. The results of Go and KEGG pathway analysis showed that the DEG in ceRNA network was mainly related to immune imbalance and inflammatory response, which may be the main pathogenesis of BD, providing an important reference for the follow-up study of BD pathogenesis.
Some strengths and limitations of this study should be acknowledged. First, we identified DE-circRNAs and DE-mRNAs in BD patients versus control subjects and comprehensively analyzed a ceRNA network involved in BD. Second, we used bioinformatics technology, which can aid in understanding the genes and networks that are responsible for diseases. Due to the limitation of bioinformatics analysis and the small number of samples, especially in BD group, the statistical ability of microarray analysis is limited. Therefore, more in vivo and in vitro experiments are needed to verify the results of bioinformatics.

CONCLUSION
Our study evaluated circRNA expression in the peripheral blood of BD patients and assessed the potential functions of DE-circRNAs via bioinformatics. Moreover, we created two DE-circRNA related regulatory networks based on the ceRNA theory. In the networks, one circRNA can indirectly regulate multiple mRNAs expression via binding to multiple miRNAs, suggesting that circRNAs may be related to the intricate mechanism of the development of BD.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found here: NCBI SRA PRJNA735880.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Jiangxi Mental Hospital. The patients/participants provided their written informed consent to participate in this study.