Five miRNAs-mediated PIEZO2 downregulation, accompanied with activation of Hedgehog signaling pathway, predicts poor prognosis of breast cancer

Expression profile of PIEZO2 in human normal and cancer tissues

A high and detectable expression level of a gene is one of the most important traits for being a promising diagnostic or prognostic biomarker. Therefore, in the first place, we determined the expression of PIEZO2 in different normal tissues using the Human Protein Atlas (HPA) database. The results demonstrated that lung, gallbladder, urinary bladder, esophagus, cerebral cortex, prostate, spleen, seminal vesicle, smooth muscle and breast were the top ten normal tissues according to expression values of PIEZO2 mRNA (Figure 1A). The top ten tissues, sorted by expression levels of PIEZO2 protein, were adrenal gland, gallbladder, pancreas, stomach, small intestine, breast, parathyroid gland, appendix, lymph node and tonsil (Figure 1B). Figure 1A and Figure 1B together told us that gallbladder and breast were the two proper candidates for further investigation. The expression of PIEZO2 mRNA and protein in different types of cancer was successively analyzed using the HPA database (Figure 1C-F). Among all types of cancer, breast cancer presented as the highest expression value of PIEZO2 in both mRNA and protein levels. Taken these findings together, breast was selected for further investigation.

Figure 1. Expression of PIEZO2 in normal and cancer tissues from the HPA database. (A) PIEZO2 mRNA expression in different normal tissues; (B) PIEZO2 protein expression in different normal tissues; (C) PIEZO2 mRNA expression in different cancer tissues; (D) PIEZO2 protein expression in different cancer tissues (HPA031974); (E) PIEZO2 protein expression in different cancer tissues (HPA040616); (F) PIEZO2 protein expression in different cancer tissues (HPA015986).

PIEZO2 expression is frequently decreased in breast cancer and inversely correlates with progression

Next, we determined the expression of PIEZO2 in breast cancer cell lines and clinical samples compared with normal breast cell line and matched non-cancerous samples. Figure 2A showed that PIEZO2 expression in four breast cancer cell lines (MCF-7, Bcap37, MDA-MB-468 and MDA-MB-231) was significantly lower than that in normal breast cell line (HBL-100). Moreover, we found lower expression of PIEZO2 in high malignant cells (MDA-MB-468 and MDA-MB-231) compared with low malignant cell (MCF-7). Then, we compared the expression of PIEZO2 in breast cancer tissues with matched adjacent normal breast tissues and suggested that cancer tissues showed a significantly decreased PIEZO2 expression (Figure 2B). To further confirm under-expression of PIEZO2 in breast cancer, corresponding breast cancer expression data from TCGA were analyzed using UALCAN database (Figure 2C). Intriguingly, a similar result was observed. Furthermore, the protein level of PIEZO2 was also decreased in breast cancer tissues compared with normal breast tissues (Figure 2D).

Figure 2. Expression of PIEZO2 in breast cancer. (A) PIEZO2 expression in breast cancer cell lines (MCF-7, Bcap37, MDA-MB-468 and MDA-MB-231) compared with that in normal breast cell line (HBL-100); (B) PIEZO2 expression in clinical breast cancer tissues compared with that in matched adjacent normal tissues (n=16); (C) expression of PIEZO2 (also known as FAM38B) in breast cancer compared with normal controls by analyzing UALCAN database; (D) PIEZO2 protein expression level in breast cancer tissue and normal breast tissue was analyzed using immunohistochemical staining from HumanProteinAtlas database. *P<0.05. Errors bars indicate respective standard deviations.

