Clinicopathological-Associated Regulatory Network of Deregulated circRNAs in Hepatocellular Carcinoma

Simple Summary Here, we present a novel strategy to identify key signatures of clinically-relevant co-expressed circRNA-mRNA networks in pertinent cancer-pathways that modulate the prognosis of HCC patients, by integrating clinicopathological features, circRNA and mRNA expression profiles. Five master circRNAs were identified and experimentally demonstrated to upregulate proliferate and promote transformation. Through further integration with miRNA-expression profiles, clinically-relevant competing-endogenous-RNA (ceRNA) networks of circRNA-miRNA-mRNAs were constructed. The most up-regulated nodal-circRNA, circGPC3 was experimentally demonstrated to up-regulate cell-cycle, migration and invasion. circGPC3 was found to act as a sponge of miR-378a-3p to regulate ASPM expression and modulate cell transformation. These 5 nodal circRNAs has potential to be good prognostic biomarkers with good prognostic performance. circGPC3 has great potential to be a promising non-invasive prognostic biomarker for early HCC. We have thus demonstrated the robustness of bioinformatically-predicted master circRNAs in clinically-relevant, circRNA-mRNA networks, underscoring the important roles that these identified deregulated key/master circRNAs play in HCC. Abstract Hepatocellular carcinoma (HCC) is one of the most common and lethal cancers worldwide. Here, we present a novel strategy to identify key circRNA signatures of clinically relevant co-expressed circRNA-mRNA networks in pertinent cancer-pathways that modulate prognosis of HCC patients, by integrating clinic-pathological features, circRNA and mRNA expression profiles. Through further integration with miRNA expression profiles, clinically relevant competing-endogenous-RNA (ceRNA) networks of circRNA-miRNA-mRNAs were constructed. At least five clinically relevant nodal-circRNAs, co-expressed with numerous genes, were identified from the circRNA-mRNA networks. These nodal circRNAs upregulated proliferation (except circRaly) and transformation in cells. The most upregulated nodal-circRNA, circGPC3, associated with higher-grade tumors and co-expressed with 33 genes, competes with 11 mRNAs for two shared miRNAs. circGPC3 was experimentally demonstrated to upregulate cell-cycle and migration/invasion in both transformed and non-transformed liver cell-lines. circGPC3 was further shown to act as a sponge of miR-378a-3p to regulate APSM (Abnormal spindle-like microcephaly associated) expression and modulate cell transformation. This study identifies 5 key nodal master circRNAs in a clinically relevant circRNA-centric network that are significantly associated with poorer prognosis of HCC patients and promotes tumorigenesis in cell-lines. The identification and characterization of these key circRNAs in clinically relevant circRNA-mRNA and ceRNA networks may facilitate the design of novel strategies targeting these important regulators for better HCC prognosis.


Introduction
Hepatocellular carcinoma (HCC) is the commonest liver cancer, representing 70-85% of all cases [1] and is the fifth commonest cancer type with a dismal prognosis of~18% 5-year overall survival [2], mainly due to late diagnosis, high rates of postoperative recurrence, metastasis and paucity of early diagnosis biomarkers. Chronic hepatitis-B/C viruses (HBV/HCV) are major risk factors for HCC [3].
Recent high-throughput RNA sequencing reveals that more than 70% of the genome is transcribed to ncRNAs [4,5]. circRNAs are a class of ncRNAs produced by a back-splice event from pre-mRNA [6] and can be classified into five groups according to the relationship with coding RNA in the transcript: "exonic", "intronic", "intergenic", "sense" and "antisense" [7]. While circRNAs can function in different ways at the molecular level, current evidence suggests that circRNAs, particularly exonic circRNAs, mainly act as sponges for miRNA harboring miRNA response elements (MREs) to affect the activity of miRNA-target mRNA interaction [8]. For instance, circMTO1 suppresses HCC progression by sponging oncogenic miR-9 to enhance p21 expression, attenuating HCC cell proliferation and invasion [9]. circRNAs can also act as protein sponges [10] modulating protein function [11], or act as scaffolds for protein [12] recruiting proteins to specific locations [13], or undergo cap-independent translation [14]. As circRNAs are conserved, stable, abundant and diverse [15], they have potential to serve as diagnostic/prognostic biomarkers.
Current studies mainly focus on single circRNA inferred from host genes or selected from deregulation profiles before their tumorigenic roles and/or clinical association were investigated [9,16]. A few studies identified/catalogued lists of circRNAs that were deregulated in only a few HCC patients, identifying 26-87 differentially expressed circRNAs in HCC patients [17,18]. A database of circRNAs was also collated from sequencing data of 5 HCC patients [19]. In addition, circRNA-based competing RNA (ceRNA) networks from limited publicly available datasets comprising between 3 and 7 HCC patients were also constructed [20,21].
As several circRNAs may work together to modulate tumorigenesis, metastasis or even patients' prognosis, it is thus pertinent to identify and characterize networks of cir-cRNA comprehensively as these may function synergistically. Here, we employed a novel network-based strategy to identify clinically relevant, differentially expressed key/master cirRNA regulators through global integration of deregulated cirRNAs, miRNAs and mR-NAs with clinical characteristics and pertinent cancer pathways.

