A Quest for miRNA Bio-Marker: A Track Back Approach from Gingivo Buccal Cancer to Two Different Types of Precancers

Deregulation of miRNA expression may contribute to tumorigenesis and other patho-physiology associated with cancer. Using TLDA, expression of 762 miRNAs was checked in 18 pairs of gingivo buccal cancer-adjacent control tissues. Expression of significantly deregulated miRNAs was further validated in cancer and examined in two types of precancer (leukoplakia and lichen planus) tissues by primer-specific TaqMan assays. Biological implications of these miRNAs were assessed bioinformatically. Expression of hsa-miR-1293, hsa-miR-31, hsa-miR-31* and hsa-miR-7 were significantly up-regulated and those of hsa-miR-206, hsa-miR-204 and hsa-miR-133a were significantly down-regulated in all cancer samples. Expression of only hsa-miR-31 was significantly up-regulated in leukoplakia but none in lichen planus samples. Analysis of expression heterogeneity divided 18 cancer samples into clusters of 13 and 5 samples and revealed that expression of 30 miRNAs (including the above-mentioned 7 miRNAs), was significantly deregulated in the cluster of 13 samples. From database mining and pathway analysis it was observed that these miRNAs can significantly target many of the genes present in different cancer related pathways such as “proteoglycans in cancer”, PI3K-AKT etc. which play important roles in expression of different molecular features of cancer. Expression of hsa-miR-31 was significantly up-regulated in both cancer and leukoplakia tissues and, thus, may be one of the molecular markers of leukoplakia which may progress to gingivo-buccal cancer.


Introduction
Gingivo buccal squamous cell carcinoma (GBSCC) is one of the most prevalent (,60%) cancers in oral cavity, especially among the tobacco users in India. Entire set of head and neck squamous cell carcinoma stands as the fifth most common malignancy worldwide [1]. But head and neck cancer comprises ,24% of total cancers in India as recorded at a tertiary hospital, Tata Memorial Centre, Mumbai and about ,13.5% of them are from the oral cavity [2]. Five years survival rate (,50%) of patients suffering from head and neck squamous cell carcinoma has not improved much even after intense research during last 15 years, so, early detection is still a key issue for better survival [3].
Since 1993, miRNA has emerged to be one of the most prominent biological regulators, which play important roles in controlling and fine tuning its target's (mRNA) expression [4][5][6][7][8][9]. In recent years, it has been demonstrated that microRNAs (miRNAs) are also involved in human tumorigenesis and could act either as tumorigenic/oncogenic or anti-tumorigenic molecules. Thus, it is making a new layer in the molecular events in human cancer. Gene expression studies revealed that many miRNAs are deregulated in different cancer types and functional studies clarified that miRNAs are involved in several molecular and biological processes that drive tumorigenesis [10][11][12]. Thus, in addition to different scales of variability in terms of new mutations in genome or genetic background of the patients (i.e. germ line mutation), variability in expression of different miRNAs is also a major aspect for investigation. With this belief, we studied expression deregulation of 762 miRNAs in GBSCC and also checked whether similar kind of expression deregulation could be observed in precancerous leukoplakia and lichen planus tissues from oral cavity. Efforts have been made to understand how these miRNAs may be involved in cancer using pathway analysis and information from literatures and databases.

Samples
This study was approved by ''Review committee for protection of research risk to humans, Indian Statistical Institute''. All individuals in this study have provided written informed consents to publish case details. All participants signed a questionnaire containing demography, tobacco habit and a statement describing that he/she has no objection for use of blood and tissue samples in this study and is participating in this study voluntarily. Unrelated patients suffering from cancer (n = 18), lichen planus (n = 12) and leukoplakia (n = 18) were considered for this study from Guru Nanak Institute of Dental Science and Research, Panihati, Kolkata, (Table 1, Table 2, Table 3). One tissue punch from cancer/precancer site and another punch from adjacent ''clinically'' normal site (at least 1.5 inch away from the border of the lesion) were biopsied by oral pathologist. One portion of biopsy tissues was used for histopathological confirmation and remaining part of tissues were kept separately in ''RNA Later'' at 280uC and used within 2 months.
RNA isolation and TLDA assay RNA isolation was done using mirVana kit (Life Technologies, USA). RNA yield quantification was done by Nano Drop and integrity was checked with Agilient's Bioanalyzer (Agilent RNA6000 Nano Kit), respectively. RIN number cut-off was chosen as $6.9. Parallel expression assay of 762 miRNAs was performed using TLDA-A (V2) and TLDA-B (V3) card in 7900HT FAST Real Time PCR system (Applied Biosystems, USA) using TLDA flat block [13]. Primary analysis was done using SDS and Data Assist (Life Technologies, USA) software packages to get expression in terms of Ct, DCt and DDCt values where Ct = Cycles at which the PCR product quantity reaches a defined threshold, DCt = Ct of a miRNA in cancer tissue -Ct of geometric mean of expression of 3 most stable endogenous control miRNAs in that tissue and DDCt = DCt of a miRNA in cancer tissue -DCt of that miRNA in control tissue . Out of the assayed 762 miRNAs, five were candidates for endogenous controls. However on the basis of their expression stability across all samples, RNU-48, RNU-44 and U6/mmu-6 were selected as endogenous controls. To get DCt as the normalized measure of expression of a miRNA in both cancer/ precancer and control tissues, expression of all miRNAs in tissues was normalized independently using geometric mean of expression of selected 3 endogenous control miRNAs in the same tissue. Ct value less than 40 were only considered for further analysis.

