Identification of potential genes correlated with breast cancer metastasis and prognosis

Abstract

We aimed to screen, via bioinformatics analyses, potential gene markers related to BC metastasis and prognosis. A gene expression dataset (GSE102818) and relevant clinical information were downloaded from the GEO database and TCGA, respectively. Then, differentially expressed genes (DEGs) were identified between primary and metastatic BC samples, followed by survival analyses. Subsequently, functional enrichment and protein–protein interaction (PPI) analyses were performed. Finally, transcriptional regulation analysis was performed to extract the pivotal genes. A total of 146 DEGs were screened. Of these, 56 downregulated DEGs were closely associated with the survival outcomes of BC. Furthermore, functional enrichment analysis showed that these survival-related genes, including CCR7, CD27, and IL2RB, played notable roles in pathways such as cytokine-cytokine receptor interaction and immune responses. Furthermore, PPI network and enriched sub-module network analyses revealed that CCR7, CD27, and IL2RB served as hub genes with a high degree of association. Accordingly, they were key players in BC-associated transcriptional regulatory networks. Notably, these genes exhibited strong relationships with miR-1908 and miR-663. CCR7, CD27, IL2RB, miR-1908, and miR-663 may participate in the molecular mechanisms underlying BC metastasis, probably by influencing cytokine-related pathways and immune responses; thus, they are potential, novel therapeutic targets for BC management.

KEYWORDS:

• Breast cancer
• survival analysis
• differential expressed genes
• functional analysis
• transcriptional regulatory analysis

Introduction

Breast cancer (BC), one of the most common malignant diseases, is a leading cause of mortality among women globally, especially, in developing nations (Bonilla et al. 2017). Investigators report that the incidence of BC will continuously increase and reach a value of approximately 85 in 100,000 women by 2021 (Akram et al. 2017). Notably, a large number of BC patients suffer from metastasis and recurrence (Siegel et al. 2017). Fortunately, advanced diagnostic techniques, sophisticated surgical operations, and improved pharmacological treatments, have considerably enhanced the survival rate of BC patients over the past few years (Peart 2015; Castaneda and Strasser 2017). Recently, multiple studies have focused on clarifying the molecular mechanisms underlying the etiology of diseases using a wide range of analytical approaches, including conventional experimental techniques, association analysis, microarrays, and high-throughput sequencing (Selvaraj et al. 2021) (Nangraj et al. 2020) (Huang et al. 2020). There have also been many reports on the molecular mechanism of BC (Marotti et al. 2017). For example, Segovia-Mendoza et al identified numerous components in immune cells that are involved in the breast tumor microenvironment and contribute to the development of BC (Segovia-Mendoza and Morales-Montor 2019). Moreover, Fatima et al identified a 32-bp deletion in chemokine CC receptor 5 that is closely correlated with the risk of BC and claimed that this signal may be responsible for the activation and invasion of immune cells at the site of tumorigenesis (Fatima et al. 2019). Qi et al suggested that some vital genes (CD27 and CR2) and pathways, such as immune responses, are possibly involved in the pathogenesis of triple-negative BC, according to microarray analysis data (Qi et al. 2019). Yang et al also identified several potential genetic markers (CCNB2, BUB1, and NDC80) associated with the occurrence and progression of basal-like BC using bioinformatics analyses based on two microarrays (Chen et al. 2019).

Additionally, the invasion and metastasis of BC has been reported to appear in a very short time-frame after the occurrence of the primary tumor (Weigelt et al. 2005). Using cDNA microarray analysis, Li et al previously found that AF1Q may be a potential genetic marker for predicting BC metastasis. Zheng et al identified several differentially expressed genes (DEGs) between normal and invasive ductal carcinoma breast samples and found that multiple genes, such as CAPG, TP53INP1, and SPI1 are possibly associated with BC metastasis (Zheng et al. 2017). Although existing evidence has explored the potential molecular markers for the development and metastasis of BC, there is a need for a deeper understanding of the mechanisms underlying, and progression of, BC.

