Pannexin1 Is Associated with Enhanced Epithelial-To-Mesenchymal Transition in Human Patient Breast Cancer Tissues and in Breast Cancer Cell Lines

Loss of connexin-mediated cell-cell communication is a hallmark of breast cancer progression. Pannexin1 (PANX1), a glycoprotein that shares structural and functional features with connexins and engages in cell communication with its environment, is highly expressed in breast cancer metastatic foci; however, PANX1 contribution to metastatic progression is still obscure. Here we report elevated expression of PANX1 in different breast cancer (BRCA) subtypes using RNA-seq data from The Cancer Genome Atlas (TCGA). The elevated PANX1 expression correlated with poorer outcomes in TCGA BRCA patients. In addition, gene set enrichment analysis (GSEA) revealed that epithelial-to-mesenchymal transition (EMT) pathway genes correlated positively with PANX1 expression. Pharmacological inhibition of PANX1, in MDA-MB-231 and MCF-7 breast cancer cells, or genetic ablation of PANX1, in MDA-MB-231 cells, reverted the EMT phenotype, as evidenced by decreased expression of EMT markers. In addition, PANX1 inhibition or genetic ablation decreased the invasiveness of MDA-MB-231 cells. Our results suggest PANX1 overexpression in breast cancer is associated with a shift towards an EMT phenotype, in silico and in vitro, attributing to it a tumor-promoting effect, with poorer clinical outcomes in breast cancer patients. This association offers a novel target for breast cancer therapy.


PANX1 Over-Expression Is Correlated With Poorer OS in Breast Cancer Patients
To investigate the role that PANX1 plays during cancer progression, we compared PANX1 expression between primary breast carcinoma (BRCA) using the RNA-seq dataset obtained from The Cancer Genome Atlas (TCGA). The TCGA data set contained 1180 breast tissue samples in total; 1170 breast cancer samples and 109 samples from breast tissue adjacent to the primary tumor that were considered normal in this study. PANX1 mRNA levels are significantly higher in breast cancer patient tissues than in the normal non-cancerous adjacent tissues ( Figure 1A).
Significantly higher PANX1 mRNA levels were seen in all of the intrinsic breast cancer subtypes when compared to normal breast cancer tissues of the TCGA data set ( Figure 1B). Compared to Luminal A (ER + PR + HER2 − ) breast cancer subtype, Luminal B (ER + PR + HER2 + ), TNBC and HER2-enriched subtypes showed significantly higher expression of PANX1. In fact, PANX1 was elevated in the different breast cancer subtypes not only at the transcriptional levels but also at the protein levels, as determined by Proteomics analysis of PANX1 protein levels in the intrinsic breast cancer subtypes ( Figure 1C). At the protein level, PANX1 had higher levels in HER2-enriched, TNBC, and Luminal B compared to Luminal A, which had the lowest PANX1 protein levels (p < 0.05 and p < 0.01) ( Figure 1C, upper panel). In addition, the levels of PANX1 protein and mRNA were correlated in the different intrinsic breast cancer subtypes (R = 0.34, p = 0.004) ( Figure 1C, lower panel).
Using qRT-PCR, we also investigated the expression of PANX1 in primary breast cancer tissues from a local cohort of archived breast cancer patients' samples. PANX1 mRNA levels were up-regulated in basal-like TNBC tissues (N = 11) and in HER2 − (N = 15) and HER2 + (N = 11) breast cancer subtypes, as compared to normal breast tissue obtained from subjects who underwent reduction mammoplasty; though statistical significance was only reached in the HER2subtype with p < 0.05 ( Figure 1D). These data indicate that PANX1 is upregulated, yet differentially in the different subtypes of breast cancer.
The elevated PANX1 expression in TCGA breast cancer tissues is correlated with clinical outcomes. Protein expression is presented as log2 Ratio while mRNA abundance is presented as log2 expression level. Linear regression alongside is 95% confidence interval (gray area). Pearson's correlation coefficient (R) and its p-value are denoted in the bottom-left corner. (D) qRT-PCR of PANX1 mRNA levels in a cohort of breast carcinoma patients of different subtypes: TNBC tissues N = 11; ER + PR − HER2 + N = 11; ER + PR + HER2 − N = 15. Patients were females with no prior therapy, selected according to the immune-histochemical tumor expression profile of ER, PR, and HER2. Normal breast tissues were obtained from breast tissue of patients who underwent reduction mammoplasty. (E) OS Kaplan Meier plots of the BRCA TCGA (left) and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC, right) breast cancer patients. The TCGA (N = 1068) and METABRIC (N = 1904) BRCA samples were divided into Low, Intermediate, or High PANX1 expression groups based on the 25th and 75th percentiles of PANX1 expression. Kaplan Meier plots were used to compare OS of High/Intermediate versus Low PANX1 expression groups. * p < 0.05, ** p < 0.01, and *** p < 0.001.

