Large-scale quantitative genomics analyzes the circRNA expression profile and identifies the key circRNA in regulating cell proliferation during the proliferation phase of rat LR

Abstract Researchers have been exploring the genetic mechanisms underlying the control of liver regeneration (LR). However, an integrated analysis of circRNAs expression of rat regenerating livers during the proliferation phase has not been performed yet. For this purpose, circRNAs expression profile was globally analyzed by high-throughput sequencing. It showed that 10,003 circRNAs were detected, and 164 circRNAs were differentially expressed. Subsequently, 27 circRNAs were predicted to bind to 58 candidate miRNAs and compete for miRNA-binding sites with 2195 mRNAs. By applying GO and KEGG analysis, it was predicted that these circRNAs significantly participated in tissue regeneration, regulation of cell proliferation and Ras, p53, Wnt, Jak-STAT, MAPK signalling pathways. Based on the number of the corresponding miRNAs and their role enriched and reported in cell proliferation of LR or hepatocellular carcinoma, four kinds of circRNAs (circ_03848, circ_08236, circ_13398 and circ_15013) were considered as the key circRNAs. The predicted competing endogenous RNA networks and bioinformatics analysis revealed the potential role of these circRNAs in LR, which would provide useful information for understanding the mechanism of LR.


Introduction
The liver is a vital organ for detoxification, protein synthesis, digestion and other functions, and its regenerative capacity is powerful in animals. 2/3 partial hepatectomy (PH) in rats is an ideal model for liver regeneration (LR) [1]. The process of LR is generally divided into three phases: priming (0-6 h after PH), proliferation (12-72 h after PH) and termination (72-168 h after PH) [2]. Hepatic cells are activated from G 0 phase to G 1 phase in the priming phase [3]. Especially, activated hepatic cells enter S phase and perform two cell cycles in the proliferation phase. A number of studies were conducted to reveal the genetic mechanism regulating the second phase of LR, for instance, hepatocyte growth factor (HGF), mitosis-promoting factors such as epidermal growth factor (EGF) and transforming growth factor alpha play an important role in this phase [4]. However, the mechanism is not yet completely clear.
circRNAs were discovered in 1976, and are predominantly generated by back-splicing reactions, in which the 5 0 and 3 0 ends of linear RNAs are directly ligated [5], and are widely found in many kinds of tissues and cells, and own the characteristics of stable structure, conservative sequence and specificity of expression pattern. According to composition and different sources, circRNAs can be divided into the following types: exonic circRNA, intronic circRNA (circular intronic RNAs), intergenic circRNA, antisense circRNA and sense-overlapping circRNA [6].
It has been verified that circRNAs work through various ways including acting as miRNAs sponge, regulating transcription by interacting with small nuclear ribonucleoproteins (snRNP), recruiting RNA polymeraseII, interacting with proteins, and coding protein [7]. We had shown that 159 circRNAs were differentially expressed (DE) in the priming phase after PH [8]. In the present study, we aimed primarily to investigate the potential mechanism of circRNAs regulating the proliferation phase after 2/3 PH in rats.
Zhengzhou University (Zhengzhou, China). Animals were housed in standard cages for one week at least before the experiments. The animals were kept at a constant temperature of 22 ± 1 C, relative humidity of 55-65% and a 12 h light-dark cycle. 2/3 PH was performed on adult male rats (10 to 12 weeks old) as previously described by Higgins [9] using ether anaesthesia. At 0, 12,24,30,36 and 72 h after PH, the right liver lobes from six rats at each time point after PH were pooled to account for biological variation between 9:00 and 11:00 am, and stored at À80 C after liquid nitrogen freezing. Excessive inhalation of CO 2 was used as the euthanasia methods after the experiment. All experimental procedures in this study were approved by and performed according to the guidelines for the care and use of experimental animals that have been established by the Animal Care and Use Committee of the Ministry of Science of the People's Republic of China.

