miRNA and circRNA expression patterns in mouse brain during toxoplasmosis development

Background Increasing evidence has shown that circular RNAs (circRNAs) are involved in neurodegenerative disorders, but their roles in neurological toxoplasmosis are yet to know. This study examined miRNA and circRNA expressions in mouse brain following oral infection with T. gondii Pru strain. Results Total RNA extracted from acutely infected (11 days post infection (DPI)), chronically infected (35 DPI) and uninfected mouse brain samples were subjected to genome-wide small RNA sequencing. In the acutely infected mice, 9 circRNAs and 20 miRNAs were upregulated, whereas 67 circRNAs and 28 miRNAs were downregulated. In the chronically infected mice, 2 circRNAs and 42 miRNAs were upregulated, whereas 1 circRNA and 29 miRNAs were downregulated. Gene ontology analysis predicted that the host genes that produced the dysregulated circRNAs in the acutely infected brain were primarily involved in response to stimulus and ion binding activities. Furthermore, predictive interaction networks of circRNA-miRNA and miRNA-mRNA were constructed based on genome-wide transcriptome sequencing and computational analyses, which might suggest the putative functions of miRNAs and circRNAs as a large class of post-transcriptional regulators. Conclusions These findings will shed light on circRNA-miRNA interactions during the pathogenesis of toxoplasmosis, and they will lay solid foundation for studying the potential regulation roles of miRNAs and circRNAs in T. gondii induced pathogenesis.


Background
Toxoplasma gondii, the causative agent of toxoplasmosis, is a medically important parasite infecting approximately 30% of the world's population [1]. Immune-potent individuals usually do not show any clinical symptoms, and infection is latent and maintained inside the host as tissue cysts. When infected persons become immunocompromised, such as people suffering from AIDS or organ transplantation, parasite cyst reactivation can lead to acute toxoplasmosis, which includes eye disease, neurological problems and even death [2]. Additionally, infection of pregnant women can cause congenital infection, with dysplasia, hydrocephaly and chorioretinitis occurring in the newborns [3].
The fast proliferating T. gondii tachyzoite is fully controlled by the host potent immune response and then transforms into a slowly replicating stage (i.e. bradyzoite) enclosed in tissue cysts, and remains dormant within the central nervous system (CNS) and muscles [4]. Various kinds of brain cells, including microglia, astrocytes and neurons, can be infected. In the chronically infected mice, the parasite cysts were mainly found in the neurons [5]. In a study of congenital toxoplasmosis, cysts were also found in the neurons [6]. In toxoplasma encephalitis (TE), the rostral basal ganglion is the most commonly affected region, followed by the cerebellum and brain stem [7,8]. In addition, it has been shown that T. gondii infection can affect levels of some neurotransmitters, such as dopamine [9]. Positive correlations between T. gondii infections and neuropsychiatric disorders like schizophrenia [10,11], cryptogenic epilepsy [12] and Parkinson's disease (PD) [13,14] were revealed in the clinical and epidemiological investigations. However, the mechanisms underlying Toxoplasma-induced neuronal disorder in the brain are poorly known.
MicroRNAs (miRNAs) are endogenous small noncoding RNA molecules which regulate gene expression at the post-transcriptional level, and are now widely recognized as essential regulatory molecules involved in neuronal development and function. miRNAs are abundant in the nervous system and a link between miRNAs and the occurrence and development of neurological diseases, such as cerebral ischemia, stroke and neurodegenerative diseases, is becoming increasingly clear [15]. Notably, miRNAs are also found to be regulators of the host response to Toxoplasma infection, as miR-146a and miR-155 are highly induced in chronically infected mice brain, but worthy to mention, compared to miR-146a −/− mice, wild type mice show a more severe TE after challenge with genotype II Pru strain [16]. In addition to using host miRNAs for parasite persistence, T. gondii infection also modulates host cell miRNA profiles in a timely manner [17]. Another type of non-coding RNA, circular RNA (circRNA), is formed by exon-scrambling and was once largely neglected due to its rarity and lack of biological functions. With the development of highthroughput sequencing technology, thousands of new circRNAs have been identified in various organisms, from microorganisms to mammals [18]. circRNA is found to antagonize miRNA activity by a sponge-like mechanism and thus regulates gene expression at the post-transcriptional level [19]. Several pieces of evidence have shown that circRNAs are most highly and specifically expressed in neural tissues, and some circRNAs are expressed in a specific spatial and temporal pattern in the brain among species. More studies have shown that circRNAs may act as important regulatory factors in the pathogenesis of neurological diseases, such as epilepsy, PD and Alzheimer's disease (AD) [20]. However, their roles in the occurrence and development processes of neurological toxoplasmosis remain largely unexplored.
To investigate the potential roles of miRNAs and circRNAs in T. gondii infection, we performed a highthroughput sequencing analysis to identify differentially expressed miRNAs and circRNAs in mice after acute and chronic T. gondii infections, respectively. Furthermore, differentially expressed miRNAs regulated by differentially expressed circRNAs were predicted, and predictive interaction analysis of miRNA-mRNA was also performed. Thus, the present study presented the first transcriptome-wide circRNA landscape and a predictive miRNA-circRNA network in T. gondii-infected mice, which provides a valuable dataset that will help elucidate the mechanisms underlying neurological disorders during T. gondii infection.