Subsequently, we further studied the expression differences of PIEZO2 based on different clinicopathological parameters in breast cancer using bc-GenExMiner database. As shown in Figure 3A, no significant difference between the <51 years and >51 years was observed. Figure 3 also demonstrated that PIEZO2 expression was markedly upregulated in breast cancer with nodal negative (Figure 3B), ER positive (Figure 3C), PR positive (Figure 3D) and HER2 negative (Figure 3E) status. PIEZO2 expression was inversely correlated with NPI score and SBR grade as presented in Figure 3F and Figure 3G, respectively. Additionally, we also discovered that expression of PIEZO2 was significantly decreased in basal-like (Figure 3H) and triple negative (Figure 3I) breast cancer patients compared with not basal-like and not triple negative breast cancer, respectively. We further investigated the relationship between PIEZO2 expression and clinicopathological characteristics using TCGA breast cancer data. Chi-square test revealed that PIEZO2 expression was significantly associated with ER status (P < 0.001), PR status (P < 0.001), HER2 status (P < 0.001), T stage (P < 0.001), N stage (P = 0.0221) and pathologic stage (P = 0.0010) (Table 1). All these findings suggest that PIEZO2 expression in breast cancer is significantly declined and negatively correlates with progression of breast cancer.

Figure 3. Expression differences of PIEZO2 in breast cancer patients based on different clinicopathological features, including age (A), nodal status (B), ER status (C), PR status (D), HER2 status (E), NPI score (F), SBR grade (G), basal-like status (H) and triple-negative status (I).

Breast cancer patients with lower expression of PIEZO2 have poorer prognosis

Then, we explored the prognostic values of PIEZO2 in breast cancer. When we entered PIEZO2 in Kaplan Meier-Plotter database, two probes (1562488_at and 222908_at) were found. As shown in Figure 4, breast cancer patients with high expression of PIEZO2 had a significantly favorable prognosis, including overall survival, relapse free survival and distant metastasis free survival in both two probes. Regarding to post progression survival, high expression of PIEZO2 was also found to be significantly linked to a favorable prognosis in the probe 222908_at. In the probe 1562488_at, however, no statistical significance of PIEZO2 for predicting post progression survival was observed between the PIEZO2 high expression group and PIEZO2 low expression group. Furthermore, we also assessed the relationship between PIEZO2 expression and prognosis in breast cancer using PrognoScan database (Table 2). As listed in Table 2, PIEZO2 expression level was positively correlated with distant metastasis free survival, relapse free survival, disease specific survival, disease free survival and overall survival in patients with breast cancer.

Figure 4. Prognostic value of PIEZO2 (Affymetrix ID are valid: 1562488_at and 222908_at) in breast cancer patients. (A-B) Overall survival curve of PIEZO2 for breast cancer patients; (C-D) relapse free survival curve of PIEZO2 for breast cancer patients; (E-F) distant metastases free survival curve of PIEZO2 for breast cancer patients; (G-H) post progression survival curve of PIEZO2 for breast cancer patients.

The prognostic roles of PIEZO2 in breast cancer patients according to various clinicopathological features were also determined by Kaplan Meier-Plotter database. The data were presented in Figure 5 which demonstrated that PIEZO2 overexpression was significantly associated with better prognosis in ER-positive, HER2-negative, luminal A and luminal B breast cancer in the two probes. Furthermore, we also found that high expression of PIEZO2 in basal-like breast cancer patients indicated a better prognosis in the probe 1562488_at. In the probe 222908_at, a positive correlation of PIEZO2 expression and survival in PR-positive and lymph node-positive breast cancer patients but a negative association between PIEZO2 expression and survival in ER-negative, HER2-positive and pathological grade 3 breast cancer patients were observed.

Figure 5. Prognostic values of PIEZO2 in breast cancer patients based on different clinicopathological features. Green bars indicate a favorable prognosis; red bars indicate an unfavorable prognosis; black bars represent no statistical significance.

Identification of miRNAs that potentially regulate PIEZO2

