The comparison of differentially expressed microRNAs in Bag-1 deficient and wild type MCF-7 breast cancer cells by small RNA sequencing

The multifunctional BAG-1 (Bcl-2 athanogene-1) protein promotes breast cancer survival through direct or indirect interaction partners. The number of the interacting partners determines its cellular role in different conditions. As well as interaction partner variability, the amount of BAG-1 protein in the cells could cause dramatic alterations. According to previous studies, while the transient silencing of Bag-1 enhanced drug-induced apoptosis, deletion of BAG-1 could induce stemness properties and Akt-mediated actin remodeling in MCF-7 breast cancer cells. Considering the heterogeneity of breast cancer and the variability of BAG-1 -mediated cell response, it has become essential to determine microRNA (miRNA) functions in breast cancer depending on Bag-1 expression level. This study aims to compare microRNA expression levels in wt and Bag-1 knockout (KO) MCF-7 breast cancer cells. hsa-miR-429 was selected as a potential miRNA in BAG-1KO MCF-7 cells because of the downregulation both in bioinformatics and validation qRT-PCR assay. According to predicted mRNA targets and functional enrichment analysis the ten hub proteins that are phosphatidylinositol-4,5-biphosphate 3-kinase catalytic subunit alpha (PIK3CA), kinase insert domain receptor (KDR), GRB2 associated binding protein 1 (GAB1), Rac family small GTPase1 (RAC1), vascular endothelial growth factor A (VEGFA), Cbl proto-oncogene (CBL), syndecan 2 (SDC2), phospholipase C gamma 1 (PLCG1), E1A binding protein p300 (EP300), and CRK like proto-oncogene, adaptor protein (CRKL) were identified as targets of hsa-miR-429. The functional enrichment analysis showed that the most significant proteins were enriched in PI3K/Akt, focal adhesion, cytoskeleton regulation, proteoglycans in cancer, and Ras signaling pathways. It was determined that hsa-miR-429 targeted these pathways in Bag-1 deficient conditions and could be used as a potential therapeutic target in future studies.