Results
Expression of miRNAs in brains from both Toxoplasmainfected and normal mice Infected mice with T. gondii at 11 DPI showed noticeable acute toxoplasmosis signs, including anorexia, hyperpyrexia, and messy hair. All mice survived and restored their physical status at 35 DPI. Initially, we set out to profile the miRNA expression across samples that represent typical infection states. miRNA libraries of brain samples were successfully constructed and summarized in Additional file 1: Table S1. Genome-wide sequencing identified more than 12 million clean reads in each sample. As shown in Fig. 1, the dominant small RNAs (sRNAs) were 20-24 nt in length. More than 93% sRNA sequences were mapped to the reference genome, and the majority of these mapped sRNA sequences were in the forward orientation (Additional file 2: Table S2). Distribution of sRNAs sequences showed that the number of miRNA sequences was higher in the control group (Additional file 3: Table S3). The mapped reads were divided into different categories, and the majority of small RNA reads were belonged to the annotated known miRNAs.
We then compared acute infection samples against normal samples and identified 20 over-expressed and 28 under-expressed miRNAs (adjusted P value < 0.05) (Fig. 2a). Comparison between chronically infected and control samples revealed more dys-regulated miRNAs, including 42 over-expressed and 29 under-expressed miRNAs. As shown in Fig. 2b, 24 miRNAs were shared between AI vs. Con and CI vs. Con. Worthy to mention, all these 24 miRNAs showed the same variation tendency during parasite infection ( Table 1). The clustering of miRNA expression profiles derived from 228 different miRNAs in control and infected samples is shown in Fig. 2c. The tree shows a very good separation between normal and infected samples. However, PCA scores plots did not clearly discriminate chronically infected mice from acutely infected group (Additional file 4: Figure  S1A). A recent study showed that infection with T. gondii oocysts caused the miRNA perturbation in mice brains [21]. The high-throughput data are available in the NCBI database with the accession number PRJNA418218. As shown in Table 2, the differentially expressed miRNAs shared by both cyst and oocyst infections are listed. It is noteworthy that, mmu-miR-155-5p were up-regulated during the whole infection course in both cyst and oocyst infections, which was further confirmed by the qPCR (Additional file 5: Figure S2).

Mapping circRNAs from total RNA-seq data
To investigate whether T. gondii infection affects the expression of the host circRNAs, we profiled the brain circRNAs during the toxoplasmosis progression. To obtain sufficient RNA-Seq data, non-polyAcontaining RNA, linear RNA and rRNA were removed prior to cDNA library preparation. Tested samples were sequenced on an Illumina HiSeq 4000 with an average data volume of 9 Gb, thereby yielding a total of more than 557 million paired-end reads with a size of 150 bp (Additional file 6: Table S4). Clean reads were mapped to the mouse reference genome (GRCm38/p6) by using Bowtie, and about 61% of these reads consisted of protein coding sequences, whereas smaller fractions aligned with ribozyme and long noncoding RNAs (Fig. 3a). circRNAs were identified and filtered by using both the find_circ and CIRI2 software as previously described [22]. A total of 16,543 circRNAs with a minimum of two reads spanning the back-splicing junction were identified. Fig. 1 The average length distribution of the sRNA sequences identified in each group. Sample groups including acutely infected, chronically infected and control are labeled as AI, CI, and Con, respectively Annotations for all circRNAs identified in this study are shown in Additional file 7: Table S5. We classified these circRNAs into three groups, namely, exonic circRNAs, intergenic circRNAs and intronic circRNAs. As shown in Fig. 3b, in all of the 9 sample data sets, exonic circRNAs predominated. The size of these circRNA candidates ranged from less than 50 nt to greater than 1500 nt, most of which ranged from 100 nt to 650 nt (Fig. 3c). As shown in Fig. 3d, most genes generated one or two circRNAs, but some genes yielded multiple circRNAs.