Analysis
Data Pruning. If expression of any particular miRNA was observed to be present in at-least 9 out of 18 cancer-normal paired tissues, then only, it has been considered for further downstream analysis. Annotation of TLDA assay V2 (A-card) and V3 (B-card) follows miRBase release-14 [13][14]. Now, expression data of those miRNAs were considered for further analysis whose annotation is still valid in miRBase release-19 [15][16]. In this way expression of 531 miRNAs was considered for subsequent down-stream analysis.
Statistical Analysis.  [17] and corrected cut off p-value was 0.00065.
Cluster Analysis of miRNAs expressed in cancer tissues. K-median clustering method was used for the entire pruned data set to see whether genome wide miRNA expression variations between individuals were large enough to reliably divide the samples into a number of sub clusters. For this clustering, chosen distance metric was Euclidean. Since, expression of many miRNAs was absent in our data set, it was creating ''Absent data'' situation. This situation is known to bias the clustering if the most popular clustering method-K-means was adopted. So, the choice of clustering was K-median. The number of reliable clusters, the data would form, was determined using ''the elbow'' method.
TaqMan Assay. Expression of significantly de-regulated miRNAs, detected in cancer tissues by TLDA method, were also validated in same cancer tissues and examined in leukoplakia and lichen planus tissues by TaqMan assay (7900HT Fast Real Time PCR system, Applied Biosystems, USA). Probes and primers were supplied by Invitrogen India Ltd and data were retrieved as ''fold change'' compared to adjacent controls. Normalization of expression of each gene in each sample was done using geometric

Expression profile in cancer tissues by TLDA
Expression data of 762 miRNAs were assayed by this method but after data pruning, expression data of 531 miRNAs (including 5 endogenous controls) was used for other down-stream analysis (File S1). During statistical test, at least 2-fold expression change in cancer tissue compared to its paired-control tissue was considered to be the bench-mark of expression deregulation. Thus, expression of 7 miRNAs was found to be significantly deregulated after multiple test correction and all of these seven miRNAs had .4fold average expression deregulation (i.e. either up-or downregulation). Validation of expression by miRNA specific TaqMan assay also reconfirmed expression deregulation in cancer tissues although fold-changes were diminished compared to TLDA outcome ( Table 4). The difference in sensitivity of these two RT-PCR based (TLDA and miR-specific TaqMan Assay) experimental methods, coupled with the fact that these two sets of experiments were performed at two different time points, might be the influencing factors for difference in degrees of expression deregulation of miRNAs. Relative locations of these 7 miRNAs along with other miRNA genes across different chromosomes revealed that hsa-miR-133a and hsa-mir-7 are located on two (Chr 18 and Chr 20) and three chromosomes (Chr 9, Chr 15 and Chr 19), respectively ( Figure 1a). Here, assays were performed for only mature miRNAs, thus, expression of hsa-miR-133a and hsamir-7 might be cumulative sum of all mature forms of respective miRNAs. Among these 7 miRNAs, expression of 4 miRNAs viz. hsa-miR-1293, hsa-miR-31, hsa-miR-31* and hsa-miR-7 were significantly up-regulated and those of 3 miRNAs viz. hsa-miR-206, hsa-miR-204 and hsa-miR-133a were significantly down regulated in cancer samples (Table 4) Normalized expression (DCt) of these 7 miRNAs in cancer and control tissues were mostly non-overlapping and DCt values of different miRNAs across the control tissues also showed quite a wide range of variation ( Figure 2). So to get correct relative expression, it is important to compare expression of miRNAs in cancer tissue with those of control tissue from the same individual. In this figure, miRNAs were placed according to increasing pvalues (from top to bottom) vertically. The DCt values of 8 miRNAs in different tumor and control samples had been plotted on the horizontal axis to show distribution of DCt values across all the samples. It is evident that more the overlap between ranges of expression of miRNAs in cancer and control tissues, less is the level of significance. Here, hsa-mir1293 with lowest p value 0.000028 had been positioned on the top and hsa-mir-1 with p value of 0.00094 (just below the corrected level of significance) was placed at the bottom on the vertical axis ( Figure 2).

