Jian Qi1,2,3†
Yu Liu
Zhen Dou
Wulin Yang- 1Anhui Province Key Laboratory of Medical Physics and Technology, Institute of Health and Medical Technology, Hefei Institutes of Physical Science, Chinese Academy of Sciences, Hefei, China
- 2Scinece Island Branch, Graduate School of USTC, Hefei, China
- 3Hefei Cancer Hospital, Chinese Academy of Sciences, Hefei, China
- 4Department of Neurosurgery, The Shenzhen People’s Hospital (The Second Clinical Medical College of Jinan University), Shenzhen, China
- 5Department of Anatomy, Shanxi Medical University, Taiyuan, China
- 6Hefei National Science Center for Physical Sciences at the Microscale, University of Science and Technology of China, Hefei, China
Although research into immunotherapy is growing, its use in the treatment of breast cancer remains limited. Thus, identification and evaluation of prognostic biomarkers of tissue microenvironments will reveal new immune-based therapeutic strategies for breast cancer. Using an in silico bioinformatic approach, we investigated the tumor microenvironmental and genetic factors related to breast cancer. We calculated the Immune score, Stromal score, Estimate score, Tumor purity, TMB (Tumor mutation burden), and MATH (Mutant-allele tumor heterogeneity) of Breast cancer patients from the Cancer Genome Atlas (TCGA) using the ESTIMATE algorithm and Maftools. Significant correlations between Immune/Stromal scores with breast cancer subtypes and tumor stages were established. Importantly, we found that the Immune score, but not the Stromal score, was significantly related to the patient's prognosis. Weighted correlation network analysis (WGCNA) identified a pattern of gene function associated with Immune score, and that almost all of these genes (388 genes) are significantly upregulated in the higher Immune score group. Protein-protein interaction (PPI) network analysis revealed the enrichment of immune checkpoint genes, predicting a good prognosis for breast cancer. Among all the upregulated genes, FPR3, a G protein-coupled receptor essential for neutrophil activation, is the sole factor that predicts poor prognosis. Gene set enrichment analysis analysis showed FRP3 upregulation synergizes with the activation of many pathways involved in carcinogenesis. In summary, this study identified FPR3 as a key immune-related biomarker predicting a poor prognosis for breast cancer, revealing it as a promising intervention target for immunotherapy.
Introduction
Breast cancer, the most common gynecological cancer worldwide (Ghoncheh et al., 2016), is becoming a majorpublic health crisis, and the number of new cases diagnosed each year is ever-increasing (Cheng et al., 2018). Breast cancer is considered less immunogenic than melanoma or renal cell carcinoma, and the results of adoptive immunotherapy (interleukin-2, interferon) have been relatively disappointing (Burugu et al., 2017). Over the past decade, with the increased understanding of the immune microenvironment of breast cancer tissues, immune escape has been considered an important feature of breast cancer development (Romaniuk and Lyndin, 2015; Takada et al., 2018). Targeting the tumor immune microenvironment in breast cancer is of high therapeutic interest (Zhao et al., 2017). However, the therapeutic effects of immune checkpoint inhibition may be limited. For example, the evaluation of avelumab, an anti-PDL1 antibody, in various subtypes of breast cancer showed that the overall response rate (ORR) for the entire cohort was 4.8% (Emens, 2018), far from achieving the intended effect.
The tumor microenvironment is composed of a variety of immune cells and stromal cells, endothelial cells along with inflammatory mediators, and extracellular matrix (ECM) molecules (Cooper et al., 2012; Straussman et al., 2012). It plays a key role in altering the tumor response to treatment (Hanahan and Weinberg, 2000; Hanahan and Coussens, 2012). Previous studies have shown that high levels of immune cell infiltration are associated with better prognosis for diseases such as (Tas and Erturk, 2017) prostate cancer (Donovan et al., 2018), cutaneous melanoma (Yang et al., 2020), and breast cancer (Papatestas et al., 1976; Manuel et al., 2012). Also, high immune infiltration is associated with neoadjuvant chemotherapy and an enhanced response to adjuvant chemotherapy (Pruneri et al., 2018). Hence, assessing TME heterogeneity and reshaping the immune microenvironment may hold promise for cancer treatments in the near future. ESTIMATE (Estimate of Stromal and Immune Cells in Malignant Tumor Tissues from Expression Data) is a newly developed algorithm that applies gene expression data to predict the fraction of stromal and immune cells in tumor samples (Yoshihara et al., 2013). Weighted gene co-expression network analysis (WGCNA) is an effective tool to establish correlation patterns between genes to identify cancer-related modules and central genes (Langfelder and Horvath, 2008). Combining multiple bioinformatic approaches, we explored the microenvironment and genetic factors associated with breast cancer to determine the prognostic biomarkers for breast cancer.
Materials and Methods
Gene Expression Datasets
The level 3 gene expression profile (level 3 data) for breast cancer was obtained from the UCSC Xena (https://xenabrowser.net/datapages/). Clinical data, such as gender, race, age, histological type, survival, and outcome, were downloaded from the TCGA data portal. Count data were used to quantitate a total of 19,986 protein-coding genes that had been annotated in the Ensembl database (http://asia.ensembl.org/index.html). The ESTIMATE algorithm was applied to the normalized expression matrix to determine the Immune/Stromal scores for each breast cancer sample. Immune score and Stromal score were calculated by applying the ESTIMATE algorithm to the downloaded database (Yoshihara et al., 2013). To verify the association between FPR3 mRNA expression and survival, a dataset, GSE11121 (contains 200 samples), was obtained from the Gene Expression Omnibus (GEO) database for testing.
Correlations Between Prognosis and Immune/Stromal Score
The overall survival rate was taken as the main prognostic indicator. We divided the relevant patients into two groups based on the Immune/Stromal scores of each breast cancer sample. Kaplan-Meier plots were generated to illustrate the relationship between patients’ overall survival and gene expression levels of differentially expressed genes (DEGs). The relationship was tested by the log-rank test.
Identification of Differentially Expressed Genes (DEGs)
According to the ESTIMATE results, all samples were divided into high/low immune-score groups and high/low stromal-score groups to select the intersection genes. p Value <0.05, Fold change >2 were set as the cutoffs to identify significantly differentially expressed genes. Heatmaps were generated using the pheatmap package in R software (Galili et al., 2018) (cran.r-project.org/web/package/pheatmap/index.html), and the limma package (Ritchie et al., 2015) (limma package; www.r-project) was used to separate the upregulated and downregulated genes in the high score. Venn diagrams were drawn up with web tools (http://bioinformatics.psb.ugent.be/webtools/Venn/). The WGCNA method (Langfelder and Horvath, 2008) was used to construct the co-expression network of the genes in the test samples of breast cancer.
Gene Ontology and KEGG Pathway Enrichment Analysis
In our study, Gene Ontology (GO) analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were employed to understand the potential function of genes using the clusterProfiler package in R software (Yu et al., 2012). Moreover, to analyze the interaction genes, the protein-protein interaction (PPI) network was built using STRING (Szklarczyk et al., 2019) and reconstructed in Cytoscape v3.6, as previously stated (Shannon et al., 2003; Liu et al., 2020). Finally, Molecular Complex Detection (MCODE) in Cytoscape was utilized to obtain clusters based on the topology to localize densely connected regions.
Gene Expression Analysis in GEPIA
The online database Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/index.html) (Tang et al., 2017) is an interactive website that includes 9,736 tumors and 8,587 normal samples from the TCGA portal and the GTEx project. GEPIA is based on gene expression with the log-rank test and the MantelCox test in 33 different types of cancer. Gene expression correlation analysis was performed for a given set of TCGA expression data. The Spearman method was used to determine the correlation coefficient.
Gene Set Enrichment Analysis
GSEA is a calculation method used in determining whether a set of basically defined gene sets exhibit statistically significant differences between two biological states (Subramanian et al., 2005). In this study, GSEA first generated an ordered list of all genes based on their correlation with FPR3 expression. Gene set permutations were performed 1,000 times for each analysis. The expression level of FPR3 is deemed to be a phenotypic marker. The nominal p value and normalized enrichment score (NES) served to sort the pathways enriched in each phenotype.
Statistical Analysis
All data are expressed as mean ± SD (standard deviation). All analyses were performed with R version 3.5.3 (http://www.R-project.org). Immune and Stromal scores were calculated by using the ESTIMATE package. Data were analyzed with standard statistical tests where appropriate. Multiple testing was adjusted by using the FDR method.
Results
Immune and Stromal Score Correlate With Breast Cancer Subtypes
We found significant correlations between Immune score and Stromal score with breast cancer subtypes. In this study, we obtained information about 1,040 breast cancer patients from the TCGA database, including their gene expression signatures and clinical profiles. The average age of the patients was 58.34 years. Using the PAM50 classification method, we used the genefu package (Gendoo et al., 2016) in R language to divide the breast cancer subtypes into luminal A (LumA), luminal B (LumB), Her2+, and Basal. The number of patients with luminal A was 286, the number of patients with subtype luminal B was 427, the number of patients with Her2+ was 128, and the number of patients with Basal was 199. According to the calculation results of the ESTIMATE and maftools (Mayakonda et al., 2018) algorithm, as shown in Figures 1A–F, the Immune score, Stromal score, Estimate score, Tumor purity, TMB, and MATH are significantly correlated with different breast cancer subtypes. Among the Immune score, the basal subtype had the highest average Immune score, followed by the Her2+ subtype, LumA third, and LumB the lowest (p = 2.7e − 12); Similarly, in breast cancer subtypes, the order of Stromal score from high to low is LumA > Her2+ > LumB > Basal (p < 2.2e − 16). In Estimate score, the score from high to low is HER2+ > Basal > LumA > LumB (p = 4.2e − 10); in Tumor purity, the score from high to low is LumB > Basal > LumA > Her2+ (p = 4.2e − 10). In TMB, the score from high to low is Basal > Her2+ > LumB > LumA (p < 2.2e − 16), while in MATH, the score from high to low is Basal > LumB > Her2+ > LumA (p = 2.4e − 15).
FIGURE 1. Tumor microenvironmental factors were tightly associated with breast cancer subtypes and Immune score (A) Box-plot shows a significant association between breast cancer subtypes and the level of Immune score (n = 1,040, p = 2.7e − 12) (B) Box-plot reveals the significant association between breast cancer subtypes and the level of Stromal score (n = 1,040, p < 2.2e − 16) (C–F) Box-plot shows a significant correlation between breast cancer subtype and Estimate score (n = 1,040, p = 4.2e − 10), Tumor purity (n = 1,040, p = 4.2e − 10), TMB (n = 1,040, p < 2.2e − 16), MATH (n = 1,040, p = 2.4e − 15) (G) Breast cancer cases were divided into two groups based on their average expression of Immune scores. A high Immune score predicts a favorable prognosis for overall survival (OS). As indicated by the log-rank test, p = 0.0011 (H) Similarly, breast cancer cases were parted into two groups based on their average expression of Stromal scores. The survival time for OS of the high Stromal score group is longer than that of the low Stromal score group, although the difference was not statistically significant (p = 0.42) (I) A similar grouping based on Immune scores was used for RFS (relapse-free survival) survival analysis (J) A similar grouping based on Stromal scores was used for RFS survival analysis.
Next, to explore the correlations between Immune/Stromal scores and breast cancer prognosis, we constructed survival curves of patients by classifying them into high and low score groups based on their gene expression profiles. We found that patients with higher Immune scores have longer survival rates for OS (Overall survival) than those with lower Immune scores (p < 0.0011). Similarly, patients with high Stromal scores have better life expectancy than patients with low Stromal scores, although the data is not statistically significant (Figures 1G,H). In addition, we also calculated survival curves for RFS (relapse-free survival), DFS (Disease-free survival), and PFS (progression-free survival). We can see that the higher the Immune score, the better the survival rates for RFS, DFS, and PFS. The results are shown in Figures 1I,J (RFS) and supplementary Figure S1 (DFS and PFS). Although the p value is not significant, it is close to 0.05.
The Significant Correlation Between Immune/Stromal Scores and T Stage of Breast Cancer
The TNM staging system is a globally recognized standard for classifying the extent of the spread of cancer into nearby tissue. Given the correlations of both Immune- and Stromal scores with breast cancer subtypes, we studied their connection to TNM stages of breast cancer. We noted that the Immune score, Stromal score, and Estimate score decrease following the T stage progression (Figures 2A–C). In line with our observation in Figure 1G, the higher the Immune score, the better the prognosis. Consistently, the Tumor purity, TMB, and MATH increased with the progression of T stages, correlating well with the degree of malignancy of breast cancer (Figures 2D–F). The M and N stages of the tumor are only correlated with TMB, but not related to the Immune score, Stromal score, Estimate score, tumor purity, or MATH (Supplementary Figure S2).
FIGURE 2. Immune and Stromal scores were significantly associated with the T stages of breast cancer. As the T stage progressed, the Immune score, Stromal score, and Estimate score decreased, and the scores of Tumor purity, TMB, and MATH increased. Pathologic_T, pathological T staging of tumor, represented by T1-T4 in turn. Kruskal_Wallis, the kruskal wallis test was used to compare multi-independent samples.
Analysis of Differentially Expressed Genes
To determine the correlation between the overall gene expression profile and the Immune/Stromal scores, we evaluated the sequencing data of 1,040 breast cancer patients in the TCGA database. In this analysis, immune-related genes were divided into the high-score and the low-score group. We then used the limma package of R language to analyze the differentially expressed genes. Heatmaps in Figures 3A,B show distinct gene expression profiles of DEGs corresponding to the high vs. low Immune score/Stromal score groups. A total of 859 genes were upregulated and 40 genes were downregulated in the group of high Immune score (FC > 2, p value < 0.05). Similarly, there were 1,011 upregulated genes and five downregulated genes in the high Stromal score group (FC > 2, p value < 0.05) (Figures 3C,D).
FIGURE 3. Comparison of gene expression profiles in high vs. low Immune or Stromal scores (A) Heatmap of the DEGs in the Immune score of high vs. low. Genes upregulated are displayed in red; genes downregulated are shown in green. Genes with no change are in black (p < 0.05, FC > 2). FC, fold change (B) Heatmap of the DEGs in Stromal score of high vs. low (p < 0.05, FC > 2) (C–D) Volcano plot of gene expression profiles in high vs low Immune or Stromal scores. A total of 859 genes were upregulated and 40 genes were downregulated in the group of high Immune score (FC > 2, p < 0.05). Similarly, there were 1,011 upregulated genes and five downregulated genes in the high Stromal score group (FC > 2, p < 0.05) (E–F) Venn diagrams display the number of commonly upregulated (C) or downregulated (D) DEGs in the high Immune or Stromal score groups.
We found that 280 genes were upregulated in the groups with both high Immune score and high Stromal score (Figure 3E) and only one common gene was downregulated (Figure 3F). On the other hand, the Stromal score was not significantly related to the prognosis of breast cancer patients (Figures 1G,H). Therefore we decided to exclude the Stromal score from our analysis when verifying gene modules associated with Immune score.
Functional Enrichment Analysis Reveals FPR3 as a Key Immune-Related Gene
Using WGCNA, we sought to identify gene modules that are associated with the Immune score. Among the eight modules, the MEblue was highly correlated with the Immune score (R2 = 0.975), Stromal score (R2 = 0.538), Estimate score (R2 = 0.885), and TMB (R2 = 0.043). The MEgreen module showed a higher correlation with the Stromal score (R2 = 0.887, Figure 4A). Since the Immune score impacts the survival time of patients (Figure 1G), we identified 388 genes that were overlapped with WGCNA analysis and the high Immune score group (Figure 4B; Supplementary Table S1). To postulate the underlying mechanism for the upregulation of these identified genes, we first used the STRING tool to perform a functional analysis of the PPI network and applied the MCODE plugin in Cytoscape to obtain the two most significant modules. The first MCODE module contains 53 genes, including almost all the established immune checkpoint genes, for example, CD274, PDCD1, CTLA4, and LAG3. (Figure 4C). Intriguingly, many of the immune checkpoint genes predict a better prognosis for breast cancer (Supplementary Figure S3), providing a molecular explanation for the limited effects of targeting immune checkpoints in the treatment of breast cancer.
FIGURE 4. Network visualization plots by WGCNA (A) Module-trait associations. Each row corresponds to a gene module, and each column corresponds to a trait. Each cell was labeled by the corresponding correlation coefficient and p-value (B) Venn diagram shows the genes obtained by the intersection of genes in the WGCNA Meblue module and genes upregulated in high Immune score group (C–D) Cytoscape analysis of 388 intersecting genes yielded two most significantly related MCODE modules. Some known immune checkpoints and FPR3 are marked yellow (E) Multivariate Cox analysis was utilized to analyze the hazard ratio (HR) of FPR3 and some known immune checkpoints.
To identify more attractive prognostic targets in the immune microenvironment, we performed a survival analysis for these 388 genes and found that a total of 101 genes showed a significant difference in patient survival (Supplementary Figure S4). We found that only one particular gene, FPR3, predicts poor prognosis under higher expression, which was included in MCODE module 2 (Figure 4D). To evaluate whether FPR3 is an independent prognostic factor in the immune microenvironment, multivariate Cox analysis was used to analyze the hazard ratio (HR) of FPR3 and some known immune checkpoints (Figure 4E). The six known immune checkpoints do not influence the survival rate of patients (Figure 4E). Only FPR3 falls to the right of the invalid line (hazard ratio >1, p <0.001), which means that the higher the FPR3 expression, the worse the patient's prognosis. In addition, the GSE11121 dataset from the GEO database was used to verify the COX analysis. The results also show that FPR3 is an independent prognostic factor (hazard ratio >1, p <0.0298), which is more significant than other immune checkpoints (Supplementary Figure S5A).
Furthermore, we analyzed the functions of genes in module two using a gene enrichment approach. A total of 46 genes are involved in leukocyte differentiation, immune response-activating cell surface receptor signaling pathway, and positive regulation of cell activation (Figure 5A), with their molecular functions relating to protein tyrosine kinase activity, peptide binding, and amide binding (Figure 5B). These genes are predicted to localize to the plasma membrane (Figure 5C) and are required for Th1 and Th2 cell differentiation, Th17 cell differentiation, and Epstein-Barr virus infection, as indicated by the KEGG pathway analysis (Figure 5D). These results indicate that FPR3, as a component of module 2, may cooperate with other partners in immune-based biological pathways.
FIGURE 5. Functional enrichment analysis of FPR3-related gene module (A) Biological process (B) Molecular function (C) Cellular component, and (D) KEGG pathways are depicted. The color of the dots demonstrates -Log10(FDR). Red dots indicate smaller FDRs than blue dots. The size of the dots indicates the number of genes enriched in each analysis.
FPR3 Expression in Breast Cancer and the Pathways Enrichment Determined by GSEA
We used TIMER to evaluate the expression of FPR3 in pan-cancer. As displayed in Figure 6A, FPR3 is highly expressed in a variety of cancers including different breast cancer subtypes, COAD cancer, and HNSC cancer. As shown in both Oncomine and GEPIA databases, the expression of FPR3 in breast cancer is much higher than that in adjacent normal tissue (Figures 6B,C). The Kaplan-Meier survival curve (Data comes from TCGA) shows that the expression of FPR3 in breast cancer is negatively related to the survival of patients (Figure 6D). Similarly, we also used the GEO dataset to verify the relationship between FPR3 expression and survival. The result is consistent with TCGA data indicating an unfavorable prognosis (Supplementary Figure S5B). The above results suggest that inhibiting FPR3 expression may be useful as an intervention strategy.
FIGURE 6. FPR3 expression in breast cancer (A) TIMER showed the expression level of FPR3 in pan-cancer, and the expression of FPR3 in total breast cancer cases and different subtypes was higher than that in adjacent normal tissues (B) The transcription levels of FPR3 in different types of cancers (Oncomine database). The graph shows the number of data sets in which the mRNA expression of the target gene was significantly up-regulated (red) or down-regulated (blue). The threshold was defined by the following parameters: P < 1E-4 and FC > 2 (C) To further examine the expression level of FPR3 between breast cancer and normal tissues, FPR3 was examined using the GEPIA web-based tool (*, p < 0.05). Red color indicates tumor tissue and gray color indicates normal tissues (D) Kaplan-Meier survival curves were generated for FPR3. The higher the expression of FPR3 (yellow line), the shorter the survival period. p value in the log-rank test. OS, overall survival.
To reveal the functional role of FPR3, gene set enrichment analysis (GSEA) was used to analyze the gene expression matrix acquired from TCGA. The samples were divided into high and low expression groups according to the median expression level of FPR3. The top two enriched pathways in the high expression group of FPR3 were “pathways in cancer” (Figure 7A) and “cytokine-cytokine receptor interaction” (Figure 7B), implying that FPR3 upregulation may lead to alterations in cancer-related pathways and cytokine-based immune regulations. We also performed a correlation analysis of FPR3 expression on the GSEA-enriched pathways. In the “pathways in cancer”, PIK3R5, SPI1, and CSF1R are most relevant to FPR3 expression (Figures 7A1–A3). In the “cytokine-cytokine receptor interaction pathway”, the most relevant genes to FPR3 expression are CCR1, IL10, and IL10RA (Figures.7B1–B3). FPR3 may directly or indirectly cooperate with these genes in promoting tumorigenesis. For instance, the high correlation of FPR3 with PIK3R5 implies that FPR3 may promote tumorigenesis through the PIK3R5-mediated G-protein coupled receptor activation (Barberis and Hirsch, 2008). As annotated in the “pathways in cancer”, high FPR3 expression associates with many pathways that promote tumorigenesis and development, such as the VEGF signaling pathway, MAPK signaling pathway, and PI3K-AKT signaling pathway. Because FPRs belong to the classic chemotactic GPCR subfamily, we propose that FPR3 is involved in the G protein receptor coupled pathway to promote carcinogenesis, as shown in the schematic diagram (Figure 7C). Through unknown ligands, FPR3 may regulate breast tumorigenesis through G-protein coupled PI3K or MAPK signaling cascades, participating in a series of carcinogenic processes, such as enhanced proliferation, sustained angiogenesis, and apoptosis evasion.
FIGURE 7. Enrichment plots by gene set enrichment analysis (GSEA). GSEA results showed that higher FPR3 expression was associated with activation of “pathways in cancer” (A) (p = 0.0011) and “cytokine-cytokine-receptor-interaction” (B) (p = 0.0011) (A1-A3) shows the co-expression correlation of FPR3 with top genes in the “pathways in cancer”. (B1-B3) shows the co-expression correlation of FPR3 with top genes in the “cytokine-cytokine-receptor-interaction” pathway (C) A model diagram was proposed that FPR3 may be involved in the occurrence and progression of cancer via G protein-coupled signaling.
Discussion
The tissue microenvironment influences the extent of tumor initiation and development (Quail and Joyce, 2013; Denkert et al., 2018). Four clinically relevant molecular subtypes of breast cancer, luminal A, luminal B, Her2+ type, and basal type, are different in terms of their morbidity, survival rate, prognosis, and tumor biological characteristics. Such patient stratification provides clinical and economic value in breast cancer treatment (Blok et al., 2018). The tumor microenvironment is where the immune system interacts with the tumor (Merlano et al., 2019) and any changes in the composition of the tumor microenvironment may affect the fate of malignant tumors. Therefore, important components of the tumor microenvironment, including stromal- and immune cells, play key roles in cancer progression (Xiong et al., 2018). In the current study, we attempted to identify genes in the tumor microenvironment that are relevant to the survival of breast cancer patients from the TCGA database. Several reports have demonstrated the successful application of estimation algorithms in exploring the gene expression signatures associated with cancers (Jia et al., 2018; Chen et al., 2019; Li and Xu, 2019). We used this algorithm to obtain the Immune score, Stromal score, Estimate score, and tumor purity for breast cancer. We found that these scores significantly correlate with breast cancer subtypes. We used Maftools to establish the correlation between TMB/MATH and breast cancer subtypes. TNM reflects an important aspect of the clinical characteristics of breast cancer. We showed that T stage progression is inversely related to the Immune score, suggesting an important role for microenvironmental factors in tumor progression.
Importantly, we found that overall survival was positively correlated with the Immune score, but not significantly with the Stromal score. This finding suggests that immune factors in the tumor microenvironment play a more important role in determining patient outcomes. WGCNA analysis identified the MEblue module that most relates to the Immune score. Key genes in the MEblue module, including most of all known immune checkpoints, were almost totally upregulated in the high Immune score group. Surprisingly, survival analysis showed that all of these immune checkpoints either had no effect on patient survival or were protective. The FPR3 gene is the only gene that is adverse to the survival of patients. We further found that FPR3 is an independent hazard factor in the immune microenvironment of breast cancer, and its expression is accompanied by the activation of many pathways in cancer, highlighting the essential role of FPR3 in cancer progression.
Immune checkpoint antagonists have altered the way cancer is treated and produced a more durable response, and the FDA has approved a number of them to treat many cancers, but not breast cancer (Lipson et al., 2015). The majority of breast cancer patients respond poorly to checkpoint blockade (Emens, 2018). Moreover, in the clinical evaluation of the safety of PD1 inhibitors, some patients suffer myalgia, fatigue, joint pain, and nausea (Dirix et al., 2018). Although the FDA approved atezolizumab in combination with nab-paclitaxel as a first-line treatment for TNBC in 2019 (Schmid et al., 2020), the efficacy of a single drug was quite poor. Our results also suggest that the expression levels of the known immune checkpoint molecules do not necessarily affect overall patient survival. Therefore, it is of great importance to find novel prognostic molecular markers in the immune microenvironment of breast cancer.
Formyl-peptide receptors (FPRs) belong to the classical GPCR subfamily, and three FPRs have been identified in humans: FPR1-FPR3. Activation of FPR1 and FPR2 by chemotactic agonists elicits a cascade of signaling events that leads to myeloid cell migration, mediator release, increased phagocytosis, and new gene transcription (Cattaneo et al., 2013). FPR1 and FPR2 have been reported to be abnormally expressed in various tumors (Prevete et al., 2015). The expression of FPR1 in gastric cancer tissue is higher than that in normal tissue and is closely related to the survival time of the patient (Yang et al., 2011), and FPR2 is also highly expressed in endometrial cancer and colon cancer (Cocco et al., 2010). However, few studies have been performed on FPR3. FPR3 is expressed in eosinophils, monocytes, macrophages, and dendritic cells, but its function is unclear (Dorward et al., 2015; Nawaz et al., 2020). Several ligands for FPR3 have been identified, including F2L, an acetylated N-terminal fragment of human heme-binding protein (Migeotte et al., 2005), and the neuroprotective peptide humanin (Harada et al., 2004). Interestingly, FPR3 does not interact with formylated chemoattract peptides or ligands for FPR1 or FPR2. Therefore, FPR3 may have a unique functional role (Rabiet et al., 2011). We proposed that FPR3 may regulate the behaviors of tumor cells or immune cells through unspecified ligands in the breast tumor microenvironment. Tumor microenvironment refers to an acidic extracellular fluid system consisting of tumor cells and non-malignant stromal cells (including fibroblast and immune cells) as well as an extracellular matrix, containing a large number of growth factors, peptides, and enzymes. The tumor microenvironment can provide ligands to activate receptors expressed on tumor cells or immune cells. Ligand-receptor binding on tumor cells might activate cancer-related pathways to promote tumorigenesis and development. Thus, activation of FPR3 may regulate various malignant behaviors of cancer cells, including proliferation, avoidance of apoptosis, and sustained angiogenesis, via its G-protein-coupled signaling cascades. On the other hand, ligand-receptor binding on immune cells may elicit signaling cascades to cause new gene transcription, mediator release, or cell migration. Therefore, it is necessary to understand the specific properties of the unknown ligand for FPR3, and to make an in-depth analysis of the signaling pathways activated in tumor cells and immune cells. Current work demonstrated that FPR3 is a prognostic marker in breast cancer progression. The detailed mechanisms of FPR3 in breast cancer still require further investigation.
Overall, this study identifies FPR3 as a key immune-related intervention target, improving our understanding of the complex roles of the immune microenvironment in breast cancer progression.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
WY and HW conceived of and supervised the study. WY, HW, JQ, and YL designed the experimental methods and analysis protocol. JQ, YL, JH, ZD, and LL conducted data collation and analysis. WY, HW, JQ, YL, ZD, LL, and HD discussed, explained the results and wrote the paper. All authors approved the submission of this manuscript.
Funding
This study was supported by the National Natural Science Foundation of China (No. 81872276), the Research Project Supported by Shanxi Scholarship Council of China (No. 2020-085), the Innovative Program of Development Foundation of Hefei Center for Physical Science and Technology (No. 2018CXFX006), and the foundation of Anhui Province Key Laboratory of Medical Physics and Technology (No. LMPT201908).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest
Acknowledgments
We thank Lim Chun-Yan from the University of California, Berkeley for his constructive comments and text editing of the manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2020.593247/full#supplementary-material.
References
Barberis, L., and Hirsch, E. (2008). Targeting phosphoinositide 3-kinase gamma to fight inflammation and more. Thromb. Haemost. 99 (2), 279–285. doi:10.1160/TH07-10-0632
Blok, E. J., Bastiaannet, E., van den Hout, W. B., Liefers, G. J., Smit, V. T. H. B. M., Kroep, J. R., et al. (2018). Systematic review of the clinical and economic value of gene expression profiles for invasive early breast cancer available in Europe. Cancer Treat Rev. 62, 74–90. doi:10.1016/j.ctrv.2017.10.012
Burugu, S., Asleh-Aburaya, K., and Nielsen, T. O. (2017). Immune infiltrates in the breast cancer microenvironment: detection, characterization and clinical implication. Breast Cancer 24 (1), 3–15. doi:10.1007/s12282-016-0698-z
Cattaneo, F., Parisi, M., and Ammendola, R. (2013). Distinct signaling cascades elicited by different formyl peptide receptor 2 (FPR2) agonists. Int. J. Mol. Sci. 14 (4), 7193–7230. doi:10.3390/ijms14047193
Chen, B., Chen, W., Jin, J., Wang, X., Cao, Y., and He, Y. (2019). Data mining of prognostic microenvironment-related genes in clear cell renal cell carcinoma: a study with TCGA database. Dis. Markers 2019, 8901649. doi:10.1155/2019/8901649
Cheng, Y., Yan, Y., Gong, J., Yang, N., and Nie, S. (2018). Trends in incidence and mortality of female breast cancer during transition in Central China. Cancer. Manag. Res. 10, 6247–6255. doi:10.2147/CMAR.S182510
Cocco, E., Bellone, S., El-Sahwi, K., Cargnelutti, M., Buza, N., Tavassoli, F. A., et al. (2010). Serum amyloid A: a novel biomarker for endometrial cancer. Cancer 116 (4), 843–851. doi:10.1002/cncr.24838
Cooper, L. A., Gutman, D. A., Chisolm, C., Appin, C., Kong, J., Rong, Y., et al. (2012). The tumor microenvironment strongly impacts master transcriptional regulators and gene expression class of glioblastoma. Am. J. Pathol. 180 (5), 2108–2119. doi:10.1016/j.ajpath.2012.01.040
Denkert, C., von Minckwitz, G., Darb-Esfahani, S., Lederer, B., Heppner, B. I., Weber, K. E., et al. (2018). Tumour-infiltrating lymphocytes and prognosis in different subtypes of breast cancer: a pooled analysis of 3771 patients treated with neoadjuvant therapy. Lancet Oncol. 19 (1), 40–50. doi:10.1016/S1470-2045(17)30904-X
Dirix, L. Y., Takacs, I., Jerusalem, G., Nikolinakos, P., Arkenau, H. T., Forero-Torres, A., et al. (2018). Avelumab, an anti-PD-L1 antibody, in patients with locally advanced or metastatic breast cancer: a phase 1b JAVELIN Solid Tumor study. Breast Cancer. Res. Treat. 167 (3), 671–686. doi:10.1007/s10549-017-4537-5
Donovan, M. J., Fernandez, G., Scott, R., Khan, F. M., Zeineh, J., Koll, G., et al. (2018). Development and validation of a novel automated Gleason grade and molecular profile that define a highly predictive prostate cancer progression algorithm-based test. Prostate Cancer Prostatic Dis. 21 (4), 594–603. doi:10.1038/s41391-018-0067-4
Dorward, D. A., Lucas, C. D., Chapman, G. B., Haslett, C., Dhaliwal, K., and Rossi, A. G. (2015). The role of formylated peptides and formyl peptide receptor 1 in governing neutrophil function during acute inflammation. Am. J. Pathol. 185 (5), 1172–1184. doi:10.1016/j.ajpath.2015.01.020
Emens, L. A. (2018). Breast cancer immunotherapy: facts and hopes. Clin. Cancer Res. 24 (3), 511–520. doi:10.1158/1078-0432.CCR-16-3001
Galili, T., O'Callaghan, A., Sidi, J., and Sievert, C. (2018). heatmaply: an R package for creating interactive cluster heatmaps for online publishing. Bioinformatics 34 (9), 1600–1602. doi:10.1093/bioinformatics/btx657
Gendoo, D. M., Ratanasirigulchai, N., Schröder, M. S., Paré, L., Parker, J. S., Prat, A., et al. (2016). Genefu: an R/Bioconductor package for computation of gene expression-based signatures in breast cancer. Bioinformatics 32 (7), 1097–1099. doi:10.1093/bioinformatics/btv693
Ghoncheh, M., Pournamdar, Z., and Salehiniya, H. (2016). Incidence and mortality and epidemiology of breast cancer in the world. Asian Pac. J. Cancer Prev. APJCP 17 (S3), 43–46. doi:10.7314/apjcp.2016.17.s3.43
Hanahan, D., and Coussens, L. M. (2012). Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell 21 (3), 309–322. doi:10.1016/j.ccr.2012.02.022
Hanahan, D., and Weinberg, R. A. (2000). The hallmarks of cancer. Cell 100 (1), 57–70. doi:10.1016/S0092-8674(00)81683-9
Harada, M., Habata, Y., Hosoya, M., Nishi, K., Fujii, R., Kobayashi, M., et al. (2004). N-Formylated humanin activates both formyl peptide receptor-like 1 and 2. Biochem. Biophys. Res. Commun. 324 (1), 255–261. doi:10.1016/j.bbrc.2004.09.046
Jia, D., Li, S., Li, D., Xue, H., Yang, D., and Liu, Y. (2018). Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging 10 (4), 592–605. doi:10.18632/aging.101415
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559. doi:10.1186/1471-2105-9-559
Li, S., and Xu, W. (2019). Mining TCGA database for screening and identification of hub genes in kidney renal clear cell carcinoma microenvironment. J. Cell. Biochem., 121, 3952. doi:10.1002/jcb.29511
Lipson, E. J., Forde, P. M., Hammers, H. J., Emens, L. A., Taube, J. M., and Topalian, S. L. (2015). Antagonists of PD-1 and PD-L1 in cancer treatment. Semin. Oncol. 42 (4), 587–600. doi:10.1053/j.seminoncol.2015.05.013
Liu, Y., Qi, J., Dou, Z., Hu, J., Lu, L., Dai, H., et al. (2020). Systematic expression analysis of WEE family kinases reveals the importance of PKMYT1 in breast carcinogenesis. Cell Prolif. 53 (2), e12741. doi:10.1111/cpr.12741
Manuel, M., Tredan, O., Bachelot, T., Clapisson, G., Courtier, A., Parmentier, G., et al. (2012). Lymphopenia combined with low TCR diversity (divpenia) predicts poor overall survival in metastatic breast cancer patients. OncoImmunology 1 (4), 432–440. doi:10.4161/onci.19545
Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118
Merlano, M. C., Abbona, A., Denaro, N., and Garrone, O. (2019). Knowing the tumour microenvironment to optimise immunotherapy. Acta Otorhinolaryngol Ital. 39 (1), 2–8. doi:10.14639/0392-100X-2481
Migeotte, I., Riboldi, E., Franssen, J. D., Grégoire, F., Loison, C., Wittamer, V., et al. (2005). Identification and characterization of an endogenous chemotactic ligand specific for FPRL2. J. Exp. Med. 201 (1), 83–93. doi:10.1084/jem.20041277
Nawaz, M. I., Rezzola, S., Tobia, C., Coltrini, D., Belleri, M., Mitola, S., et al. (2020). D-Peptide analogues of Boc-Phe-Leu-Phe-Leu-Phe-COOH induce neovascularization via endothelial N-formyl peptide receptor 3. Angiogenesis 23 (3), 357–369. doi:10.1007/s10456-020-09714-0
Papatestas, A. E., Lesnick, G. J., Genkins, G., and Aufses, A. H. (1976). The prognostic significance of peripheral lymphocyte counts in patients with breast carcinoma. Cancer 37 (1), 164–168. doi:10.1002/1097-0142(197601)37:1<164::aid-cncr2820370123>3.0.co;2-h
Prevete, N., Liotti, F., Visciano, C., Marone, G., Melillo, R. M., and de Paulis, A. (2015). The formyl peptide receptor 1 exerts a tumor suppressor function in human gastric cancer by inhibiting angiogenesis. Oncogene 34 (29), 3826–3838. doi:10.1038/onc.2014.309
Pruneri, G., Vingiani, A., and Denkert, C. (2018). Tumor infiltrating lymphocytes in early breast cancer. Breast 37, 207–214. doi:10.1016/j.breast.2017.03.010
Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19 (11), 1423–1437. doi:10.1038/nm.3394
Rabiet, M. J., Macari, L., Dahlgren, C., and Boulay, F. (2011). N-formyl peptide receptor 3 (FPR3) departs from the homologous FPR2/ALX receptor with regard to the major processes governing chemoattractant receptor regulation, expression at the cell surface, and phosphorylation. J. Biol. Chem. 286 (30), 26718–26731. doi:10.1074/jbc.M111.244590
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Romaniuk, A., and Lуndіn, M. (2015). Immune microenvironment as a factor of breast cancer progression. Diagn. Pathol. 10, 79. doi:10.1186/s13000-015-0316-y
Schmid, P., Rugo, H. S., Adams, S., Schneeweiss, A., Barrios, C. H., Iwata, H., et al. (2020). Atezolizumab plus nab-paclitaxel as first-line treatment for unresectable, locally advanced or metastatic triple-negative breast cancer (IMpassion130): updated efficacy results from a randomised, double-blind, placebo-controlled, phase 3 trial. Lancet Oncol. 21 (1), 44–59. doi:10.1016/S1470-2045(19)30689-8
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303
Straussman, R., Morikawa, T., Shee, K., Barzily-Rokni, M., Qian, Z. R., Du, J., et al. (2012). Tumour micro-environment elicits innate resistance to RAF inhibitors through HGF secretion. Nature 487 (7408), 500–504. doi:10.1038/nature11183
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47 (D1), D607–D613. doi:10.1093/nar/gky1131
Takada, K., Kashiwagi, S., Goto, W., Asano, Y., Takahashi, K., Hatano, T., et al. (2018). Significance of re-biopsy for recurrent breast cancer in the immune tumour microenvironment. Br. J. Cancer 119 (5), 572–579. doi:10.1038/s41416-018-0197-4
Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45 (W1), W98–W102. doi:10.1093/nar/gkx247
Tas, F., and Erturk, K. (2017). Tumor infiltrating lymphocytes (TILs) may be only an independent predictor of nodal involvement but not for recurrence and survival in cutaneous melanoma patients. Cancer Invest. 35 (8), 501–505. doi:10.1080/07357907.2017.1351984
Xiong, Y., Wang, K., Zhou, H., Peng, L., You, W., and Fu, Z. (2018). Profiles of immune infiltration in colorectal cancer and their clinical significant: a gene expression-based study. Cancer Med. 7 (9), 4496–4508. doi:10.1002/cam4.1745
Yang, S., Liu, T., Nan, H., Wang, Y., Chen, H., Zhang, X., et al. (2020). Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma. J. Cell. Physiol. 235 (2), 1025–1035. doi:10.1002/jcp.29018
Yang, Y., Liu, Y., Yao, X., Ping, Y., Jiang, T., Liu, Q., et al. (2011). Annexin 1 released by necrotic human glioblastoma cells stimulates tumor cell growth through the formyl peptide receptor 1. Am. J. Pathol. 179 (3), 1504–1512. doi:10.1016/j.ajpath.2011.05.059
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16 (5), 284–287. doi:10.1089/omi.2011.0118
Keywords: breast cancer, immune microenvironment, ESTIMATE algorithm, immune checkpoint, FPR3
Citation: Qi J, Liu Y, Hu J, Lu L, Dou Z, Dai H, Wang H and Yang W (2021) Identification of FPR3 as a Unique Biomarker for Targeted Therapy in the Immune Microenvironment of Breast Cancer. Front. Pharmacol. 11:593247. doi: 10.3389/fphar.2020.593247
Received: 18 September 2020; Accepted: 30 December 2020;
Published: 11 February 2021.
Edited by:
Feng-Yao Tang, China Medical University, TaiwanReviewed by:
Shao-Chih Chiu, China Medical University Hospital, TaiwanCharles C. N. Wang, Asia University, Taiwan
Copyright © 2021 Qi, Liu, Hu, Lu, Dou, Dai, Wang and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Wulin Yang, yangw@cmpt.ac.cn; Hongzhi Wang, wanghz@hfcas.ac.cn
†These authors have contributed equally to this work