Differentially expressed circRNAs in mice brain during T. gondii infection
We then performed expression profiling to show circRNA variations in mice brain during parasite infection. Expression abundance of circRNAs from nine samples was measured based on TPM, and no abnormal expression was found (data not shown). As a result, compared to control group, 76 circRNAs in acute infection group were detected to be differentially regulated by |log 2 fold change| ≥0.5, adjusted P < 0.05, among which 9 circRNAs were up-regulated while 67 circRNAs were down-regulated. In the chronic infection group, only 3 differentially expressed circRNAs were selected with 2 up-regulated and 1 down-regulated (Fig. 4a). Among these dys-regulated circRNAs, only one was common between AI vs. Con and CI vs. Con, namely novel_circ_ 0057684 (Fig. 4b). Four circRNAs were selected randomly for validations. As shown in Additional file 8: Figure S3, the results are consistent with predictions. Hierarchical clustering showed that circRNA expression pattern in acute infection group was distinguishable compared to the other two groups (Fig. 4c). Meanwhile, PCA scores plots did not clearly differentiate chronically infected mice from control group (Additional file 4: Figure S1B). These data suggested that the expression of circRNAs in acutely infected brain was severely affected.

Functional enrichment analysis of differentially expressed circRNAs
Comparison of differentially regulated circRNAs and their corresponding host genes between the acutely   infected and normal brains are shown in Additional file 9: Table S6. We performed Gene Ontology (GO) analysis on host genes. GO analysis consists of three different aspects, namely biological process, cellular component and molecular function (Fig. 5). Prediction terms with P-value less than 0.05 were selected and ranked by enrichment score (−log 2 (P-value)), and the top 10 generally affected GO terms in each categories are listed. The most enriched biological process terms were related to cellular process, signal transduction and response to stimulus, such as "response to stimulus (GO:0050896)" and "single-organism cellular process (GO:0044763)". The most enriched molecular function terms were mostly about binding activity and kinase activity, such as "ion binding (GO: 0043167)" and "transmembrane receptor protein kinase activity  Analyses of of circRNA-miRNA and miRNA-mRNA interactions circRNAs are a unique class of endogenous non-coding RNAs which act as miRNA decoys or sponges to regulate gene expression. Therefore, an integrated analysis of the expression profile from circRNA-miRNA predictive interactions was performed based on the high-throughput RNA sequencing results. First, using miRanda software, the binding site analysis and predictive interaction analysis of circRNAs-miRNAs identified in the comparison between acutely infected and normal brains were performed. As shown in Fig. 6, 38 dys-regulated miRNAs and 50 dys-regulated circRNAs were involved in this network. The data used to create Fig. 6 are listed in Additional file 10: Table S7. Down-regulated mmu-miR-214-3p is the most frequently targeted miRNA by 11 circRNAs, followed by up-regulated mmu-miR-21a-3p, down-regulated mmu-miR-455-3p and mmu-miR-497a-5p. It is worth noting that the novel_circ_0048152 showed the largest interaction network. All the above five RNAs were selected to validate the expression profiles obtained by RNA-Seq, and the qPCR results were consistent with the RNA-Seq findings (Additional file 11: Figure S4). Next, predictive miRNA-mRNA interacting pairs were ranked by the P value of the hypergeometric distribution. As shown in Table 3, three relationships between miRNA and targeted mRNA (|fold change|≧1.5, P < 0.05) in the acute infection stage and five relationships in chronic infection stage were revealed, respectively. As shown in Additional file 12: Figure S5, the expression levels of these differentially expressed mRNAs were consistent with the high-throughput sequencing data.

