Construction and Comprehensive Analysis of a circRNA-miRNA-mRNA Regulatory Network to Reveal the Pathogenesis of Hepatocellular Carcinoma

Background: Circular RNAs (circRNAs) have been demonstrated to be closely related to the carcinogenesis of human cancer in recent years. However, the molecular mechanism of circRNAs in the pathogenesis of hepatocellular carcinoma (HCC) has not been fully elucidated. We aimed to identify critical circRNAs and explore their potential regulatory network in HCC. Methods: The robust rank aggregation (RRA) algorithm and weighted gene co-expression network analysis (WGCNA) were conducted to unearth the differentially expressed circRNAs (DEcircRNAs) in HCC. The expression levels of DEcircRNAs were validated by quantitative real-time polymerase chain reaction (qRT-PCR). A circRNA-miRNA-mRNA regulatory network was constructed by computational biology, and protein-protein interaction (PPI) network, functional enrichment analysis, survival analysis, and infiltrating immune cells analysis were performed to uncover the potential regulatory mechanisms of the network. Results: A total of 22 DEcircRNAs were screened out from four microarray datasets (GSE94508, GSE97332, GSE155949, and GSE164803) utilizing the RRA algorithm. Meanwhile, an HCC-related module containing 404 circRNAs was identified by WGCNA analysis. After intersection, only four circRNAs were recognized in both algorithms. Following qRT-PCR validation, three circRNAs (hsa_circRNA_091581, hsa_circRNA_066568, and hsa_circRNA_105031) were chosen for further analysis. As a result, a circRNA-miRNA-mRNA network containing three circRNAs, 17 miRNAs, and 222 mRNAs was established. Seven core genes (ESR1, BUB1, PRC1, LOX, CCT5, YWHAZ, and DDX39B) were determined from the PPI network of 222 mRNAs, and a circRNA-miRNA-hubgene network was also constructed. Functional enrichment analysis suggested that these seven hub genes were closely correlated with several cancer related pathways. Survival analysis revealed that the expression levels of the seven core genes were significantly associated with the prognosis of HCC patients. In addition, we also found that these seven hub genes were remarkably related to the infiltrating levels of immune cells. Conclusion: Our research identified three pivotal HCC-related circRNAs and provided novel insights into the underlying mechanisms of the circRNA-miRNA-mRNA regulatory network in HCC.


INTRODUCTION
Primary liver cancer is a common malignant tumour and the third leading cause of cancer-related deaths worldwide, threatening our health and life. Hepatocellular carcinoma (HCC) is the most common type of liver cancer, accounting for 75-80% of cases (Sung et al., 2021). Although great improvements have been achieved in the diagnosis and treatment of HCC, most patients diagnosed with advanced HCC have limited treatment options and their 5-years survival rate is less than 15% (Bertuccio et al., 2017). Besides, HCC patients with surgical resection are still suffering from high risk of recurrence and lung metastasis (Kulik et al., 2018). Given the great perniciousness of HCC, there is an urgent need to explore the molecular cascades of the development and progression of HCC.
Circular RNAs (circRNAs), a class of endogenous non-coding RNAs, were first discovered in viroids in 1976 (Sanger et al., 1976). Unlike traditional linear RNAs, circRNAs are characterized by their covalently closed-loop structures with neither the 5′ end cap nor 3′ poly (A) tails, enabling them to be resistant to the digestion of RNase R (Suzuki et al., 2006). In the past, circRNAs were generally considered to be products of aberrant RNA splicing due to their low abundance and lack of poly (A) tails. However, with the recent developments in highthroughput sequencing technology, circRNAs have been found to be closely related to many human diseases including cancer, cardiovascular diseases, autoimmune diseases, Alzheimer's disease, and so on (Meng et al., 2017;Altesha et al., 2019;Zhou et al., 2019;Huang J.-L. et al., 2020).
A large body of evidence have showed that circRNAs function through a variety of ways, such as by regulating transcription, splicing and chromatin interactions, acting as microRNA (miRNA) decoys, serving as scaffolds for circRNA-protein complexes, being translated to proteins, and so on (Chen, 2020;Shi et al., 2020;Zhou et al., 2020). Emerging evidence revealed that abundant conserved miRNA response elements (MREs) exist in circRNAs so that the circRNAs can function as competing endogenous RNAs (ceRNAs) in the regulation of human cancer progression (Cui et al., 2018;Zhong et al., 2018). For example, circRNF20 harbours miR-487a to regulate the expression of HIF-1α and HK2, thereby promoting breast cancer tumorigenesis (Cao et al., 2020). CircNRIP1 can promote gastric cancer and cervical cancer progression by sponging miR-149-5p (Zhang et al., 2019) and miR-629-3p , respectively. Research have also shown that circHIPK3 regulates cancer cell growth by sponging multiple miRNAs, such as miR-7 (Zeng et al., 2018), miR-485-3p (Lai et al., 2020), miR-107 (Wei J. et al., 2020), and so on. Furthermore, numerous circRNAs have been demonstrated to play crucial roles in HCC progression through acting as ceRNAs (Han et al., 2020;Xu G. et al., 2020). However, the circRNA-related ceRNA network in HCC have not been completely elucidated and further studies are needed.
In the present study, we mined the pivotal circRNAs and constructed the circRNA-related ceRNA regulatory network in HCC using the Gene Expression Omnibus (GEO) database, the Cancer Genome Atlas (TCGA) database, and through bioinformatics methods. The flow chart of our study is exhibited in Figure 1. Firstly, robust rank aggregation (RRA) (Kolde et al., 2012) and weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008) algorithms were applied to find out the HCC-related circRNAs based on GEO datasets. Then, the expression levels of selected circRNAs were verified by quantitative real-time polymerase chain reaction (qRT-PCR). Subsequently, computational biology was utilized to predict the potential circRNA-miRNA-mRNA interactions. A protein-protein interaction (PPI) network was constructed to identify the hub genes of the ceRNA network. Finally, functional enrichment analysis, survival analysis, and immune cell infiltration analysis of hub genes were performed to thoroughly figure out the potential mechanisms of the ceRNA network.

