MicroRNA-21 regulates prostaglandin E2 signaling pathway by targeting 15-hydroxyprostaglandin dehydrogenase in tongue squamous cell carcinoma

Background Oral tongue squamous cell carcinoma (OTSCC) is one of the most aggressive forms of head and neck/oral cancer (HNOC), and is a complex disease with extensive genetic and epigenetic defects, including microRNA deregulation. Identifying the deregulation of microRNA-mRNA regulatory modules (MRMs) is crucial for understanding the role of microRNA in OTSCC. Methods A comprehensive bioinformatics analysis was performed to identify MRMs in HNOC by examining the correlation among differentially expressed microRNA and mRNA profiling datasets and integrating with 12 different sequence-based microRNA target prediction algorithms. Confirmation experiments were performed to further assess the correlation among MRMs using OTSCC patient samples and HNOC cell lines. Functional analyses were performed to validate one of the identified MRMs: miR-21-15-Hydroxyprostaglandin Dehydrogenase (HPGD) regulatory module. Results Our bioinformatics analysis revealed 53 MRMs that are deregulated in HNOC. Four high confidence MRMs were further defined by confirmation experiments using OTSCC patient samples and HNOC cell lines, including miR-21-HPGD regulatory module. HPGD is a known anti-tumorigenic effecter, and it regulates the tumorigenic actions of Prostaglandin E2 (PGE2) by converts PGE2 to its biologically inactive metabolite. Ectopic transfection of miR-21 reduced the expression of HPGD in OTSCC cell lines, and the direct targeting of the miR-21 to the HPGD mRNA was confirmed using a luciferase reporter gene assay. The PGE2-mediated upregulation of miR-21 was also confirmed which suggested the existence of a positive feed-forward loop that involves miR-21, HPGD and PGE2 in OTSCC cells that contribute to tumorigenesis. Conclusions We identified a number of high-confidence MRMs in OTSCC, including miR-21-HPGD regulatory module, which may play an important role in the miR-21-HPGD-PGE2 feed-forward loop that contributes to tumorigenesis. Electronic supplementary material The online version of this article (doi:10.1186/s12885-016-2716-0) contains supplementary material, which is available to authorized users.


Background
Head and neck/oral cancer (HNOC) is a commonly encountered malignancy. Head and neck squamous cell carcinoma (HNSCC), which arises from the epithelium lining of this region, makes up the majority (over 90 %) of HNOC. Oral tongue squamous cell carcinoma (OTSCC) is one of the most aggressive form of HNSCCs, which exhibits a propensity for rapid local invasion and spread [1], has a distinct nodal metastasis pattern [2,3]. OTSCC patients also suffer from a high recurrence rate [4]. OTSCC is a complex disease with extensive genetic and epigenetic defects, including microRNA deregulation. MicroRNAs are pivotal regulators of physiological and disease processes through their control of diverse cellular processes. Several microRNAs have been functionally classified as oncogenes or tumor suppressors, and the aberrant expression of microRNA has been observed in almost all cancer types including OTSCC [5][6][7][8]. Deregulation of these cancer-associated microRNAs can significantly impact tumor initiation and progression by activating pathways promoting uncontrolled proliferation, favoring survival, inhibiting differentiation, and promoting invasion [9,10]. MicroRNAs are not directly involved in protein coding, but are able to control the expression of their target genes at post-transcriptional levels by facilitating mRNA degradation and/or repressing translation. As such, the identification and detection of functional microRNA-mRNA regulatory modules (MRMs) are crucial components for understanding of micro-RNA functions.
MicroRNAs are a class of small non-coding RNAs of approximately 22 nucleotides in length that are endogenously expressed in mammalian cells. They are related to, but distinct from, siRNAs. A key difference between siRNA and microRNA is that siRNA requires almost complete complementary to its targeting sequence for it to exert the silencing function, whereas microRNA usually binds to its target genes through partial complementary. While numerous sequence-based bioinformatics methods for microRNA target prediction have been developed, these methods often lead to high false discovery rates [11]. In order to minimize false positives and to detect the functional microRNA targets under a specific biological condition, recent approaches often integrate the microRNA and mRNA profiling analysis in conjunction with the sequence-based target prediction. Two types of experiments are common: 1) differential mRNA profiling experiment on a microRNA transfected cell line and its negative control, and 2) simultaneous microRNA and mRNA profiling analysis on samples of different phenotypes (e.g., normal vs. tumor). The first approach has been used by many groups, including us, to define the functional micro-RNA targets when a specific microRNA is over-or under-expressed [12][13][14]. The second approach aims to discover microRNA with altered expression related to different phenotypes and to uncover their targets mRNAs. This approach is based on the simple principle that inverse relationships in their expression profiles should be held between a specific microRNA and its functional target genes. When integrated with the sequence-based bioinformatics target prediction, this approach is believed to lead to the identification of high confidence microRNA targets.
Our group and several others have recently undertaken extensive RNA-based surveys to identify gene expression and microRNA abnormalities in OTSCC. In this study, we utilized our existing transcription profiling dataset [15], and a meta-analysis of 13 published microRNA profiling studies [16], and integrate them with a collection of 12 sequence-based bioinformatics tools to define the deregulation of functional MRMs in OTSCC. We then evaluated these MRMs in 2 OTSCC patient cohorts and a panel of HNSCC cell lines. With our comprehensive approach, we identified a panel of high confidence microRNA-mRNA regulatory modules in OTSCC, including miR-21-15-Hydroxyprostaglandin Dehydrogenase (HPGD) regulatory module. We also confirmed the positive feed-forward loop that involves miR-21, HPGD and Prostaglandin E2 (PGE2) in HNOC cells that contribute to tumorigenesis.

