Next Article in Journal
Mechanobiology in Metabolic Dysfunction-Associated Steatotic Liver Disease and Obesity
Next Article in Special Issue
MicroRNAs Associated with IgLON Cell Adhesion Molecule Expression
Previous Article in Journal
MicroRNAs as Regulators of Radiation-Induced Oxidative Stress
Previous Article in Special Issue
Mechanism of DAPK1 for Regulating Cancer Stem Cells in Thyroid Cancer
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
CIMB
Volume 46
Issue 7
10.3390/cimb46070424
cimb-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Secondary Analysis of Human Bulk RNA-Seq Dataset Suggests Potential Mechanisms for Letrozole Resistance in Estrogen-Positive (ER+) Breast Cancer

by
Lincoln Sutherland
1,
Jacob Lang
1,
Norberto Gonzalez-Juarbe
2,3 and
Brett E. Pickett
1,*
1
Department of Microbiology and Molecular Biology, Brigham Young University, Provo, UT 84602, USA
2
J. Craig Venter Institute, Rockville, MD 20850, USA
3
Department of Cell Biology and Molecular Genetics, University of Maryland, College Park, MD 20742, USA
*
Author to whom correspondence should be addressed.
Curr. Issues Mol. Biol. 2024, 46(7), 7114-7133; https://doi.org/10.3390/cimb46070424
Submission received: 9 June 2024 / Revised: 26 June 2024 / Accepted: 3 July 2024 / Published: 6 July 2024
(This article belongs to the Special Issue The Significance of Transcription Factors, miRNAs, and lncRNAs in Anticancer Drug Development)
Browse Figures
Versions Notes

Abstract

:
Estrogen receptor-positive (ER+) breast cancer is common among postmenopausal women and is frequently treated with Letrozole, which inhibits aromatase from synthesizing estrogen from androgens. Decreased estrogen slows the growth of tumors and can be an effective treatment. The increase in Letrozole resistance poses a unique problem for patients. To better understand the underlying molecular mechanism(s) of Letrozole resistance, we reanalyzed transcriptomic data by comparing individuals who responded to Letrozole therapy (responders) to those who were resistant to treatment (non-responders). We identified SOX11 and S100A9 as two significant differentially expressed genes (DEGs) between these patient cohorts, with “PLK1 signaling events” being the most significant signaling pathway. We also identified PRDX4 and E2F8 gene products as being the top mechanistic transcriptional markers for ER+ treatment resistance. Many of the significant DEGs that we identified play a known role in ER+ breast cancer or other types of cancer, which partially validate our results. Several of the gene products we identified are novel in the context of ER+ breast cancer. Many of the genes that we identified warrant further research to elucidate the more specific molecular mechanisms of Letrozole resistance in this patient population and could potentially be used as prognostic markers with further wet lab validation. We anticipate that these findings could contribute to improved detection and therapeutic outcomes in aromatase-resistant ER+ breast cancer patients.

1. Introduction

Breast cancer is the second leading cause of cancer death in women [1]. Approximately 13% of women will develop invasive breast cancer, with 3% of women dying from it [1]. Roughly 83% of invasive breast cancer cases are diagnosed in post-menopausal women (50 years or older), with 91% of the deaths occurring in the same age group [1,2]. Given that those who are 65 years of age or older are the largest growing subpopulation in the United States (16.8% of the total population in 2020) [3], the overall number of new breast cancer cases is expected to rise accordingly. The American Cancer Society predicts that over 310,000 women in the United States will be diagnosed with breast cancer in 2024 [4], which is expected to rise to ~440,000 by 2030 [5].
Breast cancer can be classified as a heterogenous disease. Among the various types of breast cancer, some types have an affinity for particular hormones that contribute to rapid tumor growth. Three important markers that are generally used to classify breast cancer involve the positive or negative presence of the estrogen receptor (ER), progesterone receptor (PR), and/or the human epidermal growth factor receptor 2 (HER2) protein on the cell surface [6]. The estrogen and progesterone receptors both belong to the steroid receptor subgroup of ligand-activated transcription factors within the nuclear receptor superfamily [7]. Their presence is a strong predictor of the efficacy of adjuvant and therapeutic hormonal therapy [7]. Amplification of the ERBB2 gene in tumor cells leads to increased presence of the protein at the cellular membrane, leading to potent proliferative and anti-cell death signaling [8,9]. Taken together, the ER, PR, and HER2 play a central role in the detection, diagnosis, and therapy of breast cancer [10,11,12,13].
Given that estrogen contributes to the promotion and progression of breast cancer, some therapeutic strategies specifically target estrogen synthesis via the ER and its associated intracellular signaling pathways [14]. The aromatase enzyme then catalyzes the final conversion of androgens into estrogens [14,15]. Thus, aromatase was identified as a potential therapeutic target in the late 1970s with the synthesis of aminoglutethimide as the first aromatase inhibitor (AI) [16]. The low specificity and selectivity of this early therapeutic led to the synthesis of subsequent aromatase inhibitors, with the most recent being Letrozole [16]. The increased affinity of Letrozole to the catalytic site of the aromatase enzyme competitively inhibits the binding of the precursor reactant to the aromatase enzyme [14]. Moreover, Letrozole has greater potency than other Ais, including anastrozole, exemestane, formestane, and aminoglutethimide [14]. Seventy-seven percent of all breast cancers are ER-positive (ER+); thus, treatment with aromatase inhibitors is currently the standard treatment for all postmenopausal women [7]. Although Letrozole is one of the most commonly used treatments for ER+ breast cancer, it also has limitations [17,18]. A prior study showed that the overall response rate (ORR) among a large multi-national study of postmenopausal women was 55% [17]. The ORR quantitatively represents the percentage of patients that achieve a favorable response, either total (complete disappearance of cancerous tissue) or partial (defined as a decrease in tumor size by at least 30%) [19]. Response rates are determined by either a mammogram or an ultrasound [17]. To our knowledge, the ability of existing prognostic methods to accurately identify treatment resistance prior to treatment is lacking.
A prior study examining the intracellular transcriptional response to Letrozole treatment used RNA-sequencing data to identify Proline-Rich Protein 11 (PRR11) as the only significantly overexpressed gene among 51 genes in chromosome arm 17q23. Focusing their analysis and validation work on this single PRR11-containing locus in treatment-sensitive vs. treatment-resistant ER+ breast cancer patients was useful [20], and amplification of this locus has been shown to correlate with poor clinical outcomes in breast cancer [21].
The aim of the current study is to perform a secondary analysis of the complete transcriptome from the same original RNA-sequencing study in treatment-sensitive vs. treatment-resistant ER+ breast cancer patients. We expect that re-analyzing this valuable dataset will identify genes, signaling pathways, and transcriptional mechanisms across the whole human transcriptome that are significantly associated with resistance to letrozole treatment in patients with ER+ cancer.

2. Materials and Methods

2.1. Retrieving Fastq Files, Preprocessing, and Enriching RNA-Sequencing Data

The sratools software (version 2.10.8) package was used to download the RNA-sequencing fastq files from the desired study (Gene Expression identifier: GSE145325) from the Sequence Read Archive (SRA) database at NCBI [20,22]. The fastq files were then pre-processed and analyzed using the snakemake-based Automated Reproducible Modular Workflow for Preprocessing and Differential Analysis of RNA-seq Data (ARMOR; version 1.5.7) [23,24]. This preprocessing workflow applies the following methods using established algorithms: (1) quality control with fastqc [25] using default parameters; (2) trimming with TrimGalore! [26] using default parameters, including a minimum quality score of 20, a phred scale of 33, and a minimum length of 20 bases; (3) pseudomapping and quantification to the human transcriptome (GRCh38 release 98) using Salmon (version 1.3.0) [27] with a k-mer setting of 31, a read length of 63, and the following flags: seqBias,gcBias, fldMean=250, and fldSD=25; (4) read count normalization and differential gene expression using edgeR (version 3.36.0) [27,28] with default parameters; and (5) functional enrichment using CAMERA (version 3.50.3) [29], with default parameter enriching against the H, C2, and C5 Molecular Signature database gene sets. Although the trimming and/or differential expression algorithms may change, each of these steps is commonly applied to preprocessing RNA-sequencing data.
The differential gene expression results were then used as an input to the Signaling Pathway Impact Analysis (SPIA; version 2.50.0) algorithm [30], with default parameters other than using 2000 replicates and a Bonferroni-adjusted p-value cutoff of 0.05. This algorithm uses bootstrapping and a null distribution to calculate significantly perturbed pathways that are enriched in statistically significant DEGs, regardless of the fold-change value or direction [31,32,33]. The SPIA method was specifically chosen because of its ability to take into account the topology of gene products that contribute to a pathway rather than simply calculating whether there is a significant overlap between the significant genes and the pathway, which decreases the number of false-positive results. Another benefit is that this algorithm uses bootstrapping and permutation to calculate a custom null distribution for each pathway. This null distribution is then used to determine a specific and unique p-value for the differentially expressed genes in the dataset for each pathway. This approach generally produces fewer false-positive results when compared to other hypergeometric-based methods and enables it to predict whether any significant pathway is activated or inhibited without needing to perform separate analyses of upregulated and downregulated DEGs [30,31,32,33].

2.2. Target and Mechanistic Marker Prediction

The Pathway2Targets algorithm (version 2.2) was then used to predict and prioritize relevant therapeutic targets from the transcriptomics signatures that could be repurposed in the case of Letrozole resistance, similar to prior work [34,35,36,37]. Briefly, this algorithm identifies known drug targets within each of the significant signaling pathways. It then uses multiple criteria from the OpenTargets database [38], as well as the signaling pathway information, to rank the relevant therapeutics and their targets according to a customizable weighting scheme. The justification for this approach is that a gene product does not need to be differentially expressed to be a relevant therapeutic target that is capable of reducing the signs and/or symptoms of a given condition, particularly if the target is upstream of one or more significant DEGs.
Separately, the Salmon read counts for each sample were combined into a tabular format and labeled as “resistant” or “sensitive” to treatment with Letrozole. This table was then used as an input to the XGboost algorithm [39], which implements a tree-based method to train a model from 80% of the dataset and then quantifies its performance using the remaining 20% of the data [40,41], which minimizes model overfit. For the initial analysis, the gain metric was calculated from the read counts for all detected genes across all samples. The results were sorted in descending order, such that the genes/features with the highest gain value were listed first. Given the gain metrics from the whole transcriptome, the number of genes/features being evaluated was reduced to the best two mechanistic markers from the original analysis since this is a number that is easily accommodated by qRT-PCR (or similar) molecular methods. This approach has been successfully applied previously with acceptable performance and accuracy [40,41,42]. The XGboost algorithm was selected since prior work has shown that tree-based classifiers are faster and more accurate than other machine learning-based methods such as support vector machines, neural networks, and Bayesian approaches [43]. The XGboost parameters that were changed from the default included: constructing 10,000 parallel trees to ensure appropriate coverage of the available tree space, a subsample of 0.5, using a binary–logistic objective, and an area under the curve (AUC) evaluation metric.

