By 2018, 1,762,450 new tumor patients are been made definite diagnosed and 606,880 among them died in the United States. Of them, breast carcinoma makes up 30% (268,600) of all new female tumor cases and 15% (41,760) of all tumor-related deaths1. Although we have made great progress in prediagnosis and treatment, medicine resistance and molecular heterogeneity have influenced the active treatment of drugs. Every year, the increasing prevalence of breast carcinoma in the United States and the relevant treatment of malignant breast carcinoma bring huge pressure to the national economy and health care industry2. According to recent treatment views, like radiotherapy and adjuvant chemotherapy, have promoted the overall survival of patients, particularly in early-stage breast tumor of women who were diagnosed timely3. However, among these treatments, side effects like psychological effects, gastrointestinal discomfort, adverse physical, and psychological effects on sexuality often exist which influences the quality of life of patients. Besides, they can also take the possibility of mortality like bone loss, increased risk for secondary cancers and cardiac toxicity4,5. One of the significant advantages for new, molecularly targeted therapies is to lower the possibility of morbidity related to tumor also with mortality6. Approximately 70% of patients with breast carcinoma express the estrogen receptor (ER) and endocrine therapies are better for the treatment of ER+ breast carcinoma7. Recently, CDK4/6 inhibitors like palbociclib, ribociclib, and abemaciclib have led to a better disease-free survival compared with the antiestrogen alone in patients with advanced ER+ breast cancer (BC), and FDA has thus approved these drugs8,9. Even though plenty of patients with breast carcinoma treated with antiestrogens and CDK4/6 inhibitors gain benefits from this treatment plan, it is also important to investigate the resistance-associated gene. In the current study, microarray data for GSE117743 facilitated the investigation of differentially expressed genes (DEGs) in palbociclib-sensitive and palbociclib-resistant BC. The functions of the DEGs were assessed by using gene ontology (GO) annotation, Kyoto Encyclopedia of Genes and Genomes (KEGG), protein-protein interaction (PPI) network, and the relationship between palbociclib-resistant genes and overall survival of patients with breast carcinoma. In conclusion, we conducted this research to discover the key genes of medicine resistance and to find potential new cancer therapy targets to reduce palbociclib resistance.
The Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/gds) database is a public functional genomics dataset that helps users download statistics and achieve gene expression introduction. In our research, Gene expression profile data (GSE117743) was obtained from GEO. Three palbociclib-resistant BC cells and 3 palbociclib-sensitive cell lines were included. The array data were acquired from the Affymetrix Human Genome U133A 2.0 Array [GPL571; transcript (gene) version].
R software was utilized to compare 2 groups of tissues to identify different genes under identical experimental factors: different expression genes in resistant and sensitive cell lines. Besides, we used |log2FC|≥1 and P<0.05 as a cut-off criterion and an obvious statistical difference would be considered if the statistics met our standards10. In every experiment, the null hypothesis of the rank way was that each gene was assorted randomly; The gene would have topper ranking if the P-value was smaller. More importantly, the pheatmap package was used to draw the different genes.
Functional and pathway enrichment analysis
The KEGG pathway analysis and GO functional analysis of the DEGs we identified were conducted by taken advantage of the Database of Annotation Visualization and Integrated Discovery (DAVID) and R software. GO analysis was divided into the cellular component, biological process (BP), and molecular function, and a P-value of <0.05 was considered to indicate a statistical difference. In terms of the KEGG pathway analysis, they were identified according to the P-value of <0.05. DAVID includes many functional annotation programs to analyze plenty of biological information of genes and data that existed in DAVID was vital for the achievement of high-throughput gene functional analysis11. We also used R software to analyze datasets.
PPI network analysis
PPI network analysis was conducted by utilizing the Search Tool for the Retrieval of Interacting Genes (STRING). STRING supplies messages related to proven and predicted interactions among huge numbers of proteins. DEGs which were selected would then input into STRING. We would consider genes to be meaningful which had a score of >0.912. The PPI network was built using Cytoscape and Metascape. Besides, degrees >10 was set as the cut-off criterion. Molecular Complex Detection (MCODE) was then used to analyze the submodules of the PPI network, and the criteria we set was the number of nodes >6 and MCODE score >5.
Analysis of the hub genes and their relationship with breast carcinoma prognosis
The Kaplan-Meier Plotter (http://kmplot.com/analysis/) is a nonparametric statistic that is utilized to evaluate the overall survival statistics from recorded data. Also, every hub gene which was selected would be put into the online tool to evaluate the overall survival and progression-free survival of patients with BC for the Kaplan-Meier curve. This tool was built using gene expression and survival data of patients with BC, which were downloaded from GEO and The Cancer Genome Atlas (Affymetrix HG-U133A, HG-U133A 2.0 and HG-U133 plus 2.0 microarrays)13. As every gene mapped to a different Affy ID, each Affy ID was entered into the online tool to obtain the survival curves.
Identification of the genes between palbociclib-sensitive and palbociclib-resistant BC cells
The “t test” option in R software was used to research the gene expression profiles from the GSE117743 public microarray dataset. It highlighted the DEGs between GSM3308126, GSM3308127, GSM3308128 (palbociclib-sensitive), GSM3308129, GSM3308130, and GSM3308131 (acquired palbociclib-resistant) BC cells. According to the cut-off criteria (P<0.05 and |log2FC|≥2), 447 DEGs were identified, which were consisted of 380 downregulated and 67 upregulated genes. The top 20 upregulated and downregulated DEGs are shown in Figure 1.
GO analysis of DEGs
In the cellular component ontology, most of the enriched categories were associated with cell membrane components, including “cytoplasm” (145 genes), “plasma membrane” (121 genes), and “cytosol” (105 genes). Also, the other highly enriched terms included “extracellular exosome” (95 genes) and “the extracellular region” (70 genes). In the BP ontology, the highest number of DEGs were enriched in regulation-associated terms, including “signal transduction” (52 genes), “apoptotic process” (31 genes), “negative regulation of transcription from RNA polymerase II promoter” (27 genes), “immune response” (26 genes), and “inflammatory response” (26 genes). In the molecular function ontology, regulator-associated and binding-associated components accounted for the majority of the enriched GO categories, including “protein binding” (233 genes), “calcium ion binding” (32 genes), “protein homodimerization activity” (31 genes), “identical protein binding” (25 genes), and “receptor binding” (22 genes) (Table 1, Fig. 2).
KEGG pathways of DEGs
The KEGG pathway analysis identified many significant enriched pathways, including “IL-17 signaling pathways” (19 genes), “Herpes simplex infection” (15 genes), “Influenza A” (14 genes), “Hepatitis C” (14 genes), “NOD-like receptor signaling pathway” (14 genes), “Legionellosis” (8 genes), “Epstein-Barr virus infection” (8 genes), and “Steroid biosynthesis” (4 genes; Table 2). We also used barplot (Fig. 3A) and cnetplot (Fig. 3B) to show the enriched pathways.
Construction of a PPI network
The 447 DEGs were entered in the String website to obtain their PPI data. Then, if the combined score was ≥0.9, the PPIs were chosen to construct a PPI network (Fig. 4A). In the PPI network, 6 node proteins, including ADRA2A (α-2A adrenergic receptor), CYP2R1 (cytochrome P450 subfamily IIR polypeptide 1), CBS (Cystathionine β-synthase), NOD2 (nucleotide-binding oligomerization domain containing 2), EphA2 (erythropoietin-producing hepatocellular receptor A2), and ADM (adrenomedullin) had a high association with other node proteins (>7), which demonstrated their increased hub degrees (Fig. 4B). The hub genes may, therefore, serve a crucial role in palbociclib resistance within BC.
Relationship between hub genes and BC prognosis
The Kaplan-Meier Plotter was used to research the prognosis of patients with BC associated with genes we identified. Then, 4 hub genes were found to have a significant relationship with BC prognosis (Fig. 5). Kaplan-Meier survival analysis indicated that overexpression of the hub genes ADRA2A (Fig. 5B), CYP2R1 (Fig. 5C), and EphA2 (Fig. 5D) were related with poorer overall survival in patients with BC. However, overexpression of ADM (Fig. 5A) seemed was associated with better overall survival and the expression level of NOD2 (Fig. 5E) and CBS (Fig. 5F) may have no significant relationship with the overall survival of patients with BC. This suggested that the selected genes may be potential targets for palbociclib resistance in BC.
In this present research, gene expression data of GSE117743 were searched from the GEO database. Three palbociclib-resistant and 3 palbociclib-sensitive BC cell lines were identified for analysis. A total of 380 downregulated and 67 upregulated DEGs were identified among palbociclib-resistant and palbociclib-sensitive BC. Function annotation indicated that these DEGs were primarily related to cell membrane components, regulation-associated terms, regulator-associated and binding-associated components. KEGG pathway analysis of DEGs revealed involvement in “IL-17 signaling pathways,” “Herpes simplex infection,” “Influenza A,” “Hepatitis C,” and “NOD-like receptor signaling pathway.” According to previous studies, IL-17 signaling pathways promoted resistance to paclitaxel in breast tumors through activation of the ERK1/2 pathway14. Also, in terms of the “NOD-like receptor signaling pathway,” a previous study reported that it may be related to a higher risk of certain cancer types like colorectal cancer, gastric cancer, lymphoma, and as well as BC15. Therefore, the above KEGG pathways mentioned in this study may assist us to find new diagnostic and therapeutic ways in palbociclib-resistant BC. Also, PPI network analysis indicated that 6 hub genes may be used as new targets in palbociclib-resistant BC which had higher degrees of interaction.
The first gene ADRA2A (α-2A adrenergic receptor), encodes a protein receptor of 450 amino acid residues16 and locates in chromosome 10q25.2. Almost all tissues and organs in our body express it. It is also involved in some functions in the platelet aggregation, cardiovascular system, central nervous system, lipolysis, blood pressure, insulin secretion, and neurotransmitter release17. Besides, patients with pancreatic cancer who have a lower expression of this receptor-related with higher insulin secretion18. it also regulates adrenergic suppression of insulin secretion and associated with poor prognosis of breast carcinoma: advanced TNM stages, larger tumor size, positive lymph node involvement, higher Her2, and nottingham prognostic index positive status by learning from previous studies19,20.
CBS (Cystathionine β-synthase), the second gene, encodes a protein of 551 amino acids which is almost expressed in the liver, brain, kidney, and pancreas. Overexpression of CBS in cancer tissues or cell lines has been found in ovarian21, colon22,23, prostate24, and BC25, in contrast with normal tissues or nontransformed cell lines. On the one hand, CBS promotes cancer progression and growth, on the other hand, it can also initiate cancer formation. More than that, lower expression of CBS takes down antioxidant capacity and promotes the sensitivity of tumor cells to chemotherapeutic drugs. The cytoprotective influence of CBS is also related to the regulation of p53 apoptosis-related and NF-κB signaling pathways26. With the identification of the function of CBS in the tumor, it will become more significant to take advantage of CBS as a predictive and prognostic biomarker.
Also, the third gene CYP2R1 (cytochrome P450 subfamily IIR polypeptide 1), is a key enzyme involved in vitamin D 25 hydroxylation27. It encodes microsomal vitamin D 25-hydroxylase and catalyzes the first hydroxylation reaction in the biosynthesis of vitamin D28. According to further studies, it was noticed by a protective effect of this single nucleotide polymorphisms against more advanced disease, but we need more researches to find the mechanism in palbociclib resistance cell lines29.
EphA2 (erythropoietin-producing hepatocellular receptor A2), the fourth gene, plays a role in critical BPs, like cancer progression, metastasis, and angiogenesis. Many further studies researching EphA2 targets were applied to the treatment of tumors. According to recent researches, we identified the 64CU–DOTA-1C1 antibody in colorectal cancer, EA2, B233, 3F2-WT antibody in BC and EA5 antibody in ovarian cancer30–32. In clinical researches, EphA2 was identified to be overexpressed in many malignant carcinomas and might be used as a prognostic marker of poor prognosis in patients with tumors, except for lung cancer33. In terms of researches on drug resistance, one study believed that the EphA2-to-YAP pathway promoted gastric cancer growth and therapy resistance34, and another study found that EphA2 inhibition combined with eicosapentaenoic acid has great efforts in preclinical models of triple-negative breast tumor by disrupting cellular cholesterol efflux35.
Another huge gene, NOD2 (nucleotide-binding oligomerization domain containing 2), is a member of evolutionarily conserved Nod-like receptors family which share a tripartite structure of a terminal sensor domain (leucine-rich repeats, LRRs), a central nucleotide-binding oligomerization domain (NOD)36. NOD2 has already been certified to have a relationship with a high risk of breast tumor37–39, reduces estrogen-induced proliferative responses and promotes apoptosis in BC40. The previous analysis indicated that overexpression of NOD2 disrupts immune-related pathways, especially those signaling through NF-κB which are activated via TNF signaling pathways and antigen presentation signals (HLA class I histocompatibility proteins)41. Also, another study reported that lower expression of NOD2 promotes both cellular adhesion/migration and proliferation via MAPK signaling pathways and PI3K/Akt/mTOR signaling pathways42.
Furthermore, the relationship between the mentioned hub genes and the overall survival curves for patients with BC revealed that overexpression levels of the hub genes ADRA2A, CYP2R1, and EphA2 were related to worse prognosis in patients with BC. Besides, overexpression of ADM seemed was associated with better overall survival. However, the expression level of NOD2 and CBS may have no significant relationship with overall survival according to the above database.
These findings approved that these hub genes may play a vital role in palbociclib resistance to BC. However, the molecular reasons for palbociclib resistance, and whether these huge genes can be used as targets for the treatment of palbociclib-resistant BC remain unknown. Since all our data were achieved from the GEO database by bioinformatics tools, as well as the limited number of relevant BC cell lines and lack of multiple database analysis, more data analysis, and clinical experiments should be further performed to certify its feasibility.
Our study concluded certain mechanisms of palbociclib resistance in breast carcinoma. Bioinformatics methods were utilized to identify the DEGs in palbociclib-resistant BC cells. Also, 6 hub genes were identified as potential targets for the treatment of palbociclib resistance. These conclusions from the present research may offer a significant understanding of the mechanisms of palbociclib resistance in BC. However, these findings may be considered for future investigation.
Sources of funding
G.G. and X.S. contributed equally to this work.
Conflict of interest disclosure
The authors declare that they have no financial conflict of interest with regard to the content of this report.
Research registration unique identifying number (UIN)
1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin 2019;69:7–34.
2. Montero AJ, Eapen S, Gorin B, et al. The economic burden of metastatic breast cancer: a U.S. managed care perspective. Breast Cancer Res Treat 2012;134:815–22.
3. Zhang J, Wang Q, Wang Q, et al. Mechanisms of resistance to estrogen receptor modulators in ER+/HER2− advanced breast cancer. Cell Mol Life Sci 2019. Doi: 10.1007/s00018-019-03281-4.
4. Qiu N, He YF, Zhang SM, et al. Cullin7 enhances resistance to trastuzumab therapy in Her2 positive breast cancer via degrading IRS-1 and downregulating IGFBP-3 to activate the PI3K/AKT pathway. Cancer Lett 2019;464:25–36.
5. Ronkina N, Shushakova N, Tiedje C, et al. The role of TTP phosphorylation in the regulation of inflammatory cytokine production by MK2/3. J Immunol 2019. Doi: 10.4049/jimmunol.1801221.
6. da Silva Correia J, Miranda Y, Leonard N, et al. Regulation of Nod1-mediated signaling pathways. Cell Death Differ 2007;14:830–9.
7. Santoni M, Occhipinti G, Romagnoli E, et al. Different cardiotoxicity of palbociclib and ribociclib in breast cancer: gene expression and pharmacological data analyses, biological basis, and therapeutic implications. BioDrugs 2019. Doi: 10.1007/s40259-019-00382-1.
8. Gelmon KA, Cristofanilli M, Rugo HS, et al. Efficacy and safety of palbociclib plus endocrine therapy in North American women with hormone receptor-positive/human epidermal growth factor receptor 2-negative metastatic breast cancer. Breast J 2019. Doi: 10.1111/tbj.13516.
9. Yu E, Xu Y, Shi Y, et al. Discovery of novel natural compound inhibitors targeting estrogen receptor alpha by an integrated virtual screening strategy. J Mol Model 2019;25:278.
10. Li L, Wang G, Li N, et al. Identification of key genes and pathways associated with obesity in children. Exp Ther Med 2017;14:1065–73.
11. Xu Z, Zhou Y, Cao Y, et al. Identification of candidate biomarkers and analysis of prognostic values in ovarian cancer by integrated bioinformatics
analysis. Med Oncol 2016;33:130.
12. Zhao B, Baloch Z, Ma Y, et al. Identification of potential key genes and pathways in early-onset colorectal cancer through bioinformatics
analysis. Cancer Control 2019;26:1073274819831260.
13. Ikeda N, Nakajima Y, Tokuhara T, et al. Clinical significance of aminopeptidase N/CD13 expression in human pancreatic carcinoma. Clin Cancer Res 2003;9:1503–8.
14. Laprevotte E, Cochaud S, du Manoir S, et al. The IL-17B-IL-17 receptor B pathway promotes resistance to paclitaxel in breast tumors through activation of the ERK1/2 pathway. Oncotarget 2017;8:113360–72.
15. Fan W, Gao X, Ding C, et al. Estrogen receptors participate in carcinogenesis signaling pathways by directly regulating NOD-like receptors. Biochem Biophys Res Commun 2019;511:468–75.
16. Kobilka BK, Matsui H, Kobilka TS, et al. Cloning, sequencing, and expression of the gene coding for the human platelet alpha 2-adrenergic receptor. Science 1987;238:650–66.
17. Affolder T, Akimoto H, Akopian A, et al. Search for neutral supersymmetric Higgs bosons in pp collisions at sqrt[s] = 1.8TeV. Phys Rev Lett 2001;86:4472–8.
18. Rosengren AH, Jokubka R, Tojjar D, et al. Overexpression of alpha2A-adrenergic receptors contributes to type 2 diabetes. Science 2010;327:217–220.
19. Kaabi B, Belaaloui G, Benbrahim W, et al. ADRA2A germline gene polymorphism is associated to the severity, but not to the risk, of breast cancer. Pathol Oncol Res 2016;22:357–65.
20. Perez Pinero C, Bruzzone A, Sarappa MG, et al. Involvement of alpha2- and beta2-adrenoceptors on breast cancer cell proliferation and tumour growth regulation. Br J Pharmacol 2012;166:721–36.
21. Bhattacharyya S, Saha S, Giri K, et al. Cystathionine beta-synthase (CBS) contributes to advanced ovarian cancer progression and drug resistance. PLoS One 2013;8:e79167.
22. Phillips CM, Zatarain JR, Nicholls ME, et al. Upregulation of cystathionine-beta-synthase in colonic epithelia reprograms metabolism and promotes carcinogenesis. Cancer Res 2017;77:5741–54.
23. Szabo C, Coletta C, Chao C, et al. Tumor-derived hydrogen sulfide, produced by cystathionine-beta-synthase, stimulates bioenergetics, cell proliferation, and angiogenesis in colon cancer. Proc Natl Acad Sci U S A 2013;110:12474–9.
24. Guo H, Gai JW, Wang Y, et al. Characterization of hydrogen sulfide and its synthases, cystathionine beta-synthase and cystathionine gamma-lyase, in human prostatic tissue and cells. Urology 2012;79:483 e481–483 e485.
25. Sen S, Kawahara B, Gupta D, et al. Role of cystathionine beta-synthase in human breast cancer. Free Radic Biol Med 2015;86:228–38.
26. Ferraro E, Trapani D, Marrucci E, et al. Evaluating triptorelin as a treatment option for breast cancer. Expert Opin Pharmacother 2019;20:1809–118.
27. Lasky-Su J, Lange N, Brehm JM, et al. Genome-wide association analysis of circulating vitamin D levels in children with asthma. Hum Genet 2012;131:1495–505.
28. Wang TJ, Zhang F, Richards JB, et al. Common genetic determinants of vitamin D insufficiency: a genome-wide association study. Lancet 2010;376:180–8.
29. Carvalho IS, Goncalves CI, Almeida JT, et al. Association of vitamin D pathway genetic variation and thyroid cancer. Genes (Basel) 2019:10.
30. Cai W, Ebrahimnejad A, Chen K, et al. Quantitative radioimmunoPET imaging of EphA2 in tumor-bearing mice. Eur J Nucl Med Mol Imaging 2007;34:2024–36.
31. Coffman KT, Hu M, Carles-Kinch K, et al. Differential EphA2 epitope display on normal versus malignant cells. Cancer Res 2003;63:7907–12.
32. Landen CN Jr, Lu C, Han LY, et al. Efficacy and antivascular effects of EphA2 reduction with an agonistic antibody in ovarian cancer. J Natl Cancer Inst 2006;98:1558–70.
33. Shen W, Xi H, Zhang K, et al. Prognostic role of EphA2 in various human carcinomas: a meta-analysis of 23 related studies. Growth Factors 2014;32:247–53.
34. Huang C, Yuan W, Lai C, et al. EphA2-to-YAP pathway drives gastric cancer growth and therapy resistance. Int J Cancer 2019. Doi: 10.1002/ijc.32609.
35. Torres-Adorno AM, Vitrac H, Qi Y, et al. Eicosapentaenoic acid in combination with EPHA2 inhibition shows efficacy in preclinical models of triple-negative breast cancer by disrupting cellular cholesterol efflux. Oncogene 2019;38:2135–50.
36. Deng W, Xie J. NOD2 signaling and role in pathogenic mycobacterium recognition, infection and immunity. Cell Physiol Biochem 2012;30:953–63.
37. Krishna NR, Krushna NS, Narayanan RB, et al. Expression, purification and characterization of refolded rBm-33 (pepsin inhibitor homolog) from Brugia malayi: a human Lymphatic Filarial parasite. Protein Expr Purif 2011;79:245–50.
38. Liu J, He C, Xu Q, et al. NOD2 polymorphisms associated with cancer risk: a meta-analysis. PLoS One 2014;9:e89340.
39. Song F, Gilbody S. Bias in meta-analysis detected by a simple, graphical test. Increase in studies of publication bias coincided with increasing use of meta-analysis. BMJ 1998;316:471.
40. Worapamorn W, Haase HR, Li H, et al. Growth factors and cytokines modulate gene expression of cell-surface proteoglycans in human periodontal ligament cells. J Cell Physiol 2001;186:448–56.
41. Keestra-Gounder AM, Tsolis RM. NOD1 and NOD2: beyond peptidoglycan sensing. Trends Immunol 2017;38:758–67.
42. Velloso FJ, Campos AR, Sogayar MC, et al. Proteome profiling of triple negative breast cancer cells overexpressing NOD1 and NOD2 receptors unveils molecular signatures of malignant cell proliferation. BMC Genomics 2019;20:152.