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

BACKGROUND: The microRNA (miRNA) miR-133a-1 has been identified as a tumor suppressor in breast cancer. However, the underlying mechanisms of miR-133a-1 in breast cancer have not been fully elucidated. This study aimed to explore the targets of miR-133a-1 in breast cancer using an integrated bioinformatics approach.
METHODS: Human SKBR3 breast cancer cells were transfected with miR-133a-1 or a miRNA negative control (miRNA-NC) for 48 hours. The RNA-seq sequencing technique was performed to identify the differential expression of genes induced by miR-133a-1 overexpression. Functional enrichment analysis was conducted to determine the target genes and pathways involved in breast cancer.
RESULTS: Breast cancer patients with high levels of miR-133a-1 expression commonly showed longer overall survival compared to patients with a low level of miR-133a-1 expression. Using Cuffdiff, we identified 1,216 differentially expressed genes induced by miR-133a-1 overexpression, including 653 upregulated and 563 downregulated genes. MOCS3 was the most upregulated gene and KRT14 was the most downregulated gene. The top 10 pathways related to the differentially expressed genes were identified through Gene Ontology (GO) enrichment analysis. Sex-determining region Y-box 9 (SOX9) demonstrated the highest semantic similarities among the differentially expressed genes. Since SOX9 and CD44 were hub nodes in the protein-protein interaction network, the SOX9 gene may be a target of miR-133a-1 in breast cancer.
CONCLUSIONS: This report provides useful insights for understanding the underlying mechanisms in the pathogenesis of breast cancer. 2021 Translational Cancer Research. All rights reserved.

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, ). 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 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 (). 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 (). 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 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 (). 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).

Table 1

Top ten upregulated and downregulated genes in samples with overexpression of miR-133a-1

Gene name	Locus	Log2 fold change	t-test	P value	q value
Top ten upregulated genes
MOCS3	chr20:49575350-49578399	3.63486	3.89033	0.0002	0.00067064
CCR1	chr3:46243199-46249832	3.21679	4.80733	0.00005	0.00019055
ZNF12	chr7:6728063-6746566	3.01868	3.92279	0.00005	0.00019055
CDH5	chr16:66400524-66438689	2.8876	4.76724	0.00005	0.00019055
CPA4	chr7:129932973-129964020	2.86613	5.38712	0.00005	0.00019055
MBTD1	chr17:49254785-49337427	2.67093	4.77854	0.00005	0.00019055
IFT80	chr3:159974773-160152749	2.6679	2.63765	0.00045	0.00138782
RBM12B	chr8:94743730-94753224	2.56549	4.43156	0.00005	0.00019055
EDIL3	chr5:83236413-83680685	2.56458	8.55926	0.00005	0.00019055
ZNF253	chr19:19976713-20004293	2.5604	3.05985	0.0006	0.00178308
Top ten downregulated genes
KRT14	chr17:39738530-39743147	−6.60964	−6.2399	0.0005	0.00152285
RASD1	chr17:17397752-17399709	−4.81953	−4.59614	0.00005	0.00019055
NCF2	chr1:183524696-183560056	−4.65199	−8.12236	0.00005	0.00019055
KRT34	chr17:39533920-39538636	−4.55889	−4.97154	0.0001	0.00036075
ZFP36	chr19:39897486-39900052	−4.55068	−15.4423	0.00005	0.00019055
C11orf96	chr11:43964105-43965433	−4.42812	−4.74994	0.00015	0.00051915
CLDN4	chr7:73245192-73247023	−4.33659	−9.81406	0.00005	0.00019055
OAS2	chr12:113416273-113449528	−4.24137	−6.99281	0.00005	0.00019055
CCL5	chr17:34198495-34207377	−4.18885	−15.8991	0.00005	0.00019055
DHRS2	chr14:24105572-24114848	−4.13384	−8.46067	0.00005	0.00019055

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 ().

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 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 . 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 . 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 (). 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 , 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

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.