2.3. Protein Network Analysis

The public protein–protein interactions in the STRING database were used to construct a PPI network for all significant DEGs, the top five targets, and the top five mechanistic transcriptional markers [44]. Cytoscape was then used to visualize the network, while the 12 algorithms incorporated into the cytoHubba app were used to compute multiple node-based metrics from the Cytoscape network [45,46]. These algorithms were specifically selected to predict key nodes from a network while also minimizing bias due to either high-degree or low-degree proteins [46].

3. Results

3.1. The Majority of Differentially Expressed Genes in ER+ Treatment Resistance Are Upregulated

We began by retrieving, preprocessing, and analyzing an existing dataset, which consisted of 58 ER+ stage II or stage III primary breast tumor samples from newly diagnosed female patients who had been treated with Letrozole for a median of 7.2 months prior to surgery. These patients were at least 60 years old when the original study started, had a mix of progesterone receptor presence, were primarily HER2-negative, and consisted mostly of ductal carcinoma. No mention of prior treatment with other therapeutics was associated with the original study.
The RNA sequencing data from these clinical samples are publicly available in the NCBI Gene Expression Omnibus (GEO) database. Our analysis detected a total of 18,735 genes between responders (treatment-sensitive patients) and non-responders (treatment-resistant patients), with 105 of those genes having significant differential expression (FDR p-value < 0.05; Table 1; Supplementary Table S1). This significance cutoff was used since the FDR method applies a multiple hypothesis correction to reduce false positives in the results when many statistical tests are performed.
We observed that among these 105 Differentially Expressed Genes (DEGs) (Figure 1), the top 20 (Table 2) included 17 protein-coding genes that consisted of two S100 calcium-binding proteins, various transcription factors, and matrix metallopeptidase 7 (Figure 1 and Table 2). Three of the top twenty DEGs coded for immunoglobulin variable regions, while one was a putatively transcribed unprocessed pseudogene, Ovostatin II. All but one of these 20 DEGs were upregulated in the group that failed to respond to Letrozole treatment. The sole exception to this was the downregulation of the gene that encodes the growth arrest and DNA-damage inducible gene 45 gamma expression (GADD45G), which is involved in the regulation of cell growth and apoptosis [47] and showed a log2 fold-change (log2FC) value of −1.62.
The original study identified PRR11 as one of the few statistically significant differentially expressed genes among the 51 genes in chromosome arm 17q23 that were included and validated in their analysis [20]. However, our analysis identified PRR11 to be the 1211th differentially expressed gene across the whole transcriptome, which did not surpass the statistical (FDR-adjusted p-value < 0.05) threshold for significance.

3.2. The Signaling Pathway Impact Analysis Identified Four Significantly Affected Pathways

We then wanted to determine which intracellular signaling pathways were significantly enriched with the detected DEGs using the Signaling Pathway Impact Analysis (SPIA) algorithm (Bonferroni-adjusted p-value < 0.05). Briefly, this robust pathway enrichment algorithm uses a bootstrap-based approach to generate a null distribution for each pathway and calculate a p-value. The results from our pathway enrichment analysis included four pathways that were predicted to differ between responders and non-responders to Letrozole treatment (Table 3). These pathways included those associated with activated “PLK1 signaling events and anti-tumoral activity” as well as the activated “FOXM1 proliferation-associated transcription factor network”.

3.3. Targets Prioritized for Repurposing from the Identified Pathways

We then wanted to use the pathway results to predict any existing therapeutic targets that could be repurposed to treat the resistance phenotype with small molecules, monoclonal antibodies, peptides, or other modalities. We consequently applied the Pathway2Targets algorithm to perform this target prioritization analysis, which predicted 60 therapeutic targets (Table 4 and Supplementary Table S2). Briefly, this algorithm calculates a weighted score for each target based on the sum of values assigned to 26 attributes. Notably, the predicted targets for the non-response (resistance) phenotype included Vascular endothelial growth factor A (VEGFA), a current target for solid tumors, as well as Estrogen Receptor 1 (ESR1), Nitric Oxide Synthase 2 (NOS2), and various matrix metalloproteinases (MMP9 and MMP2).

3.4. Machine Learning Predicted Two Robust Mechanistic Transcriptional Markers

We next wanted to predict a subset of mechanistic transcriptional markers that could potentially be used to better identify patients who will not respond to treatment with Letrozole and better understand their underlying resistance mechanism(s). To perform this analysis, we applied a decision tree-based machine learning approach to predict features (i.e., expressed genes) that most accurately classify Letrozole treatment responders vs. non-responders. We used the read counts for all samples as the input for this analysis to identify and rank 278 transcriptional biomarkers by the gain metric (Table 5). The performance metrics for all 278 transcriptional biomarkers yielded an area under the receiver-operator characteristic (AUROC) curve of 0.972 (97.2%; Figure 2), indicating an exceptional ability to classify Letrozole responders vs. non-responders based on global gene expression. The sensitivity for this complete set of mechanistic transcriptional markers was predicted to be 100%, with a predicted specificity of 94%.
Given the resources required to generate whole transcriptome data using RNA-sequencing, we repeated our machine learning analysis to predict the mechanistic transcriptional markers that are most associated with treatment nonresponse in this patient population. To achieve this, we used the top two mechanistic markers from our first analysis as an input for a subsequent analysis. This more focused second analysis identified PRDX4 and E2F8 transcripts as markers with potential as diagnostic targets, with an AUROC for only these two gene products of 0.854 and an overall accuracy of 88.2% (Figure 3).

3.5. Protein–Protein Interactions Reveal a Potential Treatment Resistance Network

We then constructed a protein–protein interaction (PPI) network to predict the protein(s) that potentially played a role in Letrozole resistance for these patients (Figure 4). To achieve this, we retrieved the known PPIs from the public STRING database for the 105 DEGs, as well as the top five drug targets and the top five transcriptional biomarkers. Interestingly, this analysis revealed a relatively well-connected local network, which indicates that the majority (87/115) of the input gene products directly interacted with each other and at least partially supports our earlier predictions. We applied all 12 available algorithms within the CytoHubba app to this network to identify the key proteins with minimal bias to any individual approach. We then calculated the rank determined by each algorithm and then averaged them together. This analysis identified the top 10 “key” proteins in the network as CDK1, CCNB1, CCNA2, BIRC5, AKT1, ESR1, MMP9, CDCA8, CCND1, and PLK1. We found the ESR1, PLK1, MMP9, and AKT1 proteins to be of particular interest since they were identified previously as the targets of Letrozole, a significant signaling pathway, and two therapeutic targets. This profile at least partially represents the intracellular transcriptional response to Letrozole therapy.

4. Discussion

The primary aim of the current study was to perform a secondary analysis of a publicly available RNAseq dataset produced from 58 samples from patients with ER+ breast cancer following Letrozole treatment. In addition to finding PRR11-related mechanisms, the original study reported a cluster of E2F4-regulated genes that were sensitive to subsequent CDK4/6 inhibitor (palbociclib) treatment in a separate cohort of ER+ breast cancer patients [48].
Our analysis identified 105 statistically significant DEGs, with the top 20 gene products consisting of 17 protein-coding genes, including three encoding immunoglobulin variable regions and one pseudogene. We also identified four significant signaling pathways, twenty therapeutic targets that could be repurposed to reduce Letrozole resistance and/or target resistance cells, and two potential mechanistic transcriptional markers that demonstrate high combined specificity and sensitivity in identifying Letrozole non-responders. Below, we describe our primary findings in the context of breast, and other cancer-related studies in an effort to use prior observations to at least partially validate our mechanistic findings.

4.1. Differentially Expressed Genes

In contrast to our findings, which identified 105 statistically significant DEGs, the original study that utilized the same data identified PRR11 as the only differentially expressed gene among the 51 genes in chromosome arm 17q23 [20]. We believe that the mostly likely contributor to this seeming discrepancy was the focus of the original study on the particular locus in 17q23 that contains 51 genes. In contrast, our analysis evaluated all of the detected mRNAs in the cells, which increased the scope and the statistical stringency needed to obtain significant results after incorporating multiple hypothesis corrections. PRR11 underwent substantial differential gene expression in our analysis, although our thresholds did not categorize it as significant when the whole transcriptome was taken into account.
Transcription of the SRY-Box Transcription Factor 11 (SOX11) gene (with a log2FC value of 3.87) was our most statistically significant finding. Briefly, SOX11 is part of the Sox gene family, comprising 20 transcription factors, which can further be divided into nine subgroups based on functional similarity [49]. SOX11 belongs to the SoxC group (composed of Sox4, Sox11, and Sox12), which are critical in the development of the nervous system and show substantial transactivation potential [50], while also being required for the generation/proliferation of immature neuron progenitors [51]. SOX11 has recently been shown to behave as an oncogene and is a critical regulator of basal-like and luminal B breast cancers and their metastasis to the brain [52,53,54].
The upregulated S100 calcium-binding Protein A8 and A9 (S100A8 and S100A9) genes encode a heterodimer of two calcium-binding proteins, referred to as calprotectin [55]. This protein is secreted by neutrophils and is an important proinflammatory mediator in both chronic and acute inflammation [55]. Calprotectin has recently been discovered to play an important role in tumorigenesis, cell proliferation, and resistance to traditional cancer therapies [56,57]. S100 proteins lack the signal peptide required for secretion via the classical endoplasmic reticulum/Golgi route, and its secretion mechanism is not well understood [56,58]. Recent discoveries implicate BRCA1, a known tumor suppressor, as essential for regulating the levels of S100A8 and S100A9 in the body [56,59]. These S100 gene products have also been shown to be secreted due to IL-22 and IL-17 signaling [60], have been associated with poor outcomes in ductal carcinomas [61], and have been affected after AI treatment in an animal model [62]. While S100 proteins have been shown to be host markers of tumor development and progression, their role as modulators of Letrozole resistance has not been characterized previously. Thus, we believe that further experiments will help to elucidate the underlying mechanism(s) of S100A8 and S100A9 proteins in treatment-resistant ER+ breast cancer.
The Immunoglobulin Lambda Variable 3–25 (IGLV3–25) gene codes for light chain variable regions involved with antigen recognition [63,64] and has previously been implicated in various cancers [63,65,66]. Our analysis also found Matrix Metalloproteinase 7 (MMP7) to be upregulated in this dataset. This gene product encodes a proteolytic enzyme that secretes zinc and calcium endopeptidases, which are normally involved in wound healing and bone growth [67]. It has also been implicated in cancer progression, proliferation, differentiation, and/or apoptosis [68,69,70]. Basal-like breast cancer has upregulated MMP7 expression, which has been linked to DKK1 knockdown, a known tumor suppressor in breast cancer [71]. Interestingly, zoledronate, a common drug used to treat bone damage and osteoporosis [72,73,74], has been shown to decrease the expression of MMP7 in breast cancer cell lines [75].
Engrailed Homeobox 1 (EN1) was another upregulated gene that we identified, which is a protein-coding gene that helps control development [76] and has been shown to affect bone metastases from breast cancer tumors [77]. It has been implicated in quadruple-negative breast cancer and basal-like breast cancer [78,79,80]. Our findings could potentially expand its scope of relevance to treatment-resistant ER+ breast cancer. We believe that future validation experiments should be performed to confirm these differential expression results.