To identify PIEZO2 that can be modulated by miRNAs, starBase database was utilized to predict upstream miRNAs of PIEZO2. In total, 182 miRNAs could potentially regulate PIEZO2 as shown in Table S1. It has been widely acknowledged that there exists an inverse association between miRNA and mRNA expression [12–15]. Therefore, we further assessed expression correlations between each predicted miRNA and PIEZO2 using TCGA breast cancer data. We observed that 42 miRNAs expression were statistically inversely associated with PIEZO2 expression (Table 3). Next, the prognostic values of the 42 miRNAs in breast cancer were determined by Kaplan Meier-Plotter database. The results demonstrated that breast cancer patients with high expression of 18 miRNAs (let-7d-5p, miR-33a-5p, miR-98-5p, miR-130b-3p, miR-137, miR-138-5p, miR-186-5p, miR-193a-3p, miR-196a-5p, miR-301a-3p, miR-421, miR-429, miR-454-3p, miR-579-3p, miR-580-3p, miR-671-5p, miR-934 and miR-2355-5p) indicated a poor prognosis (Figure 6). Moreover, we also detected the 18 miRNAs expression in breast cancer using OncomiR database. As listed in Table 4, only 11 miRNAs (let-7d-5p, miR-33a-5p, miR-98-5p, miR-130b-3p, miR-137, miR-193a-3p, miR-196a-5p, miR-301a-3p, miR-421, miR-454-3p and miR-671-5p) were significantly upregulated in breast cancer when compared with normal controls. Subsequently, by text mining for the 11 miRNAs, we found that 5 miRNAs (miR-130b-3p, miR-196a-5p, miR-301a-3p, miR-421 and miR-454-3p) were reported to function as oncogenes in breast cancer (Table 5). In order to further validate the inverse regulation of the five miRNAs in PIEZO2, we determined the mRNA and protein expression levels of PIEZO2 after knockdown of the five miRNAs. As shown in Figure S1, PIEZO2 was significantly upregulated in miRNAs-knockdown groups compared with negative control group. Taken together, the 5 miRNAs, with the features of upregulation in breast cancer, indicating poor prognosis, possessing negative associations with PIEZO2 expression and being reported to act as oncogenes, were the most potential upstream miRNAs that could inversely regulate PIEZO2 in breast cancer.

Figure 6. Prognostic values of 18 potential miRNAs in breast cancer analyzed by Kaplan Meier-plotter database.

GO functional annotation and pathway enrichment analysis of co-expressed genes of PIEZO2

Co-expression of PIEZO2 was analyzed using three databases, namely GEPIA, UALCAN and cBioPortal. As shown in Figure 7A, 109 co-expressed genes of PIEZO2 were commonly appeared in all the three databases. To better understand these genes, GO functional annotation and pathway enrichment analysis were conducted using Enrichr database. Three GO categories, containing biological process, cellular component and molecular function, were included in the functional annotation. For pathway enrichment, Reactome’s cell signaling pathway, KEGG’s cell signaling pathway and cell signaling pathway from NCI-Nature were analyzed.

Figure 7. GO functional annotation and pathway enrichment analysis for the co-expressed genes of PIEZO2. (A) The Veen diagram of PIEZO2’s co-expressed genes from GEPIA, UALCAN and cBioPortal databases; (B-D) GO functional annotation (biological process, cellular component and molecular function) for the 109 co-expressed genes of PIEZO2; (E-G) pathway (Reactome, KEGG and NCI-Nature) enrichment analysis for the 109 co-expressed genes of PIEZO2; (H) the Veen diagram of the genes enriched in Hedgehog signaling pathway from Reactome, KEGG and NCI-Nature databases.

The top ten enriched GO terms were shown in Figure 7B-D, including intraciliary retrograde transport, cilium assembly and plasma membrane bounded cell projection assembly in the biological process category, microtubule, polymetric cytoskeletal fiber and pre-autophagosomal structure in the cellular component category and ATP-dependent microtubule motor activity plus-end-directed, tRNA binding and ATP-dependent microtubule motor activity in the molecular function category. The top ten enriched pathways in three pathway categories were presented in Figure 7E-G. Intriguingly, we found that Hedgehog signaling pathway was the most top enriched pathway in all the three pathway databases. The corresponding gene counts of Hedgehog signaling pathway in Reactome, KEGG and NCI-Nature pathways were displayed in Figure 7H. A total of seven genes (IFT88, INTU, WDR19, WDR35, GLI3, CDON and CSNK1G3) were enriched in Hedgehog signaling pathway.

