Novel Predictors of Breast Cancer Survival Derived from miRNA Activity Analysis.

The large public cancer dataset, the Cancer Genome Atlas (TCGA), contains multiple types of data (e.g. DNA copy number arrays, DNA methylation, exome sequencing, messenger RNA array, microRNA sequencing and reverse-phase protein arrays) performed on the same set of samples and provides comprehensive molecular portraits of several subtypes of breast cancers ( 26 ). We previously reported that high activities of let-7d and miR-18a were associated with better survival outcomes ER−/HER2− breast cancer ( 10 ). However, the majority of breast tumors in the United States are hormone receptor positive. We previously showed that the miRNA-mRNA correlation structure for ER + and ER - /HER2 - is markedly different, suggesting miRNA regulatory mechanisms are cancer-subtype specific ( 10 ). We subsequently evaluated relationships between miRNA and mRNAs expression levels in ER+ breast cancer using public databases. For each miRNA-mRNA pair, we calculated Pearson correlations between miRNA and mRNA expression levels. We inferred miRNA activity based on ActMiR method that we developed ( 10 ). miR-500a emerged as a key prognostic miRNAs for ER + breast cancer. In the TCGA dataset, high “activity” of miR-500a was strongly associated with increased breast cancer mortality ( Fig. 1A ); a similar association was observed using data from Dvinge et al. ( 14 ). It is worth noting that in comparison with the “activity”, the “expression” level of miR-500a showed no associations with breast cancer survival ( Fig. 1B ). Given that the activity of miR-500a can only be inferred by its associated genes, we profiled 195 candidate genes (161 were miR-500a associated genes) in tumors from 606 LIBCSP patients. As shown in Table 1 , among the profiled patients, 70% are postmenopausal, 20% have first degree relatives with breast cancer, and 75% and 62% are ER+ and PR+, respectively; these distributions are all similar to the overall patient cohort ( 11 ). The only notable difference is that there were fewer in situ (7% vs. 16%) and small (≤2 cm) tumors (72% vs. 80%) represented in the profiled patients compared to the overall cohort, likely due to the limited RNA yield in small and in situ tumors. From the 195 genes profiled, 98% (n=186) were detectable in more than 50% of the samples ( Fig. S3 ). Among the detectable genes, median raw counts varied from 10 ( KAT8 probe) to 12547 ( FTH1P11 probe) for candidate genes and the 6 housekeeping genes ranged from 66 ( HPRT1 probe) to 1956 ( GAPDH probe). Using Cox proportional hazards models, we calculated HRs for each of miR-500a associated genes in relation to breast cancer-specific mortality, adjusted for age and tumor stage, among the LIBCSP population-based cohort of women with breast cancer. Fig. 2 displays log-transformed p-values (y-axis) and HRs (x-axis) for these associations. Detailed results are reported in Table S3 . From 161 miR-500a associated genes, 74 displayed significant (FDR<0.05) associations with breast cancer-specific mortality. With increasing gene expression levels, 31 were associated with reduced breast cancer-specific mortality (HR1). Many well-known breast cancer-related genes (e.g. ESR1, AR ) were among those demonstrating significant association with mortality. The strongest (by p-value) inverse association for breast cancer-specific mortality was for SUSD3 (HR=0.3, 95% CI 0.2 - 0.4, p=5.2e-10), and the strongest positive association was for BIRC5 (HR=3.3, 95% CI 2.3 - 4.8, p=5.7e-10). These 74 genes also demonstrated high degree of correlation with each other; the correlation matrix ( Fig. 3 ) indicates two distinct expression clusters with opposing association with breast cancer mortality; such clustering was also supported by principal component analysis ( Fig. S4 ): PC1 showed a positive association with breast cancer mortality (HR = 4.2, 95% CI 2.9-6.2, p=4.8e-14), whereas PC2 showed the opposite (HR = 0.5, 95% CI 0.4-0.8, p=5.3e-4). We selected the genes to be profiled based on their association with miR-500a “activity”, both positive as well as negative associations. While a negative association suggests “direct targets” which represents the conventional mechanism of miRNA function, a positive association, which was also frequently observed, implicates mechanisms of indirect target or other modes of action. Fig. 4 depicts associations between HRs of each tested gene and its correlation with miR-500a activity by ActMiR algorithm for TCGA dataset. Based on our a priori hypothesis that high miR-500a activity confers increased breast cancer mortality, miR-500 associated genes should reside in Quadrant II or III, which is largely consistent with our findings ( Fig. 4 ). There are a small number of genes resided in Quadrant I and IV implicating unknown modes of action. It is important to point out that 22 (including BIRC5) of the 161 miR-500a genes are also part of PAM50 gene panel (pink dots on Fig. 4 ); 16 showed significant association with breast cancer mortality (FDR<0.05) Given that hormone receptor status and PAM50-derived risk scores are established factors associated with breast cancer survival, we set out to investigate whether miR-500a associated genes can predict survival, once we adjusted for these predictors. Among the 606 tumors profiled, 78% were ER+ and 67% were PR+ resulting in 80% ER+ or PR+ (ER+/PR+) and 20% ER− and PR− (ER−PR−) based on gene expression-based categorization. These numbers are similar to what have been reported by other studies of Caucasian Americans based on IHC status ( 27 – 29 ) (range from 81%−86% for ER+/PR+). We first re-evaluated associations between miR-500a associated genes and breast cancer-specific mortality by adjusting for ER/PR status; 43 out of 73 genes initially associated with breast cancer mortality remained significant (FDR<0.05). When we performed stratified analyses with respect to ER/PR status, 61 genes were significantly associated with breast cancer-specific mortality in ER+/PR+ group ( BIRC5 and SUSD3 remained the strongest predictors of higher and lower mortality, respectively) while only 2 genes were significant for ER−PR− group ( TMPRSS2 : HR=3.4, 95% CI 1.8 - 6.5; CLDN7 : HR=3.3, CI 1.7 - 6.3) ( Table S3 ). We then investigated whether miR-500 associated genes were associated with breast cancer mortality, once we adjusted/stratified for PAM50 risk scores. All PAM50 genes were at measurable levels in all 606 LIBCSP tumor samples. Risk-Of-Recurrence (ROR) scores (high, medium and low) were computed using the algorithm ( 15 ) that takes into account expression of PAM50 genes as well as the patient’s age and tumor stage. The associations of ROR scores and breast cancer-specific mortality are shown in Fig. S5 . The survival curves were similar among patients with low and medium ROR scores so we combined them as reference group ( Fig. S5 ). As expected, high ROR was associated with increased breast cancer-specific mortality (HR=3.6, 95% CI 2.4 - 5.2, p=1.2e-10) compared to the reference group. When we evaluated ROR-adjusted HRs for all 73 genes initially associated with breast cancer mortality, 31 remained significant (FDR < 0.05). We then stratified by ROR scores, high vs. low/medium combined. In the high-risk group (249 patients, 81 deaths), 10 genes were significantly associated with breast cancer mortality (FDR<0.05), with AKT1 (HR=3.0, 95% CI 1.9-4.8) and CROT (HR=0.3, 95% CI 0.1-0.6) showing the strongest positive and negative correlations with breast cancer mortality, respectively. In the low/medium risk group (357 participants, 38 death), TMPRSS2 (HR=4.0, 95% CI 2.1-7.8) and SUSD3 (HR=0.3, 95% CI: 0.1-0.5) demonstrated significant associations with survival. Additional evaluation of the predictive value of miR-500a-related genes was performed using a calculation of receiving operator characteristics (ROC analysis.). One of the advantages of this approach is that it is not based on the “optimal threshold” but takes the whole range of predictor’s values. As a reference, we first calculated AUC (area under curve) for the risk-of-recurrence score. Without adjustment for the survival length, AUC for the PAM-50 derived ROR-P score was very similar to the AUC of our most significantly associated gene SUSD3 (0.70 vs 0.69, Fig. 5A ). However when we applied time-dependent modeling ( 25 ) using survivalROC package to consider survival length (within one year after diagnosis), SUSD3 alone outperformed ROR-P score (AUC 0.86 vs 0.62, Fig. 5B ). Lastly, we validated our findings in two public datasets, i.e. METABRIC dataset ( 30 ) which includes 1980 breast cancer cases (77% ER+), and a combined dataset of KMPlot service kmplot.com ( 31 ) with 2862 cases (72% ER+). Most significant survival-associated genes from our study displayed the same direction of survival association in those datasets. ( Fig. S6 , Fig. S7 ). As an example, Fig. S8 shows data for our top gene SUSD3 : it is consistent among all 3 datasets, and for both ER+ and ER− subsets.

