In summary, we used the RNA-seq data (total RNA with ribosomal RNA being depleted) from six adult and fetal normal tissues (colon, heart, kidney, liver, lung, and stomach) and four gland adult normal tissues (adrenal gland, mammary gland, pancreas, and thyroid gland) to identify the circRNAs using the tool UROBORUS.
CircRNA abundance in human adult normal tissues
In human adult normal tissues, 22,260 unique circRNAs were identified, among which 8120 circRNA candidates were supported by at least two reads spanning a head-to-tail splice junction (Fig. 1a). When we searched the previously published circRNAs deposited in database circBase [31], we found that 9332 circRNAs have been included in the circBase while a total of 12,928 were novel circRNAs (Additional file 1: Table S2-S3). These novel circRNAs included 3092 that were supported by at least two junction reads, that is, circRNAs with high confidence. We noticed that approximately 88% circRNAs were supported by less than 10 reads. Meanwhile, there were 10 circRNA candidates supported by more than 150 reads in all adult samples. Among them, circSLC8A1 (hsa_circ_0000994) had the largest number of supporting reads (823 reads).
We then examined the number of exons for the circRNAs having at least 2 supporting junction reads. We found that majority of the circRNAs (77.74%) had less than 6 exons, indicating that only a small proportion of circRNAs consist of multiple exons by back-splicing (Fig. 1b). This result is consistent with the previous report that circRNAs biogenesis may follow two mechanisms: “direct back-splicing” and “exon skipping” [1, 32, 33]. We manually checked these identified circRNAs’ parental gene using Integrative Genomics Viewer (IGV). The parental genes for those circRNAs with two or more exons did not have corresponding exon-skipping phenomenon. However, the results showed that the derived gene of 35 circRNAs (35/495, 7.07%) had corresponding “exon skipping” transcripts; these circRNAs had only one exon. These results implied that these circRNAs with only one exon might be formed through the “exon skipping” mechanism, while the other circRNAs might be formed by “direct back-splicing” or other mechanism.
For all the adult normal tissues except stomach, the detected circRNAs were distributed on all the human chromosomes (autosomes and sex chromosomes). In stomach tissue, circRNAs were distributed on all the chromosomes but not on Y chromosome (Fig. 1c). In each tissue, our results showed that the number of circRNAs distributed on chromosome 1 was larger than any other chromosomes, while the number of circRNAs on Y chromosome was the smallest among all the chromosomes. This observation is consistent with the chromosome length – Y chromosome is the shortest and has the smallest number of genes in the human genome. In our comparison of the number of circRNAs versus chromosome length (in Megabase, Mb), we found that the circRNA density on chromosome 19 (4.69 per Mb) stood out, which had the highest density among all the chromosomes (average 2.55 per Mb). These results are interesting because the gene density of chromosome 19 is also the highest among chromosomes [34]. However, the distribution of circRNAs on autosomes and X chromosome is not uniform in each normal tissue. For example, the number of circRNAs on chromosome 6 (98) and chromosome 12 (83) in liver was higher than that of other tissues (mean values: 65.2 and 52.4). The number of circRNAs on chromosome 11 in stomach (95) was also higher than that of other tissues.
More than one thousand circRNAs (with at least two supporting junction reads) were identified in each human adult normal tissue (Fig. 2). Each tissue had a good number of unique circRNAs, for example, 441 (36.97%) in colon, 602 (50.04%) in heart, 474 (40.86%) in kidney, 668 (49.34%) in liver, 452 (36.92%) in lung, and 528 (40.80%) in stomach. Moreover, there were 141 circRNAs shared in all normal tissues, even though the expression level of circRNAs was different among these tissues. The result supported that some circRNAs are expressed in all tissues (e.g. in “housekeeping-like” role). These 141 circRNAs may have same or similar function in different adult normal tissues.
Although we found that the number of circRNAs detected in each tissue type was nearly equivalent, the expression of the top ten circRNAs differed in tissue types, as illustrated in Additional file 2: Figure S1. Specifically, circHIPK3 (hsa_circ_0000284) and circUBXN7 (hsa_circ_0001380) were the top two ranked based on the expression in each tissue. The circHIPK3 expression, measured in RPM (junction reads per million mapped reads), was the highest in colon, kidney, lung and stomach tissues, and the average expression level of circSLC8A1 in heart tissue and the circZKSCAN1 (hsa_circ_0001727) in liver tissue was higher than other tissues, suggesting that the circSLC8A1 and circZKSCAN1 may respectively have special function in heart and liver than any other tissues.
There were a total of 33 circRNAs ubiquitously expressed in each sample of all adult normal tissues (Additional file 1: Table S4-S5). Figure 3 displays the heatmap of commonly expressed circRNAs (in RPM) and their parental mRNA expression level (FPKM – fragments per kilobase of exon per million mapped reads). Of note, the expression of these circRNAs in stomach_#4 sample (average 0.04) and stomach_#5 sample (average 0.08) was lower than other tissues, and the circRNAs expressed in liver_#2 sample (average 0.18) was higher than other tissues. Moreover, these 33 circRNAs expressed in different samples of the same tissue were quite different, but the expression of the corresponding parental mRNAs in different individuals of tissue was pretty similar (Pearson correlation coefficient (PCC) = 0.80–0.95), except in lung tissue (PCC = 0.35).
As shown in the results above, these 33 circRNAs exert interesting features. Considering that circRNAs can play a role as miRNA ‘sponges’ to regulate gene expression, we selected 2588 published, high confidence human miRNAs from miRBase, and used the miRanda pipeline to test which miRNA can interact with these 33 circRNAs (more details were included in Materials and methods). Most circRNAs could bind multiple miRNAs and different circRNAs might have same potential miRNA targets. Therefore, we built a circRNA-miRNA interaction network, which contains 33 circRNAs and 158 miRNAs to reflect its complex interaction and regulation (Fig. 4). Next, we employed the dataset from miRTarBase [35], which contained the experimentally validated miRNA-target interactions to predict the potential role for circRNAs in molecular regulation. Then, we used network visualization software Cytoscape (version 3.3.0) [36] to display the network, and through the network analysis, we found a total of 17 circRNAs, 22 miRNAs and 90 mRNAs in circRNA-miRNA-mRNA interaction network. In the network, each gene corresponds to a node, and 2 genes are connected by an edge, indicating a tight correlation between those genes and a potential regulatory relationship (Fig. 5).
To further explore the potential function of these 17 circRNAs in organ development or differentiation, we analyzed the GO enrichment of these 90 mRNAs-derived genes in network (Additional file 3: Figure S2). We found the genes were enriched in the regulation of cell death or apoptosis. In GO molecular function category, these circRNAs were mostly enriched in transcription activity. There is no tissue-specific function detected in the cluster. Therefore, these results might explain why these circRNAs could be shared among the tissues we examined.
CircRNAs expression profile in human fetal normal tissues
In this work, we systematically analyzed 12 samples of human fetal normal tissues, including colon, heart, kidney, liver, lung, and stomach. We found a total of 47,278 circRNAs in human fetal tissues, among which 25,933 circRNAs were supported by at least two junction reads (Additional file 1: Table S6). Compared with the previously published databases obtained from circBase, 31,158 circRNAs were novel circRNA, including 14,241 circRNAs were supported by at least two reads. As shown in Fig. 6a, we identified 4537 circRNAs in colon, 7199 circRNAs in heart, 7962 circRNAs in kidney, 5760 circRNAs in liver, 9698 in lung, and 2584 circRNAs in stomach. Among these circRNAs, 637 circRNAs shared in all human fetal normal tissues. Compared to the corresponding human adult tissues, the circRNAs (average 6290) were substantially highly abundant in fetal tissues than that in adult tissues (average 1238) (Fig. 6b). For example, 9698 circRNA candidates were detected in fetal lung tissues, which is nearly an 8-fold increase to the number of circRNAs identified in adult lung tissues.
Compared with the circRNAs distributed in adult normal tissues, there were 1916 circRNAs (more than 10 supporting reads) only expressed in fetal human tissues. We analyzed the GO enrichment of potential target genes derived from the circRNA-miRNA-mRNA network. A total of 487 GO function enrichment terms, among which 56 GO terms (11.4%) are related to the development (Additional file 1: Table S7). For one example, there were 105 genes involved in GO term ‘positive regulation of developmental process’ (GO ID: 0051094), resulting in a statistically significant enrichment (Fisher’s Exact test, p = 2.31E-36).
Furthermore, we examined the circRNA expression and their parental mRNA expression in human fetal normal tissues. As shown in Additional file 4, 5, 6, 7, 8 and 9: Figures S3-S8, each sample had abundant high-expression specific circRNAs when compared to others, which might lead to individual specificity. We further found that there were a small number of circRNAs with high-expression, but the corresponding parental mRNA had low expression or even no expression; the results revealed that the formation mechanism of circRNAs might be different from mRNA.
The expression level of most circRNAs in each fetal tissue was generally higher than the corresponding circRNAs expressed in adult tissue of colon, liver, heart, kidney and lung, respectively. The same is true for the corresponding parental mRNA expression of each sample in different tissues. The heatmap in Additional file 9: Figure S8 showed that the circRNAs expressed in fetal stomach tissue was higher than that expressed in adult stomach tissue. However, their parental mRNA expression in adult stomach samples was higher than the circRNAs’ parental mRNA expression in fetal stomach samples.
Most of the circRNAs expression and their parental mRNA expression followed the same trend, that is, the expressed level was higher in fetal tissues than adult tissues. However, this is not always true for the samples from stomach tissue. The potential function of circRNAs in stomach tissue might be involved with different mechanism compared to other tissues. And this needs further investigation.
CircRNAs expression in human gland normal tissues
Gland is an organ that synthesizes substances for release into the bloodstream or into cavities inside the body or its outer surface. Thus, it plays an important role in the growth and development of human body. Accordingly, we selected four kinds of gland tissues: adrenal gland, mammary gland, pancreas, thyroid gland, to explore the distribution of circRNAs in human.
We found a total of 23,613 circRNAs in these human gland normal tissues, and 14,433 circRNAs were supported by at least two junction reads (Additional file 1: Table S8). As the Fig. 7 illustrates, there are 2311 circRNAs in adrenal gland, 9665 circRNAs in mammary gland, 1791 circRNAs in pancreas, and 3777 circRNAs in thyroid gland. These circRNAs were identified by UROBORUS (Fig. 7a). Among these circRNAs, we found 321 circRNAs shared among the four glands. The number of circRNAs in pancreas tissues was the least among these four gland tissues, but even so, it was still more abundant than that in any human adult tissues (Fig. 2). The number of circRNA candidates recognized in mammary gland tissues was higher than almost all the human adult and fetal tissues, except for fetal lung tissue. On the other hand, the expression level of circRNAs in mammary gland was also higher than that in other gland tissues (Fig. 7b). Previous studies indicated an abundance number of circRNAs expressed in brain, due to the neuronal activity in accordance with the biogenesis mechanism of circRNAs [37]. Thus, more circRNAs expressed in mammary gland than other tissues might connect with the rich innervation of mammary gland [38]. Other factors might contribute to this abundance of circRNAs in mammary gland, and future investigation is needed.
Furthermore, our results revealed 56 circRNAs had shared expression in different gland normal tissues. Most of circRNAs’ expression level was lower than their parental mRNAs’ expression, but six circRNAs (circFKBP8, circATP5C1, circFGFR2, circRAB3IP, circCORO1C, and circHIPK3) presented the opposite expression with their corresponding parental genes. CircFTO and circPICALM present same trend in adrenal gland, mammary gland and thyroid gland, but high expression of their parental mRNA was observed in pancreas (Fig. 8). In addition, there were 14 circRNAs identified in each sample of human adult tissues. Compared with circRNAs expression in adult tissues, these 14 circRNAs in gland tissues expressed higher than that in adult tissues (Fig. 9). Previous work has revealed that circHIPK3 can regulate cell growth by sponging multiple miRNAs [15]. Our finding also revealed that circHIPK3 expressed higher than their parental linear mRNA in each gland sample. These results might be attributed to the special function of glands in human body like secreting hormonal substances to regulate cell growth.