4.2. Intracellular Signaling Pathways

Our analysis found four notable intracellular pathways that were significantly modified among non-responders. This analysis predicted the “PLK1 signaling events”, “Syndecan-1 mediated signaling events”, and “FOXM1 proliferation transcription factor network” pathways to be activated; the “HIF-1-alpha transcription factor network” was predicted to be inhibited. Although the number of pathways is relatively small compared to other algorithms, we believe that applying a Bonferroni p-value correction minimized the number of false results. It is possible that this approach could be complemented by other approaches in future studies.
PLK1-related signaling was our most significantly impacted intracellular signaling pathway (activated among non-responders vs. responders) and is primarily involved in cell mitosis. PLK1 contributes to many functions, including controlling the G2/M checkpoint, coordinating the cell cycle, and serving as a regulator of cell division in eukaryotes, as well as in both DNA replication and chromosome segregation [81,82]. In general, the overexpression of PLK1 has been associated with various cancer types, including ER+ breast cancer [82,83]. PLK1 inhibition has been shown to exhibit anti-tumoral activity by blocking mitosis in cells with higher PLK1 expression, making it an attractive target for cancer therapy [84,85]. PLK1 knockdown in hormone-independent ER+ breast cancer has been associated with decreased cell viability and sensitization to radiation treatment [86]. These past studies serve to partially validate our significant pathway findings and lend credence to ongoing research into PLK1 inhibition as a potential treatment for patients with ER+ breast cancer.
The pathway for Syndecan-1 (SDC1, CD138) signaling was activated in non-responders. SDC1 is a transmembrane heparan sulfate proteoglycan responsible for maintaining typical cell morphology and is found abundantly in epithelial cells [87]. The SDC1 pathway is generally involved in cell proliferation, cell–matrix interactions, growth factor signaling, and angiogenesis [88,89]. SDC1 overexpression tends to correlate with more aggressive tumors, and it has an inverse correlation with ER expression [90]. This inverse correlation between ER expression and SDC1 signaling could help explain why Letrozole was ineffective for at least a subset of non-responders and merits further experimentation. Our analysis also showed MMP-7 as a highly significant DEG, which has been shown to anchor itself to syndecans and results in the shedding of important membrane-bound ligands, such as EGF and TGF-a, leading to cell invasion [69]. Future follow-up experiments to validate these findings could include using CRISPR-i (or similar) systems and/or drugs to affect these pathways.

4.3. Target Prioritization and Repurposing

Our approach to predicting potential therapeutic targets by combining public data with analysis-specific pathway information and a customizable weighted scoring method is unique. While this approach has been useful in the past [32,33], other approaches exist that could complement our findings. Specifically, since the target data comprise more antagonists than agonists, this could potentially be a limitation. Vascular endothelial growth factor A (VEGFA) was the top potential target identified by our prioritization algorithm, which has already been approved for certain indications, including breast cancer. This growth factor belongs to a ligand family composed of six related proteins and plays an important role in breast cancer and other forms of cancer [91]. The primary role of this protein is to stimulate angiogenesis, which facilitates tumor growth [91,92]. Our computational prediction of VEGFA as a relevant target warrants future experimentation and serves as at least a partial validation of our approach.
Estrogen Receptor 1 (ESR1) encodes the estrogen receptor and was identified as a therapeutic target. Of note, ESR1 has been found to encode gain-of-function mutations that promote tumor metastasis and resistance to endocrine therapy [93,94,95] and can be detected with a liquid biopsy [94]. The detection of ESR1 as a therapeutic target is logical given that resistance to AI treatment enables the continued activation of signaling pathways that are modulated by ESR1.
Two additional highly ranked targets identified by our algorithm were Matrix Metallopeptidase 9 (MMP9) and Matrix Metallopeptidase 2 (MMP2), both members of the extracellular matrix metalloproteinase family (MMPs). MMPs assist in extracellular matrix remodeling through the activation of substrates via enzymatic cleavage [96]. They are involved in a variety of normal processes in the body, and their dysregulation plays a vital role in processes such as aging. Their dysregulation also plays a role in several abnormal conditions, including preeclampsia among pregnant women and breast cancer [96,97,98]. Higher levels of MMP9 have been correlated with higher tumor grades and resistance to endocrine therapy [96]. MMP9 was further characterized as a vital component of the metastatic cascade early in tumorigenesis that promotes colonization in the lungs [99]. Other studies have implicated MMP2 with tumor metastasis to other organs, including the brain [100,101]. Overall, both MMP9 and MMP2 have been identified as promising markers for predicting the prognosis of patients with breast cancer [102,103]. Interestingly, a phase III clinical trial targeting MMP proteins with marimastat sought to reduce tumor blood flow in patients with metastatic breast cancer [104]. Although MMP7 was identified as a significant DEG, it was not predicted as a useful target in this analysis. However, MMP7 is one of the several MMPs that are inhibited by doxycycline, which is a commonly used treatment for various forms of cancer.
Fibroblast growth receptor 3 (FGFR3) belongs to a family of four highly conserved, transmembrane receptor tyrosine kinases. Once activated, they initiate intracellular cascades that carry out a variety of functions, including promoting cellular proliferation and survival [105]. This gene has been reported in other studies to be significantly upregulated in breast cancer [106]. FGFR3 is a known target of microRNA-593-3p, and the overexpression of miR-593-3p down-regulates FGFR3 expression, which slows breast cancer progression [107]. Various inhibitors of FGFR3 are at various stages in phase I and phase II clinical trials for breast cancer.
Another potential target we identified was AKT Serine/Threonine Kinase 1 (AKT1). AKT1, AKT2, and AKT3 are all isoforms of protein kinase B (AKT), an essential member of the PI3K/AKT signaling pathway [108]. Both AKT1 and AKT2 have been implicated in breast cancer previously, with some mixed results regarding their specific function. However, it is generally accepted that AKT1 plays an important role in the early stages of tumor development (tumor initiation and proliferation), while AKT2 primarily assists in tumor metastasis [108,109]. The original study for this dataset identified PI3K as an important player in AI treatment [20], which is supported in other studies [110]. Regardless of its specific function, AKT1 has a well-documented correlation with breast cancer and is an appropriate target for pharmaceutical intervention. Additional work at the discovery, pre-clinical, and clinical stages is needed to determine whether such treatments could reduce the effect of Letrozole resistance in breast cancer patients.

4.4. Mechanistic Transcriptional Marker Analysis

Our DEG analysis identified expressed genes that significantly differed between the two groups of samples (resistant vs. sensitive). In contrast, the aim of our mechanistic marker analysis was to determine the expressed genes that most consistently differed between the two groups of samples. Although the input for both analyses was the read counts, we employed a tree-based machine learning model for the mechanistic marker analysis that does not incorporate traditional statistics. Consequently, we did not expect a complete overlap between the DEGs and our predicted biomarkers.
Peroxiredoxin-4 (PRDX4), which was an up-regulated gene among non-responders, was the top transcriptional mechanistic marker identified in our analysis. This gene product belongs to a family of six small antioxidant isozymes and is generally localized in the endoplasmic reticulum [111]. PRDX4 plays an essential role in maintaining localized redox homeostasis and is upregulated when cells are under oxidative stress [112]. PRDX4 has been described as having a tumor-promoting effect that is well-documented among various forms of cancer, including lung and renal cancer, leukemia, and glioblastoma [113,114]. Prior work has shown the overexpression of PRDX4 to be associated with poor overall survival rates, shorter relapse-free survival among breast cancer patients [115,116], and metastasis. PRDX4 has also been shown to be associated with more advanced breast cancer tumors [117], which could also contribute to treatment resistance. Interestingly, PRDX4 was also identified as a potential therapeutic target for prostate cancer, suggesting a possible shared mechanism for androgens in both prostate and breast tissues [118,119]. Thus, the role of PRDX4 as a potential mechanistic marker and potential target in treatment-resistant breast cancer deserves future exploration.
Early Region 2 Binding Factor (E2F8), an up-regulated DEG from our study, was another potential mechanistic marker identified by our machine learning analysis. E2F8 is a member of the E2F family, which assists in assembling the core transcriptional machinery, making it crucial for cell division [120]. E2F8 has also been shown to promote angiogenesis, which helps establish the tumor microenvironment and is positively correlated with tumor malignancies [121,122]. E2F8 has been associated with CDK4/6 inhibitor resistance [123], is overexpressed in patients with breast cancer, and has been significantly correlated with poor patient survival and cancer progression [124,125]. E2F8 has been suggested as a marker of breast relapse-free survival and distant metastasis-free survival in breast cancer patients [126]. Interestingly, a prior study in HER2- has identified a different E2F Transcription Factor (E2F4) as having a role in treatment resistance in ER+ breast cancer [127] and also in the original analysis of our same ER+ cohort [48]. As such, we believe that the role of E2F8 as a mechanistic marker, specifically in aromatase-resistant breast cancer, warrants further exploration.
We acknowledge inherent deficiencies in training and applying a machine learning approach to a relatively small dataset. Unfortunately, although the GEO database had other bulk RNA-sequencing datasets that characterized HER2+ samples, there were no other public ER+ Letrozole-treated clinical samples. Consequently, augmenting our model with additional data was not possible. Additional wet lab validation experiments using qRT-PCR or similar transcriptional assays will be required in one or more distinct patient populations to better quantify the accuracy and performance of our predicted mechanistic markers.

4.5. Potential Treatment Resistance Mechanisms

Unfortunately, the original dataset only included samples from responsive and non-responsive tumors. As such, it was not possible to directly and accurately analyze and assess the pre-treatment state of the tumor. Similarly, we believe it would not be justified to incorporate additional comparisons from other datasets that evaluated triple-negative or ER- breast cancers given the distinct patient populations, enrollment criteria, and other experimental variables. Even so, our protein network analysis did identify a set of gene products that may contribute to an increased understanding of the intracellular transcriptional response to Letrozole treatment.
We recognize that a subset of the genes we identified in our analysis of this dataset are opposite of what has been reported in prior work, specifically MMP9 and ESR1, as being upregulated in Letrozole-resistant breast cancer tumors, which agrees with our findings [128]. Another study found that estrogen is known to attract and activate neutrophils, which consequently upregulates the expression of both S100A8 and S100A9, in agreement with our findings [129]. Multiple immunoglobulin lambda genes, which were upregulated in our analysis, agree with prior work in humans and mice that showed a similar response in resistant samples [130,131]. SOX11 has been reported to be downregulated due to epigenetic effects in some breast cancer tumors after treatment with Letrozole [52]; however, SOX11 was found to be upregulated in our analysis. The fact that SOX11 was included in our protein network suggests that changes in its bulk RNA-seq expression could be driven by Letrozole treatment, though we are unsure whether this upregulation primarily takes place within tumor cells or in the tumor-associated immune cells [129]. We anticipate that additional experiments in cell culture and/or clinical cohorts will be necessary to validate these findings.