Decreased expression of PIEZO2 correlates with dysregulation of Hedgehog signaling pathway

To preliminarily explore the role of Hedgehog signaling pathway in PIEZO2-mediated progression of breast cancer, we first determined the expression of the seven genes enriched in Hedgehog signaling pathway in breast cancer using UALCAN database as shown in Figure 8A-G. Among the seven genes, an obvious upregulation of IFT88 (Figure 8A) in breast cancer was found whereas INTU (Figure 8B), WDR19 (Figure 8C), WDR35 (Figure 8D) and CDON (Figure 8F) expression levels in breast cancer samples were significantly lower than that in normal breast samples. GLI3 (Figure 8E) and CSNK1G3 (Figure 8G) showed no significant differences between cancer tissues and normal tissues. Subsequently, we further evaluated the prognostic values of the seven genes in breast cancer patients by Kaplan Meier-Plotter database (Figure 8H-N). The results showed that breast cancer patients with higher expression of CDON indicated a better prognosis. However, for the other six genes, no statistical significance was found. By combination of expression and prognostic roles of these genes, CDON was thought to closely correlate with PIEZO2 in breast cancer. The positive correlation of CDON expression with PIEZO2 expression in breast cancer was further determined by four databases, namely GEPIA, UALCAN, bc-GenExMiner and cBioPortal (Figure S2). Moreover, we also employed Oncomine database and Human Protein Atlas database to further testify the expression of CDON in breast cancer. As shown in Figure S3A-E, CDON mRNA expression was significantly lower in invasive breast carcinoma, invasive ductal breast carcinoma and invasive lobular breast carcinoma. Figure S3F demonstrated that expression level of CDON protein was also markedly decreased in breast cancer. Subsequently, CDON expression among groups of patients, based on various ER, PR, HER2, nodal, basal-like and triple negative status were determined using bc-GenExMiner. The analytic data were presented in Table 6. In breast cancer patients with positive ER, positive PR, negative HER2, negative nodal, negative basal-like and negative triple-negative status, expression of CDON was significantly upregulated when compared with corresponding counterparts. Furthermore, we also experimentally demonstrated the expression of CDON in breast cancer cell lines (Figure 8O) and clinical breast cancer samples (Figure 8P). The results also revealed that CDON expression was significantly downregulated in breast cancer cell lines and clinical cancer samples when compared with their counterparts. To preliminarily explore the upstream and downstream association of CDON and PIEZO2 in breast cancer, we detected the expression change of CDON or PIEZO2 after silencing expression of PIEZO2 or CDON using siRNA-PIEZO2 or siRNA-CDON, respectively. In this study, MCF-7 was chosen as the representative breast cancer cell line as its relatively high expression of PIEZO2 and CDON compared with other cell lines (Figure 2A and Figure 8O). The knockdown effects of siRNA-PIEZO2 and siRNA-CDON were presented in Figure 8Q and Figure 8S, respectively. Figure 8R showed a significant reduction of CDON expression after knockdown of PIEZO2. However, PIEZO2 expression was not statistically changed after downregulation of CDON (Figure 8T). All these findings indicate that PIEZO2 might be associated with Hedgehog signaling pathway by regulating CDON in breast cancer.

Figure 8. Identification of potential downstream of PIEZO2 in breast cancer. (A-G) Expression of IFT88, INTU, WDR19, WDR35, GLI3, CDON and CSNK1G3 in breast cancer analyzed using UALCAN; (H-N) prognostic roles of IFT88, INTU, WDR19, WDR35, GLI3, CDON and CSNK1G3 in breast cancer analyzed using Kaplan Meier-Plotter; (O) CDON expression in human breast cancer cell lines (MCF-7, Bcap37, MDA-MB-468 and MDA-MB-231) compared with that in normal breast cell line (HBL-100); (P) CDON expression in clinical breast cancer tissues compared with that in matched adjacent normal tissues (n=16); (Q) knockdown effect of siRNA-PIEZO2 in MCF-7 cell line; (R) expression change of CDON after silencing expression of PIEZO2 in MCF-7 cell line; (S) knockdown effect of siRNA-CDON in MCF-7 cell line; (T) expression change of PIEZO2 after silencing expression of CDON in MCF-7 cell line. *P<0.05; **P<0.01; “NS” represents no statistical significance.