Human Tissue Specimens
In our study, 30 pairs of HCC and their corresponding paracancerous tissues were collected from the Affiliated Cancer Hospital of Guangxi Medical University. All patients were diagnosed with HCC based on pathological examination and did not receive any preoperative anticancer treatment. The tissues were snap frozen in liquid nitrogen after isolation then stored in liquid nitrogen until usage. All patients signed informed consent forms and our study was permitted by the Ethics Committee of Guangxi Medical University.

Expression Data Collected From the GEO and TCGA Database
Microarray datasets of circRNA expression profile in HCC tissues were retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/ geo/) using the following keywords: (circular RNA OR circRNA) AND (hepatocellular carcinoma OR HCC OR liver cancer). The datasets selected met the following filter criteria of 1) The circRNA expression profile was derived from cancer tissues and normal tissues in HCC patients; 2) The sample size was at least five for each group. Eventually, four microarray datasets (GSE94508, GSE97332, GSE155949, and GSE164803) were enrolled and further analyzed.
The RNA sequencing (RNA-seq) data of HCC tissues and normal tissues was extracted from the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) project (https:// cancergenome.nih.gov/). The expression data were obtained and converted into transcripts per million (TPM), and then log2 transformation was used for normalization. Finally, mRNA expression profile of 371 HCC tissues and 50 adjacent normal tissues were utilized for differential expression analysis.

Differential Expression and RRA Analysis
First, the microarray dataset (GSE94508) with raw counts were normalized by the "edgeR" package (version: 3.34.0). Then, the "limma" package (version: 3.48.1) was implemented to screen the differentially expressed circRNAs (DEcircRNAs) and mRNAs (DEmRNAs) of the four GEO datasets and TCGA database, respectively. DEcircRNAs and DEmRNAs were defined as those with a p-value < 0.05 and the |log2 (fold change)| > 1.
Volcano plots of the differential expression results were generated using the "ggplot2" package (version: 3.3.5). In order to integrate and rank all the DEcircRNAs of these four circRNA microarray datasets, the "RobustRankAggreg" package (version: 1.1) was conducted for RRA analysis (Kolde et al., 2012). The heatmap for the result of RRA analysis was created by the "pheatmap" package (version: 1.0.12) in R.

WGCNA Analysis
Considering that the GSE155949 dataset has the largest sample size (49 pairs of HCC tissues and matched adjacent normal tissues) and the largest number of detected circRNAs (n 10592) simultaneously, GSE155949 was selected for WGCNA analysis. The "WGCNA" package (version: 1.70-3) was used to investigate the relationships between gene modules and clinical traits (Langfelder and Horvath, 2008). In the present study, the optimum soft threshold power was six when we set the scale free R 2 > 0.8. After the adjacency matrix was converted into a Topological Overlap Matrix (TOM), the circRNAs were clustered into different module eigengenes (MEs) when the minimal module size was set to 30. Then, the correlations between MEs and HCC were calculated using the Pearson correlation analysis and the module with the highest correlation coefficient was selected for further analysis.

qRT-PCR, TA Cloning and Sequencing
Total RNA was extracted from tissues using the Trizol Reagent (Invitrogent, USA) and reverse transcribed to complementary DNA using the PrimeScript ™ RT reagent Kit with gDNA Eraser (Takara, Japan) following the manufacturer's instructions. The relative RNA expression level was evaluated by the TB GreenTM Premix Ex TaqTM II Kit (Takara, Japan) in the real-time PCR system (Applied Biosystems StepOnePlus, USA). GAPDH was employed as the internal control. The primer sequences are listed in Table 1.
To verify whether the PCR products of the four circRNAs contain their back-splice sites, the amplified products were cloned into a TA vector (Zero TOPO-TA Cloning Kit, Sangon Biotech) and sequenced by Sangon Biotech (Shanghai, China) Co., Ltd. The relative expression levels of the four circRNAs were calculated using the 2 −△△Ct method.

RNase R Treatment
Total RNA (5 µg) was incubated at 37°C for 15 min with or without 4U/µg of RNase R (Epicentre Biotechnologies, USA), followed by RNase R inactivation at 70°C for 15min. After that, reverse transcription and qRT-PCR were conducted to detect the expression levels of the four circRNAs and their corresponding mRNAs. GAPDH expression in the RNase R (-) group was used as reference, and the expression levels of target RNAs were presented as mean ± standard deviation.

Establishment of the PPI Network and Identification of Hub Genes
To elucidate the interaction among the candidate interest genes, a PPI network was constructed by the STRING (https://www. string-db.org/) (Szklarczyk et al., 2019) online tool and visualized by the Cytoscape software. The size of the network nodes indicated the concatenation degree: the larger the nodes, the higher the connection degree. The edges' size demonstrated the combined score of two interacting genes, with a thicker line representing a greater combined score. To identify the hub genes  GAPDH  AGCCACATCGCTCAGACAC  GCCCAATACGACCAAATCC  hsa_circRNA_091581  GACCACCACTAGGCCTTTGAA  TGGAGTCAGGCTTGGGTAGT  hsa_circRNA_066568  TGACTGTTCAAGATACGAGGAGG  CTCTCACGTTCCCCAACCAT  hsa_circRNA_105031  GATGCACGGTGCTACACCTA  GGAAAAGAACCCTGATTGCCTG  hsa_circRNA_104515  TGAACCTGACACCATCAGCAA  TGAGTCTGAGGCCACTCCTT  GPC3  CAGTCAGCAGGCAACTCCGAAG  AAGAAGAAGCACACCACCGAGATG  ROBO1  TGTTGCTTTGGGACGGACTGTAAC  GGCTGGATGACTGTGGTGGTTG  MBNL3  CTGATTCCTGGAAACCCACCTCTTG  CACGCTGAAATTCTCGGCAAACC  AGK  AGTGTCTCCAAGCCAGCCAGTG  CTCCACAGGCATCGCTTCATACTC qRT-PCR, quantitative real-time polymerase chain reaction.
Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 9 | Article 801478 of the PPI network, the degree of mRNAs was calculated by the Cytoscape plugin "cytoHubba", and mRNAs with a degree ≥10 were regarded as hub genes.

Survival Analysis of Hub Genes
To further explore the influence of hub genes on the prognosis, GEPIA2 (http://gepia2.cancer-pku.cn/#index)  was implemented to determine the relationships between hub genes and the overall survival (OS) and disease-free survival (DFS) of HCC patients. In the present study, Kaplan-Meier plots were drawn and hub genes with log-rank p < 0.05 were considered statistically significant. The hazard ratio (HR) of high expression group compared with low expression group could also be obtained from the survival analysis.
Association Between Core Genes and Infiltrating Immune Cells in HCC TIMER2.0 (http://timer.cistrome.org/) (Li T. et al., 2020) is a comprehensive resource for systematic analysis of immune infiltrates across diverse cancer types. We used the TIMER web tool to investigate the correlations between the hub genes and infiltrating immune cells, including B cells, CD8 + T cells, CD4 + T cells, neutrophils, macrophages, and myeloid dendritic cells. p-value less than 0.05 was considered statistically significant.

Identification of 22 DEcircRNAs in HCC Based on the RRA Method
The detailed information of four included circRNA microarray datasets (GSE94508, GSE97332, GSE155949, and GSE164803) is listed in Table 2. In brief, 523 DEcircRNAs were identified in the GSE94508 dataset, including 96 upregulated circRNAs and 427 downregulated circRNAs (Figure 2A). A total of 892 DEcircRNAs with 453 upregulated circRNAs and 439 downregulated circRNAs were found in the GSE97332 dataset ( Figure 2B). Meanwhile, 50 DEcircRNAs consisting of 30 upregulated circRNAs and 20 downregulated circRNAs were determined in the GSE155949 dataset ( Figure 2C). As for the GSE164803 dataset, 514 DEcircRNAs, of which 193 were upregulated and 321 were downregulated, were recognized ( Figure 2D). After that, we integrated and ranked the DEcircRNAs from the four microarray chips using the RRA algorithm. Finally, 22 DEcircRNAs, including nine upregulated circRNAs and thirteen downregulated circRNAs, were identified to be statistically significant with a p-value < 0.05 ( Figure 3).

Identification of Key circRNA Co-Expression Module Based on WGCNA Analysis
To identify the key module of circRNA in HCC, WGCNA analysis was conducted in the GSE155949 dataset. When setting the scale-free R 2 > 0.8 and mean connectivity of all circRNAs <100, six was chosen as the best soft threshold (power) in our study ( Figure 4A). A total of 17 coexpression modules were classified using the cluster analysis and the cluster dendrogram is presented in Figure 4B. Subsequently, the relationships between the modules and HCC were calculated and is shown in Figure 4C. The pink module containing 404 circRNAs had the strongest positive correlation (r 0.67, P 7e-14) with the HCC tissues. Besides, the scatter plot of the module membership vs. gene significance in the pink module also showed a significant correlation (r 0.43, p 1.3e-19; Figure 4D). Then, these 404 circRNAs of the pink module overlapped with 22 DEcircRNAs identified by the RRA method, and only four circRNAs (hsa_circRNA_091581, hsa_circRNA_066568, hsa_circRNA_105031, and hsa_circRNA_104515) were recognized in both algorithms ( Figure 4E).

Validation of the Expression of Four circRNAs Using qRT-PCR
The detailed information of the four circRNAs is listed in Table 3 and visualized in Figure 5. Furthermore, the result of TA cloning and sequencing confirmed the back-splicing junction sites of the PCR products of the four circRNAs ( Figure 5). To verify the resistance of circRNAs to RNase R digestion, qRT-PCR was performed to evaluate the relative expression of the four circRNAs and their host gene mRNAs under the RNase R treatment or not. As shown in Figure 6, the expression levels of the four mRNAs were significantly reduced in the RNase R (+) group compared to the RNase R (-) group. However, the expression levels of the four circRNAs did not show any obviously downregulated tendency under the RNase R treatment, implying the existence of these four circRNAs. Eventually, the expression levels of the four circRNAs in HCC were examined in 30 pairs of HCC tissues and their corresponding adjacent normal tissues. The results indicated that circRNAs, hsa_circRNA_091581 (p < 0.001; Figure 7A), and hsa_circRNA_066568 (p 0.001; Figure 7B) were upregulated and hsa_circRNA_105031 was downregulated in HCC tissues (p 0.007; Figure 7C), which showed statistical significance. However, there was no significant difference in the expression level of hsa_circRNA_104515 between HCC tissues and paracancer tissues (p 0.159; Figure 7D).

Construction of the Preliminary circRNA-miRNA-mRNA Network
Based on the qRT-PCR results, differentially expressed circRNAs, hsa_circRNA_091581, hsa_circRNA_066568, and hsa_circRNA_105031, were chosen for further analysis. Firstly, we obtained and overlapped the potential target miRNAs of the three circRNAs from two online tools, circBank and CircInteractome. As a result, 17 circRNA-miRNA pairs consisting of three circRNAs and 17 miRNAs were recognized. KEGG analysis of the 17 miRNAs in DIANA-miRPath suggested that these miRNAs were all significantly associated with some cancer-related pathways (Supplementary Figure S1). Next, a total of 1799 predicted mRNAs of these 17 miRNAs were obtained from four miRNA-prediction databases (TargetScan, miRDB, miRTarBase, and miRWalk 3.0). Besides, 2863 DEmRNAs including 2396 upregulated mRNAs and 467 downregulated mRNAs in HCC were found based on the TCGA database (Supplementary Figure S2A). By overlapping the predicted target mRNAs with the DEmRNAs, 222 mRNAs that may exert pivotal effect on HCC were selected for further study (Supplementary Figure S2B). Finally, a preliminary circRNA-miRNA-mRNA regulatory network, which contains three

Identification of Seven Core Genes Based on the Connection Degree of PPI Network
A PPI network was established to determine the interactions among 222 mRNAs of interest. After removing disconnected nodes, 169 nodes and 269 edges remained in the PPI network and was visualized by the Cytoscape software as shown in Figure 9A. To find out the hub genes involved in the PPI network, the connection degree of each mRNA was calculated by the cytoHubba plugin in Cytoscape. As shown in Figure 9B, seven genes with degree values ≥10 were identified as core genes of the PPI network: ESR1 (estrogen receptor 1), BUB1 (BUB1 mitotic checkpoint serine/threonine kinase), PRC1 (protein regulator of cytokinesis 1), LOX (lysyl oxidase), CCT5 (chaperonin containing TCP1 subunit 5), YWHAZ (tyrosine 3-monooxygenase/tryptophan 5monooxygenase activation protein zeta), and DDX39B (DExD-box helicase 39B). Then, a circRNA-miRNA-hubgene regulatory network was reconstructed and visualized in Figure 9C. In summary, ten ceRNA regulatory axes, consisting of two circRNAs (hsa_circRNA_105031 and hsa_circRNA_091581), six miRNAs (has-miR-1827, has-miR-548k, has-miR-671-5p, has-miR-767-5p, has-miR-595, and has-miR-661), and seven hub genes (ESR1, BUB1, PRC1, LOX, CCT5, YWHAZ, and DDX39B) were recognized.

GO Annotation, KEGG, and Reactome Pathway Analyses of the Seven Core Genes
To better understand the potential molecular function of the seven hub genes, enrichment analyses were carried out in the KOBAS-i database and the results are exhibited in Table 4. GO enrichment analysis indicated that the seven hub genes mainly participated in the processes of "protein kinase binding", "lung development", "nucleoplasm", and "protein binding" (corrected p < 0.001). KEGG pathway enrichment analysis revealed that the seven hub genes were significantly involved in the "Cell cycle" and "Oocyte meiosis" pathways (corrected p < 0.001). As for the Reactome analysis, the seven core genes were enriched in the pathways of "RHO GTPase Effectors" and "Signaling by Rho GTPases" (corrected p < 0.001).

Survival Analysis of the Seven Hub Genes in HCC
The GEPIA2 was utilized to explore the prognostic significance of the seven hub genes in HCC patients and the results are shown in Figure 10. Notably, all the seven core genes were related to the OS of HCC patients based on TCGA database. The high expression of ESR1 (HR 0.55) was significantly associated with longer OS of HCC patients (log-rank p < 0.05). On the contrary, the high expression of BUB1 (HR 1.8), PRC1 (HR 1.9), LOX (HR 1.5), CCT5 (HR 1.6), YWHAZ (HR 1.5), and DDX39B (HR 1.6) were remarkably correlated with worsened OS of HCC patients (log-rank p < 0.05). As for the DFS analysis, the increased expression of ESR1 (HR 0.74) was significantly associated with better DFS of HCC patients. The high expression of BUB1 (HR 1.6), PRC1 (HR 1.7), CCT5 (HR 1.5), and DDX39B (HR 1.4) were remarkably correlated with poorer DFS of HCC patients (logrank p < 0.05). However, there was no significant relationship between the expression levels of two core genes, LOX (HR 1.3) and YWHAZ (HR 1.3), and the DFS time of HCC patients (log-rank p > 0.05).

Correlation Between the Hub Genes and Immune Cells in HCC
To further investigate whether the hub genes were related to the infiltrating levels of immune cells in HCC, the TIMER database was applied in our study. As shown in Supplementary Figure S3, the expression levels of BUB1, PRC1, and LOX was significantly associated with the infiltrating levels of six immune cells, including B cells,    Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 9 | Article 801478 9 CD8 + T cells, CD4 + T cells, neutrophils, macrophages, and myeloid dendritic cells in HCC. The CCT5 expression level was closely related to the infiltrating levels of all six immune cells except for CD8 + T cells in HCC. The ESR1 expression level was significantly associated with the infiltrating levels of only B cells, CD8 + T cells, CD4 + T cells, and myeloid dendritic cells in HCC. For YWHAZ, the significant positive correlation was observed in B cells, neutrophils, macrophages, and myeloid dendritic cells infiltration in HCC. However, the relationships between DDX39B and the infiltrating levels of immune cells in HCC was not explored because DDX39B was not available in the TIMER database.

DISCUSSION
Tumorigenesis is a complex and multi-stage process involving changes in the expression levels of genes and alterations of the microenvironment. CircRNA, a novel type of non-coding RNA, has increasingly become a hot spot of cancer research due to its special covalently closed-loop structure. In this study, four DEcircRNAs in HCC were identified by RRA rank analysis of four GEO datasets (GSE94508, GSE97332, GSE155949, and GSE164803) and WGCNA analysis based on the GSE155949 dataset. Subsequently, the expression levels of these four circRNAs were verified by qRT-PCR and three DEcircRNAs (hsa_circRNA_091581, hsa_circRNA_066568, and hsa_circRNA_105031) were chosen for further analysis. To determine whether the above three circRNAs function as ceRNAs in HCC, a circRNA-miRNA-mRNA regulatory network was constructed based on the predicted targets from online tools and the DEmRNAs from the TCGA database. Seven hub genes (ESR1, BUB1, PRC1, LOX, CCT5, YWHAZ, and DDX39B) were recognized in the PPI network and functional enrichment analysis suggested that the core genes were mainly involved in "protein kinase binding", "Cell cycle", and "RHO GTPase Effectors" pathways. More importantly, all these seven hub genes were significantly associated with the prognosis of HCC patients and the infiltrating levels of immune cells in HCC. Collectively, the above results indicated that the ceRNA network constructed in this study may play important roles in the tumorigenesis of HCC. Among the three DEcircRNAs validated by qRT-PCR, hsa_circRNA_091581 was the most studied. Wei et al. found that hsa_circRNA_091581 had an increased expression in HCC tissues, and its high expression was remarkably associated with poor OS and DFS of HCC patients. Mechanism research showed that hsa_circ_091581 can promote HCC cell proliferation by acting as a sponge of miR-526b to block the degradation of c-MYC mRNA . Ji et al. also revealed that hsa_circ_091581 promoted HCC progression by targeting the miR-591/FOSL2 axis (Ji et al., 2021). Besides, as the most upregulated circRNA in the GSE155949 dataset, hsa_circRNA_091581 was experimentally demonstrated to promote the tumorigenesis of HCC by targeting miR-378a-3p to accelerate the cell cycle progression (Han et al., 2021). On the contrary, hsa_circRNA_105031 was found to be downregulated in HCC tissues and suppressed HCC progression by sponging miR-1307 to regulate the expression of ISM1 . The above research indicated that the two circRNAs played vital roles in the development and progression of HCC. However, we have not found any relevant studies on hsa_circRNA_066568 yet. Therefore, more studies on the mechanisms of these three circRNAs in HCC are urgently needed, especially for hsa_circRNA_066568.
To further investigate the ceRNA-related mechanism of the three circRNAs in HCC, a circRNA-miRNA-mRNA regulatory network was constructed based on bioinformatics analysis. As a result, three circRNAs, 17 miRNAs, and 222 mRNAs were involved in the ceRNA network. KEGG enrichment analysis of the 17 miRNAs revealed that these miRNAs were involved in multiple pathways related to HCC, such as the Hippo (Liu et al., 2020b), FoxO (Hou et al., 2016), TGF-beta (Mazzocca et al., 2012), PI3K-Akt (Zhou, 2011), and Wnt (Pez et al., 2013) signaling pathways. A PPI network of 222 mRNAs identified seven core genes (ESR1, BUB1, PRC1, LOX, CCT5, YWHAZ, and DDX39B), and the enriched GO terms for these hub genes were closely corelated to carcinogenesis processes such as protein kinase binding and protein binding. Cell cycle FIGURE 8 | Preliminary circRNA-miRNA-mRNA regulatory network in HCC. The network, consisting of three circRNAs, 17 miRNAs, and 222 mRNAs, was visualized by the Cytoscape software (3.8.3). HCC, hepatocellular carcinoma. pathway, the well-known relevant pathway to the aberrant cell proliferation of cancers (Williams and Stoeber, 2012), was the most significant KEGG enrichment pathway of the seven hub genes. For the Reactome pathway analysis, "RHO GTPase Effectors" and "Signaling by Rho GTPases" pathways were also enriched. Previous studies revealed that the RhoGTPases/Rho-effector signaling cascade played important roles in mediating HCC metastasis (Wong et al., 2010). Taken together, the circRNA-miRNA-mRNA regulatory network that we constructed might help to elucidate the regulatory mechanisms of circRNA-related ceRNA in HCC.
In this study, we also constructed the circRNA-miRNA-hubgene network. As a consequence, ten ceRNA regulatory axes (hsa_circRNA_091581/miR-661/CCT5, hsa_circRNA_105031/ miR-1827/ESR1, hsa_circRNA_105031/miR-1827/YWHAZ, hsa_circRNA_105031/miR-548k/YWHAZ, hsa_circRNA_105031/ miR-671-5p/PRC1, hsa_circRNA_105031/miR-671-5p/YWHAZ, hsa_circRNA_105031/miR-671-5p/DDX39B, hsa_circRNA_105031/miR-595/BUB1, hsa_circRNA_105031/miR-595/LOX, and hsa_circRNA_105031/miR-767-5p/LOX), were determined. In addition, the survival analysis based on the TCGA database showed that the expression levels of these seven hub genes were significantly associated with the prognosis of HCC patients,  suggesting that the seven hub genes in our network might be promising prognostic indicators for HCC patients. Another interesting finding in our study is that the expression levels of the seven hub genes were closely related to the infiltrating levels of immune cells. Immune cell infiltration, as an important constituent of tumour microenvironment, has been widely considered to be involved in tumour growth and metastasis (Hinshaw and Shevde, 2019). Our study showed that five upregulated hub genes, BUB1, PRC1, LOX, CCT5, and YWHAZ, had positive relationships with immune cells infiltration levels. In contrast, the downregulated gene, ESR1, mainly exhibited negative correlations with immune cells infiltration levels. Given the above evidence, the circRNA-miRNAhubgene regulatory network may play vital roles in HCC and is worthy of further study. Existing research have also shown that the seven core genes played important roles in the development and progression of HCC. ESR1 was identified as a tumour suppressor gene in HCC and its genetic polymorphism was significantly associated with susceptibility Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 9 | Article 801478 13 to HCC in Chinese hepatitis B virus carriers (Zhai et al., 2006;Hishida et al., 2013). BUB1, as a cell cycle related gene, might contribute to the occurrence and development of HCC when there is aberrant increased expression (Yan et al., 2017;Zhang et al., 2020). Chen et al. reported that PRC1 promoted early HCC recurrence by enhancing the Wnt/β-catenin signaling pathway (Chen et al., 2016). Besides, high PRC1 expression enhanced HCC cells' resistance to chemotherapy and was associated with unfavourable survival for HCC patients (Y. Wang et al., 2017). Evidence showed that increased LOX could promote angiogenesis, induce epithelial-mesenchymal transition, and is related to early recurrence of HCC (Umezaki et al., 2019;Yang et al., 2019). Furthermore, accumulated evidence has demonstrated that multiple molecules mediated the tumorigenesis and metastasis of HCC by regulating the YWHAZ-associated axis (Zhao et al., 2018;Wei et al., 2019;Liu et al., 2020a;Shen et al., 2020). Bioinformatic analysis also revealed the miR-139-5p/CCT5 axis might exert a momentous effect on HCC progression, but in-depth validated studies have not been carried out (Xu et al., 2021). As for DDX39B, its functions in HCC progression have not been studied yet. However, recent studies have revealed that increased DDX39B could promote the tumour metastasis of colorectal cancer and enhance chemotherapy resistance in BRCA1-proficient ovarian cancers He et al., 2021). In our study, we also found that DDX39B was upregulated in HCC tissues and high DDX39B expression was correlated to poor OS and DFS in HCC patients based on the TCGA database. These findings suggest that DDX39B may function as a tumour promoter in HCC.
Previous studies have also explored the potential circRNA-miRNA-mRNA regulatory networks in HCC. For example, Yang et al. screened out DEcircRNAs based on the GSE97332 dataset and selected the top 10 upregulated and downregulated circRNAs to construct the circRNA-miRNA-mRNA network (Yang et al., 2020). Another study obtained 26 overlapping DEcircRNAs from the GSE97332 and GSE94508 datasets. After being verified to be upregulated in the SMMC-7721 cell line, six DEcircRNAs were selected to construct a circRNA-miRNA-mRNA network . Besides, Xiong et al. identified six DEcircRNAs by integrating three microarray chips (GSE78520, GSE94508, and GSE97332) using the RRA algorithm. The qRT-PCR results of 16 pairs of HCC tissues and adjacent normal tissues revealed that only has_circRNA_102166 was downregulated in HCC tissues, with a p-value less than 0.05 . Compared with existing related studies, our research has some advantages. Firstly, we excluded microarray chips with sample sizes less than five pairs and used the RRA method to integrate the four included datasets (N 67 pairs) to eliminate the effect of small sample size. Secondly, we comprehensively utilized the RRA algorithm, WGCNA analysis, and qRT-PCR validation to screen out critical circRNAs that might play important roles in the progression of HCC. Last but not the least, integrated analyses including PPI network, functional enrichment analysis, survival analysis, and infiltrating immune cells analysis were performed to uncover the potential regulatory mechanisms of HCC progression. However, given that our circRNA-miRNA-mRNA network is constructed based on bioinformatic analysis, further in-depth studies are urgently needed. For example, dualluciferase reporter and RNA immune co-precipitation (RIP) assays can be performed to verify the interaction between miRNA and circRNA, and between miRNA and mRNA in these predicted regulatory axes. At the same time, more experiments are needed for hsa_circRNA_066568 that has not been reported. Cell counting kit 8, colony formation assay, Transwell migration and invasion assays, cell cycle assay, apoptotic assay, and xenograft tumor assay should be conducted to evaluate the effects of hsa_circRNA_066568 on HCC cells both in vitro and in vivo. Besides, mechanism researches are also needed to reveal the regulatory axes of hsa_circRNA_066568 in the progression of HCC.

CONCLUSION
In summary, our research identified three HCC-related circRNAs (hsa_circRNA_091581, hsa_circRNA_066568, and hsa_circRNA_105031) based on integrated analysis by using the RRA algorithm, WGCNA analysis, and qRT-PCR validation. A circRNA-miRNA-mRNA regulatory network was constructed by computational biology and seven hub genes were identified through the PPI analysis. We also found that the seven core genes might influence the prognosis of HCC patients through regulating the cell cycle and other pathways. Moreover, these hub genes were significantly correlated with the infiltrating levels of immune cells in HCC. Our study provides new insights into the underlying mechanisms of the circRNA-miRNA-mRNA regulatory network in HCC.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: [http://www.ncbi.nlm.nih.gov/geo// GSE94508, GSE97332, GSE155949, and GSE164803] and [https://cancergenome.nih.gov//LIHC project]. The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

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

AUTHOR CONTRIBUTIONS
MM designed the research, performed data analysis and experiments, and wrote the manuscript. BL screened for the eligible datasets, assisted with data analysis, and performed the experiments. YL screened for the eligible datasets, assisted with data analysis and manuscript preparation. JT assisted with manuscript preparation and revised the manuscript for language errors. XZ helped to perform the experiments. XYZ and DH assisted with clinical sample collection. CL assisted with data analysis. SL and XQ designed the research and provided final approval of the manuscript. All the authors read and approved the final manuscript.