In this study, the dataset GSE102818 from the Gene Expression Omnibus (GEO) was downloaded, including the gene expression data from primary tumors and metastatic lesions of BC cases. Next, the DEGs between primary and metastatic patients were identified. Meanwhile, the relevant clinical data from The Cancer Genome Atlas (TCGA) were also obtained to screen for survival-associated genes. Moreover, functional annotation, protein interaction, and transcriptional regulatory analyses of these survival-related genes were also performed. This work will contribute towards the screening of promising novel therapeutic targets for controlling the development and metastasis of BC and facilitate the development of new strategies to improve the prognosis of BC patients.

Materials and methods

Data acquisition

Gene expression profile data (from the GSE102818 dataset) were retrieved and downloaded from the National Center for Biotechnology Information GEO (Barrett et al. 2007) repository (http://www.ncbi.nlm.nih.gov/geo/). The GSE102818 dataset contains 31 primary BC samples and 17 corresponding metastatic BC samples. These 48 samples corresponded to 36 patients (10 ER-negative and 26 ER-positive), some of which came from paired primary and metastatic tissues of 10 patients (for one patient, tissues from two different metastases were obtained, and for another, two different biopsies from one primary tumor was performed). All samples were subjected to analysis using the GPL23249 Nanostring nCounter Gene Expression Assay (Pan Cancer Immune Panel). In addition, the RNA sequencing data of BC specimens and the relevant clinical data were obtained from the TCGA database, which included 1,212 BC samples. After sample examination, survival data from 1,100 BC samples were successfully obtained.

Data pre-processing and identification of dEgs

The raw microarray data was pre-processed using the scale method in R software. The affy package in R was used to normalize the data, which included format conversion, supplementing missing values, analyzing the background correction, and log2 transformation. The Limma package and t-test in R were used to extract the DEGs the between primary and metastatic BC samples. Moreover, all p values were computed and subsequently adjusted (adj.p.value) using the Benjamini & Hochberg method. Of note, the DEGs in this research were defined based on the following thresholds: adj.p.value < 0.05 and | log fold change (FC) | > 1.

Survival analysis of DEGs

TCGA has effectively enhanced the prevention and diagnosis of a wide variety of cancers, owing to its abundance of genetic and clinical data. Here, the clinical information of about 1,100 BC patients, including overall survival (OS) and OS status, was collected, and the screened DEGs were further matched with relevant clinicopathologic data. Briefly, the expression values of these DEGs were calculated, and accordingly, the specimens were divided into the high-expression and low-expression groups. Subsequently, a Log-rank test was implemented (p value < 0.05) to identify the DEGs related to the survival of BC patients. Finally, the survival curve was obtained using the Kaplan–Meier method.

Gene ontology (GO) and pathway enrichment analyses

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang et al. 2009), a free web-accessible bioinformatics resource, provided integrated functional annotation tool suites. Therefore, for candidate DEGs associated with BC survival, GO (Ashburner et al. 2000) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto 2000) pathway enrichment analyses were performed using DAVID. A p value < 0.05 and count > 2 were determined as the cutoffs for the functional analyses.

