Transcriptional Profiling of Tumorspheres Reveals TRPM4 as a Novel Stemness Regulator in Breast Cancer

Cancer stem cells (CSCs) have been implicated in the development of chemoresistance, tumor recurrence and metastasis in breast cancer, thus emerging as a promising target for novel therapies. To identify novel stemness regulators that could potentially be targeted in luminal ER+ tumors, we performed RNA-sequencing (RNA-seq) in MCF-7 adherent monolayer cells and tumorspheres enriched in breast CSCs (bCSCs). We identified 1421 differentially expressed genes (DEGs), with 923 of them being upregulated and 498 downregulated in tumorspheres. Gene ontology and pathway enrichment analyses revealed that distinct gene networks underlie the biology of the two cell systems. We selected the transient receptor potential cation channel subfamily M member 4 (TRPM4) gene that had not been associated with cancer stemness before for further investigation. We confirmed that TRPM4 was overexpressed in tumorspheres and showed that its knock-down affected the stemness properties of bCSCs in vitro. TRPM4 inhibition revealed potential anti-tumor effects by directly targeting the bCSC subpopulation. We suggest that TRPM4 plays a key role in stemness mediation, and its inhibition may represent a novel therapeutic modality against bCSCs contributing in the improvement of breast cancer treatments.


Introduction
Breast cancer is the most prevalent cancer worldwide, with 2.3 million women having been diagnosed with the disease and 685,000 patients having died of it in 2020 [1]. Conventional chemotherapy is still at the forefront of treatment regimens, being effective in inciting tumor remission, yet a considerable fraction of patients ends up developing resistance and suffering relapse overtime [2]. This occurs even in the least aggressive subtype of BC, the luminal ER + , which is treated with chemotherapy in advanced cases, and has a good prognosis in the initial 5 years after diagnosis, but has an increased chronic annual risk of recurrence thereafter [3]. More than half of the metastases of ER + tumors occur 5 years or longer after diagnosis with some patients suffering recurrence after more than 20 years [4].
A possible explanation to account for this observation is based on the existence of rare tumor subpopulations with intrinsic resistance that survive treatment and drive tumor regrowth leading to patient relapse. Such tumor cells possess oftentimes stem-like characteristics, such as self-renewal and the capacity for differentiation, and are referred to as breast cancer stem cells (bCSCs). In breast, CSCs were first isolated from breast tumor specimens by Al-Hajj and colleagues [5]. The researchers determined that only a small fraction of cancer cells, those with the phenotype CD44 + CD24 -/low /lineage -, Biomedicines 2021, 9,1368 2 of 18 could initiate tumor formation, when transplanted to the mammary fat pad of severe combined immunodeficiency disease mice [5]. Later studies further confirmed that cells with this phenotype isolated from tumor samples or breast cancer cell lines possessed tumor-initiating properties and displayed increased resistance to conventional therapeutic schemes compared to the bulk of tumor cells [6,7]. As such, bCSCs have been implicated in the etiology of tumor recurrence and metastases [8] and their specific targeting has emerged as a promising approach complementary to the standard course of therapy, in order to improve clinical outcome in breast cancer [9].
To identify pharmacological targets in bCSCs, an in-depth characterization of the molecular mechanisms underpinning their biology is needed; however, the rarity of this tumor subpopulation has been a major hurdle for their study. Several methods of isolation and/or enrichment of bCSCs have been described in the literature [10]; the most widespread in vitro system for their study and manipulation is the generation of mammospheres, which are 3D-spheroid structures first developed by Dontu and colleagues for the propagation of human mammary epithelial cells under non-adherent, non-differentiated culture conditions [11]. These discrete cell clusters were enriched in progenitor cells capable of differentiating towards distinct mammary tissue lineages [11]. A later study showed that mammospheres cultured from breast cancer cell lines and primary tumors, also known as tumorspheres, were highly enriched in CD44 + CD24 -/low cells [12]. Since then, numerous reports have meticulously characterized mammospheres/tumorspheres generated from various cancer cell lines and patient tumor samples confirming that each one is derived directly from the proliferative clonal expansion of a single tumor initiating cell and is enriched in cells that retain stem-like or progenitor cell properties [6,13,14]; thus, the 3D tumorsphere system represents a useful and reliable tool to study bCSCs in vitro.
In this study, we used next generation sequencing (NGS) to characterize the transcriptome of MCF-7-derived mammospheres and parental cells in an effort to identify novel regulators of stemness that could be potential targets in luminal ER + tumors. Our bioinformatic analysis identified numerous differentially expressed genes including several known stemness markers validating our data, as well as novel genes that may play an important role in bCSCs. One of the latter ones was TRPM4 a member of the TRPM family of ion channels that are overexpressed in various types of cancer [15]. TRPM4 itself was recently shown to be upregulated in breast cancer [16], where it was associated with epithelial-to-mesenchymal transition (EMT) gene sets [16] spurring our interest to further investigate its role in bCSCs. Our results suggest that TRPM4 is a novel bCSC-regulator, and its targeting may lead to a significant reduction of this aggressive subpopulation and may enhance therapeutic results in breast cancer.

