Regulation of Pharmacogene Expression by microRNA in The Cancer Genome Atlas (TCGA) Research Network.

Individual differences in drug responses are associated with genetic and epigenetic variability of pharmacogene expression. We aimed to identify the relevant miRNAs which regulate pharmacogenes associated with drug responses. The miRNA and mRNA expression profiles derived from data for normal and solid tumor tissues in The Cancer Genome Atlas (TCGA) Research Network. Predicted miRNAs targeted to pharmacogenes were identified using publicly available databases. A total of 95 pharmacogenes were selected from cholangiocarcinoma and colon adenocarcinoma, as well as kidney renal clear cell, liver hepatocellular, and lung squamous cell carcinomas. Through the integration analyses of miRNA and mRNA, 35 miRNAs were found to negatively correlate with mRNA expression levels of 16 pharmacogenes in normal bile duct, liver, colon, and lung tissues (p<0.05). Additionally, 36 miRNAs were related to differential expression of 32 pharmacogene mRNAs in those normal and tumorigenic tissues (p<0.05). These results indicate that changes in expression levels of miRNAs targeted to pharmacogenes in normal and tumor tissues may play a role in determining individual variations in drug response.

INTRODUCTION

Pharmacogenomics focuses on how individual genetic variations influence drug responses, and is helping to develop safer and more effective treatments for patients (Relling and Evans, 2015). Many pharmacogenomic studies have concerned single nucleotide polymorphisms (SNPs) that affect drug responses, and several SNPs have been reported (Georgitsi ). However, the diversity of drug responses is not explained by genetic mutation alone. As well as the genetic polymorphisms, drug response may be different due to factors that regulate gene expression.

MicroRNAs (miRNAs) are small, ∼21 nucleotide single-strand noncoding RNAs that can regulate gene expression by binding to partially complementary sites in 3′ untranslated regions (3′ UTRs) of messenger RNAs (mRNAs). This miRNA-mRNA interaction governs a variety of mechanisms that control gene expression, including mRNA degradation and translational repression (Wienholds ; Pasquinelli, 2012). The mRNAs affected by the miRNAs consequently influence susceptibility to cancer, as well as amentia, autoimmune diseases, and diabetes (Sayed and Abdellatif, 2011). Therefore, miRNAs are becoming recognized as important mediators that affect drug responses, without affecting the genomic sequence. An increasing number of studies on pharmacoepigenetics and pharamcoepigenomics support a role for miRNA in regulating expression of genes encoding proteins involved in drug absorption, distribution, metabolism, and excretion (ADME) (Shomron, 2010; Rukov and Shomron, 2011), as well as pharmacodynamics (Yu ). One miRNA can regulate various ADME genes via direct and/or indirect targeting, or one ADME gene may be modulated by multiple miRNAs (Yu and Pan, 2012). However, our current understanding of miRNA action was mainly obtained from in vitro cell culture systems and ex vivo systems (Rukov and Shomron, 2011). Moreover, prediction and identification of miRNAs target genes is a time-consuming, labor-intensive, and errorprone process (Huang ).

This epigenetic regulation of miRNAs in drug transporters or enzymes has a greater impact on drug responses. The influence of the epigenetic changes in cancer diseases can be expected to be even greater. Thus, we hypothesized that the drug response may be affected by expression changes of pharmacogenes in patients with cancer, in special, in organs involved in drug metabolism. The Cancer Genome Atlas (TCGA) Research Network has profiled and analyzed large numbers of human tumors to discover molecular aberrations at the DNA, RNA, and protein level, and also examined epigenetic changes, including those related to miRNA (Weinstein ). Because the TCGA also contains a significant collection of normal tissue samples, it would be an appropriate resource for pharmacogenomic miRNA studies. Tumor-induced miRNA changes are also important in drug responses and toxicity, especially responses to chemotherapy (Zheng ).