EMT Pathway Correlates Positively with PANX1 Expression
To gain a mechanistic insight into the effect of PANX1 overexpression in BRCA tissues, GSEA based on PANX1 expression in BRCA patients was run on the KEGG database and the gene ontology (GO) database. Three cell adhesion-related pathways, including adhaerens junction, focal adhesion, and gap junctions gene set, were among the highly enriched pathways in the KEGG database analysis (data not shown). GSEA analysis of the GO database revealed that the EMT pathway was one of the top enriched GO pathways, based on PANX1 expression (Figure 2A). Figure 2A also shows 16 highly enriched EMT genes that form the leading edge of the enrichment plot. In addition to their high protein and mRNA levels in the 72 samples of intrinsic breast cancer subtypes. Protein expression is presented as log2 Ratio while mRNA abundance is presented as log2 expression level. Linear regression alongside is 95% confidence interval (gray area). Pearson's correlation coefficient (R) and its p-value are denoted in the bottom-left corner. (D) qRT-PCR of PANX1 mRNA levels in a cohort of breast carcinoma patients of different subtypes: TNBC tissues N = 11; ER + PR − HER2 + N = 11; ER + PR + HER2 − N = 15. Patients were females with no prior therapy, selected according to the immune-histochemical tumor expression profile of ER, PR, and HER2. Normal breast tissues were obtained from breast tissue of patients who underwent reduction mammoplasty. (E) OS Kaplan Meier plots of the BRCA TCGA (left) and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC, right) breast cancer patients. The TCGA (N = 1068) and METABRIC (N = 1904) BRCA samples were divided into Low, Intermediate, or High PANX1 expression groups based on the 25th and 75th percentiles of PANX1 expression. Kaplan Meier plots were used to compare OS of High/Intermediate versus Low PANX1 expression groups. * p < 0.05, ** p < 0.01, and *** p < 0.001.

EMT Pathway Correlates Positively with PANX1 Expression
To gain a mechanistic insight into the effect of PANX1 overexpression in BRCA tissues, GSEA based on PANX1 expression in BRCA patients was run on the KEGG database and the gene ontology (GO) database. Three cell adhesion-related pathways, including adhaerens junction, focal adhesion, and gap junctions gene set, were among the highly enriched pathways in the KEGG database analysis (data not shown). GSEA analysis of the GO database revealed that the EMT pathway was one of the top enriched GO pathways, based on PANX1 expression (Figure 2A). Figure 2A also shows 16 highly enriched EMT Cancers 2019, 11, 1967 5 of 22 genes that form the leading edge of the enrichment plot. In addition to their high correlation with PANX1 expression, the 16 EMT genes of the leading edge are also highly inter-correlated ( Figure 2B). Therefore, PANX1 upregulation in the BRCA TCGA clinical samples correlated with alterations in the EMT pathway genes.
Cancers 2019, 11, x 5 of 23 correlation with PANX1 expression, the 16 EMT genes of the leading edge are also highly intercorrelated ( Figure 2B). Therefore, PANX1 upregulation in the BRCA TCGA clinical samples correlated with alterations in the EMT pathway genes. Notably, β-catenin (CTNNB1), known to regulate the EMT pathway, is among the genes highly correlated with PANX1 upregulation (Figure 2). This result shows that PANX1 overexpression correlates with upregulation of β-catenin expression.

PANX1 Channel Permeability Inhibition Reduces Cell Viability and Induces Cell Cycle Arrest in Breast Cancer Cell Lines
In order to validate the in silico data, we conducted in vitro assays, using luminal type MCF-7 and the highly invasive TNBC MDA-MB-231 cell lines. qRT-PCR and Western blotting analysis showed that PANX1 mRNA and protein are expressed in both MCF-7 ( Figure 3A,B) and MDA-MB-231 cell lines ( Figure 3E,F).
To assess the effect of PANX1 activity inhibition on PANX1 levels in MCF-7 and MDA-MB-231 cells, cells were treated with probenecid (PBN), a PANX1 channel inhibitor. Probenecid treatment of MCF-7 cells caused a significant decrease in the mRNA ( Figure 3A, p < 0.001) and protein levels ( Figure 3B, p < 0.01) of PANX1. Similar trends were observed in MDA-MB-231 cells ( Figure 3E,F, p < 0.05 and p < 0.01).
In addition, PANX1 channel blockade by PBN induced a significant dose-dependent decrease in cell proliferation, in both cell lines, as compared to vehicle-treated control cells ( Figure 3C,G), as early as 24 h following PBN treatment in MCF-7 and MDA-MB-231 cells. The reduction in proliferation reached a significant 20% decrease at 48 and 72 h following treatment with 0.5 mM PBN (P < 0.01) in both cell lines. A more pronounced decrease of 40% in cell proliferation was detected at 48 and 72 hours following treatment with 1 mM PBN, as compared to control cells in both cell lines ( Figure 3C, p < 0.001 and Figure 3G, p < 0.01).
In order to evaluate how PBN inhibits breast cancer cell proliferation, cell cycle analysis was performed at 72 hours following PBN treatment. Cell cycle distribution of MCF-7 and MDA-MB-231 Notably, β-catenin (CTNNB1), known to regulate the EMT pathway, is among the genes highly correlated with PANX1 upregulation (Figure 2). This result shows that PANX1 overexpression correlates with upregulation of β-catenin expression.