Cell Lines and Pharmacological Agents
The MCF-7 human breast cancer cell line was purchased from ATCC (LCG standards, Middlesex, UK) and was cultured in Dulbecco's modified Eagle's medium (DMEM, high glucose) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin in a humidified atmosphere of 5% CO 2 at 37 • C. Cells were routinely passaged every 2 or 3 days and tested for mycoplasma. The TRPM4 inhibitor used was 9-phenantrhol (Sigma-Aldrich St. Louis, MO, USA, cat. No. 211281). The chemotherapeutic drug used was doxorubicin (Adriblastina Hydrochloride 10 mg/5 mL VIAL Pfizer, New York, NY, USA).

TRPM4 Knock-Down
A small interfering RNA (siRNA) previously published [17] was used for transient TRPM4 knock-down and was transfected into cells using Lipofectamine™ RNAiMAX Transfection Reagent (ThermoFisher Scientific, Waltham, MA, USA) according to manufacturer's instructions. A scrambled siRNA was used as a negative control.

Pharmacological Treatment
Dose-response experiments with the TRPM4 inhibitor 9-phenanthrol were performed with MCF-7 tumorspheres. The MFE was calculated at different time points (3-10 days) and FACS analysis was also performed to monitor the CD44 + CD24 −/low bCSC-subpopulation. The most pronounced reduction in the MFE and in the number of bCSCs was achieved at 50 µM of 9-phenanthrol after treatment for 7 days. For the combination treatment of tumorspheres with 9-phenanthrol and doxorubicin, preliminary experiments were performed to determine the administration scheme for the two reagents that yielded the maximum response at the lowest concentrations (data not shown). The final protocol applied involved pretreatment of tumorspheres with 50 µM of 9-phenanthrol for 5 days and addition of 2.5 µM doxorubicin for 2 more days.

Immunoblot Analysis
For Western blot analysis, total protein quantitation was performed using the BCA Protein Assay (ThermoScientific, Waltham, MA, USA). The antibodies used were anti-

Immunofluorescence Analysis
MCF-7 breast cancer cells were seeded in coverslips and the next morning they were washed with PBS and fixed with 4% PFA. The cells were incubated with the primary antibody against TRPM4 (1:50) (cat. No. TA500381, Origene, Rockville, MD, USA) for 1 h at room temperature (RT) and with the secondary antibody (Goat anti-Mouse IgG Alexa Fluor 488, 1:400) (cat. No. A-1100, ThermoFisher Scientific, Waltham, MA, USA)) for 45 min at RT. For mammosphere staining, first generation MCF-7 mammospheres were collected by centrifugation (800 rpm, 3 min, RT), washed in PBS (800 rpm, 3 min, RT) and fixed in 4% PFA for 30 min. The spheres were incubated with the primary antibody overnight at 4 • C and with the secondary antibody at RT for 45 min. In both cases, TOPRO-3 (cat. No. T3605, Invitrogen, ThermoFisher Scientific, Waltham, MA, USA) was used for DNA staining. A Leica SP5 confocal microscope was used to image the specimens. For comparison reasons, the same parameters (PMTs and offset) were employed for imaging, and the images obtained were processed in the same manner.

