MicroRNA deregulation and pathway alterations in nasopharyngeal carcinoma

1MicroRNAs (miRNAs) are short non-coding RNA molecules involved in post-transcriptional gene regulation. In animals, miRNAs control the expression of target genes by inhibiting translation or degradating target mRNAs through binding to their 3′UTR. MicroRNAs have been found to regulate genes involved in diverse biological functions, including development, differentiation, proliferation, and apoptosis.

2Cumulative evidence suggests that miRNAs have major functions in the pathogenesis of tumour. Approximately 50% of miRNAs are localised in cancer-associated genome regions. Biological characterisation also identified several miRNAs function as tumour suppressors or oncogenes. Large scale profiling revealed a global alteration of miRNA expression patterns in human cancers. Recently, distinct miRNA expression signatures have been proposed as diagnostic and prognostic markers for various types of human cancer.

3Although the biological effects of many miRNAs have been characterised individually, the impact of multiple miRNA dysregulation on cellular functions and their roles in tumour progression remain largely unknown. Proteomic and microarray data reveal that although each miRNA may regulate up to hundreds of genes, their effect on individual gene is moderate at best. Recent studies suggest that multiple miRNAs may work in concert to regulate related targets in a common pathway. Therefore, pathway analysis, rather than individual target gene characterisation, may provide a better solution to evaluate the biological consequences of global miRNA dysregulation.

4To link the miRNA profiling data with biological functions, has developed a strategy, which calculates the functions and pathways collectively regulated by the coexpressed miRNA on the basis of computationally predicted targets. As all target prediction algorithms generate certain fraction of false positives, these false positives may significantly reduce the data reliability when targets predicted for multiple miRNAs were combined to calculate the functional pathways. Recent microarray and proteomic data have provided valuable insights on target prediction. With these refinements implemented in miRNA target prediction, the results of pathway enrichment analysis on the basis of the coregulated targets should provide a good insight into the functional role of dysregulated miRNAs.

5In this study, we used both experimental and computational approaches to assess the functional impact of miRNA dysregulation in nasopharyngeal carcinoma (NPC). Differentially expressed miRNAs were identified by profiling 270 human miRNAs in clinical samples. Target prediction followed by pathway enrichment analysis was conducted to identify the functional pathways specifically regulated by the down-modulated miRNAs. Two modifications were introduced to increase the confidence for target prediction and enhance the specificity for pathway analysis. The identified specific pathways were validated by the expression data and in good agreement with earlier reported pathogenesis in NPC.