PANX1 Channel Permeability Inhibition Reduces Cell Viability and Induces Cell Cycle Arrest in Breast Cancer Cell Lines
In order to validate the in silico data, we conducted in vitro assays, using luminal type MCF-7 and the highly invasive TNBC MDA-MB-231 cell lines. qRT-PCR and Western blotting analysis showed that PANX1 mRNA and protein are expressed in both MCF-7 ( Figure 3A,B) and MDA-MB-231 cell lines ( Figure 3E,F).
To assess the effect of PANX1 activity inhibition on PANX1 levels in MCF-7 and MDA-MB-231 cells, cells were treated with probenecid (PBN), a PANX1 channel inhibitor. Probenecid treatment of MCF-7 cells caused a significant decrease in the mRNA ( Figure 3A, p < 0.001) and protein levels ( Figure 3B, p < 0.01) of PANX1. Similar trends were observed in MDA-MB-231 cells ( Figure 3E,F, p < 0.05 and p < 0.01).
In addition, PANX1 channel blockade by PBN induced a significant dose-dependent decrease in cell proliferation, in both cell lines, as compared to vehicle-treated control cells ( Figure 3C,G), as early as 24 h following PBN treatment in MCF-7 and MDA-MB-231 cells. The reduction in proliferation reached a significant 20% decrease at 48 and 72 h following treatment with 0.5 mM PBN (P < 0.01) in both cell lines. A more pronounced decrease of 40% in cell proliferation was detected at 48 and 72 h following treatment with 1 mM PBN, as compared to control cells in both cell lines ( Figure 3C, p < 0.001 and Figure 3G, p < 0.01).
In order to evaluate how PBN inhibits breast cancer cell proliferation, cell cycle analysis was performed at 72 h following PBN treatment. Cell cycle distribution of MCF-7 and MDA-MB-231 cells showed that 0.5 or 1 mM PBN treatment induced cell cycle arrest in the G0/G1 phase ( Figure 3D,H), preventing cells from progressing towards mitosis. In fact, a higher percentage of cells were blocked in the G0/G1 phase following 72 h of PBN treatment; MCF-7 cells-39% of cells were arrested at G0/G1 when treated with 0.5 mM PBN (p < 0.01) and 41% when treated with 1 mM PBN (p < 0.001) versus 26% of control cells ( Figure 3D); MDA-MB-231 cells-53% of cells were arrested at G0/G1 when treated with 0.5 mM PBN (p < 0.05) and 57% when treated with 1 mM PBN (p < 0.05) versus 38% in control cells ( Figure 3H); in addition to a reduced percentage of cells in the G2/M phase in PBN-treated cells compared to control. cells showed that 0.5 or 1 mM PBN treatment induced cell cycle arrest in the G0/G1 phase ( Figure  3D,H), preventing cells from progressing towards mitosis. In fact, a higher percentage of cells were blocked in the G0/G1 phase following 72 hours of PBN treatment; MCF-7 cells-39% of cells were arrested at G0/G1 when treated with 0.5 mM PBN (p < 0.01) and 41% when treated with 1 mM PBN (p < 0.001) versus 26% of control cells ( Figure 3D); MDA-MB-231 cells-53% of cells were arrested at G0/G1 when treated with 0.5 mM PBN (p < 0.05) and 57% when treated with 1 mM PBN (p < 0.05) versus 38% in control cells ( Figure 3H); in addition to a reduced percentage of cells in the G2/M phase in PBN-treated cells compared to control. These data indicate that PBN has decreased cell growth and viability of both MCF-7 and MDA-MB-231 cells by inducing a cell cycle arrest at the G0/G1 phase.

Pharmacological Inhibition or Genetic Ablation of PANX1 Channels Abrogate Ethidium Bromide (EtBr) Dye Uptake by PANX1 Channels
In addition to pharmacological blockade of PANX1 channels, genetic knock-out of PANX1 expression was also performed. CRISPR/Cas9 gene editing system was employed to generate MDA- These data indicate that PBN has decreased cell growth and viability of both MCF-7 and MDA-MB-231 cells by inducing a cell cycle arrest at the G0/G1 phase.

Pharmacological Inhibition or Genetic Ablation of PANX1 Channels Abrogate Ethidium Bromide (EtBr) Dye Uptake by PANX1 Channels
In addition to pharmacological blockade of PANX1 channels, genetic knock-out of PANX1 expression was also performed. CRISPR/Cas9 gene editing system was employed to generate MDA-MB-231 cells with knocked-out PANX1 (MDA-PANX1cells). MDA-PANX1cells had diminished levels of PANX1 mRNA, as shown by qRT-PCR ( Figure 4A), and protein, as shown by Western blotting ( Figure 4B) and immunofluorescence ( Figure 4C).
The probenecid mode of action is thought to be mediated at least in part through the inhibition of PANX1 channels [48]. The dye uptake assay was performed to evaluate the ability of cells to uptake EtBr through PANX1 channels. PBN-treated and control MDA-MB-231 cells were incubated with EtBr. Control cells were able to uptake EtBr, while PBN-treated MDA-MB-231 cells displayed a significant reduction of EtBr uptake as determined by lower EtBr mean fluorescence intensity (MFI; Figure 4D). These data confirm effective inhibition of EtBr dye uptake. The knocked-out PANX1 expression was further confirmed by the lack of EtBr uptake in MDA-PANX1 -cells as revealed by dye uptake assay ( Figure 4D). The decrease in PANX1 channel activity in MDA-PANX1 -cells was similar to that exhibited by PBN-treated MDA-MB-231 cells ( Figure 4D). The EtBr dye uptake was also evaluated in the presence of normal divalent ion solution (2 mM CaCl2 and 1 mM MgCl2), at room temperature and in the presence of 1 mM ATP to activate EtBr uptake by PANX1 [49]. PBN effectively inhibited ATP-induced EtBr dye uptake ( Figure S2).  The probenecid mode of action is thought to be mediated at least in part through the inhibition of PANX1 channels [48]. The dye uptake assay was performed to evaluate the ability of cells to uptake EtBr through PANX1 channels. PBN-treated and control MDA-MB-231 cells were incubated with EtBr. Control cells were able to uptake EtBr, while PBN-treated MDA-MB-231 cells displayed a significant reduction of EtBr uptake as determined by lower EtBr mean fluorescence intensity (MFI; Figure 4D). These data confirm effective inhibition of EtBr dye uptake. The knocked-out PANX1 expression was further confirmed by the lack of EtBr uptake in MDA-PANX1cells as revealed by dye uptake assay ( Figure 4D). The decrease in PANX1 channel activity in MDA-PANX1cells was similar to that exhibited by PBN-treated MDA-MB-231 cells ( Figure 4D). The EtBr dye uptake was also evaluated in the presence of normal divalent ion solution (2 mM CaCl 2 and 1 mM MgCl 2 ), at room temperature and in the presence of 1 mM ATP to activate EtBr uptake by PANX1 [49]. PBN effectively inhibited ATP-induced EtBr dye uptake ( Figure S2).