MicroRNA target prediction
The microRNA target prediction was performed using the comparative analysis function of the miRWalk [17], which contains a collection of 10 bioinformatics tools, including DIANAmT, miRanda, miRDB, miRWalk, RNAhybrid, PicTar (4-way), PicTar (5-way), PITA, RNA22, TargetS-can5.1. In addition, MicroCosm 5.0 and TargetScanHuman 6.2 were also used for predicting the microRNA targets. For our study, genes that were predicted by at least one method were defined as candidate microRNA targets. The base-pairing and the minimum free energy (mfe) for the binding of microRNA to its targeting sequences were predicted using the RNAhybrid program [18].

Clinical samples from OTSCC patients
We downloaded the RNASeq and miRNASeq profiling datasets on 12 OTSCC and paired normal tissue samples from The Cancer Genome Atlas (TCGA) Data Protal [tcga-data.nci.nih.gov]. The gene expression values were extracted as normalized count, and the microRNA levels were extracted as reads per million miRNA mapped from the datasets. The demographics of the patients were as follows: 6 male, 6 female and average age = 62 (range: 36-88), 1 stage T1 cases, 5 stage T2 cases, 3 stage T3 case and 3 T4 cases. Oral cytology samples were obtained from 13 patients with pathologically characterized primary OSCC of the tongue before tumor resection (including 6 stage T1 cases 6 stage T2 cases and 1 stage T3 case) as previously described [29,30]. These procedures are in compliance with the Helsinki Declaration, and was approved by the Ethical Committee of the First Affiliated Hospital, Sun Yat-Sen University (reference number: 2014-C-001). The informed consent was obtained from participants. Patients were excluded if there is a history of lung carcinoma or HNSCC elsewhere and may represent metastatic disease. The demographics of the patients were as follows: 8 male, 5 female and average age = 51.8 (range: 32-78). The total RNA was isolated using miRNeasy Mini kit (Qiagen), and quantified by a spectrophotometer or the RiboGreen RNA Quantitation Reagent (Molecular Probes).

Quantitative RT-PCR Analysis
The relative microRNA levels were determined by Taq-Man microRNA assays (Applied Biosystems) as previously described [16,31]. The relative mRNA levels were determined by quantitative two-step RT-PCR assay with pre-designed gene specific primer sets (Origene) as described before [16,31]. The relative microRNA and mRNA levels were computed using the 2 -delta delta Ct analysis method, where U6 and beta-actin were used as internal controls, respectively.

Western-blot analysis
Western blots were performed as described previously [16] using antibodies specific for HPGD (Cayman Chemical) and beta-actin (Sigma-Aldrich) and an immuno-star HRP substrate Kit (Bio-RAD).

