Construction and Investigation of circRNA-miRNA-mRNA ceRNA Regulatory Network in Breast Cancer

Background: Circular RNAs (circRNAs) have drawn lots of attention in tumorigenesis and progression. However, circRNAs as crucial regulators in multitudinous biological processes have not been systematically identied in breast cancer (BC). Our research aims to explore novel circRNAs in BC and their mechanisms of action. Methods: The circRNA expression prole data, as well as RNA-sequencing data of BC, were downloaded from public database, respectively. The differentially expressed circRNAs, miRNA, and mRNA were determined via fold change ltering. The competing endogenous RNAs (ceRNAs) network were established on the foundation of the relationship between circular RNAs, miRNAs and mRNAs. GO and KEGG analysis of the overlapped genes were performed to predict the potential functions and mechanisms of circRNAs in BC. The CytoHubba was used to determine the hub genes from the PPI regulatory network. Morever, we further used Kaplan–Meier plotter to perform survival analysis of these hub genes. Real-time PCR was used to validate the expression of the circRNAs in BC tissues. Results: A total of seven differential expressed circRNAs were screened. After the predicted target miRNA and DEmiRNA were intersected, four circRNA-miRNA interactions including three circRNAs and four miRNAs were determined. Furthermore, the Venn diagram was used to intersect the predicted target genes and the downregulated differentially expressed genes, and screened 149 overlapped genes. Moreover, we constructed a PPI network, and selecting six hub genes, including DGAT2, ACSL1, ADIPOQ, LPL, LEP, PCK1. Moreover, the survival analysis results revealed that low expression of ADIPOQ, LPL, LEP were obviously correlated with poor prognosis of BC patients. The real-time PCR results demonstrated that, the levels of circ_0028899, circ_0000375, and circ_0000376 were signicantly down-regulated in breast cancer tissues. Conclusions: Our study constructed and analyzed a circRNA-associated ceRNA regulatory network and discovered that circ_0028899, circ_0000375, and


Introduction
Breast cancer (BC) is the most common occurring neoplasm and the second leading cause of cancer mortality in females. The incidence of BC has been rising over the past decade, far exceeding the rate of increase in other cancers. Despite progress in treatment strategies, some BC patients still have poor prognosis, especially for metastatic BC patients [1]. Therefore, identifying novel biomarkers and elucidating molecular mechanisms are warranted to be discovered for the early detection, direct treatment, and judge prognosis for BC. Circular RNAs (circRNAs) are a subclass of non-coding RNAs that forms a closed loop structure by a process called "back splicing" [2,3]. Due to the lack of 3' tails and 5' caps structure, thus circRNAs appear to be more stability than their linear types [4]. Accumulating evidence suggests that circRNAs are aberrantly expressed in many different kinds of malignant tumors, and play crucial regulatory roles in tumorigenesis and development [5,6]. In molecular biology, competing endogenous RNAs (ceRNAs) modulate other RNA transcripts via competing binding to shared miRNAs. Some evidences have indicated that many circRNAs are participated in the generation of tumors through the ceRNA mechanism [7]. For instance, ciRs7 can abrogate the function of miR-7, resulting in the reduced miR-7 activity and relieved the latter's inhibition of its target gene [8][9][10]. CircHIPK3 can adsorb multiple miRNAs to regulate cancer cell growth, and these miRNAs have been reported as tumor suppressor miRNAs [11]. Analogously, circTP63 is proven to interact with FOXM1 by abrogating the function of miR-873-3p in lung squamous cell carcinoma [6]. CircRNAs as ceRNAs that mediate malignant processes have also been con rmed in BC [12,13]. However, there are still many unknown circRNAs to be further studied.
In our study, we collected circRNAs microarray data and RNA-seq data from the GEO and TCGA database, respectively. We used the limma package to identify differentially expressed circRNAs. Then, we further predicted the differentially expressed circRNAs sponge miRNAs and miRNAs target genes, and constructed a circRNA-miRNA-mRNA regulatory network. GO and KEGG on the differentially expressed mRNAs of the ceRNA regulation network were performed to assess the potential molecular mechanisms participated in tumorigenesis. Subsequently, we constructed a PPI network and picked out six hub genes. We further constructed the circRNA-miRNA-hub genes subnetwork regulation module to better study the pathogenesis of BC. Furthermore, the prognosis of patients with BC was further evaluated in the Kaplan-Meier plotter databases. The current research was conducted to shed light on underlying mechanisms of BC pathogenesis and provide novel targets for BC.