Therefore, the aim of this study was to explore miRNA expression difference in normal tissues derived from patients with five different cancer types and identify significant miRNAs regulating pharmacogene expression, using an integrated analysis of miRNA and mRNA. In addition, we purposed to assess miRNA expression difference, in special, in tumor tissues compared with normal tissue of cancer patient samples.

MATERIALS AND METHODS

miRNA data collection using TCGA datasets

The miRNA data of normal and tumor tissues was downloaded from the TCGA Research Network portal (cancergenome.nih.gov) which dataset was available as of May 2016. All data for cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), kidney renal clear cell (KIRC), and lung squamous cell carcinoma (LUSC) samples were collected in the United States, whereas liver hepatocellular carcinoma (LIHC) samples originated from patients in the United States, France, Japan, and China, considering various ethnic backgrounds. The miRNA sequencing (miRNAseq) data was gathered using an Illumina® HiSeq 2000 platform at the Michael Smith Genome Sciences Centre (GSC) of the BC Cancer Agency (Vancouver, BC, Canada). From the Illumina® HiSeq RNASeqV2 level 3 dataset, the “normalized_count” (quantile normalized relative standard error of the mean) value of each miRNA was collected. The miRNAseq data was integrated in to a matrix with log2 transformed for the downstream analysis.

Pharamcogenes selection and mRNA data collection