4.6. Study Limitations

We are aware that our analysis was not able to incorporate at least some of the relevant covariates due to the lack of publicly available metadata, which could have impacted the accuracy of our findings. A subset of these covariates could include the patient’s age, tumor size, clinical stage, ER and progesterone receptor status, and histologic subtype. In addition to these covariates, there are some questions that could not be answered because of the original study design that involved surgically collecting samples after treatment. As such, a direct analysis of samples collected before and after treatment is not possible.
We purposefully selected our cutoff criteria in the various analyses such that the statistical significance value was primarily used to determine the threshold. Our decision to not impose a fold-change cutoff enabled us to identify genes with relatively small but statistically significant changes in gene expression. This approach at least partially accounts for the cellular heterogeneity in bulk RNA-seq experiments [42,130,131]. We anticipate that adjusting our fixed threshold approach to a different set of cutoff values could alter the results.
Although our use of robust and established large-scale statistical methods should improve confidence in our results, the possibility of false-positives is inevitable. Consequently, we anticipate follow-up experiments in the wet lab to further validate our results. Additional studies that compare our results to those from triple-negative and/or ER- breast cancer patients could also be beneficial in efforts to anticipate precision medicine efforts.

5. Conclusions

Our analysis identified genes, pathways, and mechanistic transcriptional markers that could contribute to improved understanding, prognostics, and characterization of mechanisms for Letrozole resistance in ER+ breast cancer patients. Future work is needed to better characterize the identified gene products that do not display a well-understood function to better understand how the tumor microenvironment interacts with Letrozole. Further analysis of our predicted biomarkers could also yield clinical applications that enable treatment regimens specified to particular subsets of patients.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cimb46070424/s1, Table S1: Excel download for differentially expressed gene information. Table S2: Excel download for therapeutic targets.

Author Contributions

Conceptualization, B.E.P.; methodology, B.E.P.; software, B.E.P.; validation, L.S. and J.L.; formal analysis, L.S. and J.L.; investigation, L.S., J.L. and B.E.P.; data curation, L.S.; writing—original draft preparation, L.S., J.L. and B.E.P.; writing—review and editing, L.S., J.L., N.G.-J. and B.E.P.; visualization, L.S. and B.E.P.; supervision, B.E.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article and Supplementary Materials.

Acknowledgments

We acknowledge the technical expertise and the high-performance computing environment provided by the BYU Office of Research Computing.

Conflicts of Interest

B.E.P. owns equity in Pythia Biosciences.

