Identification of an lncRNA‑miRNA‑mRNA interaction mechanism in breast cancer based on bioinformatic analysis.

Non‑coding RNAs serve important roles in regulating the expression of certain genes and are involved in the principal biological processes of breast cancer. The majority of studies have focused on defining the regulatory functions of long non‑coding RNAs (lncRNAs) and microRNAs (miRNAs/miRs), and few studies have investigated how lncRNAs and miRNAs are transcriptionally regulated. In the present study, based on the breast invasive carcinoma dataset from The Cancer Genome Atlas at cBioPortal, and using a bioinformatics computational approach, an lncRNA‑miRNA‑mRNA network was constructed. The network consisted of 601 nodes and 706 edges, which represented the complex web of regulatory effects between lncRNAs, miRNAs and target genes. The results of the present study demonstrated that miR‑510 was the most potent miRNA controller and regulator of numerous target genes. In addition, it was observed that the lncRNAs PVT1, CCAT1 and linc00861 exhibited possible interactions with clinical biomarkers, including receptor tyrosine‑protein kinase erbB‑2, estrogen receptor and progesterone receptor, demonstrated using RNA‑protein interaction prediction software. The network of lncRNA‑miRNA‑mRNA interactions will facilitate further experimental studies and may be used to refine biomarker predictions for developing novel therapeutic approaches in breast cancer.

Introduction

Breast cancer was the most common malignancy in women worldwide in 2015 (1). Breast cancer is a clinically heterogeneous disease, encompassing multiple histological types, pathological characteristics and a variety of clinical behaviors, making clinical management difficult (2). Numerous studies have sought to understand the pathogenesis of breast cancer, and to determine biomarkers for use as diagnostic and prognostic tools. A primary focus of research has been to investigate the role of microRNAs (miRNAs/miRs) in breast cancer (3). The first observation of miRNA deregulation in breast cancer was by Iorio et al (4), who performed microarray analysis to evaluate the miRNA expression profiles of 76 neoplastic breast tissue and 10 healthy adjacent tissue samples, and identified 29 disordered miRNAs, including miR-10b, miR-125b, miR-145, miR-21 and miR-155, which emerged as the most consistently deregulated in breast cancer. In order to elucidate the association between miRNAs and cancer metastasis, Farazi et al (5) conducted Solexa sequencing of small RNAs from 11 healthy breast tissue samples, 17 ductal carcinoma in situ samples, 151 invasive breast carcinoma samples and 6 cell lines; 269 novel miRNAs were identified and it was demonstrated that patients with increased expression of miR-423 were more likely to develop metastasis.

In addition to miRNAs, long non-coding RNAs (lncRNAs) have emerged as important factors contributing to the development and progression of breast cancer (6). lncRNAs are a heterogeneous group of non-coding RNAs, including the newly-discovered long intervening non-coding RNAs (lincRNAs), and they are defined as larger than 200 bp in length (7). Gupta et al (8) demonstrated that lincRNAs in the HOX loci become systematically dysregulated in breast cancer, and identified the lincRNA HOTAIR to be overexpressed in primary breast tumors and metastases; therefore, HOTAIR may serve as a predictor of eventual metastasis and mortality in primary mammary tumors.

Numerous studies have been performed on breast cancer; however, further research is required to elucidate the mechanisms involved and identify novel methods of treating breast cancer. The association between lncRNAs and miRNAs in human diseases has been previously studied, including some evidence to suggest that non-coding RNAs may act as competing endogenous RNAs (ceRNAs) (9,10). Wang et al (11) provided experimental evidence that non-coding RNAs may form an lncRNA-miRNA-mRNA interaction network to regulate liver cancer; it was demonstrated that the lncRNA HULC may function as an endogenous sponge, which interacts with and downregulates miRNA-372 and thereby decreases the translational repression of its target gene.

Although a number of regulatory mechanisms influencing the association between gene expression and protein expression have been elucidated, the involvement of numerous lncRNAs and miRNAs in transcriptional regulation has not been investigated in breast cancer. In the present study, a bioinformatics approach was used to predict the regulatory mechanisms of miRNAs and lncRNAs, and to construct the lncRNA-miRNA-mRNA network in breast cancer. The results of the present study represent a view of breast cancer from a concurrent analysis of lncRNA, miRNA and mRNA.