High-throughput RNA sequencing
All harvested liver tissues were sent to OE Biotech (Shanghai, China) for the whole transcriptome sequencing and quantitative analysis. Briefly, total RNA was extracted from each sample at 0, 12, 24, 30, 36 and 72 h after PH by using TRIzol reagent (Invitrogen, USA) following the manufacturer's protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA, USA). The samples with RNA Integrity Number ! 7 were subjected to the subsequent analysis. The libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold (Illumina, San Diego, CA, USA) and Ribonuclease R according to the manufacturer's instructions. Then these libraries were sequenced on the Illumina sequencing platform (HiSeqTM 2500 or other platform) and 150/125 bp paired-end reads were generated. The reads were mapped, quantified and annotated as described by Li et al. [8]. Briefly, the clean reads were aligned to the reference genome by Bowtie2. And for unmapped reads, the junctions were picked out using backsplice algorithm. circRNAs were verified with a software developed by OE which were considered as the reference sequence for further analysis. Finally, the expression level of circRNAs was measured by RPM (mapped backsplicing junction reads per million mapped reads, RPM ¼ number of circular reads/number of totalreads (units in million)).

Identification of DE circRNAs
To assess the changes of circRNA expression during the proliferation phase, DE circRNAs were defined by two statistical criteria (jlog2Foldchangej ! 1 and p values .05). It was performed by the negative binomial distribution test based on the DESeq package in the same circRNA in each group (12,24,30,36 and 72 h after PH) compared to normal liver tissue (0 h after PH), which was used to control filtration upon the statistics of alignment quality scores (the number of junction reads counts of circRNA in each sample was standardized). The hierarchical heatmap was performed to show the circRNA expression pattern among samples by Cluster 3.0.

Target miRNAs prediction of circRNAs
More and more researcher reported that circRNAs act as miRNA sponges by directly targeting miRNA to influence related miRNA activities. The candidate target miRNAs of circRNAs were predicted by miRanda [10]. Interaction between miRNAs and circRNAs were defined with a maximum binding Max energy of -30 and a minimum Max score of 150. Free energy refers to the energy released when circRNA formed a stable structure with miRNA, the lower value of free energy indicated the stronger binding ability. Hypergeometric distribution assay was used to identify miRNAs that had a significant effect on circRNA during the proliferation phase of LR. While smaller p-values of enrichment indicated DE circRNAs was enriched in the corresponding miRNA.