Pharmacological Inhibition or Genetic Ablation of PANX1 Channels Reverse EMT in Breast Cancer Cells
To validate the EMT data obtained in the BRCA TCGA samples in silico, key EMT markers' expression levels were measured by qRT-PCR in MCF-7 and MDA-MB-231, in the presence or absence of PBN, and in the MDA-PANX1cells ( Figure 5A  In addition, mRNA levels of key regulators of EMT that included HIF-1α, Snail, Slug, TGF-1β, and β-catenin were assessed. HIF-1α mRNA levels were significantly decreased in MDA-PANX1cells (p < 0.05, Figure 5B) and in PBN-treated MCF-7 cells (p < 0.001, Figure 5A); while PBN exposure did not seem to affect HIF-1α mRNA levels in MDA-MB-231 cells ( Figure 5B). Furthermore, a significant decrease in mRNA levels of TGF-1β (p < 0.05) and Slug (p < 0.01) were exhibited by MCF-7 cells upon PBN treatment ( Figure 5A). Probenecid-treated MDA-MB-231 cells had a non-significant decrease in Snail mRNA levels, while MDA-PANX1 -cells expressed significantly lower levels of Snail In addition, mRNA levels of key regulators of EMT that included HIF-1α, Snail, Slug, TGF-1β, and β-catenin were assessed. HIF-1α mRNA levels were significantly decreased in MDA-PANX1cells (p < 0.05, Figure 5B) and in PBN-treated MCF-7 cells (p < 0.001, Figure 5A); while PBN exposure did not seem to affect HIF-1α mRNA levels in MDA-MB-231 cells ( Figure 5B). Furthermore, a significant decrease in mRNA levels of TGF-1β (p < 0.05) and Slug (emphp < 0.01) were exhibited by MCF-7 cells upon PBN treatment ( Figure 5A). Probenecid-treated MDA-MB-231 cells had a non-significant decrease in Snail mRNA levels, while MDA-PANX1cells expressed significantly lower levels of Snail (p < 0.01). Notably, β-catenin transcript levels were significantly downregulated in PBN-treated MCF-7 cells (p < 0.01), and in MDA-PANX1 -(p < 0.001). Figure 6 further supports a mesenchymal-epithelial transition (MET) induced by PANX1 downregulation. Immunostaining for E-cadherin and N-cadherin showed enhanced E-cadherin protein levels and reduced N-cadherin protein levels upon PANX1 downregulation, in accordance with the transcriptional data. Furthermore, Figure 6A (right panel) shows that PANX1 downregulation led to more prominent cytoplasmic/membranous distribution and reduced nuclear localization of β-catenin. Overall, pharmacological or genetic downregulation of PANX1 in breast cancer cells reverses the already existing EMT, i.e., promote a MET state in those cells.

Pharmacological Inhibition or Genetic Ablation of PANX1 Channels Reduce the Metastatic Potential of Breast Cancer Cells
Epithelial-to-mesenchymal transition is inherently correlated with the increased metastatic potential of cancer cells. Reversal of EMT upon PBN treatment prompted us to assess the effect of functional inhibition of PANX1 channels on the metastatic potential of breast cancer cells. Given their major role in invasion and metastasis, enzymatic activity of the matrix metalloproteinases MMP-2 and MMP-9 was measured by gelatin zymography in MCF-7 and in MDA-MB-231 cells, in the Overall, pharmacological or genetic downregulation of PANX1 in breast cancer cells reverses the already existing EMT, i.e., promote a MET state in those cells.

