Bladder | 2023 | Vol. 10 | e21200011
DOI: 10.14440/bladder.2023.852
ARTICLE

Single-cell analyses EMP1 as a marker of the ratio of M1/M2 macrophages is associated with EMT, immune infiltration, and prognosis in bladder cancer

Jinqiao Li1†, Jianyu Liu2†, Honglei Wang1†, Jinpeng Ma1, Yueze Wang1, Wanhai Xu1*
1Urology Surgery Department, Harbin Medical University Cancer Hospital, Harbin, Heilongjiang 150081, PR China.
2Breast Surgery Department, Harbin Medical University Cancer Hospital, Harbin, Heilongjiang 150081, PR China.
†These authors contributed equally to this work.
*Corresponding author: Wanhai Xu, e-mail: xuwanhai@hrbmu.edu.cn
Conflict of interest: The authors have declared that no competing interest exists.
Abbreviation used: EMP1: epithelial membrane protein 1; EMT: epithelial–mesenchymal transition; TAM: Tumor-associated macrophages; TCGA: The Cancer Genome Atlas; GSEA: Gene Set Enrichment Analysis; GSVA: Gene Set Variation Analysis; GO: Gene Ontology; ECM: Extracellular Matrix; FA: Focal Adhesion; FN1: fibronectin 1; SDC1: Syndecan 1 GAS3: growth arrest-specific 3; PMP22: peripheral myelin protein 22; GEO: Gene Expression Omnibus
Received October 12, 2022; Revision received January 5, 2023; Accepted February 17, 2023; Published December 18, 2023
Abstract
Background: Bladder cancer is among the most lethal urinary system cancers across the globe. Macrophage 1 and Macrophage 2 play an essential role in the pathogenesis of tumors. Nevertheless, prior studies failed to investigate the implication of the two cells, working in combination, in the development, growth, progression and metastasis of bladder cancer.
Methods: We computed the M1/M2 ratio of the samples retrieved from The Cancer Genome Atlas (TCGA) by using the Cibersortx algorithm and calculated the ratio in 32 patients in our series by employing flow cytometry. SurvivalRandomForest was utilized to reduce the dimension of the list of the M1/M2-related genes, with an aim to obtain the most survival-predictive gene (EMP1) encoding epithelial membrane protein 1 (EMP1). The EMP1 was biologically characterized by using Gene Set Enrichment Analysis (GSEA), Gene Set Variation Analysis (GSVA), and Gene Ontology (GO). The single-cell transcriptome (sc-RNA) analysis was then applied to further look into the function of EMP1. Finally, Cellchat was employed to examine the interaction between macrophages and epithelium cells.
Results: The results showed that higher M1/M2 ratio was found to be associated with a more favorable prognosis of bladder cancer. EMP1 was identified to be the key gene indicative of M1/M2 ratio and higher EMP1 expression was associated with poor prognosis. Further analyses showed that EMP1 might promote tumor invasion and metastasis via epithelial-mesenchymal transition (EMT) and focal adhesion (FA). Moreover, the expression level of EMP1 could serve as an indicator of immunotherapy efficacy. The scRNA-seq data indicated that EMP1 in cancer cells was strongly associated with tumor proliferation. Finally, the Cellchat results exhibited that EMP1 might promote the interaction between macrophages and cancer cells through the fibronectin 1-syndecan 1 (FN1-SDC1) pathway.
Conclusion: Our study identified EMP1, an M1/M2-related gene, the expression of which may act as a prognostic indicator for the proliferation, metastasis, and response to immunotherapy. EMP1 might be involved in the regulation on M1/M2 ratio.
Keywords: bladder cancer; EMP1; EMT; FN1-SDC1; macrophage

Introduction

Bladder cancer is the second mortal urinary system cancer across the globe[1,2]. For patients with localized and non-muscle-invasive bladder cancer (NMIBC), surgery or radiotherapy is currently the standard therapeutic approach. However, these treatments are often ineffective in preventing recurrent or distant metastatic diseases, especially in the case of bladder cancer, which typically forms microscopic metastases.
Macrophage 1 (M1) and Macrophage 2 (M2) are two types of macrophages that differentiate from monocytes under different stimuli[3]. In tumors such as breast cancer, M1 and M2 macrophages have been shown to have opposite effects on tumor growth[4]. M1 macrophages are the predominant phenotype in normal immunological responses, while M2 macrophages induce immunosuppression and tumorigenesis[5]. Although the prognostic value of M1 macrophages has been established, M2 macrophages have also been extensively investigated, but with inconsistent results. Researchers have recognized the limitation of analyzing the effect of M1 or M2 macrophages separately and have attempted to combine them as an M1/M2 macrophage ratio. In tumors such as oral squamous cell carcinoma[6], some researchers have already proven the strong prognostic value of the M1/M2 ratio. However, in bladder cancer, the value has not yet been clearly investigated.
Epithelial membrane protein 1 (EMP1) is a member of the epithelial membrane protein family (including EMP1, EMP2, and EMP3) and is encoded by the same gene family as growth arrest-specific 3 (GAS3) and peripheral myelin protein 22 kDa (PMP22). It is an important target for tumor therapy [7]. Previous research has confirmed that EMP1 is an important prognostic factor that is closely related to tumor metastasis and proliferation [8-10]. However, the relationship between EMP1 and tumors, especially bladder cancer, has not been fully studied and needs further exploration. Using Cox regression and SurvivalRandomForest to perform dimension reduction, we found that EMP1 is associated with poor prognosis in bladder cancer and further explored its relationship with the M1/M2 macrophage ratio and bladder cancer.
In this study, we identified EMP1 as the most valuable gene related to the M1/M2 macrophage ratio through several bioinformatic approaches, and we performed its clinical, biological and multi-omic characterization in both public databases and our own patient cohort. Furthermore, this study revealed the essential role of EMP1 in predicting the response to immunotherapy or serving as a biomarker for bladder cancer. Finally, we analyzed scRNA-seq data from the Gene Expression Omnibus (GEO) database, in which our finding demonstrated EMP1’s significant impact on the communication between myeloid cells and other cell types.

Materials and Methods

Data of Bladder Cancer Gathering and Processing

