Downregulation of Hsa-miR-3616-3p in Triple-Negative Breast Cancer is Associated with Cell Migration and Invasion


 Background: Micro RNA (miRNA), a class of endogenous small RNAs with a length of about 20-24 nucleotides, have been played a variety of important regulatory roles in cells and have attracted the attention of researchers because involvement of initiation and progression of diverse kinds of human diseases, especially cancer. However, the regulatory role of miRNA in triple negative breast cancer (TNBC), its clinical significance and potential mechanism are still largely unknown.Methods: In this study, the differentially expressed miRNAs were analyzed using GEO2R from GSE57897 and GSE97811 of the Gene Expression Omnibus (GEO). Then, we performed the overall survival (OS) analysis of four candidate miRNAs. And the expressions of four candidate miRNAs in TNBC and non-TNBC tissues were tested by Quantitative real-time PCR. We determined the level and prognostic values of hsa-miR-3616-3p in TNBC patients. Then, TNBC cells migration and invasion were studied in vitro with hsa-miR-3616-3p by using wounding heal assays and transwell assays. Next, target genes and transcription factors of hsa-miR-3616-3p enrichment analysis, were performed by the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.Results: Based on GSE57897 and GSE97811, we found 962 differentially expressed miRNAs including 943 up-regulated miRNAs and 19 down-regulated miRNAs. We selected top2 up-regulated miRNAs as hsa-miR-4428 and hsa-miR-575, and top2 down-regulated miRNAs as hsa-miR-3616-3p and hsa-miR-3909, and then we found that low expression of 4 candidated miRNAs had a worse prognosis in TNBC subtype. Next, we verified that the expression of hsa-miR-3616-3p was the most downregulated in TNBC tissues compared with matched surrounding tissues by qRT-PCR. Moerover, our results showed that overexpression of hsa-miR-3616-3p inhibited TNBC cells migration, and invasion in vitro.Conclusions: Our study might reveal important roles of hsa-miR-3616-3p plays in TNBC migration and invasion.