Protein–protein interaction (PPI) network construction and sub-module analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Szklarczyk et al. 2015) database is publicly available and has a powerful PPI-predictive ability. To further explore the interactions among the identified DEGs correlated with BC survival, we performed PPI analysis using a PPI score of 0.4, based on the STRING database. Next, PPI networks were constructed and visualized using Cytoscape (Shannon et al. 2003) software (v 3.6.1, http://www.cytoscape.org/). Furthermore, the MCODE (Bandettini et al. 2012) plugin in the CytoNCA software (v 1.4.2, http://apps.cytoscape.org/apps/MCODE) was utilized to analyze the topological characteristics of significantly enriched modules using the following criterion: PPI score parameter > 10.

miRNA-TF-target regulatory analysis

To gain a deeper understanding of the potential functions of the candidate DEGs, the ENCODE chip-seq dataset in the Enrichr (Chen et al. 2013) (http://amp.pharm.mssm.edu/Ehrichr/) database was used to predict the transcription factors (TFs) regulated by DEGs in the PPI module with the cutoff of p value < 0.05. Meanwhile, the miRTarBase dataset in the Enrichr database was used to predict miRNA-DEG pairs with a p value < 0.05 (used as the screening threshold). Finally, miRNA-TF-target regulatory networks were constructed using Cytoscape by integrating the TF-DEG interactive pairs with the miRNA-DEG interactions.

Results

Identification of dEgs and survival analysis

In total, 146 DEGs, including three significantly upregulated genes and 143 markedly downregulated genes, were extracted using the methods outlined above (Figure 1). Subsequently, these DEGs were subjected to survival analysis by matching them with the 1,100 BC samples from TCGA. Our results showed that 145 DEGs were further extracted via integration with relevant clinical information (Table S1). Moreover, a total of 56 downregulated DEGs were found to be correlated with BC survival characteristics, while three upregulated DEGs were not associated with BC prognosis, as shown in the Table S2.

Figure 1. Volcano plot of differentially expressed genes (DEGs). The red dots represent the upregulated DEGs and the blue dots indicate the downregulated DEGs.

Functional enrichment analysis of DEGs

The GO enrichment analysis indicated that downregulated DEGs, such as CD27, were markedly enriched in 75 GO terms correlated with biological processes such as immune response, 19 GO terms involved with molecular functions including transmembrane signaling receptor activity, and 14 GO terms related to cellular components that mainly focused on the plasma membrane (Figure 2). Additionally, the KEGG pathway analysis demonstrated that these downregulated DEGs, including CD27, C–C-chemokine receptor (CCR7), and interleukin 2 receptor- beta (IL2RB), played essential roles in cytokine-cytokine receptor interactions and the T cell receptor signaling pathway (Figure 2).

Figure 2. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of survival-related genes. The black line denotes the p value and the length of the bars shows the number of genes.

PPI network analysis

PPI network analysis was performed using the STRING database to explore the potential PPI relationships of DEGs with BC survival. A total of 56 nodes and 518 protein pairs were obtained (Figure 3A). The top 15 gene nodes are shown in Table 1; these were considered hub genes with a high degree of association. Subsequently, we extracted a significantly enriched sub-module using MCODE with a score ≥ 10; this sub-module included 22 nodes and 206 interactive pairs (Figure 3B).

Figure 3. The protein-protein interaction (PPI) network and sub-module network of survival-related genes. (A) PPI network; (B) significantly enriched sub-module network. The green circular nodes represent the downregulated genes.

Table 1. Top 15 hub genes in protein-protein interaction and sub-module network.

PPI network		Sub-module network
Gene symbol	Degree	Gene symbol	Degree
CD79B	22	CD27	21
CXCL9	22	TBX21	21
CCL5	22	CD2	21
CD27	21	CD3E	21
TBX21	21	SELL	20
CD1C	21	PRF1	20
CD2	21	CCR7	20
CD3E	21	ZAP70	20
KLRB1	20	CD247	20
SELL	20	IL2RB	20
PRF1	20	KLRK1	19
CCR7	20	CD5	19
ZAP70	20	LCK	19
CD247	20	CD8A	19
IL2RB	20	GZMA	18

miRNA- TF-target regulatory network analysis

According to the prediction analysis based on the Enrichr database, 31 miRNA-target pairs and four TF-target pairs were predicted and subsequently used to construct the TF-miRNA-target regulatory network. Moreover, the hub genes, including CCR7, CD27, and IL2RB, showed close interactions with miR-663 and miR-1908 (Table 2). Our findings also suggested that 35 miRNA-TF-target regulative pairs were achieved, containing 14 significant DEGs, six miRNAs (miR-3677, miR-663, miR-1908, miR-1296, miR-3193, and miR-492), and one TF (BCL11A_GM12878_hg19) (Figure 4).

Figure 4. Transcriptional regulatory network. The green circles represent the downregulated genes; the purple triangles represent the miRNAs; and the yellow hexagons represent the transcription factors.

Table 2. The interaction pairs among miRNAs, transcription factors and genes.

Regulators	Targets	Regulators	Targets
hsa-miR-1908	CD40LG	hsa-miR-1296	ITK
hsa-miR-1908	CD8A	hsa-miR-1296	TBX21
hsa-miR-1908	TBX21	hsa-miR-1296	CCR7
hsa-miR-1908	IL2RB	hsa-miR-1296	CD247
hsa-miR-1908	CD27	hsa-miR-1296	IL2RG
hsa-miR-1908	CCR7	hsa-miR-3193	CD247
hsa-miR-1908	PDCD1	hsa-miR-3193	PDCD1
hsa-miR-663	CD40LG	hsa-miR-3193	IL2RG
hsa-miR-663	CD8A	hsa-miR-492	ZAP70
hsa-miR-663	TBX21	hsa-miR-492	ITK
hsa-miR-663	IL2RB	hsa-miR-492	CD5
hsa-miR-663	CD27	hsa-miR-492	IL2RB
hsa-miR-663	CCR7	hsa-miR-492	PDCD1
hsa-miR-663	PDCD1	BCL11A_GM12878_hg19	LCK
hsa-miR-3677-3p	TBX21	BCL11A_GM12878_hg19	IL2RB
hsa-miR-3677-3p	IL2RB	BCL11A_GM12878_hg19	CCR7
hsa-miR-3677-3p	CCR7	BCL11A_GM12878_hg19	PDCD1
hsa-miR-3677-3p	CD3E	–	–

Discussion

BC has been an international health concern due to its relatively high metastasis and recurrence rate. Despite considerable advancements in the early detection and clinical management of this condition, it remains a feared and refractory ailment; BC is the leading cause of tumor-related deaths among women. Amjad et al. discovered that a 53-mRNA signature may be a viable gene set for identifying TNM stages, as well as formulating prospective treatments; this signature may show a high performance in predicting and controlling recurrence-free survival rates at stages I, II, III, and IV, as well as in groups of patients at all stages (Amjad et al. 2020). Liu et al. found that TF-related signatures can independently predict the prognosis of BC patients without considering clinical variables (Liu et al. 2021). However, the molecular mechanism underlying BC metastasis has not been fully understood; this hinders the improvement of therapeutic interventions for BC. Therefore, in our study, comprehensive bioinformatics analyses were performed using data from TCGA and GEO datasets to screen promising new therapeutic targets to manage the development and metastasis of BC.

In this study, we identified 146 DEGs, including three upregulated DEGs and 143 downregulated DEGs, between primary BC and metastatic BC samples. Moreover, 56 DEGs were closely correlated with BC survival outcomes. Functional enrichment analysis showed that several genes, including CCR7, CD27, and IL2RB, were markedly enriched in pathways such as cytokine-cytokine receptor interaction and immune responses. Moreover, the PPI network and sub-module analyses demonstrated that downregulated genes (CCR7, CD27, and IL2RB) acted as hub genes with higher node degrees. Of note, based on the regulatory analysis, miR-1908 and miR-663 strongly interacted with these three genes.

CCR7 has been generally recognized to be responsible for tumor proliferation, migration, and invasion (Xu et al. 2017). Recently, a dearth of evidence has also confirmed that CCR7 is implicated with the occurrence, metastasis, and invasion of BC (Mishan et al. 2015; Sonbul et al. 2017). For example, Strien et al examined the CCR7 expression levels in different BC cells and observed that CCR7 is expressed in the overwhelming majority of primary BC tumors but is more closely related to aggressive tumors (Strien et al. 2017). Vehedi et al also highlighted that CCR7 shows a correlation with the stage and grade of BC (Vahedi et al. 2018). More interestingly, emerging research has demonstrated that CCR7 is predominately associated with BC prognosis and survival outcomes (Cassier et al. 2011). Herein, we found that this gene was downregulated and validated that it was a survival-related gene in BC patients. Moreover, functional analysis revealed that CCR7 primarily participated in cytokine-cytokine receptor interaction pathways. In addition, according to our network analysis, CCR7 had strong relationships with miR-663 and miR-1908. Rawlings-Goss et al previously showed that miR-1908 is a crucial regulator involved in the onset and development of tumor metastasis, indirectly implying that the CCR7/miR-1908 axis possibly participates in the metastasis of BC (Rawlings-Goss et al. 2014). Similarly, multiple investigators have argued that miR-663 plays significant roles in different cancers, such as colorectal cancer and prostate carcinoma (Tsikitis et al. 2016; Wang et al. 2017). Taken together, we inferred that CCR7 can be considered a promising mediator of BC metastasis and prognosis.

CD27, a member of the tumor necrosis factor family, can stimulate cell apoptosis and modulate immune lymphocyte function. Via an association analysis, Xu et al showed that the aberrant expression of CD27 could disrupt immune responses and influence BC risk (Xu et al. 2012). Our findings demonstrated that CD27, as a notable BC survival-associated gene, also played a prominent node in the PPI and sub-module networks, showing a high degree. Furthermore, functional enrichment analyses revealed that CD27 was significantly enriched in cytokine-cytokine receptor interaction pathways and immune responses. Accumulating evidence has implied that cytokines are mainly produced by immune cells and are engaged in intra- and intercellular communication (Gravestein and Borst 1998). Numerous researchers have pointed out that cytokines and their corresponding receptors participate in the metastasis of various cancers, including BC (Méndez-García et al. 2019); this suggests that CD27 may be involved in BC metastasis. Notably, a growing body of research has emphasized that investigating the crosstalk between the molecular footmarks of BC and the immune system serves as a promising therapeutic approach for BC. In addition, in the present study, strong interactions between CD27 and miR-1980/miR-663 were observed. Therefore, we speculated that CD27 is closely correlated with the metastasis and prognosis of BC.

IL2RB is considered a potential marker for the susceptibility to various cancers, such as lung cancer and BC (Marra et al. 2014; Ho et al. 2018). A previous study indicated that IL2RB is broadly expressed in the majority of BC cells (Ridings et al. 1995). Later, numerous investigators demonstrated that a higher expression of IL2RB is related to the progression of BC tumors, and more importantly, this alteration is possibly involved in the malignancy of BC tumors (García-Tuñón et al. 2004). Furthermore, Hong et al previously stated that the variants in IL2RB are associated with BC risk and this variant was mainly responsible for cytokine-correlated pathways (Hong et al. 2018). Consistently, our results suggested that IL2RB was significantly enriched in the cytokine-cytokine receptor interaction pathway. Additionally, we performed a survival analysis and found that this gene was predominately linked with the survival outcomes of BC, indicating that IL2RB may serve as a vital prognostic gene with regard to BC treatment. More intriguingly, we also observed that IL2RB closely interacted with two pivotal regulators (miR-1908 and miR-663) that were reported to be implicated with tumor progression. Thus, our results provided a novel avenue to improve the understanding of the potential molecular mechanisms underlying the role of IL2RB in BC development.

Although several prominent gene targets were screened using bioinformatics methods, our investigation had a few limitations. Firstly, according to the screening criteria, we only screened the GSE102818 data set. Therefore, more rigorous studies are required to verify our findings; these studies should primarily involve exhaustive bioinformatics analyses with a larger sample size and strong experimental evidence. Secondly, we found strong interactions between several key genes and specific signaling pathway; this may represent a crucial therapeutic strategy for BC. However, the mechanisms underlying these interactions and their roles remain to be elucidated. In addition, comprehensive clinical data should also be integrated further to evaluate and validate the prognosis of BC patients.

Conclusion

In summary, we identified vital genetic markers, such as CCR7, CD27, and IL2RB, which were correlated with the mechanisms underlying the metastasis of BC. Notably, we found that these genes were possibly associated with the survival rate of BC patients. Meanwhile, numerous pathways, including cytokine-related pathways, were also observed to play essential roles in the initiation and progression of BC. Additionally, miR-1908 and miR-663 were found to closely interact with these three genes. Thus, our findings may provide novel insights to formulate effective therapeutic interventions for the prevention and treatment of BC. However, further experimental confirmations and integrated bioinformatics analyses are required to validate the findings of the present study.

Data availability

The data that support the findings of this study are publicly available in the GEO database (https://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE102818.

Disclosure statement

No potential conflict of interest was reported by the author(s).