RNA Extraction, cDNA Synthesis and q-RT-PCR
For total RNA extraction the RNeasy Kit (Qiagen) was used. For total RNA extraction from mammospheres, first generation mammospheres were dissociated and replated in mammosphere forming conditions. Four days upon replating, spheres were collected by centrifugation (800 rpm, 3 min, RT) and processed for RNA extraction. 2.9. RNA-Sequencing (RNA-seq) and Bioinformatic Analysis RNA-seq libraries were prepared with the TruSeq RNA v2 kit (Illumina, San Diego, CA, USA) using 1 µg of total RNA. The libraries were checked with the Agilent bioanalyzer (DNA1000 chip) (Agilent, Santa Clara, CA, USA), quantitated with the qubit HS spectrophotometric method and pooled in equimolar amounts for sequencing. Approximately 25 million, 75 bp long, single-end reads were generated for each sample on an Illumina NextSeq500 or Novaseq 6000 sequencer. For each sample, two biologically independent replicates were sequenced yielding highly similar results.
Quality Control was performed with the fastq raw data file using the "FASTQC" software (GPL v.3, Babraham Institute, Cambridge, UK) FastQ files were aligned to the human genome (hg19) using HISAT2 (HISAN 3N-beta release) [19]. Counts were defined using the HTSeq htseq-count command with the "intersection non empty" option [20]. Normalization was performed with the estimate size factor function followed by Differentially Expressed Genes Analysis. The count files were used as input for DESeq2 (Bioconductor version: Release (3.13) [21] for identification of DEGs between mammospheres and adherent monolayer cells with a statistically significant cut-off value of p < 0.05. Additional cut-off criteria were set as fold-change ≥ 2, p-adjust < 0.01 and number of reads > 10.
All raw and processed data files have been submitted to GEO (accession number GSE182532).
The volcano plot was designed using GraphPad Prism 8.

Gene Ontology and Gene Set Enrichment Analysis of DEGs
Gene ontology (GO) analysis for the DEGs was performed using the PANTHER classification system (version 16) [22]. The annotation data set selected was PANTHER GO-Slim Biological Process. Only categories with a p-value < 0.05 (as determined by Fisher's exact test) and an FDR < 0.05 were further analyzed. Biological processes that were underrepresented in our data, as well as the ones that had to do with neurogenesis, were excluded from further analysis.
Gene set enrichment analysis for the DEGs was also performed using the GSEA software [23,24]. The number of permutations were set to 1000, the permutation type was 'gene_set' and gene sets with <15 and >500 genes were excluded [25]. The gene set database selected was "hallmarks" that provides gene sets that are well-defined biological states or processes.

Statistical Analysis
Data presented were generated from at least three independent biological experiments, unless otherwise stated, and they are expressed as mean ± SE. An unpaired two-tailed Student's t-test was used for statistical analysis.

RNA Sequencing Data Analysis
Mammospheres or tumorspheres are 3D spheroids generated by breast cancer cells grown under non-adherent conditions and are significantly enriched in bCSCs [13,26]. We confirmed that this was also the case in our MCF-7 mammosphere system before proceeding with subsequent experiments. Indeed, these spheres exhibited an~11-fold enrichment in the CD44 + CD24 −/low bCSC-subpopulation (Supplementary Figure S1A,B) and overexpression of known stemness genes (Supplementary Figure S1C), when compared to MCF-7 cells grown in adherent monolayers. We used this system to identify novel stemness regulators in breast cancer. For this purpose, we analyzed the transcriptome profiles of MCF-7 breast cancer cells grown as 2D monolayers or 3D spheres by RNA-seq. In total, we analyzed 27,967 genes that mostly exhibited a similar expression pattern between the two cell systems (Pearson correlation coefficient r = 0.913), as expected ( Figure 1A). Mammospheres expressed 14,937 genes, with 496 of them being expressed only in this system, while adherent monolayer cells expressed 14,641 genes with 316 uniquely expressed genes (Supplementary Figure S2).

Statistical Analysis
Data presented were generated from at least three independent biological experiments, unless otherwise stated, and they are expressed as mean ± SE. An unpaired twotailed Student's t-test was used for statistical analysis.

RNA Sequencing Data Analysis
Mammospheres or tumorspheres are 3D spheroids generated by breast cancer cells grown under non-adherent conditions and are significantly enriched in bCSCs [13,26]. We confirmed that this was also the case in our MCF-7 mammosphere system before proceeding with subsequent experiments. Indeed, these spheres exhibited an ~11-fold enrichment in the CD44 + CD24 -/low bCSC-subpopulation (Supplementary Figure S1A,B) and overexpression of known stemness genes (Supplementary Figure S1C), when compared to MCF-7 cells grown in adherent monolayers. We used this system to identify novel stemness regulators in breast cancer. For this purpose, we analyzed the transcriptome profiles of MCF-7 breast cancer cells grown as 2D monolayers or 3D spheres by RNA-seq. In total, we analyzed 27,967 genes that mostly exhibited a similar expression pattern between the two cell systems (Pearson correlation coefficient r = 0.913), as expected ( Figure 1A). Mammospheres expressed 14,937 genes, with 496 of them being expressed only in this system, while adherent monolayer cells expressed 14,641 genes with 316 uniquely expressed genes (Supplementary Figure S2).
To assess the validity of our RNA-seq data, we identified the top 100 most highly expressed genes in the MCF-7 monolayer adherent cells and in the mammospheres (Supplementary Table S1) and cross-referenced some of them to published literature. In both cases the keratins KRT8 and KRT18, which are considered typical luminal markers [27,28], were very highly expressed. We also detected elevated levels of CSDE1, which encodes for an RNA-binding protein that plays an important role in several cellular processes, such as cell cycle and apoptosis, and has been shown before to be overexpressed in MCF-7 cells [29]. High expression levels of the elongation factor 2 (EEF2) have been associated with poor prognosis in hormone receptor-positive breast cancer [30]. Other well-known breast cancer associated genes that were highly expressed in both cell systems included GATA3 [31], GNAS [32], FASN [33] and NCOA3 [34]. To assess the validity of our RNA-seq data, we identified the top 100 most highly expressed genes in the MCF-7 monolayer adherent cells and in the mammospheres (Supplementary Table S1) and cross-referenced some of them to published literature. In both cases the keratins KRT8 and KRT18, which are considered typical luminal markers [27,28], were very highly expressed. We also detected elevated levels of CSDE1, which encodes for an RNA-binding protein that plays an important role in several cellular processes, such as cell cycle and apoptosis, and has been shown before to be overexpressed in MCF-7 cells [29]. High expression levels of the elongation factor 2 (EEF2) have been associated with poor prognosis in hormone receptor-positive breast cancer [30]. Other well-known breast cancer associated genes that were highly expressed in both cell systems included GATA3 [31], GNAS [32], FASN [33] and NCOA3 [34].
Technical validation of the RNA-seq data was performed by RT-qPCR experiments as described in the next section.

Identification and Validation of Differentially Expressed Genes
Next, we focused on identifying the DEGs between mammospheres and parental adherent cells, aiming to discover novel potential regulators of bCSCs. Their transcrip-Biomedicines 2021, 9, 1368 6 of 18 tomes were compared using the R package DESEQ2 and DEGs were identified by setting a 2-fold-change and a p-value ≤ 0.05 (with p-adj ≤ 0.01) as cut-off values. In this manner, we identified 1421 DEGs with 923 of them being upregulated and 498 downregulated in mammospheres ( Figure 1B). A complete list of these genes is provided in the Supplementary Table S2. Notably, among the most upregulated genes we find the well-known stemness marker CXCR4 [35,36], CYP1A1, which was reported to control bCSCs proliferation, development and self-renewal [37], and TM4SF1 that was shown to promote breast cancer stem cell traits [38]. In the downregulated genes, we find CAV1, which was recently reported to be downregulated in bCSCs and was associated with a metabolic switch from mitochondrial respiration to aerobic glycolysis [39] and ELOVL2, an enzyme involved in lipid metabolism with a reduced expression in a spheroid-induced EMT model [40].
We also validated the results of our DEG analysis by assessing the expression levels of several DEGs by RT-qPCR. A side-by-side comparison of these results with the corresponding RNA-seq data is shown in Figure 2. We selected to examine genes that showed high (55-230), medium (22)(23)(24)(25)(26)(27)(28) or low (6)(7)(8)(9)(10)(11)(12) fold-change of expression levels in the mammospheres compared to monolayer cells in the RNA-seq data. Overall, the changes in mRNA levels detected by the two techniques were similar for most genes examined with the exception of KCNE4 in the upregulated ( Figure 2A) and of the TGM2 in the downregulated genes ( Figure 2B). In both cases, the RNA-seq experiments overestimated the fold-change of gene expression in the mammospheres compared to adherent cells; however, the pattern of expression (up or down) was verified.
Technical validation of the RNA-seq data was performed by RT-qPCR experiments as described in the next section.

Identification and Validation of Differentially Expressed Genes
Next, we focused on identifying the DEGs between mammospheres and parental adherent cells, aiming to discover novel potential regulators of bCSCs. Their transcriptomes were compared using the R package DESEQ2 and DEGs were identified by setting a 2fold-change and a p-value ≤ 0.05 (with p-adj ≤ 0.01) as cut-off values. In this manner, we identified 1421 DEGs with 923 of them being upregulated and 498 downregulated in mammospheres ( Figure 1B). A complete list of these genes is provided in the Supplementary Table S2. Notably, among the most upregulated genes we find the well-known stemness marker CXCR4 [35,36], CYP1A1, which was reported to control bCSCs proliferation, development and self-renewal [37], and TM4SF1 that was shown to promote breast cancer stem cell traits [38]. In the downregulated genes, we find CAV1, which was recently reported to be downregulated in bCSCs and was associated with a metabolic switch from mitochondrial respiration to aerobic glycolysis [39] and ELOVL2, an enzyme involved in lipid metabolism with a reduced expression in a spheroid-induced EMT model [40].
We also validated the results of our DEG analysis by assessing the expression levels of several DEGs by RT-qPCR. A side-by-side comparison of these results with the corresponding RNA-seq data is shown in Figure 2. We selected to examine genes that showed high (55-230), medium (22)(23)(24)(25)(26)(27)(28) or low (6)(7)(8)(9)(10)(11)(12) fold-change of expression levels in the mammospheres compared to monolayer cells in the RNA-seq data. Overall, the changes in mRNA levels detected by the two techniques were similar for most genes examined with the exception of KCNE4 in the upregulated ( Figure 2A) and of the TGM2 in the downregulated genes ( Figure 2B). In both cases, the RNA-seq experiments overestimated the foldchange of gene expression in the mammospheres compared to adherent cells; however, the pattern of expression (up or down) was verified.

Gene Ontology and Gene Set Enrichment Analysis for DEGs
To further investigate the biological significance of our data, we performed GO analysis for the DEGs, using the PANTHER classification system [22]. PANTHER GO-Slim biological process analysis was performed separately for the up-and down-regulated genes in mammospheres, using the statistical overexpression test.

Gene Ontology and Gene Set Enrichment Analysis for DEGs
To further investigate the biological significance of our data, we performed GO analysis for the DEGs, using the PANTHER classification system [22]. PANTHER GO-Slim biological process analysis was performed separately for the up-and down-regulated genes in mammospheres, using the statistical overexpression test.
We determined a total of 50 overrepresented biological processes in the upregulated genes and 33 in the downregulated genes in mammospheres, clustered in 8 and 5 major groups, respectively (Supplementary Table S3). Each group contained similar (or identical) biological processes, as recognized by the corresponding GO terms and GO terms grouped together included mostly common genes (data not shown). For simplicity and clarity reasons, we chose to illustrate one GO term per group, the one that contained the largest number of genes (Supplementary Figure S3). In the upregulated DEGs in MCF-7 mammospheres, the most numerous groups were "signaling" and "developmental process" (Supplementary Figure S3A), while in the downregulated DEGs we had "cellular component organization or biogenesis" and "developmental process" (Supplementary Figure S3B).
For a more informative and biologically relevant depiction of our data, we also ranked the biological processes by fold-enrichment with a cut-off value of 2. Figure 3 illustrates the GO term/biological pathway that had the highest fold-enrichment score in each group. The most overrepresented biological pathway in the upregulated genes was the "unsaturated fatty acid metabolic process" (Figure 3A). The altered metabolism is one of the hallmarks of cancer. It has been reported that CSCs have increased needs in fatty acids and that unsaturated fatty acids maintain cancer cell stemness [41]. The Wnt pathway, an established CSC-promoting signaling cascade [42], was also highly enriched in the upregulated DEGs ( Figure 3A). "Regulation of cell migration" was another overexpressed pathway in our analysis, a property that in CSCs has been correlated with metastasis [43]. Migration of cells is affected by the extracellular matrix (ECM) distribution and it has been reported that ECM of CSCs differs in composition compared to normal cells, playing a role in cancer stemness [44]. Divalent metal ions have a major role in tumor progression through modulation of the structure and function of various components of the ECM [45]. Consequently, the associated GO terms/biological pathways that we found enriched in the upregulated DEGs ( Figure 3A) are in accordance with previous studies, further validating our results.
In the downregulated DEGs, the most overrepresented biological pathway was the "transforming growth factor beta receptor signaling pathway" ( Figure 3B). Transforming have shown that TGFβ suppresses breast cancer tumorigenesis by reducing the CSC pools or by promoting the differentiation, whereas other studies have reported to promote or sustain stemness of the pool of CSCs in breast cancer [46]. In addition, "stem cell differentiation", "cell growth", "regulation of cell cycle process" and "negative regulation of locomotion" pathways were overrepresented in the downregulated DEGs, as one would expect to find in a cell population enriched for slowly proliferating, metastasis-prone CSCs ( Figure 3B).
Next, we performed hallmark analysis, using the GSEA software [23,24], implementing separate analysis for the up-and downregulated genes in mammospheres. The analysis revealed 11 enriched gene sets in the upregulated DEGs and 2 enriched gene sets in the downregulated DEGs in mammospheres ( Figure 3C and Supplementary Table S4). It is noteworthy that EMT and hypoxia, two of the most widely recognized CSC hallmarks, were among the most represented ones in the upregulated DEGs. Taken together the results described above support the validity of our RNA-seq data and particularly the use of mammospheres as a discovery tool for novel breast cancer stemness regulators.

TRPM4 Is Overexpressed in MCF-7 Mammospheres
The GO analysis of the upregulated genes in the MCF-7 mammospheres confirmed overrepresentation of several well-known stemness associated pathways, but also identified one that, to our knowledge, has not been extensively studied in bCSCs before. This was "inorganic ion homeostasis" (Figure 3A), which was clustered with similar biological processes under the more general GO term "ion transport" (Supplementary Table S3). Manual inspection of the genes in this group revealed several genes that had not been associated with cancer stemness before, one of them being the TRPM4 gene, a member of the TRPM family of ion channels [15]. TRPM4 is known to be overexpressed in breast [16] and other types of cancer [47]. We confirmed these data by interrogating the Oncomine data base [48]. We found that TRPM4 was overexpressed in aggressive breast tumors in two separate patient data sets (Supplementary Figure S4). These findings prompted us to investigate further its role in bCSCs. associated with cancer stemness before, one of them being the TRPM4 gene, a member of the TRPM family of ion channels [15]. TRPM4 is known to be overexpressed in breast [16] and other types of cancer [47]. We confirmed these data by interrogating the Oncomine data base [48]. We found that TRPM4 was overexpressed in aggressive breast tumors in two separate patient data sets (Supplementary Figure S4). These findings prompted us to investigate further its role in bCSCs. Firstly, we sought to confirm our RNA-seq data, which showed that TRPM4 mRNA expression levels were 2.4-fold higher in MCF-7 mammospheres compared to the adherent cells grown in monolayers ( Figure 4A, blue bars). Indeed, RT-qPCR experiments validated these findings, showing a ~3-fold overexpression of TRPM4 in mammospheres (Figure 4A, yellow bars). Western blot experiments showed that TRPM4 protein levels were ~2.7-fold higher in the MCF-7-derived mammospheres compared to the corresponding Firstly, we sought to confirm our RNA-seq data, which showed that TRPM4 mRNA expression levels were 2.4-fold higher in MCF-7 mammospheres compared to the adherent cells grown in monolayers ( Figure 4A, blue bars). Indeed, RT-qPCR experiments validated these findings, showing a~3-fold overexpression of TRPM4 in mammospheres ( Figure 4A, yellow bars). Western blot experiments showed that TRPM4 protein levels were~2.7-fold higher in the MCF-7-derived mammospheres compared to the corresponding attached cells ( Figure 4B). Immunostaining of cells grown in monolayers or mammospheres with an antibody against TRPM4 further corroborated its overexpression in the latter ones and showed that it is mostly localized on the cell membrane as expected ( Figure 4C).
Overall, these data demonstrate convincingly that TRPM4 is overexpressed in MCF-7 mammospheres both on the mRNA and protein level. attached cells ( Figure 4B). Immunostaining of cells grown in monolayers or mammospheres with an antibody against TRPM4 further corroborated its overexpression in the latter ones and showed that it is mostly localized on the cell membrane as expected ( Figure 4C).
Overall, these data demonstrate convincingly that TRPM4 is overexpressed in MCF-7 mammospheres both on the mRNA and protein level.

TRPM4 Regulates the Stemness Properties of bCSCs In Vitro
The upregulation of TRPM4 expression in the MCF-7 mammospheres, known to be highly enriched in bCSCs, urged us to investigate its potential role in the regulation of this subpopulation.
To this end, we knocked down TRPM4 expression in MCF-7 cells (Supplementary Figure S5) using a previously published siRNA [17] and examined its effects on the stemness properties of bCSCs by applying the mammosphere-forming assay ( Figure 5A,B). The results of these experiments revealed that TRPM4 knock-down resulted in the formation of a significantly smaller number of mammospheres compared to control ( Figure 5A). Specifically, the TRPM4 siRNA-treated cells had 32% lower efficiency in forming mammospheres compared to the scramble siRNA-treated ones ( Figure 5B). To further show that the decreased mammosphere forming efficiency was a specific effect of the TRPM4 knockdown on the bCSCs, we performed FACS analysis to monitor the CD44 + CD24 −/low CSC subpopulation. The last day of the experiment the mammospheres in each condition were collected, dissociated and stained against the surface markers CD44 and CD24 ( Figure 5C). FACS analysis showed 40% reduction in the percentage of the CD44 + CD24 −/low CSC subpopulation upon TRPM4 knock-down compared to control cells ( Figure 5D).
The upregulation of TRPM4 expression in the MCF-7 mammospheres, known to be highly enriched in bCSCs, urged us to investigate its potential role in the regulation of this subpopulation.
To this end, we knocked down TRPM4 expression in MCF-7 cells (Supplementary Figure S5) using a previously published siRNA [17] and examined its effects on the stemness properties of bCSCs by applying the mammosphere-forming assay (Figures 5A,B). The results of these experiments revealed that TRPM4 knock-down resulted in the formation of a significantly smaller number of mammospheres compared to control ( Figure  5A). Specifically, the TRPM4 siRNA-treated cells had 32% lower efficiency in forming mammospheres compared to the scramble siRNA-treated ones ( Figure 5B). To further show that the decreased mammosphere forming efficiency was a specific effect of the TRPM4 knock-down on the bCSCs, we performed FACS analysis to monitor the CD44 + CD24 -/low CSC subpopulation. The last day of the experiment the mammospheres in each condition were collected, dissociated and stained against the surface markers CD44 and CD24 ( Figure 5C). FACS analysis showed 40% reduction in the percentage of the CD44 + CD24 -/low CSC subpopulation upon TRPM4 knock-down compared to control cells ( Figure 5D). The above data strongly support the role of TRPM4 as a key regulator of the stemness properties of bCSCs in vitro.
Stemness-associated signaling cascades, such as the Wnt, Hedgehog and Notch pathways, regulate the maintenance of bCSCs (reviewed in [49]). Some of them are also involved in the induction of EMT that has been associated with the generation of bCSCs (reviewed in [50]). To gain some insight into the mechanisms that TRPM4 employs to The above data strongly support the role of TRPM4 as a key regulator of the stemness properties of bCSCs in vitro.
Stemness-associated signaling cascades, such as the Wnt, Hedgehog and Notch pathways, regulate the maintenance of bCSCs (reviewed in [49]). Some of them are also involved in the induction of EMT that has been associated with the generation of bCSCs (reviewed in [50]). To gain some insight into the mechanisms that TRPM4 employs to regulate breast cancer stemness, we selected several molecules that participate in these pathways and showed considerable expression in our system and interrogated their mRNA levels in TRPM4 k/d MCF-7 tumorspheres by RT-qPCR. The results are presented in . Among the EMT-associated molecules examined (blue in Figure S6), we found that the knock-down of TRPM4 led to significantly reduced mRNA levels of vimentin, an EMT marker. However, we did not see any changes in the levels of CDH1, which encodes for e-cadherin, or in the transcription factors SNAI1 and TWIST. Among the Wnt-associated molecules examined (green in Figure S6), PRICKLE1 was significantly downregulated in the TRPM4 k/d MCF-7 tumorspheres. This gene is a member of the non-canonical Wnt/planar cell polarity (PCP) pathway and was recently described to be involved in breast cancer cell dissemination [51]. We did not observe any changes in the expression levels of HHIP (yellow in Figure S6) or NOTCH1 (orange in Figure S6) involved in the Hedgehog and Notch signaling pathways, respectively.

TRPM4 Is a Potential Druggable Target in bCSCs
As mentioned above, several studies have reported on the overexpression of TRPM4 in different types of cancer [52][53][54][55], and even though, its precise mechanism of action in tumors is currently unknown, it is already portrayed as a promising anticancer drug target [47].
Bearing this in mind and based on our data described above, we hypothesized that pharmacological inhibition of TRPM4 could, potentially, have anti-tumor effects through targeting the bCSC-subpopulation.
To this end, we examined whether 9-phenanthrol, a TRPM4 channel blocker, could affect bCSCs. MCF-7 mammospheres were dissociated to single cells, re-plated under mammosphere-forming conditions and treated with 50 µM of the inhibitor for 7 days; at the end of treatment, the MFE was calculated ( Figure 6A,B). Inhibition of TRPM4 resulted in 45% decrease in the number of mammospheres compared to vehicle-treated one ( Figure 6A,B). At the same time, FACS analysis of the inhibitor-treated mammospheres showed an approximately 50% decrease in the CD44 + CD24 −/low CSC-subpopulation compared to control mammospheres ( Figure 6C,D).
These results strongly suggest that TRPM4 inhibition diminishes the stemness potential of bCSCs and targets directly this subpopulation in MCF-7 mammospheres. Notably, TRPM4 pharmacological inhibition phenocopied the effects of the siRNA-mediated knockdown, as described previously ( Figure 5). Therefore, blocking of this ion channel may serve as a novel therapeutic modality for targeting the associated stemness pathways in bCSCs.
As it has been suggested before [56] and we have also shown [18], both a CSC-targeted agent and a conventional cytotoxic drug are needed to achieve significant tumor reduction and, presumably, improve cancer treatment. Therefore, we sought to investigate whether TRPM4 inhibition in combination with doxorubicin treatment could have such an effect. MCF-7-derived tumorspheres were treated for 5 days with 9-phenanthrol (50 µM) and for an additional two days with doxorubicin (2.5 µM). On the last day of treatment, the number of tumorspheres was counted ( Figure 6E,F). Monotreatment with the chemotherapeutic agent led only to a 20% reduction in the number of spheres; on the other hand, pretreatment with the TRPM4 inhibitor followed by doxorubicin administration led to a 70% decrease in the number of spheres compared to the vehicle-treated ones ( Figure 6E,F). These results show that combination of a chemotherapeutic drug and a TRPM4 inhibitor is more successful in eliminating breast tumorspheres, probably because both stem/progenitor and differentiated cancer cells are targeted effectively. This was further confirmed by FACS experiments that were carried out in tumorspheres during the last day of treatment in order to monitor the CD44 + CD24 −/low CSC-subpopulation ( Figure 6G,H). Doxorubicin treatment alone could not target the CSC-pool and it rather enriched the culture in bCSCs. However, its combination with the TRPM4 inhibitor was efficient in reducing the size of the CD44 + CD24 −/low CSC subpopulation by 35% ( Figure 6G,H).
In conclusion, our data suggest that TRPM4 is a mediator of stemness in breast cancer in vitro and its inhibition may hold promise as an efficient complementary therapeutic approach along with standard treatments.
Biomedicines 2021, 9, x FOR PEER REVIEW 12 of 18 In conclusion, our data suggest that TRPM4 is a mediator of stemness in breast cancer in vitro and its inhibition may hold promise as an efficient complementary therapeutic approach along with standard treatments.

Discussion
Current anti-cancer therapeutic schemes have been mostly designed based on the prevalent tumor lesions identified through the employment of genome-wide sequencing techniques that analyze the bulk of the tumor. This strategy, while useful, ignores the molecular aberrations of small but, potentially, aggressive tumor subpopulations that may have different drug sensitivities [57]. Such subpopulations are often endowed with stemness properties and are referred to as cancer stem cells. More specifically in breast cancer, CSCs have been described as manifesting remarkable resistance to standard treatments [58,59], making them prime candidates for the development of novel, targeted therapies that could be administered in combination with commonly used drugs, aiming to achieve complete tumor elimination [60,61].
To develop such therapies, the first crucial step is to isolate bCSCs for further characterization and manipulation in vitro. Among the methods that have been developed to enrich for bCSCs, the mammosphere or tumorsphere culture system has been established as a particularly advantageous model for their investigation [10]. In this study, we used the widely applied and extensively characterized 3D mammosphere culture system developed from the luminal ER + breast cancer MCF-7 cells. Notably, this is one of the few cell lines where the CD44 + C24 −/low bCSC-subpopulation has been isolated, transplanted into mice and shown to possess stemness properties [62]. We used this system to employ a nonbiased genome-wide approach for the identification of the transcriptional networks that are distinct between cells grown in 2D monolayers and 3D spheres and are potentially essential for the maintenance of the CSC-pool in tumors. The quality of our sequencing data was extensively validated by RT-qPCR experiments; both methodologies yielded similar results regarding the fold-change in mRNA levels between monolayers and mammospheres for most genes examined. We also performed a literature search for some of the genes that were found to be overexpressed in both cell systems and we confirmed that they were breast cancer-associated ones according to several published reports [27][28][29][30][31][32][33][34]. We also performed DEG analysis. Among the most highly upregulated genes in mammospheres we detected several known stemness-associated genes [35][36][37][38]. Other highly upregulated genes in mammospheres included the ion channel KCNE4 and the monoamine oxidase MOAB, which have not been associated with cancer stemness before. Interestingly, both genes have been found to be overexpressed in glioblastoma [63,64] a tumor type highly enriched in CSCs [65].
Our gene ontology analysis offered further validation of our data, but also provided us with new insights into the biology of bCSCs and enriched our record of potential specific targets. A stark example was the "unsaturated fatty acid metabolism", which emerged as the most overrepresented GO term in the upregulated genes in mammospheres. This is in accordance with recent literature that shows that CSCs require more unsaturated fatty acids than their non-stem counterparts (reviewed in [66]) and thus, they overexpress the related metabolic pathways. However, not all key genes involved have been identified and studied yet. Several such candidate genes are listed in our supplementary data.
Probably the least studied pathways in CSCs are the ones regarding ion transport and ion homeostasis, even though recent studies have demonstrated that ion channels are present in various CSCs (reviewed in [67]). Thus, we turned our interest to these gene lists in an effort to discover novel CSC-regulators. We hypothesized that by identifying and modulating such genes, we might be able to suppress the bCSC-population and highlight a new therapeutic avenue in BC. We identified TRPM4 as an overexpressed ion channel that had not been associated with cancer stem cells before. This gene is one of the 8 members of the TRPM (transient receptor potential melastatin) family, which belongs to the superfamily of TRP cation channels [15]. The TRPM family members are widely expressed transmembrane proteins that contribute to intracellular calcium homeostasis by promoting Ca 2+ entry into the cytosol in response to various stimuli with the exception of TRPM4 and 5 that are Ca 2 + impermeable [15]. The TRPM4 channel is only permeable to monovalent cations, mostly to Na + and K + , upon activation by an increase in intracellular Ca 2 + [15]. Most interestingly to us, TRPM4 protein levels are increased in a variety of tumor types compared to healthy tissue, where it contributes to important cancer features, such as increased proliferation and migration and cell cycle shift (reviewed in [15,47]). Particularly in breast cancer, two recent studies showed that both the mRNA and protein levels of TRPM4 were upregulated in specimens from breast cancer patients [16,55]. Increased TRPM4 staining intensity in the membrane and cytosol of tumor cells was significantly correlated with a worse patient prognosis [55]. These data were also confirmed by our own Oncomine analysis that showed that TRPM4 overexpression is associated with aggressive breast tumors in the public datasets examined.
Taking into account these reports as well as the pathway analysis of our RNA-seq data, we hypothesized that TRPM4 might play an important role in the regulation of breast CSCs. After confirming the overexpression of TRPM4 in MCF-7 mammospheres, we set out to examine its functional effects on the same model of CSC-enrichment by knockingdown its gene expression. TRPM4 silencing reduced the CSC fraction and impaired its stemness properties, suggesting that this ion channel is critical for CSC maintenance and activity in breast cancer. These results were also reproduced, when an TRPM4 inhibitor, 9-phenanthrol, was used. Our experiments with the combined use of 9-phenanthrol and doxorubicin in an in vitro tumorsphere system showed that this administration scheme was more efficient in reducing the number of spheres than the drug alone, which could not target the bCSC-subpopulation. This is in agreement with our previous results [18] and the view expressed by others [47] that both a CSC-specific agent and a chemotherapeutic drug are needed for complete tumor elimination.
It should be mentioned that 9-phenanthrol also inhibits TMEM16A/ANO1 [68], a calcium-activated chloride channel that is expressed in very low levels in our cells according to our RNA-seq data (data not shown) and probably not on the protein level. Hence, it is safe to assume that the effects observed in our cell system after 9-phenanthrol treatment are due solely to inhibition of the TRPM4 channel.
The precise mechanisms of action of TRPM4 in cancer are currently under investigation [47]. Our results suggest that it may function through maintaining the CSCcompartment in tumors. The pathways that mediate TRPM4 activity in cancer stemness are yet unknown; our preliminary data suggest that this ion-channel may be involved in the EMT process. In support of this notion, recent studies have highlighted the connection of TRPM4 expression with activation of the EMT in prostate [54] and breast [16] cancer cells, a process that has long been associated to cancer stemness [69]. TRPM4 may also be involved in the regulation of the Wnt/PCP signaling pathway, which is often aberrantly expressed in cancer [70]. However, these findings need further investigation. As we could examine only a small set of stemness-associated genes by RT-qPCR, it is plausible that other signaling cascades also mediate TRPM4 actions.
Our results also suggest that TRPM4 can be used as a druggable target in breast cancer. Indeed, the cell membrane localization of the protein renders it an ideal candidate for novel anti-cancer drugs that could either fall into the small-molecule inhibitors category [71] or they could be specific antibodies that block the protein [72].
Future studies including in vivo experiments will clarify the role of TRPM4 in cancer, will elucidate the regulatory mechanisms it employs to sustain CSCs and will evaluate its potential as a therapeutic target in breast tumors.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/biomedicines9101368/s1: Figure S1. MCF-7 mammospheres are highly enriched in bCSCs, Figure S2. Expressed genes in MCF-7 cells grown as adherent monolayers or mammospheres, Figure S3. Gene ontology analysis for DEGs, Figure S4. Oncomine analyses of breast cancer databases, Figure S5. Western blot analysis for TRPM4 siRNA knock-down in MCF-7 cells, Figure S6. mRNA expression levels of genes associated with stemness pathways in TRPM4 k/d MCF-7 mammospheres. Table S1. The top 200 genes expressed in MCF-7 adherent monolayer cells and mammospheres. Table S2. The 1421 differentially expressed genes (DEGs) between adherent monolayer cells and mammospheres identified by RNA-sequencing, Table S3. The genes of each biological pathway overrepresented in DEGs, Table S4. The genes of the enriched hallmarks in DEGs. Funding: We acknowledge support of this work by the project "Advanced Research Activities in Biomedical and Agro alimentary Technologies" (MIS 5002469), which is implemented under the "Action for the Strategic Development on the Research and Technological Sector", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All raw and processed data files have been submitted to GEO (accession number GSE182532).