Materials And Methods
CircRNA Microarry and RNA-sequencing data Collection The circRNA expression pro le data for GSE101123 was downloaded from the GEO database and included 8 BC sample and 3 mammary gland sample. The RNA-sequencing data were obtained from TCGA database.
The miRNA expression pro le data contained 1103 BC samples and 104 normal samples, and the mRNA expression pro le data included 1109 BC samples and 113 normal samples.

Differential expression analysis of circRNAs, miRNAs and mRNAs
Using P-adjust < 0.05 and |log 2 FC|≥ 1 as a threshold, we performed differential expressed circRNAs analysis in GSE101123 by using limma package in R. Then, we screened differential expressed miRNAs by using edgeR package (FDR< 0.05 and |log 2 FC|≥ 1). In addition, the Deseq package were used to select differential expressed mRNAs with criteria of FDR < 0.05 and |log2FC|≥ 2.5.
Prediction and screening of miRNA binding sites Prediction of interactions between circRNAs and miRNAs using a cancer-speci c CircRNA (CSCD) database.
We then used the Venn diagram to intersect the target miRNAs and the differentially expressed miRNAs (DEmiRNAs) from TCGA database.

Prediction and screening of miRNA-Targeted mRNAs
The miRNA-targeted mRNAs were predicted based on the miRWalk 3.0 software. Furthermore, we intersected the target genes predicted by miRWalk software and the differential genes in the TCGA database.

Construction of ceRNAs regulatory network
The ceRNA network was established on the foundation of the interactions between differentially expressed circRNAs and DEmiRNAs and between DEmiRNAs and differentially expressed mRNAs. CytoCope 3.7.1 software was used to visualize the established network.