HCC Patients and Clinical Samples
Tumorous tissue, adjacent non-tumorous liver tissue and blood of HCC patients from Singapore General Hospital were obtained after informed consent and prior approval from the SingHealth Institutional Review Board (SingHealth CIRB Ref: 2018/3155) and kept in liquid nitrogen. All methods were carried out in accordance with relevant guidelines and regulations, and the tissues were anonymized to the bench researchers.
The demographics and clinical characteristics of these HCC patients are summarized in Supplementary Table S1. The circRNA, miRNA and mRNA expression profiles of 49 HCC patients' tumorous tissue and adjacent non-tumorous liver tissue were obtained via microarray or sequencing. The expression profiles for circRNAs, miRNAs and mRNAs in this study were deposited in Gene Expression Omnibus with series entries GSE155949, GSE156087 and GSE138178, respectively.

Workflow for Analyses of Clinically Relevant circRNA Networks in HCC
The strategy to identify clinically relevant circRNA-mRNA and circRNA-miRNA-mRNA competing RNA (ceRNA) networks from circRNA, miRNA and mRNA expression profiles of 49 HCC patients is presented in Figure 1A. Clinically relevant deregulated circRNA-mRNA pairs are identified as follows: (1) Significantly deregulated circRNA (|FC| > 1.5, FDR < 0.05) are associated with various clinical characteristics (Table S1) using either a Student's t-test for binary characteristics or a log-ranked test for survival analyses to identify clinically associated deregulated circRNAs ( Figure 1A). (2) The correlation between the expression of circRNAs and mRNAs is then determined using Pearson correlation to identify significant circRNA-mRNA pairs (|r| > 0.7, FDR < 0.05) ( Figure 1A). (3) The clinically associated deregulated circRNAs are then integrated with significantly correlated circRNA-mRNA pairs to identify clinically relevant deregulated circRNA-mRNA pairs ( Figure 1A).
The following steps were performed to generate a circRNA-miRNA-mRNA ceRNA network. (1) circRNA-miRNA and miRNA-mRNA pairs were first identified through integrating data from expression correlation analyses using Spearman correlation as well as prediction algorithms, namely miRanda [22] and PITA [23] to ensure that the circR-NAs/mRNAs were not only predicted to be targets of the miRNAs, but their expressions were also inversely correlated with the corresponding miRNA. (2) The significance of sharing of miRNAs by the different ceRNAs (i.e., circRNA-circRNA, circRNA-mRNA and mRNA-mRNA) was then evaluated using various statistical analyses including the Hypergeometric test [24], Pearson correlation (PC) [24,25], partial Pearson correlation (PPC) [26], sensitive partial Pearson correlation (SPPC) [26], and conditional mutual information (CMI) [27]. (3) The relative confident ceRNA pairs that were found to be significant by PC (r > 0 & p < 0.05) and PPC methods (p < 0.05) were then integrated with clinically relevant deregulated circRNA-mRNA pairs (described earlier) to generate clinically relevant deregulated circRNA-miRNA-mRNA triplets. All of the analysis above was performed by using R and/or Matlab2017a.
All the primers for RNAs, siRNAs for circRNAs and miRNAs, and probes for circR-NAs are listed in Table S2. Further details can be found in the Supplementary Materials and Methods.
In total, 74 of the 132 circRNAs were potential oncogenic circRNAs since high tumor expression was consistently associated with poorer clinical characteristics (e.g., higher tumor grade), while 35 were potential tumor suppressor circRNAs since high expression of these circRNAs in tumors was consistently associated with better clinical characteristics (e.g., lower tumor grade) ( Figure 2A).
Cytoscape software [28] was employed to better visualize the relationship of the 460 clinically associated deregulated circRNA-mRNA co-expressed pairs. In total, 17 oncogenic ( Figure 2E) and 6 tumor suppressive circRNAs ( Figure 2F) clustered separately to form independent networks. The 17 oncogenic circRNAs were significantly co-expressed (|r| > 0.7, FDR < 0.05) with 159 genes forming two networks, namely a singleton (one circRNA-one mRNA) ( Figure 2E(ii)) and a very large and complex network comprising 16 circRNAs and 158 genes ( Figure 2E(i)). Notably, only five circRNAs (circGPC3/circR091581, circW7/circR001387, circW8/circR001388, circW3/circR103583, and circRaly/circR406097) in the large oncogenic circRNA-gene network are predicted to be sufficient to modulate the expression of 93% (147/158) of the genes in the network as these are significantly co-expressed with these 147 mRNAs ( Figure 2E(i),G). In total, 89.8% of the co-expressed genes of these five nodal/master circRNAs were also associated with worse clinical characteristics, including high tumor grade or absent tumor capsule ( Figure 2E(i),G). The six tumor suppressive circRNAs were significantly co-expressed with 20 mRNAs forming six distinct networks as shown in Figure 2F, four of which comprise only one or two genes. Of the two larger networks, co-expressed genes [11] of the largest network ( Figure 2F(i)) were mainly metallothionein genes while the co-expressed genes [4] of the second largest network ( Figure 2F(ii)) were found to be involved in butanoate, fatty acid and amino acid metabolism. These observations suggest that oncogenic and tumor suppressive circRNAcentred networks modulate prognosis of HCC patients through distinct pathways. Notably, the expression of these five nodal/master circRNAs with red and bold labels were also significantly correlated with each other (Figure S1A), suggesting a potential built-in redundancy in pathway regulation. Interestingly, three of the oncogenic nodal circRNAs (circW7; circW8; circW3), together with four other circRNAs (circR103586, circR103582, circR103584, and circR103587) in the 16 circRNA-centred oncogenic networks and circR103585 were all derived from the same gene: WHSC1 ( Figure 2E(i) and Figure S2B). In fact, 48.5% of circRNAs were expressed in more than one isoform ( Figure S2A). On the other hand, expression of most of the six tumor suppressive circRNAs was less strongly correlated with each other (|r| < 0.6) ( Figure S1B) leading to less complex networks observed ( Figure 2F).

