Predictive miRNAs Patterns in Blood of Breast Cancer Patients Demonstrating Resistance Towards Neoadjuvant Chemotherapy

Introduction

Breast cancer (BC) remains a common problem affecting the health of women worldwide, and unquestionably, a leading cause of cancer-related deaths.¹ The known methods of therapy for BC include chemotherapy, surgery, endocrine therapy, molecular targeted therapy, radiation, and neoadjuvant therapy.^2,3 Neoadjuvant chemotherapy (NAC) refers to preoperative chemotherapy and has become a pivotal adjuvant treatment method for locally advanced BC.⁴ If reasonably adopted, NAC can effectively reduce the clinical staging of BC patients and increase the success rate of surgery and survival rate.^5,6 Unfortunately, drug resistance is responsible for the failure of chemotherapy in Hormone receptor-positive (HR+) BC.⁷ It’s obvious that an effective method that can predict chemotherapy resistance before chemotherapy is of great significance for HR+ BC patients with drug resistance to avoid excessive chemotherapy and to determine the therapeutic strategy after surgery.

microRNAs (miRs/miRNAs) are endogenous single-strand RNAs with 19–25 nucleotides that negatively modulate the expression of their target genes.⁸ Former investigations have revealed the differential expression pattern of miRNAs in drug-resistant BC patients and its close association with drug resistance in BC.^9,10 Cun et al¹¹ found that miR-92a-3p is upregulated in sera of BC patients and correlated with tamoxifen resistance. Serum miRNAs might be a promising biomarker in the early diagnosis, treatment of sensitivity/drug resistance, and prognostic evaluation of different cancer types including BC.^12–14 Therefore, the present study aimed to establish a drug resistance prediction model based on the differentially expressed serum miRNAs in drug-resistant BC. Recently, next-generation sequencing and microarray techniques have been increasingly adopted to measure miRNA expression levels.^15,16 Next-generation sequencing can reveal differentially expressed miRNAs in exosomes of BC patients.¹⁷ Moreover, researchers have adopted the nomogram model to predict breast pathologic complete response of NAC in BC, which helps clinicians to formulate the most beneficial individualized treatment plan for BC patients.¹⁸ However, the construction of the NAC resistance prediction model to screen out serum miRNA biomarkers in BC has not been described. By the use of the high-throughput sequencing technology to screen out differentially expressed serum miRNAs in BC patients sensitive or resistant to chemotherapy and verification of a large number of clinical samples, we analyzed the function of differentially expressed miRNAs in chemosensitivity in BC and predicted NAC resistance in BC patients and explored the clinical application value of serum miRNAs by constructing the nomogram prediction model of drug resistance based on analysis of pre-processed data, hoping to offer new ideas for prognostic evaluation of BC.

Materials and Methods

Ethics Statement

This study was a retrospective cohort study approved by the Ethics Committee of Cancer Hospital Affiliated to Xinjiang Medical University, approval number, XJZ-CR-2019-017. All participants of this study signed the informed consent before enrollment. The research is Declaration of Helsinki.

Study Subjects

The serum samples were collected from 110 chemosensitive and chemoresistant BC patients treated in the Affiliated Tumor Hospital of Xinjiang Medical University from 2017 to 2022, all of whom were diagnosed with stage IIA to IIIC locally advanced BC and received NAC before surgery, consisting of 53 HER-2 (human epidermal growth factor receptor 2) negative patients who received 6 cycles of NAC by using TAC regimen and 57 HER-2 positive patients who received 6 cycles of TcbPH regimen and the targeted therapy. Hormone receptor-positive patients, including those with positive estrogen receptor (ER) and progesterone receptor (PR), received adjuvant endocrine therapy. The inclusion criteria were as follows: a) diagnosed as BC by pathology; b) complete and clear pathological data; c) receiving NAC. The exclusion criteria were as follows: a) inflammatory BC; b) distant metastasis; c) during pregnancy or lactation period; d) chemotherapy intolerance or reneging; e) with other serious diseases. The treatment efficacy was assessed by the response evaluation criteria in solid tumors:¹⁹ complete remission (CR), no tumor found by clinical means; partial remission (PR), reduction of breast mass > 50%; stable disease (SD), reduction of breast mass < 50% or enlargement < 25%; progressive disease (PD), enlargement of breast mass > 25%. Patients with CR and PR were considered chemotherapy-sensitive (N = 68), and patients with SD and PD were defined as chemotherapy-resistant (N = 42).

