An integrated bioinformatics analysis to investigate the targets of miR-133a-1 in breast cancer

Introduction

MicroRNAs (miRNAs) are small, non-coding RNAs with an approximate length of 22 nucleotides. They are found in many organisms and play an important role in the regulation of gene expression at the post-transcriptional level (1,2) via degradation of RNA targets or inhibition of messenger RNA (mRNA) translation (3). Increasing evidence has shown that miRNAs are aberrantly expressed in multiple cancers and more than 50% of miRNA genes are located in cancer-associated genomic regions (4-6), suggesting a crucial role of miRNAs as oncogenes or tumor suppressors in the pathogenesis of cancers.

Profiling studies have reported the abnormal expression of miRNAs in human breast cancers and their important functions in tumor pathogenesis. In MCF7 breast cancer cells, miR-206 was shown to reduce the expression of the estrogen receptor (ER)-α and to suppress cell growth (7). Other studies have reported that miR-155 could promote tumor angiogenesis by targeting the von Hippel-Lindau (VHL) gene in breast cancer cells and was associated with poor prognosis (8). Owing to the large number of miRNAs and their crucial role in breast cancer, breast cancer-related miRNAs and their underlying mechanisms deserve to be further investigated.

A single miRNA could act as a regulator of hundreds of genes, and a transcript could be regulated by numerous different miRNAs. With the development of high-throughput profiling and comprehensive bioinformatics analysis, it is possible to identify multiple miRNA-associated genes using next-generation sequencing. This may help in investigating miRNA-targeted genes or cellular pathways involved in cancer. Ye et al. performed next-generation sequencing to explore the alterations in the whole transcriptome induced by overexpression of miR-145. Their study identified 49 mRNAs targeted by miR-145 in breast cancer cells (9). Considering that miR-133a can function as a tumor suppressor in breast cancer (10), this study set out to explore the multiple miR-133a-1 associated target genes and their pathways using next-generation sequencing followed by an integrated bioinformatics analysis. We present the following article in accordance with the MDAR checklist (available at: http://dx.doi.org/10.21037/tcr-20-2681).

Methods

Gene expression database

The gene expression profiles of 916 patients with breast cancer were downloaded from the Kmplotter database (www.kmplot.com). Data analyses were conducted in strict accordance with The Cancer Genome Atlas (TCGA) publication guidelines (http://cancergenome.nih.gov/publications/publicationguidelines). Kaplan-Meier analysis was performed to assess the differences in survival between patients with high and low levels of miR-133a-1 expression using the “survival” package in R (version 3.24). Statistical significance was compared using the log-rank test. A P value of <0.05 was considered to represent statistical significance. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Cells and transfection

Human miR-133a-1 and miR-negative control (NC) were purchased from Vigene Biosciences (Jinan, China; Cat. #CR207927). The human breast cancer cell line SKBR3 was purchased from ATCC (American Type Culture Collection; Manassas, VA, USA). The SKBR3 cells were seeded onto 6-well plates and cultured in ATCC-formulated McCoy’s 5a Medium Modified (Catalog No. 30-2007; Manassas, VA, USA) plus 10% fetal bovine serum (Catalog No. 30-2020; Manassas, VA, USA) overnight at 37 °C. Cells were then transfected with miR-133a-1 or mi-NC using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. After 48 hours, the cells were harvested for further experiments.

Total RNA extraction and quantitative real-time polymerase chain reaction analysis

Total RNA was isolated using the RNeasy Mini kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The RNA-sequencing (RNA-seq) library was constructed using Hieff NGSâ MaxUP II DNA Library Prep Kit for Illumina (Yeasen, Shanghai, China) following the manufacturer’s instructions. Sequencing was performed by Novogene Biotech (Beijing, China) using the HiSeq X Ten system (Illumina, San Diego, CA, USA). Quality control was carried out using Fastp, and clean reads were mapped to the human genome (hg19) using Tophat.

Western blot analysis

For immunoblotting, nitrocellulose filter membranes were incubated with the sex-determining region Y-box 9 (SOX9) antibody (cat. #82630, 1:1,000), MSH2 (cat. #2850, 1:1,000), CD44 (cat. #5640, 1:1,000), and actin (cat. #3700, 1:1,000). All antibodies were purchased from Cell Signaling Technology, Inc. (Danvers, MA, USA).

Statistical analysis

Cuffcompare was used to compare the similarity of transcripts and assess the construction of transcripts. Cuffmerge was used to combine multiple transcript sets into a single transcript set. The gene expression pattern was analyzed using the PlotPCA package, and heatmaps were drawn according to the gene expression using the Pheatmap package in the R language. Cuffdiff was used to screen the differentially expressed genes. The screening threshold for differentially expressed genes in Cuffdiff was false discovery rate (FDR) <0.05, abs (log foldchange) >1; that is, P value <0.05 and a fold change of >2 after correction. A volcano plot was drawn using Ggplot2 package. The distribution of differentially expressed genes was mapped to chromosomes using the graphics function of R language. Gene Ontology (GO) enrichment analysis of the differentially expressed genes was conducted using the Clusterprofiler package of the R language. The similarity between GO terms was calculated using GOSemSim and the GO terms cluster was labeled by Ggtree. The potential downstream genes of miR-133a-1 were predicted with TargetScan. Combined with the differentially expressed genes from the RNA-seq, the intersection was used for semantic similarity analysis using GOSemSim. The protein-protein interaction (PPI) network was analyzed using the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database (https://string-db.org/).

Results

High expression of miR-133a-1 was associated with favorable overall survival

To investigate the potential function of miR-133a-1 in breast cancer, the expression profiles of 916 patients with breast cancer were downloaded from the Kmplotter database. Patients were stratified into two sub-cohorts using the median value of miR-113a-1 expression of all the patients as the cutoff value. A total of 926 patients were included in the high miR-133a-1 expression subgroup, while 326 patients were included in the low miR-133a-1 expression subgroup. Kaplan-Meier analysis revealed that patients with low levels of miR-133a-1 expression showed superiority in median overall survival compared to patients with a high expression of miR-133a-1, with a hazard ratio (HR) of 1.4, and a 95% confidence interval (CI) of 1.1 to 1.8 (P=0.0056, Figure 1). These results indicated that miR-133a-1 may be a potential biomarker of overall survival in patients with breast cancer.

Figure 1 Kaplan-Meier analysis of overall survival in breast cancer patients with high or low expression of miR-133a-1.

MiR-133a-1 overexpression in breast cancer cells triggered variations in the mRNA expression profiles

To explore the potential target genes of miR-133a-1 in breast cancer, RNA-seq was performed on SKBR3 cells transfected with miR-133a-1 or the negative control miR-NC for 48 hours. Real-time quantitative PCR revealed that the levels of miR-133a-1 expression in cells transfected with miR-133a-1 were significantly higher than those in cells transfected with the negative control (Figure 2A).

Figure 2 Transfection of miR-133a-1 in SKBR3 cells. (A) MiR-133a-1 levels in SKBR3 cells following transfection with miR-133a-1 or miR-NC. (B) A scatter plot of the differential genomic profiles in cells overexpressing miR-133a-1 or miR-NC using principal component analysis (PCA). A PCA plot showing variations in the three replicates of miR-133a-1 or miR-NC. The percentages in axis indicate the percentages of variance explained by the two principal components. NC, negative control.

A sequencing library was constructed for further RNA-seq analysis. A principal component analysis was conducted to assess the similarity or dissimilarity of gene expression patterns between cells overexpressing miR-133a-1 and cells transfected with the miR-NC. Cells overexpressing miR-133a-1 and cells transfected with miR-NC were grouped into distinct clusters, and the first two principal components explained more than 80% of the variability between the samples (Figure 2). This finding indicated the presence of differential gene expression patterns between the two samples, and confirmed the excellent performance of the RNA-seq.

The differential expression of genes

Hierarchical cluster analysis demonstrated the comprehensive differential gene expression patterns of samples with miR-133a-1 overexpression and cells transfected with miR-NC (Figure 3A). The results showed that the gene expression profiles of cells overexpressing miR-133a-1 were markedly different from those of the negative control cells. The immense impact of miR-133a-1 on gene expression in SKBR3 cells further supported the key role of miR-133a-1 in breast cancer.

Figure 3 Differential gene expression in cells overexpressing miR-133a-1 or miR-NC. (A) A heat map demonstrating the differential expression patterns between cells overexpressing miR-133a-1 or miR-NC. The horizontal axis indicates the name of samples. The first three samples are replicates of miR-133a, and the other three samples are replicates of miR-NC. The top horizontal axis indicates the clusters of samples. Genes with high expression levels are marked in red and those with low expression levels are marked in green. (B) A volcano plot presenting the differentially expressed genes in cells with miR-133a-1 overexpression and NC cells. The x-axis represents the log2 fold change and the y-axis indicates the −log10 (q value, which equals the adjusted P value). Genes upregulated in cells with miR-133a-1 overexpression are presented as red dots and those with downregulated miR-133a-1 expression are presented as green dots. Genes that were not differentially expressed are presented as black dots. (C) The distribution of the differentially expressed genes in the whole genome. Upregulated genes are marked in red and downregulated genes are marked in green. NC, negative control.

The volcano plot in Figure 3B visualizes 1,216 differentially expressed genes between cells overexpressing miR-133a-1 and cells transfected with miR-NC, including 653 upregulated genes and 563 downregulated genes. The top 10 upregulated genes were MOCS3, CCR1, ZNF12, CDH5, CPA4, MBTD1, IFT80, RBM12B, EDIL3, and ZNF253. The top 10 downregulated genes were KRT14, RASD1, NCF2, KRT34, ZFP36, C11orf96, CLDN4, OAS2, CCL5, and DHRS2 (Table 1). BRCA1, the well-known breast cancer-related gene, was also found to be significantly downregulated in samples with miR-133a-1 overexpression (log2 fold change =1.08502; P value <0.001, q value <0.001).

To investigate the enrichment status of the genes with differential expression caused by miR-133a-1 overexpression, the chromosomal locations of the 1,216 genes were examined. The results revealed that the differentially expressed genes were located throughout the genome, with no enrichment preference for any genomic region observed (Figure 3C).

Analysis of the functions and pathways related to the differential gene expression

To investigate the pathways related to the differential expression of the 653 upregulated and 563 downregulated genes, GO enrichment analysis was performed using the Clusterprofiler package. An adjusted P value of <0.05 was considered to be statistically significant. The top 10 pathways in which the differentially expressed genes were significantly enriched included negative regulation of protein phosphorylation, regulation of the intrinsic and extrinsic apoptotic signaling pathways, embryonic organ development, gland morphogenesis, ossification, response to viruses, and positive regulation of chromosome organization (Figure 4A).

Figure 4 Functional and pathway enrichment analysis. (A) Gene Ontology (GO) terms enriched in cells with miR-133a-1 overexpression. (B) Clustering dendrogram presenting the GO-based pathway clusters which were associated with miR-133a-1 overexpression. FC, fold change.

Cluster similarity analyses were conducted on the basis of the GO results. Four GO-based pathway clusters containing more than 1 pathway were generated and are presented in the clustering dendrogram in Figure 4B. The observation that miR-133a-1 regulated genetic events across subgroups and multiple pathways further supported the hypothesis that it plays an important role in breast cancer.

The miR-133a-1 associated target genes

TargetScan was used to predict the potential target genes of miR-133a-1. The selected genes were combined with the differentially expressed genes from our sequencing data, and the intersection was used for the analysis of semantic similarity using GOSemSim. Twelve genes were selected; their semantic similarities were ranked and are presented in the box plot shown in Figure 5A. Sex determining region Y-box 9 (SOX9) showed the highest semantic similarity among these genes, followed by neural precursor cell expressed, developmentally down-regulated 9 (NEDD9), and msh homeobox 2 (MSX2). These results suggested the potential role of SOX9 in the miR-133a-1-associated gene network. In addition, from the GO database (http://www.geneontology.org), we found that SOX9 was enriched in embryonic organ development pathway and ossification pathway, which indicating the potential regulation of the two pathways by miR-133a-1.

Figure 5 In silico analysis revealed that SOX9 is a crucial target of miR-133a-1. (A) The box plot shows the semantic similarity of potential targeted genes of miR-133a-1. (B) The protein-protein interaction network analysis of the differentially expressed genes obtained from the STRING database.

To further validate the potential target genes of miR-133a, the PPI network was analyzed. The 12 genes obtained from the semantic similarity analysis were uploaded to the STRING website. The network of potential target genes of miR-133a-1 was constructed with 12 nodes and 11 edges (Figure 5B). SOX9 and CD44 were hub nodes in the network, which further supported the hypothesis that the SOX9 gene might be a target of miR-133a. Among the genes presented in Figure 5B, MSX2 demonstrated the closest relationship with SOX9, suggesting a potential role of MSX2 in the miR-133a-1-SOX9 axis.

Real-time quantitative PCR and Western blot analyses were performed to validate the expression of SOX9, MSH2, and CD44 in miR-133a-overexpressing breast cancer cells. The results revealed that both the gene and protein expressions of SOX9, MSH2, and CD44 were downregulated by miR-133a (Figure 6).

Figure 6 Quantitative real-time PCR and Western blotting both validated the downregulation of SOX9, MSH2 and CD44 in miR-133a overexpressed breast cancer cells. (A) Real-time quantitative PCR. (B) Western blotting. SKBR3 cells were transfected with miR-133a-1 or miR-NC using Lipofectamine 2000 and incubated for 48 hours before measurement. *, P<0.05. NC, negative control.

Discussion

In this study, we demonstrated that low levels of miR-133a-1 expression were associated with longer median survival in patients with breast cancer, thus supporting the potential role of miR-133a-1 as a prognostic biomarker. Integrated bioinformatics analysis identified multiple potential target genes of miR-133a-1 and associated pathways in breast cancer cells. Semantic similarity analysis demonstrated that SOX9 showed the highest semantic similarity among the differentially expressed genes. PPI analysis confirmed the key role of SOX9 in the miR-133a-1-associated PPI network.

Several miR-133a-1-upregulated and downregulated genes identified in this study have been reported to play important roles in breast cancer cells. It has been demonstrated that the expression of the CCR1 is higher in breast invasive ductal carcinoma tissue than in adjacent normal tissues, and silencing of CCR1 can weaken the invasive and metastatic capabilities of breast cancer cells (11). CDH5 has been shown to stimulate the progression of breast cancer via the TGF-β signaling pathway (12). Experiments in MCF-7 cells have shown high levels of CPA4 secretion, indicating a role of CPA4 in the progression of breast cancer (13). Upregulation of RASD1 has also been shown to result in apoptosis of MCF-7 breast cancer cells (14). All of the above reports concur with and support the findings of the current study. Moreover, we also identified several novel potential targets of miR-133a-1, including MOCS3, ZNF12, MBTD1, IFT80, KRT14, NCF2, and KRT34. Further investigations are warranted to elucidate the underlying mechanisms of miR-133a-1 and its targets in breast cancers.

GO enrichment analysis generated several pathways in which the differentially expressed genes were abundantly enriched. These pathways included the negative regulation of protein phosphorylation, the apoptotic signaling pathway, and the development of embryonic organs. The association of miR-133a-1 with these pathways have been reported in other diseases and were consistent with our results. Upregulation of miR-133a reduced phosphorylation of ERK1/2 in cardiac hypertrophy (15) and miR-133 acted as an anti-apoptotic regulator in cardiomyocytes (16). A microarray study revealed that expression of miR-133 was increased during the development of mouse embryos (17). However, the relationships between these pathways and miR-133a-1 have not been explored in breast cancers previously. Here, we reported that genes in these pathways were regulated by miR-133a-1 overexpression in breast cancer cells, and the mechanisms involved warrant further investigation.

Previous studies have reported that miR-133b regulates the pathogenesis and metastasis of breast cancer by targeting SOX9 (18). MiR-133a, an isoform of miR-133 that has only 1 base difference from miR-133b at the 3'-end (19), was reported to suppress breast tumor growth (10). However, the correlation of miR-133a-1 with SOX9 has not been investigated. Here, we observed that SOX9 had the highest semantic similarities among all the miR-133a-1-induced differentially expressed genes, and SOX9 was the hub node in the PPI network. These observations provided a basic understanding of the function of miR-133a-1 in breast cancer via targeting of SOX9.

SOX9 is a transcription factor which regulates several developmental and tumor pathways. Previous studies have described the involvement of SOX9 in tumor initiation, proliferation, migration, and metastasis (20). SOX9 was found to contribute to all the miR-133b-mediated effects including cell proliferation, colony formation, and invasion ability in vitro, while re-expression of SOX9 could reverse miR-133b-mediated metastasis suppression in vivo (18). Therefore, we speculated that SOX9 might also function downstream of miR-133a by regulating cell proliferation, colony formation, invasion, and metastasis in breast cancer cells. Targeting SOX9 may be a potential therapeutic strategy for the treatment of breast cancers. However, further experiments are required to validate the function of miR-133a and SOX9 in breast cancers.

In conclusion, we performed an integrated bioinformatics analysis to investigate the role of miR-133a-1 in breast cancer. We identified several differentially expressed genes in breast cancer that are regulated by miR-133a-1, as well as their related pathways. This comprehensive investigation may contribute to a deeper understanding of cancer-related miRNAs and their targets, and provide a basis for further exploration of the mechanisms involved in the pathogenesis of breast cancer.

Acknowledgments

Funding: This work was supported by grant from National Natural Science Foundation of China (81502682, 81902996), Natural Foundation of Jiangsu Province (BK20140735), Fudan University Shanghai Cancer Center Fund (YJ201701), Natural Foundation of Shanghai (16411964900), National Science and Technology Major Project (2020ZX09201-013), and the Pudong New District Health and Family Planning Commission Youth Science and Technology Project (PW2016B-4).

Footnote

Reporting Checklist: The authors have completed the MDAR checklist. Available at http://dx.doi.org/10.21037/tcr-20-2681

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at http://dx.doi.org/10.21037/tcr-20-2681). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Lin S, Gregory RI. MicroRNA biogenesis pathways in cancer. Nat Rev Cancer 2015;15:321-33. [Crossref] [PubMed]
2. Inui M, Martello G, Piccolo S. MicroRNA control of signal transduction. Nat Rev Mol Cell Biol 2010;11:252-63. [Crossref] [PubMed]
3. Mattick JS, Gagen MJ. The evolution of controlled multitasked gene networks: the role of introns and other noncoding RNAs in the development of complex organisms. Mol Biol Evol 2001;18:1611-30. [Crossref] [PubMed]
4. Di Leva G, Garofalo M, Croce CM. MicroRNAs in cancer. Annu Rev Pathol 2014;9:287-314. [Crossref] [PubMed]
5. Iorio MV, Croce CM. microRNA involvement in human cancer. Carcinogenesis 2012;33:1126-33. [Crossref] [PubMed]
6. Jansson MD, Lund AH. MicroRNA and cancer. Mol Oncol 2012;6:590-610. [Crossref] [PubMed]
7. Adams BD, Furneaux H, White BA. The micro-ribonucleic acid (miRNA) miR-206 targets the human estrogen receptor-alpha (ERalpha) and represses ERalpha messenger RNA and protein expression in breast cancer cell lines. Mol Endocrinol 2007;21:1132-47. [Crossref] [PubMed]
8. Kong W, He L, Richards EJ, et al. Upregulation of miRNA-155 promotes tumour angiogenesis by targeting VHL and is associated with poor prognosis and triple-negative breast cancer. Oncogene 2014;33:679-89. [Crossref] [PubMed]
9. Ye P, Shi Y, An N, et al. miR-145 overexpression triggers alteration of the whole transcriptome and inhibits breast cancer development. Biomed Pharmacother 2018;100:72-82. [Crossref] [PubMed]
10. Sui Y, Zhang X, Yang H, et al. MicroRNA-133a acts as a tumour suppressor in breast cancer through targeting LASP1. Oncol Rep 2018;39:473-82. [PubMed]
11. Shin SY, Lee DH, Lee J, et al. C-C motif chemokine receptor 1 (CCR1) is a target of the EGF-AKT-mTOR-STAT3 signaling axis in breast cancer cells. Oncotarget 2017;8:94591-605. [Crossref] [PubMed]
12. Labelle M, Schnittler HJ, Aust DE, et al. Vascular endothelial cadherin promotes breast cancer progression via transforming growth factor beta signaling. Cancer Res 2008;68:1388-97. [Crossref] [PubMed]
13. Tan AA, Phang WM, Gopinath SC, et al. Revealing Glycoproteins in the Secretome of MCF-7 Human Breast Cancer Cells. Biomed Res Int 2015;2015:453289. [Crossref] [PubMed]
14. Tian J, Duan YX, Bei CY, et al. Calycosin induces apoptosis by upregulation of RASD1 in human breast cancer cells MCF-7. Horm Metab Res 2013;45:593-8. [Crossref] [PubMed]
15. Li Y, Cai X, Guan Y, et al. Adiponectin Upregulates MiR-133a in Cardiac Hypertrophy through AMPK Activation and Reduced ERK1/2 Phosphorylation. PLoS One 2016;11:e0148482. [Crossref] [PubMed]
16. Xu C, Lu Y, Pan Z, et al. The muscle-specific microRNAs miR-1 and miR-133 produce opposing effects on apoptosis by targeting HSP60, HSP70 and caspase-9 in cardiomyocytes. J Cell Sci 2007;120:3045-52. [Crossref] [PubMed]
17. Care A, Catalucci D, Felicetti F, et al. MicroRNA-133 controls cardiac hypertrophy. Nat Med 2007;13:613-8. [Crossref] [PubMed]
18. Wang QY, Zhou CX, Zhan MN, et al. MiR-133b targets Sox9 to control pathogenesis and metastasis of breast cancer. Cell Death Dis 2018;9:752. [Crossref] [PubMed]
19. Wang Y, Li L, Moore BT, et al. MiR-133a in human circulating monocytes: a potential biomarker associated with postmenopausal osteoporosis. PLoS One 2012;7:e34641. [Crossref] [PubMed]
20. Jana S, Madhu Krishna B, Singhal J, et al. SOX9: The master regulator of cell fate in breast cancer. Biochem Pharmacol 2020;174:113789. [Crossref] [PubMed]