Fluorescent immunocytochemical analysis
Immunofluorescence analysis was performed as previously described [16]. In brief, cells were cultured on 8 chamber polypropylene vessel tissue culture treated glass slides (Millipore) fixed with cold methanol, permeabilized with 0.5 % Triton X-100/PBS, and blocked with 1% BSA in PBS. The slides were incubated with primary antibodies against HPGD (1:500, Cayman Chemical). The slides were then incubated with a FITC-conjugated anti-rabbit IgG antibody (1:50, Santa Cruz). The slides were mounted with ProLong Gold antifade reagent containing DAPI (Invitrogen) following the manufacturer's protocol. The slides were then examined with a fluorescence microscope (Carl Zeiss).

Statistical analysis
Data was analyzed using the Statistical Package for Social Science (SPSS), version 17.0. Student's t-test was used to compare differences between groups. Pearson's correlation coefficient was computed for examining the relationship between the expression of microRNA and their target genes. For all analyses, p < 0.05 was considered statistically significant.

Results
We first developed a list of putative microRNA-mRNA regulatory modules (MRMs) based on the simple principle that inverse relationships should be anticipated in the expression of a specific microRNA and its functional target gene (mRNA). We used a total of 97 differentially expressed coding genes (44 up-regulated and 53 down-regulated mRNAs, see Additional file 1: Table S1A and S1B, respectively) and 9 differentially expressed  The putative microRNA-mRNA regulatory module (MRM) was constructed based on microRNA and mRNA expression profiles of OTSCC, as we previously reported in [16] and [15], respectively b Differential expression of microRNAs and mRNAs was validated using dataset on 12 OTSCC and paired normal tissue samples that were extracted from TCGA. Genes that show statistically significant differential expression were identified with bold font c The candidate targets of a microRNA were predicted using a collection of 12 bioinformatics tools, including DIANAmT, miRanda, microCosm, miRDB, miRWalk, RNAhybrid, PicTar (4-way), PicTar (5-way), PITA, RNA22, TargetScan5, and TargetScanHuman 6.2. The number of bioinformatics tools (out of a total of 12 tools tested here) that predict a gene to be a microRNA target was presented. The gene/microRNA pairs predicted by at least 3 tools were listed in the table   d Correlations of microRNA and mRNA levels were assessed using dataset on 12 OTSCC and paired normal controls that were extracted from TCGA. Inverted correlation (negative Pearson r value) is expected for a MRM, and p value was calculated e Correlations of microRNA and mRNA levels were assessed by quantitative real-time PCR based on 13 HNSCC cell line. Inverted correlation (negative Pearson r value) is expected for a MRM, and p value was calculated f Correlations of 4 pairs of microRNA and mRNA levels were assessed by quantitative real-time PCR based on 13 OTSCC patient oral cytology samples. Inverted correlation (negative Pearson r value) is expected for a MRM, and p value was calculated microRNAs (5 up-regulated and 4 down-regulated microRNAs, see Additional file 1: Table S1C) from our previous genomic profiling studies on OTSCC [15,16] for the development of this putative MRMs list. This putative MRMs list consists of 265 putative MRMs defined by microRNA up-regulation and mRNA downregulation, and 176 putative MRMs defined by micro-RNA down-regulation and mRNA up-regulation. We then tested these putative MRMs using a panel of 12 different sequence-based microRNA target prediction algorithms (DIANAmT, miRanda, microCosm, miRDB, miRWalk, RNAhybrid, PicTar (4-way), PicTar (5-way), PITA, RNA22, TargetScan5.1, and TargetScanHuman6.2) to refine our putative MRMs list. A total of 132 candidate MRMs were identified (predicted as microRNA target by at least 1 bioinformatics algorithm, see Additional file 2: Table S2A and Additional file 3: Table S2B). As shown in Table 1, 38 potential MRMs were predicted by at least 3 bioinformatics target prediction algorithms, where the up-regulation of the microRNA contributes to the down-regulation of mRNA, and 15 potential MRMs were predicted by at least 3 bioinformatics target prediction algorithms (Table 2), where down-regulation of the microRNA contributes to the up-regulation of mRNA.
The differential expression of microRNAs and coding genes (mRNAs) involved in these 53 potential MRMs (9 microRNAs and 34 mRNAs) was then validated using dataset on 12 OTSCC and paired normal tissues (extracted from TCGA Data Portal). As shown in Additional file 4: Table S3, statistically significant differential expression was observed for 8 out of 9 microRNAs and 23 out of 34 mRNAs tested in the validation OTSCC cohort.
To further evaluate these potential MRMs, we examined the correlative relationship between the microRNA levels and the expression of their target genes in these 12 OTSCC and 12 paired normal tissues (extracted from TCGA Data Portal), as well as 13 HNSCC cell lines ( Table 1   The putative microRNA-mRNA regulatory module (MRM) was constructed based on microRNA and mRNA expression profiles of OTSCC, as we previously reported in [16] and [15], respectively b Differential expression of microRNAs and mRNAs was validated using dataset on 12 OTSCC and paired normal tissue samples that was extracted from TCGA. Genes that show statistically significant differential expression were identified with bold font c The candidate targets of a microRNA were predicted using a collection of 12 bioinformatics tools, including DIANAmT, miRanda, microCosm, miRDB, miRWalk, RNAhybrid, PicTar (4-way), PicTar (5-way), PITA, RNA22, TargetScan5, and TargetScanHuman 6.2. The number of bioinformatics tools (out of a total of 12 tools tested here) that predict a gene to be a microRNA target was presented. The gene/microRNA pairs predicted by at least 3 tools were listed in the table   d Correlations of microRNA and mRNA levels were assessed using dataset on 12 paired OTSCC and normal controls that was extracted from TCGA Data Portal. Inverted correlation (negative Pearson r value) is expected for a MRM, and p value was calculated e Correlations of microRNA and mRNA levels were assessed by quantitative real-time PCR based on 13 HNSCC cell line. Inverted correlation (negative Pearson r value) is expected for a MRM, and p value was calculated miR-21-HLF and miR-21-HPGD were also statistically significant in HNSCC cell lines.
We further explore the interaction of miR-21 and HPGD in our study. As shown in Fig. 2a, ectopic transfection of miR-21 mimic to UM1, UM2, SCC9 and Tca8113 cells led to a statistically significant reduction in HPGD mRNA level as compared to cells treated with control mimic. The miR-21 has no apparent effect on HPGD expression in HeLa cells. As shown in Fig. 2b and c, ectopic transfection of miR-21 mimic to UM1 cells led to reduced HPGD expression at protein level and reduced immunostaining of HPGD, respectively, as compared to the cells treated with control mimic. As shown in Fig. 2d, ectopic transfection of miR-21 mimic also enhanced the proliferation of OTSCC cells, which is consistent with previous observations [10,32], and confirmed the oncogenic effect of miR-21.
Bioinformatics analysis revealed that there are three miR-21 targeting sites located in the 3′-UTR of the HPGD mRNA (E1 at position 2652 to 2671, E2 at position 2880 to 2901, E3 at 2890 to 2911) and the targeting sites E2 and E3 are partially overlapped (Fig. 3a). The predicted minimum free energy (mfe) for the binding of these sites to miR-21 are −17.6, −11.4 and −16.5 kcal/ mol, respectively. To test whether the miR-21 directly interacts with these predicted targeting sites in HPGD mRNA, dual luciferase reporter assays were performed using constructs containing these targeting sites (Fig. 3b). When cells were transfected with miR-21, the luciferase activities of the construct containing targeting site E1 (pGL-E1) was significantly reduced as compared to the cells transfected with negative control. When the seed region of this targeting site was mutated (pGL-E1m), the effect of miR-21 on the luciferase activity was abolished. For sites E2 and E3, when cells were transfected with miR-21, the luciferase activities of the construct containing both targeting sites E2 and E3 (pGL-E2E3) was not changes as compared to the cells transfected with negative control. Interestingly, when the seed region of E2 was mutated (pGL-E2mE3), the miR-21-mediated down-regulation of the luciferase activity was observed. MiR-21 has no effect  GPD1L (a, b, c), HLF (d, e, f), HPGD (g, h, i), and the correlation of miR-130b level with the expression of MGLL (j, k, l) were assessed, and the Pearson's correlation coefficient (r) was calculated on constructs with E3 mutation (pGL-E2E3m) or mutations of both E2 and E3 (pGL-E2mE3m).
As shown in Fig. 4a and b, both siRNA-mediated knockdown of COX2 and treatment with COX2 inhibitor (CelecoxiB) led to down-regulation of miR-21 in UM1 cells. As shown in Fig. 4c, directly apply PGE2 to the UM1 cells led to the up-regulation of miR-21, and knockdown of HPGD (Fig. 4d) also led to the upregulation of miR-21. As anticipated, treating cells with PGE2 and CelecoxiB led to up-regulation and downregulation of cell proliferation, respectively, which is consistent with previous observations [33,34] (Fig. 4e). These results are in agreement with observation made by Lu et al. in cholangiocarcinoma [35], which confirm the PGE2-mediated miR-21 up-regulation in OTSCC and suggest a PGE2-miR-21-HPGD positive feed-forward loop that contributes to tumorigenesis (Fig. 4f).

Discussion
Despite the significant increase in the number of experimentally validated microRNA-mRNA regulatory relationships, the majority of the microRNA targeted genes remains unknown. MicroRNA usually binds to its target genes through partial complementary. While numerous sequence-based bioinformatics methods for microRNA target prediction have been developed, these methods often lead to high false discovery rates [11]. However, the integration of these bioinformatics tools with mRNA/microRNA differential expression profiles often lead to the identification of high confidence microRNA-mRNA regulatory modules. In this study, we carried out this integrated analysis to identify MRMs in two steps. First, based on the simple principle that inverse relationships should be anticipated in the expression of a specific microRNA and its functional target genes, we developed a list of putative microRNA-mRNA regulatory modules by linking each microRNAs with all inversely regulated mRNAs based on the results of our previous mRNA and microRNA profiling studies on OTSCC [15,16]. The second step is to these putative MRMs bioinformaticsly using sequence-based micro-RNA target prediction algorithm. Since there are many available sequence-based microRNA target prediction tools, and each of these tools utilizes a different model to define targeting sequences that are associated with functionality, the predictions differ when applied to the same microRNAs, with each method having different levels of coverage and false positive prediction [11]. In order to reduce the potential false positives, we used a voting scheme to combine the predictions from the 12 commonly used bioinformatics tools, including DIANAmT, miRanda, microCosm, miRDB, miRWalk, RNAhybrid, PicTar (4-way), PicTar (5-way), PITA, RNA22, TargetS-can5.1, and TargetScanHuman6.2. With this integrated approach, we developed a list of 53 potential MRMs that are differentially expressed in OTSCC.
Since the microRNA regulates its target gene mainly at post-transcriptional level, inverse correlation between Deregulations of miR-21 and miR-130b, as well as deregulation of GPD1L, HLF, HPGD and MGLL have been reported either in HNOC or other cancer types [15,16,[36][37][38][39][40][41][42], and these MRMs represent significant functional relevance in OTSCC. GPD1L has the glycerol-3-phosphate dehydrogenase enzyme activity and is a regulator of HIF-1α stability [40]. And a recent study showed that the GPD1L expression is a strong predictor for local recurrence and survival in HNSCC [39]. HLF belongs to the PAR (proline and acidic amino acid-rich) subfamily of bZIP transcription factors [43,44], and plays a role in development and circadian rhythm regulation in the mammalian. HLF fusion proteins that resulted from chromosomal translocation (e.g., E2A-HLF) are often linked to leukemia. However, the role of HLF in OTSCC is not entirely clear. MGLL is involved in Prostaglandin E2 (PGE2) production in response to inflammation and infection which leads to fever [45]. Arachidonic acid (AA), a precursor for PGE2, is typically liberated from AA-containing phospholipids by the action of phospholipases A2 (PLA2s). MGLL is a monoacylglycerol lipase which hydrolyzes 2arachidonoylglycerol (2-AG), an endocannabinoid that functions in the central nervous system, to AA and glycerol, representing an alternative AA-producing pathway. MGLL may also play a role in certain types of cancer by regulating both endocannabinoid and fatty The base-pairing and the minimum free energy (mfe) for the binding of miR-21 to the targeting sequences were predicted using the RNAhybrid program [18]. b Dual luciferase reporter assays were performed to test the interaction of miR-21 and its targeting sequences in the HPGD mRNA using constructs containing the predicted targeting sequences (pGL-E1 and pGL-E2E3) and mutated targeting sequences (pGL-E1m, pGL-E2mE3, pGL-E2E3m, pGL-E2mE3m) cloned into the 3′-UTR of the reporter gene. Data represent at least 3 independent experiments with similar results. *: p < 0.05 acid pathways, and supporting protumorigenic metabolism [46]. This appears to be contradict with the apparent down-regulation of MGLL observed in OTSCC [15], and the miR-130b-MGLL regulatory module predicted here. Nonetheless, whether MGLL plays a role in OTSCC and, if so, by what mechanism are questions that remain unanswered. HPGD is a known anti-tumorigenic effecter, and it regulates the tumorigenic actions of Prostaglandin E2 (PGE2) by converts PGE2 to its biologically inactive metabolite, and down-regulation of HPGD has been observed in many human cancer types [47][48][49][50][51][52][53]. Since miR-21 is one of the most consistently observed up-regulated microRNA in OTSCC [16,54], the miR-21-HPGD regulatory module may represents a critical mechanism of regulating PGE2 signaling.
Our functional study confirmed the effect of miR-21 on HPGD expression level, and the direct interaction of miR-21 with the HPGD mRNA in OTSCC cells. We identified three miR-21 targeting sites located in the 3′-UTR of the HPGD mRNA, including a previously reported site (E1) [35], and two partially overlapped sites (E2 and E3). While we confirmed the miR-21-mediated and E1 site-dependent target gene downregulation, E2 and E3 sites appear to have no effect. This may be because that targeting sites E2 and E3 are partially overlapped, and may interfere with the proper interaction with the RISC complex. The elimination of E2 may partially restore the capability of E3 (which has a stronger binding affinity among the two sites) to binding to the RISC complex. This is different than our previous observation where miR-138 was able to interact with multiple overlapping target sites on the FOSL1 mRNA [55]. Additional studies are needed to explore this mutual exclusive phenomenon among multiple targeting sites. The HPGD gene has 6 known transcript variants (NCBI accession: NM_000860, NM_001145816, NM_001256301, NM_001256305, NM_001256306, NM_001256307), and all 6 variants have the same 3′-UTR. As such, the interaction between miR-21 and HPGD mRNA is not likely to be affected by alternative splicing. Interestingly, we did not observe any miR-21 effect on HPGD expression in HeLa cells (a cell line that originated from a cervical cancer case). It is possible that this apparent difference in the miR-21 effect on HPGD expression may be due to differences in cancer types. It is worth noting that the effect of miR-21 on HPGD expression has also been observed in other cancer type [35]. Alternatively, this difference may be cell-line specific. HeLa cells (or the OTSCC cell lines used here) may have specific mutation(s) that dictate the miR-21 effects on HPGD. More in-depth functional analysis will be needed to fully evaluate the miR-21-HPGD regulatory module in different cancer types and in other biological systems.
The levels of COX2 and its catalytic product PGE2 are increased in a variety of malignancies, including HNOC [56][57][58][59]. The tumorigenic actions of PGE2 are attributable to its modulation of cell proliferation, survival, migration, and invasion. The level of PGE2 is controlled by the status of PGE2 synthesis and degradation. Whereas the cyclooxygenases (COX1 and COX2) are rate-limiting key enzymes that control PGE2 biosynthesis, HPGD is a key enzyme that converts PGE2 to its biologically inactive metabolite, 13,14-dihydro-15-keto-PGE2, thus leading to PGE2 inactivation [60,61]. Consistent with the antitumorigenic effect of HPGD, the down-regulation of HPGD has been observed in many human cancer types [47][48][49][50][51][52][53]. Lu et al., first reported the PGE2-mediated up-regulation of miR-21 in cholangiocarcinoma, and suggested a positive feed-forward loop that involves PGE2, miR-21 and HPGD [35]. Our results are consistent with these previous observations, and confirm the existence of a PGE2-miR-21-HPGD positive feed-forward loop in OTSCC that contributes to tumorigenesis (Fig. 4f ).