Data Collection

The general information and clinical data of study subjects were collected using patient information forms. General information included age, gender, nation, menstrual status, contact information and address. Clinical data consisted of tumor size, lymph node metastasis, number of lymph nodes, histological grade, molecular classification, vascular invasion, nerve invasion, ER, ER ratio, PR, PR ratio, HER-2, and Ki-67 expression. The BC tissue samples were collected by core needle puncture. According to the pathological diagnostic criteria of BC and ER and PR based nuclear staining, ER and PR negativity (-) were defined at the staining score of 0 or 1, and ER and PR positivity (+) were defined at the staining scores of 2 or 3. Through immunohistochemistry, as an immunohistochemical score of 3+, and HER-2 negativity (-) was defined at the immunohistochemical score of 0 or 1, and HER-2 positivity (+) was defined at the immunohistochemical scores of 2 or 3.

Collection and Treatment of Serum Samples

Firstly, 5 mL of fasting blood sample was collected from each participant by venipuncture before NAC treatment. The blood samples were then centrifuged at 1600 g for 10 min at 4 C. The upper serum was divided into 2 tubes and kept at −80 C. All these procedures were completed within 2 h of blood collection.

High-Throughput RNA Sequencing

Total RNA was extracted from sera following the protocols of TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The RNA concentration and purity were detected using Nanodrop 2000 (Thermo Scientific Pierce, Waltham, MA, USA). The RNA integrity was detected by agarose gel electrophoresis and RNA electrophoresis. The extracted RNA satisfied with the following parameters: total concentration ≥ 5 μg, concentration ≥ 200 ng/μL, 1.8 ≤ OD260/280 ≤ 2.0. The paired-end sequencing library was constructed according to TruSeqTM strand. RNA sequencing was carried out on Hiseq4000 platform (Illumina, San Diego, CA, USA) by Xinjiang Ouyi Biotechnology (Xinjiang, China).

Bioinformatics Analysis