Introduction
The multifunctional Bag-1 protein regulates several cellular pathways, including apoptosis, autophagy, proliferation, invasion, and metastasis in cancer cells (Ozfiliz et al., 2015;Kilbas et al., 2019;Kizilboga et al., 2019;Turk et al., 2021). Bag-1 has three different isoforms (Bag-1S, 33 kDa; Bag1M 46 kDa; and Bag-1L 52 kDA) that have distinct antiapoptotic effects in breast cancer. Bag-1L isoform has a role in the increased cell viability and resistance to apoptosis (Liu et al., 2009). Moreover, according to the isoformspecific interactome analysis, protein processing, and ERassociated protein degradation (ERAD) pathways were controlled by Bag-1 expression in protein homeostasis of breast cancer (Can et al., 2021). It has been well implicated that increased Bag-1 expression is associated with breast cancer progression and invasion (Brimmell et al., 1999). The clinical relevance of Bag-1 was shown as a prognostic biomarker in Oncotype DX and PAM50 assay, which highlighted its role in breast cancer (Paik et al., 2004;Bernard et al., 2009). The immunohistochemical studies showed that Bag-1 cytoplasmic positive staining scores were correlated with breast cancer progression in estrogen receptor-positive tissue samples (Turner et al., 2001). However, downregulation of Bag-1 expression sensitized tamoxifen-resistant cells to tamoxifen and enhanced the apoptotic potential of cisplatin and paclitaxel in MCF-7 cells (Kilbas et al., 2019;S. Lu et al., 2019). Moreover, stable knockout of Bag-1 enhanced the stemness properties of cells through increased expression of stem cell markers in mice and caused stress-mediated hyperactivation of Akt and actin cytoskeleton remodeling in MCF-7 cells (Tang et al., 2017) (Ozfiliz-Kilbas P et al., n.d.). These findings highlighted that Bag-1 expression is critical in cellular homeostasis.
miRNAs are 22 nucleotide length endogenous, singlestranded small noncoding RNAs and negatively regulate gene expression and affect several cellular pathways. miRNAs regulate post-transcriptional mechanisms by degrading mRNA and inhibiting target protein translation through sequence-specific binding to the 3'-untranslated region (UTR) of target mRNA or coding sequences (Ghildiyal & Zamore, 2009;Hausser et al., 2013). miRNAs have essential biological functions in cell proliferation, apoptosis and differentiation, drug resistance, and mostly in cancer (Volinia et al., 2006). The tumor-associated miRNAs are classified as a tumor suppressor (TS)-miRNAs and cancer-causing onco-miRNAs (Takahashi et al., 2013). Dysregulation on the expression profiles of miRNAs have been commonly associated with cancer development, cellular transformation, and tumorigenesis (Melo & Esteller, 2011).
The association between heterogeneous characteristics of breast cancer and differential expression of miRNAs during breast cancer progression made it important to elucidate miRNA functions and the profiling of miRNAtargeted genes in breast cancer. Although numerous studies have shown that different expression levels of BAG-1 show different functions in breast cancer, miRNA profiles of Bag-1-deficient cells associated breast cancer development remains largely unknown. The current study aimed to screen the differential expressed miRNAs and determine the potential targeted genes and pathways in Bag-1 deficient conditions in MCF-7 breast cancer cells. High-throughput small RNA sequencing, identification of differentially expressed miRNAs, and prediction of target genes and networks in wt and BAG-1KO MCF-7 cells showed that PI3K/Akt signaling, ECM-receptor interaction, proteoglycans in cancer, and TGF-β signaling pathways were highly enriched in BAG-1KO samples compared to wt samples. Experimental qRT-PCR validation of differentially expressed miRNAs confirmed that hsa-miR-429 significantly decreased in BAG-1KO samples. The protein-protein interaction (PPI) network of potential gene targets of hsa-miR-429 showed that PI3K/Akt signaling and cell adhesion pathways are highly enriched in BAG-1KO cells when hsa-miR-429 is downregulated. Alternations in target miRNA expression due to directly or indirectly affected by Bag-1 deficiency provide novel results for the relationship between different Bag-1 expressions and miRNA profiles in breast cancer. The altered miRNA-associated potential gene and pathway network, obtained from this study, can be used as preliminary data for future studies.

Cell culture of wt and BAG-1KO MCF-7 cells
The wt MCF-7 (HTB-22) human breast cancer cells were maintained in Dulbecco's Modified Eagle Medium (DMEM, PanBiotech, Germany) augmented with 10% Fetal Bovine Serum (FBS, PanBiotech, Germany), 100-units/ml penicillin, and 100-µg/ml streptomycin (Pan Biotech, Germany). Cells were kept at 37°C in a humidified 5% CO 2 incubator. BAG-1KO MCF-7 cells were generated through co-transfection of Bag-1 CRISPR/Cas9 Knockout (KO) plasmid (sc-417179) (Santa Cruz Biotechnology, USA) and Bag-1 Homology Directed Repair (HDR) plasmid (sc-417179-HDR) (Santa Cruz Biotechnology, USA) into MCF-7 cells. Control CRISPR/Cas9 Knockout plasmid (sc-418922) was transfected into MCF-7 cells for negative control. Cells were transfected with UltraCruz Transfection Reagent (Santa Cruz Biotechnology, USA) for 48 h according to the manufacturer's protocol. Following transfection cells were treated with puromycin (Santa Cruz Biotechnology, USA) for 48h and selected by single-cell dilution process. Comparison of wt and Control KO cells showed similar survival, proliferation and growth potentials. After validation and characterization of Bag-1 deficiency in MCF-7 cells, one colony showing the lack of Bag-1 expression was selected for further experiments (Ozfiliz-Kilbas P et al., n.d.). The confirmation of knockout of Bag-1 in MCF-7 cells were shown in Figure S1. BAG-1KO MCF-7 cells were also grown in the same conditions with wt MCF-7 cells.