Materials and methods

Identification of differentially expressed lncRNAs and miRNAs

In order to investigate the role of aberrant expression of lncRNAs in breast cancer, a bioinformatics approach was used to predict the differentially-expressed genes. The symbols and names of the lncRNA genes were obtained from the database of the HUGO Gene Nomenclature Committee (genenames.org). The breast invasive carcinoma dataset was obtained from The Cancer Genome Atlas (TCGA) at the cBioPortal for Cancer Genomics, consisting of 1,105 samples (cbioportal.org). By selecting the cancer study and genomic profiles, the genes of lncRNAs that have been downloaded may be entered and the information submitted. miRNAs may be selected in the same way.

lncRNA-miRNA interaction analysis

In order to identify the miRNAs which are able to target lncRNAs, the binding of lncRNAs to miRNAs (including the folded RNA predicted structure of the lncRNAs and miRNAs) was predicted using the bioinformatics tool RegRNA 2.0 (regrna2.mbc.nctu.edu.tw/detection.html). The protein sequence was obtained from the GenBank database of the National Center for Biotechnology Information (NCBI; ncbi.nlm.nih.gov/genbank). In addition, when predicting the miRNA target sites, the minimum folding free energy was set under <-25 and the system score was set to >160. An increased score indicates a stronger binding ability. Due to the large number of lncRNAs which are differentially expressed in breast cancer, the lncRNAs selected fell above 15% alteration frequency. Additionally, the lncRNA sequences associated with Homo sapiens were searched in advance using the NCBI database.

miRNA target prediction and construction of the lncRNA-miRNA-mRNA network