The Illumina 5’ and 3’ sequencing adapters were filtered using Cutadapt software (https://cutadapt.readthedocs.io/en/stable/guide.html#), and rRNA, tRNA, miscRNA, snRNA and snoRNA were removed using BLASTN software. The filtered sequences were searched again in the miRBase 21.0 database (http://mirbase.org/) to identify new miRNAs. The expression level of miRNA was calculated using DEGseq. The differentially expressed miRNAs were screened out using p < 0.05 and FC (fold change) > 2. The target genes of differentially expressed miRNAs were predicted using miRNA target software miRDB (http://www.mirdb.org/) and miRWALK (http://mirwalk.umm.uni-heidelberg.de/). GO and KEGG function classification of target genes regulated by differentially expressed miRNAs was conducted using Blast2GO software. The miRNA-target gene co-expression network was constructed using Cytoscape 3.6.1 (https://cytoscape.org/). Functional enrichment analysis of the 5 candidate serum miRNAs was performed using Reactome (https://reactome.org/).

RT-qPCR

Total RNA was extracted from serum samples using a TRIzol kit (ET111, TransGen Biotech, Beijing, China) and reversely transcribed into cDNA using a reverse transcription kit (QP013, GeneCopoeia, Rockville, MD, USA). RT-qPCR was run using ABI QuantStudioTM6 Flex (ABI, Foster City, CA, USA) under conditions of 5 s at 95 C and 40 cycles of 30s at 60 C and 2 min at 72 C. The primers were synthesized by Xinjiang Ouyi Biotechnology. The relative mRNA expression was estimated using the 2^−ΔΔCt method. U6 was utilized as a reference gene. Table 1 shows the primer sequences.

Table 1 RT-qPCR Primers Used in This Study

Statistical Analysis

Statistical analysis was carried out using R software version 4.0.1^20,21 (R development core team; http://www.r-project.org) with the caret, rms, pROC packages. The p value of 0.05 was set as a significance level, and all tests were two-sided. We distributed 70% BC patients into the development group and 30% into the validation group at random to construct and validate the nomogram model. Data characteristics between groups were described and compared using chi-square test or t-test. The variables to develop the nomogram model were selected following the Akaike information criterion (AIC) using the stepwise selection method. The nomogram model validation consists of 3 parts: identification, calibration and clinical application. The discriminative ability of the nomogram model was assessed using the area under the curve (AUC) of ROC. The closer AUC to 1, the better the discriminative ability of the prediction model. The calibration of the prediction model refers to the consistency between the predicted result of the nomogram model and the observation probability, which was evaluated using the calibration curve. The decision curve analysis (DCA) was used to evaluate the clinical application value of the nomogram model.

Results

Comparison of the Baseline Characteristics of BC Patients

Table 2 shows the general information and clinical data of all the participants receiving NAC treatment (N = 110). The development group and validation group showed similar distribution of pathological parameters and no significant differences in age, nation, hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p, hsa-miR-8485, menopausal status, tumor diameter, lymph node metastasis, number of lymph nodes, histological grade, molecular classification, vascular invasion, nerve invasion, ER, ER ratio, PR, PR ratio, HER-2, and Ki-67 (all p > 0.05).

Table 2 Analysis of Clinical Characteristics and miRNA

High-Throughput Sequencing Result of miRNAs in Chemosensitive and Chemoresistant BC Patients

To begin with, high-throughput sequencing (Illumina) of 4 serum samples from BC chemotherapy-sensitive patients and 3 serum samples from BC chemotherapy-resistant patients was conducted, with the presence of significant differences between groups (p < 0.05) and FC > 2 as screening standards. The differentially expressed miRNAs were screened out by bioinformatics analysis. Compared with chemotherapy-sensitive patients, chemotherapy-resistant patients had 44 differentially expressed genes in the sera, including 20 upregulated miRNAs and 24 downregulated miRNAs (all p < 0.05, Table 3, Figure 1).

Table 3 Analysis of Expression Levels of miRNAs in Human Sera ()

Figure 1 Analysis of differentially expressed miRNAs. (A) Volcano map of differentially expressed miRNAs, each dot representing one miRNA. Red dots represent upregulated miRNAs, blue dots represent downregulated miRNAs, and gray dots indicate no significant change. (B) Heat map of differentially expressed miRNA clusters. The abscissa shows the sample name and the ordinate shows the miRNA name. Red indicates upregulated expression and blue indicates downregulated expression. DR1-DR3 refer to chemoresistant serum sample, and S1-S4 refer to chemosensitive serum sample.

Prediction and Enrichment Analysis of miRNA Target Genes

To further investigate the regulatory network of miRNAs, we analyzed the downstream target genes of 44 differentially expressed miRNAs on miRWALK and miRDB databases. Totally 869 target genes were predicted for the 44 differentially expressed miRNAs, among which 806 target genes were predicted on miRWALK database, and 75 on miRDB database, with 12 shared target genes (Figure 2A).

Figure 2 Prediction and enrichment analysis of miRNA target genes. (A) Venn Diagram of target genes regulated by 44 differentially expressed miRNAs analyzed on miRWALK and miRDB databases; (B) GO function enrichment analysis of target genes regulated by differentially expressed miRNAs using Blast2GO software; (C) KEGG signaling pathway enrichment analysis of target genes regulated by differentially expressed miRNAs.

To understand the function of miRNAs, we conducted GO and KEGG enrichment analyses. GO enrichment analysis manifested target genes mainly participate in the positive regulation of type I interferon, amino acid transporter activity, and regulatory process of synaptic vesicle secretion (Figure 2B). KEGG analysis revealed that target genes are primarily engaged in the regulation of the thyroid hormone, MAPK, PI3K/Akt, and Rap1 signaling pathways (Figure 2C). The target genes of the differentially expressed miRNAs might be involved in the regulation of drug resistance in BC through the modulation of these functions and signaling pathways.

Expression Levels of Serum miRNAs in Chemoresistant and Chemosensitive BC Patients

To screen out miRNAs that are associated with chemotherapy resistance in BC, we identified 5 miRNAs in association with drug resistance in BC or BC progression among 40 differentially expressed miRNAs by literature review, including hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p and hsa-miR-8485. RT-qPCR revealed that hsa-miR-4741, hsa-miR-6730-3p, and hsa-miR-6831-5p were differentially expressed between chemosensitive and chemoresistant BC patients (Table 4, all p < 0.05). Hsa-miR-4741, hsa-miR-6730-3p, and hsa-miR-6831-5p might play pivotal roles in the progression and drug resistance in BC. Thereby, the 5-iRNA features might become the candidate diagnostic factors in chemoresistance in BC patients.

Table 4 General Information and Clinical Data of BC Patients

Multivariate Logistics Regression Analysis and Establishment of Nomogram Model

Hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p, hsa-miR-8485, nation, ER ratio, and PR ratio were subjected multivariate logistic regression analysis (Table 5). The p value was relaxed to p < 0.15. Hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p, hsa-miR-8485, nation, ER ratio, and PR ratio were included in the nomogram prediction model for NAC resistance in BC using software R4.0.1 (Figure 3A). The Points and Total points of these factors were assessed through nomogram model to predict whether NAC resistance occurs. The consistency, predictive effect, and clinical application value of the nomogram model were examined using calibration curve, ROC curve, and decision curve. The calibration curve of the development group and validation group showed that the predicted calibration curve was closer to the standard curve (45 C straight line), indicating that the nomogram model had good consistency (Figure 3B). The AUCs for ROC curves in the development group and validation group were 0.884 and 0.750, respectively, suggesting that the nomogram model had good prediction ability (Figure 3C). The decision curve of the development group demonstrated higher net benefit of the nomogram model in predicting NAC resistance in BC than other models, and the net benefit of the nomogram model in predicting NAC resistance in BC was higher than other models in the validation group when the threshold value > 0.1, which elicited that the established nomogram prediction model had good clinical application value (Figure 3D).

Table 5 Multivariate Logistic Regression Analysis of the Development Group

Figure 3 Establishment of nomogram model for the diagnosis of drug resistance in BC. (A) Nomogram model predicting drug resistance in BC patients; (B) calibration curve of the nomogram model in the development group (a) and validation group (b), with the X axis representing predicted probability and the Y axis representing the actual probability of drug resistance in BC patients through the nomogram model. The standard curve was a straight line with slope 1 through the origin of the axes. The closer the calibration curve to the standard curve, the better the consistency of the model; (C) ROC curve of the nomogram model in the development group (a) and validation group (b); (D) decision curve of the nomogram model in the development group (a) and validation group (b), with the Y axis representing net benefit and the X axis the threshold probability. The black solid line represents no intervention for all patients, the gray solid line represents intervention for all patients, and the dotted line represents the nomogram prediction model.

Identification of Gene Sets Enriched with 5-miRNA Signature-Based Risk Score

To explore the possible mechanism of 5-miRNA in NAC resistance in BC, Reactome and KEGG signaling pathway analyses were carried out on 5-miRNA signatures in patients with high risk of drug resistance that was assessed through the nomogram model. Reactome function enrichment analysis disclosed that 5-miRNA signatures were enriched in inflammation, circulation, and metabolism, and bubble map showed the association between 5-miRNA signature with RHOJ GTPase cycle, cGMP effects, and signaling by GSK3beta (Figure 4A). Additionally, KEGG pathway analysis indicated that most target genes of the 5-miRNA were of paramount importance to the signaling pathways regulating pluripotency of stem cells, insulin, Ras, and ErbB pathways (Figure 4B). These results suggested that 5-miRNA signatures were most relevant to inflammation and cancer metastasis. Based on the KEGG analysis result, we hypothesized that immune escape promotes NAC resistance in BC, which might lead to higher drug resistance rates in patients with high-risk scores than patients with low-risk scores. Overall, the 5-miRNA signature was defined as the diagnostic factor of drug resistance in BC in the present paper.

Figure 4 Identification of gene sets enriched with 5-miRNA signature-based risk score. (A) Reactome function enrichment analysis; (B) KEGG pathway analysis.

TCGA Database Analysis

The expression of 5 miRNAs in breast cancer was analyzed by TCGA database. Only HSA-MIR-4741 and hsa-miR-6730 were found to be highly expressed in breast cancer (p = 2.416400e-04), the expression of hsa-miR-6730 was significantly lower in breast cancer (P = 6.38199999999367e-05) (Figure 5).

Figure 5 TCGA database analysis. (A) Hsa-miR-4741 is highly expressed in breast cancer; (B) Hsa-miR-6730 is low expressed in breast cancer.

Discussion

Chemoresistance is a major contributor to the failure of anti-BC chemotherapy.²² In clinical work, patients with luminal-type breast cancer who are inconclusive in response to chemotherapy are less likely to respond to preoperative neoadjuvant therapy for a subset of patients whose tumors are 2–5 cm in diameter with or without lymph node metastasis, potentially contributing to disease progression; Whether neoadjuvant chemotherapy or surgical treatment should be given first is a difficult problem for clinicians. For some patients with HR (+) advanced breast cancer who can be treated with chemotherapy or endocrine therapy to reduce tumor burden more rapidly, there is an urgent need for accurate prediction models of chemotherapy resistance risk to assist clinical decision-making. It’s a known fact that miRNA expression levels can inform the response to chemotherapy in BC.²³ Hence, this study for the first time constructed a nomogram prediction model of NAC resistance in HR(+) BC based on multiple serum miRNAs, which could help clinicians predict whether NAC resistance occurs in HR(+) BC patients in a convenient and non-invasive manner, optimize clinical treatment plan to preclude delay in the implementation of other therapeutic approaches and condition advancement due to unnecessary preoperative NAC.

As indicated in a previous study, some miRNAs play key roles in regulating BC NAC resistance and have the potential to become NAC response predictors, prognosticators, and therapeutic targets.²⁴ First of all, we screened out 44 differentially expressed miRNAs in sera of the 3 selected chemoresistant BC patients by high-throughput sequencing and bioinformatics analysis, 20 of which were upregulated, and 24 were downregulated. Subsequently, we identified 12 shared target genes for these miRNAs. GO enrichment analysis revealed the engagement of target genes in the regulation of type I interferon, amino acid transporter, and synaptic vesicle secretion. KEGG analysis demonstrated that these target genes were enriched in the thyroid hormone pathway, MAPK pathway, PI3K/Akt pathway, and Rap1 pathway. Thereby, the target genes of the differentially expressed miRNAs might modulate chemoresistance in BC through above-mentioned functions and pathways.

Furthermore, we identified 5 miRNAs related to BC chemoresistance or development by literature review, including hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p and hsa-miR-8485. miR-3646 is a highly expressed miRNA in docetaxel-resistant BC cells.²⁵ Our results indicated that hsa-miR-4741, hsa-miR-6730-3p, and hsa-miR-6831-5p were differentially expressed between chemosensitive and chemoresistant BC patients. miR-4741 is implicated in the common miRNA signatures for colon and rectal cancer.²⁶ miR-6730-3p increases proliferation of melanoma cells.²⁷ miR-6831-5p is upregulated in bladder cancer.²⁸ We inferred that hsa-miR-4741, hsa-miR-6730-3p, and hsa-miR-6831-5p play essential roles in BC development and chemoresistance. We thereafter selected hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p and hsa-miR-8485 as study subjects.

The ethnocultural disparities might have a role to play in the screening, diagnosis, management and prognosis of BC patients.²⁹ The positive expression of ERα is in association with chemoresistance in BC.³⁰ PR expression is beneficial to the resistance of BC cells to gemcitabine.³¹ Then, we included hsa-miR-3646, hsa-miR-4741, hsa-miR-6730-3p, hsa-miR-6831-5p, hsa-miR-8485, nation, ER ratio and PR ratio for the multivariate logistic regression analysis and the development of the nomogram model. The nomogram model has potential predictive values and practical implications in BC.³² In our study, the nomogram model showed good consistency, predictive effect, and clinical application value.

Inflammation, circulation, and metabolism constitute a mechanism of chemoresistance in BC.³³ Reactome function enrichment analysis revealed that the signatures of above 5 miRNAs were enriched in inflammation, circulation, and metabolism. The regulation of Rho-GTPase can regulate cell shape in triple negative BC.³⁴ The production of cGMP is implicated in the regulation of chemosensitivity in BC cells.³⁵ GSK3β induces chemoresistance in BC by regulating Akt phosphorylation.³⁶ We found that the 5-miRNA signatures were related to RHOJ GTPase cycle, cGMP effects, and signaling by GSK3β. Most of the target genes of the 5 miRNAs play essential roles in signaling pathways regulating stem cell pluripotency, insulin, and Ras and ErbB pathways, which are important participants in inflammation and metastasis in cancers.^37–40 These results further suggested the association between 5-miRNA signatures and inflammation and metastasis. All in all, this study illuminated that the 5-miRNA signatures could be regarded as the diagnostic factor of drug resistance in BC.

The present study had two major limitations. First, the sample size we included was relatively small. Second, the downstream target genes and pathways of the 5 miRNAs and the deeper drug resistance mechanism in BC were not explored at the cellular level. In future investigations, we shall expand the sample size to further validate the clinical value and availability of the 5-miRNA drug resistance diagnostic factors and nomogram prediction model, and conduct cell experiments to further dig into the deep mechanism of NAC resistance in BC.