RNA extraction and quality control
According to the manufacturer's protocol, the total RNA of wt and BAG-1KO MCF-7 cells were isolated using SPLIT RNA Extraction Kit (Lexogen, Austria). According to the manufacturer's instructions, total RNA quality, concentration, and integrity were measured with the Agilent 2100 Bioanalyzer (Agilent Technologies, USA).

Construction of sequencing libraries
The total RNA from two biological replicates for each group was used in sequencing experiments. Small RNA libraries were constructed directly from total RNA using the Small RNA-Seq Library Prep Kit for Illumina Platforms (Lexogen, Austria) following the manufacturer's instructions. Small RNAs are ligated at 3, ' and 5' ends, respectively, and converted into cDNA.

RNA-seq and data analysis of small RNA sequencing
High-throughput RNA-seq was performed by Illumina MiSeq Kit V3 sequencing platform (Illumina, USA). The data were uploaded to the Galaxy web platform, and all downstream analyses were done at usegalaxy.org (Afgan et al., 2018). FastQC reported the quality check of raw reads with sequence quality scores, sequence-based quality scores, sequence GC and Kmer content, duplicate sequences, overrepresented sequences and overall sequence length distribution read distribution plots. The FastQC quality control reports from two biological replicates of wt and BAG-1KO MCF-7 cells were shown in Figure S2 A-D. The Illumina small RNA-Seq Library kit adapter sequence 5'-TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC-3' was trimmed, leaving a minimum overlap length=5 and maximum error rate= 0.1 by Cutadapt. The minimum length was selected 18 to discard sequences fewer than 18 base pairs in size, and the total length was selected 35 to exclude more than 35 base pairs in length. The quality of trimmed reads was controlled by FastQC and shown in Figure S3 A-D. The trimmed reads were aligned to the reference human genome GRCh37-hg19 by Bowtie2 with BAM format output. The aligned reads were annotated to mature human miRNA (http://www.mirbase.org/ftp. shtml) by featureCounts and were shown in Figure S4. The144 miRNAs were found in either in two replicated wt or BAG-1KO samples ( Figure S5), 88 miRNAs were only present in two replicated wt samples ( Figure S6), and 71 miRNAs were only present in two replicated BAG-1KO samples with at least 1 TPM (transcripts per million) ( Figure S7). The overall alterations between wt and BAG-1KO MCF-7 cells were analyzed by Rosalind OnRamp software ( Figure S8 A-B).

Identification of differentially expressed miRNAs in wt and BAG-1KO cells
To analyze differential expressions of 2576 miRNA counts obtained from featureCounts annotation, DESeq2 differential expression analysis was performed. The output of DESeq2 gives the fold changes (log2) of genes based on the factors (untreated/treated) and generates a tabular file. In our study, the fold changes of "Bag-1KO samples" against "wt" samples were compared, and the result was shown in Figure S9. The false discovery rates (FDR) ≤0.01 was accepted significantly. The fold change determined the expression ratio between wt and BAG-1KO cells. 2.6. The functional analysis of miRNAs using DIANA-miRPath v3.0 prediction algorithm DIANA-miRPath v.3.0 web platform was used to investigate the molecular function of miRNAs and enrichment of pathways. The algorithm performs a statistical model mentioned in Vlachos et al. (Vlachos et al., 2015) (http://snf-515788.vm.okeanos.grnet.gr/). Functional enrichment of target genes for downregulated and upregulated miRNAs was also performed on DIANA-miRPath v.3.0. KEGG pathway analysis was performed to predict target genes and molecular pathways associated with differentially expressed miRNAs in BAG-1KO and wt samples. GO analysis was performed to annotate the biological processes, cellular components, and molecular functions of miRNA-associated mRNA targets. Both functional enrichments analyzed FDR p-value <0.05 were enriched using Fisher's exact t-test hypergeometric distribution. 2.7. Validation of selected differentially expressed miRNAs by quantitative polymerase chain reaction (RT-qPCR) Total RNA, including smaller than ~200 nucleotides was extracted from MCF-7 wt and BAG-1KO cells using miRNeasy Mini Kit (Qiagen, USA), and the reverse transcription of total RNA containing miRNA was performed using miScript II RT Kit (Qiagen, USA) according to the manufacturer's protocol. The qRT-PCR experiment was performed using miScript SYBR Green PCR Kit (Qiagen, USA), and the cycling conditions were as follows: initial activation 95 °C for 15 min, followed by 40 cycles of denaturation 94 °C for 15 s, annealing 55 °C for 30 s, extension 70 °C for 30 s. The reactions were performed in two replicates and the selected miRNA primer sequences were listed in Table 1. RNU6 small nuclear RNA was used as a standard internal control and to normalize expression profiles. Relative expression of miRNA samples was calculated using the 2 -ΔΔCq methods. Data analyses were performed via GraphPad Prism 8.0.0.

Protein-Protein Interaction (PPI) Network Analysis
Search Tool for the Retrieval of Interacting Genes/ Proteins (STRING) database was used to predict proteinprotein interactions derived from genomic context, highthroughput lab experiments, co-expression, text mining database knowledge (Szklarczyk et al., 2019). The active interaction sources were selected from experiments and databases for reliable results. Moreover, the minimum required interaction score was set at high confidence (0.700). The mRNA regulatory network was visualized by Cytoscape open-source platform version 3. 8. 2 (Shannon et al., 2003).

Statistical analysis
The differential expression of miRNAs was performed using the Galaxy web platform, and the FDR ≤0.01 expression results were considered as significant. Then, the functional KEGG and GO enrichment analysis were performed using FDR <0.05. The qRT-PCR results were calculated using 2 -ΔΔCq and RNU6 was used as a normalization control. For the analyses of results "Two-way ANOVA followed by Sidak's multiple comparisons tests was performed using GraphPad Prism 8.0.0, Graphpad Software, San Diego, California USA, www.graphpad.com". To predict protein-protein interaction STRING was used at 0.7 high confidence ratio, and interaction network was visualized by Cytoscape web-platform using CytoHubba application in the Cytoscape database using the ˚the degree˚ calculation method. The Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics tools (https://david.ncifcrf.gov/) was performed for the functional characterization of mRNA targets.

Small RNA sequencing and determination of differentially expressed miRNAs in wt and BAG-1KO cells
High throughput next-generation sequencing was performed from two biological replicates of wt and BAG-1KO MCF-7 cells using the computational workflow in Figure 1A. Following sequencing, the quality of raw reads was controlled by FastQC. According to the basic quality report statistical report, the total number of raw sequences in wt and BAG-1KO samples was processed and the adapter sequence was removed by the Cutadapt tool. The trimmed reads with less than 18 bp were excluded for reliable analysis. The trimmed and processed reads of wt and BAG-1KO bases were sequentially aligned to the GRCh37-hg19 human genome using Bowtie2. The sequence counts for 2576 miRNAs of two biological replicates of wt and BAG-1KO samples annotated in the hg19 mature human miRNAs using the featureCounts tool. The differential expression analyses of 2576 miRNA counts were performed by DESeq2 and the mean average (MA) plot of differentially expressed miRNAs between BAG-1KO and wt MCF-7 samples were shown in Figure  1B. The Following the Benjamini-Hochberg correction method, False Discovery Rates (FDR) ≤0.01 was set to obtain significantly differentially expressed miRNAs, and with this correction method, 25 known miRNAs were found significantly differentially expressed in BAG-1KO reads comparing wt reads. According to the comparison results in BAG-1KO and wt MCF-7 cells, while 14 different miRNA targets were downregulated in BAG-1KO cells, 11 of the significantly upregulated miRNA targets were determined (Table 2).

Determination of predicted target genes of differentially expressed miRNAs in wt and BAG-1KO samples
To determine the potential functions of differentially expressed miRNAs and further explore miRNAsassociated target genes controlling molecular pathways, in silico DIANA-miRPath v3.0 miRNA target prediction algorithm was performed incorporating KEGG and GO functional enrichment analysis.
KEGG pathway enrichment analysis: We first explored the potential mRNA targets and significantly enriched pathways of all downregulated 14 miRNAs in Bag-1KO cells compared with wt MCF-7 cells. According to KEGG pathway analysis, we found that signaling pathways regulating pluripotency, PI3K/Akt signaling mechanism, Mucin type O-Glycan biosynthesis, ECMreceptor interaction, proteoglycans in cancer, and TGF-β signaling pathway were targeted in Bag-1KO MCF-7 cells ( Figure 2A). The upregulated 11 miRNAs were enriched in pathways involving morphine addiction, ECMreceptor interaction, mucin-type O-Glycan biosynthesis, and metabolism of xenobiotics by cytochrome p450, TGF-β signaling pathway, and prion diseases ( Figure 2B). Downregulated miRNA-associated 95 target genes as the highest number were involved in PI3K/Akt signaling pathway. Then, proteoglycans in cancer were found as the second enriched pathway with 57 target gene involvement.   44 genes were involved in the pluripotency of stem cell signaling pathways in Bag-1KO MCF-7 samples ( Figure  2C). However, it was observed that the number of enriched pathways of the upregulated miRNA-targeted genes was a small number compared to the downregulated miRNA-targeted genes enriched pathways. ECM-receptor interaction and TGF-β signaling pathways were the highly enriched pathways with the 18 and 16 genes associated with upregulated miRNAs ( Figure 2D). The enriched KEGG pathways, differentially expressed miRNAs, and miRNA-associated target genes of 14 downregulated and 11 upregulated miRNAs were shown in Table 3 and Table  4, respectively. GO functional annotation analysis: The significance of downregulated and upregulated miRNAs was analyzed using GO analysis in DIANA-miRPath v.3.0 by examining biological processes, molecular functions, and cellular components. The downregulated 14 miRNAs were enriched in biological processes involving gene expression, Fc-epsilon receptor signaling pathway, neurotrophin TRK receptor signaling pathway, cellular protein modification process, biosynthetic process, epidermal growth factor receptor signaling pathway, response to stress, cytoskeletal protein binding, extracellular matrix disassembly and organization and cell-cell signaling ( Figure 3A). The molecular functions of downregulated miRNAs were enriched in nucleic acid and protein binding transcription factor activity, ion binding, enzyme binding, cytoskeletal protein binding, enzyme regulator activity, RNA, and mRNA binding ( Figure 3B). The enriched cellular components include organelle, protein complex, nucleoplasm, and cytosol ( Figure 3C). The detailed GO enrichment analysis, including p-value and gene number of 14 downregulated miRNAs, was shown in Figure S10 A-C.
The biological process of upregulated 11 miRNAs was enriched in the biosynthetic process, cellular protein modification process, cellular nitrogen compound metabolic process, gene expression, phosphatidylinositolmediated signaling, mitotic cell cycle, and cell junction organization ( Figure 3D). The molecular function category of GO terms showed nucleic acid binding transcription factor activity, ion binding, protein binding transcription factor activity, and enzyme binding pathways of 11 upregulated miRNAs ( Figure 3E). The enriched cellular components include organelle, cytosol, protein   complex, and nucleoplasm ( Figure 3F). The detailed GO enrichment analysis, including p-value and gene number of 11 upregulated miRNAs, was shown in Figure S11A-S11C.

Selection of hsa-miR-429 as a potential target and analyzing the PPI Network of genes targeted by hsa-miR-429
According to the results of differentially expressed 25 miRNAs in Bag-1KO samples, hsa-miR-429 was found the third most downregulated expression profile of the 14-downregulated miRNAs after hsa-miR-21-3p and hsa-miR-21-5p (Table 2). It was also involved in 8 KEGG terms out of 13 enriched KEGG pathways of downregulated miRNA-associated target genes (Table  3), and experimentally validated with the confirmation of its downregulated expression profile (Figure 4). The results of hsa-miR-429 were found to be consistent in both sequencing and experimental PCR analysis, as they showed a similar downregulation profile. Therefore, hsa-miR-429 was selected as a potential target for Bag-1KO samples in further analyses.
The PPI network of the predicted gene targets for downregulated hsa-miR-429 was developed using the STRING database, and the visualization of the PPI network was established through Cytoscape software. The network was composed of 57 nodes and 97 edges ( Figure 5). The top ten hub proteins for the target of hsa-miR-429, which are the primary nodes of several interaction partners, were identified as phosphatidylinositol-4,5-biphosphate 3-kinase catalytic subunit alpha (PIK3CA), kinase insert domain receptor (KDR), GRB2 associated binding protein 1 (GAB1), Rac family small GTPase1 (RAC1), vascular endothelial growth factor A (VEGFA), Cbl protooncogene (CBL), syndecan 2 (SDC2), phospholipase C   gamma 1 (PLCG1), E1A binding protein p300 (EP300), and CRK like proto-oncogene, adaptor protein (CRKL). The PPI network of the top ten hub genes was shown in Figure 6 and Table 6.

Functional annotation and pathway enrichment analyses for the target gene of hsa-miR-429
The GO functional annotation and KEGG pathway enrichment analyses of the top ten hub genes included essential proteins of predicted target genes for hsa-miR-429, were performed using the DAVID database. The most significant GO terms were listed in Table  7, including positive regulation on gene expression, phosphatidylinositol-mediated signaling, angiogenesis, cell migration, and cell adhesion in the biological process category. Moreover, plasma membrane, cytoplasm, intracellular membrane, and focal adhesion were found in the cellular compartment category, protein binding, phosphatidylinositol-4,5-biphosphate3-kinase activity, protein kinase, receptor tyrosine kinase, and integrin binding were observed in the molecular function category. The KEGG pathway enrichment analysis revealed that the most significant genes targeted by hsa-miR-429 were enriched in proteoglycans in cancer, PI3K-Akt signaling, focal adhesion, cancer pathways, actin cytoskeleton regulation pathway, and Ras signaling pathway, which was listed in Table 8.

Discussion
BAG-1 protein is one of the critical targets in breast cancer cells, which functions as a co-chaperone for HSP70/HSC70 proteins, triggers the cellular survival and results in a high proliferation rate compared to normal cells (Song et al., 2001;Ozfiliz et al., 2015). Elevated protein levels of BAG-1 determined in metastatic cancer patients and associated with decreased drug-induced apoptosis in breast cancer (Kizilboga et al., 2019). On the contrary, elimination of BAG-1 expression through silencing increased the efficiency of therapeutics (Kilbas et al., 2019). Following this, BAG-1 expression levels were shown as a limiting factor in the efficiency of breast cancer therapies. miRNAs, the posttranscriptional regulators, are 18-23 nucleotide long noncoding RNAs (Castro, Moreira, Gouveia, Pozza, & Mello, 2017). They have regulatory roles controlling the development, progression, and metastasis of breast cancer at the posttranscriptional phase; however, these miRNAs' functions and expression profile differs according to the subtype of breast cancer (Kurozumi et al., 2017). Different miRNAs function independently to repress a single target, or combinatorial cause more dramatic effects in repressing a target (Linsen, Tops, & Cuppen, 2008). Although previous studies have demonstrated miRNA-mediated breast cancer development, the miRNA expression profile in Bag-1 deficient-breast cancer cells is unknown (Singh & Mo, 2013). One of the previous studies showed that the tumor suppressor hsa-miR-138 expression was decreased in the tumor tissues, and the regain of miR-138 expression enhanced apoptosis of gallbladder cancer through directly targeting the 3'-UTR sites of Bag-1. The loss and gain-of-function studies confirmed that the Bag-1 is a direct target for miR-138 . Moreover, miR-28-5p was found as a tumor suppressor directly  targeting and deregulating Bag-1 in B cells (Schneider et al., 2014).
Here, we discuss the differential expression profiles of miRNAs in Bag-1 deficiency leading to breast cancer progression in MCF-7 cells. Following small RNA sequencing analyses, 25 significant differentially expressed miRNAs were observed, 14 miRNAs were found downregulated, and 11 miRNAs were upregulated in BAG-1KO samples compared to wt MCF-7 samples. It was noteworthy that the enriched pathways associated with upregulated miRNAs were less in number compared with downregulated miRNAs enriched pathways. To eliminate the false-positive rates of prediction algorithms, the relative expression of 14 miRNAs selected out of 25 differentially expressed miRNAs to sample different log2 fold change number was experimentally validated by qRT-PCR analyses.
According to the results of downregulated miRNA profiles in BAG-1KO samples, previous studies showed that hsa-miR-21 is required for cell proliferation, PI3K/ Akt activation, and epithelial-mesenchymal transition (EMT) in MDA-MB-231 and SKM-1 cells (Li et al., 2018;Arisan et al., 2021). In addition, hsa-miR-27a-3p has been identified as an onco-miR and promotes proliferation and migration of colorectal and triple negative breast cancer cells through targeting BTG-1 and GSK3β, respectively (C. Su, Huang, Liu, Liu, & Cao, 2019). hsa-miR-342 and hsa-miR-146 are associated with estrogen levels, and have a role in negative regulation of BRCA1 in breast cancer (AI et al., 2011;Crippa et al., 2014). hsa-miR-345 has been known with its inhibitory role on tumor metastasis and induction of apoptosis in hepatocellular and pancreatic carcinoma (Bao et al., 2014;M et al., 2017). hsa-miR-151 was found to repress migration in breast cancer through  Heparan sulfate proteoglycan core protein 10 PLCG1 1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase gamma-1 10 EP300 Protein propionyltransferase p300 9 CRKL V-crk avian sarcoma virus CT10 oncogene homolog-like 9 Figure 6. PPI-network for the critical targets of hsa-miR-429. Circular yellow nodes represented the top 10 hub genes of predicted target genes for hsa-miR-429, grey nodes represented the significant targets of hsa-miR-429.
targeting TWIST1 (Yeh et al., 2016). In addition, hsa-miR-98-5p has been known for its anti-oncogenic function in breast cancer (Shi, Wang, Wang, & Gu, 2019). The miR-200 family is one of the most studied miRNAs associated with EMT and metastasis. Significant downregulation of the miR-200 family, which consists of miR-141, miR-200a, miR-200b, miR-200c, and miR-429, regulates the expressions of mesenchymal markers ZEB1, ZEB2, TWIST1, TWIST2, Snai1, and Snai2 and promotes EMT in breast cancer (Bockmeyer et al., 2011). Previous studies demonstrated that miR-200a showed decreased expression profile in EMT-induced cancer cells, but overexpression was associated with downregulated N-cadherin, ZEB1, vimentin, and upregulated E-cadherin profile (Y. Lu et al., 2014;Pichler et al., 2014). miR-429, a member of the miR-200 family, has been reported by its inhibitory role on EMT, tumor migration and metastasis in hepatocellular, nasopharyngeal and breast cancer (Guo, Zhao, Zhang, Liu, & Sun, 2018). Loss of function studies showed a decreased expression profile in breast cancer, inducing EMT and promoting bone metastasis directly targeting ZEB1 transcription factor (Ye et al., 2015). In the current study, hsa-miR-429 was consistently downregulated in BAG-1KO MCF-7 cells in both bioinformatics and qRT-PCR validation assays. Moreover, hsa-miR-429 was involved in downregulated miRNAsassociated KEGG enrichments in 8 of a total 13 pathways; therefore, we suggested that hsa-miR-429 and its gene targets can be a potential purpose to explain the survival properties of BAG-1KO cells. According to microRNA targets database (miRDB) and Targetscan, has-miR-429 significantly binds to 3'UTRs of Bag-5 and Bag-6, which Table 7. GO Functional Annotation analysis of the top ten hub genes included essential proteins of predicted target genes for hsa-miR-429. are also member of Bag-1 gene family, which possess antiapoptotic function in the cells due to presence of Bag domain to regulate Hsp70 co-chaperone activity. Although hsa-miR-429 is not directly bound to Bag-1 mRNA to regulate its function, it was observed from our results that the predicted gene targets of hsa-miR-429 were PIK3CA, KDR, GAB1, RAC1, VEGFA, CBL, SDC2, PLCG1, EP300, and CRKL, which are significantly enriched in PI3K/Akt signaling, cell migration, cell adhesion, focal adhesion, integrin binding, and actin cytoskeleton regulation. Although there is no direct evidence that miR-429 targets PIK3CA, previous studies showed that it inhibits cell proliferation and migration through targeting AKT serine/threonine kinase 1 (AKT1) involved in PI3K/ Akt pathway in renal cell carcinoma (Z. Su et al., 2020). Similarly, although direct interaction of KDR with hsa-miR-429 has not been demonstrated, it was observed that hsa-miR-429 suppresses breast cancer invasion by inhibiting the Wnt/β-catenin pathway through potentially targeting SOX2, CBL, FN1, KDR, KLF4, NR3C1, QK1, RBFOX3, SYNJ1, and VEGFA detected by CytoHubba algorithm in Cytoscape (Zhang et al., 2020). GAB1 was found to be directly targeted by hsa-miR-200a, a miRNA family including hsa-miR-429, to suppress cell invasion and migration of hepatocellular carcinoma (Wang & Sun, 2017). hsa-miR-429 decreased the expression of VEGFA in human endothelial cells by directly inhibiting Hypoxia-inducible factor 1 (HIF1) (Bartoszewska et al., 2015). In breast cancer, it was determined that hsa-miR-200bc/429 cluster directly targets PLCG1 and regulates EGF-driven invasion. (Uhlmann et al., 2010) Moreover, hsa-miR-429 also directly targets and negatively regulates (v-crk sarcoma virus CT10 oncogene homolog (avian)like) CRKL adaptor protein in breast cancer to promote cell migration and invasion (Guo et al., 2018). Although several studies determined the direct or indirect network between PIK3CA, KDR, GAB1, VEGFA, PLCG, and CRKL and hsa-miR-429 in various cancers, the functional regulation of these genes in Bag-1 deficiency conditions is still needed to elucidate. However, there is no direct or indirect interaction network between hsa-miR-429 and RAC1, CBL, SDC2 and EP300. Therefore, these preliminary screening data can be used to investigate the regulation of survival and proliferation conditions of Bag-1KO breast cancer cells.

Acknowledgment
The authors would like to thank Gizem Alkurt for valuable technical support for Illumina MiSeq's high throughput sequencing. We are thankful to DAPGenomics LLC for providing a small RNA sequencing library and supporting materials for this study.

Funding
This research was funded by İstanbul Technical University (Internal grant no: TDK-2017-40794) and DAPGenomics LLC. R&D grants.

Conflicts of interest
The authors declare no conflict of interest. Per base sequence quality Figure S2.  Per base sequence quality Figure S3.