Screening of Hub Genes related to NAC resistance in Breast Cancer
To explore the common mechanism that promotes chemoresistance in breast cancer, we generated EPI-resistant cell lines of different subtypes of breast cancer, including MCF-7/EPI, SKBR3/EPI, and MDA-MB-231/EPI. The IC50 of EPI in resistant cells was seven times higher than that of parental wild-type (WT) cells (Fig. 1A and Table 1). Furthermore, consistent with clinical experience, tumor chemoresistance showed characteristics of multidrug resistance in our experiment. Breast cancer drug-resistant cells induced by EPI (a cytotoxic agent that inhibits RNA and DNA synthesis) are also resistant to DOC (a cytotoxic agent that stabilizes the tubulin polymer bundles and therefore interferes with microtubular assembly and cell replication). We then compared the expression profiles of EPI-resistant cells with WT cells. Three hundred and two genes were commonly upregulated or downregulated at least 2-fold in MCF-7/EPI, SKBR3/EPI, and MDA-MB-231/EPI compared to that in parental MCF-7, SKBR3, and MDA-MB-231 cells (Fig. 1B and C). In GO analysis, we analyzed the top 30 GO with the highest enrichment in cellular component (CC), biological process (BP), and molecular function (MF). Upregulated DEGs were widely distributed in the intracellular, intracellular part, intracellular organelle, organelle, and nucleus of the breast cancer resistant cells, and were enriched in “regulation of cellular metabolic process,” “negative regulation of cellular process” and “negative regulation of biological process.” The molecular functions of the upregulated DEGs were “protein and transcription regulatory region DNA” and “RNA polymerase II proximal promoter sequence-specific DNA binding.” Different to upregulated DEGs, downregulated DEGs were enriched in extracellular of the breast cancer resistant cells and were enriched in “cell adhesion” and “regulation of cell motility” (Supplementary material 1). The results showed that slowing down the cell cycle and decreasing biological processes and metabolic abnormalities are important mechanisms for the survival of chemo-resistant cells. In KEGG enrichment analysis, upregulated chemoresistance cell feature genes were enriched in “DNA replication,” “HIF-1 signaling pathway,” “p53 signaling pathway, and “Pentose phosphate pathway.” Downregulated chemoresistance cell feature genes were enriched in “Glutathione metabolism” and “Chemokine signaling pathway” (Supplementary material 1). The enrichment of these functions suggested that breast cancer chemo-resistant cells may resist chemotherapy by slowing down the cell cycle and strengthening DNA repair and synthesis. In order to more comprehensively characterize the biological characteristics of breast cancer drug resistance, we analyzed the hallmarks of breast cancer chemo-resistant cells and tumor tissues of breast cancer chemo-resistant patients by GSEA. We found that E2F and MYC targets, mTORC1 signaling, P53 pathway and KRAS signal gene sets were significantly upregulated in breast cancer resistant cells (Supplementary material 1). This suggests that the chemoresistance of breast cancer cells is closely related to the regulation of the cell cycle and apoptosis.
We also introduced clinical data and transcriptome profiles of breast cancer patients undergoing NAC using the GSE25066 database. According to the GGI evaluation method, GSE25066 samples were divided into two groups: GGI-low (NAC resistance, n=157) and GGI-high (NAC sensitivity, n=335) (Fig. 1D). Principal component analysis of gene expression showed that there was a significant difference in gene expression between the two groups (Fig. 1E). There were three hundred and forty-seven DEGs between the two groups. Thirty-six genes were commonly overexpressed or downregulated in both NAC chemo-resistant tumors and cells compared to the sensitive tumors and WT cells (Fig. 1F and G).
A Risk Model with NAC response was Constructed Based on GGI Lever
To further explore the clinical value of the differentially expressed genes, we performed the LASSO regression model based on the expression and prognosis data of 246 breast cancer patients who received NAC from the GSE25066 training set (Fig. 2A and B). We obtained two gene sets: 1se (containing 12 genes) and min (containing 18 genes). The ROC analysis showed that the 1se set (AUC=0.97) and min set (AUC=0.98) of characteristic genes both have good diagnostic values for evaluating the resistance of breast cancer to neoadjuvant chemotherapy (Fig. 2C). Considering the cost of patient detection, we chose the 1se set: HJURP, IFI27, RAD51AP1, EZH2, DNMT3B, SLC7A5, DBF4, USP18, ELOVL5, PTGER3, KIAA1324, and CYBRD1 (all 12 genes mentioned below refer to these 12 genes). The results of the validation set analysis were consistent with those of the training set. (Supplementary material 2). The complete names of the 12 genes and their main functions are listed in Table 2. The mathematical model of the risk score was as follows: Risk score =h0(t)×exp (HJURP0.0887787219+IFI27×0.0047893433+ RAD51AP1×0.0554382065+EZH2×0.0351434032+DNMT3B×0.0301466701+SLC7A5×0.0511106486+DBF4×0.0100252843-USP18×0.0251921906-ELOVL5×0.0158776302-PTGER3×0.0116989697-KIAA1324×0.0002110923- CYBRD1×0.0062603332). The cutoff risk score was 0.51 (Fig. 2D). Univariate and multivariate analyses also indicated that the signature-based risk factor score was a good prognostic indicator (Fig. 2E and F).
Characteristics of Tumor Cells and Tumor Microenvironment in High-risk Patients
To explore the internal and external factors of breast cancer cell NAC resistance, we analyzed the gene expression characteristics and tumor microenvironment of NAC-resistant cancer cells. According to the GSE25066 gene expression profile, the positive genes of risk score (HJURP, IFI27, RAD51AP1, EZH2, DNMT3B, SLC7A5, DBF4, USP18, and ELOVL5) were negatively correlated with the expression of negative genes (PTGER3, KIAA1324, and CYBRD1) (Fig. 3A). Utilizing the HitPredict database, we searched for interaction partners of their protein products, which were displayed in the protein interaction network diagram (Fig. 3B). We then performed GO and KEGG enrichment analyses of these genes. GO enrichment analysis indicated that 12 genes promote breast cancer resistance to NAC by upregulating DNA repair and metabolism-related pathways and downregulating some membrane receptor signaling mechanisms (Fig. 3C and Table 3). Furthermore, KEGG pathway enrichment analysis also revealed that the G protein-coupled receptor signaling pathway and transition metal ion transport were enriched by these genes (Table 3). The chord plot showed that 12 genes had complex interactions with the enriched GO pathway (Fig. 3D). For example, EZH2 and DNMT3B are involved in transferase activity (GO0016740) and negative regulation of gene expression (GO0045892 and GO0045814). DNMT3B, DBF4, and CYBRD1 jointly regulate metal ion binding (GO0046872). The above results suggest that these 12 genes cause drug-resistant cells in dormant cells to resist the killing effect of chemotherapeutic drugs by increasing DNA repair, reducing cell macromolecule synthesis, and cell metabolism. We also analyzed the correlation between the risk score and some known molecules associated with chemotherapy resistance in breast cancer. The risk score was significantly positively correlated with genes related to cell cycle regulation, DNA repair, transport, and efflux (Fig. 3E).
Tumor-infiltrating lymphocytes are critical components of the tumor microenvironment and are important external factors of chemotherapy resistance [29,30,31]. Therefore, we analyzed the difference in the proportion of tumor-infiltrating lymphocytes between the high-risk and low-risk groups. Compared with the low-risk group, the high-risk group recruited more CD4, CD8, and NK cells (Fig. 3F and G). Interestingly, we found that there were more tumor-infiltrating cells in the high-risk group, but the survival rate was significantly lower than that in the low-risk group (Fig. 2). We further explored the correlation between the risk score and immunosuppressive molecules in the GSE25066 dataset. The risk score was positively correlated with many immunosuppressive molecules, such as CTLA4, LAG3, ICOS, IDO1, and ADORA2A (Fig. 3H). The above results indicate that a large number of tumor-infiltrating cells in the high-risk group are depleted, which may lead to the failure of NAC for breast cancer and the lower survival rate of patients.
Evaluation of Risk Score on Prognosis of Overall Breast Cancer
To confirm whether the risk score for evaluating NAC resistance can predict the prognosis of breast cancer, we conducted the following study. According to TCGA BRCA data, the expression of the positive genes of risk score (HJURP, IFI27, RAD51AP1, EZH2, DNMT3B, SLC7A5, DBF4, USP18, and ELOVL5) in tumors was higher than that in normal and tumor-adjacent tissues. The expression of the negative genes (PTGER3, KIAA1324, and CYBRD1) in tumors was lower than that in normal and tumor-adjacent tissues (Supplementary materiale 3). We also found that high expression of all the positive genes and low expression of all negative genes were associated with poor prognosis of breast cancer (Supplementary material 3). The above results suggest that the risk score composed of these 12 genes may also serve as a prognostic factor for breast cancer. The expression heatmap showed that RAD51AP1, EZH2, DNMT3B, SLC7A5, DBF4, and USP18 were highly expressed in the high-risk group, whereas ELOVL5, PTGER3, KIAA1324, and CYBRD1 were expressed at low levels in the high-risk group (Fig. 4A). Survival status and risk score showed the distribution of patients into different risk groups (Fig. 4B). Kaplan–Meier survival analysis showed that patients in the high-risk group had a worse prognosis than those in the low-risk group (Fig. 4C). To facilitate clinical application, we constructed a nomogram based on data from the GSE25066 database in order to predict patient survival probability by weighing age, ER, PR, HER2, grade, stage, lymph nodes, tumor, and the signature-based risk score (Fig. 4D and Supplementary material 3). The above results suggest that the novel signature of the 12 genes can not only evaluate NAC response but also predict prognosis in breast cancer.
Compound Screening for Reversing Breast Cancer Resistance
We further explored other potential drugs that could be used for the treatment of patients in the high-risk group. We compared the transcriptome data of samples from the high- and low-risk groups. Under the conditions of FoldChange >1, 169 genes were upregulated, and 92 genes were downregulated in the high-risk group in GSE25066. Differentially expressed genes were imported into the CMap database and the top 30 potential drugs that could be used to treat high-risk patients were screened (Fig. 5A). The top three compounds are bambuterol (bronchodilators) [32], pravastatin (lipid-lowering agents) [33], and isocarboxazid (antidepressant, a non-selective and irreversible inhibitor of monoamine oxidase) [34]. It is noteworthy that there were five anticancer drugs among the candidate compounds. They are imexon (alkylating agents) [35], temozolomide (alkylating agents) [36], axitinib (inhibition of tumor growth and phosphorylation of VEGFR-2) [37], semaxanib (VEGFR (Flk-1/KDR) inhibitor) [38] and crizotinib (ATP competitive multi-target protein kinase inhibitor that inhibits met/ALK/ROS) [39]. Furthermore, we used the previously induced chemo-resistant cell line MCF7/EPI to verify the function of the selected candidate drugs to reverse chemoresistance. Here, we selected the top three compounds (bambuterol, pravastatin, and isocarboxazid) and five antitumor drugs (imexon, temozolomide, axitinib, semaxanib, and crizotinib) as candidate compounds for verification. Pravastatin, isocarboxazid, imexon, axitinib, and crizotinib had significant cytotoxic effects on MCF7/EPI cells (Fig. 5B). We also observed the morphology of cells in each group after 72 h of drug treatment. Cells treated with pravastatin, isocarboxazid, temozolomide, crizotinib, and fiacitabine were swollen, had many protrusions, and tended to die (Fig. 5C). Finally, we detected the apoptosis level of cells in each group using a TUNEL apoptosis detection kit. This is consistent with the previous results that groups treated with pravastatin, isocarboxazid, imexon, axitinib, and crizotinib displayed numerous dead cells (Fig. 5D). This suggests that pravastatin, isocarboxazid, imexon, axitinib, and crizotinib could inhibit or kill chemo-resistant cells.