References

  1. Giaquinto, A.N.; Sung, H.; Miller, K.D.; Kramer, J.L.; Newman, L.A.; Minihan, A.; Jemal, A.; Siegel, R.L. Breast Cancer Statistics, 2022. CA A Cancer J. Clin. 2022, 72, 524–541. [Google Scholar] [CrossRef]
  2. Postmenopause. Cleveland Clinic. Available online: https://my.clevelandclinic.org/health/diseases/21837-postmenopause (accessed on 8 June 2024).
  3. US Census Bureau. The Older Population: 2020; US Census Bureau: Washington, DC, USA, 2023.
  4. Breast Cancer Statistics. Available online: https://www.cancer.org/cancer/types/breast-cancer/about/how-common-is-breast-cancer.html (accessed on 8 June 2024).
  5. Rosenberg, P.S.; Barker, K.A.; Anderson, W.F. Estrogen Receptor Status and the Future Burden of Invasive and In Situ Breast Cancers in the United States. JNCI J. Natl. Cancer Inst. 2015, 107, djv159. [Google Scholar] [CrossRef]
  6. Breast Cancer Types: What your Type Means. Mayo Clinic. 2022. Available online: https://www.mayoclinic.org/diseases-conditions/breast-cancer/in-depth/breast-cancer/art-20045654 (accessed on 8 June 2024).
  7. Korourian, S.; Kumarapeli, A.R.; Klimberg, V.S. Breast Biomarker Immunocytochemistry. In The Breast; Elsevier: Amsterdam, The Netherlands, 2018; pp. 197–206.e6. [Google Scholar]
  8. Siddiqa, A.; Long, L.M.; Li, L.; Marciniak, R.A.; Kazhdan, I. Expression of HER-2 in MCF-7 breast cancer cells modulates anti-apoptotic proteins Survivin and Bcl-2 via the extracellular signal-related kinase (ERK) and phosphoinositide-3 kinase (PI3K) signalling pathways. BMC Cancer 2008, 8, 129. [Google Scholar] [CrossRef]
  9. Zambrano, J.; Yeh, E.S. Autophagy and Apoptotic Crosstalk: Mechanism of Therapeutic Resistance in HER2-Positive Breast Cancer. Breast Cancer 2016, 10, 13–23. [Google Scholar] [CrossRef]
  10. Onitilo, A.A.; Engel, J.M.; Greenlee, R.T.; Mukesh, B.N. Breast cancer subtypes based on ER/PR and Her2 expression: Comparison of clinicopathologic features and survival. Clin. Med. Res. 2009, 7, 4–13. [Google Scholar] [CrossRef]
  11. Calhoun, B.C.; Collins, L.C. Predictive markers in breast cancer: An update on ER and HER2 testing and reporting. Semin. Diagn. Pathol. 2015, 32, 362–369. [Google Scholar] [CrossRef]
  12. Caldarella, A.; Crocetti, E.; Bianchi, S.; Vezzosi, V.; Urso, C.; Biancalani, M.; Zappa, M. Female breast cancer status according to ER, PR and HER2 expression: A population based analysis. Pathol. Oncol. Res. 2011, 17, 753–758. [Google Scholar] [CrossRef]
  13. Bagaria, S.P.; Ray, P.S.; Sim, M.-S.; Ye, X.; Shamonki, J.M.; Cui, X.; Giuliano, A.E. Personalizing breast cancer staging by the inclusion of ER, PR, and HER2. JAMA Surg. 2014, 149, 125–129. [Google Scholar] [CrossRef]
  14. Bhatnagar, A.S. The discovery and mechanism of action of letrozole. Breast Cancer Res. Treat. 2007, 105 (Suppl. 1), 7–17. [Google Scholar] [CrossRef]
  15. Sohl, C.D.; Guengerich, F.P. Kinetic analysis of the three-step steroid aromatase reaction of human cytochrome P450 19A1. J. Biol. Chem. 2010, 285, 17734–17743. [Google Scholar] [CrossRef]
  16. Mukherjee, A.G.; Wanjari, U.R.; Nagarajan, D.; Vibhaa, K.K.; Anagha, V.; Chakraborty, R.; Renu, K.; Dey, A.; Vellingiri, B.; Gopalakrishnan, A.V. Letrozole: Pharmacology, toxicity and potential therapeutic effects. Life Sci. 2022, 310, 121074. [Google Scholar] [CrossRef]
  17. Ellis, M.J.; Ma, C. Letrozole in the neoadjuvant setting the P024 trial. Breast Cancer Res. Treat. 2008, 112, 371. [Google Scholar] [CrossRef]
  18. Guarneri, V.; Generali, D.G.; Frassoldati, A.; Artioli, F.; Boni, C.; Cavanna, L.; Tagliafico, E.; Maiorana, A.; Bottini, A.; Cagossi, K.; et al. Double-blind, placebo-controlled, multicenter, randomized, phase IIb neoadjuvant study of letrozole-lapatinib in postmenopausal hormone receptor-positive, human epidermal growth factor receptor 2-negative, operable breast cancer. J. Clin. Oncol. 2014, 32, 1050–1057. [Google Scholar] [CrossRef]
  19. Sachdev, A.; Sharpe, I.; Bowman, M.; Booth, C.M.; Gyawali, B. Objective response rate of placebo in randomized controlled trials of anticancer medicines. eClinicalMedicine 2023, 55, 101753. [Google Scholar] [CrossRef]
  20. Lee, K.M.; Guerrero-Zotano, A.L.; Servetto, A.; Sudhan, D.R.; Lin, C.C.; Formisano, L.; Jansen, V.M.; González-Ericsson, P.; Sanders, M.E.; Stricker, T.P.; et al. Proline rich 11 (PRR11) overexpression amplifies PI3K signaling and promotes antiestrogen resistance in breast cancer. Nat. Commun. 2020, 11, 5488. [Google Scholar] [CrossRef]
  21. Liu, Y.; Xu, J.; Choi, H.H.; Han, C.; Fang, Y.; Li, Y.; Van der Jeught, K.; Xu, H.; Zhang, L.; Frieden, M.; et al. Targeting 17q23 amplicon to overcome the resistance to anti-HER2 therapy in HER2+ breast cancer. Nat. Commun. 2018, 9, 4718. [Google Scholar] [CrossRef]
  22. Kodama, Y.; Shumway, M.; Leinonen, R.; International Nucleotide Sequence Database Collaboration. The Sequence Read Archive: Explosive growth of sequencing data. Nucleic Acids Res. 2012, 40, D54–D56. [Google Scholar] [CrossRef]
  23. Orjuela, S.; Huang, R.; Hembach, K.M.; Robinson, M.D.; Soneson, C. ARMOR: An utomated eproducible dular Workflow for Preprocessing and Differential Analysis of NA-seq Data. G3 2019, 9, 2089–2096. [Google Scholar] [CrossRef]
  24. Köster, J.; Rahmann, S. Snakemake—A scalable bioinformatics workflow engine. Bioinformatics 2018, 34, 3600. [Google Scholar] [CrossRef]
  25. Babraham Bioinformatics—FastQC A Quality Control Tool for High Throughput Sequence Data. Available online: www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 8 June 2024).
  26. Babraham Bioinformatics—Trim Galore! Available online: https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ (accessed on 8 June 2024).
  27. Patro, R.; Duggal, G.; Love, M.I.; Irizarry, R.A.; Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 2017, 14, 417–419. [Google Scholar] [CrossRef]
  28. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef]
  29. Wu, D.; Smyth, G.K. Camera: A competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res. 2012, 40, e133. [Google Scholar] [CrossRef]
  30. Tarca, A.L.; Draghici, S.; Khatri, P.; Hassan, S.S.; Mittal, P.; Kim, J.-S.; Kim, C.J.; Kusanovic, J.P.; Romero, R. A novel signaling pathway impact analysis. Bioinformatics 2009, 25, 75–82. [Google Scholar] [CrossRef]
  31. Gifford, K.T.; Pickett, B.E. Comparative meta-analysis of host transcriptional response during Streptococcus pneumoniae carriage or infection. Microb. Pathog. 2022, 173, 105816. [Google Scholar] [CrossRef]
  32. Ferrarini, M.G.; Lal, A.; Rebollo, R.; Gruber, A.J.; Guarracino, A.; Gonzalez, I.M.; Floyd, T.; de Oliveira, D.S.; Shanklin, J.; Beausoleil, E.; et al. Genome-wide bioinformatic analyses predict key host and viral factors in SARS-CoV-2 pathogenesis. Commun. Biol. 2021, 4, 590. [Google Scholar] [CrossRef]
  33. Viedma, M.d.P.M.; Pickett, B.E. Characterizing the Different Effects of Zika Virus Infection in Placenta and Microglia Cells. Viruses 2018, 10, 649. [Google Scholar] [CrossRef]
  34. Spendlove, M.D.; Gibson, T.M.; McCain, S.; Stone, B.C.; Gill, T.; Pickett, B.E. Pathway2Targets: An open-source pathway-based approach to repurpose therapeutic drugs and prioritize human targets. PeerJ 2023, 11, e16088. [Google Scholar] [CrossRef]
  35. Scott, T.M.; Solis-Leal, A.; Lopez, J.B.; Robison, R.A.; Berges, B.K.; Pickett, B.E. Comparison of Intracellular Transcriptional Response of NHBE Cells to Infection with SARS-CoV-2 Washington and New York Strains. Front. Cell. Infect. Microbiol. 2022, 12, 1009328. [Google Scholar] [CrossRef]
  36. Gray, M.; Guerrero-Arguero, I.; Solis-Leal, A.; Robison, R.A.; Berges, B.K.; Pickett, B.E. Chikungunya virus time course infection of human macrophages reveals intracellular signaling pathways relevant to repurposed therapeutics. PeerJ 2022, 10, e13090. [Google Scholar] [CrossRef]
  37. Scott, T.M.; Jensen, S.; Pickett, B.E. A signaling pathway-driven bioinformatics pipeline for predicting therapeutics against emerging infectious diseases. F1000Research 2021, 10, 330. [Google Scholar] [CrossRef]
  38. Ochoa, D.; Hercules, A.; Carmona, M.; Suveges, D.; Gonzalez-Uriarte, A.; Malangone, C.; Miranda, A.; Fumis, L.; Carvalho-Silva, D.; Spitzer, M.; et al. Open Targets Platform: Supporting systematic drug–target identification and prioritisation. Nucleic Acids Res. 2021, 49, D1302–D1310. [Google Scholar] [CrossRef] [PubMed]
  39. GitHub—dmlc/xgboost: Scalable, Portable and Distributed Gradient Boosting (GBDT, GBRT or GBM) Library, for Python, R, Java, Scala, C++ and More. Runs on Single Machine, Hadoop, Spark, Dask, Flink and DataFlow. GitHub. Available online: https://github.com/dmlc/xgboost (accessed on 8 June 2024).
  40. Rapier-Sharman, N.; Clancy, J.; E Pickett, B. Joint Secondary Transcriptomic Analysis of Non-Hodgkin’s B-Cell Lymphomas Predicts Reliance on Pathways Associated with the Extracellular Matrix and Robust Diagnostic Biomarkers. J. Bioinform. Syst. Biol. 2022, 5, 119–135. [Google Scholar] [CrossRef] [PubMed]
  41. Moreno, C.; Bybee, E.; Tellez Freitas, C.M.; Pickett, B.E.; Weber, K.S. Meta-Analysis of Two Human RNA-seq Datasets to Determine Periodontitis Diagnostic Biomarkers and Drug Target Candidates. Int. J. Mol. Sci. 2022, 23, 5580. [Google Scholar] [CrossRef] [PubMed]
  42. Clancy, J.; Hoffmann, C.S.; Pickett, B.E. Transcriptomics secondary analysis of severe human infection with SARS-CoV-2 identifies gene expression changes and predicts three transcriptional biomarkers in leukocytes. Comput. Struct. Biotechnol. J. 2023, 21, 1403–1413. [Google Scholar] [CrossRef] [PubMed]
  43. Piccolo, S.R.; Mecham, A.; Golightly, N.P.; Johnson, J.L.; Miller, D.B. The ability to classify patients based on gene-expression data varies by algorithm and performance metric. PLoS Comput. Biol. 2022, 18, e1009926. [Google Scholar] [CrossRef]
  44. Szklarczyk, D.; Kirsch, R.; Koutrouli, M.; Nastou, K.; Mehryary, F.; Hachilif, R.; Gable, A.L.; Fang, T.; Doncheva, N.T.; Pyysalo, S.; et al. The STRING database in 2023: Protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023, 51, D638–D646. [Google Scholar] [CrossRef]
  45. Su, G.; Morris, J.H.; Demchak, B.; Bader, G.D. Biological network exploration with Cytoscape 3. Curr. Protoc. Bioinform. 2014, 47, 8.13.1–8.13.24. [Google Scholar] [CrossRef] [PubMed]
  46. Chin, C.-H.; Chen, S.-H.; Wu, H.-H.; Ho, C.-W.; Ko, M.-T.; Lin, C.-Y. cytoHubba: Identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 2014, 8 (Suppl. 4), S11. [Google Scholar] [CrossRef] [PubMed]
  47. Ying, J.; Srivastava, G.; Hsieh, W.-S.; Gao, Z.; Murray, P.; Liao, S.-K.; Ambinder, R.; Tao, Q. The Stress-Responsive Gene GADD45G Is a Functional Tumor Suppressor, with Its Response to Environmental Stresses Frequently Disrupted Epigenetically in Multiple Tumors. Clin. Cancer Res. 2005, 11, 6442–6449. [Google Scholar] [CrossRef] [PubMed]
  48. Guerrero-Zotano, A.L.; Stricker, T.P.; Formisano, L.; Hutchinson, K.E.; Stover, D.G.; Lee, K.M.; Schwarz, L.J.; Giltnane, J.M.; Estrada, M.V.; Jansen, V.M.; et al. ER Breast Cancers Resistant to Prolonged Neoadjuvant Letrozole Exhibit an E2F4 Transcriptional Program Sensitive to CDK4/6 Inhibitors. Clin. Cancer Res. 2018, 24, 2517–2529. [Google Scholar] [CrossRef] [PubMed]
  49. Dey, A.; Kundu, M.; Das, S.; Jena, B.C.; Mandal, M. Understanding the function and regulation of Sox2 for its therapeutic potential in breast cancer. Biochim. Biophys. Acta (BBA)—Rev. Cancer 2022, 1877, 188692. [Google Scholar] [CrossRef] [PubMed]
  50. She, Z.-Y.; Yang, W.-X. SOX family transcription factors involved in diverse cellular events during development. Eur. J. Cell Biol. 2015, 94, 547–563. [Google Scholar] [CrossRef]
  51. Potzner, M.R.; Tsarovina, K.; Binder, E.; Penzo-Méndez, A.; Lefebvre, V.; Rohrer, H.; Wegner, M.; Sock, E. Sequential requirement of Sox4 and Sox11 during development of the sympathetic nervous system. Development 2010, 137, 775–784. [Google Scholar] [CrossRef]
  52. Shepherd, J.H.; Uray, I.P.; Mazumdar, A.; Tsimelzon, A.; Savage, M.; Hilsenbeck, S.G.; Brown, P.H. The SOX11 transcription factor is a critical regulator of basal-like breast cancer growth, invasion, and basal-like gene expression. Oncotarget 2016, 7, 13106–13121. [Google Scholar] [CrossRef]
  53. Wang, H.; Liu, B.; Long, J.; Yu, J.; Ji, X.; Li, J.; Zhu, N.; Zhuang, X.; Li, L.; Chen, Y.; et al. Integrative analysis identifies two molecular and clinical subsets in Luminal B breast cancer. iScience 2023, 26, 107466. [Google Scholar] [CrossRef]
  54. Zvelebil, M.; Oliemuller, E.; Gao, Q.; Wansbury, O.; Mackay, A.; Kendrick, H.; Smalley, M.J.; Reis-Filho, J.S.; Howard, B.A. Embryonic mammary signature subsets are activated in Brca1-/- and basal-like breast cancers. Breast Cancer Res. 2013, 15, R25. [Google Scholar] [CrossRef] [PubMed]
  55. Gebhardt, C.; Németh, J.; Angel, P.; Hess, J. S100A8 and S100A9 in inflammation and cancer. Biochem. Pharmacol. 2006, 72, 1622–1631. [Google Scholar] [CrossRef] [PubMed]
  56. Li, J.; Shu, X.; Xu, J.; Su, S.M.; Chan, U.I.; Mo, L.; Liu, J.; Zhang, X.; Adhav, R.; Chen, Q.; et al. S100A9-CXCL12 activation in BRCA1-mutant breast cancer promotes an immunosuppressive microenvironment associated with resistance to immunotherapy. Nat. Commun. 2022, 13, 1481. [Google Scholar] [CrossRef] [PubMed]
  57. Chen, Y.; Ouyang, Y.; Li, Z.; Wang, X.; Ma, J. S100A8 and S100A9 in Cancer. Biochim. Biophys. Acta (BBA)—Rev. Cancer 2023, 1878, 188891. [Google Scholar] [CrossRef]
  58. Rammes, A.; Roth, J.; Goebeler, M.; Klempt, M.; Hartmann, M.; Sorg, C. Myeloid-related Protein (MRP) 8 and MRP14, Calcium-binding Proteins of the S100 Family, Are Secreted by Activated Monocytes via a Novel, Tubulin-dependent Pathway. J. Biol. Chem. 1997, 272, 9496–9502. [Google Scholar] [CrossRef]
  59. Wang, L.; Di, L.-J. BRCA1 and estrogen/estrogen receptor in breast cancer: Where they interact? Int. J. Biol. Sci. 2014, 10, 566–575. [Google Scholar] [CrossRef] [PubMed]
  60. Karpisheh, V.; Ahmadi, M.; Abbaszadeh-Goudarzi, K.; Saray, M.M.; Barshidi, A.; Mohammadi, H.; Yousefi, M.; Jadidi-Niaragh, F. The role of Th17 cells in the pathogenesis and treatment of breast cancer. Cancer Cell Int. 2022, 22, 108. [Google Scholar] [CrossRef] [PubMed]
  61. Arai, K.; Takano, S.; Teratani, T.; Ito, Y.; Yamada, T.; Nozawa, R. S100A8 and S100A9 overexpression is associated with poor pathological parameters in invasive ductal carcinoma of the breast. Curr. Cancer Drug Targets 2008, 8, 243–252. [Google Scholar] [CrossRef] [PubMed]
  62. Boutas, I.; Kontogeorgi, A.; Koufopoulos, N.; Dimitrakakis, C.; Dimas, D.; Kalantaridou, S.N.; Korou, L.-M.; Perrea, D. Effect of Aromatase Inhibitors on Serum Calprotectin Levels in an Animal Experimental Model: Trial. Cancer Diagn. Progn. 2022, 2, 720–730. [Google Scholar] [CrossRef]
  63. Ralli, S.; Jones, S.J.; Leach, S.; Lynch, H.T.; Brooks-Wilson, A.R. Gene and pathway based burden analyses in familial lymphoid cancer cases: Rare variants in immune pathway genes. PLoS ONE 2023, 18, e0287602. [Google Scholar] [CrossRef]
  64. Patil, S.; Linge, A.; Grosser, M.; Lohaus, F.; Gudziol, V.; Kemper, M.; Nowak, A.; Haim, D.; Tinhofer, I.; Budach, V.; et al. Development and validation of a 6-gene signature for the prognosis of loco-regional control in patients with HPV-negative locally advanced HNSCC treated by postoperative radio(chemo)therapy. Radiother. Oncol. 2022, 171, 91–100. [Google Scholar] [CrossRef] [PubMed]
  65. Dai, D.; Xie, L.; Shui, Y.; Li, J.; Wei, Q. Identification of Tumor Microenvironment-Related Prognostic Genes in Sarcoma. Front. Genet. 2021, 12, 620705. [Google Scholar] [CrossRef] [PubMed]
  66. Li, R.; Liu, X.; Zhou, X.J.; Chen, X.; Li, J.P.; Yin, Y.H.; Qu, Y.Q. Identification and validation of the prognostic value of immune-related genes in non-small cell lung cancer. Am. J. Transl. Res. 2020, 12, 5844–5865. [Google Scholar] [PubMed]
  67. Liao, H.Y.; Da, C.M.; Liao, B. Roles of matrix metalloproteinase-7 (MMP-7) in cancer. Clin. Biochem. 2021, 92, 9–18. [Google Scholar] [CrossRef]
  68. Köhrmann, A.; Kammerer, U.; Kapp, M.; Dietl, J.; Anacker, J. Expression of matrix metalloproteinases (MMPs) in primary human breast cancer and breast cancer cell lines: New findings and review of the literature. BMC Cancer 2009, 9, 188. [Google Scholar] [CrossRef] [PubMed]
  69. Bouris, P.; Skandalis, S.S.; Piperigkou, Z.; Afratis, N.; Karamanou, K.; Aletras, A.J.; Moustakas, A.; Theocharis, A.D.; Karamanos, N.K. Estrogen receptor alpha mediates epithelial to mesenchymal transition, expression of specific matrix effectors and functional properties of breast cancer cells. Matrix Biol. 2015, 43, 42–60. [Google Scholar] [CrossRef] [PubMed]
  70. Gialeli, C.; Theocharis, A.D.; Karamanos, N.K. Roles of matrix metalloproteinases in cancer progression and their pharmacological targeting. FEBS J. 2011, 278, 16–27. [Google Scholar] [CrossRef] [PubMed]
  71. Sizemore, S.T.; Keri, R.A. The Forkhead Box Transcription Factor FOXC1 Promotes Breast Cancer Invasion by Inducing Matrix Metalloprotease 7 (MMP7) Expression. J. Biol. Chem. 2012, 287, 24631–24640. [Google Scholar] [CrossRef] [PubMed]
  72. Wang, B.; Zhan, Y.; Yan, L.; Hao, D. How zoledronic acid improves osteoporosis by acting on osteoclasts. Front. Pharmacol. 2022, 13, 961941. [Google Scholar] [CrossRef]
  73. Reid, I.R.; Green, J.R.; Lyles, K.W.; Reid, D.M.; Trechsel, U.; Hosking, D.J.; Black, D.M.; Cummings, S.R.; Russell, R.G.G.; Eriksen, E.F. Zoledronate. Bone 2020, 137, 115390. [Google Scholar] [CrossRef]
  74. Black, D.M.; Delmas, P.D.; Eastell, R.; Reid, I.R.; Boonen, S.; Cauley, J.A.; Cosman, F.; Lakatos, P.; Leung, P.C.; Man, Z.; et al. Once-yearly zoledronic acid for treatment of postmenopausal osteoporosis. N. Engl. J. Med. 2007, 356, 1809–1822. [Google Scholar] [CrossRef] [PubMed]
  75. Dedes, P.; Gialeli, C.; Tsonis, A.; Kanakis, I.; Theocharis, A.; Kletsas, D.; Tzanakakis, G.; Karamanos, N. Expression of matrix macromolecules and functional properties of breast cancer cells are modulated by the bisphosphonate zoledronic acid. Biochim. Biophys. Acta (BBA)—Gen. Subj. 2012, 1820, 1926–1939. [Google Scholar] [CrossRef] [PubMed]
  76. Chang, J.; Guo, C.; Li, J.; Liang, Z.; Wang, Y.; Yu, A.; Liu, R.; Guo, Y.; Chen, J.; Huang, S. EN1 Regulates Cell Growth and Proliferation in Human Glioma Cells via Hedgehog Signaling. Int. J. Mol. Sci. 2022, 23, 1123. [Google Scholar] [CrossRef] [PubMed]
  77. Majumder, A.; Singh, M.; Tyagi, S.C. Post-menopausal breast cancer: From estrogen to androgen receptor. Oncotarget 2017, 8, 102739–102758. [Google Scholar] [CrossRef] [PubMed]
  78. Kim, Y.J.; Sung, M.; Oh, E.; Van Vrancken, M.; Song, J.-Y.; Jung, K.; Choi, Y.-L. Engrailed 1 overexpression as a potential prognostic marker in quintuple-negative breast cancer. Cancer Biol. Ther. 2018, 19, 335–345. [Google Scholar] [CrossRef] [PubMed]
  79. Beltran, A.S.; Graves, L.M.; Blancafort, P. Novel role of Engrailed 1 as a prosurvival transcription factor in basal-like breast cancer and engineering of interference peptides block its oncogenic function. Oncogene 2013, 33, 4767–4777. [Google Scholar] [CrossRef]
  80. Yang, L.; Li, X.; Luo, Y.; Yang, T.; Wang, H.; Shi, L.; Feng, M.; Xie, W. Weighted gene co expression network analysis of the association between upregulated AMD1, EN1 and VGLL1 and the progression and poor prognosis of breast cancer. Exp. Ther. Med. 2021, 22, 1030. [Google Scholar] [CrossRef] [PubMed]
  81. Liu, Z.; Sun, Q.; Wang, X. PLK1, A Potential Target for Cancer Therapy. Transl. Oncol. 2017, 10, 22–32. [Google Scholar] [CrossRef]
  82. Chiappa, M.; Petrella, S.; Damia, G.; Broggini, M.; Guffanti, F.; Ricci, F. Present and Future Perspective on PLK1 Inhibition in Cancer Treatment. Front. Oncol. 2022, 12, 903016. [Google Scholar] [CrossRef]
  83. Wierer, M.; Verde, G.; Pisano, P.; Molina, H.; Font-Mateu, J.; Di Croce, L.; Beato, M. PLK1 Signaling in Breast Cancer Cells Cooperates with Estrogen Receptor-Dependent Gene Transcription. Cell Rep. 2013, 3, 2021–2032. [Google Scholar] [CrossRef]
  84. Kahl, I.; Mense, J.; Finke, C.; Boller, A.; Lorber, C.; Győrffy, B.; Greve, B.; Götte, M.; A Espinoza-Sánchez, N. The cell cycle-related genes RHAMM, AURKA, TPX2, PLK1, and PLK4 are associated with the poor prognosis of breast cancer patients. J. Cell. Biochem. 2022, 123, 581–600. [Google Scholar] [CrossRef] [PubMed]
  85. Gutteridge, R.E.A.; Ndiaye, M.A.; Liu, X.; Ahmad, N. Plk1 Inhibitors in Cancer Therapy: From Laboratory to Clinics. Mol. Cancer Ther. 2016, 15, 1427–1435. [Google Scholar] [CrossRef] [PubMed]
  86. Bhola, N.E.; Jansen, V.M.; Bafna, S.; Giltnane, J.M.; Balko, J.M.; Estrada, M.V.; Meszoely, I.; Mayer, I.; Abramson, V.; Ye, F.; et al. Kinome-wide Functional Screen Identifies Role of PLK1 in Hormone-Independent, ER-Positive Breast Cancer. Cancer Res. 2015, 75, 405–414. [Google Scholar] [CrossRef] [PubMed]
  87. Couchman, J.R. Syndecan-1 (CD138), Carcinomas and EMT. Int. J. Mol. Sci. 2021, 22, 4227. [Google Scholar] [CrossRef] [PubMed]
  88. Yang, Z.; Chen, S.; Ying, H.; Yao, W. Targeting syndecan-1: New opportunities in cancer therapy. Am. J. Physiol. Physiol. 2022, 323, C29–C45. [Google Scholar] [CrossRef] [PubMed]
  89. Fleurot, E.; Goudin, C.; Hanoux, V.; Bonnamy, P.-J.; Levallet, J. Estrogen receptor α regulates the expression of syndecan-1 in human breast carcinoma cells. Endocr.-Relat. Cancer 2019, 26, 615–628. [Google Scholar] [CrossRef] [PubMed]
  90. Barbareschi, M.; Maisonneuve, P.; Aldovini, D.; Cangi, M.G.; Pecciarini, L.; Mauri, F.A.; Veronese, S.; Caffo, O.; Lucenti, A.; Palma, P.D.; et al. High syndecan-1 expression in breast carcinoma is related to an aggressive phenotype and to poorer prognosis. Cancer 2003, 98, 474–483. [Google Scholar] [CrossRef] [PubMed]
  91. Al Kawas, H.; Saaid, I.; Jank, P.; Westhoff, C.C.; Denkert, C.; Pross, T.; Weiler, K.B.S.; Karsten, M.M. How VEGF-A and its splice variants affect breast cancer development—Clinical implications. Cell. Oncol. 2022, 45, 227–239. [Google Scholar] [CrossRef]
  92. Claesson-Welsh, L.; Welsh, M. VEGFA and tumour angiogenesis. J. Intern. Med. 2013, 273, 114–127. [Google Scholar] [CrossRef]
  93. Jeselsohn, R.; Buchwalter, G.; De Angelis, C.; Brown, M.; Schiff, R. ESR1 mutations—A mechanism for acquired endocrine resistance in breast cancer. Nat. Rev. Clin. Oncol. 2015, 12, 573–583. [Google Scholar] [CrossRef]
  94. Dustin, D.; Gu, G.; Fuqua, S.A.W. ESR1 mutations in breast cancer. Cancer 2019, 125, 3714–3728. [Google Scholar] [CrossRef]
  95. Brett, J.O.; Spring, L.M.; Bardia, A.; Wander, S.A. ESR1 mutation as an emerging clinical biomarker in metastatic hormone receptor-positive breast cancer. Breast Cancer Res. 2021, 23, 85. [Google Scholar] [CrossRef] [PubMed]
  96. Joseph, C.; Alsaleem, M.; Orah, N.; Narasimha, P.L.; Miligy, I.M.; Kurozumi, S.; Ellis, I.O.; Mongan, N.P.; Green, A.R.; Rakha, E.A. Elevated MMP9 expression in breast cancer is a predictor of shorter patient survival. Breast Cancer Res. Treat. 2020, 182, 267–282. [Google Scholar] [CrossRef]
  97. Cancemi, P.; Aiello, A.; Accardi, G.; Caldarella, R.; Candore, G.; Caruso, C.; Ciaccio, M.; Cristaldi, L.; Di Gaudio, F.; Siino, V.; et al. The Role of Matrix Metalloproteinases (MMP-2 and MMP-9) in Ageing and Longevity: Focus on Sicilian Long-Living Individuals (LLIs). Mediat. Inflamm. 2020, 2020, 8635158. [Google Scholar] [CrossRef] [PubMed]
  98. Nikolov, A.; Popovski, N. Role of Gelatinases MMP-2 and MMP-9 in Healthy and Complicated Pregnancy and Their Future Potential as Preeclampsia Biomarkers. Diagnostics 2021, 11, 480. [Google Scholar] [CrossRef]
  99. Owyong, M.; Chou, J.; Bijgaart, R.J.v.D.; Kong, N.; Efe, G.; Maynard, C.; Talmi-Frank, D.; Solomonov, I.; Koopman, C.; Hadler-Olsen, E.; et al. MMP9 modulates the metastatic cascade and immune landscape for breast cancer anti-metastatic therapy. Life Sci. Alliance 2019, 2, e201800226. [Google Scholar] [CrossRef]
  100. Mendes, O.; Kim, H.-T.; Lungu, G.; Stoica, G. MMP2 role in breast cancer brain metastasis development and its regulation by TIMP2 and ERK1/2. Clin. Exp. Metastasis 2007, 24, 341–351. [Google Scholar] [CrossRef] [PubMed]
  101. Mendes, O.; Kim, H.-T.; Stoica, G. Expression of MMP2, MMP9 and MMP3 in Breast Cancer Brain Metastasis in a Rat Model. Clin. Exp. Metastasis 2005, 22, 237–246. [Google Scholar] [CrossRef]
  102. Jiang, H.; Li, H. Prognostic values of tumoral MMP2 and MMP9 overexpression in breast cancer: A systematic review and meta-analysis. BMC Cancer 2021, 21, 149. [Google Scholar] [CrossRef]
  103. Somiari, S.B.; Somiari, R.I.; Heckman, C.M.; Olsen, C.H.; Jordan, R.M.; Russell, S.J.; Shriver, C.D. Circulating MMP2 and MMP9 in breast cancer—Potential role in classification of patients into low risk, high risk, benign disease and breast cancer categories. Int. J. Cancer 2006, 119, 1403–1411. [Google Scholar] [CrossRef]
  104. ClinicalTrials.gov. Available online: https://www.clinicaltrials.gov/study/NCT00003010 (accessed on 8 June 2024).
  105. Chew, N.J.; Nguyen, E.V.; Su, S.-P.; Novy, K.; Chan, H.C.; Nguyen, L.K.; Luu, J.; Simpson, K.J.; Lee, R.S.; Daly, R.J. FGFR3 signaling and function in triple negative breast cancer. Cell Commun. Signal. 2020, 18, 13. [Google Scholar] [CrossRef] [PubMed]
  106. Long, X.; Shi, Y.; Ye, P.; Guo, J.; Zhou, Q.; Tang, Y. MicroRNA-99a Suppresses Breast Cancer Progression by Targeting FGFR3. Front. Oncol. 2020, 9, 1473. [Google Scholar] [CrossRef]
  107. Xie, J.; Wan, Y.; Zhang, M.; Jin, Z.; Yao, Y. Circ_0061825 Acts as a miR-593-3p Sponge to Promote Breast Cancer Progression by Regulating FGFR3 Expression. Cancer Manag. Res. 2020, 12, 11243–11255. [Google Scholar] [CrossRef] [PubMed]
  108. Hinz, N.; Jücker, M. Distinct functions of AKT isoforms in breast cancer: A comprehensive review. Cell Commun. Signal. 2019, 17, 154. [Google Scholar] [CrossRef] [PubMed]
  109. Riggio, M.; Perrone, M.C.; Polo, M.L.; Rodriguez, M.J.; May, M.; Abba, M.; Lanari, C.; Novaro, V. AKT1 and AKT2 isoforms play distinct roles during breast cancer progression through the regulation of specific downstream proteins. Sci. Rep. 2017, 7, srep44244. [Google Scholar] [CrossRef] [PubMed]
  110. Sanchez, C.G.; Ma, C.X.; Crowder, R.J.; Guintoli, T.; Phommaly, C.; Gao, F.; Lin, L.; Ellis, M.J. Preclinical modeling of combined phosphatidylinositol-3-kinase inhibition with endocrine therapy for estrogen receptor-positive breast cancer. Breast Cancer Res. 2011, 13, R21. [Google Scholar] [CrossRef] [PubMed]
  111. Jia, W.; Chen, P.; Cheng, Y. PRDX4 and Its Roles in Various Cancers. Technol. Cancer Res. Treat. 2019, 18, 1533033819864313. [Google Scholar] [CrossRef]
  112. Elko, E.A.; Manuel, A.M.; White, S.; Zito, E.; van der Vliet, A.; Anathy, V.; Janssen-Heininger, Y.M. Oxidation of peroxiredoxin-4 induces oligomerization and promotes interaction with proteins governing protein folding and endoplasmic reticulum stress. J. Biol. Chem. 2021, 296, 100665. [Google Scholar] [CrossRef] [PubMed]
  113. Park, M.H.; Jo, M.; Kim, Y.R.; Lee, C.K.; Hong, J.T. Roles of peroxiredoxins in cancer, neurodegenerative diseases and inflammatory diseases. Pharmacol. Ther. 2016, 163, 1–23. [Google Scholar] [CrossRef]
  114. Kocatürk, B. In silico analysis reveals PRDX4 as a prognostic and oncogenic marker in renal papillary cell carcinoma. Gene 2023, 859, 147201. [Google Scholar] [CrossRef]
  115. Mei, J.; Hao, L.; Liu, X.; Sun, G.; Xu, R.; Wang, H.; Liu, C. Comprehensive analysis of peroxiredoxins expression profiles and prognostic values in breast cancer. Biomark. Res. 2019, 7, 16. [Google Scholar] [CrossRef] [PubMed]
  116. Rafiei, S.; Tiedemann, K.; Tabariès, S.; Siegel, P.M.; Komarova, S.V. Peroxiredoxin 4: A novel secreted mediator of cancer induced osteoclastogenesis. Cancer Lett. 2015, 361, 262–270. [Google Scholar] [CrossRef]
  117. Wang, G.; Zhong, W.-C.; Bi, Y.-H.; Tao, S.-Y.; Zhu, H.-X.; Xu, A.-M. The Prognosis of Peroxiredoxin Family in Breast Cancer. Cancer Manag. Res. 2019, 11, 9685–9699. [Google Scholar] [CrossRef] [PubMed]
  118. Ummanni, R.; Mundt, F.; Pospisil, H.; Venz, S.; Scharf, C.; Barett, C.; Fälth, M.; Köllermann, J.; Walther, R.; Schlomm, T.; et al. Identification of clinically relevant protein targets in prostate cancer with 2D-DIGE coupled mass spectrometry and systems biology network platform. PLoS ONE 2011, 6, e16833. [Google Scholar] [CrossRef] [PubMed]
  119. Farooqi, A.A.; Hou, M.-F.; Chen, C.-C.; Wang, C.-L.; Chang, H.-W. Androgen receptor and gene network: Micromechanics reassemble the signaling machinery of TMPRSS2-ERG positive prostate cancer cells. Cancer Cell Int. 2014, 14, 34. [Google Scholar] [CrossRef] [PubMed]
  120. Lee, D.Y.; Chun, J.N.; Cho, M.; So, I.; Jeon, J.H. Emerging role of E2F8 in human cancer. Biochim. Biophys. Acta (BBA)—Mol. Basis Dis. 2023, 1869, 166745. [Google Scholar] [CrossRef]
  121. Jiang, X.; Wang, J.; Deng, X.; Xiong, F.; Zhang, S.; Gong, Z.; Li, X.; Cao, K.; Deng, H.; He, Y.; et al. The role of microenvironment in tumor angiogenesis. J. Exp. Clin. Cancer Res. 2020, 39, 204. [Google Scholar] [CrossRef] [PubMed]
  122. Weijts, B.G.M.W.; Bakker, W.J.; A Cornelissen, P.W.; Liang, K.-H.; Schaftenaar, F.H.; Westendorp, B.; Wolf, C.A.C.M.T.d.; Paciejewska, M.; Scheele, C.L.G.J.; Kent, L.; et al. E2F7 and E2F8 promote angiogenesis through transcriptional activation of VEGFA in cooperation with HIF1. EMBO J. 2012, 31, 3871–3884. [Google Scholar] [CrossRef]
  123. Whittle, J.R.; Vaillant, F.; Surgenor, E.; Policheni, A.N.; Giner, G.; Capaldo, B.D.; Chen, H.R.; Liu, H.K.; Dekkers, J.F.; Sachs, N.; et al. Dual Targeting of CDK4/6 and BCL2 Pathways Augments Tumor Response in Estrogen Receptor–Positive Breast Cancer. Clin. Cancer Res. 2020, 26, 4120–4134. [Google Scholar] [CrossRef]
  124. Ye, L.; Guo, L.; He, Z.; Wang, X.; Lin, C.; Zhang, X.; Wu, S.; Bao, Y.; Yang, Q.; Song, L.; et al. Upregulation of E2F8 promotes cell proliferation and tumorigenicity in breast cancer by modulating G1/S phase transition. Oncotarget 2016, 7, 23757–23771. [Google Scholar] [CrossRef]
  125. Li, Y.; Huang, J.; Yang, D.; Xiang, S.; Sun, J.; Li, H.; Ren, G. Expression patterns of E2F transcription factors and their potential prognostic roles in breast cancer. Oncol. Lett. 2018, 15, 9216–9230. [Google Scholar] [CrossRef] [PubMed]
  126. Oshi, M.; Takahashi, H.; Tokumaru, Y.; Yan, L.; Rashid, O.M.; Nagahashi, M.; Matsuyama, R.; Endo, I.; Takabe, K. The E2F Pathway Score as a Predictive Biomarker of Response to Neoadjuvant Therapy in ER+/HER2- Breast Cancer. Cells 2020, 9, 1643. [Google Scholar] [CrossRef] [PubMed]
  127. Hoogstraat, M.; Lips, E.H.; Mayayo-Peralta, I.; Mulder, L.; Kristel, P.; van der Heijden, I.; Annunziato, S.; van Seijen, M.; Nederlof, P.M.; Sonke, G.S.; et al. Comprehensive characterization of pre- and post-treatment samples of breast cancer reveal potential mechanisms of chemotherapy resistance. NPJ Breast Cancer 2022, 8, 60. [Google Scholar] [CrossRef]
  128. Selli, C.; Turnbull, A.K.; Pearce, D.A.; Li, A.; Fernando, A.; Wills, J.; Renshaw, L.; Thomas, J.S.; Dixon, J.M.; Sims, A.H. Molecular changes during extended neoadjuvant letrozole treatment of breast cancer: Distinguishing acquired resistance from dormant tumours. Breast Cancer Res. BCR 2019, 21, 2. [Google Scholar] [CrossRef] [PubMed]
  129. Huang, H.; Zhou, J.; Chen, H.; Li, J.; Zhang, C.; Jiang, X.; Ni, C. The immunomodulatory effects of endocrine therapy in breast cancer. J. Exp. Clin. Cancer Res. 2021, 40, 19. [Google Scholar] [CrossRef] [PubMed]
  130. Pires, B.R.B.; Panis, C.; Alves, V.D.; Herrera, A.C.S.A.; Binato, R.; Pizzatti, L.; Cecchini, R.; Abdelhay, E. Label-Free Proteomics Revealed Oxidative Stress and Inflammation as Factors That Enhance Chemoresistance in Luminal Breast Cancer. Oxidative Med. Cell. Longev. 2019, 2019, 5357649. [Google Scholar] [CrossRef] [PubMed]
  131. Dabydeen, S.A.; Kang, K.; Díaz-Cruz, E.S.; Alamri, A.; Axelrod, M.L.; Bouker, K.B.; Al-Kharboosh, R.; Clarke, R.; Hennighausen, L.; Furth, P.A. Comparison of tamoxifen and letrozole response in mammary preneoplasia of ER and aromatase overexpressing mice defines an immune-associated gene signature linked to tamoxifen resistance. Carcinog. 2015, 36, 122–132. [Google Scholar] [CrossRef] [PubMed]