6Samples of NPC tissue and adjacent normal nasopharynx tissue were obtained from patients undergoing surgery and were frozen immediately after surgical resection. Collection and distribution of tissue specimens were performed in accordance with the Institutional Regulation Board of Chang Gung Memorial Hospital, Taiwan. Total RNA was prepared using TRIzol reagent (InVitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The concentration of RNA was quantified using a NanoDrop Spectrophotometer. The RNA integrity was evaluated by Agilent 2100 BioAnalyzer (Agilent Technologies, Palo Alto, CA, USA). RNAs with an RNA integrity number (RIN) >7.5 were used for miRNA and mRNA quantifications.

7For miRNA quantification, a pulsed RT reaction, as described by, was performed to convert all miRNAs into corresponding cDNAs in one RT reaction. Briefly, 10-μl reaction mixture containing miRNA-specific stem-loop RT primers (final 2 nM each), 500 μM dNTP, 0.5 μl Superscript III (InVitrogen), and 1 μg total RNA were used for the RT reaction. The pulsed RT reaction was performed as follows: 16°C for 30 min, followed by 50 cycles at 20°C for 30 s, 42°C for 30 s, and 50°C for 1 s. To quantify mRNA transcripts, 1 μg of total RNA was converted into cDNA using a random hexamer as primer in a 10-μl reaction with the following conditions: 16°C for 30 min, followed by 42°C for 90 min. Reverse transcription products were diluted 20-folds before using for miRNA and mRNA Q-PCR reactions.

8For miRNA quantification, 1 μl of diluted RT product was used as template for a 10 μl PCR. Briefly, 1 × SYBR Master Mix (Applied Biosystem, Foster City, CA, USA), 200 nM miRNA-specific forward primer, and 200 nM universal reverse primer were used for each PCR reaction. The condition for Q-PCR is 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 63°C for 32 s, and a dissociation stage. Endpoint reaction products were analysed on a 10% polyacrylamide gel stained with ethidium bromide to discriminate the correct amplification product (57–60 bp) and the potential primer dimmers (<44 bp). For mRNA quantification, the following PCR conditions were used: 95°C for 10 min, followed by 45 cycles of 95°C for 15 s and 60°C for 1 min, and a dissociation stage. All Q-PCR reactions were performed on an ABI Prism 7500 Fast Real-Time PCR system (Foster City, CA, USA).

9The threshold cycle (Ct) is defined as the cycle number at which the change of fluorescence intensity crosses the threshold of 0.2. Of the 270 miRNAs evaluated, 47 showed the expression level below the detection limit (Ct>35) in more than 70% of samples, and were excluded from the analysis. For the remaining 223 miRNAs, raw Ct data were converted to 39 – Ct normalised by global median normalisation before further analysis. For mRNA expression, the average Ct of β-2-microglobulin, actin, and GAPDH was subtracted from the raw Ct value to obtain δ-Ct (dCt). The experimentally normalised dCt values were converted to 39 – Ct used to analyse the expression level of human mRNA transcripts.

10Statistical analyses used for miRNA and mRNA expression data including two-sample t-tests (two-tailed), paired-sample t-test (two-tailed), Mann–Whitney test, principle component analysis, Pearson's correlation analysis, and hierarchical clustering were performed with Partek Genomics Suite (version 6.3, St Louis, MO, USA). Pathway enrichment analysis was performed using the MetaCore database (GeneGO, St Joseph, MI, USA). P-values for pathway enrichment analysis were calculated using the formula for hypergeometric distribution and reflects the probability for a pathway to arise by chance.

11The roles of miRNA dysregulation in the pathogenesis of NPC have not been well explored. To identify miRNAs differentially expressed in NPC, we analysed the expression levels of 270 human miRNAs in 13 NPC and 9 normal nasopharyngeal tissues. These samples include seven pairs of NPC and their adjacent normal tissues from same patients. MicroRNA profiling was performed using a stem-loop RT–qPCR method as described by. The stem-loop RT–qPCR assay has been proven to offer high sensitivity and specificity for the quantification of mature miRNAs. The expression data of 47 miRNAs were eliminated from further analysis because of their low abundance in all samples tested (Ct>35). Expression levels of the remaining 223 miRNAs were expressed as 39 – Ct after global median normalisation.

12Unsupervised hierarchical clustering analysis using expression levels of all 223 detectable miRNAs in the seven paired normal–NPC tissues generated a tree with normal and NPC samples clearly separated into two groups (Figure 1A). To identify differentially expressed miRNAs in normal and NPC tissues, two statistical tests (two-sample t-test and Mann–Whitney rank test) were performed. Expression of 35 miRNAs was significantly altered in NPC samples (fold change: ⩾3; false discovery rate: <0.05; Figure 1B). Among them, 11 miRNAs, including miR-196b, miR-138, miR-155, miR-142-3p, and miR-18a, were significantly upmodulated and 24 miRNAs, including miR-204, miR-449a, miR-34c-3p, miR-143, and miR-145, were down-modulated in NPC samples. Complete list of differentially expressed miRNAs and their chromosome location are shown in Table 1. Similar results were obtained when the miRNA expression levels were normalised using the geometric mean of two reference miRNAs, miR-103 and miR-191, as suggested by; Supplementary Figure S1 and Supplementary Table S1).

13To test whether the differentially expressed miRNAs selected from paired samples show similar trend of alteration in unpaired samples, we included two additional normal and six additional NPC samples for the principal component analysis (PCA). As shown in Figure 1C, the PCA analysis showed a complete segregation between the 9 normal and 13 NPC samples. Unsupervised hierarchical clustering with the 35 altered miRNAs showed a clear separation between NPC and normal samples. Expression levels of four most significantly upregulated miRNAs and four most significantly down-modulated miRNAs in NPC tissues were shown in Figure 2A and B. The global alteration in miRNA expression pattern observed in our study is similar to earlier reports on other human cancers.

14The profiling analysis indicated that a large number of miRNAs are down-modulated in NPC tissues. As miRNAs are negative regulators of protein-coding genes, down-modulation of these miRNAs are expected to cause an upregulation of their target genes and alterations of the associated cellular pathways in NPC tissues. To estimate the overall effect of these down-modulated miRNAs on cellular functions, we adopted a two-stage approach as depicted in Figure 3A. The first stage was designed to identify target genes coregulated by these down-modulated miRNAs. The second stage performed a pathway enrichment analysis using the coregulated targets to identify cellular functions specifically regulated by these down-modulated miRNAs.