In order to predict which genes to target with screening miRNAs, three databases were used: TargetScan (targetscan.org); Microcosm Targets (ebi.ac.uk/Enright-srv/microcosm/htdocs/targets/v5); and PicTar (pictar.mdc-berlin.de/cgi-bin/PicTar_vertebrate.cgi). Genes which were identified using all of the three databases were used, in order to limit the number of false positive results. In addition, it was confirmed that the lncRNAs exhibited alterations >15%, in order to maximize the clarity of the network diagram. The lncRNA-miRNA-mRNA interaction network was constructed using Cytoscape software (version 3.0; http://www.cytoscape.org/download.php).

Association between clinical pathological feature and lncRNAs

In order to investigate the association between clinical pathological features and lncRNAs, computer analysis based on the RPISeq RNA-protein interaction prediction program (pridb.gdcb.iastate.edu/RPISeq), which uses random forest (RF) and support vector machine (SVM) classifiers to predict the association between lncRNAs and breast cancer-associated proteins. The protein sequences were obtained from the NCBI GenBank database.

Gene ontology (GO) analysis

In order to further investigate the biological effects of aberrantly-expressed lncRNAs and miRNAs in breast cancer, GO enrichment of the target gene were carried out using the GOrilla tool (cbl-gorilla.cs.technion.ac.il). For each GO term, a list of associated genes is returned with the most optimal at the top of the list. Each gene name is specified by the gene symbol and followed by a short description of the gene.

Results

Genetic alterations in lncRNAs and microRNAs

In the present study, among 2,772 lncRNAs analyzed, 626 lncRNAs exhibited alterations in 1–22% of breast cancer cases, with the exception of 22 invalid gene symbols. A total of 41 lncRNAs exhibited alterations in ≥10% cases, and 15 of them (CASC8, PVT1, CCAT1, CCDC26, LINC00536, LINC00861, LINC00964, PCAT1, PCAT2, LINC0097, HPYR1, BAALC-AS2, LINC00051, LINC00535) accounted for a total of 30%. The 15 identified lncRNAs with alterations in breast cancer have been well-studied previously. Zhang et al (12) evaluated the expression of CCAT1 by reverse-transcription quantitative polymerase chain reaction (RT-qPCR) analysis, demonstrating that the expression of the lncRNA CCAT1 is increased in breast cancer tissue compared with healthy adjacent tissues, and that it may be a novel biomarker of poor prognosis in patients with breast cancer. For the miRNAs, it was identified in the present study that the majority of them exhibited alterations in the 971 records of breast cancer data; specifically, 1,172/1,189 miRNAs exhibited alterations and the frequency ranged between 1 and 21%. miR-510 overexpression in breast cancer has been demonstrated to increase cell growth, migration and invasion (13). The results of the present study are presented in Fig. 1.

Figure 1.

Identification of differentially-expressed long non-coding RNAs from The Cancer Genome Atlas; alternation frequency >15%.

lncRNA-miRNA network analysis

Although it is known that miRNAs are able to target a number of protein-coding genes, little is known about whether miRNAs are able to target lncRNAs. In order to elucidate the precise mechanism underlying the role of lncRNAs and miRNAs in breast carcinoma, RegRNA 2.0 was used to analyze the functional association between lncRNAs and miRNAs. RegRNA 2.0 is an integrated web server used to compare the homologs of regulatory RNA motifs and elements against an input mRNA sequence. The predicted miRNAs from this software ought to intersect with the breast invasive carcinoma dataset from TCGA. Among the significantly differentially-expressed lncRNAs, using a threshold alteration frequency >10%, 393 miRNAs are able to exert regulatory functions on 41 lncRNAs. The lncRNAs with the highest alteration frequencies were selected, demonstrating that 13 lncRNAs, targeted by 158 miRNAs, are associated with breast cancer. It was predicted, for example, that lncRNA CCAT1 may be regulated by hsa-miR-595 and hsa-miR-345-5p, and linc00536 may be regulated by hsa-miR-93-5p and hsa-miR-214-5p. By contrast, hsa-miR-93-5p may be targeted by linc00535 and linc00964. A previous study demonstrated that miR-21 is able to regulate lncRNA expression. RT-qPCR analysis of lncRNA demonstrated that miR-21 is capable of downregulating the lncRNA GAS5 (14).

The folded RNA predicted structure of the lncRNAs and miRNAs was analyzed using RegRNA2.0 software, and the partial reliability information of pair probabilities are depicted in Fig. 2. Hairpin loops are the most frequently-predicted secondary structure, which may facilitate binding of lncRNAs to miRNAs. lncRNA CCAT1 is able to bind to hsa-miR-1178 and hsa-miR-1304-3p, and linc00861 is able to target hsa-miR-510. Additionally, previous studies demonstrated that lncRNAs harbored putative miRNA regulatory elements and served important roles in the miRNA regulatory network. Pilyugin and Irminger-Finger (15) observed that the newly-discovered lncRNA BARD1 9′L shared miR-203 and miR-101 miRNA response elements (MREs) with BARD1 mRNA in their homologous 3′ untranslated regions. The network constructed in the present study demonstrated the unknown complexity of non-coding RNA regulatory interactions and how lncRNAs may serve as important factors in the miRNA regulatory network, with the aim of elucidating the role of these interactions in disease processes.

Figure 2.

RNA fold reliability data of probable long non-coding RNA-microRNA pairs. (A) CCAT1-hsa-miR-1178. (B) PVT1-hsa-miR-5196. (C) linc00861-hsa-miR-510 (D) linc00861-hsa-miR-140-3p. linc, long intervening non-coding; miR, microRNA.

miRNA target analysis and the lncRNA-miRNA-mRNA regulatory network

A previous study demonstrated that miRNAs served a role in regulating gene expression by causing the degradation of target mRNAs. Lin et al (16) demonstrated that miR-33b inhibits breast cancer metastasis by targeting high mobility group AT-hook 2, spalt like transcription factor 4 and twist family bHLH transcription factor 1. In the present study, 35 miRNAs were identified to be differentially-expressed in breast cancer. In order to eliminate false positive rates of the target prediction, only the miRNA-mRNA pairs simultaneously predicted by ≥2 applications were taken forward. A total of 549 genes were predicted to be targets of aberrantly expressed miRNAs. The majority of the targets were associated with cancer, including cell division cycle 34 (CDC34), mitogen-activated protein kinase kinase kinase 11 (MAP3K11), AT-rich interaction domain 4B and caspase recruitment domain family member 11. The target genes exhibited functions in cell proliferation, apoptosis, invasion and metastasis. The genes were submitted to Cytoscape for visualization (Fig. 3).

Figure 3.

Interaction network of lincRNA-miRNA-mRNA in breast cancer. Diamond nodes represent lncRNAs, triangle nodes represent miRNAs and circle nodes represent mRNAs. Edges represent the possible associations between lncRNAs, miRNAs and mRNA. lincRNA, long intervening non-coding RNA; lncRNA, long non-coding RNA; miRNA, microRNA.

A miRNA-mRNA regulatory network was constructed in the present study. As presented in the network diagram (Fig. 3), a gene may be targeted by multiple miRNAs. The number of target genes of hsa-let-7e-5p is 10 (solute carrier family 35 member D2, CDC34, basic leucine zipper and W2 domains 1, adrenoceptor β-2, heart and neural crest derivatives expressed 1, glutaminyl-tRNA synthetase, adaptor related protein complex 1 sigma 1 subunit, G protein subunit γ-5, leucine rich repeats and immunoglobulin like domains 3 and coilin). Apolipoprotein C2 mRNA is able to be co-regulated by hsa-miR-637, hsa-miR-422a and hsa-miR-765. A previous study demonstrated the interaction between miR155 and tumor protein p53 inducible nuclear protein 1 (TP53INP1). TP53INP1 is a pro-apoptotic stress-induced p53 target gene. Gironella et al (17) demonstrated that TP53INP1 was a target of miR-155 using bioinformatics and cell experiments. In the present study, the regulation of miR155-TP53INP1 was demonstrated in breast cancer. Zhang et al (18) demonstrated that overexpression of miR-155 led to the downregulation of TP53INP1 and reversed the effect of TP53INP1, by promoting breast cell proliferation and suppressing cellular apoptosis. It has additionally been demonstrated that neurofibromatosis type 1 (NF1β) is targeted by miR-106a and miR-21, and that the NF1β mRNA flanking region around the core sequence for miR-21 recognition is conserved in vertebrates up to 500 bp (19). NF1β knockdown led to a phenocopy of overexpressing miR-21 in the study by Dellago et al (20). In the present study, it was hypothesized that certain common target mRNAs may be detected between different dysregulated miRNAs. From the miRNA-mRNA regulatory network constructed in the present study, the regulation of hsa-miR-199a-5p-MAP3K11 was observed. Byrnes et al (21) demonstrated that miR-199a-5p acts as a tumor suppressor in carcinogenesis, and that its downregulation contributes to enhanced cellular proliferation by targeting MAP3K11. Further details of the results of the present study are presented in Fig. 3.

The regulatory network constructed in the present study exhibited 601 nodes and 706 edges, demonstrating the complexity regulatory interactions between lncRNAs, miRNAs and target genes. However, the functional roles of a number of lncRNAs and miRNAs in the regulatory network have been rarely reported in breast cancer. Cesana et al (22) demonstrated that a muscle-specific lncRNA, linc-MD1, ‘sponges’ miR-133 and miR-135 to regulate the expression of mastermind like transcriptional coactivator 1 and myocyte enhancer factor 2C, transcription factors that activate muscle-specific gene expression. The results of this previous study indicated that there may be direct competition for miRNA binding between lncRNAs and mRNAs, and that this serves a role at the post-transcriptional level.

Analysis of clinical pathological features and lncRNA

In order to further investigate the regulatory information network associated with the clinical pathological features of breast cancer, RPISeq software was used to predict the potential target proteins of lncRNAs. Interaction probabilities generated by RPISeq range between 0 and 1. In performance evaluation experiments, predictions with probabilities >0.5 were considered ‘positive,’ indicating that the corresponding RNA and protein are likely to interact. The score of the predicted interaction possibility between lncRNAs and clinical pathological features are presented in Fig. 4. The results of the present study demonstrated that a number of lncRNAs, including PVT1, CCAT1 and linc00861, exhibited a strong possibility of interaction with receptor tyrosine protein kinase erbB-2 (HER2), estrogen receptor (ER), progesterone receptor (PR), breast cancer type 1 susceptibility protein, breast cancer type 2 susceptibility protein, cellular tumor antigen p53, phosphatase and tensin homolog, and tumor necrosis factor, with scores of RF>0.6 and SVM>0.9. It was hypothesized that lncRNAs PVT1, CCAT1 and linc00861 may govern fundamental biological processes and may be used as clinical biomarkers for diagnosis. Merry et al (23) reported that HER2 amplification affects the expression of lincRNAs, using RNA sequencing in human breast cancer. The results of the present study may provide a foundation for future experiments, leading to increased understanding of the associations between RNA-based biomarkers and clinical pathological features.

Figure 4.

Scores of the predicted interaction possibilities between long non-coding RNAs and breast cancer associated proteins, generated using RPISeq. (A) Interaction possibilities for CCAT1. (B) Interaction possibilities for linc00861. HER2, receptor tyrosine protein kinase erbB2; ER, estrogen receptor; PR, progesterone receptor; BRCA1, breast cancer type 1 susceptibility protein; BRCA2, breast cancer type 2 susceptibility protein; p53, cellular tumor antigen p53; PTEN, 4,5-trisphosphate 3-phosphatase and dual-specificity protein phosphatase PTEN; TNF, tumor necrosis factor; linc, long intervening non-coding; RF, random forest; SVM, support vector machine.

GO analysis

In order to investigate the role of lncRNAs and miRNAs in breast cancer, a systems biology approach was used to examine the functions of predicted targets. Genes (548) were recognized by gene symbol out of 549 gene terms entered. A total of 442 of these genes were associated with a GO term. The results of the enriched biological processes and molecular functions analysis, which were generated with the GOrilla tool, are presented in Fig. 5. GOrilla is a widely-used application for identifying and visualizing enriched GO terms in ranked lists of genes (24). Each differential target gene may return a P-value; if this value is <5%, it denotes an enrichment gene. Enrichment genes may contribute to multiple biological processes, including the type I interferon (IFN) signaling pathway, natural killer cell activation in the immune response, the response to double stranded RNA and B cell proliferation. The most notable biological function was the regulation of the type I IFN signaling pathway. Browne et al (25) demonstrated that type I IFN, including IFNα, inhibit the replication of human immunodeficiency virus I by upregulating the expression of important genes. An additional enrichment ranking of GO and signaling pathways indicated that cytokines may be the most enriched genes. A study of eight cytokine genes conducted by Kim et al (26) demonstrated that 74 (24%) patients were classified using eight cytokine genes as exhibiting prevalent depression, and 19 (8%) and 25 (10%) patients were classified with persistent and incident depression, respectively. These previous results supported the role of cytokines in the etiology of depression associated with breast cancer.

Figure 5.

Biological categories of microRNA target genes. Graphical representation of gene ontology-term enrichment of target genes. The diagram is colour coded according to the degree of enrichment; a deeper orange color indicates a greater degree of enrichment.

The functional analysis conducted in the present study was supported by the above previous studies, and demonstrated using a variety of bioinformatic approaches, that lncRNAs and miRNAs exerted regulatory effects on breast cancer activation, by affecting molecular functions and signaling pathways.

Discussion

The present study synthetically analyzed the effects of non-coding RNAs to provide a framework for elucidating the mechanism of regulation by an lncRNA-miRNA-mRNA network in breast cancer. Comprehensive analysis methods were used to investigate aberrantly expressed genes, miRNAs and lncRNAs in pathological processes or diseases states. High throughput screening of the interactions between lncRNAs, miRNAs and target genes in breast cancer has rarely been performed. In the present study, expression data was used to identify abnormally expressed lncRNAs and miRNAs, and an lncRNA-miRNA-mRNA regulatory network in breast cancer was constructed.

The differentially-expressed lncRNAs and miRNAs associated with breast cancer were obtained from the database of TCGA, in order to identify genes which are likely to be associated with oncogenesis. Following future experimental studies, the function of non-coding RNAs in breast cancer will be elucidated more completely. lncRNAs were initially hypothesized to be ‘transcription noise’ and termed the ‘dark genome’; however, lncRNAs are now hypothesized to serve important roles in healthy cellular functioning and tumorigenesis through a variety of mechanisms (27). miRNAs have been demonstrated to serve as potential biomarkers for cancer diagnosis and prognosis. The amplification of PVT1, for example, contributes to breast cancer pathophysiology. The dynamic expression of miR-21 may serve multiple biological roles in tumorigenesis. Si et al (28) observed that miR-21 may function as an oncogene and modulate tumorigenesis through the regulation of genes, including apoptosis regulator Bcl-2.

The interactions between lncRNAs and miRNAs were predicted using RegRNA software. A total of 549 genes were targeted by miRNAs, according to the results of three different algorithms. miR-510 is the core element in the network constructed in the present study, as the series of target genes were controlled and regulated by miR-510. In order to further investigate the biological effects of aberrantly-expressed lncRNAs and miRNAs in breast cancer, the GOrilla tool was used to analyze GO enrichment and pathways of the targets. The prediction results indicated that these target genes may be involved in a number of biological processes, including B-cell proliferation, natural killer cell activation and the adaptive immune response.

A network was constructed, depicting the association between lncRNAs, miRNAs and mRNAs in the occurrence or development of breast cancer. From the results of the present study, a connection was identified between the non-coding RNAs and mRNAs, demonstrating that non-coding RNAs possibly exert a regulatory effect on the mRNAs and vice versa. Further functional research is required to improve the clinical treatment of breast cancer.

As a number of miRNAs were associated with the clinical pathological features of breast cancer, it was hypothesized that lncRNAs may additionally be associated with clinical pathological features. The results of the present study demonstrated that a number of lncRNAs, including PVT1, CCAT1 and linc00861, exhibited the possibility of interaction with clinical biomarkers, including HER2, ER and PR, using RPISeq software. Previous studies have demonstrated that a number of miRNAs were associated with clinical biomarkers. It is possible that lncRNAs and miRNAs may co-regulate HER2, ER and PR in breast cancer. However, the functional roles of the majority of lncRNAs in cancer pathology remain to be elucidated, although certain examples, including HOTAIR and H19, have been well-studied. Liu et al (29) demonstrated that HOTAIR was a target of miR-331-3P and may act as a ceRNA, becoming a sink for miR-331-3p, and thereby modulating the de-repression of HER2 and imposing an additional level of post-transcriptional regulation. lncRNAs and miRNAs have been demonstrated to be involved in various pathophysiological processes in human disease. However, the dynamics and corresponding functions remain to be elucidated. At present, two theories have been put forward about the association between lncRNAs and miRNAs. Previous studies demonstrated that some lncRNAs, which may serve as endogenous miRNA sponges, prevent themselves from binding to mRNAs based on the ceRNA hypothesis. Salmena et al (30) outlined the ceRNA hypothesis, and suggested that mRNAs, pseudogenes and lncRNAs may crosstalk with each other by competing for miRNAs through shared MREs, thereby acting as ceRNAs. Further studies provided evidence for this hypothesis. Johnsson et al (31) observed that the lncRNA PTENP1 regulated PTEN transcription and mRNA stability by acting as a miRNA sponge. A previous study demonstrated that lncRNAs served as precursors for miRNAs. miRNAs were able to be processed from lncRNA via sequential processing by the RNase III enzymes Drosha and Dicer (32). Cai and Cullen (33) demonstrated that the lncRNA H19 was able to function as a primary miR-675 precursor. Non-coding RNAs are, therefore, not independently regulators of oncogenesis in breast cancer; lncRNAs, miRNAs and mRNAs construct an interaction network to co-regulate gene expression.

The results of the present study presented the lncRNAs and miRNAs that were possibly involved in lncRNA-miRNA-mRNA regulatory network in breast cancer. The present study provided a novel insight into the molecular mechanism of breast cancer. However, the present study exhibited a number of limitations. Certain miRNAs cannot be simultaneously recorded in the three prediction software algorithms, which may have influenced the selection of miRNAs. For example, a number of miRNAs are unable to be selected from PicTar or Microcosm Targets. The abundance of differentially-expressed genes means that the predictions in the present study may be refined further. In addition, the lncRNA-miRNA-mRNA regulatory mechanism analyzed in the present study was predicted using a bioinformatics approach; therefore, further in vivo and in vitro experiments are required.