Abstract
Background. Oral squamous cell carcinoma (OSCC) is one of the most common head and neck squamous cell tumors. MicroRNAs and DNA methylation, as epigenetic mechanisms, regulate the expression of oncogenes and tumor suppressor genes, contributing to the carcinogenic development. However, the current knowledge on the genetic and epigenetic landscape of OSCC is still limited.
Objectives. To assess the transcriptomic impact of microRNAs found to be methylated through Infinium genome-wide methylation profiling of archived OSCC tissues, and to analyze their biological role using gene network analysis.
Materials and methods. We used the Infinium array-based methylation assay to assess the genome-wide methylation status at the single-CpG-site level of DNA purified from archived OSCC tissue samples. After quality control, filtering out poorly performing probes and normalization of data, we identified the differentially methylated microRNA loci. We performed a literature-based analysis of OSCC transcriptomic data to identify the predicted target genes for each microRNA, followed by individual network and pathway enrichment analyses.
Results. The analysis of Infinium methylation array data revealed 1469 differentially hypomethylated loci, 4 of which were of interest, namely hsa-microRNA-124-3, hsa-microRNA-24-1, hsa-microRNA-769, and hsa-microRNA-4500. Network and pathway enrichment analyses revealed multiple pathways modulated through DNA methylation-microRNA expression axes.
Conclusions. We describe the transcriptomic impact of 4 differentially methylated microRNAs in OSCC tissues samples and discuss their role in the pathology of OSCC. These results may contribute to a better understanding of how epigenetic mechanisms such as DNA methylation and microRNAs cooperate to impact the development of OSCC.
Key words: DNA methylation, microRNAs, network analysis, oral squamous cell carcinoma
Background
Oral squamous cell carcinoma (OSCC) is one of the most common head and neck squamous cell carcinoma (HNSCC) tumors, with over 350,000 new cases annually.1 Despite recent therapeutic progress, the cure rate of OSCC is relatively low and closely correlated with the tumor stage at diagnosis. Moreover, in recent years, there has been a change in the incidence of OSCC, with an increased occurrence in young people.2
Oral cancer is less studied than other cancers, and there is currently a lack of bona fide diagnostic and prognostic biomarkers. Similarly to other solid cancers,3, 4, 5 OSCC has a multifactorial and multigenic character. Its development is driven by various intricate genetic and epigenetic mechanisms which play important roles in the initiation and development of these tumors.6
The DNA methylation and non-coding RNAs are 2 of the most studied epigenetic mechanisms impacting global gene expression.7 The DNA methylation is essential for several important cellular mechanisms, including regulation of gene expression, differentiation, genome imprinting, and stability.8 The accumulation of alterations to genome methylation leads to the destabilization of cellular homeostasis. The aberrant DNA hypo- and hypermethylation patterns are critical players in tumorigenesis by promoting the expression of oncogenes and silencing tumor suppressor genes.9, 10, 11
MicroRNAs are small non-coding RNAs involved in the post-transcriptional regulation of approx. 60% of human genes. Multiple microRNAs have been associated with various aspects of oral cancer biology, including onset, metastasis, recurrence, prognosis, and survival.12, 13, 14 The interest for understanding the regulation of microRNA expression in the context of other epigenetic processes has been growing steadily, and there is a need to understand the global impact of microRNAs on cancer cell transcriptomics.
Network analysis approaches, such as protein–protein interaction networks, have the potential to add biological meaning to large-scale gene datasets such as those obtained from high-throughput RNA sequencing, DNA methylation assays or microRNA target gene analysis.15 Combining both biological networks and pathway enrichment is meaningful for predicting complex disease mechanisms, disease biomarkers, novel gene interactions, and drug design.16 An example of this is the integrative analysis of DNA methylation and microRNA expression data that has been conducted for hepatocellular carcinoma.17 However, a network analysis that integrates these 2 epigenetic mechanisms in the context of OSCC is lacking. Therefore, we investigated the potential interaction between these 2 epigenetic mechanisms, shedding light on unknown mechanisms that might contribute to the development of OSCC.
Objectives
The aim of the current study was to evaluate the impact of differentially methylated microRNAs in the OSCC transcriptome. Therefore, we used Infinium array-based methylation assay to analyze the genome-wide methylation status of microRNA loci in a set of OSCC-archived samples and correlated the results with publicly available microRNA expression data. We retrieved the OSCC transcriptome data from the Gene Expression Omnibus (GEO) database and used miRWalk v. 3.0 prediction algorithms (http://mirwalk.umm.uni-heidelberg.de/) to predict the gene targets of the differentially methylated microRNAs. The commonly predicted gene sets were analyzed further using network clustering coupled to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, in order to predict the biological impact of the differentially methylated microRNAs.
Materials and methods
Sample collection
The present case-control study was reviewed and the analysis of the provided formalin-fixed paraffin-embedded (FFPE) tissue samples was approved by the Ethics Committee of Victor Babeș University of Medicine and Pharmacy, Timișoara, Romania (approval No. SCUM TM 9414/30.05.2018). Due to the retrospective design of the current study and patient anonymity, the ethical review board determined that informed consent was not required. We used 13 FFPE samples from male OSCC patients (demographic characteristics are presented in Table 1); all samples were assessed and graded by 2 independent, trained pathologists. Adjacent tissues, both cancerous and normal, were dissected out of every FFPE block and pooled into 4 tumor (T) and 4 normal (N) samples for downstream methylation analysis.
DNA purification
Genomic DNA was extracted from the FFPE tissue samples using the QIAamp DNA FFPE Tissue Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instructions. Paraffin was dissolved in 1 mL of xylene and removed. Afterwards, the samples were lysed using 180 µL of Buffer ATL under denaturing conditions with 20 µL of proteinase K digestion. The samples were then incubated at 90°C to reverse the formalin crosslinking. To obtain RNA-free genomic DNA, 2 µL of RNAse A were added. The genomic DNA was bonded to the membrane using a mixture of Buffer AL and ethanol. In the washing steps, the residual contaminants were washed away using freshly reconstituted Buffer AW1 and AW2, according to the manufacturer’s instructions. Finally, the DNA was eluted in 50 µL of nuclease-free water. The quality and purity of the extracted DNA were measured using a spectrophotometer (NanoDrop™; Thermo Fisher Scientific, Wilmington, USA). The DNA purity of the samples, as measured using the 260/280 ratio, had an average of 1.89 ±0.18.
Genome-wide DNA methylation in the well-differentiated OSCC samples
We used the Infinium 450K array-based methylation assay (Illumina Inc., San Diego, USA) to quantitatively assess the genome-wide methylation status at single-CpG-site level, using the DNA purified from the pooled FFPE samples. The 1-μg genomic DNA samples were incubated in sodium bisulfite buffer (Sigma-Aldrich, St. Louis, USA) in a SimpliAmp Thermal Cycler (Applied Biosystems, Waltham, USA) with alternating cycles of 95°C and 60°C. The resultant DNA was bound again to a membrane, washed, incubated in a desulphonation buffer (Zymo Research, Orange, USA), and eluted with 20 µL of DNA elution buffer (Zymo Research). The quality and purity of the treated DNA was measured using a spectrophotometer (Nanodrop Technologies, Wilmington, USA). The yields were in the range of 50–100 ng/µL.
After bisulfite conversion, whole-genome methylation patterns were assayed using the Infinium® HumanMethylation450 BeadChip Kit (Illumina Inc.), following the manufacturer’s protocol. Briefly, the raw data were obtained in intensity data file format (*.idat) by scanning the processed chips using an iScan reader (Illumina Inc.). Each point of methylation data was presented with fluorescent signals from the methylated (M) and unmethylated (U) alleles. The methylation status of each CpG site was measured based on the ratio of fluorescent signals from one allele relative to the sum of both methylated and unmethylated alleles (β value) using GenomicRatioSet (Illumina Inc.). Probes with a detection p-value cutoff >0.05 were excluded from the subsequent analysis.
To minimize the variation between and within the samples, we applied the preprocess quantile normalization methon, and then performed a PCA analysis and generated multidimensional scaling (MDS) plots to visualize and explore the differences between OSCC and normal tissues.
The differentially methylated microRNA gene loci with an adjusted p-value < 0.05 were selected for further analysis.
Statistical analyses
Illumina methylation *.idat files were imported into R programminh language (R Foundation for Statistical Computing, Vienna, Austria) using the minfi function. The quality control of the imported samples was evaluated by calculating the detection p-value for every CpG island in each sample, which is indicative for the quality of the signal. Since the samples from our dataset were all derived from male donors, it was not necessary to filter the sites from X and Y chromosomes. Furthermore, we have removed the probes where common single nucleotide polymorphisms (SNPs) affected the CpG islands. Once the data were normalized and filtered, we visualized the results using MDS plots. In the final step, we calculated the β (percentage of methylation) and M (log2 ratios of the intensities of methylated compared to unmethylated probes) values.
The differentially methylated microRNAs loci were identified by a two-tailed, heteroscedastic Student’s t-test, being ranked according to their Bonferroni adjusted p-values (p = 0.05 as a cutoff).
Cross validation with literature analysis
To determine whether the differentially methylated microRNAs retrieved from our genome-wide DNA methylation analysis had been previously validated through experimental studies, we chose to perform a literature search in PubMed. In contrast to Web of Science and Google Scholar, PubMed is a human-curated database, that retrieves only peer-reviewed journal articles that fit into the searchable subject matter, with reliable sorting, reproducibility and reportable results.18 The key words combinations used in the PubMed literature search were: microRNAs AND oral cancer, microRNAs AND oral squamous cell carcinoma. We filtered out reviews, meta-analyses and articles that did not relate to OSCC and microRNAs.
Network analysis
Target predictions (3ʹUTR, 5ʹUTR and coding sequence (CDS) region) for hsa-miR-24-1, hsa-miR-124-3, hsa-miR-769, and hsa-miR-4500 were computed with the miRWalk v. 3.019 machine learning algorithm (Bonferroni adjusted p-value ≤0.05 as a cutoff), using the following criteria: 1) microRNA binding probability >0.99 and 2) microRNAs concurrently predicted by at least 1 other algorithm (miRDB or TargetScan). The microRNA target prediction datasets were cross-referenced with the differentially expressed genes (DEGs) retrieved after GEO2R analysis (Benjamini–Hochberg adjusted false discovery rate (FDR) <0.05) of the GSE138206 OSCC transcriptome dataset. The common genes of microRNA targets that were differentially expressed in OSCC transcriptomes were then subjected to KEGG pathway enrichment analysis and cluster analysis on the STRING platform20 using the following criteria: adjusted FDR value of 0.05 as cutoff and high confidence interaction score set at 0.700.
Prediction of CpG islands
MethPrimer platform (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi) was used for designing the bisulfite sequencing polymerase chain reaction (BSP PCR) primers and for predicting the CpG islands associated with the 4 microRNAs, using the 2-kb DNA sequence upstream of the microRNA loci as input data. For BSP primer selection and CpG island prediction, the general parameters have been set according to the recommendations by Li and Dahiya.21
Results
The β values represent the ratio between the methylated probe intensity and the overall intensity. The β-value statistics result in a number between 0 and 1, where 0 indicates that all copies of the CpG islands from analyzed samples are completely unmethylated, and a value of 1 indicates that every copy of CpG sites is methylated.22 In our case, the β-values results are situated close to 0 and are have higher density for tumor samples in contrast to normal samples. Conversely, the M-values are calculated as the log2 ratio of methylated compared to unmethylated probes. The M-values close to 0 indicate that CpG sites are approximately half-methylated, negative values mean that there are more unmethylated CpG sites than methylated ones, and positive M-values mean the opposite.22 In our case, the highest density plot of M-values close to 0 was created for the tumor group in contrast to the normal group. Therefore, the analysis of the Infinium methylation array β-values and M-values indicated that there is a higher degree of hypomethylation globally in the tumor group than in the control group (Figure 1).
To find the statistically significant genes that encode for microRNAs, we filtered the 440,364 gene loci assayed on the Infinium methylation array based on an adjusted p-value < 0.05. A total of 1469 loci were obtained, of which only 4 proved to encode microRNAs, namely hsa-miR-124-3, hsa-miR-769, hsa-miR-24-1, and hsa-miR-4500. Compared to controls, the microRNA loci from OSCC samples are hypermethylated, with log fold change (logFC) values between 3.21 and 4.99. The genomic locations of all 4 microRNAs are presented in Table 2.
To determine whether our DNA methylation array results are consistent with other published findings, we searched the GEO database the key words “OSCC” and “DNA methylation”. We retrieved 9 datasets, of which GSE123781 and GSE87053 contained the 4 loci of interest and were further selected for GEO2R analysis. We found that 3 out of the 4 gene loci were significantly hypermethylated, namely cg19267861 (the locus entailing hsa-miR-124-3), cg02065637 (the locus containing hsa-miR-124-3) and cg26110710 (the locus coding for SLITRK5 and hsa-miR-4500) (Table 3).
We used MethPrimer to design BSP PCR primers and predict the presence of CpG islands upstream of genes encoding the 4 microRNAs (Figure 2, Table 4). The CpG island prediction identified the presence of CpG islands upstream of hsa-miR-124-3 (Figure 2A) and hsa-miR-769 (Figure 2B), but not for hsa-miR-4500 and hsa-miR-24-1.
Five sets of BSP primers with their sequences were retrieved for all 4 microRNA molecules, each of them consisting of 2 pairs of forward and reverse primers (Supplementary File 1).
Our PubMed search for microRNAs associated with OSCC led to a total of 189 microRNAs, of which 101 are downregulated, and 88 are upregulated in cancer tissues compared to normal tissues (Supplementary File 2,3). Of note, hsa-miR-24-1 was found to be upregulated,23 while hsa-miR-76924 and hsa-miR-124-325 were downregulated in OSCC samples compared to controls. To date, no published data link hsa-miR-4500 to oral cancer.
Next, we used miRWalk v. 3.0 to predict the targets of the 4 microRNAs. By filtering out the genes with a binding probability below 0.99, we retrieved 4147 unique predicted targets for hsa-miR-24-1, 6084 predicted targets for hsa-miR-124-3, 743 targets for hsa-miR-4500, and 7256 targets for hsa-miR-769.
To understand the biological significance of these predicted targets in the specific context of oral cancer, we cross-referenced them with the differentially expressed genes in the GSE138206 OSCC dataset (after GEO2R analysis, adjusted p-value < 0.05). The GEO2R analysis of the GSE138206 dataset retrieved 2289 genes with dysregulated expression, of which 1490 were downregulated (FC < −0.538) and 799 were upregulated (FC > 1.0). Since hsa-miR-24-1 has a documented upregulated expression in OSCC,23 we chose to compare its predicted target genes with the GSE138206 list of downregulated genes, and found 261 putative common targets. The hsa-miR-76924 and hsa-miR-124-325 have both previously been reported to be downregulated in OSCC. Therefore, we compared their predicted targets with the list of upregulated target genes from the GSE138206 dataset and found 274 putative common targets for hsa-miR-124-3 and 298 for hsa-miR-769. To date, there are no published data linking hsa-miR-4500 to OSCC; therefore, we compared its predicted targets with both the upregulated genes (and found 22 common genes) and downregulated genes (and found 51 common genes).
The STRING functional analysis of the hsa-miR-24-1 putative targets led to a protein–protein interaction (PPI) network characterized by 259 nodes, 68 edges and a PPI enrichment p-value of 3.75e−12, for which no significant signaling pathway could be found after KEGG functional enrichment analysis.
The STRING functional analysis of the hsa-miR-124-3 putative target genes led to a PPI network characterized by 266 nodes, 562 edges and a PPI enrichment p-value below 1.0e−16. The average local clustering coefficient of the network is 0.359 with an average node degree of 4.23. The cluster analysis using the Markov Cluster Algorithm (MCL) method (inflation parameter of 4) led to the identification of 3 clusters functionally enriched through KEGG pathway analysis (referred to as Cluster 1, 2 and 3) containing 23, 12 and 7 genes, respectively. (Figure 3, Table 5). Cluster 1 (23 genes) is enriched with genes involved in viral infections, and the local clustering coefficient is high with a value of 0.84 and an average node degree of 10.8. Cluster 2 (local clustering coefficient = 0.769, average node degree = 7.83) includes 12 genes related to protein digestion and absorption, ECM-receptor interaction pathway and focal adhesion. Cluster 3 contains 7 genes significantly associated with the p53 signaling pathway, with the local clustering coefficient of 0.867 and an average node degree of 2. The high values of the local clustering coefficients suggest that members within these clusters are involved in similar pathways and functions.
The STRING functional analysis of the hsa-miR-769 targets identified a PPI network characterized by 290 nodes, 581 edges, and a PPI enrichment p-value below 1.0e−16. The average local clustering coefficient of the network was 0.357 and the average node degree was 4.01. The MCL cluster analysis of the network also resulted in 3 clusters functionally enriched in KEGG pathways, containing 25, 13 and 8 genes, respectively (Figure 4, Table 6). Cluster 1 contains 25 genes that were also enriched in viral infection pathways, with a local clustering coefficient of 0.809 and an average node degree of 10.6. Cluster 2 (local clustering coefficient = 0.79, average node degree = 7.54) includes 13 genes which are most significantly enriched in protein digestion and absorption, ECM receptor interaction and AGE-RAGE signaling pathways. Cluster 3 contains 8 genes significantly enriched in AGE-RAGE signaling pathway, ECM-receptor interaction and cancer-related pathways. The local clustering coefficient for Cluster 3 was 0.893, with an average node degree of 2.5.
Due to the limited number of hsa-miR-4500 target genes found in the OSCC transcriptome, the resulting PPI network was very poor and lacked statistical significance. Therefore, this was not pursued further.
Discussion
Oral squamous cell carcinoma (OSCC) is one of the most common types of HNSCC. Therefore, it is critical to gain a better understanding of the epigenetic mechanisms that underlie OSCC initiation, development, evolution, and response to therapy.26 In patients with OSCC, increased methylation status of microRNA loci has been observed compared to controls, suggesting that this parameter may serve as a predictor for malignant transformations.27 In the present study, we used archived FFPE samples to explore the OSCC genome-wide methylation profile of microRNA loci.
Our analysis led to the identification of 4 microRNA loci significantly hypermethylated in OSCC tissue compared to normal adjacent tissue, namely hsa-miR-24, hsa-miR-769, hsa-miR-124, and hsa-miR-4500.
The hsa-miR-24 is one of the most heavily studied microRNAs in cancer, being associated with tumor initiation, progression, metastasis, and response to therapy in a wide array of malignant pathologies.28 In our study, the hsa-miR-24 loci were significantly hypermethylated both at the 3ʹUTR site and 200 bp upstream of the transcription start site (TSS). There is evidence that if DNA methylation occurs within the 3ʹUTR region, it can both positively and negatively impact the transcription.29 Therefore, the hypermethylation status of the hsa-miR-24-1 locus may still correlate with the reported upregulated expression of miR-24-1/-3. Of note, in OSCC, the mature miR-24-1/-3 is overexpressed not only in malignant tissues23 (where it correlates with the staging of the disease) but also in the patients’ saliva, serum and plasma, qualifying it as a possible diagnostic and predictive biomarker.30, 31, 32
We also found that the hsa-miR-124-3 locus is methylated both distally (1500 bp upstream, TSS1500) and in close proximity (200 bp upstream, TSS200) to the TSS, which suggests gene silencing. This is also consistent with reports showing the downregulation of hsa-miR-124-3 in an OSCC cell line25 and animal models of oral cancer.33 The KEGG pathway enrichment analysis results for the hsa-miR-124-3 gene targets in OSCC were related to viral infections in the most populated gene cluster. Hepatitis C, Epstein–Barr virus, herpes simplex 1 and human papillomavirus are among the most significant viruses associated with this set of genes. The links between the etiologic causes of OSCC and viruses have long been investigated.34, 35, 36 Hepatitis C virus was demonstrated to be involved in the development of oral lichen planus, a chronic inflammatory condition that increases the risk of OSCC.37 Epstein–Barr virus,38 herpes simplex virus 139 and human papillomavirus40 have also been investigated in the promotion of the malignant transformation of cells in OSCC. In Cluster 2, pathways such as protein digestion and absorption, and ECM-receptor interaction pathway proved to be the most relevant ones. The presence of ECM-receptor interaction pathway associated with the Cluster 2 gene set confirms previous literature findings reporting the involvement of hsa-miR-124-3 in the ECM receptor interaction pathway through its direct repression of ITGB1 in OSCC cell lines.25
Regarding hsa-miR-769, we found the hypermethylation of both 1500 bp upstream of the TSS and within the miR-locus. Our data are consistent with the documented downregulation of hsa-miR-769 expression in OSCC tissues, since hypermethylation of its TSS may block the transcription.24 The KEGG pathway enrichment analysis of Cluster 1 identified Epstein–Barr virus infection, hepatitis C and influenza A as the most significant pathways associated with hsa-miR-769 target genes that were differentially expressed in OSCC. In this case, the ECM-receptor interaction was also retrieved from among the significant pathways in Clusters 2 and 3. Interestingly, this pathway has previously been predicted to be associated with OSCC.41 A new entry for this set of genes is the AGE-RAGE signaling pathway retrieved for Clusters 2 and 3. This pathway was initially reported to be responsible for the initiation of diabetic complications.42 The receptor for advanced glycation end products (RAGE) molecule, a constituent of this pathway, has been reported to be a key factor that accelerates tumor progression and metastasis in various malignancies.43, 44, 45 Regarding OSCC, there have been several investigations that aimed to evaluate the role and expression of RAGE receptor.46, 47, 48 In a study that evaluated 74 OSCC samples, the RAGE receptor was reported to be a prognostic factor, and was significantly associated with the depth of invasion and local recurrence.46 Moreover, a particular role of RAGE in OSCC pathogenesis was confirmed in an in vitro study, where this receptor proved to be activated in a paracrine manner by the extracellular molecule high-mobility group box 1 (HMGB1) that was secreted by the OSCC cells. This activation promotes tumor progression and stimulates bone destruction.47 Furthermore, RAGE proved to be significantly decreased by evodiamine, a novel anti-tumor drug, inhibiting proliferation, invasion and angiogenesis of OSCC both in vitro and in vivo.48 Finally, the gene that encodes the RAGE receptor is predicted to be one of the hsa-miR-769-5p targets both in CDS and 5ʹUTR positions, according to the miRWalk v. 3.0 database.19
Finally, for hsa-miR-4500, there is no evidence so far concerning its expression in OSCC. When performing the KEGG pathway functional enrichment analysis for the commonly predicted target genes of this microRNA, there was no significant pathway associated with it due to a limited number of resulting genes.
Limitations
Due to financial limitations, we were not able to validate our methylation results in additional independent cohorts. Caution must be taken given the small sample size and the fact that the findings might not be entirely transferable to a larger population. Further investigations are needed to validate the correlation between the expression and methylation status of hsa-miR-24-1, hsa-miR-124-3 and hsa-miR-769 in OSCC samples through quantitative real-time (qRT)-PCR.
Conclusions
Returning to the study’s central question, we suggest that, in stark contrast with the overall genome DNA hypomethylation, microRNA loci could show increased methylation status, which may translate into differential expression of their mature forms. Of note, one of the identified microRNAs has been proposed as an OSCC biomarker, given its presence in multiple biological fluids. Our bioinformatics analysis also suggests that methylation-dependent modulation of microRNA expression could significantly impact the response pathways to viral infections, the ECM-receptor interaction pathway and AGE-RAGE signaling pathway. These may be relevant for the development, progression and pathogenesis of OSCC.
Supplementary data
The supplementary Tables are available at https://doi.org/10.5281/zenodo.6759729.
The package contains the following files:
Supplementary File 1. Primer picking results for bisulfite sequencing PCR (BSP);
Supplementary File 2. Downregulated microRNAs curated from literature analysis;
Supplementary File 3. Upregulated microRNAs curated from literature analysis.
The references used for this supplementary data are listed at the end of the files.
