Abstract
Objective: To determine the influence of the single nucleotide polymorphism (SNP) rs4245739 on the binding and expression of microRNAs and subsequent MDM4 expression and the correlation of these factors with clinical determinants of ER-negative breast cancers.
Methods: FindTar and miRanda were used to detect the manner in which potential microRNAs are affected by the SNP rs4245739-flanking sequence. RNA sequencing data for ER-negative breast cancer from The Cancer Genome Atlas (TCGA) were used to compare the expression of miR-184, miR-191, miR-193a, miR-378, and MDM4 in different rs4245739 genotypes.
Results: Comparison of ER-negative cancer patients with and without the expression of miR-191 as well as profile microRNAs (miR-184, miR-191, miR-193a and miR-378 altogether) can differentiate the expression of MDM4 among different rs4245739 genotypes. Although simple genotyping alone did not reveal significant clinical relationships, the combination of genotyping and microRNA profiles was able to significantly differentiate individuals with larger tumor size and lower number of involved lymph nodes (P < 0.05) in the risk group (A allele).
Conclusions: We present two novel methods to analyze SNPs within 3'UTRs that use: (i) a single miRNA marker expression and (ii) an expression profile of miRNAs predicted to bind to the SNP region. We demonstrate that the application of these two methods, in particular the miRNA profile approach, permits detection of new molecular and clinical features related to the rs4245739 variant in ER-negative breast cancer.
keywords
Introduction
Genetic variations, including single nucleotide polymorphisms (SNPs), have emerged as an important resource in the understanding of carcinogenesis1–3. Recently, several studies have reported the role of SNPs in carcinogenesis through their associations with the binding and expression of microRNA (miRNAs)4–6. MiRNAs are short non-coding RNAs (~18–25 nucleotides) that negatively regulate gene expression post-transcriptionally by directly binding to mRNA targets7,8. MiRNAs have been associated with the pathogenesis of almost every human cancer through regulation of several important biological pathways, including cell proliferation, migration, and apoptosis8–11. Unique patterns of miRNA expression have been proposed as clinical biomarkers and potential therapeutic targets in several cancers12–15.
One SNP, rs4245739, has been identified as a risk factor for several cancers, including esophageal cancer16, non-Hodgkin’s lymphoma17, prostate cancer18, lung cancer19, and ovarian cancer20. Rs4245739(A > C) has also been associated with susceptibility to breast cancer21, the most common cancer among women and the fifth leading cause of cancer mortality worldwide22,23. In a study involving a total of more than 10,000 breast cancers negative for estrogen receptor (ER) expression and 75,000 control cases from patients of Caucasian ancestry, pooled results showed that rs4245739 was associated with the risk of ER-negative breast cancer (P < 0.01) but not the risk of ER-positive breast cancer, with an odds ratio (OR) for the A risk allele of 1.14 (95% CI = 1.10–1.18)22. Subsequent studies from China (1,100 breast cancer cases/1,400 controls) and Norway (1,717 breast cancer patients/1,870 controls) have also documented rs4245739 as a risk variant for breast cancer21,24.
Rs4245739(A > C) is located in the three prime untranslated region (3’UTR) of MDM425, an oncogene that negatively regulates p53 expression by directly binding and masking the transactivation domain in response to DNA damage26. This SNP creates a target binding site for the miRNA miR-191-5p, which has been linked to several cancer types18,25. A functional study revealed miRNA-mediated fine-tuning of MDM4 expression as a result of miR-191-5p binding to the C-allele but not the A-allele of rs424573925. Therefore, in A-allele genotypes, MDM4 is overexpressed and associated with a higher risk of breast cancer25. In ovarian cancer, patients with the AA genotype showed significantly increased risk of recurrence and mortality25. Upregulation of miR-191-5p has been reported in breast cancer, although its specific role in carcinogenesis or progression in this cancer type is not yet known27.
In this study, we used a combination of genetic variant data, in silico prediction, and miRNA expression profiles to demonstrate processes by which molecular mechanisms involved in the SNP located at the 3’UTR of the MDM4 gene can alter the expression of several miRNAs associated with specific clinical features, i.e., tumor size and lymph node infiltration in breast cancer.
Materials and methods
Tumor array and sequencing data
Affymetrix Genome-wide SNP 6.0 microarray (Affymetrix Inc., Santa Clara, CA, USA), gene-level miRNA, and gene-level mRNA sequencing data were retrieved for The Cancer Genome Atlas (TCGA) breast cancer cohort28. RNA sequencing data determined with RNA-Seq by Expectation-Maximization (RSEM) software were used. ER-negativity was determined from immunohistochemical records. Altogether, SNP array and expression data were available for 189 primary ER-negative breast cancers. For these cancers, we were able to find corresponding clinical stage, tumor size, and lymph node status.
SNP survey
SNPs listed on dbSNP were extracted, and those with a minor allele frequency (MAF) of less than 0.1 were removed to yield a total of 11,052,394 SNPs29. The genomic location of SNPs was determined with reference to the human genome assembly GRCh37.
SNP genotyping
Genotypes for the 906,600 SNP probes of the Affymetrix SNP 6.0 platform (Affymetrix, Inc.) were determined by first mapping probe set IDs to dbSNP rs IDs. The R package genomewidesnp6Crlmm was then used for genotype calls, with calls having a confidence level of less than 90% filtered out. All 189 samples were retained for rs4245739.
MiRNA-mRNA binding alteration prediction
We used the algorithms FindTar (http://bio.sz.tsinghua.edu.cn/findtar)30 and miRanda (http://www.microrna.org/microrna/home.do) to predict the mRNA targets of miRNAs.
Bioinformatic and statistical analysis
Bioinformatic and statistical analyses were conducted with the R statistical package version 3.1.2. Cut-offs for the different methods were determined as described previously31. P values less than 0.05 were considered statistically significant unless otherwise stated.
Results
Screen of SNPs located at UTRs and miRNA-mRNA binding alterations due to rs4245739
We analyzed the role of SNPs present in the 3’UTRs of the human genome. After discarding those with a MAF of less than 0.1, we obtained 11,052,394 SNPs, the majority of which were inter-genic. Of the SNPs located in genes, more than 100,000 were in 3’UTRs, a quantity almost two and four times the amounts in the coding regions (non-synonymous only) and 5’UTRs, respectively. We next calculated the probability of finding SNPs in these regions calibrated by their total base length in the genome (Figure 1A, left). This probability was higher for 3’UTR regions compared to the coding regions and 5’UTRs, suggesting a selective pressure for this occurrence. We therefore predicted that some of the 3’UTR-located SNPs may alter the binding ability of regulatory miRNAs (Figure 1A, right), which may in turn mediate allele-dependent mRNA post-transcriptional stability. We assessed this phenomenon by analyzing rs4245739, which is present in the 3’UTR of the MDM4 gene and for which the A allele predisposes individuals to ER-negative breast cancer (Figure 1B).22 We used the FindTar330 and miRanda32 algorithms to predict miRNA-mRNA binding alterations between risk and wild type alleles of 3’UTRs produced by rs4245739 (Figure 1C), which showed not only miR-191 but other miRNAs to be affected by this SNP.
MiRNA binding alterations due to SNP rs4245739. (A) Distribution of SNPs located at the 3’UTR and 5’UTR per 100,000 bases in comparison to coding regions (NS) in the human genome (left panel). Schematic description of how SNP located at the 3’UTR might influence binding sites of miRNA (right panel). (B) Schematic representation of SNP rs4245739, which is located at the 3’UTR of the MDM4 gene. (C) Potential binding alterations of several miRNAs due to rs4245739 that affect not only miR-191 but also 16 additional miRNAs.
Expression of miRNAs affected by SNP rs4245739 genotypes in patients with ER-negative breast cancer
We used an ER-negative breast cancer data set retrieved from the TCGA to study the manner in which the potential alterations in miRNA binding result in phenotypic changes. First, we obtained the rs4245739 genotype of individuals from this data set and observed the homozygotic genotype (CC) and heterozygotic genotype (AC) in approximately 20% of primary tumors (Figure 2A). Interestingly, there were very few genotypic changes between germline and tumor tissue, with only conversion from the heterozygote to either the homozygote C-allele or homozygote A-allele being observed. We next focused on the expression of miRNAs with binding to MDM4 that could be altered in an allele-dependent manner (Figure 2B). Only miR-184, miR-191, miR-193a, and miR-378 had an appreciable average expression, albeit with a large variance. Notably, two of these miRNAs (miR-193a and miR-378) were predicted to bind to the A-allele MDM4 3’UTR only, whereas the other two (miR-184 and miR-191) were predicted to bind to it only in the presence of the C-allele. Examination of whether the expression of these miRNAs differed between genotypes did not yield significant results (Figure 2C).
Expression of miRNAs according to different rs4245739 genotypes in ER-negative breast cancer. (A) Distribution of ER-negative breast cancer patients according to rs4245739 genotypes. (B) Expression of 17 miRNAs with binding to MDM4 mRNAs that is altered owing to rs4245739. (C) Expression of miR-184, miR-191, miR-193a, and miR-378 in different rs4245739 genotypes.
Characterization of rs4245739-affected miRNA expression and 3’UTR SNP functions
Reasoning that the expression level of these miRNAs may complement our ability to understand 3’UTR SNP function and its effect on disease, we approached the characterization of rs4245739 in three different ways (Figure 3A). In the first approach, we simply grouped tumors according to rs4245739 genotype (Figure 3A, left). In the second method, we additionally stratified patients according to expression or lack of expression of miRNA-191 (Figure 3A, center), which would enable us to assert with more confidence whether allele-dependent miRNA modulation of MDM4 was occurring. Our final approach for classifying samples consisted of generating an miRNA profile of each individual (Figure 3A, right). In the case of rs4245739, we retrieved from the TCGA data sets the expression levels of the four miRNAs predicted to affect MDM4 expression (miR-184, miR-191, miR-193a, and miR-378). Given that miR-193a and miR-378 would bind to the A-allele-containing 3’UTR, and miR-184 and miR-191 would bind to the C-allele-containing 3’UTR, we divided the individuals into an A-allele SNP+/miR-193a-/miR-378- group and a C-allele SNP+/miR-184+/miR-191+ group. These two groups would represent A-genotype individuals without miRNA-mediated MDM4 alterations, and C-genotype individuals with miRNA-mediated MDM4 alterations. We compared these two novel approaches to 3’UTR SNP analysis with the genotype-only classification by performing differential MDM4 gene expression analysis between the groups with each approach. In order to categorize tumors into miR-191-expressing and non-miR-191-expressing cases, we compared the results obtained with different thresholds against the results from the simple genotyping approach (Figure 3B, left). We did not observe significant changes in MDM4 expression among any of the groups generated by simple genotyping or the miR-191-expressing approach. Nevertheless, when we applied the miRNA-profile approach with different miRNA expression thresholds (differences in Figure 3B, right, and P values in Figure S1A, left), we observed significant changes in MDM4 expression. Unexpectedly, we found overexpression of MDM4 in the case of miRNA binding to its 3’UTR. In light of these results, we analyzed the expression of other genes located near (< 2Mb) the rs4245739 SNP using the miRNA profile approach. Three other genes (MFSD4, SLC26A9, and SLC41A1) were significantly repressed in the miRNA-bound groups (Figure S1A, right).
Characterization of rs4245739-affected miRNA expression using several methods to show 3’UTR function. (A) Three approaches for characterization of rs4245739 at the 3’UTR of the MDM4: grouping patients according to rs4245739 genotype only, stratifying patients according to presence or absence of miR-191 expression, and classifying patients according to miRNA expression profile. Round object represents A-allele, square object represents C-allele, red color represents MDM4 overexpression, green color represents MDM4 underexpression. Left panel: using rs4245739 genotype only, we could not differentiate tumors with over- or underexpression of MDM4. Previous research demonstrated that A-allele was associated with MDM4 overexpression and C-allele was associated with miR-191 binding and MDM4 underexpression25. Middle panel: using additional determinant of presence or absence of miR-191 expression, we could better discriminate patients with A- or C-allele and status of MDM4 expression. Right panel: using determinant of microRNA profile, we could precicely discriminate patients with A-allele with high MDM4 expression and patients with C-allele with low MDM4 expression. (B) Comparison of MDM4 expression between groups of specific 3’UTR genotypes after stratifying patients using presence of miR-191 expression (left panel) and miRNA expression profile (right panel). Left panel: as miR-191 binds only to C-allele to mediate MDM4 degradation, we divided individuals into A- and C-allele groups to compare MDM4 expression (X-axis) according to miR-191 expression (Y-axis). We observed inverse MDM4 expression at a certain miR-191 threshold that was not statistically significant. Right panel: we compared MDM4 expression (X-axis) in individuals with A- and C-alleles according to miR profile (Y axis) in which miR-193a and miR-378 would bind to the A-allele- and miR-184 and miR-191 would bind to the C-allele; therefore, we divided the individuals into an A-allele SNP+/miR-193a-/miR-378- group and a C-allele SNP+/miR-184+/miR-191+ group. We observed significant inverse MDM4 expression compared to miR profile expression. (C) We categorized individuals into groups with A-allele/presence of miR-191 expression and C-allele/absence of miR-191 expression (left panel) as well as an A-allele SNP+/miR-193a-/miR-378- group and a C-allele SNP+/miR-184+/miR-191+ group (right panel) and then compared with the tumor size (P value, X-axis) across miR expression thresholds (Y-axis). We demonstrated significantly different tumor size at higher thresholds, indicating that A-risk allele correlated with larger tumor size at higher miR expression threshold. (D) We categorized individuals into A-allele/presence of miR-191 expression and C-allele/absence of miR-191 expression groups (left panel), as well as an A-allele SNP+/miR-193a-/miR-378- group and a C-allele SNP+/miR-184+/miR-191+ group (right panel), and then compared with the number of lymph nodes infiltrated (P value, X-axis) across miR expression thresholds (Y-axis). We demonstrated significant differences in the number of lymph node infiltrated miR profile but not miR-191 expression, indicating that A-risk allele correlated with higher number of lymph nodes infiltrated using miR profile.
Activated biological pathways in patients with ER-negative cancer with different rs4245739 genotypes and miRNA expression profiles
We next analyzed the full transcriptome to detect whether gene expression differences were conserved among the different methods. We observed 100% overlap between the sets of differentially expressed genes identified by these two techniques up to a threshold of 300 units of miR-191 expression, at which point the overlap rapidly fell below 50% (Figure S1B, S1C, left). At the same time, the number of genes detected by the marker method remained relatively stable across different expression thresholds (Figure S1B, S1C, right). Finally, when we analyzed the results obtained with the third method, which combined the expression levels of miRNA predicted to bind to a given allele to give an integrated score, we found that the differentially expressed gene overlap between the combined approach and simple genotyping approaches rapidly decreased (Figure 1C, left). As with the previous method comparison, differences were not related to the number of genes obtained with each approach (Figure 1C, right). To analyze the biological importance of the genes obtained, we used the 95th percentile of the entire miRNA expression data set in the case of the miR-191 marker. As a result, we obtained two groups: one comprising individuals with the C-genotype and miR-191 expression above the aforementioned cut-off (n = 37), and another comprising all other individuals (n = 157). In the case of the miRNA profile method, we restricted our samples to the 70% of individuals following the miRNA profile results. The two groups comprised 75 individuals with the A-genotype, and 45 individuals with the C-allele. We performed pathway analysis of these results using the DAVID webserver (Tables 1-3). We found processes such as “immune response” and “cell adhesion” to be overrepresented in all three methods, but each method also detected specific pathways. For example, a protocadherin cluster was found only when using the miRNA profile.
Clinical relevance of characterization of SNP rs4245739 genotypes and miRNA profiles
Finally, we assessed the strength of association between the results of these new methods and clinical parameters. The use of simple genotyping alone failed to reveal any significant relationships (Figure 3C and D and Figure S2, P values marked by dotted lines). When we examined tumor size categories of samples, we found significant differences among groups using both methods at higher thresholds (P values in Figure 3C, differences in Figure S2A). In addition, the number of lymph nodes affected in individuals was found to be significantly reduced in the C-allele group obtained by the miRNA profile method (P values in Figure 3D, right; differences in Supplementary Figure 2B, right), using several different thresholds. Results with the marker method, similar to those using the simple genotyping approach, did not reach significance in this association (Figure 3D, left; Supplementary Figure 2B, left). Associations with other clinical parameters, including breast cancer mortality (Supplementary Figure 2C) and clinical stage (Supplementary Figure 2D), were not statistically significant. Our novel methods thus improve the prediction of these associations compared to mere use of the presence or absence of a SNP variant.
DAVID results for rs4245739 analysis using method 1 (only SNP presence/absence)
DAVID results for rs4245739 analysis using method 2 (SNP and miR-191 expression)
DAVID results for rs4245739 analysis using method 3 (SNP and miR profile)
Discussion
We demonstrated the use of a single miRNA expression marker and expression profile of miRNAs associated with rs4245739, a SNP located at the 3’UTR of the MDM4 gene, to identify tumor features linked to this variant. Our methods, in particular, use of the expression profile, were superior compared to classification based on a single SNP in predicting disease characteristics and outcomes among patients with ER-negative breast cancer. In addition, we found overexpression of MDM4 in the case of miRNA binding to its 3’UTR. Previous research has shown that rs4245739 affects binding of miR-191 and regulates expression of MDM4 mRNA in ovarian cancer25. These findings demonstrate the complexity of the processes governing miRNA-mRNA interactions18,33. We combined SNP genotyping data with miRNA as well as MDM4 expression in ER-negative breast cancer. Single miRNA expression in the presence of the SNP variant did not influence MDM4 expression (Figure 3B, left). However, using the miRNA profile approach, we were able to observe differential changes in MDM4 expression (Figure 3B, right; Figure S1A, left). Our results confirm the findings of a previous study of ovarian cancer by Wynandaele and colleagues25.
More importantly, by combining SNP variant presence and the miRNA profile approach, we were able to differentiate subgroups of patients at risk with greater tumor size and number of lymph nodes infiltrated, indicating the potential for the use of SNPs in this way as diagnostic and prognostic markers in ER-negative breast cancer. Tumor size is used together with nodal and metastasis status to determine clinical breast cancer staging and suitable breast cancer management. In addition, around 30% of breast cancers are ER-negative tumors, with a higher proportion in young breast cancer patients, and these tumors constitute a subgroup with relatively poorer prognosis and fewer treatment options34–36. Our enrichment results showed a common theme of immune response across the three different methods; these findings that are consistent with previous results highlighting the importance of immune pathways in ER-negative breast cancer37–39.
We found differential MDM4 expression using microRNA profiles in relation to different genotypes. Our findings support the previous notion that SNPs located at 3’UTRs may lead to dysregulation of target mRNAs and are potentially implicated in the development of cancer and other chronic diseases25. Indeed, many studies have shown that SNP variants in the 3’UTR elevate the risk of developing several cancer types. Saetrom and associates40 showed that rs1434536 located at the bone morphogenetic receptor type 1B gene (BMPR1B) 3’UTR influences the binding of miR-125b and correlates with breast cancer risk. Located at the estrogen receptor-1 gene (ESR1) 3’UTR, rs9341070 is a risk variant for breast cancer. The SNP affects interaction with miR-206 and influences the expression of ESR141. As a risk variant for squamous cell carcinoma of the head and neck, SNP rs8126 (T > C) is located at the TNFAIP2 3’UTR and influences the expression of TNFAIP242. In addition, rs3134615, which is a susceptibility variant for small-cell lung cancer, is located at the MYCL1 3’UTR and influences the affinity for miR-182743. Rs34764978, which is mapped at the 3’UTR of human dihydrofolate reductase (DHFR), affects the binding and expression of miR-24 and correlates significantly with resistance to the chemotherapeutic agent methotrexate44. These data show that SNPs of specific genes located at the 3’UTR play an important role in the pathogenesis of several diseases through modulating the binding of particular miRNAs.
MiR-19125 as well as miR-184, miR-193, and miR378 interact and modulate MDM4 expression to certain rs4245739 genotypes. MDM4 is an oncogene that inhibits p53 tumor suppressor proteins during regulation of cell proliferation and response to DNA damage25. Nagpal et al27. have reported that miR-191 is an estrogen-sensitive microRNA that acts in regulation of cell proliferation through inhibition of BDNF, CDK6, and SATB1. MiR-184 regulates the PI3K/AKT/mTOR pathway during mammary development and tumor progression45 as well as p53/p21 activity during the cell cycle and cell migration46. MiR-193a inhibits cell growth and migration in breast cancer cells by directly targeting NLN, CCND1, PLAU, and SEPN147. MiR-378 induces a metabolic shift from the glycolytic to the oxidative pathway through PGC-1β/ERRγ48, and modulates therapeutic responses to tamoxifen by targeting GOLT1A49. In addition, gene enrichment pathway analysis revealed that biological processes such as “immune response” and “cell adhesion” were overrepresented in all three methods we used, suggesting that these microRNAs play roles in immune response and cell migration. As its binding affinity for the MDM4 3’UTR is affected by SNP rs4245739, miR-92b is an important candidate for further functional studies because the precursor sequence forms a stem loop for miR-17-92 cluster, which plays a significant role in the carcinogenesis of many cancers, including breast cancer50. MiR-17-92 cluster is known as oncomir-150,51. Further investigations could involve experiments to confirm functional rs4245739, miR-92b, and MDM4 or TP53 expression with the use of luciferase reporter assays and overexpression of miR-92 mimics in combination with reverse transcription-PCR in cell lines.
Acknowledgements
Wahyu Wulaningsih was supported by grants from UK Medical Research Councils (Grant No. MC_UU_12019/2 and MC_UU_12019/4). Sumadi Lukman Anwar received a grant from PTUPT (Grant No. Ristekdikti 09_18).
Footnotes
Conflict of interest statement No potential conflicts of interest are disclosed.
- Received April 28, 2017.
- Accepted August 16, 2017.
- Copyright: © 2017, Cancer Biology & Medicine
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY) 4.0, which permits unrestricted use, distribution and reproduction in any medium, provided the original author and source are credited.