15Most computational algorithms predict miRNA targets on the basis of sequence complementarity and/or thermostability. However, recent studies indicated that many additional factors, such as local AU content and site position, can significantly affect the target site efficacy. The overall target efficacy can be expressed by a context score. The correlation between context score and site efficacy has been validated in both mRNA and protein levels. Therefore, the concept of context score was implemented to select high confidence targets for the down-modulated miRNAs. The context score threshold for high efficacy target was set at −0.2.

16To identify coregulated targets for the down-modulated miRNAs, we retrieved all target genes listed in the TargetScan database (http://www.targetscan.org/; ), which includes targets for 456 miRNA families predicted by seed sequence complementarity. Low confidence targets were eliminated by filtering out targets with total context score greater than −0.2. Two miRNAs, miR-34c-3p, and miR-30a* were not present in the TargetScan database, and therefore were excluded from further analysis. The 22 down-modulated miRNAs (‘down-miR’) were used for target prediction and pathway analysis. To enhance the specificity of the target and pathway analysis, a control set of 22 miRNAs (‘ctrl-miR’) was included in the target prediction and pathway analysis. These 22 miRNAs showed constant expression level in normal and NPC tissues and shared no sequence homology in their seed regions with the 22 down-modulated miRNAs. The seed sequence and the number of unique target predicted for individual miRNA were listed in Table 2.

17The number of total targets predicted for down-miR and ctrl-miR was 8014 and 8039, respectively. Distribution of predicted targets vs coregulating miRNAs for down-miR and ctrl-miR was shown in Figure 3B. Within each target group, approximately half of the predicted targets were recognised by single miRNA (44% of down-miR targets and 50% of the ctrl-miR targets). A majority of targets were recognised by three or fewer miRNAs (85% of down-miR targets and 90% of ctrl-miR targets). Between the two target groups, more than 55% of predicted targets were shared by both groups when the number of coregulating miRNA was equal or less than three. However, when the number of coregulating miRNA was equal or more than four, the fraction of overlapping targets decreased to less than 20%. These results suggested that genes coregulated by four or more miRNAs provided a better distinction between the two target groups and were more suitable for pathway-enrichment analysis.

18The number of targets regulated by four or more miRNAs was 1223 and 838 for down-miR and ctrl-miR, respectively. These two target groups were selected and uploaded into MetaCore for pathway enrichment analysis. Statistically enriched pathways were identified using a threshold P-value of 0.001. Sixty-five pathways were found to be statistically enriched with the down-miR targets, whereas only four pathways were enriched with the ctrl-miR targets. For each pathway, we calculated the ratio of P-values between down-miR and ctrl-miR to determine the specificity. With a P-value ratio >1000 as the criteria, we identified 45 pathways specifically enriched by down-miR targets. These specifically targeted pathways include many pathways involved in tumour pathogenesis, such as TGF-Wnt pathways, G1-S cell cycle progression, VEGF signalling pathway, apoptosis and survival pathways, and IP3 signalling pathways. The five most specifically targeted pathways were shown in Figure 3D. Table 3 listed the specifically targeted pathways grouped by their cellular functions.

19Regulation of G1/S transition and cross-talk between VEGF and angiopoietin-1 signalling are the two most significantly pathways targeted by down-miR. To validate the effect of down-miR on these two pathways in NPC, we selected three predicted target genes from each pathway and compared their expression levels in normal and NPC samples (Table 4). Real-time RT–PCR analysis revealed that the expression levels of all six predicted targets, including cyclin D2 (CCND2), cyclin E2 (CCNE2), CDC25A, VEGFA, phospholipase C-γ1 (PLCG1), and AKT3, were significantly elevated in NPC tissues (Figure 4A and B). We also examined the protein levels of CCNE2 in NPC tissue sections by immunohistochemical staining to identify which cell type in the tumour mass express cyclin E2 protein. Two representative cases of cyclin E2 staining are shown in Figure 4C. Strong positive cyclin E2 staining was detected in tumour cells (left panels), but not in paired normal nasopharyngeal epithelium (right panels). This result confirmed that the predicted targets of down-miR are upregulated in NPC.

20The 3′UTR of cyclin E2 contains binding sites for many different miRNAs. Six of them (miR-34c-5p, miR-449, miR-200a, miR-200b, miR-195, and miR-9) were down-modulated in NPC (Figure 5A and B). We further analysed the expression correlation between three of the down-modulated miRNAs and cyclin E2 in 16 samples, including 10 NPC and 6 normal nasopharyngeal tissues. Expression of all three miRNAs showed significant inverse correlation (Pearson's correlation coefficient <−0.5, P<0.05) with the expression level of cyclin E2 (Figure 5C). These results indicated that the down-modulated miRNAs cooperatively upregulated critical targets in the G1/S transition pathway. To determine whether the expression level of CCNE2 is indeed regulated by down-modulated miRNAs, we examined the mRNA and protein level of CCNE2 in HK1 cells overexpressing one of the coregulating miRNA, miR-9 (Figure 5D). As shown in Figure 5E and F, HK1/miR-9 cells showed a significant decrease in CCNE2 mRNA and protein levels, indicating that CCNE2 is a direct target of miR-9.

21This study used a stem-loop RT–PCR method to quantify the expression levels of 270 miRNAs in NPC tissues. Using a three-fold change as the cutoff, we identified 35 miRNAs whose expression levels were significantly altered in NPC tissues. Unsupervised hierarchical clustering using these differentially expressed miRNAs clearly segregated NPC tissues from the normal samples. Among the differentially expressed miRNAs identified in this study, many are known oncogenic or tumour-suppressive miRNAs previously shown to be dysregulated in other cancer types. Our results provided further support for using miRNA expression patterns as diagnostic markers.

22The profiling analysis identified a group of miRNAs whose expression was significantly down-modulated in NPC. To assess the global impact of these down-modulated miRNAs, we conducted a high stringency target prediction to identify potential target genes, which were coregulated by multiple down-modulated miRNAs, and then analysed the cellular pathways, which were specifically enriched with these coregulated genes. Results from the pathway enrichment analysis suggested that these down-modulated miRNAs selectively target signalling cascades involved in cell cycle regulation, cell survival and apoptosis, and cytoskeleton remodelling. Quantification of target genes in clinical samples confirmed the elevated expression of several coregulated targets, and revealed a significant inverse correlation between the down-modulated miRNAs and their targets. These results indicated that the target prediction combined with pathway enrichment analysis is an effective approach in investigating the influence of multiple miRNAs.

23Recently, identified eight miRNAs differentially expressed in NPC. Three of the most significantly altered miRNAs, including miR-29c, miR-34b, and miR-34c, also showed significantly reduced expression in our NPC samples. However, the remaining five miRNAs identified in their study did not show significant alteration in the current study. The discrepancy could be due to the difference in sample preparation or difference between microarray and RT–PCR assay.

24Expression of many oncogenic miRNAs, including miR-155, miR-106a, and three miRNAs from the miR-17-92 cluster, were found significantly elevated in NPC. The miR-155 upregulates the NF-κB signalling pathway and is highly expressed in many haematological malignancies and solid tumours. MicroRNAs from the miR-17-92 cluster suppress multiple proliferation inhibitors as well as proapoptotic targets, and promote tumour formation. The upregulation of the miR-17-92 cluster is noted in leukaemia, thyroid cancer, and lung cancer.

25The down-modulated mRNAs include miR-34b/miR-34c cluster, the miR-143/miR-145 cluster, and the miR-195/miR-497 cluster. Both miR-34b and miR-34c are downstream components of the p53 tumour suppressor network. Induction of miR-34b and miR-34c leads to apoptosis or cellular senescence, whereas reduced miR-34b/c expression attenuates p53-mediated cell death. Reduced expression of miR-143 and miR-145 has been documented in cervical cancer, colon cancer, and B-cell malignancies. Forced expression of miR-143 and miR-145 inhibit cell growth in cervical cancer cells. Although the biological function of miR-195 and miR-497 has not been characterised, reduced expression of miR-195 and miR-497 on the 17p13 locus has been documented in primary peritoneal carcinoma.

26In addition to the commonly observed miRNAs, this study also identified several dysregulated miRNAs whose aberrant expression has not been reported before. For example, miR-25* was significantly elevated in NPC and miR-532, and miR-199b-5p were down-modulated. MiR-449 shared the same seed sequence with miR-34c and was found down-modulated in our study. Similarly, miR-152 shared the same seed sequence with miR-148a, and both miRNAs were down-modulated in NPC.

27Recently, two studies using microarray and proteomic approach clearly showed that a single miRNA can regulate the expression of hundreds of target genes. The widespread repression can be detected in both mRNA and protein levels, but the degree of repression for individual target gene is usually very mild. On the other hand, miRNA profiling analysis typically uncovered aberrant expression of multiple miRNAs in the diseased tissues. How to assess the biological consequence of miRNA dysregulation has become a major challenge. To overcome this difficulty, used a computational approach to calculate the global pattern of cellular functions and pathways affected by the coexpressed miRNAs.

28In this study, we adopted the systems biology-based approach to investigate the global impact of differentially expressed miRNAs in the pathogenesis of NPC. Similar to the strategy described by, we analysed the major signalling pathways that are most likely to be affected by the coexpressed miRNAs. Two modifications were introduced to enhance our confidence in target prediction and pathway analysis.

29To enhance the accuracy for miRNA target prediction, we implemented the ‘context score’ concept during the target prediction process. Only genes with a total context score less than −0.2 were included in our target pool. The concept of context score was first proposed by to evaluate the site efficacy of miRNA target using microarray data. The ‘context score’ combined the complementarity of seed sequence, local AU content, site position, and additional base pairing at the 3′ site of miRNA to evaluate the site efficacy between miRNAs and their targets. This concept was further supported by the data from recent proteomic studies.

30To enhance the specificity of pathway enrichment analysis, we selected a set of control miRNAs, which were not altered in NPC, and conducted the target prediction and pathway enrichment analysis in parallel. The probability of a gene targeted by multiple miRNAs increases as the length of its 3′UTR increases. This potential bias introduced by genes with longer 3′UTR can be minimised by comparing it with the targets coregulated by the control miRNAs. The pathway specificity was calculated by comparing the P-values between these two target sets. This strategy effectively addressed the complexity that arise from the coexpressed miRNAs and coregulated targets.

31Our studies revealed that the targets coregulated by the down-modulated miRNA are specifically enriched with genes involved in ‘TGF-WNT and cytoskeletal remodelling’ pathways. Predicted targets include several key regulators of the Wnt pathway, such as WNT1, WNT2, frizzled homologue 3, 7, and 8, AKT, and NF-κB1, as well as multiple extracellular matrix proteins, such as collagen 4A1, 4A4, and fibronectin 1. Studies by, using profiling and immunohistochemical stains, shows that the protein levels of these coregulated targets were significantly elevated in NPC tissues, whereas elevated expression of extracellular matrix proteins was shown by using Q-PCR. The analysis is consistent with the observation that WNT signalling pathway is abnormally regulated in NPC.

32Another pathway preferentially targeted by the down-modulated miRNAs involved growth factor-regulated G1-S cell cycle progression. Several cell cycle regulators, including cyclin D2, cyclin E2, and CDC25A are coregulated by the downregulated miRNAs. Reverse transcriptase PCR analysis confirmed the elevated expression of these target genes in NPC tissues. Correlation analysis further revealed a significant inverse correlation between the level of cyclin E2 and three of its cognate miRNAs, miR-34c, miR-200a, and miR-9. These results suggest that the aberrant cell cycle control in NPC tissues may be tightly linked to the dysregulated expression of miRNAs.

33The third signalling pathway, significantly targeted by down-modulated miRNAs, is the VEGF-regulated angiogenesis. Predicted targets include VEGF-A, PLC-γ, and AKT. Earlier studies have shown that VEGF is widely expressed in NPC tissues, and elevated VEGF expression is associated with local recurrence and distal metastasis. Current studies revealed that several key regulators on the VEGF-regulated signalling pathway are cotargeted by multiple miRNAs, including miR-135a, miR29c, miR-195, and miR-497. Simultaneous down-modulation of these miRNAs may be responsible for the elevated VEGF signalling in NPC.

34Recent studies using an miRNA target reverse screening method showed that miR-16 family triggers a cell cycle arrest by silencing multiple cell cycle regulatory genes simultaneously, rather than the individual target. Similarly, miR-17-5p regulates the cell cycle progression by coordinately suppressing more than 20 genes involved in the G1/S transition. It is evident that a better connection between cellular function and miRNA expression requires both computational as well as experimental approach. The context score filter and the use of control gene list significantly increased our confidence for target prediction and pathway analysis. Experimental validation and literature mining confirmed that the expression data detected in clinical samples are consistent with the computational prediction. This study has shown the systems biology approach as a valuable tool to uncover global functional change associated with miRNA dysregulation. Similar approach can be used to assist the functional interpretation of miRNA expression data in other human diseases.