Pharmacological Inhibition or Genetic Ablation of PANX1 Channels Reduce the Metastatic Potential of Breast Cancer Cells
Epithelial-to-mesenchymal transition is inherently correlated with the increased metastatic potential of cancer cells. Reversal of EMT upon PBN treatment prompted us to assess the effect of functional inhibition of PANX1 channels on the metastatic potential of breast cancer cells. Given their major role in invasion and metastasis, enzymatic activity of the matrix metalloproteinases MMP-2 and MMP-9 was measured by gelatin zymography in MCF-7 and in MDA-MB-231 cells, in the presence or absence of PBN treatment and in the MDA-PANX1 − cells. Data showed a prominent and significant decrease in both proMMP9 (p < 0.05 and p < 0.01) and MMP-9 (p < 0.05 and p < 0.01) activity in MDA-MB-231 and in MDA-PANX1cells ( Figure 7C,D). MMP2, on the other hand, showed a slight significant activity decrease in the PBN-treated MDA-MB-231 cell, only (p < 0.05). MCF-7 cells, treated with PBN, showed a significant decrease in the activity of MMP2, only ( Figure 7A,B, p < 0.05).  Decreased metastatic potential of breast cancer cells with diminished PANX1 activity prompted us to assess, using Real-Time Cell Analysis (RTCA), the invasive capacity of MDA-MB-231 cells, upon PBN treatment, and that of MDA-PANX1 -cells. MDA-MB-231 cells treated with PBN showed a significant decrease in their invasiveness (80% decrease compared to control, p < 0.01) and their proliferation rate (50% decrease compared to control, p < 0.01) ( Figure 7E). These findings suggest that PANX1 channel permeability blockade or gene deletion decrease the metastatic potential of breast cancer cells MDA-MB-231 cells by decreasing their extracellular matrix penetrating ability and their overall invasive potential.

Discussion
Connexins mediate direct cell-cell communication [50][51][52] and maintain cell and tissue homeostasis [53,54]. Loss of connexins is a hallmark of cancer progression, including breast cancer [55][56][57]. Pannexins, a vertebrate glycoprotein family [17], share functional and structural features with connexins and also mediate communication between the cell and its microenvironment [20,25]. PANX1, the best characterized member of the pannexin family, plays important roles in tissue Decreased metastatic potential of breast cancer cells with diminished PANX1 activity prompted us to assess, using Real-Time Cell Analysis (RTCA), the invasive capacity of MDA-MB-231 cells, upon PBN treatment, and that of MDA-PANX1cells. MDA-MB-231 cells treated with PBN showed a significant decrease in their invasiveness (80% decrease compared to control, p < 0.01) and their proliferation rate (50% decrease compared to control, p < 0.01) ( Figure 7E). These findings suggest that PANX1 channel permeability blockade or gene deletion decrease the metastatic potential of breast cancer cells MDA-MB-231 cells by decreasing their extracellular matrix penetrating ability and their overall invasive potential.