Human Protein Atlas database analysis

Expression of PIEZO2 mRNA and protein in different human normal and cancer tissues were determined using the Human Protein Atlas (HPA) database [23]. CDON protein expression level in breast cancer tissues and normal breast tissues was also analyzed using HPA database.

UALCAN database analysis

UALCAN database, a user-friendly and interactive web resource, provides easy access to publicly available cancer transcriptome data from The Cancer Genome Atlas (TCGA) [24]. In this study, it was utilized to analyze gene expression and assess the correlation between two genes. Statistical analysis was calculated and displayed on the webpage and logrank P-value < 0.05 was considered as statistically significant.

Breast cancer gene expression miner

Breast cancer gene expression miner (bc-GenExMiner) is an easy-to-use online platform for analyzing gene expression, prognosis and correlation in breast cancer [25,26]. Bc-GenExMiner was introduced to determine PIEZO2 and CDON expression in breast cancer based on different clinicopathological features. Moreover, correlation of PIEZO2 and CDON was also assessed by bc-GenExMiner. P-value < 0.05 was considered as statistically significant.

Kaplan Meier-Plotter database analysis

Kaplan Meier-Plotter database is established using gene expression data and survival information of cancer patients downloaded from the Gene Expression Omnibus database [27]. The database was used to analyze associations between PIEZO2 expression and overall survival, relapse-free survival, distant metastases-free survival or post-progression survival in breast cancer. In this study, briefly, PIEZO2 was firstly entered into the database to obtain Kaplan Meier survival plots. PIEZO2 expression above or below the median classified these cases into a low expression group and a high expression group. These cohorts were then compared with a Kaplan-Meier survival plot, and hazard ratio (HR), 95% confidence interval (CI), and logrank P-value were determined and displayed on the webpage. In addition, the prognostic values of predicted miRNAs, IFT88, INTU, WDR19, WDR35, GLI3, CDON and CSNK1G3 in breast cancer were also evaluated using Kaplan Meier-Plotter database. A logrank P-value < 0.05 was considered as statistically significant.

PrognoScan database analysis

The correlation between the expression of PIEZO2 and survival in breast cancer was also determined using the PrognoScan database which is a database for meta-analysis of the prognostic values of genes in human cancers, including bladder cancer, blood cancer, brain cancer, breast cancer, colorectal cancer, esophagus cancer, head and neck cancer, lung cancer, ovarian cancer, prostate cancer, renal cell carcinoma, skin cancer and soft tissue cancer [28,29]. A Cox P-value < 0.05 was considered as statistically significant.

starBase database

The upstream miRNAs of PIEZO2 were predicted using starBase database, which is an open-source platform for studying the miRNA-ncRNA, miRNA-mRNA, ncRNA-RNA, RNA-RNA, RBP-ncRNA and RBP-mRNA interactions from CLIP-seq, degradome-seq and RNA-RNA interactome data [30,31]. starBase database was also employed to assess the inverse correlations of each miRNA expression and PIEZO2 expression. R < -0.1 and P-value < 0.05 were set as the thresholds for further identifying potential miRNAs that could target PIEZO2.

OncomiR database analysis

OncomiR database is an online resource for exploring pan-cancer microRNA dysregulation in cancer [32]. miRNA expression in breast cancer was determined using OncomiR database. P-value < 0.05 was considered as statistically significant.

GEPIA database analysis

GEPIA, a newly developed interactive web server for analyzing the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx projects, was used to identify the co-expressed genes of PIEZO2 in breast cancer [33]. The similar genes listed on the webpage were directly downloaded.

cBioPortal database analysis

cBioPortal is an online database for integrative analysis of complex cancer genomics and clinical profiles, which currently provides access to data from more than 48,668 tumor samples from 172 cancer studies in the TCGA pipeline [34,35]. It was also used to obtain the co-expressed genes of PIEZO2 in breast cancer. Only correlated genes with PCC > 0.3 were selected for subsequent investigation.