Discussion
A previous study indicated that T. gondii infection could modulate the behavior of the intermediate hosts with the presence of tissue cysts in the brain, and the infected rodents showed impaired learning and memory [23]. In humans, increased rates of psychiatric disorders and suicide, and decreased psychomotor performance have been associated with persistent T. gondii infection [24,25]. In order to reveal the mechanisms underlying neuronal disorder, we need to better understand brain gene expression and regulation induced by T. gondii infection. Previous studies have shown that more genes were differentially expressed during chronic infection compared to acute infection after type II strain infection [26,27]. Among these differentially expressed genes, more genes were up-regulated and were primarily involved in immune regulation and cell activation. The fast development of RNA-seq technology offers unique opportunities to discover novel regulating factors in toxoplasmosis pathogenesis. In this study, we comprehensively profiled the miRNA and circRNA expressions of mice brains in order to further our understanding of the disease pathogenesis. Using RNA-seq data, we revealed a global bias for miRNA and circRNA accumulation in T. gondii infected mouse brains. miRNAs, widely spread in animal cells, can bind to the complementary site on the 3′ untranslated region (UTR) of the targeting mRNAs and thus facilitate mRNA degradation or translation inhibition on a post-transcriptional level. In a previous study, 637 miRNAs were identified in Kunming mouse brains at 14 and 21 days post infection with T. gondii cysts, which is a phase of rehabilitation [17].
However, the miRNAs expression has never been profiled in the acute and chronic infection stages after T. gondii infection. In this study, we identified 1314 known miRNAs and 91 novel miRNAs, of which 48 and 71 differentially expressed miRNAs were found in acute and chronic infection, respectively. Mmu-miR-155-5p is the most up-regulated miRNA during the whole infection course, which also showed significantly perturbed expression in oocyst-induced toxoplasmosis [21]. mmu-miR-155-5p showed an elevated level in chronic infection compared to acute infection, which might be resulted from the increased high chronic cyst burdens. Several studies provide evidences that miR-155 (the precursor of miR-155-5p) is involved in the innate and acquired immune responses, hematopoiesis and autoimmune disorders [28][29][30]. More evidences show that elevated levels of miR-155 occur in the formation and development of several tumors, such as leukemias and breast cancer. Both the levels of mmu-miR-142a-3p and mmu-miR-142a-5p were up-regulated in this study. miR-142a-3p has been identified as an essential player in the formation and differentiation of hematopoietic stem cells, in which miR-142a-3p directly targets transcription factor p53 which is involved in the biological process of hematopoietic stem and progenitor cells [31]. Transplantation of mesenchymal stem cells (MSCs) shows good therapeutic effects on acute lung injury, during which Beclin-1 protein is up-regulated. miR-142a-5p is found to negatively regulate Beclin-1 production, and its expression is inhibited during the transplantation of MSCs [32].
Over-expression of the miR-223 significantly reduced Plasmodium falciparum infection during the intraerythrocytic life cycle [33]. In this study, both mmu-miR-223-3p and mmu-miR-223-5p were up-regulated, which might be involved in host defense against T. gondii infection. Mmu-miR-203-3p was also constantly up-regulated during the whole infection course. Interestingly, the miR-203 functions as a tumor suppressor and is down-regulated in many kinds of cancers, such as primary prostatic tumors, cervical cancer and rhabdomyosarcoma [34][35][36]. In the present study, mmu-miR-185-3p was the most down-regulated miRNA. More evidences showed that the expression level of miR-185 decreased in carcinogenesis, in which miR-185 suppresses tumor proliferation by directly targeting DNA methyltransferase 1 (DNMT1) [37,38]. circRNAs are produced by RNA back splicing and highly prevalent in the eukaryotic transcriptome [39]. While the functions of circRNAs remain poorly known, some circRNAs regulate gene expression by acting as miRNA sponges. For instance, a heart-related circRNA (HRCR) was found to function as an endogenous miR-223 sponge, which might be involved in the protection of the heart from pathological hypertrophy and heart failure [40]. In the present study, a large number of circRNAs were identified in both infected and normal mouse brains. Most circRNAs have low abundance, and certain circRNAs are expressed in one gene locus. Interestingly, 76 and 3 differentially expressed circRNAs were identified in the acute and chronic infection stages, respectively, which suggests that these circRNAs may play key roles in toxoplasmosis pathogenesis. In order to explore the roles of these differential circRNAs in mouse toxoplasmosis, GO analysis was used to annotate the biological functions of host linear transcripts in the acute infection stage, and the results showed that a significant amount of GO terms were related with the cellular process and binding activity. circRNA-miRNA predictive interaction networks were also established to predict the relationships between circRNAs and miR-NAs. Generally, a circRNA harbors more binding sites for different miRNAs, which means that the circRNA may play multiple functional roles in the regulation of miRNA target gene expression. In this study, it was found that the novel_circ_0048152 showed the largest interaction network and might potentially interact with seven miRNAs including mmu-miR-206-3p, mmu-miR-1a-3p, mmu-miR-598-3p, mmu-miR-147-3p, mmu-miR-484, mmu-miR-34a-5p and mmu-miR-21a-3p. Further biological roles of the novel_circ_0048152 still need to be explored. Moreover, we found that mmu-miR-214-3p, mmu-miR-21a-3p, mmu-miR-455-3p and mmu-miR-497a-5p were potentially co-expressed with multiple circRNAs. However, their functional roles during toxoplasmosis have never been reported.