Discussion
Connexins mediate direct cell-cell communication [50][51][52] and maintain cell and tissue homeostasis [53,54]. Loss of connexins is a hallmark of cancer progression, including breast cancer [55][56][57]. Pannexins, a vertebrate glycoprotein family [17], share functional and structural features with connexins and also mediate communication between the cell and its microenvironment [20,25]. PANX1, the best characterized member of the pannexin family, plays important roles in tissue development and differentiation [17,[58][59][60][61] and is implicated in several human diseases [25]. PANX1 expression has been reported during mouse mammary gland development and in the adult mouse mammary gland [31]; however, its role in breast cancer is still evolving [62].
Our results uncovered a role for PANX1 in breast cancer metastasis through EMT regulation. We provide evidence that PANX1 upregulation is correlated with poor prognosis in breast cancer patients and that it is differentially expressed in the different "intrinsic" human breast cancer subtypes. This observation is supported in a previous report by Stewart et al. who correlated PANX1 with negative clinical outcomes leading to poor overall survival and distant metastasis free survival in patients with breast cancer using in silico arrays [31]. Furthermore, Furlow et al. have shown that PANX1 is enriched at metastatic breast cancer foci without any elucidation into the underlying mechanism [30]. Here, we have provided data on pathways and biological processes that are enriched upon PANX1 increased expression. EMT pathway genes, were among the highly enriched gene ontology (GO) pathways following GSEA.
To support our in silico findings we used two human breast cancer cell lines, MCF-7 and MDA-MB-231, these are widely used cell line models to study metastasis [63]. We assessed the effect of PANX1 channel on breast cancer progression via pharmacological inhibition of PANX1 in these cells. Probenecid (PBN), a drug that has been used for gout treatment, is a potent PANX1 inhibitor [48,64,65] and is among several compounds that have been reported to effectively block PANX1 channel activity [66]. Treating MCF-7 and MDA-MB-231 cells with PBN not only decreased PANX1 activity and expression but also significantly reduced the proliferation of these cells by inducing a cell cycle arrest at the G0/G1 transition. Effectively, PBN reduced the number of cycling cells progressing through DNA replication and mitosis. Previously, similar effects on cell viability and PANX1 expression were evidenced in MDA-MB-231 cells [67] and in hippocampal cells from aged rats [68], respectively. In the hippocampal cells, there was a marked reduction in PANX1 protein levels following PBN treatment, similar to the results obtained in this study. The PBN-induced cell cycle arrest can explain its recent use as an effective adjuvant therapy to sensitize breast cancer cells and enhance the efficacy of bisphosphonate chemotherapy [67].
Our data show that the inhibition of PANX1 channel activity, using PBN, or the deletion of the PANX1 gene in breast cancer cells, did decrease the mRNA levels of several EMT genes, thus confirming the GSEA in silico data. Importantly, inhibiting PANX1 channels or PANX1 deletion upregulated E-cadherin (an epithelial marker) and down-regulated N-cadherin (a mesenchymal marker). PANX1 upregulation supports an EMT phenotype, while its inhibition or downregulation promotes a mesenchymal to epithelial transition. EMT [34] is crucial in the metastatic process [36] and is associated with tumor cell migration and invasion into the surrounding stroma and dissemination into secondary organs [35,37]. EMT is, in part, mediated through the cell adhesion molecule E-cadherin [38,39]. E-cadherin is crucial in maintaining the epithelial cell phenotype and its downregulation has been correlated with increased invasiveness of breast cancer cells [69]. In addition, the expression of the non-epithelial and mesenchymal-associated molecule N-cadherin contributes to increased invasiveness and motility of breast cancer cells [69,70]. A set of pleiotropically acting transcription factors, including Snail, Slug, and other regulators, orchestrate EMT. Our study revealed a decrease in Snail and Slug mRNA levels in TCGA samples and upon PANX1 downregulation in MCF-7 and in MDA-MB-231 cells. These transcription factors regulate EMT in different types of malignant tumors and their invasiveness. This was confirmed by studies that revealed that EMT transcription factors such as Snail could elicit metastasis when ectopically overexpressed [70][71][72]. We also reported a significant decrease in TGF-1β expression in MCF-7 cells. TGF-1β, known for its anti-proliferative effects and immune evasion by cancer cells [73,74], activates EMT [75,76]. In addition, PANX1 has a role in tumor hypoxia [77]. PANX1 downregulation led to a reduction in the mRNA levels of the hypoxia transcription factor and marker, HIF-1α in breast cancer cells. Hypoxia stabilizes several transcription factors that promote EMT progression and cancer cell metastasis [78,79]. Taken together, we show that ablating the PANX1 gene or inhibiting its activity reverts EMT. This result confirms the GSEA data which demonstrated that PANX1 upregulation correlates with an enhanced EMT state. Similar results have been recently reported in testicular cancer cells, where downregulation of PANX1 or its inhibition downregulated vimentin protein levels and upregulated E-cadherin protein levels, through signal-regulated kinase (ERK) [80].
GSEA showed that PANX1 upregulation correlated with the expression of several cell adhesion molecules including β-catenin. Malignant and invasive breast tumors are associated with mutations and over-expression of β-catenin. During EMT, nuclear localization of β-catenin activates target genes, a group of which are EMT genes [41,44]. The β-catenin/TCF/LEF transcription complex directly upregulates genes associated with EMT, such as EMT transcription factors Snail1 and Zeb1 [41,44,45]. Snail1 and Zeb1 act as transcription repressors of E-cadherin, leading to down-regulation of E-cadherin expression during EMT [46,47]. Of these established EMT markers, our study revealed that Snail and Slug transcription factors highly correlate with PANX1 expression. PANX1 could possibly act indirectly by regulating genes such as β-catenin which then regulate EMT. On the other hand, a direct role for PANX1 as a transcription regulator of EMT in the nucleus is plausible. This direct role for PANX1 is not far-fetched since such a role has been reported for Connexin 43 [81], a protein that shares structural and functional similarities with PANX1. However, this mechanism requires further scrutiny.
Interestingly, the decreased expression of EMT key markers, upon PANX1 inhibition or genetic deletion, was also correlated with decreased matrix metalloproteinases (MMPs) activity in MCF-7 and MDA-MB-231 cells, and decreased MDA-MB-231 cells invasion potential. MMPs play a vital role in the cancer microenvironment and are markers of breast cancer patients' poor prognosis [82]. Reduced MMP activity and invasiveness are common characteristics of less metastatic forms of cancer [82]. Overall, a pathway of PANX1 action could be envisaged to explain why PANX1 overexpression correlates with poor prognosis in breast cancer patients. In this pathway, upregulated PANX1 enhances breast cancer cells EMT phenotype, the enhanced EMT phenotype enhances breast cancer cells metastatic potential and invasiveness, for example by increasing the levels or activity of MMPs and other metastatic genes. This scenario is supported by data from recent studies on testicular cancer cells, where PANX1 inhibition downregulates EMT genes, upregulates E-cadherin, downregulates MMP-9, and decreases cell invasion [80]. Freeman et al. showed that decreasing PANX1 protein levels, in melanoma cells, using shRNA or blocking PANX1 channel using PBN or carbenoxolone decreased cell growth, migration, and invasiveness, probably through the Wnt/β-catenin pathway [83].
Together, our findings identify a role for PANX1 in breast cancer progression and that PANX1 mediates this tumor-promoting role by modification of the EMT pathway.
Probenecid (PBN, Santa Cruz Biotechnology, Dallas, TX, USA), also known as p-(dipropylsulfamoyl) benzoic acid, is a PANX1 channel inhibitor [48]. PBN was dissolved in DMSO and DMSO-treated cells were used as vehicle control.

CRISPR/Cas9 Mediated Targeting of Human PANX1
To knock out PANX1 gene from MDA-MB-231 cells, gRNA oligo pairs targeting exon 2 of human PANX1 were designed with high specificity and low off-targets using the CRISPOR software (www.crispor.tefor.net) and inserted into the plasmid pX330-puro-Rosa26-H3F3B, (Addgene plasmid # 73131; a kind gift from Dr. Agnel Sfeir [86]). MDA-MB-231 cells at 70-80% confluence were transfected with the vector containing the PANX1 gRNA, using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). Limiting dilution cloning was performed by trypsinizing the cells and seeding them at a density of 0.8 cell per well of a 96-well plate. After 10 days of selection in media containing 1 µg/mL puromycin, growing colonies were isolated and transferred to larger culture plates for further expansion. The selected and expanded colonies were further grown and harvested for RNA and protein analysis. Knock-out of PANX1 expression was confirmed by qRT-PCR, Western blotting, and immunofluorescence. MDA-MB-231 cells with knocked-out PANX1 are designated as MDA-PANX1cells. Mock-transfected MDA-MB-231 cells were generated by transfection of an empty pX330-puro-Rosa26-H3F3B plasmid without any targeting gRNAs into MDA-MB-231 cells followed by selection in puromycin as described above and these cells were used as a control to MDA-PANX1cells.

Transcriptomic Datasets
BRCA RNAseq data were acquired from The Cancer Genome Atlas (TCGA) database [87]. Gene expression values were represented as Fragment Per Kilobase Million mapped reads (FPKM). Prior to analysis, FPKM values were transformed to log2 scale after the addition of one (log2(x+1). Clinical annotations were obtained in R environment using cbioportal's cgdsr package [88,89]. The dataset contained 1071 BRCA and 113 normal breast tissue samples. Among the normal tissues, a total of 109 samples had paired BRCA tissues while 4 remained unmatched and were, thereby, excluded. ERα, PR and HER2 status was determined by immunohistochemistry in the initial TCGA publication [87]. Normalized microarray data from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) [90] and corresponding clinical annotations were downloaded from cbioportal using the cgdsr package.

Gene Set Enrichment Analysis (GSEA)
GSEA was performed to uncover Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) biological processes that are positively correlated with PANX1 expression. KEGG and GO biological process gene sets were downloaded from the Molecular Signatures Database (MSigDB) (http://www.broadinstitute.org/gsea/). Analysis was done in the desktop GSEA software developed by the Broad Institute [91,92]. Genes were ranked based on their Pearson correlation with PANX1 expression level. Enrichment analysis was performed using 1000 phenotype permutations, and gene sets with nominal P value < 0.05 and false discovery rate (FDR) < 0.25 were considered significant.

Proteomic Analysis
Isobaric Tags for Relative and Absolute Quantification (iTRAQ) proteomics data of 105 TCGA BRCA samples were collected from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) data portal [93]. Analyzed samples covered 4 BRCA subtypes: Luminal A, Luminal B, Basal-like, and HER2-enriched. Protein expression was calculated as log2 Ratio to a pool of internal reference comprising a mixture of 40 samples with equal representation of the 4 breast subtypes. Among the 105 available samples, only 72 samples possessed PANX1 protein expression data.

Survival Analysis
Overall survival (OS) was defined as the duration from the time of diagnosis (in months) until death. Patients of TCGA or METABRIC data sets were segregated into High, Intermediate or Low PANX1 expression score groups based on the 25th and 75th percentiles of PANX1 expression. PANX1 High, Intermediate, and Low expression score groups were plotted using Kaplan-Meier survival curves. Kaplan-Meier plots were used to compare OS of High/Intermediate versus Low PANX1 expression groups.

Patients and Specimens
Breast carcinoma samples were retrieved from the pathology departments at the American University of Beirut Medical Center and Hammoud Hospital University Medical Center, after obtaining approval from the ethics committee (IRB reference number: PALM. FB. 01). A total of 37 patients with matched age (over 50 years old) and no treatment history were classified into 3 groups. Group I TNBC (11 specimens of triple negative breast cancer [TNBC]), group II HER2 + (11 specimens of estrogen receptor [ER] − /progesterone receptor [PR] -/human epidermal growth factor receptor 2 [HER2] + breast cancer), and group III HER2 -(15 specimens of ER + PR + HER2breast cancer). All groups were evaluated for histologic subtype, grade, central fibrosis, and tumor necrosis. A representative formalin-fixed paraffin-embedded (FFPE) tissue block from each case was obtained for molecular analyses. Negative controls were obtained from breast tissue of patients who underwent reduction mammoplasty [94]. The protocol of this study was approved by the American University of Beirut under the permission of the pathology departments at the American University of Beirut Medical Center and Hammoud Hospital University Medical Center (IRB PALM. FB. 01). All patients provided written informed consent to the surgical procedures and gave permission for the use of resected tissue specimens.

Gene Expression Analysis by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Gene expression levels of PANX1 were determined in archived breast cancer tissues. mRNA levels of PANX1 along with EMT markers were also determined in MCF-7, MDA-MB-231 and MDA-PANX1cells by qRT-PCR. Real time-PCR primers are listed in Table 1. Briefly, extraction of total RNA from cells was performed using RNeasy ® Plus Mini Kit (QIAGEN, Hilden, Germany) while extraction of RNA from mouse or human tissues was performed using TRIzol ® Reagent (Invitrogen), following the manufacturer's protocols and as previously described [95]. One µg of total RNA was reverse-transcribed to a single stranded complementary DNA (cDNA) using Revert Aid first strand cDNA synthesis kit (Thermo, Waltham, MA, USA). qRT-PCR was performed using iQ SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) in a CFX96™ Real-Time PCR Detection System (Bio-Rad). PCR amplification steps were as follows-95 • C for 10 s annealing temperature of the target gene for 30 s and then 72 • C for 30 s. The fluorescence threshold cycle value was obtained for each gene. ∆∆Cq method was used to calculate the relative fold change in gene expression after normalization to the housekeeping gene, GAPDH.

Protein Expression Analysis by Western Blotting
Proteins from cultured cells were extracted using a lysis buffer (126 mM Tris/HCl, 20% glycerol (v/v), 40 mg/mL of sodium dodecyl sulfate [SDS]) containing protease and phosphatase inhibitors (Roche, Penzberg, Germany). Protein concentrations were determined using DC Protein Assay II kit (Bio-Rad). Total protein lysates (100 µg) were resolved by SDS-PAGE and then transferred to PVDF membranes (Bio-Rad). The PVDF membranes were blocked with 5% skimmed milk for 1 hour before incubation with primary antibodies at a concentration of 1 µg/mL: anti-PANX1 (cat #710184, Invitrogen, Carlsbad, CA, USA) and anti-GAPDH (cat #G8795, Sigma). The membranes were then washed in PBS and incubated with appropriate horseradish peroxidase-conjugated IgG secondary antibody (Santa Cruz Biotechnology). Protein bands were visualized using chemiluminescence, and their intensity was quantified by densitometry and normalized to GAPDH using the Image J software (https://imagej.nih.gov/ij/).

Invasion and Proliferation Assays
Real-Time Cell Analysis (RTCA) was used to study invasion and proliferation of MDA-MB-231 and MDA-PANX1cells. Cells were seeded on a cellular invasion plate (CIM-plate 16) with micro-electronic sensors on the underside of an 8 µm microporous polyethylene terephthalate membrane [18,[96][97][98] of a Boyden-like upper chamber, and real-time evaluation of cell impedance was performed using xCELLigence RTCA [A2] DP instrument (Roche). For invasion assays, 30 µL of growth factor-reduced Matrigel TM (BD Biosciences, San Jose, CA, USA) diluted in serum-free medium at a ratio of 1:20 were used to coat the upper surface of the membrane, followed by incubation at 37 • C in 5% CO 2 for 4 h in a humidified incubator and then washed with PBS. A volume of 160 µL of RPMI-1640 containing 10% FBS was added to the lower chamber of each well and 30 µL to the upper chamber. The plate was incubated for 1 hour at 37 • C before cell seeding. Cells were seeded in the upper chamber at a density of 2 × 10 4 cells per 100 µL of serum-free media. MDA-MB-231 cells were allowed to adhere before PBN (1 mM) treatment was applied. For cell proliferation assay, cells were seeded in an E-plate as described above, at a density of 1 × 10 4 cells with an additional 150 µL of media containing 10% FBS. Both invasion and proliferation were monitored by recording cell impedance every 15 m for a minimum of 18 h.

Gelatin Zymography
Cell culture conditioned media were collected 48 h following incubation of MCF-7, PBN-treated MCF-7, MDA-MB-231, PBN-treated MDA-MB-231 and MDA-PANX1in serum-free media, and the enzymatic activity of metalloproteinase (MMP)-2 and MMP-9 was assessed by gelatin zymography. Concentrated protein supernatants (70 µg) were run on an SDS-PAGE gel, containing gelatin as a substrate. Gels were then stained with Coomassie TM brilliant blue R-250 (Bio-Rad) for one hour at room temperature and de-stained with a solution of methanol, acetic acid, and water. Band staining intensity was determined by densitometry, using Image J software.  [99]. Measurement of EtBr uptake by the cells was performed by recording EtBr fluorescence by live imaging for 15 m at 37 • C using a Laser Scanning Confocal Microscope (LSM710, Carl Zeiss, Oberkochen, Germany). The assay was also performed in normal divalent solution (145 mM NaCl, 5 mM KCl, 2 mM CaCl 2 , 1 mM MgCl 2 , 13 mM glucose, 10 mM HEPES, pH 7.3) and at room temperature, with the addition of 1 mM ATP to stimulate PANX1 channel opening, as previously documented [49]. Results from this assay are displayed in Figure S2.

Statistical Analysis
For survival analyses, the Log-rank test was used to estimate significance and hazard ratio (HR). Univariate Cox regressions were performed using dichotomic PANX1 score as variable. Prognosis significance was estimated using Wald's p value. Analyses were performed in R environment using survival and survminer packages (https://CRAN.R-project.org/package = survival; https://CRAN.Rproject.org/package = survminer).
Histograms, boxplots, and correlation plots were generated in R using ggpubr (https://CRAN. R-project.org/package = ggpubr). Pearson's method was used for gene correlation analysis in R as implemented in the corrplot package (https://github.com/taiyun/corrplot). Adjustment for multiple comparisons was done using the Benjamini and Hochberg method [100].
Microsoft Excel and GraphPad Prism software were also used to perform statistical analysis. Results are expressed as individual data or as mean ± standard deviation. Student's t-test was used to compare various groups. Statistical significance was determined by unpaired Student's t-test.

Conclusions
This study supports a role for PANX1 in cancer progression and metastasis, where its expression is upregulated in all breast cancer subtypes and upregulated PANX1 levels in breast cancer tissues are correlated with poor clinical outcomes. Furthermore, this study uncovers a role for PANX1 in EMT regulation in breast cancer. Overall, the ablating PANX1 function, pharmacologically or genetically, may constitute a novel multi-pronged therapeutic target for metastatic breast cancer by reverting EMT and decreasing the metastatic ability and invasiveness of breast cancer cells.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6694/11/12/1967/s1, Figure S1: PANX1 expression does not correlate with age. Figure S2: PBN attenuates ATP-induced EtBr dye uptake by PANX1 channels; Figure S3: Western blot images captured on the Chemidoc MP system; Table S1: Intensity ratios obtained by densitometry of PANX1 bands, normalized to densitometry of GAPDH bands, using Image Lab software. Funding: This work was funded and supported by grants from the Lebanese National Council for Scientific Research (LNCSR) to M.E.S., the American University of Beirut, and MPP and URB grants to M.E.S.