Introduction
Breast cancer (BC) is one of the most popular malignant tumors among women and the second major cause of cancer-related deaths in women worldwide [1]. TNBC accounts for 15% ~ 20% of all the pathological types of breast cancer without expression of estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor2 (HER2) [2]. It has special biological behavior and clinicopathological characteristics. Compared with other breast cancer subtypes, TNBC occurs more frequently in young women, usually with a high histologic grade and poor prognosis [3]. The most important characteristics of TNBC are high proliferation and frequent metastasis to the lung and brain [4]. Although decades of research have provided considerable insight into the multistep transfer process, little is known about the molecular basis of TNBC transfer. The main therapeutic strategies are conventional chemotherapy and radiation. To date, no targeted therapy for TNBC has been approved by Food and Drug Administration (FDA) [5,6]. Therefore, it is urgent to nd out prognosis predictive factors and potential therapeutic targets in BC diagnosis and treatment. And, further elucidation of the molecular mechanisms of TNBC metastasis is essential for the development of new treatments and the successful treatment of TNBC patients.
MiRNA is a highly conserved non-coding small RNA, of 21-23 nucleotides that can regulate gene expression at the post-transcriptional level by inhibiting protein translation or degradation of target mRNA [7]. As we all know, microRNAs can serve as either an oncogene or a tumor suppressor in the occurrence and development of various cancers, including breast cancer [8,9]. Some miRNAs were also found to regulate cancer cell growth and metastasis in lung cancer, liver cancer, colorectal cancer, gastric cancer, renal cancer, prostate cancer,breast cancer and other types of cancer [10][11][12][13][14][15][16]. To date, the functions of only a few miRNAs have been well characterized, and the biological functions of most miRNAs remain unknown.
In our study, the differentially expressed miRNAs were analyzed using GEO2R from GSE57897 and GSE97811 of the Gene Expression Omnibus (GEO). We found 4 candidate miRNAs and performed the overall survival (OS) analysis. The results showed that a high expression of hsa-miR-4428, hsa-miR-575, hsa-miR-3909 and hsa-miR-3616-3p were associated with poor prognosis in the BC, but great prognosis in the TNBC. Therefore, we further veri ed the relationship between 4 candidate miRNAs and TNBC, we detected the expression of four candidate miRNAs in 15 TNBC tissues and 15 non-TNBC tissues by using real-time qPCR. Our qPCR data showed that the expressions of hsa-miR-575, hsa-miR-3909 and hsa-miR-3616-3p were consistent with the result of the microarray data and hsa-miR-3616-3p was the most signi cant and stable down-regulation in TNBC tissues. Hsa-miR-3616-3p overexpression inhibited migration and invasion of TNBC cells in vitro. We found that target genes and transcription factors of hsa-miR-3616-3p were associated with TNBC. Our work may reveal the important role of hsa-miR-3616-3p in the occurrence and development of TNBC, and hsa-miR-3616-3p may be clinically valuable biomarkers and novel therapeutic targets. The Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) is a public genome database [17]. In this study, two Noncoding RNA pro les GSE57897 and GSE97811 were downloaded from the GEO. The criteria for both Non-coding RNA pro les were (a) the samples included two groups of breast cancer tissue samples and normal breast tissue samples; (b) the sample size of each Non-coding RNA pro le was greater than 60; (c) they were recently updated in the last two years (2019-2020). The GSE57897 and GSE97811 were derived from GPL18722 and GPL21263 platform, respectively. The GSE57897 data pro le contained 422 breast cancer tissue samples and 31 normal breast tissue samples. The GSE97811 data pro le contained 45 breast cancer tissue samples and 16 normal breast tissue samples.
2.2 Identi cation of differentially expressed miRNAs GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) online analysis software was used to analyze the differentially expressed miRNAs in breast cancer tissue samples and normal breast tissue samples. They were analyzed by Limma rapid differential analysis tool in sanger box (http://sangerbox.com/). p<0.05 was used as the cut-off criteria of differentially expressed miRNAs.

Survival analysis of candidate miRNAs
The OS analysis of 4 candidate miRNAs was performed using the online Kaplan-Meier Plotter survival analysis server (https://kmplot.com/analysis/). The p<0.05 was considered to indicate statistical signi cance, and screened out promising miRNAs with the prognostic value.

Tissue samples
The samples of 15 TNBC tissues and 15 non-TNBC tissues were from patients in the First A liated Hospital of Chongqing Medical University (Chongqing, China). Tissue samples of patients were treated without chemotherapy or radiotherapy before surgical resection and all patients gave informed consent. The tissue samples were quickly cryopreserved with liquid nitrogen tanks for further analysis. All biopsy tissue samples of patients were diagnosed as TNBC or non-TNBC by immunohistochemical staining. The study was conducted in accordance with the Helsinki Declaration and has been approved by the Ethics Committee of Chongqing Medical University.

RNA extraction and the qRT-PCR analysis of miRNAs expression
Total RNA of TNBC or non-TNBC tissue samples and TNBC cells was extracted by RNAiso Reagent (TaKaRa, China) under the manufacturer's instructions. To approve the bioinformatics results, the extracted RNA was reverse-transcribed into cDNA by using the PrimeScript RT Reagent Kit (TaKaRa, China), and then detect the fold changes of miRNAs by using the qRT-PCR. The qRT-PCR response setting procedures were as follows: 3 min for 95℃, 40 cycles at (95℃ for 5 s, 58℃ for 35s), 72℃ for 60s ), 65℃ for 5 s, 95℃ for 0.5s.
All experiments were conducted more than three times. The expression levels of miRNAs were de ned based on the cycle threshold (Ct), and relative expression levels were calculated using the 2-ΔΔCt method. The expression levels of U6 was used as reference miRNA. The primers used for the experimental validation are listed in Table 1.
qRT-PCR was used to measure transfection e ciency.

Wound healing assays
Cell migration was mensurated by wound healing assays. MDA-MB-231 cells were transfected with mimic NC or hsa-miR-3616-3p mimic and inhibitor NC or hsa-miR-3616-3p inhibitor, and when the cells reached 90% con uence, the monolayer of the cells was scratched with the tip of the 200 microliter pipette to produce wound stripes. Cells migration into the wound was observed through ve randomly selected microscopic horizons at two predetermined time points (0 and 24 h), respectively. The area of the remaining gaps after cell growth after 0 and 24 h intervals were measured, and the area of the gap in NC group was set to 0 h as 100% and relative cell migration was calculated. A larger gap is considered to be lower cell migration. All experiments were conducted at least three times.

Transwell assays
Transwell assays were performed to evaluate cell migration and invasion capability. Transfected MDA-MB-231 cells were starved in serum-free medium for two hours, detached, and resuspended in medium containing serum-free. For migration assays, 1.5× TNBC cells were suspended in 200 µL serum-free medium and inoculated into the upper chambers (BD BioCoat, MA, USA), and then 500 µL complete medium was added into the bottom chambers. For invasion assays, 3×TNBC cells were suspended in 200 µL serum-free medium and inoculated into the upper chambers coated with matrigel (BD Biosciences, NJ, USA), and then 500 µL complete medium was added into the bottom chambers. After 24 h, the cells on the upper chambers were removed and cells on the lower compartment were xed with ethanol and stained by crystal violet, then photographed and counted with a microscope (Leica, Wetzlar, Germany).

GO and KEGG enrichment analysis of target genes of hsa-miR-3616-3p
The DAVID (https://david.ncifcrf.gov/) is a database that provides systematic and comprehensive annotation information of biological functions for genes and proteins [18]. Gene ontology (GO) is an international standard classi cation system of gene function [19]. It includes biological process (BP), cellular component (CC), and molecular function (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database for systematic analysis of gene function and genomic information [20]. We used DAVID to carry out GO functional annotation and KEGG pathway enrichment analysis on the target genes of hsa-miR-3616-3p, p < 0.05 were considered statistically signi cant [21].
The function of TFs were conducted by KEGG enrichment analysis. The expression of TFs in normal tissues and TNBC tissue samples from TCGA were analyzed by the Ualcan database (http://ualcan.path.uab.edu/).

Statistical analysis
Statistical analyses were performed in GraphPad Prism 8.0 software. The values of different groups were represented by the mean±SD. The comparison of expression levels of TNBC and non-TNBC tumor tissues were analyzed by paired t test. p<0.05 was considered to represent a statistically signi cant difference.

Identi cation of differentially expressed Non-coding RNAs
To explore the role of systems biology in the pathogenesis of BC, we analyzed two chip data of GSE57897 and GSE97811 by GEO2R. They were normalized by Limma rapid differential analysis tool in sanger box (http://sangerbox.com/) to eliminate batch effects. Based on adjust p value < 0.05 and |log2 Fold Change| ≥ 1 were used as the cut-off criteria, the differentially expressed miRNAs were derived 962, included 943 upregulated miRNAs and 19 downregulated miRNAs between breast tissue samples and normal breast tissue samples ( Figure 1A). We selected 4 candidate miRNAs for further analysis: top2 up-regulated miRNAs as hsa-miR-4428 and hsa-miR-575, and top2 down-regulated miRNAs as hsa-miR-3616-3p and hsa-miR-3909 ( Figure 1B and Table 2).

OS analysis of candidate miRNAs
OS is considered to be the best end point in tumor clinical trials. When the survival time of patients can be fully evaluated, it is the rst end point. The prognostic value of 4 candidate miRNAs was measured by the online Kaplan-Meier Plotter survival analysis server in three subtypes of BC, namely Luminal (positive for ER and/or PR), Her-2 overexpression (negative for ER and PR, positive for Her-2), and TNBC. Obtain OS of patients based on the expression of 4 candidate miRNAs in breast cancer and its three subtypes. Based on the screen criteria p<0.05, OS analysis of 4 candidate miRNAs in BC patients and luminal subtype of BC patients, the results showed the patients with high expression of 4 candidate miRNAs were related to a poor prognosis (Figure 2A and Figure 2B). As for Her-2 overexpression subtype, the results showed that the high expression of hsa-miR-575 was related to a poor prognosis, and 3 candidate miRNAs (hsa-miR-4428, hsa-miR-3909 and hsa-miR-3616-3p) were not associated with OS (p>0.05) ( Figure 2C). In TNBC subtype, the low expression of 4 candidate miRNAs exhibited shorter OS time than those with high expression ( Figure 2D). Highly aggressive TNBC was considered to have less favorable prognosis compared to other BC types. Thus, we further explored the relationship between 4 candidate miRNAs and TNBC.
Moreover, hsa-miR-3616-3p was the most signi cant and stable down-regulation, and the low hsa-miR-3616-3p expression was related to a poor prognosis in TNBC. All of the above results suggested that hsa-miR-3616-3p may be correlated with TNBC. Therefore, we focused on hsa-miR-3616-3p for further study.

Down-regulation of hsa-miR-3616-3p enhances migration and invasion in TNBC cells
To further explore whether the expression of hsa-miR-3616-3p could regulate cell migration and invasion in TNBC, MDA-MB-231 cells were transfected with mimic NC or hsa-miR-3616-3p mimic, inhibitor NC or hsa-miR-3616-3p inhibitor, respectively. The transfection e ciency was measured by qRT-PCR ( Figure 4A). Next, wound healing assays and transwell assays were performed to examine the effect of hsa-miR-3616-3p with migration and invasion in MDA-MB-231 cell ( Figure 4B -4D). All of the above results showed that up-regulation of hsa-miR-3616-3p inhibited the migratory and invasive capabilities of MDA-MB-231 cells, while down-regulation of hsa-miR-3616-3p revealed an opposite effect. These data suggested that hsa-miR-3616-3p should exert an anti-oncogenic role in the TNBC, and the downregulation of hsa-miR-3616-3p induced migratory and invasive in TNBC cells.

GO and KEGG Enrichment Analysis of Target genes of hsa-miR-3616-3p
In order to better understand the biological function of hsa-miR-3616-3p, we found 4,218 target genes of hsa-miR-3616-3p by using the TargetScan database (http://www.targetscan.org/) and chosed 1654 target genes based on the criteria: Cumulative weighted context++ score<-2. Next, the functions of 1654 target genes were performed by GO and KEGG enrichment analysis. GO results showed that the target genes of hsa-miR-3616-3p signi cantly enriched in protein binding, integral component of plasma membrane and transcription factor activity, and so on (Figure5A). Moreover, KEGG analysis showed that the target genes of hsa-miR-3616-3p were enriched in focal adhesion, cell adhesion molecules (CAMs) and regulation of action cytoskeleton, and so on (Figure5B).

Bioinformatics Analysis of transcription factors of hsa-miR-3616-3p
Transcription factors (TFs) could bind to speci c sequences of DNA to exert transcription effect on one or more genes. The dysregulated miRNAs were controlled by transcription factors. By using the Transmir database (http://www.cuilab.cn/transmir), we found that 24 TFs were predicted to bind on the promoter domain of hsa-miR-3616-3p ( Figure 6A). KEGG enrichment analysis of 24 TFs was operated by DAVID and the results showed that these TFs of hsa-miR-3616-3p were enriched in transcriptional misregulation in cancer, TGF-beta signaling pathway and pathways in cancer, and so on ( Figure 6B). Moreover, the Transmir database revealed that 7 of the 24 TFs were associated with the breast tissues (Table 3). We validated and visualized the the expression of 7 TFs in Ualcan database (http://ualcan.path.uab.edu/) with TNBC tissue samples from TCGA. The results con rmed that ERG, MYC and NRF1 were downregulated (p<0.05) and GRHL2 and KDM5B were upregulated (p<0.05) in TNBC tissues compared with normal tissues, and CTCF and TRIM25 were no statistically signi cant (p>0.05) in TNBC tissues compared with normal tissues (Figure 6C-6I). These data demonstrated that hsa-miR-3616-3p could be regulated by these TFs (ERG, MYC,NRF1,,GRHL2 and KDM5B) through activating or inhibiting promoter domain.

Discussion
Worldwide, breast cancer is the most common cancer affecting women, and its incidence and mortality rates are expected to increase signi cantly the next years [22]. This highlights the pressing need for new biomarkers and novel mechanism for BC prediction and therapy. More and more studies have con rmed that miRNAs have already been the new biomarkers for diagnosis, prognosis, therapy prediction and therapeutic tools for breast cancer [23][24][25]. Meanwhile, miRNAs are key regulators in biological processes of BC, such as the occurrence and development of BC, proliferation signal transduction, resistance to cell death, activation of cell migration and invasion [26][27][28]. However, only a few miRNAs have been well characterized and little is known about the biological functions of miRNAs in BC. Excitingly, the high-throughput technology, public databases and bioinformatics methods greatly promote the understanding of miRNAs in BC. In this sduty, integrating multiple sets of BC miRNA microarrays to explore key miRNAs via bioinformatics analysis and public databases is a promising approach to identify the potential miRNA and provide new biomarkers for BC.
In our study, we analyzed two chip data from GSE57897 and GSE97811, and harvested 962 differentially expressed miRNAs between BC tumor and normal breast tissue samples, including 943 up-regulated miRNAs and 19 down-regulated miRNAs. The top 2 most signi cantly up-regulated miRNAs (hsa-miR-4428 and hsa-miR-575) and down-regulated miRNAs (hsa-miR-3909 and hsa-miR-3616-3p) were chosed to further experiments. Most of them have been reported that were associated with other cancers and diseases. Silencing of hsa-miR-4428 weaken the proliferative and migratory abilities of lung adenocarcinoma cells [29]. Hsa-miR-575 targeted BLID to promote growth and invasion of non-small cell lung cancer cells [30], moreover hsa-miR-575 exerted oncogenic function in hepatocellular carcinoma by increasing hepatocellular carcinoma proliferation and metastasis [31]. Downregulation of hsa-miR-3909 enhanced the in ammation level in the progression of rheumatic heart disease by acting on it's target gene IL1R1 [32]. However, hsa-miR-3616-3p hasn't been reported. Overall, it suggests that these 4 miRNAs in BC should be the important molecules.
Gene expression pro ling studies have shown that breast cancer is not caused by a single gene. Molecular subtyping of breast cancer can provide important prognostic information and provide a reference for the choice of treatment strategies [33]. We have mainly explored OS analysis of the 4 candidated miRNAs in the luminal, Her-2 overexpression and TNBC subtypes. TNBC constitutes the most aggressive molecular subtype among BC, usually appears in the form of high-grade invasive ductal carcinoma, with a higher rate of early recurrences, distant metastases, and poorer prognosis compared to other BC subtypes [34]. Our results showed that high expression of 4 candidated miRNAs had a worse prognosis in both BC and luminal subtype. In Her-2 overexpression subtype, the OS of hsa-miR-575 was only statistically signi cant, the main reason may be that the number of Her-2 overexpression samples was relatively small. In contrast, low expression of 4 candidated miRNAs had a worse prognosis in TNBC subtype. We further focused on the relationship between 4 candidated miRNAs and TNBC. Subsequently, we tested the expression of 4 candidated miRNAs by qRT-PCR on 15 TNBC and 15 non-TNBC tissue samples and found that the expression level of hsa-miR-3616-3p in TNBC tissue samples was consistent and stable with the microarray data. Further, it was showed that the down-expression hsa-miR-3616-3p could enhanced the metastasis and invasion ability of TNBC cells by wounding heal assays and transwell assays. In a word, hsa-miR-3616-3p would be a new diagnostic and prognostic marker or therapeutical target for TNBC patients.
MicroRNA (miRNAs) are small non-coding RNA (ncRNA) molecules that bind messenger RNA (mRNA) complementary sequences and generally direct post-transcriptional silencing in the 3′-untranslated region (UTR) of target genes [35]. To understand the relationship between hsa-miR-3616-3p and target genes, we screened out 1654 target genes based on the predictions and analyses of TargetScan database, and then we conducted functions of 1654 target genes by GO and KEGG enrichment analyses. The function of these target genes were enriched in focal adhesion, cell adhesion molecules (CAMs) and regulation of action cytoskeleton, and so on. Studies had shown that high focal adhesion kinase expression in breast carcinoma was associated with triple-negative phenotype [36] and focal adhesion kinase (FAK) inhibition prevented the migration of TNBC [37]. Previous researches had indicated that cell adhesion was an important pathway leading to TNBC poor prognosis [38,39]. And regulation of action cytoskeleton inhibited migration, invasion and adhesion of TNBC cells [40]. These results suggested that hsa-miR-3616-3p may bind target gene to affect biological behavior of TNBC via the above enrichment pathways.However, the target genes of hsa-miR-3616-3p were showing up in greater numbers, and it dose need further research to understand.
Studies had indicated that the expression of miRNA was regulated by transcription factors [41,42], so we predicted and analyzed transcription factors of hsa-miR-3616-3p by Transmir database. Transmir is a database of regulatory relationships between transcription factors and miRNA collated from the literature [43]. We found 24 transcription factors were related with hsa-miR-3616-3p, and the results of KEGG enrichment analysis of 24 TFs showed that these TFs of hsa-miR-3616-3p were enriched in transcriptional misregulation in cancer, TGF-beta signaling pathway and pathways in cancer, and so on, which were related to occurrence and development of BC. Moreover, 7 of the 24 TFs were associated with the breast tissues. The expression of these TFs was evaluated by Ualcan database. In TNBC tissues compared with normal tissues, the expression of ERG, MYC and NRF1 were downregulated (p<0.05), GRHL2 and KDM5B were upregulated (p<0.05). However, CTCF and TRIM25 were no statistically signi cant (p>0.05).. Studies had shown that high expressed ERG was relevant to prostate cancer including invasion and metastasis, EMT, epigenetic reprogramming, differentiation and in ammation [44]. MYC, high expressor in TNBC, which could inhibited miR-200b to promote the epithelial to mesenchymal transition and mammosphere formation of TNBC cells [45]. High expressed NRF1 in BC, showed unfavorable prognosis with a high risk of breast cancer death in white women [46]. It was reported that high GRHL2 expression was highly correlated with worse survival in TNBC, with regulation of EMT [47]. In previous studies, high expressed KDM5B enhanced TNBC cells metastasis and progression by regulating hsa-miR-448[48]. Together, it is possible that the expression of hsa-miR-3616-3p may be regulated by transcription factors.
The above conclusions indicated that target genes and TFs may regulate hsa-miR-3616-3p via these pathways affecting the progression of TNBC, which were worth further study. Together, unreported hsa-miR-3616-3p may play a key role in TNBC migration and invasion and may provide new ideas for the diagnosis, prognosis, and therapeutic targeting of TNBC.

Conclusion
In this study, based on database and bioinformatical analysis, we found out that down-expressed hsa-miR-3616-3p have a worse prognosis in TNBC. Moreover, down-expressed hsa-miR-3616-3p regulated the metastasis and invasion ability of TNBC cells by wounding heal assays and transwell assays. These results suggest that hsa-miR-3616-3p may be potential biomarker to diagnosis, prognosis, and therapy of TNBC.    Validation of four candidate miRNAs through qRT-PCR (A) The relative expression of hsa-miR-4428 was lower in TNBC tissues than that in non-TNBC tissues. (p=0.0371) (B) The relative expression of hsa-miR-575 was higher in TNBC tissues than that in non-TNBC tissues.
(p=0.0059) (C) The relative expression of hsa-miR-3909 was lower in TNBC tissues than that in non-TNBC tissues. (p=0.7897) (D) The relative expression of hsa-miR-3616-3p was lower in TNBC tissues than that in non-TNBC tissues. (p=0.0346) (E) Comparison of the expression of 4 candidate miRNAs identi ed by microarray and qRT-PCR. Our qRT-PCR datas showed that the expressions of hsa-miR-575, hsa-miR-3909 and hsa-miR-3616-3p were consistent with the microarray data.