Expression of miRNAs in precancerous leukoplakia and lichen planus samples
Among 7 miRNAs which was observed to be significantly deregulated in cancer samples (Table 4), expression of only miR-31 was significantly up-regulated in leukoplakia tissues where as expression of none of these 7 miRNAs was deregulated significantly in lichen planus tissues. Expression of miR-31* and miR-204 was also up-(4.75 folds) and down-regulated (1.99 folds), respectively, in leukoplakia tissues but not significantly different.

Database mining and bioinformatics analyses
Published reports on these 7 miRNAs had also shown similar direction of expression deregulation in some cancers including head and neck [10,[18][19][20]. A total of 561 unique targets were identified when these 7 miRNAs were used to search targets using miRWalk [21] and further cross-validated from Pubmed (http:// www.ncbi.nlm.nih.gov/pubmed). The hsa-miR-31* has validated target, RhoA, which is reportedly implicated in mouth neoplasm [18]. The hsa-miR-1293 till now is known to target GCN1L1 (Table 5) and hsa-miR-1293 mediated down regulation of this tumor antigen gene (i.e. GCN1L1) could contribute to poor prognosis of the tumor [22]. IPA tool was used for ''disease term search'' (IngenuityH Systems, www.ingenuity.com) and most significant ''disease term'' in cancer category of IPA was ''head and neck'' cancer. IPA could not provide any hit by hsa-miR-1293 and considered hsa-miR-31 and hsa-miR-31* as synonymous so the output was from 5 miRNAs. It was also noticed that IPA considered hsa-miR-206 synonymous to hsa-miR-1 since they have identical seed sequence [23]. It was also observed that expression deregulation status of hsa-miR-1 and hsa-miR-206 across the samples were very similar to each other (r = 0.94) (Figure 1c). Inclusion of hsa-miR-1 in target search list increased the number of targets to 702. In KEGG pathway mapping portal [24][25] validated targets of these 8 miRNAs have been used. Top most pathways in the mapped list were ''microRNAs in cancer'' followed by ''proteoglycans in cancer'' and ''global cancer pathway'' (Table S1 in File S2). Other top relevant pathway was PI3K-AKT which is one of the most common pathways implicated in cancer. Other most significant changes that could have happened due to these miRNAs are disruption of actin cytoskeleton maintenance and focal adhesion. Other probable implicated signaling pathways were RAS, RAP1, MAPK, HIF-1, FOXO, TNF, ErBB, apoptosis etc (Table S1 in File S2). ''GO'' term enrichment search was performed using input of 702 validated targets of these eight miRNAs (data not shown) and it was observed that most prominent disrupted biological processes are primarily related to cell migration, loss of apoptosis, cell proliferation etc. But DAVID based ''GO'' term enrichment for biological process showed most prominent biological process is ''Apoptosis'' (Table S2 in File S2) [26][27][28]. All these observations support that these 8 miRNAs may play important functional roles in gingivo buccal cancer. Our RNA-Seq data shows that expression of FN1, MSN and MMP9, which belongs to ''Proteoglycans in cancer'' gene list and are targets of downregulated miRNAs, was up-regulated in 10 of 13 tissue samples (on average, 6.83, 2.43 and 21.89 folds, respectively, compared to adjacent normal). Similarly predicted up-regulation of LAMC and down regulation of PAI, which belong to PI3K-AKT pathway, was also validated in RNA-Seq data in a sub set of these tissue samples (unpublished data). These observations also support our predictive pathway analysis regarding involved biological process/ pathways that could be targeted by this set of 8 miRNAs.
Selected samples had little variation in terms of its differentiation status or clinical staging. Cluster analysis was performed to check whether genome wide miRNA profile could explain such characteristic variation. It was performed using expression deregulation data (DDCt) of 531 miRNAs from 18 cancer samples using K-median method. Optimally, only two distinct clusters were obtained; one consisted of 13 paired tumor-normal samples and other consisted of 5 paired tumor-normal samples. It was evident that these two clusters of samples were not formed on the basis of their differentiation status or clinical stages. Expression of 30 miRNAs was significantly deregulated in the cluster of 13 samples (p-value cut off 0.00298 after Benjamini-Hochberg Correction, Figure 3a, Figure 3b) but none of the miRNAs were significantly deregulated in the cluster of 5 samples (data not shown). Fold expression (up-or down-regulated) of a miRNA in a tumor tissue compared to its adjacent control and number of samples providing expression data of a miRNA could be observed in Heat map diagrams of cluster of 13 samples (Figure 3b). It showed that expression of all miRNAs was not available from all samples. Expression of some miRNAs was up-regulated in most of the samples (e.g. expression of miR-31*, miR-7, miR-21 shown at the bottom of the figure) and expression of some miRNAs was down-regulated in most of the samples (e.g. miR-206, miR-1, miR-133a shown at the middle of the figure) (Figure 3b).
Out of these 30 miRNAs (Table 6), expression of 28 miRNAs showed similar expression deregulation pattern as it was already reported in earlier studies on cancer [10,[29][30]. Of the remaining two miRNAs, one miRNA (hsa-miR-770-5p) has not been associated yet with any cancer. Expression of remaining one miRNA (hsa-miR-211) has been deregulated in opposite direction in this study compared to previous reports on colorectal cancer [29][30]. Actually, 7 miRNAs, whose expression was significantly deregulated in the analysis of 18 samples, are a subset of these 30 miRNAs. However, expression of hsa-miR-411* was observed to be significantly deregulated in this study, but no previous report had shown its involvement with cancer. Again, it is known that hsa-miR-411 and miR-411* are originated from the same precursor miRNA and hsa-miR-411 has a strong association with cancer [31]. When we checked expression of mature hsa-miR-411 and hsa-miR-411*, a strong positive correlation was observed (r = 0.95) (Figure 3c) although expression of hsa-miR-411 was not found to be statistically significant. Similarly, we have observed significant expression deregulation of hsa-miR-135b* and hsa-miR-99a* and reports of association of cancer with hsa-miR-135b and hsa-miR-99a [32]. Searching in similar way, 1207 unique target genes were obtained for these 30 miRNAs. IPA and KEGG mapping were performed using these 30 miRNAs and their known 1207 target genes. In IPA analysis, three cancers, which were found to be associated with expression deregulation of a major subset of these miRNAs, were hypo pharyngeal, esophageal and head and neck cancer (p-values 1.71610 207 , 2.20610 207 and 1.02610 207 respectively). Interestingly, same biological signaling pathways (relevant to cancer) were reported to be involved with these set of 30 miRNAs as it was observed with 8 miRNAs earlier.
But the difference lies in the repertoire and number of targeted nodes for all these pathways (Table S1 in File S2, Table S2 in File S2, Table S3 in File S2 and Table S4 in File S2).

Discussion and Conclusions
Expressions of 7 miRNAs were significantly deregulated in 18 cancer samples and significant up-regulation of hsa-miR-1293 is being reported for the first time in gingivo buccal cancer ( Table 4). As of now, functional role of hsa-miR-1293 in cancer is very limited. According to miRWalk [21] and StarBase [33], one of the predicted targets of hsa-miR-1293 is MAPK14 (p38) which has been shown to be associated with tumor's sensitivity to cisplatinum treatment [34]. If this predicted relationship could be validated, then expression of this miRNA may be useful in prediction of patient's sensitivity to cis-platinum treatment. Two miRNAs, hsa-miR-206 and hsa-miR-1, are known to play antitumorigenic role [35][36][37] and hsa-miR-206 may indirectly activate 7 miRNAs whose expression was significantly deregulated samples. Each row represents expression of a miRNA and each column represents a sample. Sky-blue colored cells stand for failed assay i.e. no data in those cells. Red and green colors signify up-and down-regulation, respectively. Dendogram along the vertical axis represents hierarchical classification of miRNAs on the basis of expressions similarity. Distance metrics of hierarchical clustering was Euclidean distance. Heat map was constructed using Heatmap 2 of R's ''gplot'' package. C: Highly correlated expression of miR-206 and miR-1. doi:10.1371/journal.pone.0104839.g001 Table 4. Expression profile of 7 miRNAs in precancers and cancer samples.  apoptosis, inhibition of cell migration and focus formation [38]. Down regulation of expression of hsa-miR-1 and hsa-miR-133a, which has been observed in this study (Table 5), has already been reported in an earlier study with oral squamous cell carcinoma [39]. In fact, hsa-miR-133a targets several oncogenes and is reported to be commonly down regulated in a number of other oral cancer studies [35][36][37][38][39][40]. Expression of both hsa-miR-31 and hsa-miR-31* was reported in some earlier study on cancer [41][42]. Study on OSCC cell line showed that exogenous delivery of pre-mir-31, which boosts up quantity of mature hsa-miR-31 and hsa-miR-31*, enhanced OSCC oncogenicity [41][42]. Hence, our observation of up-regulation of hsa-miR-31 and hsa-miR-31*,   corroborates with existing reports on oral and other cancers. Similar to a previous report on head and neck cancer, here also, expression of hsa-miR-204 was also found to be significantly down regulated [43]. Expression of hsa-miR-7, a known OncomiR, is reportedly up-regulated in OSCC [44]. It targets primarily tumor suppressor transcripts from RECK [16]. In this study, expression of hsa-miR-7 was also significantly up-regulated and thus corroborates with previous findings related to oral cancer. Efforts have been made to understand possible biological implication by mining different databases. Pathway analysis revealed that most   disrupted biological processes would be cell migration, apoptosis and proliferation (Table S2 in File S2). Most disrupted pathways were predicted to be ''proteoglycans in cancer'' and PI3K-AKT (Table S1 in File S2) which are also supported by our RNA-Seq data on expression deregulation of some proteoglycan genes and LAMC and PAI which belong to PI3K-AKT pathway (unpublished data). These observations also support our predictive pathway analysis regarding involved biological process/pathways that could be targeted by this set of 8 miRNAs.
Cluster analysis of 18 samples had shown that expression of 30 miRNAs was found to be significantly deregulated in the cluster of 13 samples (Figure 3b). Literature search and database mining with these 30 miRNAs showed relevance of these genes with oral and other cancers (Table 6). Although, pathway analysis using 8 miRNAs (identified from expression data of 18 samples) and 30 miRNAs (identified from expression data of cluster of 13 samples) narrowed down to similar signaling pathways but the number and repertoire of nodes targeted by these set of miRNAs are quite different (data not shown). This cluster analysis is actually revealing existence of molecular heterogeneity that may chew up resolution of finding finer molecular marker. Interestingly, observed molecular heterogeneity pattern in no way related to differentiation subtype that exists within the tissues or clinical stages of samples.
Among two pre-cancers, leukoplakia is known to be associated with tobacco habit, whereas lichen planus is an auto immune disease. It is reported that all oral cancers are preceded by precancerous lesions. So, we checked whether, similar to cancer samples, expression deregulation of miRNAs could be observed in these two precancerous lesions. Here, TaqMan data showed that expression of only hsa-mir-31 was significantly up-regulated (4.55 fold more than adjacent control tissue) in leukoplakia samples (Table 4) but not in lichen planus tissue. Though, expression of another miRNA, hsa-mir-31*, was also up-regulated by 4.75 folds but significantly not different due to greater standard deviation among samples. So, expression of these two miRNAs needs to be checked in more leukoplakia samples. This once again reiterates about molecular proximity of leukoplakia with GBSCC than other pre-cancer, lichen planus. So, expression of hsa-mir-31 and hsa-mir-31* in leukoplakia tissues could be potential risk-markers for progression of precancer to GBSCC.
Small number and mixed stage of tumor samples and expression assay of only 762 miRNAs by TLDA (available at the time of our experiment) limited this study to infer that expression of only 7 miRNAs were significantly deregulated in GBSCC tissues compared to adjacent control tissues. There is high chance that number of deregulated miRNAs will increase if we could assay expression of ,2578 miRNAs presently known in human tissue as mentioned in miRBase-20 [45]. As a result, expression of more miRNAs could have been checked in cancer and leukoplakia tissues to infer more about molecular markers. More importantly, role of these miRNAs in carcinogenesis is to be validated by functional study.

Supporting Information
File S1 Genome wide miRNA expression profile for 528 miRNAs (excluding 3 endogenous controls) in GBSCC. Each column from S1 to S24 represents data for 18 samples, each row represents data for DDCt of a specific miRNA. N/A represents no data.