Conclusions
The present study revealed unique sets of miRNAs and circRNAs along with their expression profiles. Their potential roles were predicted by bioinformatics softwares. Predictive interaction networks were constructed for both miRNA-mRNA and circRNA-miRNA. Our data lay a solid foundation for studying the functions of miRNAs and circRNAs in toxoplasmosis pathogenesis, which will be performed in further experimental studies.

Mice and parasite
Six to eight-week-old female BALB/c mice were purchased from the Experimental Animal Center of Lanzhou Veterinary Research Institute, Chinese Academy of Agriculture Sciences, PR China. Strict animal handling guidelines were followed in this study. All animals had free access to food and water. They were kept at 22°C in a 12/12 h light/dark cycle.T. gondii type II Pru strain was used in this study, and was maintained throughout the study via oral inoculation of cysts in mice. Cysts were obtained from brain tissues of infected Kunming mice after 60 days post infection (DPI).

Mouse infection and sample collection
Nine BALB/c mice were randomly assigned to three experimental groups and each group contains three individuals. Six mice were infected orally with~10 freshly prepared T. gondii cysts by gavage. Meanwhile, three mock-infected (control) mice received 100 μL saline alone. All infected mice were checked daily for the development of clinical signs of toxoplasmosis. At 11 and 35 DPI with T. gondii cysts, infected (n = 3 at each time point) and control mice (n = 3) were anesthetized with 5% isoflurane gas at around 10 a.m., and then sacrificed by cervical dislocation. The brain tissues were collected and immediately snap-frozen in liquid nitrogen and kept at − 80°C until use.

RNA extraction
Total RNA was extracted from brain tissues using Trizol method (Invitrogen, Carlsbad, CA, USA). The purity of RNA was determined using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). Integrity of RNA was confirmed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The concentration of RNA was quantified using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA samples were stored at − 80°C until used.

RNA sequencing
The detailed method for mRNA sequencing was described previously [41]. Three micrograms of the total RNA per sample was used for the construction of the small RNA library. Sequencing libraries were constructed using NEB-Next® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) according to the manufacturer's recommendations. The quality of sequencing library was evaluated on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. Small RNA sequencing libraries were constructed with HiSeq 2500 sequencing platform and 50 bp single-end reads were generated. Raw data of fastq format were firstly processed through in-house perlscripts. Q20, Q30, and GC-content of the clean data were calculated. In order to analyze the small RNA expression and distribution on the genome, small RNA tags were mapped to reference sequence by Bowtie [42]. MiRBase 20.0 was used to align small RNA to the miRNA precursor to obtain the known miRNA. By exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated small RNA tags, softwares miREvo and mirdeep2 were integrated to predict novel miRNAs [43,44].
For the construction of circRNA library, 5 μg RNA per sample was used as input material. Ribosomal RNA was removed by Epicentre Ribozero™ rRNA Removal Kit (Epicentre, USA) followed by the linear RNA digestion with RNase R (Epicentre, USA). The sequencing libraries were constructed by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's instructions. According to the manufacturer's instructions, circRNA sequencing libraries were constructed with HiSeq 4000 sequencing platform and 150 bp paired-end reads were generated. Reference mouse genome was downloaded from the NCBI genome website directly. Index of the reference genome was built using bowtie2 v2.2.8 and paired-end clean reads were aligned to the reference genome using Bowtie. The mapped reads were removed. The unmapped reads, including splice junctions-aligned reads were further processed by find_circ and CIRI2 [45,46]. At last, the circRNAs that contain splice junctions-aligned reads were identified. The intersection between the two algorithms for circRNA prediction was selected.