Bladder cancer samples with adequate clinical information were gathered from The Cancer Genome Atlas (TCGA) database and Gene-Expression Omnibus (GEO). In terms of the TCGA cohort (403 bladder cancer samples), the RNA-seq data and accompanying clinical information were extracted from the TCGA database (http://cancergenome.nih.gov/), which was then transformed into transcripts per kilobase million (TPM). Four GEO microarray cohorts were employed in this study, and the expression and survival information were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), with background adjustment, and were normalized using the RMA algorithm. Before further analysis, all gene expression data were log2 converted and quantile normalized utilizing the ‘normalizeBetweenArrays’ technique in the R package limma 3.46.0. For the analysis using merged gene expression data from different data sets, we removed the batch effects using the R package sva 3.36.0.

Inference of the Immune Infiltrating Microenvironment

To evaluate the abundance of M1 and M2 cells, we employed the Cibersortx algorithm (https://cibersortx.stanford.edu/), an analytical tool to impute gene expression profiles and provide an estimation of the abundances of member cell types in a mixed cell population using gene expression data. Other algorithms, including TIMER, CIBERSORT, MCPCounter, XCell, and EPIC deconvolution algorithms, were also implemented in immunedeconv R package 0.35.

Unsupervised Clustering for M1M2 ratio-related genes

With the assistance of the R package ConsensuClusterPlus 1.54, it was possible to find robust M1M2 ratio-related clusters in TCGA patients by using a consensus clustering technique of partition (based on the Euclidean distance and Ward's linkage) based on the M1M2 related genes. The cumulative distribution function (CDF) and consensus heatmap were used to determine the optimal K value. This method was performed 1000 times to ensure the stability of the stratification process.

Dimension reduction for discovering core gene

The TCGA data were subjected to univariate Cox regression analysis to identify the DEGs associated with any prognoses with a P-value < 0.01. SurvivalRandomForest with 1000 trees was used to validate the results and rank the importance of the 11 genes obtained by Lasso regression in the R package SurvivalRandomForest version 3.6.4.

Gene Set Variation Analysis (GSVA)

We downloaded canonical pathways gene sets from the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes (KEGG) databases. The GSVA score was then measured for each bladder cancer sample using the GSVA algorithm in the R package GSVA 1.36.2.

Gene Set Enrichment Analysis (GSEA)

The GSEA algorithm analyzed the enriched biological processes between different groups. The data in TCGA were first Transformed to be prepared for Linear Modelling by voom in R package limma 3.46.0. The differential genes between the two groups were calculated with the R package limma 3.46.0. Subsequently, they were pre-ranked by log2 fold-change and delivered to the R package clusterProfiler 3.18.1 for GSEA analysis. The results with a p.adjust < 0.05 were recognized as statistically significant.

Multi-Omics Features of EMP1 high and low groups

Somatic mutations and somatic copy number alternations (CNAs), which corresponded to the cases with RNA-seq data, were downloaded from the TCGA database. GISTIC analysis was performed online (https://www.genepattern.org/) to determine the genomic event enrichment. CNAs associated with the two groups and the threshold copy number at alteration peaks were obtained using GISTIC 2.0 analysis.

Prediction of Immunotherapy Response

The IMvigor210 cohort, a urothelial carcinoma cohort treated with the anti‐PD‐L1 antibody atezolizumab, was used to predict patients’ responses to immunotherapy. Based on the Creative Commons 3.0 License, complete expression and clinical data were downloaded from http://research-pub. Gene.com/IMvigor210CoreBiologies. Raw data were then normalized using the limma R package, and the count value was transformed into TPM value.

Single-cell transcriptomic analysis

Single-cell RNA sequencing (scRNA-seq) gene expression data from 11 samples (GSM4006644, GSM4006645, GSM4006646, GSM4006647, GSM4006648, GSM4751267, GSM4751268, GSM4307111, GSM3729178, GSM3729179, GSM5655079) were downloaded from published study via Gene Expression Omnibus, through accession number GSE135337, GSE145137, GSE130001 and GSE186520. we filtered cells with a unique feature count of > 2500 or < 200 and cells with a mitochondrial count of > 5%. Then, the default parameters of the “NormalizeData” function of Seurat were used to normalize the feature-expression measurements for every cell by the total expression. Finally, 14348 cells were introduced into a combined Seurat object via the “FindIntegrationAnchors” and “IntegrateData” functions. Then, variable genes were carried forward into scaling and principal component (PC) analysis. Significant PCs (top 15) were used for t-SNE analysis and clustering through “RunTSNE” (perplexity=100, max iteration=3000) and “FindClusters” functions (resolution = 0.5). To identify cell types, we checked whether the well-studied marker genes were the top differentially expressed genes (DEGs) and annotated the most probable identity for each cell cluster. The remaining cell types were then identified manually from the cell marker database (http://biocc.hrbmu.edu.cn/CellMarker/).

Differential Expression and Analysis of Signaling Pathways

DEGs among clusters were detected by the Seurat function “FindAllMarkers”. Volcano plots were applied to visualize the genes with upregulated or downregulated expression. Genes with Benjamini–Hochbergadjusted P < 0.05 and absolute log2 fold change between two groups > two were used to analyze functional enrichment using the GO and KEGG databases (clusterProfiler 3.18.1). Additionally, ranked gene set enrichment analysis (GSEA) was undertaken. Genes were ranked based on their phenotypes, and the GSEA algorithm proposed by Subramanian et al. [11] was used to compute the enrichment score of each gene.

Ligand–Receptor Expression and Cell Interactions

Cell-to-cell communication (CellChat 0.0.2; R package) was ascertained by evaluating the expression of pairs of ligands and receptors within CellChatDB. human. We examined the interaction between different cell types, and communication with many cells of less than ten was filtered.

Statistical Analysis

Unpaired Student t-tests and Wilcoxon rank-sum tests were applied to evaluate the statistical significance of normally-distributed and non-normally distributed variables, respectively when comparing the two groups. Kruskal-Wallis tests and One-way analysis of variance (ANOVA) were also used to conduct difference comparisons of more than two groups. Spearman and distance correlational analyses were conducted using the R package Hmisc 4.4.1. Objects with a Spearman correlation coefficient greater than 0.4 were strongly correlated. The hazard ratios (HR) of all prognostic factors were computed using a univariate Cox proportional hazards regression model. The “surv_cutpoint” function of the R package survminer 0.4.8 was used to estimate the optimal cut-off point for prognostic factors according to their relationship with the patients’ survival probability with the maximum rank statistic. Next, survival curves were drawn using the Kaplan–Meier method. The statistical significance of the difference in survival probability was estimated using the log-rank test. The R package forestplot 1.10 was employed to demonstrate the univariate prognosis analyses of different groups of prognostic factors. The networks of canonical signaling pathways were constructed using Cytoscape 3.7.2, and the hub pathway was estimated by the Cytoscape plug-in CytoHubba. The R package ComplexHeatmap 2.4.3 generated all the heat maps. A Waterfall Chart was applied to exhibit the overview of the gene mutation landscape, which was generated using the R package maftools 2.4.12. The statistical difference of CNVs between the EMP1 low and EMP1 high groups was determined using the R package cnvaq 0.1.3. Then, the IGV 2.8.2 software was employed to visualize the CNV landscape of the two groups. All the above analyses were performed using the R 4.0.0 software. All the statistical analyses were two-sided, and a P < 0.05 was considered statistically significant.

Specimens and cell culture

This research was approved by the Ethics Committee. Specimens (N=32) were histopathologically confirmed by three independent pathologists and derived with informed consent. We recruited 32 biopsied bladder cancer patients and their tumor excisions were placed in sterile 50 mL centrifuge tubes containing RPMI1640 media (Invitrogen/GIBCO, Carlsbad, CA), 10% heated-inactivated fetal bovine serum (FBS, Gibco), and 1% antibiotics (100 U/mL penicillin, 100 µg/mL streptomycin). The centrifuge tubes were transported from the operating room to the laboratory or stored in liquid nitrogen.
HT1197 and RT112 cell lines were obtained from the American Type Culture Collection (ATCC, Manassas, VA). Cells were cultured in RPMI1640 media (Invitrogen/GIBCO, Carlsbad, CA) supplemented with 10% heated-inactivated fetal bovine serum (FBS, Gibco) and 1% antibiotics (100 U/mL penicillin, 100 µg/mL streptomycin). All cell cultures were incubated at 37℃ in a humidified environment of 5% CO2 in the air.

Real-time quantitative RT-PCR

The total RNA of patients was extracted by using a RNeasy kit (Qiagen Sciences, Hilden, Germany) and using a PrimeScript II 1st-strand cDNA synthesis kit (TaKaRa, Dalian, China) to obtain cDNA. qPCR was performed using a LightCycler 480 real-time PCR machine (Roche, Basel, Switzerland) with SYBR and GAPDH (glyceraldehyde-3-phosphate dehydrogenase) as an internal control. The resulting data from the PCR machine were all analyzed based on the cycle threshold (DDCT) method, and all primers were listed in supplementary Table S8.

Western blots

Cells in six-well plates were lysed with NP-40 buffer (Beyotime, Shanghai, China) with 0.1% mM phenylmethylsulfonyl fluoride (PMSF) protease inhibitor cocktail (Roche Molecular Biochemicals) for 30 min in ice. SDS was introduced to the supernatant and boiled for 10 min after centrifugation at 12,000 ×g for 10 min at 4°C. Before Western blotting, protein concentrations were determined by BCA assay. Cell lysates were separated by SDS-PAGE and transferred to nitrocellulose filter (NC) membranes (Millipore, Billerica, MA), then blocked with 5% nonfat dry milk in TBST (20 mM Tris, pH 7.4, 150 mM NaCl, 0.1% Tween 20) for 2 h at room temperature, and then incubated with anti-EMP1 (ab230445, Abcam, Cambridge, MA) antibody or anti-GAPDH (LF205, Epizyme, Shanghai, China) antibody at 4°C overnight. After three rounds of washing with TBST, the membrane was incubated with a secondary antibody in PBS for 1 h at room temperature in the dark. After using TBST for washing three times, the membrane was scanned and analyzed in an Odyssey infrared imaging system (Li-Cor Biosciences, Lincoln, Nebraska USA).

Flow cytometry

Tumor excisions from patients were cut into small pieces and homogenized in a cold staining buffer to make single-cell suspensions in the presence of collagenase. The cells were stained with fluorescence-labeled FITC-CD11b antibody, APC-CD68 antibody, PE/Cy7-CD80 antibody, and PE/Cy7-CD206 antibody (Biolegend, San Diego, CA). Stained cells were measured on an FC500 flow cytometer (Beckman, Indianapolis, Indiana USA) and analyzed by Kaluza software (Beckman).

RESULTS

The predictive value of M1M2 macrophages and the M1M2 ratio correlated genes

To evaluate the prognostic value of the M1M2 ratio in bladder cancer, we retrieved the TCGA data from the GDC data portal (https://portal.gdc.cancer.gov/) and excluded samples with purity less than 60%. The Cibersortx algorithm was utilized online (https://cibersortx.stanford.edu/) to estimate the cell infiltration (Table S1). We then filtered out those samples in which the quantity of Macrophage 1 or 2 was zero. After that, we observed that the survival curves were separated by the optimal cutoff value of the M1M2 ratio (Fig. 1A). Subsequently, experiments were performed to validate the association between the M1M2 ratio and survival in patients with bladder cancer. Tumor resections were collected from 32 patients diagnosed with bladder cancer, and their clinical information was listed (Table S2). All of the patients underwent surgery, and there was no distant metastasis. The median follow-up time was 1173 days. Tumor tissues were stained for CD11b, CD68, CD80 (M1), and CD206 (M2) with flow cytometric analysis to determine the M1M2 ratio (Fig. 1B). Combining tumor resection flow cytometry results and survival data from 32 patients, we found that patients with a higher M1M2 ratio separated by the best cut-off value had better survival than patients with a lower ratio (Fig. 1C) which was consistent with TCGA database analysis result. To investigate the inner significance of the M1M2 ratio, we identified genes that were highly correlated to the M1M2 ratio (|r|>0.3, Spearman correlation, Table S3). These genes clustered the samples into groups through an unsupervised clustering algorithm. Based on the results of a consensus matrixes heatmap (Fig. 1D, S1A‒C), the Delta area curve, tracking plot, and cumulative distribution function (CDF, Fig. S1D‒F), the optimal K value (the number of groups) was 3, and the assignment of each sample was listed (Table S4). A total of 1389 intersected differentially expressed genes (DEGs) among these three clusters were detected by R package limma 3.46.0 (Figure 1E, Table S5). As illustrated in the networks of canonical signaling pathways visualized by clueGO (a Cytoscape plug-in), these DEGs were highly enriched in Myeloid leukocyte migration (Fig. 1F).
Figure 1. The association between M1M2 ratio and bladder cancer prognosis. (A) The survival curve of 271 patients from TCGA was divided into two groups with high and low ratios based on their optimal cut-off value. (B) Results from three representatives of 32 bladder cancer patients depicted in cytometry plots and the relative quantification of M2-like macrophages (CD206) and M1-like macrophages (CD80) gating on CD11b+CD68+ cells. (C) Survival data of 32 patients divided into two groups with high and low ratios using the best cut-off value. (D) Consensus matrixes of TCGA cohort for k = 3. (E) Venn plot generated to show DEGs among these three clusters. (F) Cytoscape plug-in ClueGo exhibiting the GO clustering result. (G) The relation between error rate and the number of trees running the SurvivalRandomForest algorithm. (H) The relative importance of 17 genes selected by dimension reduction. (I) Survival curve of 403 TCGA patients categorized into two groups with high and low ratios based on the median level of EMP1. (J, K) Forest plots depicting the result of univariate and multivariate Cox regression analysis of these eight genes.
To identify the most prognostic gene, we first performed univariate COX regression analysis with the survival status and time of patients as the dependent variables, and filtered out 253 genes (P < 0.01, Table S6). Then, we applied SurvivalRandomForest to further reduce the dimensionality based on the survival of patients. We set the total number of trees to 1000 and obtained the lowest error rate when there were 194 trees (Fig. 1G, S1H). This resulted in 17 genes that were the most predictive of patients’ survival (Table S7). The actual variable importance of these genes is shown in Figure 1H. We selected the top eight genes for further investigation (SERPINB7, ANXA1, CCL26, MYO5A, DUSP14, SERPINB2, CFL2, and EMP1). Survival analysis revealed a significant difference in long-term survival between patients stratified by the median level of these eight genes (Fig. 1I, S2A‒D). The hazard ratio (HR) for each of these genes was higher than 1.0, with a P-value lower than 0.05 by univariate COX regression analysis (Fig. 1J). However, multi-COX regression analysis indicated that only EMP1 was associated with survival when combining all the genes (Fig. 1K). Therefore, we chose EMP1 as the key gene for subsequent analysis.

The biological value of EMP1

We first validated the correlation between the number of macrophages and EMP1. EMP1 was moderately correlated with the M1M2 ratio and M2, and highly related with M1 (Fig. 2A‒C), while there was no correlation between EMP1 and CD8+T cells (Fig. S3A). The correlation between EMP1 and the markers level of M1 and M2 was also presented (Fig. S3B). GSVA was applied to score each sample's enrichment of macrophage-related pathways. EMP1 was highly correlated with macrophage differentiation and myeloid differentiation (Fig. 2D). Regarding the significant impact of EMP1 on M1 and M2, six algorithms, including TIMER, CIBERSORT, MCPCounter, XCell, EPIC, and ssGSEA, were implemented to comprehensively compute the relation between EMP1 expression and immune cell infiltration via the R package immunedeconv 0.35 (Fig. 2E). It was found that increased cells such as fibroblasts, DCs, Tregs, and Monocytes were related to EMP1, with worse survival probability. Then GSEA analysis was performed in the Hallmark, and KEGG gene was set to investigate the biological function of EMP1 (Fig. 2F‒G). Consequently, EMP1 was testified to positively correlate with epithelial-mesenchymal transition (EMT) and focal adhesion (FA), which can be an essential reason for EMP1 inducing metastasis of bladder cancer. Previous research has confirmed that EMT and focal adhesion may be related to tumor metastasis [12-15].
The results of qPCR on tumor tissues and adjacent normal tissues of 32 patients showed that the expression of EMP1 in normal bladder tissues was higher than that in tumor tissues (Fig. 2H). The same results were also observed in cell experiments. The results based on Western blotting showed that the expression of EMP1 in the normal bladder cell line (SV-HUC-1) was higher than that of seven bladder tumor cell lines (HT-1197, RT-112, 647-V, J82, EJ, T24 and UM-UC-3) (Fig. 2I). This is an interesting phenomenon, so next we continue to analyze the clinical value of EMP1.

The clinical value and multi-omics characterization of EMP1

To evaluate the clinical value of EMP1, we collected five cohorts and generated a forest plot. EMP1 showed great predictive value in three out of five cohorts, with P-values lower than 0.05 (Fig. 3A). For GSE32894, the largest cohort, the HR value was more than 2, and the Kaplan-Meier curve revealed a significant difference between the high and low EMP1 groups, separated by the median level of EMP1 (Fig. 3B). To validate EMP1’s clinical value, we collected frozen bladder cancer tissues from 32 patients with complete clinical information (Table S3). We also combined the RT-PCR results with patient survival data and found that patients with high EMP1 expression had a significantly lower long-term survival rate than those in the low group, separated by the median level of EMP1 (Fig. S3C). To exclude other factors affecting survival, we applied multivariate Cox regression in the merged cohort and our validation cohort (Fig. 3C-D). We concluded that EMP1 could be an independent risk factor for bladder cancer.
To further investigate the genomic heterogeneity of patients with different EMP1 levels, we compared the tumor mutation burden (TMB) and the mutational signatures of the EMP1 high and low groups, separated by the median level of EMP1. We found no significant difference in TMB between the two groups (Fig. 4A). We identified four mutational signatures that were highly correlated with bladder cancer, namely, SBS1 (age-related), SBS2 and SBS13 (APOBEC activity-related), and SBS5 (ERCC2 mutation-related). The EMP1 low group had significantly more mutations in APOBEC-related signatures (SBS2 + SBS13), while the EMP1 high group had more SBS1 and SBS5 signatures (Fig. 4A). Among the frequently mutated genes, the EMP1 low group harbored significantly more mutations of KDM6A, FGFR3, and SYNE1, while the PIK3CA mutation was more enriched in the EMP1 high group (Fig. 4A). Human chromosome 17p11.2, the region of the EGFR gene, showed more amplification in the EMP1 low group. Moreover, the EMP1 high group had more EMP1 copy number gain, and the EMP1 low group had more EMP1 copy number loss, corresponding to the EMP1 expression level. The copy number variation of macrophage markers such as IFNG, CD68, and CD163M also exhibited a significant difference between the EMP1 high and low groups (Fig. 4A).
Finally, we performed GSVA analysis in the meta-cohort to score the enrichment of pathways from the GO and KEGG databases for each sample. The volcano plot shows the most differentially enriched pathways between the EMP1 low and high groups (Fig. 4B), including CXCR chemokine receptor binding, antigen processing and presentation, and MHC II protein complex, which are essential for M1M2 function and chemotaxis. The heatmap shows that the EMP1 high group is active in myeloid cell differentiation, proliferation, and regulation (Fig. 4C).
Figure 2. EMP1 deficiency inhibits bladder cancer cell migration and invasion, whereas overexpression promotes this process. (A‒C) Correlation between EMP1 and levels of M1M2 ratio, M1, and M2, respectively. (D) Correlation matrix exhibiting the association between EMP1 and Myeloid cell development pathway. (E) Heatmap for the infiltration of immune cells with rows representing different kinds of immune cells and columns representing different samples with ECMI increasing from left to right. (F, G) GSEA running plot demonstrating the biological process of hallmark and KEGG enriched in EMP1 high and low. (H) qPCR of tumor tissues and adjacent normal tissues of 32 patients. (I) Western blot analysis of EMP1 expression in the normal bladder epithelial cell line SV-HUC-1 and seven bladder cancer cell lines (including HT-1197, RT-112, 647-V, J82, EJ, T24, and UM-UC-3). * P < 0.05.
Figure 3. EMP1 is the most prominent gene associated with prognosis. (A) Forestplot was plotted to present the predictive value of EMP1 in each cohort. (B) Kaplan-Meier survival analysis of EMP1 high group and EMP1 low group in GSE32894 cohort. (C, D) Forestplot showing the independent predictive value of EMP1 when combing clinical information in Merge-cohort and our validation cohort.

The efficacy of EMP1 in predicting the response of immunotherapy

Since EMP1 is related to immune infiltration and immune-relevant pathways, we investigated its potential value in predicting the response to immunotherapy. We collected data from IMvigor210, in which mice received anti-PD-1. There was a significant difference in EMP1 expression between the anti-PD-1 response and non-response groups, with a Wilcoxon test P-value less than 0.05 (Fig. 4D). More mice in the EMP1 high group achieved a good effect of immunotherapy than in the EMP1 low group (Fig. 4E). Moreover, the mice in the EMP1 high group had a significantly longer survival time than those in the EMP1 low group (Fig. 4F). This phenomenon may be attributed to the different expression of immune checkpoint molecules. As shown in Figure 4G, most of the immune-suppressive IC molecules, such as CTLA-4, PD-1, and ICOS, were more expressed in the EMP1 high group. However, the trend of TNFRSF14 was the opposite, and further research is needed to investigate the underlying mechanism.

Single-cell analysis of the role of EMP1 in bladder cancer

To investigate the role of EMP1 at a higher resolution level, we collected 11 samples (GSM4006644, GSM4006645, GSM4006646, GSM4006647, GSM4006648, GSM4751267, GSM4751268, GSM4307111, GSM3729178, GSM3729179, GSM5655079) that underwent 10×genomics standard single-cell transcriptome analysis. We filtered out cells with a unique feature count of > 2500 or < 200 and cells with a mitochondrial count of > 5%. After normalization and scaling, we combined 14348 cells into a Seurat object. We then projected the single-cell data into two-dimensional space using uniform manifold approximation and projection, including initial PCA analysis (Fig. S4A), and labeled 13 clusters with different cell types. We extracted the DEGs between each cluster with p-value-adjusted < 0.05 and log2FC > 0.5 (Fig. 5A). To identify cell types, we checked whether the well-studied marker genes were among the top DEGs and annotated the most probable identity for each cell cluster. We identified the remaining cell types by manually retrieving them from the cell marker database (http://biocc.hrbmu.edu.cn/CellMarker/, Fig. 5B). The tSNE plot showed the distribution of cells from each GEO data set (Fig. 5C). We presented the abundance of M1, M2, and Myeloid cells and annotated them according to their surface markers (CD14, CD274, CD163, MbSR1, etc., Fig. 5D-E). Interestingly, EMP1 was moderately expressed in every cell type, without a significant difference (Fig. S4B).
Figure 4. Genomic alteration landscape according to EMP1 and the predictive value of EMP1 for immunotherapy response. (A) from the top to the bottom panels showing tumor mutation burden (TMB), the relative contribution of four mutational signatures, selected differentially mutated genes (> 5%) and broad-level copy number alterations (> 20%), and specific genes. The proportion of patients in each alteration is displayed in the right bar charts. (B) Volcano plot exhibiting the differentially enriched GO pathways between EMP1 high and low groups in merged-cohort. (C) Heatmap illustrating the GSVA score of macrophage-related pathways between EMP1 high and low groups. (D) Violin plot presents the levels of EMP1 between EMP1 high and low groups. (E) Stacked bar chart showing the proportion of immunotherapeutic response in different EMP1 groups. (F) Kaplan-Meier survival analysis of EMP1 high group and EMP1 low group in IMvigor210 cohort. (G) The boxplot of the ten immune checkpoint levels in EMP1 high and low groups.
Figure 5. Single-cell transcriptomic analysis reveals the transcriptome of cells in the tumor microenvironment of bladder cancer samples. (A) The differentially expressed genes among each cluster. (B–C) tSNE visualization of reference cell types and 11 GEO data sets. (D) tSNE plot exhibiting the expression of macrophage-related markers on myeloid linage cells. (E) The composition of myeloid linage cells shown on tSNE graph. (F–G) GSEA and GO analysis elucidating EMP1 relevant pathways.
Since EMP1 is commonly expressed in most cells (Fig. S4C), we focused on EMP1 in epithelial cells[16]. We performed GSEA analysis to explore the EMP1 relevant pathways and GO analysis to extract the enriched pathways of DEGs between EMP1 high and low cells. We found that EMP1 was not only highly associated with extracellular matrix (ECM) receptor interaction and focal adhesion (Fig. 5F), which corresponded to our results in bulk-RNA-seq, but also enriched in epithelial cell proliferation-related pathways (Fig. 5G).

Cellchat reveals the inner cell communication pattern

Since EMP1 is closely related to the M1M2 ratio, we focused on the cell communication pattern among various types of cells. The GSEA results showed that EMP1 was closely correlated with EMT and other biological processes related to ECM components (Fig. 5F-G). Therefore, we explored cell communication based on ECM signaling pathways containing Collagen, FN1, LAMININ, THBS, and VTN networks. We used the R package CellChat 0.0.2 to generate the general interaction number and strength, showing that Stem cells and Myeloid cells were actively communicating with other cells (Fig. 6A-B), and the interaction between M1 and M2 with Stem cells was also strong (Fig. 6C-D). Next, we depicted the distribution of specific signaling axes from the five ECM signaling pathways among different cell-cell communications (Fig. 6E). Surprisingly, in the interactions involving Myeloid lineage cells, the FN1-SDC1 axis was highly active, indicating that Myeloid cells may have an impact on cancer cells. We plotted the cell-cell communication via FN1-SDC1 in a chord graph (Fig. 6F). In the FN1-SDC1 signaling network, we found that the signal sent from epithelial cells lacked interaction among tumor cells, while immune cells significantly affected various epithelial cells (Fig. 6G). Correspondingly, we presented the contribution of each L-R pair in FN1, Collagen pathway, and LAMININ (Fig. 6H, S5A-B). A violin plot showed the expression status of all the receptors and ligands involved in the FN1 signaling pathways. FN1 and CD44 were highly expressed in myeloid cells, while SDC-1 was mainly expressed in epithelial cells (Fig. 6I). To further investigate the ECM communication pattern of different types of cells, we performed hierarchical clustering. The most suitable number of groups was three, based on Silhouette metrics (Fig. 6J). Pattern 3, in which VTN and FN1 pathways played a dominant role, was more active in B cells, Myeloid cells, and Macrophage cells (Fig. 6K).

The distinct cell communication between EMP1 high and low group

To investigate the impact of EMP1 on cell-cell interaction, we divided samples into EMP1 high and low groups, based on the average EMP1 level of epithelial cells. The EMP1 low group had lower interaction number and strength than the EMP1 high group (Fig. 7A). The EMP1 high group enhanced the internal communication of epithelial cells and weakened the communication between Myeloid lineage cells (Fig. 7A). Among the various types of ECM signaling pathways, all pathways except LAMININ were more active in the EMP1 high group than in the EMP1 low group (Fig. 7B) [17]. The circle plot shows the significantly lower interaction number and strength among cells in the EMP1 low group (Fig. 7C-D). Considering the importance of the FN1 pathway, as previously mentioned, we focused on the difference in FN1 activity between the EMP1 high and low groups. The EMP1 low group had a much weaker interaction among cells than the EMP1 high group (Fig. 7E). Moreover, the FN1 signaling was significantly more active in the EMP1 high group, which indicates that EMP1 may increase the efficacy of M2 through the FN axis (Fig. 7F). We also compared the activity of the other four pathways between the EMP1 high and low groups, and the results are presented in supplementary figures (Fig. S6A-D).

DISCUSSION

We innovatively started our research from the M1/M2 ratio, which reflects the balance of macrophages in the tumor microenvironment (TME). Macrophages can be divided into M1 type, which promotes tumor regression, and M2 type, which promotes tumorigenesis[18]. Previous studies have confirmed that M1 or M2 macrophages play a role in the occurrence and prognosis of bladder cancer[19-21], but no research has combined the two to further explore the mechanism in bladder cancer. In our study, we first combined M1 and M2 macrophages in a public database and patient flow cytometry results to explore their association with bladder cancer. To gain insights into the mechanisms of bladder cancer progression and how M1M2 macrophages function within them, we derived the most survival-predictive gene associated with M1M2 macrophages, EMP1, through dimensionality reduction using SurvivalRandomForest. We found that EMP1 was associated with epithelial-mesenchymal transition (EMT) and focal adhesion by bioinformatics analysis, suggesting that EMP1 is related to tumor metastasis and invasion[12-15]. We further investigated the role of EMP1 in epithelial cells and its biological function at a higher resolution level using single-cell transcriptome analysis. We found that EMP1 was associated with extracellular matrix (ECM) receptor interaction and focal adhesion [14,15,22,23] and was also enriched in epithelial cell proliferation-related pathways. This corresponded to our results in bulk-RNA-seq, demonstrating the role of EMP1 in metastasis and proliferation. We also explored the relationship between EMP1 and M1M2 cells using Cellchat and found that myeloid cells mainly interacted with epithelial cells through the FN1-SDC1 pathway, which has been confirmed to be related to tumors in other types of cancers[24-28]. Cells with high expression of EMP1 had significantly more active FN1-SDC1, and more active macrophage and epithelial cells, indicating that EMP1 may regulate the biological process of epithelial cells through FN1-SDC1.
One limitation of this study is the discrepancy between the WB results and the qPCR results and survival analysis. The WB and qPCR results showed that the normal bladder cells had higher EMP1 content than the bladder cancer cells such as RT-112, which seemed to indicate that EMP1 might have a negative effect on cancer occurrence. However, the survival analysis showed that patients with high EMP1 content had a worse prognosis. This may be due to the fact that some genes play different roles in different stages of cancer development, or in tumorigenesis and metastasis, such as TGF-β. TGF-β acts as a tumor suppressor early in the carcinogenesis process but turns into a tumor promoter at a later stage [29]. We speculate that EMP1 may have a similar dual role, showing a negative effect during the transformation of normal cells into tumors, and a positive effect on the metastasis and proliferation of bladder cancer after tumor formation. To address this problem, our follow-up research will continue to explore the specific mechanism of EMP1.
EMP1 may have potential value in clinical prediction, for example, patients with high EMP1 in tumor tissue may have more metastases. In addition, EMP1 may also have great potential in clinical treatment. Immunotherapy is one of the hopes of curing cancer, but it also faces many challenges, such as personalized precision treatment[30]. Studies have confirmed that macrophages can be used as therapeutic targets and play a role in tumor immunotherapy [31]. Currently, the approved immune checkpoint inhibitors (ICIs) in bladder cancer are mainly PD-1/PD-L1 related drugs, such as atezolizumab, pembrolizumab, nivolumab, etc. [32]. By comparing the response of EMP1 high and low groups to immunotherapy, we found that patients with high EMP1 expression had a poor immunotherapy effect and survival. Interestingly, the EMP1 high expression group had a high expression of immune checkpoint proteins such as CTLA4, PD-1, and ICOS, except for the TNF receptor family member TNFRSF14. Therefore, EMP1 may become a potential biomarker for checkpoints and personalized precision therapy strategies such as anti-PD-1, CTLA4, and ICOS, while a low expression of EMP1 may guide anti-TNF therapy.
Figure 6. CellChat algorithm revealing macrophage communication characteristics in bladder cancer. (A–B) Interaction number and strength among all types of cells. (C–D) Interaction number among macrophages and the other cells. (E) Summary of ligand-receptor interactions between different cell clusters via ECM-related signaling pathways. (F) Chord plot demonstrating the interaction between different cell types via FN1 signaling network. (G) Connection plot illustrating the interaction between cells in TME. (H) Bar plot exhibiting the relative contribution ligand-receptor in FN1 regarding pathways. (I) Violin plots presenting FN1-related protein levels on different cell types. (J) Line chart displaying the measure score in a different number of communication patterns. (K) Heatmap demonstrating the different cell communication patterns in different cells and the pathway composition of each pattern.
Figure 7. Cell communications difference between EMP1 high and low groups. (A) Bar plot depicting different interaction numbers and strengths between EMP1 high and low groups. (B) Bar plot illustrating the relative and absolute pathway activity comparison between EMP1 high and low groups. (C-E) Circle plot presenting the interaction difference between EMP1 high and low groups in the whole picture or focusing FN1 pathway. (F) Heatmap displaying the cell-cell communication difference via FN1 signaling pathway.
The association between EMP1 and the FN1-SDC1 axis was revealed by single-cell analysis and Cellchat. The underlying mechanism of this interaction remains to be elucidated. Therefore, the next step of this study is to conduct animal in vivo experiments to validate the findings from the in vitro assays. Furthermore, the clinical significance of EMP1 in the precision and personalized treatment of bladder cancer will be explored.
In conclusion, the M1/M2 ratio was identified as a negative prognostic marker of bladder cancer through TCGA and flow cytometry analysis. To investigate the underlying mechanism, the most relevant and significant gene related to the M1/M2 ratio, EMP1, was discovered. EMP1 was found to be highly associated with EMT, FA, and proliferation through TCGA and experimental data. Moreover, EMP1 could also serve as a predictor of immunotherapy response. Furthermore, scRNA revealed a widespread expression of EMP1 in different cell types. Focusing on the epithelial cells, GSEA analysis demonstrated a positive correlation between EMP1 and bladder cancer, EMT, and ECM. The Cellchat result indicated that M1/M2 mainly interacted with epithelial cells through the FN1-SDC1 axis, which was more active in the EMP1 high group. Taken together, these data suggested that EMP1 might be a potential prognostic marker for predicting bladder cancer proliferation, metastasis, and response to immunotherapy, and that EMP1 might modulate tumor proliferation and metastasis and improve the prognosis of cancer patients by regulating the M1/M2 ratio.

Ethics approval and informed consent

This study was approved by the Investigation and Ethics Committee at Harbin Medical University. Patients gave their informed consent to be subjected to the protocol.

Funding

This research was supported by National Natural Science Foundation of China (U20A20385).

Acknowledgments

We would like to express our sincere appreciation for the constructive suggestions provided by our colleagues in Harbin Medical University Cancer Hospital. Their feedback helped us improve the quality of our research. We are also grateful to the data analysis platform Bigpt (biogpt.org.cn) for their data support, which made our research more comprehensive and accurate.

Author contributions

Study design: Wanhai Xu, Jinqiao Li
Biology experiment: Jinqiao Li, Honglei Wang
Statistical analysis: Jinqiao Li, Jianyu Liu
Data collection: Jinpeng Ma, Yueze Wang
Drafting and editing of the manuscript: Jinqiao Li, Jianyu Liu, Honglei Wang

References

  1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics. CA Cancer J Clin. 2022;72(1):7-33. Epub. 2022;20220112: https://doi.org/10.3322/caac.21708.
  2. Global Burden of Disease Cancer C, Kocarnik JM, Compton K, Dean FE, Fu W, Gaw BL, et al. Cancer Incidence, Mortality, Years of Life Lost, Years Lived With Disability, and Disability-Adjusted Life Years for 29 Cancer Groups From 2010 to 2019: A Systematic Analysis for the Global Burden of Disease Study 2019. JAMA Oncol. 2022;8(3):420-44. https://doi.org/10.1001/jamaoncol.2021.698 PMID: 34967848
  3. Yang Q, Guo N, Zhou Y, Chen J, Wei Q, Han M. The role of tumor-associated macrophages (TAMs) in tumor progression and relevant advance in targeted therapy. Acta Pharm Sin B. 2020;10(11):2156-70. Epub 20200419. https://doi.org/10.1016/j.apsb.2020.04.004 PMID: 33304783
  4. De Palma M, Biziato D, Petrova TV. Microenvironmental regulation of tumour angiogenesis. Nat Rev Cancer. 2017;17(8):457-74. Epub 20170714. https://doi.org/10.1038/nrc.2017.51 PMID: 28706266
  5. Mohapatra S, Pioppini C, Ozpolat B, Calin GA. Non-coding RNAs regulation of macrophage polarization in cancer. Mol Cancer. 2021;20(1):24. Epub 20210201. https://doi.org/10.1186/s12943-021-01313-x PMID: 33522932
  6. Dan H, Liu S, Liu J, Liu D, Yin F, Wei Z, et al. RACK1 promotes cancer progression by increasing the M2/M1 macrophage ratio via the NF-kappaB pathway in oral squamous cell carcinoma. Mol Oncol. 2020;14(4):795-807. Epub 20200220. https://doi.org/10.1002/1878-0261.12644 PMID: 31997535
  7. Wang YW, Cheng HL, Ding YR, Chou LH, Chow NH. EMP1, EMP 2, and EMP3 as novel therapeutic targets in human cancer. Biochim Biophys Acta Rev Cancer. 2017;1868(1):199-211. Epub 20170410. https://doi.org/10.1016/j.bbcan.2017.04.004 PMID: 28408326
  8. Ahmat Amin MKB, Shimizu A, Zankov DP, Sato A, Kurita S, Ito M, et al. Epithelial membrane protein 1 promotes tumor metastasis by enhancing cell migration via copine-III and Rac1. Oncogene. 2018;37(40):5416-34. Epub 20180604. https://doi.org/10.1038/s41388-018-0286-0 PMID: 29867202
  9. Wang M, Liu T, Hu X, Yin A, Liu J, Wang X. EMP1 promotes the malignant progression of osteosarcoma through the IRX2/MMP9 axis. Panminerva Med. 2020 Sep;62(3):150–4.
  10. Liu Y, Ding Y, Nie Y, Yang M. EMP1 Promotes the Proliferation and Invasion of Ovarian Cancer Cells Through Activating the MAPK Pathway. Onco Targets Ther. 2020;13:2047-55. Epub 20200309. https://doi.org/10.2147/OTT.S240028 PMID: 32210572
  11. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545-50. Epub 20050930. https://doi.org/10.1073/pnas.0506580102 PMID: 16199517
  12. Pastushenko I, Brisebarre A, Sifrim A, Fioramonti M, Revenco T, Boumahdi S, et al. Identification of the tumour transition states occurring during EMT. Nature. 2018;556(7702):463-8. Epub 20180418. https://doi.org/10.1038/s41586-018-0040-3 PMID: 29670281
  13. Bakir B, Chiarella AM, Pitarresi JR, Rustgi AK. EMT, MET, Plasticity, and Tumor Metastasis. Trends Cell Biol. 2020;30(10):764-76. Epub 20200813. https://doi.org/10.1016/j.tcb.2020.07.003 PMID: 32800658
  14. Paluch EK, Aspalter IM, Sixt M. Focal Adhesion-Independent Cell Migration. Annu Rev Cell Dev Biol. 2016;32:469-90. Epub 20160804. https://doi.org/10.1146/annurev-cellbio-111315-125341 PMID: 27501447
  15. Sulzmaier FJ, Jean C, Schlaepfer DD. FAK in cancer: mechanistic findings and clinical applications. Nat Rev Cancer. 2014;14(9):598-610. Epub 20140807. https://doi.org/10.1038/nrc3792 PMID: 25098269
  16. Czerniak B, Dinney C, McConkey D. Origins of Bladder Cancer. Annu Rev Pathol. 2016;11:149-74. Epub 20160222. https://doi.org/10.1146/annurev-pathol-012513-104703 PMID: 26907529
  17. Qin Y, Rodin S, Simonson OE, Hollande F. Laminins and cancer stem cells: Partners in crime? Semin Cancer Biol. 2017;45:3-12. Epub 20160801. https://doi.org/10.1016/j.semcancer.2016.07.004 PMID: 27491691
  18. Ngambenjawong C, Gustafson HH, Pun SH. Progress in tumor-associated macrophage (TAM)-targeted therapeutics. Adv Drug Deliv Rev. 2017;114:206-21. Epub 20170425. https://doi.org/10.1016/j.addr.2017.04.010 PMID: 28449873
  19. Liu C, He D, Zhang S, Chen H, Zhao J, Li X, et al. Homogeneous Polyporus Polysaccharide Inhibit Bladder Cancer by Resetting Tumor-Associated Macrophages Toward M1 Through NF-kappaB/NLRP3 Signaling. Front Immunol. 2022;13:839460. Epub 20220504. https://doi.org/10.3389/fimmu.2022.839460 PMID: 35603205
  20. Martinez VG, Rubio C, Martinez-Fernandez M, Segovia C, Lopez-Calderon F, Garin MI, et al. BMP4 Induces M2 Macrophage Polarization and Favors Tumor Progression in Bladder Cancer. Clin Cancer Res. 2017;23(23):7388-99. Epub 20170919. https://doi.org/10.1158/1078-0432.CCR-17-1004 PMID: 28928159
  21. Kobatake K, Ikeda KI, Nakata Y, Yamasaki N, Ueda T, Kanai A, et al. Kdm6a Deficiency Activates Inflammatory Pathways, Promotes M2 Macrophage Polarization, and Causes Bladder Cancer in Cooperation with p53 Dysfunction. Clin Cancer Res. 2020;26(8):2065-79. Epub 20200211. https://doi.org/10.1158/1078-0432.CCR-19-2230 PMID: 32047002
  22. Gilkes DM, Semenza GL, Wirtz D. Hypoxia and the extracellular matrix: drivers of tumour metastasis. Nat Rev Cancer. 2014;14(6):430-9. Epub 20140515. https://doi.org/10.1038/nrc3726 PMID: 24827502
  23. Yuzhalin AE, Lim SY, Kutikhin AG, Gordon-Weeks AN. Dynamic matrisome: ECM remodeling factors licensing cancer progression and metastasis. Biochim Biophys Acta Rev Cancer. 2018;1870(2):207-28. Epub 20181012. https://doi.org/10.1016/j.bbcan.2018.09.002 PMID: 30316942
  24. Kumra H, Reinhardt DP. Fibronectin-targeted drug delivery in cancer. Adv Drug Deliv Rev. 2016;97:101-10. Epub 20151127. https://doi.org/10.1016/j.addr.2015.11.014 PMID: 26639577
  25. Wu T, Zhang DL, Wang JM, Jiang JY, Du X, Zeng XY, et al. TRIM29 inhibits miR-873-5P biogenesis via CYTOR to upregulate fibronectin 1 and promotes invasion of papillary thyroid cancer cells. Cell Death Dis. 2020;11(9):813. Epub 20200929. https://doi.org/10.1038/s41419-020-03018-3 PMID: 32994394
  26. Ruiz-Garcia E, Scott V, Machavoine C, Bidart JM, Lacroix L, Delaloge S, et al. Gene expression profiling identifies Fibronectin 1 and CXCL9 as candidate biomarkers for breast cancer screening. Br J Cancer. 2010;102(3):462-8. Epub 20100112. https://doi.org/10.1038/sj.bjc.6605511 PMID: 20068563
  27. Gharbaran R. Advances in the molecular functions of syndecan-1 (SDC1/CD138) in the pathogenesis of malignancies. Crit Rev Oncol Hematol. 2015;94(1):1-17. Epub 20141218. https://doi.org/10.1016/j.critrevonc.2014.12.003 PMID: 25563413
  28. Ibrahim SA, Gadalla R, El-Ghonaimy EA, Samir O, Mohamed HT, Hassan H, et al. Syndecan-1 is a novel molecular marker for triple negative inflammatory breast cancer and modulates the cancer stem cell phenotype via the IL-6/STAT3, Notch and EGFR signaling pathways. Mol Cancer. 2017;16(1):57. Epub 20170307. https://doi.org/10.1186/s12943-017-0621-z PMID: 28270211
  29. Morikawa M, Derynck R, Miyazono K. TGF-beta and the TGF-beta Family: Context-Dependent Roles in Cell and Tissue Physiology. Cold Spring Harb Perspect Biol. 2016;8(5). Epub 20160502. https://doi.org/10.1101/cshperspect.a021873 PMID: 27141051
  30. Hegde PS, Chen DS. Top 10 Challenges in Cancer Immunotherapy. Immunity. 2020 Jan;52(1):17–35.
  31. Mantovani A, Marchesi F, Malesci A, Laghi L, Allavena P. Tumour-associated macrophages as treatment targets in oncology. Nat Rev Clin Oncol. 2017;14(7):399-416. Epub 20170124. https://doi.org/10.1038/nrclinonc.2016.217 PMID: 28117416
  32. Patel VG, Oh WK, Galsky MD. Treatment of muscle-invasive and advanced bladder cancer in 2020. CA Cancer J Clin. 2020;70(5):404-23. Epub 20200807. https://doi.org/10.3322/caac.21631 PMID: 32767764