GO and KEGG enrichment analyses
GO enrichment analysis of the overlapped genes was implemented by clusterPro ler R package. GO terms (biological processes, molecular function, and cellular components) with P-adjust< 0.05 were considered to be signi cantly enriched by the overlapped genes. KEGG pathway analysis of the overlapped genes was accomplished by using DAVID bioinformatics tool (https://david.ncifcrf.gov/). KEGG pathways analysis with P < 0.05 were evaluated for signi cance.

Formation of PPI regulatory network and recognition of hub genes
To further obtain the interactions between the overlapped genes, we constructed the PPI network by using the STRING (V 11.0, https://string-db.org/), and the con dence score≥0.4 was selected. Additionally, screen the hub genes from the PPI network by applying the CytoHubba plugin in cytoscape. The top 6 genes were screened as hub genes from the PPI network by using MCC method.

Patients and tissue specimens
Thirty paired fresh frozen BC tissues and adjacent tissues were obtained from patients with breast cancer who underwent surgery at the Department of Breast center of the Fourth Hospital of Hebei Medical University between January 2017 and February 2018. None of these patients received any radiotherapy or chemotherapy prior to surgical resection. This study was consented by the Ethics Committee of the Fourth Hospital of Hebei Medical University.

Real-time PCR validation
Total RNA from each tissue was extracted using Trizol reagent (Invitrogen, USA) and converted to complementary DNA (cDNA) using the GoScriptTM Reverse Transcription System (Promega, USA). Real-time PCR was performed using the SYBR Green Real-Time PCR Reagent Kit (Promega, USA). The expression level of circRNA was detected by using speci c primers (RiboBio, China) and normalized against U6 small nucleolar RNA.

Survival Analysis for hub genes
To assess the prognostic signi cance of these hub genes on survival, the Kaplan-Meier plotter database was applied to evaluate the correlation between six hub genes and overall survival. Dividing patients into two groups according to the median expression of the gene, and p less than 0.05 was considered statistically signi cant. Real-time PCR data was presented as mean ±SD and analyzed by Studen's t test, and p less than 0.05 was considered statistically signi cant.

Identi cation of Differentially Expressed circRNAs in BC
Raw data from GSE101123 were downloaded from GEO database and then performed differential gene expression analysis by limma package in R. Using P-adjust < 0.05 and |log 2 FC|≥ 1 as a threshold, 7 differential expressed circRNAs were screened out, of which, 4 were upregulated and 3 were downregulated. The differential expressed circRNAs were displayed in the volcano plot (Fig. 1a), and also exhibited in heat map (Fig. 1b). The essential characteristics and basic structural of the seven circRNAs were displayed in Table 1 and Fig. 2.

Prediction Of Circrna-bound Targeted Mirna
Since many circRNAs included the binding site of miRNAs, it could exert sponge activity by binding to a speci c miRNA, and therefore could suppress miRNA function. We rst predicted circRNA-bound miRNAs using CSCD miRNA target online prediction software, and total 441 target miRNAs of the seven circRNAs were predicted. Furthermore, we downloaded raw sequencing data from the o cial TCGA website. Using FDR < 0.05 and |log 2 FC|≥ 1 as the threshold, we identi ed a total of 251 differentially expressed miRNAs, among them 64 were downregulated and 187 were upregulated (Fig. 3a). We further used the Venn diagram to intersect the target miRNAs and DEmiRNAs, and screened 5 target miRNAs (miR-934, miR-7760, miR-760, miR-4443, miR-4784) that play important roles in BC (Fig. 3b). We then analyzed the correlation between the expression of ve miRNAs and the differential expressed circRNAs. Finally, four circRNA-miRNA interactions including 3 circRNAs (circ_0028899, circ_0000375, circ_0000376) and 4 miRNAs (miR-934, miR-7760, miR-4443, miR-4784) were determined.

Prediction Of Mirna-targeted Mrnas
We mapped the abovementioned 4 DEmiRNAs into the miRwalk target gene prediction software and searched for the target gene. Total 7168 target genes of the four miRNAs were predicted. In addition, 553 downregulated differentially expressed genes in BC were obtained from TCGA database (FDR < 0.05, |log 2 FC| ≥ 2.5) (Fig. 4a). We then used the Venn diagram to intersect the target genes and the downregulated differentially expressed genes, and screened 149 target genes (Fig. 4b).
CircRNA-miRNA-mRNA network construction In order to better comprehend the interaction of circRNA, miRNA and mRNA associated with BC, we established a circRNA-miRNA-mRNA network. As shown in Fig. 5, this network was composed of 3 circRNA nodes, 4 miRNA nodes, 149 mRNA nodes and 462 edges. This network showed the possible ceRNA relationships among the three differential expressed circRNAs (circ_0028899, circ_0000375, circ_0000376), the four miRNAs (miR-934, miR-7760, miR-4443, miR-4784), and the 149 differentially expressed genes.
GO and KEGG enrichment analysis of the overlapped genes The 149 overlapped genes were annotated, and these genes were selected for GO functional enrichment analysis. When considering Biological Process (BP), the overlapped genes were primarily enriched in response to acid chemical, regulation of membrane potential, and regulation of lipid metabolic process. In terms of cellular components (CC), the top three enriched items were plasma membrane protein complex, receptor complex, membrane raft. With regards to molecular function (MF), cation transmembrane transporter activity, peptide binding, and amide binding were enriched in the top three places (Fig. 6a). As for the pathways, the overlapped genes were enriched in several KEGG pathways, for instance PPAR signaling pathway, PI3K − Akt signaling pathway, and AMPK signaling pathway (Fig. 6b).

Formation of PPI network and hub gene selection of hub genes
We constructed the PPI network by using the string database and visualized it with Cytoscape software to explore the interactions among the 149 overlapped genes. As shown in Fig. 7a, a total of 91 nodes and 159 edges were ltered following removal of the unconnected nodes in this network. Furthermore, we screened the hub genes from the PPI regulatory network by applying the CytoHubba plugin in cytoscape. The top 6 hub genes (DGAT2, ACSL1, ADIPOQ, LPL, LEP, PCK1) were screened by using MCC method, and one subnetwork was constructed with 6 nodes and 14 edges (Fig. 7b).

CircRNA-miRNA-hub genes network construction
Then, we established a circRNA-miRNA-hub gene network to show the link between circRNAs, miRNA and hub genes. As shown in Fig. 8

Assessing the overall survival in patients with BC
We further used Kaplan-Meier plotter to perform prognostic analysis of the six hub genes, the survival analysis results showed that low expression of ADIPOQ, LPL, LEP were associated with poor overall survival in patients with BC ( Fig. 9), while the other had no statistically signi cant effect on overall survival in patients with BC.

Validation of Differentially Expressed circRNAs in BC tissue
To verify the reliability of bioinformatics analysis, the expression of the three circRNAs were detected in 30 pairs breast cancer tissues and adjacent tissues. The results indicated that circ_0028899, circ_0000375, and circ_0000376 were signi cantly down-regulated in breast cancer tissues (P < 0.05, Fig. 10), which was consistent with those of bioinformatics analysis.

Discussion
During the past few decades, a large number of circRNAs have been discovered in different kinds of diseases by analyzing high-throughput sequencing data [14]. It has been determined that circRNA has a variety of important biological functions, such as acting as a miRNA sponge element (MRE), forming a RNA-protein complex (RBP) or modulating targeted gene transcription [7,15,16]. Increasing study has indicated the vital role circRNAs play in tumor progression, and because of the relative tolerance of these circRNAs to exonuclease, it may be a better biomarker to judge the prognosis of patients [17]. Additionally, more and more circRNAs, for instance hsa_circ_0025202 [18], circAGFG1 [12], circGFRA1 [19], circ_ANKS1B [20], have been pointed out to play important roles in tumorigenesis, chemoresistance, apoptosis or metastasis of BC.
Nevertheless, many of the functions of circRNAs in BC remain to be explored.
In the present study, the circRNA expression pro le data of GSE101123 were downloaded from GEO database, and seven circRNAs (circ_0028899, circ_0000376, circ_0000519, circ_0000375, circ_0000517, circ_0003645, circ_0000516) were identi ed by limma package in R. Since many circRNAs have a large number of miRNA binding sites, it can exert sponge activity by binding to a speci c miRNA, and therefore regulate gene expression through the ceRNA mechanism [21,22]. In order to con rmed whether the abovementioned circRNAs function as ceRNAs in BC, we predicted circRNA-bound miRNAs using CSCD miRNA target online prediction software, and identi ed DEmiRNAs from TCGA databases. We further used the Venn diagram to intersect the target miRNAs and DEmiRNAs, and screened 5 target miRNAs that play important roles in BC. Subsequently, we analyzed the correlation between the expression of ve miRNAs and the differential expressed circRNAs. The expression of hsa_circ_0000519 and miR-760 was elevated in BC, which did not conform to the mechanism of ceRNA. Therefore, four circRNA-miRNA interactions including three circRNAs (circ_0028899, circ_0000375, circ_0000376) and 4 miRNAs (miR-934, miR-7760, miR-4443, miR-4784) were determined.
Furthermore, we used the Venn diagram to intersect the target genes of the four miRNAs and the downregulated differentially expressed genes in BC, and screened 149 overlapped genes that play a crucial role in the carcinogenesis and development of BC. To better understand the interaction of circRNA, miRNA and mRNA associated with BC, we established a circRNA-miRNA-mRNA network. We observed that circ_0028899, circ_0000375, circ_0000376 perhaps serving as ceRNAs to sponge miR-934, miR-7760, miR-4443, or miR-4784, and then regulate the 149 overlapped genes expression. These results provided strong evidence for the ceRNA regulation mechanisms of the three circRNAs in BC. Besides, the GO and KEGG analysis indicated that the 149 overlapped genes were participated in many key biological functions associated with tumors and metabolic pathways, for instance PPAR signaling pathway, AMPK signaling pathway, PI3K − Akt signaling pathway, Jak-STAT signaling pathway, and Ras signaling pathway. To further elucidate the functional mechanism of the ceRNA network, we established a PPI regulation network, and selecting six hub genes, containing DGAT2, ACSL1, ADIPOQ, LPL, LEP, PCK1. Previous studies have shown that ve genes (DGAT2, ACSL1, ADIPOQ, LPL, and LEP) play important roles in the carcinogenesis of BC [23][24][25][26][27]. However, the relationship between PCK1 and BC has not been investigated, and the correlation between these hub genes and circRNAs in BC has not been reported. Herein, we discovered nine circRNA-miRNA-hub genes regulatory modules, indicating competitive regulatory associations of the three circRNAs with the six genes in BC. The survival analysis results showed that low expression of ADIPOQ, LPL, LEP were associated with poor overall survival in patients with BC, while the other had no statistically signi cant effect on overall survival in patients with BC. In addition, the expression of the three circRNAs were detected in 30 pairs breast cancer tissues and adjacent tissues. The results indicated that circ_0028899, circ_0000375, and circ_0000376 were signi cantly down-regulated in breast cancer tissues, which was consistent with those of bioinformatics analysis.
Despite the ndings in clinical implications, several limitations should be noted in our study. First, the sample size of breast cancer was relatively small for the pooled analysis, and a larger sample sizes are required to validate the ndings of our study. Second, the molecular functions of key circRNAs, miRNAs, and mRNAs in BC remained unclear. Therefore, more experiments, such as the loss or gain of function assays, should be conducted to con rm the function of these key genes in BC cells.

Conclusion
Our study constructed and analyzed a circRNA-associated ceRNA regulatory network and discovered that circ_0028899, circ_0000375, and circ_0000376 may function as ceRNAs to serve key roles in BC. This also provided new insights into the pathogenetic mechanisms of breast cancer and may provide potential therapeutic targets for it.     The ceRNA regulatory network of circRNA-miRNA-mRNA in BC. Hexagon indicate circRNAs, triangles represent miRNAs, and ovals represent miRNAs. CircRNA-miRNA-hub genes network. The network consists of three circRNAs, four miRNAs, and six hub genes. Hexagon indicate circRNAs, triangles represent miRNAs, and ovals represent miRNAs.