Downstream mRNA prediction of miRNAs
As a small endogenous non-coding RNAs, miRNAs mediate post-transcriptional gene silencing by guiding Argonaute (AGO) proteins to RNA targets [11]. For high-throughput experimental studies, it is a popular way that miRNA function was defined by its predicted downstream mRNA using TargetScan (http://www.targetscan.org) according to the threshold of total context þþ score À0.6 [12].

GO and KEGG enrichment analysis of the downstream mRNAs
To reveal the potential role of DE circRNAs, their downstream mRNAs of the candidate target miRNAs were submitted to the DAVID website (http://david.abcc.ncifcrf.gov) for the GO and KEGG enrichment analysis [13]. The enriched GO terms were presented by enrichment score. KEGG enrichment analysis was used to determine the involvement of downstream mRNAs in different biological pathways. The -log10 (p-values) yielded an enrichment score indicating the significance of GO term and KEGG pathway correlations.

Construction of circRNA-miRNA-mRNA networks
In recent years, there is an increasing amount of evidence that circRNAs act as miRNA sponges [14], and miRNAs cause post-transcriptional gene silencing by base pairing with mRNA for degradation and/or translational repression [15]. The miRNA expression profiles in liver injury were reported [16,17]. It has been demonstrated that miRNAs played an important role in LR [18]. One circRNA can sponge a number of target miRNA, and a target miRNA can be regulated by multiple circRNAs. Similarly, a miRNA can silence a number of target genes, and a target gene can be regulated by multiple miRNAs. Thus, to well elucidate the possible function of DE circRNAs during the proliferation phase, the network of circRNA-miRNA-mRNA is necessary. Complying with the above criterion, circRNA-miRNA-mRNA networks were displayed by Cytoscape 3.2.

Screening of key circRNAs
A combination of deep sequencing and bioinformatics analysis is a more effective way to identify key circRNAs and their potential function in the rat liver. The selection of key circRNAs among DE circRNAs was carried out according to the following principles. The number of nodes in the circRNA-miRNA network was not less than 3. Moreover, give priority to circRNAs that target the miRNAs associated with cell proliferation in rat LR or hepatocellular carcinoma (HCC) that were confirmed by experiments.

Validation of candidate key circRNA by qRT-PCR
To validate the data from high-throughput sequencing, the candidate key circRNAs were amplified with specific divergent primer in the PH groups and the control groups. Each pair divergent primers were designed according to the corresponding linear transcript sequence and synthesized in Sango Biotech (Shanghai, China). Specific primers designed for circRNAs are listed in Table 1. Total RNAs from regenerating liver and normal tissues were extracted using the TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. RNA quality was estimated by spectrophotometry and agarose gel electrophoresis. Then the total RNAs was incubated 1 h at 37 C with 20 units of RNase R (Epicentre Biotechnologies). RNA was subsequently purified by phenolchloroform extraction. 1 lg of RNA per sample was reverse transcribed into cDNA using a cDNA Reverse Transcription Kit (Takara, Tokyo, Japan). The PCR program was as follows: 95 C for 1 min, followed by 40 cycles of 95 C for 30 s, and 60 C for 30 s and 72 C for 45 s. Expression of circRNAs was normalized to GAPDH (internal standard control) and calculated using the 2 ÀDDCt method. All qRT-PCRs were performed in triplicate.

Statistical analysis
Statistical analysis was performed with software SPSS 17.0. The data were presented as mean ± SEM. Group comparisons were evaluated by one-way ANOVA to determine statistical significance. Differences were statistically significant when p values .05.

Genomewide identification of circRNAs during the proliferation phase of rat LR
To investigate circRNAs at a genome-wide level, the total RNAs was isolated from six rat livers. After eliminating ribosome RNAs (rRNAs) and treating with RNase R, it was utilized to construct libraries for deep sequencing. The sequencing data reached 18 circRNAs expression pattern during the proliferation phase of rat LR The whole transcriptome sequencing was used to detect circRNA expression profile using high-throughput sequencing and the Circbase database [8]. 14.41% of circRNAs had the predicted spliced length of more than 2000 nt. circRNAs with the length of less than 500 nt were 52.14%, and the length from 501 nt to 2000 nt were 33. 45% (Figure 1(B)). Expression of circRNAs in all samples was measured based on RPM (mapped back-splicing junction reads per million mapped reads), which indicated that there was no abnormal expression in the six groups ( Figure 1(C)). circRNAs were classified into five categories according to their relation with protein-coding genes: 1.30% circRNAs were antisense, 3.92% were exonic, 6.32% were intergenic, 0. 78% were intronic, and 87.68% were sense overlapping ( Figure 1(D)). They were widely distributed in all chromosomes (Figure 1(E)).

Identification of DE circRNAs during the proliferation phase of rat LR
In this study, the following criteria: jlog2Foldchangej !1 and p-values .05 were used to identify the DE circRNAs, which performed by the negative binomial distribution test. As a result, 164 DE circRNAs were found between the normal (0 h) and regenerating livers (Additional file S3), which suggested that multiple circRNAs were involved in LR. Venn diagram displayed the DE circRNAs with statistical significance between regenerating and normal livers (Figure 2(A)). Hierarchical clustering exhibited the circRNA expression pattern in regenerating livers (Figure 2(B)). 126 out of 164 DE circRNAs were upregulated while 38 were downregulated. Thereinto, the number of upregulated circRNAs at 12, 24, 30, 36 and 72 h after PH was 55, 28, 38, 51 and 7. And the number of downregulated circRNAs in the same order was 6, 12, 8, 2 and 23, respectively (Figure 2(C)). According to different composition and sources, the upregulated circRNAs were divided into 1 intergenic, 112 sense overlapping, 1 intronic, and 11 exonic. On the other hand, the downregulated circRNAs were divided into 34 sense overlapping and 4 intronic (Figure 2(D)).
Bioinformatics analysis of competing endogenous RNA (ceRNA) networks during the proliferation phase of rat LR circRNAs can affect gene expression at post-transcriptional level by acting as miRNAs sponge [14]. According to the threshold of max energy À30 and max score ! 150, 35 out of 164 DE circRNAs were predicted to target 92 miRNAs. Furthermore, a total of 58 out of above 92 miRNAs were predicted to target 2195 downstream genes according to the threshold of total context þþ score À0.6 using TargetScan (Additional file S4).
According to the threshold of total context þþ score ( À0.6), 27 out of 35 DE circRNAs might indirectly regulate 2195 downstream mRNAs by their corresponding target miRNAs. To elucidate the potential function of 27 circRNAs during the proliferation phase of LR, we furtherly combined 2195 predicted downstream mRNAs of 58 candidate miRNAs. Finally, we established and displayed the circRNA-miRNA-mRNA networks using Cytoscape 3.2 (Additional file S5, Figure 3).

Enrichment of DE circRNAs during the proliferation phase of rat LR
To explore the potential functions of above 27 DE circRNAs, GO and KEGG enrichment analyses were performed to classify the downstream mRNAs. The downstream mRNAs mainly encoded cellular components of the plasma membrane, organelle membrane, mitochondrial membrane, etc. (Figure  4(A)) and associated with transcription regulator activity, GTP binding, DNA binding, etc. molecular functions (Figure 4(B)). Simultaneously, they were enriched in regulation of cell proliferation, small GTPase mediated signal transduction, intracellular signalling cascade, etc. biological process ( Figure  4(C)). Moreover, KEGG showed that they were significantly (p .05) enriched in p53, MAPK, Wnt and Jak-STAT signalling pathway (Figure 4(D)). It predicted that most circRNAs might play a vital role in acting as a "sponge" or activating factor in organelle membrane to regulate cell proliferation through p53, MAPK, Wnt and Jak-STAT signalling pathway during the proliferation phase of rat LR.

Screening of candidate key circRNAs during the proliferation phase of rat LR
In this study, the criteria for screening the candidate key circRNAs were based on the number of the corresponding miRNA and their role enriched and reported in cell proliferation of LR or HCC ( Figure 5(A)). Additionally, the expression profiles of miRNA and mRNA in the regenerating livers of rats during the proliferation phase of LR were reported in a published paper of our group (DOI:10.1002/jcp.28529). Combined with previous studies and as shown in Figure 5(B), 4 out of above 27 DE circRNAs, including circ_03848, circ_08236, circ_13398 and circ_15013 were considered as the candidate key circRNAs. For example, circ_03848 was predicted to bind to 4 miRNAs (rno-miR-207, rno-miR-30c-1-3p, rno-miR-532-3p and miR-150-5p) which were reported involved in LR. While the downstream mRNAs of these miRNAs were associated with transcription regulator activity, protein transport, Ras signal transduction and so on. Furthermore, the network of the candidate key circRNAs-miRNAs-mRNAs, which with a correlated expression pattern according Figure 5(B) and associated with cell proliferation was constructed by Cytoscape 3.2 ( Figure 5(C)).

Validation of the candidate key circRNA
To test whether above candidate key circRNAs were indeed circular, PCR was performed using divergent specific primers with cDNA or genomic DNA (gDNA) firstly. It showed that they could be amplified using cDNA but not gDNA ( Figure 6(A)). The PCR products of these amplified circRNAs were identified directly by sequencing. The analysis of sequencing data confirmed that head-to-tail splicing existed in these circRNAs. To further verify the circular characteristics and the expression changes of these circRNAs during the proliferation phase of rat LR, they were also amplified and quantified by qRT-PCR after RNase R digestion. The expression tendencies of them were consistent with the highthroughput sequencing data ( Figure 6(B-E)).

Discussion
circRNAs were a class of molecules with important biological functions, which were first reported more than 40 years ago. Recently, a large number of circRNAs have been successfully identified in different species such as archaea [19], drosophila [20], rat [8], sheep [21], pig [22] and human [23]. Among them, some circRNAs were found to play an important role in cell proliferation. For example, Du et al. reported that circ-Foxo3 retarded cell cycle by forming ternary complexes with p21 and CDK2 [24]. It was known that the recovery of liver mass mainly resulted from hepatocyte proliferation after 2/ 3 PH in rats [25]. Our team studied the expression changes of circRNAs during the priming phase (0-6 h after PH) of rat LR and found that 159 circRNAs changed significantly [8]. After this phase, hepatocytes entered the cell cycle in the proliferation phase (12-72 h after PH). But how circRNAs work at this phase is not clear. The present study found that 164 dysregulated circRNAs could be involved in this critical step. Among them, 164 DE circRNAs originated from 150 host linear transcripts. Six circRNAs were intergenic indicating that they did not have the corresponding linear transcripts. And seven genes could produce two or three circRNAs during the proliferation phase of rat LR. To explore the role of circRNAs in the proliferation phase of rat LR, the ceRNAs interaction network of DE circRNAs were drawn, and GO enrichment analysis of the downstream target mRNAs were performed. As a result, the biological processes were closed with tissue regeneration, regulation of cell cycle and regulation of cell proliferation, etc. It was consistent with the previous research results that quiescent hepatocytes undergo one or two rounds of replication to restore the liver mass after PH [2,26]. Meanwhile, KEGG analysis showed that downstream mRNAs of DE circRNAs were connected with Ras, p53, Wnt, Jak-STAT and MAPK signalling pathways, which had been acknowledged as the important signalling pathways in the proliferation phase of rat LR. For example, Kunimoto et al. reported that the novel IQGAP3/ Ras/ERK-related signalling was involved in hepatocyte proliferation of LR [27], indicating the critical role of Ras signalling pathway in LR. Moreover, Stepniak et al. reported that c-Jun/ AP-1 promoted G 1 /S transition in mice hepatocytes through a p53-dependent pathway after PH [28]. In addition, Wnt/ b-catenin, JAK-STAT and MAPK signalling pathways have been proved to be a major regulator of regeneration [29][30][31]. Overview, we hypothesized that DE circRNAs might regulate the expression of downstream mRNAs through their target miRNAs, thereby controlling cell proliferation through above signalling pathways in the proliferation phase of rats LR.
After having validated their expression pattern and compared the expression trends between high-throughput sequencing and qRT-PCR, a combination of target miRNA function of circRNAs must be considered, we gave the highest priority to the candidate key circRNAs, including circ_03848, circ_08236, circ_13398 and circ_15013. Because their candidate target miRNAs were involved in LR or HCC. For instance, Chaveles et al. described that miR-207 and miR-532 were significantly downregulated in rat liver after PH during the proliferation phase [32,33], indicating they might play a crucial role in LR. With a negative correlated expression pattern, we found that circ_03848 was upregulated sharply, and it was predicated to inhibit the activity of rno-miR-207 and rno-miR-532 to regulate cell proliferation.
It was known that circRNA played a role in post-transcriptional regulation by several ways such as miRNA sponges or protein templates [14]. After predicting circRNA-miRNA pairs, the predicted target miRNAs of circ_03848, including rno-miR-150-5p and rno-miR-378a-5p were reported to have a correlated expression pattern with it and associated with LR [34,35], suggesting that they could be directly bound to circ_03848 to affect rat LR. In detail, it was found that circ_03848 exhibited an upregulated pattern at 30 h after PH in this paper, and it might target rno-miR-150-5p and rno-miR-378a-5p. Coincidentally, Yu et al. reported that the expression of miR-150 exhibited the sharpest downregulation from 12 to 48 h after PH during LR [18]. In addition, miR-378 expression declined, accompanied by remarkably increased expression of its target ornithine decarboxylase 1 (ODC1) from 0 to 18 h after PH [36], and ODC1 was found that involved in positive regulation of cell population proliferation [37]. Our results were consistent with those reports by others about the expression of ceRNAs was negatively correlated [38], showing that circ_03848 might target these miRNAs and affect downstream gene expression in regulating cell proliferation. Additionally, circ_08236 exhibited an upregulated pattern at 36 h after PH and predicted to target rno-miR-143-5p in this paper. With a negative correlated expression pattern, rno-miR-143-5p was reported to exhibit the sharpest downregulation from 12 to 48 h after PH during LR [39]. Moreover, both of circ_13398 and circ_15013 were predicted to bind to rno-miR-484, showing this miRNA could be regulated by them. Interestingly, miR-484 was highly expressed in tumor liver compared to normal liver [40]. And Yang et al. demonstrated that miR-DEN-induced precancerous lesions and HCC were dramatically impaired in miR-484 À/À mice [41]. Coincidentally, we also found that rno-miR-484 was upregulated in the proliferation phase of rat LR, suggesting that it might play an important role. Therefore, these candidate key circRNAs were very likely to contribute to cell proliferation by sponging their target miRNAs during the proliferation phase of rat LR.
The ceRNAs network and bioinformatics analysis preliminarily revealed the role of circRNAs in regulating cell proliferation of LR in this paper. For example, during the proliferation phase of rat LR, circ_03848 might sponge rno-miR-378a-5p and/or rno-miR-532-3p and promote the expression of their target genes ODC1 and/or RAD51, in turn, accelerate cell proliferation. In addition, circ_03848 might also be combined with rno-miR-324-3p and promote the expression of APEX1 to regulating the G 1 /S transformation and thus promote the process of cell proliferation. It was consistent with the results of He et al. [33,36,42]. Additionally, Assy et al. found that the vascular endothelial growth factor A (VEGFa) played an important role in rat LR which induced by PH [43]. Schug and Yu et al. found that rno-miR-143-5p and miR-150 were significantly downregulated in regenerating liver at 12-48 h after PH, and both of them regulated LR through targeting VEGFa [18,39]. It was consistent with the results in this study, and presumably circ_03848 and/or circ_08236 might regulate rat LR process through rno-miR-143-5p/VEGFa and/or rno-miR-150/VEGFa axis. In summary, we made scientific speculation that circ_03848, circ_08236, circ_13398 and circ_15013 regulated hepatocytes proliferation through the ceRNAs network of circRNA-miRNA-mRNA axis in the proliferation phase, thereby regulating the process of LR.
In conclusion, 10,003 circRNAs were detected and 164 of them were preliminarily determined to change remarkably in the proliferation phase of rat LR using high-throughput RNA sequencing. The ceRNAs network and bioinformatics analysis preliminarily revealed the potential role of DE circRNAs in LR, suggesting that DE circRNAs might regulate tissue regeneration and cell proliferation through Ras, p53, Wnt, Jak-STAT and MAPK signalling pathways during the proliferation phase of LR. And then, four kinds of circRNAs, including circ_03848, circ_08236, circ_13398 and circ_15013 were considered as the candidate key circRNAs and speculated to regulate cell proliferation at 12-72 h after PH by sponging miRNAs and modulating the downstream mRNAs, which based on the number of the corresponding miRNAs and their role enriched and reported in cell proliferation of LR or HCC. It will provide the novel information for LR. In the further studies, using gene interference or additive techniques to verify our conclusion that circRNAs regulate the process of LR are necessary.

Disclosure statement
The authors report no conflict of interest  qRT-PCR, high-throughput sequencing. The infinite was represented by 9 and the infinitesimal was represented by 0.3 for the high-throughput sequencing data.