Important pharmacogenomic-related genes were searched on the Pharmacogenomics Knowledge Base (Klein ). Additional pharmacogenetic genes, derived from the U.S. Food and Drug Administration (FDA) Table of Pharmacogenomic Biomarkers in Drug Labels (http://www.fda.gov/drugs/scienceresearch/researchareas/pharmacogenetics/ucm083378.html/), were included. The final phamacogenes for analysis were selected by eliminating duplicates.

The public sequencing data of mRNA, associated with selected pharmacogenes, was also collected from the TCGA Research Network portal. RNA sequencing (RNASeq) data were produced by the University of North Carolina (Chapel Hill, NC, USA) using an Illumina® HiSeq 2000 platform. An mRNAseq matrix with log-2 transformation was made for downstream analysis.

Comparison of miRNA expression in normal and tumor tissues

All normal and tumor tissues samples were clustered using a hierarchical method. The clustering dendrogram was drawn using the Ward linkage method. To plot miRNA expression data in a heat map, we selected miRNAs that had >0.1 deviations in expression levels across samples. In addition, a distance matrix for miRNA expression variables in normal tissue samples was constructed using the Euclidean distance and was visualized by multidimensional scaling (MDS). This step was implemented using cmdscale in the R statistics software package.

Correlation analysis of miRNAs and gene expression

We selected only paired data in sold primary tumors and normal tissues to compare the difference in expression of miRNA. We analyzed the correlation between the expression levels of miRNA and mRNA in normal tissues of cancer samples and found a significant negative correlation. In addition, the Pearson’s correlation analysis was performed to identify in tumor specific downregulated miRNA by analyzing the significant association between miRNA and mRNA expression, and correlation coefficients were calculated with adjustment for cancer types.

miRNA target prediction

We next matched the significant correlations with target information using TargetScan (Agarwal ), miRANDA (Betel ), miRDB (Wong and Wang, 2015), Diana Tools (Paraskevopoulou ), miRMap (Vejnar and Zdobnov, 2012), and miRNAMap (Hsu ) as appropriate. Given that no program was consistently superior to the others, and that we aimed to minimize the probability of introducing false positives and/or negatives, we selected genes that were identified by at least three databases as potential targets (Dai and Zhou, 2010). Data extraction and analyses were performed using Python version 3.4 (http://www.python.org/).

Evaluation using GEO dataset

For evaluate with our founding, we collected expression datasets of miRNA and mRNA for tumor and non-tumor tissues derived from colonic adenocarcinoma (GSE29623) (Chen ) and intrahepatic cholangiocarcinoma and hepatocellular carcinoma patients (GSE57555) (Murakami ).

Statistical analysis

Differences between the number of miRNAs and mRNA expression in each cancer patient were analyzed by Student’s t-test. Pairwise comparisons of miRNA expression levels in normal tissues were analyzed with a paired t-test. Regression analysis tested whether changes in miRNA expression correlated with mRNA expression after adjusted by tissue types. All statistical tests were performed in R Statistics version 3.3.2 (http://www.r-project.org/). Statistical significance was defined as a p-value of less than 0.05. Multiple testing correction was performed by controlling the false discovery rate (Benjamini and Hochberg, 1995) at α=0.05.

RESULTS

Pharmacogenes selection

Through searching database, 63 genes were selected and 31 genes were added from FDA table. After adding cytochrome P450 oxidoreductase (POR), a total of 95 genes, including 30 drug-metabolizing enzymes and 12 transporter genes, are listed in Table 1.

Table 1.

List of 95 pharmacogenes in this study

Classification	Gene
Metabolizing enzymes	ADH1A, ADH1B, ADH1C, ALDH1A1, COMT, CYP1A2, CYP2A6, CYP2B6, CYP2C19, CYP2C8, CYP2C9, CYP2D6, CYP2E1, CYP2J2, CYP3A4, CYP3A5, CYP4F2, DPYD, G6PD, GSTP1, GSTT1, NAT1, NAT2, POR, SULT1A1, TPMT, UGT1A1
Transporters	ABCB1, SLC19A1, SLC22A1, SLCO1B1
Targets/Pathway	ABL1, ABL2, ACE, ADRB1, ADRB2, ALK, ALOX5, ASL, ASS1, BCR, BRAF, BRCA1, CFTR, CPS1, CYB5R1, CYB5R2, CYB5R3, CYB5R4, DRD2, EGFR, ERBB2, F2, F5, FIP1L1, HMGCR, HPRT1, IL28B, IL2RA, KCNH2, KCNJ11, KIT, KRAS, LDLR, MS4A1, MTHFR, NAGS, NQO1, NRAS, OTC, P2RY12, P2RY1, PDGFRA, PDGFRB, PGR, PROC, PROS1, PTGIS, PTGS2, SCN5A, SERPINC1, TYMS, VKORC1
gDNA repair	POLG
Transcription factor	AHR, ESR1, NR1I2, PML, RARA, RYR1, VDR
Miscellaneous	HLA-A, HLA-B, HLA-DQA1, HLA-DRB1

A total of 1,448 samples were downloaded from the TCGA portal (36 CHOL, 458 COAD, 244 KIRC, 373 LIHC, and 337 LUSC samples). After excluding unpaired data, 1,870 miRNAs remained in 9 CHOL, 8 COAD, 67 KIRC, 48 LIHC, and 43 LUSC primary tumor and paired normal tissue samples. Through Ward linkage analysis, the samples were clustered into one of four major groups that each represented a human tissue (Fig. 1). The number of mRNAs expressed at lower levels in primary solid tumors was higher than that seen for normal solid tissues (Table 2). Meanwhile, for KIRC and LUSC the number of miRNAs expressed at higher levels in primary solid tumors was lower than that seen for normal solid tissues. The number of miRNAs having lower expression levels in primary tumor tissues was lower than that for normal tissues in patients with CHOL, COAD, and LIHC.

Fig. 1.

Heat map representing miRNA levels of normal tissues derived from colon, kidney, liver, and lung cancer patients. The 55 miRNAs have standard deviations >0.1 across all samples. Each row and column represents a marker and sample, respectively. The clustering dendrogram was drawn using the Ward linkage method.

Table 2.

Comparisons of miRNA and mRNA expression levels between tumors and normal solid tissues derived from cancer patients

Cancer	Number of patients	Number of miRNAs		Number of mRNAs
		Increased in tumors	Decreased in tumors	Increased in tumors	Decreased in tumors
Cholangiocarcinoma	9	120	82	21	48
Colon adenocarinoma	8	255	120	8	26
Kidney renal clear cell	67	182	270	36	45
Liver hepatocellular carcinoma	48	212	168	19	53
Lung squamous cell carcinoma†	43†	157	441	26	42

Significantly differently expressed miRNAs or mRNAs between tumor tissues and normal solid tissues were determined by paired t-test, respectively (p<0.05).

The 36 samples had mRNA expression data.

Based on assessment of miRNA relationships among the 95 pharmacogenes in different tissues, the overall pattern of the MDS plot separated the colon, kidney, liver, and lung into four discrete identities, while bile duct tissues were included with the liver (Fig. 2). A pairwise comparison of miRNA profiles between tissues showed that the profile for normal kidney tissues was closer to that seen for normal colon and lung tissues (Fig. 3). miRNA expression profiles for bile duct tissues were most similar to those for the liver, which differed most significantly from those seen for the kidney. Of the 1,870 miRNAs analyzed, miR-122 exhibited the greatest differences in comparisons between KIRC-LIHC, COAD-LIHC, LIHC-LUSC, CHOL-KIRC, CHOL-LUSC, and CHOL-COAD (p=2.08×10−111, p=6.98×10−79, p=1.03×10−62, p=3.52×10−44, p=3.17×10−43, and p=2.04×10−15, respectively). Similarly, miR-450b, miR-375, miR-590, and miR-26b levels significantly differed among COAD-KIRC (p=1.79×10−60), KIRC-LUSC (p=9.88×10−55), COAD-LUSC (p=1.03×10−32), and CHOL-LIHC (p=2.14×10−11) comparisons.

Fig. 2.

Multidimensional scaling analysis plot of normal tissues based on miRNA distance.

Fig. 3.

Pairwise comparison of miRNA expression levels in normal tissues. CHOL: cholangiocarcinoma, LIHC: liver hepatocellular cell carcinoma, COAD: colon adenocarcinoma, LUSC: lung squamous cell carcinoma, KIRC: kidney renal clear cell carcinoma.

Correlation of miRNA and mRNA expression in normal and tumor tissues

The correlation analysis results showed that 23 miRNAs showed a negative correlation between miRNA and mRNA expression for 14 pharmacogenes (Table 3, Fig. 4), resulting in 33 combinations of miRNAs and mRNAs. HsamiR-429 decreased 3 mRNA expression levels, including ADH1B (p=2.48×10−24), AHR (p=1.63×10−2), and ALDH1A1 (p=1.44×10−3). Meanwhile, hsa-miR-181d decreased the expression levels of AHR (p=2.88×10−3), BCR (p=6.25×10−3), and CYB5R4 (p=6.30×10−3), whereas hsa-miR-152 decreased the expression levels of ABL2 (p=1.48×10−47), AHR (p=7.44×10−8), and CYB5R4 (p=2.24×10−45). Hsa-miR-98 decreased the expression levels of ADRB2 (p=9.13×10−13).

Table 3.

miRNA expression negatively correlated with pharmacogene expression in different normal solid tissues derived from cancer patients (r2>0.3)

Classification	Gene	miRNA	FDR adjusted p-value	Adjusted Pearson correlation coefficient (r2)
Metabolizing enzymes	ADH1B	hsa-miR-429	2.48e−24	0.538
		hsa-miR-577	2.15e−20	0.468
	CYB5R4	hsa-miR-152	2.24e−45	0.812
		hsa-miR-758	1.21e−07	0.437
		hsa-miR-181d	6.30e−03	0.361
Receptors	ADRB1	hsa-miR-let-7c	2.39e−05	0.538
	ADRB2	hsa-miR-98	9.13e−13	0.738
Targets	ABL1	hsa-miR-378g	1.08e−05	0.377
	ABL2	hsa-miR-152	1.48e−47	0.800
		hsa-miR-107	6.30e−12	0.452
		hsa-miR-217	1.96e−09	0.323
		hsa-miR-410	4.19e−04	0.317
	ALOX5	hsa-miR-134	1.93e−02	0.636
Transcription factors	ACE	hsa-miR-511	3.92e−09	0.529
	AHR	hsa-miR-152	7.44e−08	0.597
		hsa-miR-181d	2.88e−03	0.561
		hsa-miR-429	1.63e−24	0.643
		hsa-miR-520b	1.22e−03	0.558
		hsa-miR-653	6.27e−13	0.630

FDR: false discovery rate.

Fig. 4.

Correlation of RNA expression and miRNA changes across normal colon, bile duct, kidney, liver, and lung tissues derived from cancer patients. Line is fitted to the points. Open circle, bile duct; closed circle, kidney; open square, liver; closed square, lung; open triangle, colon (A) correlation of hsa-miR-152 with ABL2 (p=1.48e−47); (B) correlation of hsa-miR-429 with ADH1B (p=2.48e−24); (C) correlation of hsa-miR-98 with ADRB2 (p=9.13e−13).

The correlation analysis results showed that 19 miRNAs had a negative correlation between miRNA and expression levels of 15 pharmacogene mRNAs (Table 4, Fig. 5) to yield 24 combinations between miRNAs and mRNAs. Hsa-miR-520b (1.59×10−3) decreased ADRB1 mRNA expression levels, whereas hsa-miR-152 decreased the expression levels of ABL2 (p=1.49×10−43), AHR (p=4.06×10−19), and CYB5R4 mRNA (p=1.42×10−49). Hsa-miR-98 decreased the expression levels of ADRB2 mRNA (p=1.23×10−41).

Table 4.

miRNA expression negatively correlated with pharmacogene expression in different normal and tumor solid tissues derived from cancer patients (r2>0.3)

Classification	Gene	miRNA	FDR adjusted p-value	Adjusted Pearson correlation coefficient (r2)
Receptors	ADRB1	hsa-miR-520b	1.59e−03	0.450
	ADRB2	hsa-miR-98	1.23e−41	0.450
Targets	ABL2	hsa-miR-152	1.49e−43	0.482
Metabolizing enzymes	CYB5R4	hsa-miR-152	1.42e−49	0.510

FDR, false discovery rate.

Fig. 5.

Correlation of RNA expression and miRNA changes across normal and tumor colon, kidney, liver and lung tissues derived from cancer patients. Line is fitted to the points. Open circle, bile duct; closed circle, kidney; open square, liver; closed square, lung; open triangle, colon (A) correlation of hsa-miR-152 with CYB5R4 (p=1.42e−49); (B) correlation of hsa-miR-98 with ADRB2 (p=1.23e−41) (C) correlation of hsa-miR-152 with ABL2 (p=1.49e−43).

Evaluation using GEO datasets

Through evaluation using GSE29623 and GSE57555 datasets, Hsa-miR-520b decreased mRNA expression of ADRB1, while hsa-miR-98 decreased mRNA expression of ADRB2 (p<0.05). Additionally, hsa-miR-152 decreased mRNA expression levels of ABL2 and CYB5R4 (p<0.05).

DISCUSSION

In the present study, we used the integrative analysis to identify miRNAs that contribute to altered expression of pharmacogenes in different tissues and tumors. The integrative analysis of mRNA and miRNA expressions is a powerful tool for identifying individual genes and genetic or epigenetic mechanisms of gene expression, as well as a means to understand the relationship between target genes and downstream regulation by miRNA (Yang ; Ye ). miRNA and mRNA pharmacogene expression was analyzed in paired normal and tumorigenic samples derived from CHOL, COAD, KIRC, LIHC, and LUSC patients using TCGA data.

The data included 95 pharmacogenes that were selected for analysis in our study. For LIHC, drug-metabolizing enzymes and transporters are abundantly expressed in both the liver and bile duct. The colon, kidneys, and lungs are the main organs involved in the elimination of chemotherapeutic drugs. Since lung and colorectal cancer are the first and second leading causes of cancer-related deaths worldwide, respectively (World Health Organization, 2014), patients with these types of cancer may receive chemotherapy despite the stage-dependence of these drugs.

The United States has announced a research initiative that aims to accelerate progress toward a new era of precision medicine that is tailored to individuals (http://www.whitehouse.gov/precisionmedicine/). Genetic variations and epigenetic changes between individuals may be related to differences in drug responses (Dluzen and Lazarus, 2015). Most previous studies of miRNA in pharmacogenes examined only a limited number of genes with small sample sizes using traditional methods (Rieger ), such that few global miRNA analyses of pharmacogene expression have been performed (Kim ). Our results showed that the number of mRNAs expressed at lower levels in primary solid tumors was higher than that seen for normal solid tissues, while the number of miRNA expression levels of pharmacogenes varied in tumor tissues compared to normal tissues. These results indicate that there are considerable differences in the level and distribution of miRNAs across normal and tumorigenic tissues. However, as expected, our results showed that miRNA and mRNA expression levels were similar between liver and bile duct tissues.

The expression of several drug-metabolizing enzymes and transporter genes was regulated by miRNAs. For example, miR-27a and miR-548a repressed mRNA expression levels of ABCB1 and CYP3A4, respectively (Wei ; Messingerova ). Although we found negative correlations of the expression of these miRNAs and mRNAs in our study, they were excluded because their relationships did not occur in more than three miRNA target prediction databases.

Nevertheless, we could use the integrative analysis of massive miRNA-mRNA expression data to identify new various miRNAs for various drug-metabolizing enzyme (ADH1B, CYB5R4), receptor (ADRB1, ADRB2), target (ABL1, ABL2, ALOX5) genes, and transcription factor (ACE, AHR) that contribute to their differential expression in bile duct, colon, kidney, liver, and lung tissues. Expression of hsa-miR-148 and hsa-miR-152 was reported to be downregulated in gastrointestinal cancer tissues, suggesting that these two miRNAs may be involved during the early stage of gastric carcinogenesis (Chen ). The hsa-miR-520 was also decreased in n colorectal carcinoma when compared with normal colorectal tissues (Bahar ). Associations between these miRNAs and pharmacogenes have not been previously reported. let-7 family members such as let-7, let-7a, let-7b, let-7c, let-7d, let-7e, let-7f, let-7g, let-7i, and miR-98 were previously shown to target ADRB2 (Wang ), but to our knowledge this is the first study to show that hsa-miR-98 can also regulate ADRB2 expression.

Even with targeted therapy, the response to cancer drugs is not solely dependent on tumor epigenetics (Nasr ). Moreover, germ line epigenetics can play a role in drug effects. Therefore, understanding and considering the contribution of both somatic and germ line epigenetics is important when predicting drug response and toxicity.

Recently, there has been a rapid increase in knowledge of how pharmacogenes are regulated by epigenetic mechanisms and methods to analyze this regulation (Koturbash ). Although we examined a limited set of genes known to be involved in drug responses, the methodology described herein can be easily applied to future studies. One limitation of our study is that we did not stratify the data for age, gender, or racial/ethnic backgrounds, although miRNAs have been shown to exhibit differences related to these parameters (Huang ; Kwekel ). miRNAs regulate gene expression by repressing translation and/or by mRNA deadenylation and decay (Djuranovic ). Several groups demonstrated that protein repression can occur in the absence of mRNA degradation (Wilczynska and Bushell, 2015), but we did not analyze protein expression levels of the pharmacogenes targeted in our study. Although there are further challenges to defining the role of miRNA in drug responses, here we identified miRNA-mediated changes in pharmacogene expression that may influence therapeutic responses.

In conclusion, epigenomic changes, including miRNA-induced regulation of expression of genes encoding drug-metabolizing enzymes, transporters, or targets, can potentially lead to changes in drug activity that may contribute to drug sensitivity, resistance, and toxicity. Here we investigated miRNA using publicly available epigenomic and transcriptomic databases in an effort to advance pharmacogenomics research. We believe the current analysis will lead to more rapid identification of functional miRNAs that are relevant to understanding variability in drug responses of cancer patients.