We utilized resources from the follow-up component of the Long Island Breast Cancer Study Project (LIBCSP), which includes a population-based cohort of women newly diagnosed with first primary breast cancer who participated in the baseline interview within two to three months of diagnosis ( 11 ) and were subsequently re-interviewed about five years later and followed for vital status ( 12 ). Participants were 1508 women newly diagnosed with first primary breast cancer between August 1, 1996 and July 31, 1997, who resided in Nassau and Suffolk counties of Long Island, NY, predominantly (90%) Caucasian, and aged 25-98 (median age 58). Potentially eligible women were identified through daily or weekly contact with pathology/cytology departments of 33 collaborating institutions, and confirmed by medical record. The National Death Index (NDI) was used to ascertain all-cause and breast cancer-specific mortality among study participants, who were followed from diagnosis until December 31, 2014, for an average of 162 months (range 2.7 – 224 months, median 213 months). We obtained sufficient RNA for profiling from paraffin-embedded tumor tissue from 616 of the 1508 patients (for the other LIBCSP patients, tumor slides were not made available by the diagnosing hospitals (n=763) or the extracted RNA was of insufficient/low quality (n=139); of these, 606 further passed QC. Demographic description of both the overall patient cohort and the subset with sufficient tissue RNA is presented in Table S1 . Among the 606 women with successfully profiled tumor tissue, 275 deaths occurred by the end of the follow-up period, of which 119 (43%) were related to breast cancer (ICD codes 174.9 and C-50.9 listed as a primary or secondary code on the death certificate). The study protocol was approved by the Institutional Review Boards of the collaborating institutions. We recently developed a computational method, ActMiR, which infers miRNA activity based on expression levels of miRNAs and their potential target genes ( 10 ). In brief, three pieces of information are used in this method: (i) a list of predicted targets of a given miRNA; (ii) expression levels of this miRNA samples of interest; (iii) mRNA expression levels of predicted targets in the same samples. This approach infers “functional activity” of a given miRNA using a regression-based model based on the expression levels of all miR predicted target genes. For the predicted target list of miRNAs, we used a collection of predicted target genes for 1537 unique mature miRNAs from TARGETSCAN ( www.targetscan.org ) that considers all conserved miRNA binding sites inherited from 23-way alignments of UTR sequences ( 13 ). In order to obtain robust results, we filtered out miRNAs whose number of target genes was smaller than 10. High miRNA activity corresponds to the high degradation effect on their target genes. To identify microRNAs of which activity was associated with breast cancer survival, we used “The Cancer Genome Atlas” database ( tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp ) and data described in Dvinge et al. ( 14 ). We previously published a novel computational method, ActMiR , for inferring microRNA activity ( 10 ). High miR-500a “activity” was predicted to be associated with increased ER+ breast cancer mortality (see Supplementary Methods ). Because the “activity” of miR-500 can only be inferred by expression of its mediated genes, we selected miR-500a mediated genes to be profiled in the LIBCSP tumor samples. They included genes that were significantly correlated with imputed miR-500a activity and also those differentially expressed after modulation of miR-500a in breast cancer cells. A total of 161 miR-500a mediated genes were selected. We also included all genes of the PAM50 panel by NanoString ( 15 ) in order to take into account of known risk-of-recurrence factors; 22 genes of the PAM50 set overlapped with the miR-500a-mediated gene set. Lastly, 6 genes were used as housekeeping controls and for normalization: CLTC , GAPDH , GUSB , HPRT1 , PGK1 and TUBB . These bring the total number of genes profiled to 195; the complete list is shown in Table S2 . Formalin-fixed paraffin-embedded (FFPE) tumor sections from the population-based cohort were histopathologically reviewed by a trained pathologist and the cancer tissue was separated using manual microdissection ( 16 ). Total RNA was extracted from tumor tissues using the Qiagen miRNeasy FFPE kit (Qiagen, MD). RNA concentration and quality were determined with a NanoDrop spectrophotometer (Thermo Scientific, MA). In total, 606 samples with sufficient RNA concentration were used in our analyses. Gene expression analysis was performed with the NanoString nCounter platform (Seattle, WA), using a custom-designed codeset and following the manufacturer’s instructions. Briefly, 100 ng RNA was incubated with reporter and capture probes overnight at 65°C. Following hybridization, unbound probes were removed, and the purified complexes were aligned and immobilized on imaging cartridges using an nCounter Prep station. Code count detection was carried out by scanning cartridges in an nCounter Digital Analyzer to determine gene expression levels. As a quality control filter, we excluded all samples having less than 50000 total counts for positive probes, or having less than 1.5*blank counts for housekeeping and endogenous genes. Raw counts were imported from RCC files with the NanoStringQCPro R package (version 1.6.0), background was subtracted as mean of negative controls counts, geometrical mean of housekeeping genes counts was used for normalization, and values were log-transformed for further analysis. Coexpression of profiled genes was analyzed in R by building a covariance matrix; principal component analysis was performed with “ pca” function of “ pcaMethods” package ( 17 ). We were able to obtain hormone receptor (HR) status determined by immunohistochemistry (IHC) from medical records for 475 out of 606 profiled samples, among which 75% and 62% were positive for estrogen receptor (ER+) and progesterone receptor (PR+), respectively. In addition, mRNA levels corresponding ESR1 and PGR were measured by NanoString in all tumor samples. For analyses in this study, we redefined ER/PR status based on ESR1 / PGR mRNA expression levels (see Supplementary Methods ). The expression-based method resulted in prevalence of 78% for ER+ (kappa=0.77) and 67% for PR+ (kappa=0.66) breast tumors ( Fig. S1 , Fig. S2 ). Recently, more advanced prognostic tools based on larger gene panels, such as PAM50-derived “Prosigna ™ ” , have been approved for clinical use and shown to improve breast cancer survival ( 18 ). In this study, we profiled expression of PAM50 ( 15 ) genes and determined molecular subtype and risk-of-recurrence (ROR) score using R code of Joel S. Parker from https://genome.unc.edu/pubsup/breastGEO/PAM50.zip archive ( 15 ). Risk-of-recurrence score was calculated in the “Subtype + Proliferation” variant. Cox proportional hazards models ( 19 ) were used for estimating the multivariate-adjusted hazard ratios (HRs) and 95% confidence intervals (CIs) for breast cancer-specific mortality associated with gene expression among the 606 women in our population-based cohort. Survival analysis was performed using “ survival” package in R ( 20 ). For each gene, range of cutoff thresholds (from second to eighth decile) was tested and the threshold with the lowest p-value was used ( 21 ). For each gene, a Cox model was built using coxph function with age and tumor stage (invasive or in situ ) included as cofactors. Multiple comparison adjustment was applied, where indicated, with the false discovery rate (FDR) threshold at 0.05. All models were initially adjusted for age at diagnosis (continuous) and tumor stage. We further evaluated potential confounding using the methods described by Rothman and Greenland ( 22 ) starting with a full multivariate model and using backward elimination. Potential confounder considered included: menopausal status (pre-/post-menopausal), family history of breast cancer in a first-degree relative, body mass index (BMI) at diagnosis, tumor size and education. If eliminating a covariate from the full Cox regression model changed the effect estimate by 10% or more, the covariate was considered a confounder and kept in the model. Otherwise that covariate was dropped from the multivariate model. None of the covariates tested met such criterion, thus, only results adjusted for age and tumor stage were presented. We also performed stratified analyses by ER/PR status as well as ROR scores discussed above. Kaplan-Meier survival ( 23 ) plots were built for visualization of unadjusted associations between gene expression and breast cancer-specific mortality. To estimate receiving operator characteristics (ROC), we first used R package pROC ( 24 ) . This package allows to perform a standard ROC analysis and calculates the Area Under the Curve (AUC) for binary outcome, i.e. censored (death by breast cancer) vs. uncensored (survived at the end of followup). However, this is not the optimal estimator for such time-dependent variable as patients survival, therefore we further re-analyzed our data using the survivalROC package ( 25 ) . This package takes into account the survival duration and is considered more appropriate for this type of data ( 25 ) . We used Nearest Neighbor Estimation (NNE) variant of the survivalROC function, with time point of 365 days and span parameter of 0.07.

Despite the fact that overall breast cancer mortality rates declined 36% from 1989 to 2012 (2015-2016 Report from American Cancer Society), breast cancer remains the second leading cause of cancer death among US women. The decline in breast cancer mortality has been partially attributed to improvements of treatment, in which molecular phenotyping has played a significant role. Hormone receptor-positive breast cancer represents more than 80% of patients among US women ( 32 ). This category of patients is typically treated with adjuvant endocrine therapy, with or without adjuvant chemotherapy ( 4 ). Whether the patient should receive chemotherapy is a difficult choice between risk of recurrence and adverse effects of the therapy: while treatment is usually toxic and decreases patients’ quality-of-life, it is often used to prevent relapse of the tumor. Early approaches of estimating a high-risk subgroup were based on “clinical” characteristics of disease including tumor size, presence/number of lymph node metastases, and histology of tumor tissue ( 33 ), then came the inclusion of immunohistological staining for hormone receptor proteins (estrogen and progesterone receptors) ( 33 ). The simplest assay, IHC4, is based on just four immunohistochemical markers (ER, PR, HER2 and Ki67) determined in an FFPE tumor tissue sample, and was shown to have good performance in a centralized study when combined with clinical characteristics ( 34 ). However, usage of this assay in decentralized conditions is much less robust ( 33 ). All IHC-based assays share similar technical limitations: staining is highly dependent on the quality of sample, variability of antibodies used, and subjective nature of scoring of the intensity of staining. Gene expression-based assays minimize these limitations and allow measurement of multiple targets in one assay. As a result, almost all recently developed assays are based on measuring mRNA instead of protein levels. Currently there are several approved gene expression panels for clinical use, three of which have the majority of the market: 21-gene Oncotype DX ™ , 50-gene Prosigna ™ and 70-gene Mammaprint ™ ( 35 ). MicroRNAs represent a special type of molecular mechanism of gene regulation, and have the capacity to repress almost any gene having a corresponding target sequence. Their involvement in cancer-related pathways has been clearly demonstrated, and including miRNA-related scores in disease prediction panels is a logical next step ( 7 ). It should be noted, however, that not all platforms designed for mRNA level detection are compatible with miRNA measurement, mainly due to shorter length of miRNAs. More importantly, miRNA “expression” level is not equivalent to its regulatory “activity” ( 8 , 36 ), because activities of miRNAs are modulated through the accessibility of proteins or RNAs that is essential for miRNA functions and the relative abundance of miRNA to its target genes ( 37 – 39 ). Thus, miRNA expression levels might not be the best functional indicator. For our cases, miR-500a expression levels were measured by qPCR on a small subset of paired tumor/adjacent premalignant tissues and no significant difference was observed ( Fig. S9 ); therefore we did not measure miR-500a expression in the entire cohort. In our study, we reported that miR-500a activity (but not expression) was associated with survival of ER-positive patients with early stage breast cancer in two public databases. We profiled a set of miR-500a associated genes, including potential direct and indirect targets, and aimed to investigate their prognostic significance in the well-characterized population-based cohort of breast cancer patients with long-term follow-up. Our study revealed a set of genes strongly associated with the survival of breast cancer patients, in agreement with what is predicted by ActMiR ( 10 ). Moreover, most of these associations were validated in two independent populations, including a large dataset of METABRIC ( 30 ) and the combined dataset provided by KMplot.com portal ( 31 ). At the same time, gene expression data from the set of clinical samples used in this study are not sufficient to claim that ActMiR-calculated `miR-500a activity` matches real biological activity: more specific in vitro experiments need to be performed in order to justify that type of equivalence. However, the direction of associations observed in this study was largely consistent with what was predicted by ActMiR (see Fig. 4 ) which strongly suggests that the panel of genes we assayed was indeed related to miR-500a activity. Although many of the predictor genes revealed in our study are novel for breast cancer prognosis, it does not mean they are better predictors of breast cancer survival because they can be highly correlated with established prognostic markers. To address this issue, we re-iterated survival analysis by taking into account two well-established prognostic markers (ER/PR status and PAM50-derived ROR score). Large proportion of genes that showed significant associations in overall breast cancer patients remained significant, even after adjusting for ER/PR status (43 out of 73) or PAM50-derived ROR scores (31 out of 73). Among the high-risk (ER−PR− or high ROR) or low-risk (ER+/PR+, or low/medium ROR) group, we were able to identify a set of genes that show significant associations with survival and differentiate the group with different survival profiles. This is an important finding because it might potentially help to better classify patients and improve precision of treatment or general care for these patients. Although miR-500a has been identified as one of key players in breast cancer mortality, its exact role remains unknown. TargetScan ( 40 ) predicts that as many as 3816 possible target transcripts for this microRNA but none of them are experimentally validated. There is only a handful of publications which suggest that on miR-500a is a potential circulating biomarker of such pathologies as autism ( 41 ), lung cancer ( 42 ), endometriosis ( 43 ), and hepatocellular carcinoma ( 44 ). Our study discovered a large number of miR-500a genes that demonstrated strong associations with breast cancer survival. Many of these associated with poor survival in our population-based cohort of women with breast cancer are already part of PAM50 panel (such as BIRC5, CEP55, CENPF – all known to be markers of poor survival). However, there were many others where we are the first to be report that they are predictive of breast cancer survival, even once we adjusted for the PAM50 ROR score. For example, high expression of SUSD3 (Sushi Domain-Containing Protein 3) displayed significantly reduced mortality for both HR-positive and HR-negative breast cancer cases in our population, and this finding was validated in both the METABRIC and combined KMplot datasets ( Fig. S6 ). Its association remained significant after adjustment for ROR (HR 0.4, 95% CI 0.3 - 0.7). More importantly, SUSD3 level effectively stratified low-ROR group (HR 0.3, 95% CI 0.1 - 0.5) identifying the patients with poor survival within this group. According to AUC values, in our study the SUSD3 alone performed equaly well with ROR score in its predictability of breast cancer survival and performed even better when adjusted for survival length within one-year window ( Fig. 5 ). This gene was shown to be a survival biomarker based on METABRIC data ( 45 ), and increased expression of corresponding protein was reported in breast cancer tumor tissue ( 46 ). At the same time, there was only one previous study focused on deciphering molecular mechanism of this gene. Moy et al. ( 47 ) showed that SUSD3 was a direct transcriptional target of estrogen receptor in MCF7 breast cancer cells, and its expression was induced after stimulation with estradiol. SUSD3 product was shown to be localized at plasma membrane, and its knockdown lead to decreased cell motility and proliferation ( 47 ). These results cannot easily fit with the inverse association with breast cancer mortality in our study. The molecular mechanism of SUSD3 deserves to be further investigated experimentally. Another gene highly correlated with breast cancer survival after adjustment for PAM50-derived ROR scores (HR 0.4, 95% CI 0.2 - 0.6) was CROT (Carnitine O-Octanoyltransferase). Again, its association was highly significant not only in our dataset but also by KMplot data (METABRIC data not available for CROT: gene symbol not found in the list). This gene is known to be involved in lipid acids metabolism, but was not previously linked to cancer. Similar inverse association was shown for ARMCX1 (armadillo repeat containing, X-linked 1) gene in patients from all three datasets, including ROR-adjusted LIBCSP patients (HR 0.5, 95% CI 0.4 - 0.8), with other studies having confirmed it as potential tumor suppressor gene ( 48 ). While our study has many strength including a well-characterized population-based cohort with long follow-up coupled with a novel systems biology approach, we acknowledge limitations of the study. Part of those limitations come from the lack of detailed clinical information. For example, we did not have integrated staging information (I – IIA – IIB etc.), and only used invasiveness status as an available stage. In addition, we did not perform in-depth analysis on treatment, which may influence breast cancer survival. Such exclusion stems from several reasons. First, only ~66% participants completed the follow-up survey, which may introduce bias to the cohort. Second, only information on the first course of treatment was collected while future therapy could be different depending on disease progression. Third, multi-center nature of this study (more than 30 participating hospitals) assumes there could be certain heterogeneity in choice of the treatment strategies and the way it was reported. Finally and importantly, by definition, a confounder (i.e. treatment) has to be related to both the independent variable (i.e. gene expression) and outcome (i.e. survival). In our case, all tumor samples were collected prior to treatment other than surgery, the treatment and molecular profile of the tumor should be considered independent. Also, treatment is part of the causal pathway (being impacted by tumor stage or hormone receptor status at diagnosis) and including a causal intermediate in the model can bias the effect estimates ( 49 ). In summary, by leveraging public databases and our own epidemiologic studies and cooperating systems biology approach, we identified a key microRNA, miR-500a, that appeared to be strongly linked with breast cancer mortality across platforms. More importantly, we identified a miR-500a associated gene panel that was not only strongly associated with breast cancer mortality but also can differentiate survival profiles with existing prognostic factors such as hormone receptor status and PAM50 derived risk of recurrence. These results suggest miR-500a associated genes may improve prognostic efficiency of current molecular phenotyping and ultimately improve breast cancer survival. Our study on these genes sheds light on breast cancer progression mechanisms and may lead to development of new targets for cancer therapy.

Recent advances in molecular phenotyping of breast cancer have significantly advanced our capabilities to treat breast cancer based on molecular features of the tumor and ultimately to improve survival ( 1 ). Different gene expression panels, such as Oncotype DX ™ , MammaPrint ™ and Prosigna ™ , have been developed and approved for clinical use for predicting risk of recurrence for early stage breast cancer ( 2 ). Using different platforms and detection techniques, these panels allow relatively robust prediction of the risk of breast cancer relapse in the short term (0-5 years), but lower prediction power for long-term follow-up ( 3 ). Development of a risk score plays an important role in choosing treatment modality: high-risk patients typically undergo more aggressive chemotherapy, while low-risk patients fare well with hormonal treatment alone ( 4 ). While the treatment regimen for high- and low-risk groups is chosen with a high degree of certainty, the choice of treatment for the patients with intermediate risk is more difficult; heterogeneity of outcome within this group means a certain number of patients will be misclassified ( 5 ). Thus there is urgent need to identify additional molecular markers that can refine the risk prediction for long-term survival. Given that currently approved diagnostic gene panels involve between 21 and 70 genes, it would be beneficial to develop a simpler prognostic panel with a smaller number of genes to better characterize breast cancer mortality risk. MicroRNAs are known as “universal regulators” of gene expression and are involved in post-transcriptional regulation of multiple fundamental cellular processes in cancer pathophysiology including cell proliferation, differentiation, apoptosis and the stemness of the cancer cells. These processes in turn, can be linked to tumor initiation, progression, and metastasis ( 6 ). Because their mechanism of action involves limited complementarity, each miRNA can potentially repress hundreds of mRNA targets. It is now clear that miRNAs are involved in virtually all signaling pathways including those crucial for cell transformation and cancer progression ( 6 ). It is thus not surprising that recent studies have attempted to include miRNAs as potential biomarkers of progression ( 7 ). However, miRNA expression level is not equivalent to its functional activity ( 8 ). Instead, the relative abundance of miRNAs to its target genes determines the functional activity levels of miRNAs ( 9 ). We have recently reported a novel computational approach, ActMiR , for identifying active miRNAs and miRNA-mediated regulatory mechanism and have shown stronger associations between the “activity” of miRNAs with clinical outcomes (including breast cancer survival) compared to the level of expression ( 10 ). For example, we identified that high activities of let-7d and miR-18a were associated with better survival outcomes for estrogen receptor negative human epidermal growth factor receptor 2 negative (ER−/HER2−) breast cancer. As we report in the current study, the same ActMiR algorithm also suggests that high miR-500a activity would be associated with poor survival among ER+ breast cancer patients. Because miR-500a activity cannot be directly measured but only inferred through expression of its associated genes, we designed the translational study reported here to investigate the association of miR-500a mediated gene expression with long-term breast cancer survival.