Cimb 46 00424 g001
Figure 1. A volcano plot of all differentially expressed genes (DEGs). Each dot represents a gene product, with blue dots representing those that were significantly downregulated (log2FC < 1; FDR-corrected p-values < 0.05) and red representing DEGs that were significantly upregulated (log2FC > 1.5; FDR-corrected p-value < 0.05, equivalent to −log10FDR of 1.3). The 15 most statistically significant DEGs (11 up-regulated and 4 down-regulated) are labeled.
Figure 1. A volcano plot of all differentially expressed genes (DEGs). Each dot represents a gene product, with blue dots representing those that were significantly downregulated (log2FC < 1; FDR-corrected p-values < 0.05) and red representing DEGs that were significantly upregulated (log2FC > 1.5; FDR-corrected p-value < 0.05, equivalent to −log10FDR of 1.3). The 15 most statistically significant DEGs (11 up-regulated and 4 down-regulated) are labeled.
Cimb 46 00424 g001
Cimb 46 00424 g002
Figure 2. The area under the receiver operator characteristic (AUROC) curve shows high predicted sensitivity vs. specificity using all transcripts. An initial machine learning analysis using all detected gene products as inputs showed high overall specificity (x-axis) and sensitivity (y-axis) for detecting Letrozole resistance in patients, with an AUC value of 97.2%.
Figure 2. The area under the receiver operator characteristic (AUROC) curve shows high predicted sensitivity vs. specificity using all transcripts. An initial machine learning analysis using all detected gene products as inputs showed high overall specificity (x-axis) and sensitivity (y-axis) for detecting Letrozole resistance in patients, with an AUC value of 97.2%.
Cimb 46 00424 g002
Cimb 46 00424 g003
Figure 3. The area under the receiver operator characteristic curve shows acceptable performance for the top two transcripts. A secondary computational analysis using the highest ranked mechanistic markers from the whole transcriptome prediction (PRDX4 and E2F8). This analysis showed the combination of these two gene products having acceptable specificity (x-axis) and sensitivity (y-axis) for detecting Letrozole resistance in patients, with an AUC value of 85.4%.
Figure 3. The area under the receiver operator characteristic curve shows acceptable performance for the top two transcripts. A secondary computational analysis using the highest ranked mechanistic markers from the whole transcriptome prediction (PRDX4 and E2F8). This analysis showed the combination of these two gene products having acceptable specificity (x-axis) and sensitivity (y-axis) for detecting Letrozole resistance in patients, with an AUC value of 85.4%.
Cimb 46 00424 g003
Cimb 46 00424 g004
Figure 4. Protein–protein interaction network of DEGs, predicted targets, and top predicted markers. A graphical representation showing the interactions (edges) for each of the proteins (nodes) from the STRING-DB database. The 87 gene products that directly interacted with each other were included in a single network. The input consisted of the 105 DEGs (blue-colored nodes), top five predicted therapeutic targets (salmon-colored nodes), and top five predicted transcriptional biomarkers (orange-colored nodes).
Figure 4. Protein–protein interaction network of DEGs, predicted targets, and top predicted markers. A graphical representation showing the interactions (edges) for each of the proteins (nodes) from the STRING-DB database. The 87 gene products that directly interacted with each other were included in a single network. The input consisted of the 105 DEGs (blue-colored nodes), top five predicted therapeutic targets (salmon-colored nodes), and top five predicted transcriptional biomarkers (orange-colored nodes).
Cimb 46 00424 g004
Table 1. Metrics for the original 58 ER+ breast cancer samples treated with Letrozole.
Table 1. Metrics for the original 58 ER+ breast cancer samples treated with Letrozole.
MetricsStudy Information
Gene Expression Omnibus (GEO) Study IdentifierGSE145325
TitleRNA sequencing of ER+ breast tumor treated with Letrozole
PlatformIllumina HiSeq 3000
Diagnostic criteriaResistance to estrogen suppression was defined by a preoperative endocrine prognostic index (PEPI) ≥ 4 and/or evidence of cancer relapse after a median follow-up of 5 years
Number of responders vs. non-responders36 vs. 22
Table 2. Metrics and information for the top five statistically significant differentially expressed genes.
Table 2. Metrics and information for the top five statistically significant differentially expressed genes.
Gene SymbolLog2 Fold-ChangeLog Counts Per Million (CPM)FDR p-ValueGene Name
SOX113.873.748.07 × 10−6SRY-box transcription factor 11
S100A94.535.061.77 × 10−5S100 calcium-binding protein A9
S100A84.72.961.77 × 10−5S100 calcium-binding protein A8
IGLV3–255.156.270.000184Immunoglobulin lambda variable 3–25
MMP74.114.840.000665Matrix metallopeptidase 7
Table 3. Metrics for all significantly affected intracellular signaling pathways and their predicted modulation.
Table 3. Metrics for all significantly affected intracellular signaling pathways and their predicted modulation.
RankPathway NameTotal # of Pathway Members# DEGs in PathwayBonferroni p-ValuePredicted Modulation
1PLK1 signaling events/anti-tumoral activity44235.84 × 10−10Activated
2Syndecan-1-mediated signaling events/leukocyte adhesion3080.00018429Activated
3FOXM1 proliferation-associated transcription factor network36140.0105863Activated
4HIF-1-alpha transcription factor network/oxygen homeostasis60170.01923722Inhibited
Table 4. Ranked list of top five repurposed therapeutic targets predicted to aid with the resistance phenotype.
Table 4. Ranked list of top five repurposed therapeutic targets predicted to aid with the resistance phenotype.
Target SymbolTarget NameWeighted Score (Higher Is Better)
VEGFAVascular endothelial growth factor A2220.5
ESR1Estrogen Receptor 11771
MMP9Matrix Metallopeptidase 91699
FGFR3Fibroblast Growth Factor Receptor 91589
AKT1AKT serine/threonine kinase 11571
Table 5. Description and metrics for the top five predicted mechanistic transcriptional markers from all detected genes, ranked in decreasing order by gain.
Table 5. Description and metrics for the top five predicted mechanistic transcriptional markers from all detected genes, ranked in decreasing order by gain.
FeatureGainCoverFrequencyAntibody AvailableLog2FC
(FDR p-Value)		Location
PRDX40.0320.0270.023Monoclonal0.807
(0.0326)		Secretory granules; the ER; and exosomes
E2F80.0310.0250.021Monoclonal1.4
(0.266)		Cytosol and nucleus
IQGAP30.0240.0200.016Monoclonal1.2
(0.266)		Cytosol and plasma membrane
ATP6V1C20.0210.0180.015Monoclonal2.44
(0.0263)		Cytosol; lysosomes; and exosomes
CDCA80.0170.0140.011Monoclonal1.37
(0.101)		Cytosol and nucleus
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Sutherland, L.; Lang, J.; Gonzalez-Juarbe, N.; Pickett, B.E. Secondary Analysis of Human Bulk RNA-Seq Dataset Suggests Potential Mechanisms for Letrozole Resistance in Estrogen-Positive (ER+) Breast Cancer. Curr. Issues Mol. Biol. 2024, 46, 7114-7133. https://doi.org/10.3390/cimb46070424