Figure 2.
Prognostic circRNA-centric network modulating cell cycle and metabolism: (A) Prognostic RNA (e.g., circRNA) is defined as RNA that is differentially expressed in tumorous tissue vs. adjacent non-tumor, and is associated with worse/better clinical characteristics. Worse prognostic RNA is RNA whose expression is higher in tumors and the higher tumor expression is associated with clinical outcomes (e.g., higher grade tumors, absence of tumor capsule, worse overall survival) while better prognostic RNA is one whose expression is lower in tumors and the low tumor expression is associated with worse clinical outcomes. Examples of worse (left) and better prognostic cirRNAs (right) associated with tumor grade To determine whether circRNA functions in cis through the parental mRNA/transcript or in trans through modulating other mRNAs, we evaluated whether the expression changes of the circRNAs were consistent with the host mRNAs. The expression changes of circR-NAs and host genes between tumorous versus non-tumorous tissue were overall positively correlated (r = 0.277) ( Figure S3A, quadrant (ii) and (iii)), although~2.5% of these circRNAs and their parental mRNA/transcripts have differential changes (Supplementary Figure S3A, quadrant (i) and (iv)), as evident from the examples given in the figure ( Figure S3B).

Oncogenic Roles of the Five Nodal/Master circRNAs Associated with Worse Clinical Characteristics
The oncogenic roles of the five nodal/master circRNAs (circGPC3, circW7, circW8, circW3, and circRaly) were experimentally evaluated since these key circRNAs are associated with clinical characteristics related to poorer patient prognosis and predicted to modulate the expression of~93% (147/158) of the genes associated with the cell cycle pathway. Expression of all five nodal circRNAs were significantly associated (p < 0.05) with tumor grade with lowest expression in the non-tumorous tissue and highest expression in the higher grade tumors ( Figure 3A). Expression of three of the five nodal circRNAs (circW7, circW8, and circW3) were also significantly associated (p < 0.05) with tumor capsulation, with non-tumorous tissue exhibiting the lowest expression, while the tumors with incomplete capsule, indicative of poorer prognosis, exhibited the highest expression ( Figure 3A). Higher expression of circRaly was also significantly associated with poorer overall survival ( Figure 3A). The higher expression of these nodal circRNAs in the tumors was validated using RT-qPCR of samples in independent cohort of HCC patients ( Figure 3B).
To further characterize the roles of these nodal circRNAs in HCC, we evaluated their expression in a non-transformed liver cell-line, LO2, and a transformed cell-line, Huh7. All five circRNAs were expressed at higher levels in Huh7 compared to LO2 ( Figure 3C). The junction of all five nodal circRNAs was identified through sequencing of the amplified products using the divergent primers ( Figure S4A). All five nodal circRNAs are resistant to Rnase R digestion ( Figure S4B), suggesting that the expression observed is circular RNA and not its linear host gene. All five nodal circRNAs were found to be primarily localized to cytoplasm as evident from sub-cellular fractionation ( Figure S4C). RNA immunoprecipitation (RIP) with Ago2-antibodies on Huh7 cells, which expresses high levels of all five nodal circRNAs, revealed that all five circRNAs can interact with Ago2 ( Figure S4D), suggesting that these may be silenced by miRNAs since Ago2 protein are essential components of the RNA-induced silencing complex (RISC) and play key roles in RNA silencing. This may also be evidence of their sponging ability since Ago2 is also involved in miRNA sponging. Data are from three independent experiments. Mean ± SEM (n.s.: not significant; * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001 by two-tailed Student's t-test).

The Five Nodal circRNAs Modulate Tumorigenesis by Promoting Proliferation and Anchorage Independent Growth
Introduction of circGPC3, circW3, circW7 and circRaly constructs into LO2 (the nontransformed cell-line) led to increased soft-agar colony formation, while attenuating the expression of circW8 in Huh7 led to significantly less colony formation ( Figure 3D,E), confirming the oncogenic role of these five nodal circRNAs. Overexpressed circW3, circW7 and circRaly were also found to promote soft-agar colony formation in Huh7 ( Figure S4E).
Overexpressed nodal circGPC3, circW3 and circW7 in LO2 were found to significantly increase the cell proliferation, leading to shortened doubling time ( Figure 3D,F), while knockdown of circW8 in Huh7 significantly inhibited cell proliferation, increasing their doubling time ( Figure 3D,F). A similar trend was observed for circRaly, although it did not reach statistical significance ( Figure 3D,F and Figure S4F).
RNA sequencing was performed for the five key nodal circRNAs that were overexpressed or inhibited in non-transformed LO2 and/or transformed Huh7. In total, 98 upregulated and 96 down regulated genes were co-regulated by all five key nodal circRNAs (Table S4). These were predicted to mainly co-upregulate signal transduction/immune system pathways ( Figure S5A), and co-downregulate organelle biogenesis/maintenance and metabolism pathways ( Figure S5B). Furthermore, the key nodal circRNAs were also found to upregulate genes involved in the deadenylation of mRNA ( Figure S5A), which is implicated in miRNA function [34], affirming their roles as sponges for miRNAs.
3.6. Clinically Relevant Differentially Expressed circRNA-miRNA-mRNA ceRNA Networks circRNAs can modulate gene expression through various modalities, and one of the major ways they do so is by acting as competing endogenous RNAs (ceRNAs) or miRNA sponges/inhibitors to modulate the expression of mRNAs by competing for shared miRNAs. ceRNA networks were thus computationally inferred from circRNA, miRNA and mRNA expression profiles. Correlation (Spearman) analyses coupled with prediction (miRanda and PITA) algorithms ( Figure 1A) were employed to identify 773 shared circRNA-miRNA pairs and 22,233 shared miRNA-mRNA pairs ( Figure 4A, red circle). The significance of sharing the same miRNA response element (MRE) by circRNA/mRNA for each ceRNA pair was evaluated using the hypergeometric test [24]. A total of 784, 838,878 and 27,107 significant circRNA-circRNA, mRNA-mRNA and circRNA-mRNA pairs, respectively, were identified. In total, 0, 1679 and 41 highly confident circRNA-circRNA, mRNA-mRNA and circRNA-mRNA ceRNAs respectively were identified through integrating PC [24,25], PPC [26], SPPC [26] and CMI [27] statistical analyses ( Figure 4B, red circle).
The 41 circRNA-mRNA ceRNA pairs generated 406 potential circRNA-miRNA-mRNA triplets (Table S5). In total, 12 of the genes in these triplets which were upregulated in the tumors of HCC patients were mainly enriched in the cell-cycle, signaling pathways and carbohydrate metabolism, while 26 genes in these triplets were downregulated and enriched in genes associated with platelet activation, signaling, aggregation and degranulation ( Figure 4C). None of the 5 nodal circRNAs in the key oncogenic clinically relevant circRNA-mRNA network ( Figure 2E), or any of the other circRNAs in the same network, were amongst the 41 highly confident circRNA-mRNA ceRNA pairs. We thus further evaluated if any of the five nodal circRNAs are predicted to compete with their co-expressed mRNA for miRNA binding at a lower stringency, by identifying moderately confident ceRNAs. The moderately confident ceRNAs were identified through the integration of only the PC and PPC statistical analyses, and 114,419 moderately confident circRNA-miRNA-mRNA triplets were derived. The 114,419 moderately confident circRNA-miRNA-mRNA triplets were intersected with the 460 clinically relevant circRNA-mRNA pairs (|r| > 0.7) ( Figure 2B and Table S3) to yield 34 circRNA-miRNA-mRNA triplets in the clinically relevant circRNA-miRNA-mRNA ceRNA network ( Figure 4D,E and Table S6). Notably, one of the nodal circRNAs in the oncogenic clinically relevant circRNA-mRNA network, circGPC3 (circR091581) (Figure 4Dii, red circle), was found to compete with 11 mRNA/genes (ASPM, CENPW, KIF14, NEK2, POLQ, TOP2A, DBF4, ERCC6L, E2F7, GPC3, and MMS22L) involved in cell-cycle for two miRNAs (miR-378c and miR-378a-3p) ( Figure 4F).

Role of circGPC3 in Modulating Tumorigenesis through Sponging miR-378a-3p to Upregulate ASPM
Of the five nodal circRNAs, circGPC3 was predicted in silico to reside in a clinicallyrelevant circRNA-miRNA-gene ceRNA network and compete with several mRNAs for binding to two miRNAs (Figure 4Dii). Its host gene, GPC3 (Glypican 3), was implicated in tumorigenesis through the Wnt signaling pathway [35]. Hence, we further characterized cirGPC3 to understand how circGPC3 modulates cancer phenotypes leading to higher tumor grade and poorer cancer prognosis. circGPC3 is derived from exon 3 of GPC3 (chrX:132887508-132888203). Expression of circGPC3 and its host gene, GPC3, were upregulated in the tumor tissue compared to the adjacent non-tumorous liver tissue in both the 47 discovery phase HCC samples and 15 independent cohort HCC samples ( Figure S6A,B and Figure 3B).
Overexpression of circGPC3 in LO2 cells and inhibition of circGPC3 in Huh7 cells upregulated the expression of 467 genes in the Toll-like receptors (TLRs) and transcription pathways (Table S7, Figure S7A), while downregulating the expression of 409 genes in the transcription and DNA repair pathways (Table S7, Figure S7B).
Interestingly, modulating the expression of circGPC3 altered the protein ( Figure S6D) but not the transcript ( Figure S6C) expression of its host gene GPC3. Inhibiting circGPC3 expression not only inhibited GPC3 protein expression, but it also inhibited β-catenin expression ( Figure S6D), which is consistent with previous observation of β-catenin acting downstream of GPC3 [36]. circGPC3 is predicted in silico to have potential to be translated ( Figure S6E) as a peptide, raising the possibility that circGPC3 may modulate phenotypes as a protein.
circGPC3 can be amplified using divergent primers from cDNA but not genomic DNA, suggesting that it is in the RNA form ( Figure 5A). Fluorescence in-situ hybridization (FISH) revealed that circGPC3 mainly resides in the cytoplasm ( Figure 5B). Live cell imaging showed that inhibition of circGPC3 expression in Huh7 significantly inhibited cell proliferation increasing their doubling time, while overexpression of circGPC3 in LO2 significantly increased the cell proliferation ( Figures 3F and 5C,D). These observations are consistent with the predicted oncogenic role of circGPC3. Cell-cycle analyses revealed that introduction of circGPC3 in LO2 led to more cells entering S/G2 phases and fewer in G0/G1 phases (although not significant), while the reverse is observed in Huh7 cells whose circGPC3 expression is inhibited with siRNAs ( Figure 5E). Notably, attenuating the expression of circGPC3 in Huh7 cells led to significantly less soft-agar colony formation ( Figure 5F). Inhibition of circGPC3 expression in Huh7 led to reduced invasive ability ( Figure 5G) and slower cell migration ( Figure 5H), while over-expression of circGPC3 in non-transformed LO2 cells led to enhanced invasive capability ( Figure 5G) and more rapid cell movement ( Figure 5H). These observations in Huh7/LO2 cells were also validated in HepG2 and SNU449 cells ( Figure S8A-F), except that the attenuation of circGPC3 expression in HepG2 cells led to faster wound healing ( Figure S8F). This may be due to inhibition of circGPC3 expression upregulating genes in the VEGFA pathway in HepG2 but downregulating different genes in the same VEGFA pathway in Huh7 cells (Table S8). Nonetheless, as circGPC3 also plays a role in cell proliferation, the observation may also be due to the effect of cell proliferation. Further experiments with mitomycin C to eliminate the effect of cell proliferation on wound healing will clarify this.
As circGPC3 was predicted in silico to reside in a clinically relevant circRNA-miRNAgene ceRNA network and compete with 11 mRNAs for binding to miR-378a-3p (Figure 4Dii), we thus evaluated if the expression of these genes is affected by circGPC3. Using quantitative real-time RT-PCR (qRT-PCR), the trend of a few mRNAs is consistent with the observation in patients ( Figure 6A) where their expressions are upregulated in LO2 cells transfected with circGPC3, and downregulated in Huh7 cells transfected with siRNA against circGPC3. As the ASPM gene showed the highest upregulation in LO2 and deepest downregulation in Huh7 cells, we further characterized the relationship between ASPM and circGPC3 experimentally. We hypothesized that circGPC3 competes with ASPM for binding to miR-378a-3p, leading to the modulation of cellular cancer phenotypes ( Figure 6A).
We initially examined if miRNA-378a-3p is a target of both circGPC3 and ASPM by genetically engineering luciferase reporter constructs carrying the following: wild-type circGPC3 (circGPC3 WT); mutant circGPC3 (circGPC3 Mut) with the single predicted miRNA-378a-3p binding site altered ( Figure 6B, top panel); wild-type ASPM (ASPM WT); and mutant ASPM (ASPM Mut) with all three predicted miRNA-378a-3p binding sites in the protein coding sequence (CDS) altered in the same construct ( Figure 6B, top panel).
These constructs were then co-transfected with either control miRNA or miRNA-378a-3p into HEK293K cells and luciferase activity was determined. As evident in Figure 6B (bottom panel), miRNA-378a-3p inhibited luciferase activity of cells carrying the WT forms of circGPC3 and ASPM but not the mutant forms of both circGPC3 and ASPM where the binding to miR-378a-3p was abrogated. Hence, miRNA-378a-3p targets both circGPC3 and ASPM to inhibit their expression.
As Ago2 is an essential component of the RNA-induced silencing complex (RISC) and plays key roles in RNA silencing, we then determined if both circGPC3 and miR-378a-3p can bind to Ago2 as part of this complex. RNA immunoprecipitation (RIP) was performed on Huh7 cells, which expresses high levels of circGPC3. Endogenous circGPC3 and miR-378a-3p were immunoprecipitated from Huh7 cells with Ago2-antibodies (Sigma Aldrich, Saint Louis, MO, USA) and their expression was analysed using qRT-PCR. Figure 6C shows that both circGPC3 and miR-378a-3p were precipitated with the Ago2-antibodies, suggesting that they are in the RISC complex. We then proceeded to determine whether circGPC3 directly binds miR-378a-3p. As shown in Figure 6D, higher enrichment of circGPC3 and miR-378a-3p transcripts was observed, when Huh7 and circGPC3-expressing LO2 celllysates were probed with 3 -terminal-biotinylated-circGPC3 compared to control probes (Huh7 cells) or vector control (circGPC3-expressing LO2 cells), suggesting that circGPC3 interacted directly with miR-378a-3p. We then evaluated whether circGPC3 modulates the hallmarks of cancer in LO2 cells, through sponging miR-378a-3p, leading to the up-regulation of ASPM expression. As shown in Figure 6E (top panel), the transformation potential of circGPC3-expression LO2 was significantly attenuated in the presence miR-378a-3p (histogram bars 1 versus 3) but not the control miRNA (histogram bars 1 versus 4). Similar trends were observed with cell proliferation ( Figure S9A) and cell migration ( Figure S9B) experiments in LO2 cells. Cells expressing circGPC3 were found to express higher levels of ASPM at the transcript levels ( Figure 6E, middle panel, histogram bars 1 versus 2) as well as the protein levels ( Figure 6F, right panel) in LO2 cells. On the other hand, inhibiting circGPC3 in Huh7 cells with siRNA against circGPC3 inhibited ASPM protein expression ( Figure 6F, left panel). These data are consistent with previous reports that ASPM is associated with vascular invasion, early recurrence, and poor prognosis in HCC [37] and suggests that circGPC3 may modulate cancer phenotypes through the miR-378a-3p-ASPM pathway. Figure 6. CircGPC3-centred ceRNA regulation through circGPC3/miR-378a-3p/ASPM axis: (A) Right: Validation of genes predicted to correlate with circGPC3 in clinically relevant circRNA-miRNA-mRNA ceRNA network (see Figure 5D). Left: Hypothesis of the role of circGPC3 in modulating various oncogenic phenotypes (e.g., cell proliferation, cell cycle and cell transformation) through the circGPC3/miR-378a-3p/ASPM axis. (B) Top: Sequence of Mutation (in red) generated for circGPC3 Mut and ASPM Mut (at the miRNA binding site within the protein coding sequence). Bottom: Luciferase assays with reporter constructs containing the wild-type or mutant circGPC3/ASPM downstream of a luciferase gene were performed after co-transfection with miR-378a-3p in HEK293T cells. (C) RNA Immunoprecipitation (RIP) experiments were performed using an antibody against Ago2 on extracts from Huh7 cells. Expression of circGPC3 and miR-378a-3p was detected. (D) Huh7 and circGPC3-expressing LO2 cell-lysates were probed with 3 -terminal-biotinylated-circGPC3 compared to control probes (Huh7 cells) or vector control (circGPC3-expressing LO2 cells). (E) circGPC3 promoted the cell transformation via the miR-378a-3p-ASPM pathway. (F) Protein level of ASPM after circGPC3 KD or OE. Data are from three independent experiments. Mean ± SEM (n.s.: not significant; * p < 0.05; ** p < 0.01; *** p < 0.001; by two-tailed Student's t-test).

Significance of Clinically Relevant circRNA-Centric Regulatory Network in HCC
In summary, we propose a model for a clinically relevant circRNA-centric regulatory network in HCC ( Figure 7A) by integrating clinical characteristics with circRNA, miRNA and mRNA expression computationally and validating experimentally. From the receiver operating characteristic (ROC) analysis, all five key nodes showed excellent performance for distinguishing tumor grade (3,4) groups from adjacent nontumorous tissue (p < 0.001, AUC > 0.9) ( Figure 7B). circW8 and circRaly showed good performance for distinguishing tumor grade (3,4) from tumor grade (1,2) groups (p < 0.001, AUC > 0.8), while circW3 and circW7 showed fair performance (p < 0.01, AUC > 0.7) ( Figure 7B). Although circGPC3 showed poor performance for distinguishing tumor grade (3,4) groups from tumor grade (1,2) groups (p < 0.05, AUC = 0.69), it was able to reasonably distinguish tumor grade (1,2) groups from adjacent non-tumorous tissue (p < 0.001, AUC > 0.8), similarly to circW3 and circW7 ( Figure 7B). Logistic regression analysis showed that combination of circGPC3, circW3, circW7 and circRaly can excellently distinguish tumor grade (1,2) groups from adjacent non-tumorous tissue (p < 0.001, AUC = 0.92) with a high sensitivity of 0.82 and a high specificity of 0.95 ( Figure 7C, bottom left), which suggested the potential roles of the key nodes in the early detection of HCC. Additionally, the combination of circGPC3, circW7, circW8 and circRaly had good performance for distinguishing tumor grade (3,4) groups from tumor grade (1,2) groups (p < 0.001, AUC = 0.88) ( Figure 7C, bottom right). Hence, these nodal circRNAs have potential to be good prognostic biomarkers. Notably, the key nodal circGPC3 circRNAs were found to be a promising non-invasive biomarker for detecting HCC or even early HCC, as its expression is highest in the plasma of patients with higher grade tumors (3,4), compared to those with lower grade tumors (1,2), and lowest in healthy individuals ( Figure 7D).

Discussion
Current limitations in HCC diagnosis/therapy necessitates the identification of novel molecules through better elucidation of the molecular network mechanisms of HCC. Here, based on circRNA, mRNA and miRNA expression profiles of HCC patients, clinically relevant circRNA-centric networks were computationally predicted. As illustrated in Figure 7A, the prognostic significance of these predicted circRNA-centric networks was evident from key cancer pathways molecules in these networks are predicted to reside in; the experimental demonstration of tumorigenic/metastatic phenotypes of the master/key nodal circRNAs; as well as the relevant cancer-related mechanism, through which the nodal circGPC3 cirRNA acts as a sponge of miR-378a-3p to alter the expression of ASPM (a cancer molecule [37]), leading to the modulation of tumorigenesis.
Seventeen circRNAs residing in two circRNA-mRNA networks were associated with clinical characteristics that relate to worse prognosis in HCC patients (higher tumor grade and less tumor encapsulation). Six circRNAs, residing in six separate circRNA-mRNA networks, were associated with clinical characteristics of better prognosis (smaller tumor size, absent tumor capsule and better overall survival). Encapsulated HCC tumors are less aggressive, with lesser invasion events [38]. Tumor grade was reported to be associated with deregulation of proliferation, cell cycle, metastasis and invasion-associated genes in HCC [39][40][41][42].
Five key/master nodal circRNAs (circGPC3, circW3, circW7, circW8 and circRaly), accounting for~93% (147/158) of the co-expressed genes within the large 16-member clinically relevant worse prognosis circRNA-mRNA network, were identified. Notably, all five master circRNAs are significantly correlated with one another (Figure S1A), suggesting a potential built-in redundancy of circRNAs in pertinent pathway regulation. Interestingly, three of these (circW3, circW7 and circW8) are derived from the same parental gene: WHSC1 ( Figure S2B). In fact, another circRNA (circ-NSD2), also derived from WHSC1, was reported to promote migration, invasion and metastasis of colorectal cancer (CRC) cells in a mouse model of CRC liver metastasis [43], highlighting that this family of circRNAs derived from WHSC1 may be worth further investigation.
To facilitate the design of therapeutic strategies for HCC based on circRNAs in key cancer-related circRNA-mRNA networks, it is pertinent to elucidate the molecular mechanisms by which circRNAs in these cirRNA-mRNA networks modulate tumorigenesis. While circRNAs can modulate gene expression through diverse ways, one of the major ways is by acting as ceRNAs/miRNA sponges to influence gene expression through competition for shared miRNAs [44]. Here again, through computational analyses, we identified a confident clinically relevant circRNA-miRNA-mRNA network that hosts circGPC3, one of the five key nodal circRNAs.
circGPC3 is the most highly upregulated circRNA in the tumors of HCC patients, and its high expression is significantly associated with a high tumor grade. High circGPC3 expression promoted cell proliferation, cell cycle, cell transformation, migration and invasion in both transformed and non-transformed liver cell lines ( Figures 3E and 5C-H). circGPC3 was computationally predicted to modulate expression of 11 genes by competing with two miRNAs, eight of which were experimentally validated ( Figure 6A). Notably, we were able to experimentally demonstrate the circGPC3/miR-378a-3p/ASPM axis through luciferase assays, Ago2, RIP and circRNA pull-downs, as well as cancer phenotype assays ( Figure 6). Our data are consistent with previous reports that both miR-378a-3p [45] and ASPM [37,46,47] can modulate cellular proliferation, migration and invasion. In fact, ASPM was previously reported to be over-expressed in HCC, and was highlighted as a biomarker associated with poorer prognosis in HCC [37,48] and breast cancer [49]. This study thus highlights circGPC3 as the link between miR-378a-3p and ASPM in their role in tumorigenesis and modulating HCC prognosis.
The ability of circRNAs (e.g., circGPC3) to sponge miRNAs (e.g., miR-378a-3p), to modulate tumorigenesis, suggests their potential as prognostic biomarkers and therapeutic targets [50]. In fact, the nodal circRNAs, circGPC3, circW3 and circW7 displayed good performance for distinguishing low tumor grades from adjacent non-tumorous tissue (p < 0.001, AUC > 0.8) ( Figure 7B), suggesting the potential roles of these key nodal cir-cRNAs in the early detection of HCC. circRNAs represent attractive minimally invasive biomarkers as they are present in various bodily fluids (including blood and saliva [51]) and are insensitive to ribonuclease, hence they are more stable than other linear RNA molecules due to the special circular closed structure [52]. Notably, plasma levels of cir-cGPC3 were lowest in healthy individuals and highest in HCC patients with higher grade tumors ( Figure 7D), suggesting their potential as promising minimally invasive prognostic biomarkers for HCC, although validation in a larger cohort of samples is necessary.
circRNAs also represent attractive targets for therapy, as targeting the unique backsplice junction of oncogenic circRNAs led to higher specificity with less off-target effects [44]. Interestingly, circGPC3 upregulated the protein but not the mRNA expression of its host gene GPC3 ( Figure S6C,D). circGPC3 was predicted by circBank [53] to have potential to be translated into protein ( Figure S6E), suggesting the circGPC3 may be bifunctional. On the one hand, circGPC3 can act at the RNA level to modulate the expression of its co-expressed genes (e.g., ASPM), perhaps as a miRNA sponge ( Figure 6). On the other hand, we hypothesize circGPC3 may be translated into protein and influence the protein expression ( Figure S6D) of its host GPC3 gene, perhaps through affecting protein stability. In fact, high GPC3 protein expression was reported to be associated with poor prognosis of HCC patients [54]. circGPC3 was also found to alter the protein expression of β-catenin ( Figure S6D), which is consistent with reports that GPC3 can activate the targeted therapeutic Wnt signaling [55]. GPC3 is an established oncofetal glycoprotein that is highly overexpressed in HCC. It is hailed not only as a diagnostic [56] and prognostic [57,58] biomarker, but also as a novel, attractive therapeutic target that is currently in clinical trials [56,59]. It is thus tempting to propose circGPC3, which acts upstream GPC3, as another attractive prognostic biomarker which can also be developed as another novel therapeutic target. Further studies are necessary to validate if circGPC3 can act as a protein to modulate its host GPC3 protein expression as well as the Wnt pathway.

Conclusions
In summary, these data suggest the robustness of bioinformatically-predicted clinically relevant circRNA-mRNA networks and key/master nodal circRNAs, and underscore the important roles that these identified deregulated key/master circRNAs play in HCC. It also highlights the robustness of our computational strategy to identify clinically relevant circRNA-miRNA-mRNA ceRNA networks, which can be experimentally validated, e.g., circGPC3/miR-378a-3p/ASPM regulated axis. Indeed, circGPC3 has great potential to be a promising non-invasive prognostic biomarker for early HCC.

Data Availability Statement:
The expression profiles for circRNAs, miRNAs and mRNAs in this study are deposited in Gene Expression Omnibus with series entry GSE155949, GSE156087 and GSE138178, respectively.

Conflicts of Interest:
The authors declare no conflict of interest.