RNA quantification and differential expression analysis
The FPKM (fragments per kilo-base of exon per million fragments) of the coding genes was calculated using Cuffdiff (v.2.1.1) to evaluate the mRNA expression level in each sample [47]. The expression levels of miRNAs and circRNAs were estimated by TPM (transcripts per million reads) through the following criteria: Normalized expression = actual miRNA readcount or circRNA junction readcount× 1,000,000/libsize (libsize is the sum of miRNA or circRNA readcount) [48].
Differential expression analysis of two groups was conducted using the DESeq 2(1.14.1). The P-values were adjusted using the Benjamini & Hochberg method. Corrected P-value of 0.05 was set as the threshold for significantly differential expression by default [49]. By using metaboanalyst 4.0 online software, these differentially expressed RNAs were selected to perform hierarchical cluster analysis and principal component analysis (PCA) [50].

Quantitative PCR
Total RNA was extracted using Trizol reagent (Life Technologies). To quantify the amount of mature miRNA, we used the miRcute Plus miRNA qPCR Detection Kit (SYBR Green) (FP411, Tiangen Biotech, Beijing, China) and small nuclear U6 as an internal standard. The reactions were performed at 95°C for 15 min, 94°C for 20 s, 60°C for 30 s, 72°C for 34 s. All reactions were run in triplicate. To quantify the amount of circRNAs, cDNA was synthesized with the HiScript II Q RT Super-Mix (Vazyme, Nanjing, China). Quantitative real-time PCR (qRT-PCR) was performed with a ChamQ SYBR qPCR Master Mix kit (Vazyme, Nanjing, China) and GAPDH was used as a normalization control. The reactive cycle consisted of 95°C for 30 s, then 40 cycles of 95°C for 5 s and 60°C for 34 s. A set of divergent primers (the forward primer being located downstream of the reverse primer) were designed and used for backspliced circRNA detection. All PCR primers were designed using oligo7 software and are listed in Additional file 13: Table S8. The qRT-PCR was run on a Bio-Rad CFX connect Real-Time PCR Detection System (BIORAD, USA) and the qRT-PCR data were analyzed using the comparative 2 −ΔΔCT method.

miRNA target gene prediction and functional analysis
The miRanda, PITA and RNAhybrid softwares were used to predict potential target genes of all the differentially expressed miRNAs. An intersection among the three combined predictions were selected. Threshold parameters for PITA were set as follows: max target length: 50,000; energy cutoff: − 10; p-value cutoff: 0.05. Parameters for RNAhybrid were set as follows: utr: 3utr.fa; mir: mature.fa. Functional annotation of the predicted microRNA targets were conducted using the GO database. GOseq package was downloaded (http:// bioinf.wehi.edu.au/software/goseq/) and used in this study [51]. GOseq based on Wallenius non-central hyper-geometric distribution was implemented for GO enrichment analysis.
Prediction of circRNA-miRNA interactions circRNAs function as miRNA sponges and affect miRNA-mediated regulation of target gene expression via natural sequestration or competitive suppression of miRNA activity [52]. circRNA-miRNA interactions were explored with these obtained differentially expressed miRNAs and circRNAs. By using miRanda software, circRNA-miRNA predictive interaction network was built based on the prediction and analysis of miRNA binding sites and the correlations between circRNA and miRNA that were ranked as previously described [53]. In this study, we set the Smith-Waterman hybridization alignments match score higher than 140 and the minimum free energy of the duplex structure less than − 10 to guarantee the reliability of our results. Other setting parameters were as follows: gap open penalty:-9; gap extend penalty: − 4; scaling parameter: 4. Cytoscape 3.71 software was used to diagram the map of the circRNA-miRNA interaction network [54].

Statistical analysis
Statistical analyses were performed using SPSS 20.0 software. All data were expressed as the mean ± SEM. P value less than 0.05 was considered statistically significant. Student's t test was used to compare the qPCR results.