Enrichr database analysis

Enrichr, a comprehensive gene set enrichment analysis web server, was employed to conduct Gene Ontology (GO) functional annotation and pathway enrichment analysis for these commonly appeared co-expressed genes of PIEZO2 [36,37]. The top ten enriched GO items and pathways were displayed on the webpage and directly downloaded.

Oncomine analysis

Oncomine database is a web-based data mining platform for cancer research [38]. The expression of CDON in breast cancer was evaluated by Oncomine analysis of Cancer vs. Normal and meta-analysis. P-value < 0.0001 and fold change > 1.5 were set as the thresholds.

Cell culture

All cell lines (HBL-100, MCF-7, Bcap37, MDA-MB-468 and MDA-MB-231) used in this study were purchased from the cell bank of Chinese Scientific Academy (Shanghai, China). HBL-100, MDA-MB-468 and MDA-MB-231 were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Gibco, 12430047) supplemented with 10% fetal bovine serum (FBS; Biological Industries, 04-0101-1, Cromwell, CT, USA) and MCF-7 and Bcap37 were maintained in Roswell Park Memorial Institute (RPMI) 1640 medium (Gibco, 31800105, Life Technologies, Carlsbad) containing 10% FBS under a humidified atmosphere of 5% CO₂ at 37°C.

Patient and sample collection

16 clinical breast cancer tissues and matched adjacent normal tissues were collected from breast cancer patients who underwent surgery at the Zhejiang Cancer Hospital (Hangzhou, China) between 2016 and 2017. All procedures performed in this study involving human participants were conducted in accordance with the ethical standards of the Zhejiang Cancer Hospital and written informed consent from every participant was obtained.

Cell transfection

Cell transfection was performed as we previously described [12,13,39]. miRNA inhibitors, siRNAs and their negative control oligonucleotides (NC) were purchased from RiboBio Co. Ltd (Guangzhou, China). 2 x 10⁵ of MCF-7 cells were seeded onto six-well plates, and cultured for 12 hours under a humidified atmosphere of 5% CO₂ at 37°C. Subsequently, these cells were transfected with 50 nM of these oligonucleotides using Liopfectamine 3000 reagent (Invitrogen, Shanghai, China) according to the manufacturer’s instructions. The sequences of siRNA used in this work were listed in Table S2.

RNA isolation and quantitative real-time polymerase chain reaction

RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR) were performed as described previously [12]. RNAiso plus Reagent (TaKaRa biotechnology, 9109, Kusatsu, Japan) was utilized to extract total RNA from cell lines and clinical samples. Then, total RNA was reversely transcribed into complementary DNA (cDNA) by the PrimeScript^TM RT Reagent Kit (TaKaRa biotechnology, RR037A). Subsequently, the cDNA was used for real-time polymerase chain reaction (RT-PCR) analysis with gene-specific primers. Real-time PCR was performed in a Roche LightCycle480 II Real-Time PCR Detection System through SYBR Premix Ex Taq (TaKaRa biotechnology, RR420A). Glyceralddehyde-2-phosphate dehydrogenase (GAPDH) and U6 was used as the internal control for genes and miRNAs, respectively. The primers used in this study were listed in Table S2. The expression of PIEZO2 or miRNAs was normalized to GAPDH or U6, and calculated using the comparative threshold method (2^−ΔΔCT).

Western blot

Western blot assay was performed to determine PIEZO2 protein level as we previously described [40]. β-Actin was used for normalization of protein loading. Experiments were repeated at least three times.

Statistical analysis

All experiments were performed in triplicates. GraphPad prism 7 software (GraphPad Software, Inc., LaToIIa, CA, USA) was used to analyze expression data for statistical significance. The results were shown as mean ± SD. Differences between two groups were determined using Student’s t-test. The Chi-Square test was applied to the examination of relationship between PIEZO2 levels and clinicopathological features. A P-value < 0.05 was considered as statistically significant.