Abstract
Objective. To identify trastuzumab-resistant genes predicting drug response and poor prognosis in human epidermal growth factor receptor 2 positive (HER2+) breast cancer. Methods. Gene expression profiles from the GEO (Gene Expression Omnibus) database were obtained and analyzed. Differentially expressed genes (DEGs) between the pathological complete response (pCR) group and non-pCR group in a trastuzumab neoadjuvant therapy cohort and DEGs between Herceptin-resistant and wild-type cell lines were detected and evaluated. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses were performed to select the functional hub genes. The hub genes’ prognostic power was validated by another trastuzumab adjuvant treatment cohort. Results. Fifty upregulated overlapping DEGs were identified by analyzing two trastuzumab resistance-related GEO databases. Functional analysis picked out ten hub genes enriched in mitochondrial function and metabolism pathways: ASCL1, CPT2, DLD, ELVOL7, GAMT, NQO1, SLC23A1, SPR, UQCRB, and UQCRQ. These hub genes could distinguish patients with trastuzumab resistance from the sensitive ones. Further survival analysis of hub genes showed that DLD overexpression was significantly associated with an unfavorable prognosis in HER2+ breast cancer patients. Conclusion. Ten novel trastuzumab resistance-related genes were discovered, of which DLD could be used for trastuzumab response prediction and prognostic prediction in HER2+ breast cancer.
1. Introduction
HER2 (human epidermal growth factor receptor 2) positive breast cancer accounts for 20-30% of the total primary breast cancer population [1]. HER2 positive (HER2+) patients were traditionally correlated with the poor outcome before the discovery of trastuzumab, a monoclonal antibody targeting to suppress HER2 activity [2]. Until now, trastuzumab is still the first-line therapeutic drug for treating HER2+ breast cancer patients [3]. Trastuzumab can immensely increase the clinical prognosis of primary HER2+ breast cancer and metastatic HER2+ breast cancer [4]. However, the resistance of trastuzumab is a major impediment to impairing its efficacy for patients. So, further identification of genomic alterations of trastuzumab insensitivity is crucial to provide potential biomarkers for precise treatment of trastuzumab and to dig novel therapeutic strategies.
Gene profiling and signatures are utilized to quickly detect differentially expressed genes (DEGs), which have significantly promoted the progress of tumor research during the past decades. A few studies have analyzed the predictive and prognostic value of transcript sequencing-based genes from breast cancer patients. Through analysis of genome-wide expression profiling, Sotiriou et al. found out that topoisomerase II-alpha (TOP2A) expression could effectively predict the response of anthracycline (epirubicin) for estrogen receptor– (ER–) negative breast cancer patients who received neoadjuvant chemotherapy (NCT) [5]. Additionally, a 4-microRNAs-based signature was found to be of significant value in predicting the pathological complete response (pCR) for basal-like triple negative breast cancer (TNBC) breast cancer patients [6], and a 17-gene signature was found to have an excellent predictive and prognostic role in TNBC patients who receive anthracycline-based NCT [7]. However, a few studies were conducted to identify gene signatures predicting trastuzumab sensitivity for HER2+ patients. Hence, we performed a profound analysis using integrated bioinformatic methods for massive public data. Ultimately, we identified essential genes that could predict the therapeutic sensitivity of trastuzumab in HER2+ patients, which is of great clinical significance.
With the widespread application of high-throughput techniques, GEO (Gene Expression Omnibus) provides a platform that includes millions of datasets, tissues, and cell samples to analyze public gene expression. This allowed gene analysis from discovering molecular biomarkers and classifying disease by comparing phenotypes. Our study first calculated the DEGs by drug response in a trastuzumab neoadjuvant treatment cohort in GEO databases. We also found another DEG cluster between Herceptin-resistant and wild-type cell lines. We identified an upregulated gene set with DEGs based on the overlapping analysis. Several functional evaluations revealed that 10 DEGs were related to mitochondria, identified as hub genes. The survival analysis of the hub genes in another trastuzumab adjuvant cohort demonstrated that only DLD overexpression was significantly associated with poor outcomes. Further survival and immune studies explored that DLD may be discovered as a biomarker for trastuzumab response prediction and HER2+ breast cancer prognosis assessment.
2. Methods
2.1. Data Collection
The NCBI- (National Center for Biotechnology Information-) GEO datasets GSE62327, GSE15043, and GSE58984 were searched using “trastuzumab” and “breast cancer” as keywords. GSE62327 contains the RNA-seq data and clinical information, especially drug responses to different neoadjuvant therapy. We chose 6 patients with pCR and 18 patients with non-pCR in patients who used trastuzumab as neoadjuvant therapy. GSE15043 comprises the RNA-seq data of 8 Herceptin-resistant breast cancer cell lines and 2 wild-type breast cancer cell lines. GSE58984 has the RNA-seq data and clinical information from a cohort of 94 patients who used trastuzumab as adjuvant therapy. All microarray data were obtained from the GEO database: https://www.ncbi.nlm.nih.gov/geo/. The raw data were downloaded as MINiML files.
The RNA-sequencing data of all upregulated genes and corresponding clinical information were downloaded from the cancer genome atlas (TCGA) dataset (https://portal.gdc.com) in HER2+ breast cancer despite the status of hormone receptors. Finally, 185 samples were chosen. Counting data is converted to transcripts per million (TPM) and normalized log2 (TPM+1) while keeping clinical information intact.
2.2. Identification of DEGs (Differentially Expressed Genes)
The “limma” R package was used to compare the genome expression profile of GSE62327 and GSE15043. In GSE62327, DEGs between the pCR and non-pCR group with a false discovery rate (FDR) <0.05 were picked. In GSE15043, DEGs between wild-type and Herceptin-resistant breast cancer cell lines with an FDR <0.05 were chosen. DEGs associated with overall survival (OS) were assessed using a univariate Cox proportional hazard regression model analysis. Volcano plots and forest plots were constructed using ggplot2 packages. The overlapping DEGs were drawn in a Venn diagram using the “VennDiagram” R package.
2.3. Functional Enrichment Analysis
The “clusterProfiler” R package was used to conduct Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis, and GO and KEGG pathways’ network figure. The network’s top nodes with the most connectivity lines were identified as hub genes.
2.4. Construction of PPI (Protein-Protein Interaction) Network
The PPI networks for overlapping DEGs were conducted by the STRING database (version 11.5) and Cytoscape software (version 3.9.1). The top nodes with mitochondrial function were identified as hub genes.
2.5. Mitochondrial Proteins and Pathways Annotation
The annotations of mitochondrial-related genes were interpreted by MitoCarta3.0 downloaded from https://pubs.broadinstitute.org/mitocarta/mitocarta30-inventory-mammalian-mitochondrial-proteins-and-pathways.
2.6. Validation and Screening of Hub Genes
To confirm the prognostic effect of hub genes in patients using trastuzumab as adjuvant therapy, survival analyses were carried out in GSE58984 using the hub genes’ profile. The distant disease-free survival (DDFS) of each hub gene was calculated by Kaplan-Meier analysis with the log-rank test. The Kaplan-Meier curves were drawn using the “surviminer” R package.
2.7. Estimate the Proportion of Immune and Cancer Cells (EPIC)
Estimation of immune fractions of DLD was determined through HER2+ breast tissue in TCGA database using EPIC by R package “GfellerLab/EPIC.” The student -test was calculated and visualized using R software.
2.8. Statistical Analysis
Wilcoxon rank-sum test was used to compare the gene expression between two groups. Kruskal-Wallis test was used to compare the gene expression among different groups. The OS, disease-specific survival (DDS), and progress-free interval (PFI) between other groups were calculated by Kaplan-Meier analysis with the log-rank test. Spearman’s correlation analysis was used to describe the correlation between quantitative variables. All statistical analyses were performed with R software (Version 4.1.3). All values were two-tailed unless otherwise specified, and values less than 0.05 were considered statistically significant.
3. Results
3.1. Identify DEGs in Trastuzumab Neoadjuvant Therapy Cohort and Herceptin-Resistant Cells
Twenty-four patients from GSE62327 were assigned to the “pCR” group () and “non pCR” group () by the response to trastuzumab monotherapy as neoadjuvant therapy at the pathological level (Figure 1). The detailed clinical information of these patients was summarized (Supplementary Table 1). In another queue, 10 cell lines from GSE15043 were divided into the “wild-type” group () and the “Herceptin-resistant” group () (Figure 1). We compared the whole-genome expression profile between groups in each cohort. A total of 399 DEGs from GSE62327 and 2250 DEGs from GSE15043 were detected, respectively (Table 1). The volcano plots of the DEGs were presented in Figure 2(a) (). The Venn diagrams (Figure 2(b)) obtained 50 upregulated DEGs and 38 downregulated DEGs in common. Most upregulated overlapping DEGs (31/50, 62%) were overexpressed in tumor tissues compared to adjacent nontumorous tissues in TCGA datasets (Figure 2(c)). In turn, the majority of the downregulated overlapping DEGs (31/38, 81.56%) had a lower expression level in tumor tissues (Figure 2(c)). The univariate Cox proportional hazard regression analysis of all overlapping DEGs was performed (Figures S1(a)–S1(b)). The trastuzumab resistance-related DEGs did not provide a unified prognostic value in each group. There were more significant results in the upregulated group than downregulated ones (Figure 2(d)), so we focused our attention on upregulated genes.
(a)
(b)
(c)
(d)
3.2. Conduct Functional Analysis of Upregulated DEGs and Find out Hub Genes
GO enrichment analysis and KEGG pathway analysis was carried out to elucidate the biological functions and pathways of the upregulated DEGs. The GO analysis demonstrated that upregulated DEGs were significantly enriched in mitochondrial-related terms (Figure 3(a)), such as mitochondrial electron transport, respiratory chain complex III, and electron transfer activity. According to the GO analysis, the KEGG pathway analysis of upregulated DEGs (Figure 3(b)) also proved that they participated in many energy metabolism pathways. The protein-protein interaction (PPI) network among the DEGs also indicated that mitochondrial-related genes play hub gene roles (Figures S2(a)–S2(b)). To further investigate which part of genes plays the most critical role, the GO and KEGG network with nodes and connectivity lines was performed (Figure 3(c)). The top-ranked linker nodes were ASCL1, CPT2, DLD, ELVOL7, GAMT, NQO1, SLC23A1, SPR, UQCRB, and UQCRQ, which were defined as hub genes cluster (Figure S2C). The enrichment ID table of hub genes also demonstrated the recurring keyword “mitochondrial” (Table 2). To explore the mitochondrial-related upregulated DEGs’ function, we searched the hub genes in MitoCarto 3.0 database (Supplementary Table 2). Interestingly, we found that these genes are widely related to metabolism, mitochondrial central dogma, mitochondrial dynamics and surveillance, OXPHOS, etc. During all related genes, DLD had been mostly labeled (Figure S2(d)).
(a)
(b)
(c)
3.3. Validate Prognostic Effect of the Hub Genes in Trastuzumab Adjuvant Therapy Cohort
We further tested their prognostic robustness in an adjuvant treatment cohort to determine if the hub genes could predict the prognosis. In GSE58984, 94 patients received trastuzumab as adjuvant treatment (Supplementary Table 2). The cohort also offered DDFS information, allowing us to test hub genes’ prognostic value (Figures 4(a)–4(j)). Among ten hub genes, high expression of ACSL1, DLD, ELVOL7, and SPR increased the risk of distant disease relapse. To our interest, only patients with DLD overexpression showed a significantly reduced DDFS compared to the low expression group (, ).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.4. Explore DLD in Breast Cancer and HER2 Overexpression Subtype
To make in-depth knowledge of DLD, we further compared breast cancer tissues and adjacent normal tissues in TCGA breast cancer patients (Figure 5(a)). The DLD expression exceeded in the adjacent normal tissue than in tumorous tissues. We also made correlation analysis stratified by DLD expression level in the T stage, N stage, M stage, pathological stage, histological type, ER status, HER2 status, PAM50 subtypes, menopause status, and age (Table 3). As a result, DLD expression level was significantly associated with M stage, age, and histological type (Figures 5(b)–5(d) and S3(a)–S3(e)). In survival analysis of the TCGA cohort, the high expression of DLD was associated with an increased risk of death and relapse (Figures 5(e) and S3(f)). We also performed a survival analysis of DLD, especially in those with HER2 overexpression (Figures 5(f) and S3(g)). In particular, the survival analyses showed that DLD was a more sensitive index in the HER2 overexpression subtype. Moreover, we also evaluated the immune landscape of DLD in the EPIC algorithm. EPIC estimated scores revealed that NK cells had the most significant correlation with DLD expression in HER2+ BRCA (Figure 5(g)). DLD had a negative correlation with NK cells (, ), while CD4+ T cells (, ) and CD8+ T cells (, ) had a positive correlation (Figure 5(h)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4. Discussion
The widespread application of trastuzumab for HER2+ breast cancer patients significantly increased patients’ metastatic status and prognosis. However, a considerable number of people will develop trastuzumab insensitivity or resistance. Thus, our study first discovered a 10-gene signature that could predict the pCR rate of patients who received trastuzumab as a neoadjuvant therapeutic drug. Moreover, we also found that one of 10 genes, DLD, could be the predictive factor for the DDFS of HER2+ patients who received trastuzumab as postoperative adjuvant therapy. Totally, our study is the first one to uncover the trastuzumab-related gene both for neoadjuvant and adjuvant trastuzumab therapy.
Our study analyzed one GEO database based on the HER2+ breast cancer tissues of patients who received single trastuzumab therapy and one GEO database based on HER2+ cell lines, which are trastuzumab-resistant. Ultimately, we identified 86 DEGs with 50 upregulated genes and 38 downregulated genes in all three GEO databases. Then, we performed GO enrichment analysis and KEGG pathway analysis using genes upregulated. The pathways shown in GO enrichment were oxidoreductase complex, mitochondrial ATP synthesis coupled electron transport, respiratory chain complex, and mitochondrial respirasome pathways. And KEGG pathway analysis involved the metabolic pathway, fatty acid pathway, and oxidative phosphorylation (OXPHOS). Tan et al. reported that glycolysis is negatively linked to trastuzumab sensitivity [8]. Yan et al. found that OXPHOS enhancement is associated with trastuzumab resistance [9]. Our work disclosed fatty acid pathway is involved in trastuzumab resistance. Other studies also reported that fatty acid metabolism-associated proteins overexpressed in HER2-positive breast cancer cell lines and tumor sample [8]. Still, the exact function and potential molecular mechanisms of fatty acid contributing to trastuzumab resistance are worth and necessary for profound investigation.
Additionally, 10 hub DEGs were further determined. DLD (dihydrolipoamide dehydrogenase) was negatively associated with pCR rate in HER2+ patients receiving trastuzumab neoadjuvant therapy and DDFS in HER2+ patients receiving trastuzumab as postoperative adjuvant therapy. So, DLD might be a trastuzumab resistance-related gene, and its function and underlying mechanism in trastuzumab resistance need further exploration. Besides, high DLD expression could predict the shorter overall survival in HER2+ breast cancer patients. The immune analysis also revealed that the immune microenvironment is relevant to DLD expression.
Based on these results, we speculate that DLD possibly is the key regulated gene in HER2+ breast cancer patients. DLD, a class-I pyridine nucleotide-disulfide oxidoreductase family member, is a mitochondrial enzyme. It is responsible for decarboxylating pyruvate to form acetyl-CoA during glucose metabolism and the production of mitochondrial adenine triphosphate (ATP) [10]. In eukaryotes, DLD plays an essential role in energy metabolism. DLD variants could cause dihydrolipoamide dehydrogenase deficiency (DLDD). People who suffer from DLDD have lactic acidosis and neurologic deterioration due to oxidative metabolism defects [11]. There only exist a few studies about the superficial effect of DLD in cancers. For example, downregulation of DLD in melanoma could inhibit cell proliferation by regulating energy metabolism [12], while DLD overexpression in head and neck cancer (HNC) cells could induce ferroptosis [13]. However, in DLD, as a mitochondrial-related gene, its roles in breast cancer are poorly understood. Prior studies have noted that mitochondrial dysfunction and dynamics are closely correlated with breast cancer progression and chemosensitivity [14–17]. A strong relationship between ATP synthase and resistance to HER2-targeted antibody therapies has been reported [9]. However, there is no research exploring the role of DLD in breast cancer, let alone HER2+ breast cancer. Therefore, exploring how DLD regulates mitochondrial energy metabolism to mediate trastuzumab resistance in HER2+ patients is meaningful.
To interpret our results, two limitations still need to be considered. First, our study only enrolled one GEO dataset, which includes HER2+ breast cancer patients receiving trastuzumab as neoadjuvant therapy, and one GEO dataset, which includes trastuzumab-resistant cells and normal cells. We need to enroll more GEO datasets to dig trastuzumab-resistant genes. Second, patients in the GEO dataset who received postoperative adjuvant therapy did not receive trastuzumab as a single therapeutic drug and lacked more survival-related information.
Taken together, our present work elucidates a novel gene signature that can predict the pCR rate of HER2+ breast in neoadjuvant therapy and the prognosis of HER2+ breast cancer in adjuvant treatment simultaneously. These hub genes may influence trastuzumab resistance by regulating mitochondrial-related metabolism. Moreover, our study highlighted the gene DLD as a diagnostic and prognostic factor and a potential target for HER2+ breast cancer treatments.
Data Availability
The data used to support the findings of this study are included within the article and supplementary files.
Conflicts of Interest
The authors declare that they have no competing interests.
Authors’ Contributions
(I) Conception and design were contributed by Yueyao Du, Zhenfeng Yu, and Huijuan Dai. (II) Administrative support was contributed by Yueyao Du. (III) Data analysis and interpretation were conducted by Xinrui Dong. (IV) Manuscript writing was contributed by Xinrui Dong and Huijuan Dai. (V) All authors did the final approval of the manuscript. (VI) Xinrui Dong and Huijuan Dai contributed equally to this research. Xinrui Dong, Huijuan Dai, and Aijun Sun contributed equally to this work.
Acknowledgments
The conduct of this study was sponsored and funded by the National Natural Science Foundation of China (grant number 82002777) and the Multidisciplinary Cross Research Foundation of Shanghai Jiao Tong University (Grant numbers YG2019QNA26) of Yueyao Du.
Supplementary Materials
Figure S1. Forest plots to demonstrate the results of the univariate Cox regression analysis between DEGs expression and OS. Figure S2. PPI network of DEGs. Figure S3. Correlation between clinicopathological characteristics and the expression of DLD and survival analysis of DLD in TCGA. Supplementary Table 1. Clinical characteristics of the patients in GSE62327. Supplementary Table 2. MitoPathway analysis of mito-related genes in upregulated DEGs. Supplementary Table 3. Clinical characteristics of the patients in GSE58984. (Supplementary Materials)