AMA Style

Sutherland L, Lang J, Gonzalez-Juarbe N, Pickett BE. Secondary Analysis of Human Bulk RNA-Seq Dataset Suggests Potential Mechanisms for Letrozole Resistance in Estrogen-Positive (ER+) Breast Cancer. Current Issues in Molecular Biology. 2024; 46(7):7114-7133. https://doi.org/10.3390/cimb46070424

Chicago/Turabian Style

Sutherland, Lincoln, Jacob Lang, Norberto Gonzalez-Juarbe, and Brett E. Pickett. 2024. "Secondary Analysis of Human Bulk RNA-Seq Dataset Suggests Potential Mechanisms for Letrozole Resistance in Estrogen-Positive (ER+) Breast Cancer" Current Issues in Molecular Biology 46, no. 7: 7114-7133. https://doi.org/10.3390/cimb46070424

APA Style

Sutherland, L., Lang, J., Gonzalez-Juarbe, N., & Pickett, B. E. (2024). Secondary Analysis of Human Bulk RNA-Seq Dataset Suggests Potential Mechanisms for Letrozole Resistance in Estrogen-Positive (ER+) Breast Cancer. Current Issues in Molecular Biology, 46(7), 7114-7133. https://doi.org/10.3390/cimb46070424

Article Metrics

Back to TopTop