Prognostic Value of MiR-21: An Updated Meta-Analysis in Head and Neck Squamous Cell Carcinoma (HNSCC)

Head and neck squamous cell carcinoma (HNSCC) is a group of malignancies with serious impact on patient quality of life due to a reduced rate of response to chemotherapy or radiation therapy. MiR-21 has been identified as one of the most common proto-oncogenes. It is hypothesized that upregulated miR-21 could serve as a potential biomarker for human cancer diagnosis. Considering the target genes identified for miR-21 in HNSCC, this transcript is an important player in several cellular processes that control carcinogenesis. The abnormal expression of miR-21 in this group of pathologies has been assessed in several publications, but given the heterogeneity of the published results, a meta-analysis and proper bioinformatics analysis of expression databases are needed to correctly establish the prognostic potential of this molecule. The present meta-analysis comprises the published survival data on HNSCC patients, reported as HR and 95% CI, in association with the expression levels of miR-21. Our investigation revealed that miR-21 could be used successfully as a prognostic biomarker in HNSCC patients, confirming its oncogenic potential. Specifically, the upregulation of miR-21 in these patients predicts a worse outcome in terms of survival rate.


Introduction
Head and neck cancer (HNC) is considered a pathology with severe consequences on the quality of life of patients, with a reduced rate of response to treatment [1,2]. HNC is a highly heterogeneous group of tumors occurring on the mucosal surfaces of the nasal and oral cavity, oropharynx, and larynx [1]. From this heterogeneous group, over 90% of diagnosed cases are head and neck squamous cell carcinomas (HNSCCs), statistically being the sixth most widespread cancer worldwide [1].
The survival rate remains low, since up to 25% of the patients progress to metastatic disease. The presence of the lymph node metastases (LNM) represents an important determinant factor in HNSCC, considering the positive correlation between the lymphatic metastatic spread and the survival rate [3]. The main risk factors are related to smoking and alcohol consumption [4,5]. Also, infection with high-risk human papillomaviruses (HPVs) has recently been associated with the pathogenesis of HNSCC [6].
Considering the complexity of the HNSCC molecular mechanism, most studies have focused on coding and noncoding genes. MicroRNAs (miRNAs) are noncoding transcripts, 18-22 nucleotide long RNAs, being vital posttranscriptional regulators of gene expression [7]; at the cellular level they might act as oncogenes or tumor suppressor genes affecting cell proliferation, invasion, migration, self-renewal capacity, resistance to therapy, and autophagy [8,9].
This noncoding RNA is upregulated in various human cancers [10,11], including HNSCC [12]. MiR-21 regulates the expression of multiple target genes with fundamental implications in cancer progression [13]. In oral squamous cell carcinoma, this microRNA is part of the mechanism that regulates sensitivity to cisplatin, by targeting phosphatase and tensin homolog (PTEN), as well as programmed cell death protein 4 (PDCD4) genes. Moreover, this feature is transmissible through the presence of exosomes carrying miR-21 to other cells [14], hence supporting previous findings [15]. Comparable results were reported by Yan et al., 2018, in salivary adenoid cystic cancer cell lines, with miR-21 being involved in tumor cell growth via the regulation of PDCD4 expression and B-cell lymphoma 2 (Bcl-2) genes. In addition, miR-21 has been shown to target the PTEN gene, being involved also in the ability of salivary adenoid cystic carcinoma cells to metastasize [16]. Another study observed that the overexpression of miR-21 and cyclin-dependent kinase 5 (CDK5) were associated with epithelial-to-mesenchymal transition (EMT) and lymph node metastasis (LNM). The results suggested that miR-21/CDK5 interaction impacted EMT by targeting signal transducers and activators of transcription 3 (STAT3) in HNSCC cells [12]. In tongue squamous cell carcinoma, upregulation of mir-21 could independently predict shorter survival of patients by modulating the expression of the tropomyosin 1 (TPM1) gene and inhibiting the apoptosis of tumor cells [17]. Furthermore, there is evidence regarding the modulation of radiosensitivity through miR-21; it was observed that the downregulation of miR-21 increased the radiosensitivity of nasopharyngeal cell lines [18]. All these investigations emphasize the role of miR-21 in the progression of head and neck pathologies, but also in resistance to treatment, given that radiotherapy is used as the primary treatment along with surgery, as it is confined to local effects [19].
Several studies have focused on evaluating the expression levels of this transcript; however, inconsistencies or heterogeneities concerning the prognostic value of miR-21 remain. Our work focused on performing a meta-analysis in order to assess the expression levels of miR-21 and to accurately evaluate its prognostic value in HNSCC. We also did a statistical analysis of various databases to gain a more in depth understanding of the miR-21 expression across various experiments and of the significance of its inhibitory activity.

Search Strategy and Eligibility Criteria for Meta-Analysis
Following the guidelines and criteria of Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA), we did a systematic search directed by interrogating the PubMed and Embase databases. Our search approach assumed the following terms: ('otolaryngolog' OR 'head and neck' OR 'buccal' OR 'mouth' OR 'oral cavity' OR 'larynx' OR 'lip' OR 'hypopharynx' OR 'oropharynx' OR 'tongue') AND 'patients' AND 'survival' AND ('miR-21' OR 'miRNA 21' OR 'microRNA 21'). This strategy was applied for both databases. Studies that investigated the prognostic value of miR-21 expression in the tissue/serum/plasma of HNSCC patients and that were published before May 2019 were assessed for eligibility. We chose only publications that reported hazard ratios (HR) and confidence intervals (CI) for the detailed investigation.
We excluded articles that met any of the following criteria: (1) conference/meeting abstracts; (2) not published in English; (3) full text not available; (4) written as review/meta-analysis; (5) the expression of miR-21 was not evaluated in human samples; (6) did not offer sufficient data regarding the HR and 95% CI for the survival analysis.

Search Strategy in Open Databases
We searched the ArrayExpress database with the terms: "miR-21, head and neck cancer" or "miRNA head and neck cancer". The study selection was done accordingly: (1) inclusion of human subjects; (2) expression analysis done on normal and tumor tissue; (3) the localization of tumors must include the head and neck area; (4) the expression analysis should include miR-21 expression analysis; (5) exclusion of unprocessed datasets; (6) datasets with at least 10 samples. Addition datasets were included based on a similar study done on miR-375-5p in head and neck cancer [20]. We also searched the TCGA database and retrieved the miRNA Expression Quantification (n = 569) GDC Hub.
In order to assess the expression level of these genes, especially the genes targeted by miR-21, we retrieved data from the GDC and TCGA Head and Neck Cancer (HNSC) datasets, which included HTSeq-Counts (n = 546).

Quality Appraisal
The eligibility criteria were separately employed by two individual reviewers who then proceeded to data extraction (C.P.-B. and R.C.). Studies were selected for analysis if they met the following benchmarks: (1) the expression of miR-21 was determined in patients with HNSCC; (2) the assay used to quantify the expression was clearly defined; (3) the outcome of patients was expressed as overall survival (OS), relapse-free survival (RFS), cancer-specific survival (CSS), disease-free survival (DFS) or 5-year survival; (4) the HRs of the obtained miR-21 are clearly explained as high vs. low/low vs. high.

Data Extraction
The information retrieved from the selected articles was organized in an excel worksheet that included data about publication (DOI, author's name), population (country, smoking status, age), sample type and count, tumor data (anatomical site, staging, lymph node invasion, human papillomavirus (HPV) status, experiment (assay, miRNA expression), and statistical analysis (HRs, 95% CI, p-value, survival type, direction of result).

Data Analysis
Analyses were conducted by using Comprehensive Meta-Analysis software, version 2.2.050 (Biostat Inc., Englewood, NJ, USA). As an indicator of effect sizes, the hazard ratio was used. Given the nonsignificant heterogeneity of the studies, all analyses were based on a fixed effects model. Publication bias analysis was performed using two methods. First, we used the Classic Fail-safe N of Rosenthal, in order to evaluate the number of studies with null effects that are required to turn the effect size to zero. Second, we calculated the Begg and Mazumdar's rank correlation test. This test computes the rank order correlation (Kendall's tau b) between the effect size and the standard error (which is driven primarily by sample size). This determines whether large studies tend to be included in the analysis, regardless of their effect size, whereas small studies are more likely to be included when they show a relatively large effect size. For the moderation analysis, we performed statistical comparisons between the categories of each moderator (the case of categorical moderators) and meta-regressions for continuous moderators.
The miRNA expression data retrieved from the TCGA and GEO databases were analyzed in GraphPad Prism version 8.2.1(GraphPad Software, San Diego, USA). We first separated the normal tissue from the tumor tissue, and then applied an unpaired two-way T-test for each dataset.
The miRNA targets were discovered by searching the TargetScan online tool, which gave us 384 results. At the same time, we searched the GEPIA2 web server that works with standard processed data from the TCGA and GTEx projects. For the differential expression analysis, the Log 2 FC was set at 1 and the q-value cutoff at 0.01. The chosen differential expression method was ANOVA. We then found the common genes for both groups and constructed an interaction network in the String database (https://string-db.org/). The minimum required interaction score was set at low confidence (0.150) and clustered with the MCL clustering with an inflation parameter of 2. The gene functional profiling was done in g:Profiler (https://biit.cs.ut.ee/gprofiler/gost). The Kaplan Meier survival plots for gene expression were calculated for the overall survival Cutoff-High (%) of 75% and Cutoff-Low (%) of 25%, and calculation of hazard ratio based on Cox PH Model. The Pearson correlation between genes and miR-21 was calculated in GraphPad Prism.

Database Screening and Characteristics of the Included Studies
The identification and selection trial are briefly illustrated in Figure 1. The online search in PubMed and Embase retrieved 84 articles for miR-21 and the HNSCC group of cancers. After the removal of duplicates, we manually excluded conference abstracts and articles with no full-text available. From the remaining 59 articles, we eliminated those that were irrelevant for the analysis, meaning reviews/meta-analysis publications, experiments that were not conducted on patients, and those that were not available in English. Eventually, we removed 11 more articles, as survival analyses was not reported as HR and 95% CI. The current meta-analysis collected seven articles that investigated the prognostic value of miR-21 in head and neck cancers [21][22][23][24][25][26][27].
The basic aspects in the selected studies are represented in Table 1. The present meta-analysis included a total number of 757 patients from six countries (Brazil, Denmark, Korea, Serbia, Taiwan, and USA). The investigations were focused on HNSCC localized in the oropharynx, hypopharynx, larynx, oral cavity, and tongue. All studies used tissue samples to assess the expression of miR-21, and the expression quantification assays included quantitative reverse-transcription PCR (qRT-PCR) and microarray or in situ hybridization (ISH). The involvement of miR-21 in the outcome of HNSCC patients was examined as DFS, OS, RFS, CSS, and 5-year survival.   Figure 2 describes the forest plot for miR-21 with the hazard ratio for each study included in the meta-analysis. As shown, from a total of seven studies included in the analysis, five revealed significant hazard ratios [22][23][24][25][26] indicating high expressions of miR-21 associated with poor survival, while two studies obtained nonsignificant hazard ratios [21,27]. The overall corrected hazard ratio is HR = 1.719, 95% CI = (1.402, 2.109) and statistically significant (Z = 5.198, p < 0.001), indicating that a high expression of miR-21 is associated with poor survival for patients diagnosed with HNSCC.

Data Analysis
The heterogeneity analysis performed for the distribution of effect sizes in our meta-analysis indicated nonsignificant heterogeneity, Q (6) = 3.078, p = 0.799. As a consequence, all the data analysis was performed using a fixed-effects model.

Classic Fail-safe N
In order to test the publication bias, first, we used the Classic Fail-safe N of Rosenthal. The concern of publication bias is that some nonsignificant studies are missing from the analysis, and that these studies, if included, would nullify the observed effect. If this number is relatively small, then there are reasons for concern. However, if this number is large, it generates confidence that the effect size, while possibly inflated by the exclusion of some studies, is nevertheless not null.
Our meta-analysis incorporates data from seven studies, which yields a z-value of 5.015 and a corresponding 2-tailed p-value of 0.0001. The fail-safe N is 39. This means that we would need to locate and include 39 'null' studies in order for the combined two-tailed p-value to exceed 0.05. In other words, there would have to be 5.6 missing studies with a null effect for every observed study for our effect to be nullified.

Begg and Mazumdar's rank correlation
Our second publication bias analysis was performed using the Begg and Mazumdar's rank correlation test. This test is concerned with the potential relationship between the size of the studies and the effect size obtained by each one. The results revealed a nonsignificant Kendall's tau b of 0.095, with a two-tailed p-value of 0.763 (based on continuity-corrected normal approximation), suggesting there is no tendency for studies which are more precise (and implicitly larger) to generate larger effect sizes.

Moderators' Analysis
In our moderation analysis, we tested two categorical moderators (type of sample: FFPE tissue vs tissue, technique: ISH vs. Microarray/qRT-PCR) and five continuous moderators (percentage of women, percentage of smokers, percentage of HPV-16 positive patients, percentage of stage III/IV patients, and percentage of patients who have undergone therapy in the follow-up time, counted as data reported in the studies analyzed).

Categorical moderators Sample Type
As shown in Table 2, the analysis of hazard ratios as a function of sample type revealed significant HR for both types of samples analyzed, HR = 1.771, 95% CI = (1.313, 2.391), for FFPE tissue and HR = 1.675, 95% CI = (1.267, 2.215) for tissue, with no significant differences between those two categories, Q (1) = 0.071, p = 0.789.

Continuous moderators
In order to test continuous moderators, we performed a metaregression in which the percentage of women, smokers, HPV-16 positive patients, stage III/IV patients, and patients who had undergone treatment after the sampling, were treated as predictors of hazard ratios of miR-21.

Percentage of Women
Among the studies analyzed, only six reported the proportion of women. We performed a meta-regression analysis in order to test the predictive value of the proportion of women in the samples upon the hazard ratios. The results revealed that this variable did not have a significant predictive value, B = −0.004, p = 0.674.

MiR-21 Expression Across Various Datasets
In our final analysis, we included a total of six independent datasets: five retrieved from GEO, which included microarray analyses done on Agilent or the Affimetrix platform, and one from TCGA, which was done with the help of the miRNASeq technique. The studies included a total number of 969 tumor samples and 101 normal tissue samples. In all cases, the unpaired two-tailed T-test showed p < 0.001. In five of the six included databases, miR-21-5p was upregulated in tumor tissue in comparison with normal tissue. In one dataset, however, GEOD-32960, miR-21-5p was down-regulated in tumor tissue. This might be caused by a very small number of normal tissues provided (Figure 3).

miR-21 HNSCC targets are involved in angiogenesis and RUNX expression
miR-21-5p was discovered to have 384 target genes; a differential gene analysis done in GEPIA2 revealed that there are 2077 under-expressed genes in HNSCC tumor tissue versus normal tissue. From these two lists of genes, 28 genes were found to be common ( Figure 4A). The list of common genes is presented in Table 3. The STRING protein network of interactions between the 28 common genes; the minimum required interaction score was set at low confidence (0.150) and clustered with the MCL (Markov Cluster Algorithm) clustering with an inflation parameter of 2. (C) Gene annotation of the miR-21 potential targets in HNSCC together with associated biological function with the g:Profiler online tool. The genes are involved in vasculature development, cardiovascular system, cell surface receptor signaling pathway, blood vessel development, tube development, circulatory system development, regulation of locomotion, angiogenesis, and anatomical structure formation. Some of these genes also regulate commissural axon pathfinding, and RUNX1 (Runt-related transcription factor 1) expression and activity.

Targeted gene expression correlation with overall survival
From the 28 common genes, only 5 showed a statistically-significant correlation with the overall survival rate of HNSCC patients (Log rank p < 0.05). These genes are NETO2, CDK6, OLR1, and NTF3, TGFBI ( Figure 5A). HNSCC patients with low expression of NETO2, CDK6, OLR1, and TGFBI have a better overall survival rate than HNSCC patients with high expression of these genes, while patients with low expression of NTF3 have a worse overall survival rate than those with high expression of NTF3. This means that even though TCGA data shows that NETO2, CDK6, OLR1, and TGFBI are underexpressed in HNSCC tissue compared to normal adjacent one, patients who show a low expression of this gene have a better disease outcome with regard to overall survival rate. Therefore, the data highlights a more complex signaling mechanism, which is not limited only to the interaction between miR-21 and the predicted target genes in HNSCC. The Pearson correlation analysis between the expression of all 28 genes from Table 3 and miR-21 expression shows that within the five genes with an impact on the overall survival of HNSCC patients, four have positive correlations with miR-21 (ORL1, TGFBI, NETO2, and CDK6) (see Figure 5B). The data might suggest that these are not true miR-21-5p targets in HNSC, because the expression of these genes has a dependent increase with the increase of miR-21-5p expression. The AIF1L gene showed the most negative r value of correlation with miR-21; however, this gene does not have a significant impact upon the overall survival of HNSCC patients. Only one gene passed all the evaluations: NTF3. This gene is underexpressed in HNSCC tumor tissue versus normal adjacent tissue, it is targeted by miR-21-5p, according to the TargetScan online tool, its low expression in tumor tissue reduces the overall survival rate of HNSCC patients, and it has a negative correlation with miR-21 expression.

Discussion
Accumulating data suggest that aberrantly-expressed miR-21 is a mark of malignant phenotype in several cancers [28], including HNSCC [29], where, even if the pathologies forming this group are very different in terms of histology or anatomic site, this microRNA is upregulated in the vast majority of cases. There are no disputes regarding the importance of miR-21 as a cancer-related transcript. Furthermore, there are studies showing that a combination of miR-21 and miR-375 expression profiles has a great potential to discriminate patients with oral tongue squamous cell carcinoma from normal controls [30]. The present meta-analysis brings direct confirmation of the oncogenic role of miR-21 in head and neck cancers. Using fixed-effects model, the investigation showed that miR-21 upregulation predicts shorter survival in HNSCC patients. We also found that the expression of miR-21 in the TCGA and GEO datasets is upregulated in the majority of cases, and that the underexpression of some of its targets is associated with worse outcome in terms of survival rates. Our results confirm the oncogenic potential of miR-21 in head and neck cancers, given the fact that in this category of cancers, the anatomic sites include many different histological types, and the upregulation of miR-21 remains significant and allows the separation of these patients' outcomes. A comparison between reported continuous variables that could influence the predictive value of miR-21 in these patients was not possible, as the selected studies did not report separate HRs for each group, i.e., males vs. females, smokers vs. nonsmokers, stage I/II vs. III/IV, HPV-16+ vs. HPV-16-, and treated vs. untreated patients; therefore, we investigated the moderator effect-evaluated as proportions in the cohort-of these variables on the predictive value of miR-21 in HNSCC patients. When continuous moderators were assessed, our analysis revealed that the prognostic value of miR-21 expression remains stable, that it is not moderated by the percentage of women, smokers, stage III/IV patients, and HPV-16 status, and that patients who underwent conventional chemo-or/and radiotherapy in the sample did not influence the predictive value of miR-21. Given that the studies reported the proportion of patients that followed a treatment course after the moment of the surgical resection of the tissue sample, these results should be interpreted with caution. Therefore, all these investigations confirm the reliable nature of miR-21 in assessing the outcome of HNSCC patients.
The upregulation of miR-21 was evaluated in HNSCC lines, revealing its involvement in cellular growth and cell cycle regulation, and its function as a regulatory molecule being exerted by complementary binding to the 3' UTR of PDCD4, with similar observations in HPV+ and HPV-cell lines [31]. It has been shown that in oral cancer, miR-21 is involved in the regulation of proliferation in cancer cells by inhibition of tumor necrosis factor alpha (TNF-α), and in the upregulation of miR-21, resulting in increased proliferation and decreased TNF-α expression [32]. The oncogenic potential of miR-21 in oral squamous cell carcinoma (OSCC) cell lines was underlined in relation with PTEN protein levels. Specifically, the induced underexpression of miR-21 correlated with the upregulation of PTEN gene expression and inhibition of cell growth, with these results showing miR-21 to be a potential therapeutic target in OSCC [33]. MiR-21 was also reported as a regulator of the immune response in tongue cancer, playing a role in expressly activating the nuclear factor NF-κB pathway in CAL27 cell line exosomes [34]. In salivary adenoid cystic carcinoma, by targeting genes such as PDCD4, PTEN, and BCL-2, miR-21 regulates the invasion and metastasis potential of the SACC-LM cell line [16]. In nasopharyngeal carcinoma, in vitro studies pointed out that the mechanism by which miR-21 controls the apoptosis, invasion, and cellular growth involves the targeting of the PTEN gene, thus interfering with the Akt pathway and BCL-2 expression [35,36]. The studies also revealed details related to the mechanism of miR-21 regulation. Precisely, the promoter of miR-21 requires the direct binding of STAT3 protein in order to initiate the transcription of this molecule [35]. Resistance to therapy is also mediated by miR-21, as, in tongue cancer, it was reported that in the Tca8113 and SCC-25 cell lines, transfection with miR-21 mimic resulted in enhanced chemoresistance capacity coupled with reduced apoptosis. MiR-21 interferes with resistance to chemotherapy by binding to cell adhesion molecule-1 (CADM1) mRNA [37]. Figure 6 summarizes the mechanisms behind the oncogenic role of miR-21 in malignancies from the head and neck area, according to the data presented in the scientific literature. Our bioinformatics analysis supports the oncogenic nature of miR-21 via the targeting of NTF3 and the association with a decreased survival rate in HNSCC patients. The association between human papillomaviruses (HPV), and particularly HPV-16 and NCSCC, is widely acknowledged, its carcinogenicity being proven previously [38,39], but also recognized by the International Agency for Research on Cancer (IARC) in 2007 as a risk factor for the development of oropharyngeal carcinoma [40]. There are reports indicating that the prognostic role of miR-21 is stronger in HNSCC patients with HPV-negative status tumors. Keeping in mind that miR-21 targets molecules that interfere with p53 pathways, this could be explained by the fact that in HPV-positive tumors, the p53 protein expression is not altered, being more sensitive to chemo-or/and radiotherapy [24,41]. Our study revealed that HPV status did not moderate the prognostic value of the expression levels of miR-21, indicating that this molecule may act as a stable prognostic marker for the survival of patients diagnosed with HNSCC. Therefore, more details regarding the characteristics of these patients should be provided in each study to facilitate data accessibility in order to obtain reliable results in this respect.

Conclusions
In conclusion, using fixed effects model, our meta-analysis was successfully underlined the potential of miR-21 and its stability as prognostic biomarker in head and neck squamous cell carcinoma (HNSCC) patients. Our results confirm the oncogenic nature of this transcript in this group of malignancies. Taking into consideration the diversity of these cancers, the oncogenic effect of miR-21 is undeniable, with a bioinformatics analysis confirming the interaction with genes with known involvement in carcinogenesis. The current meta-analysis was successful at confirming the prognostic biomarker potential of miR-21 in HNSCC patients, with highlights regarding the oncogenic role of this transcript observed in HNSCC cell lines. However, in order to transfer its utility in clinical trials with the ultimate purpose of developing personalized treatments for these patients, standardization of the used methodologies and of the clinical data that should be published is required, together with deciphering the molecular processes that lead to the abnormal expression of miR-21 and its downstream consequences in this group of cancers.

Conflicts of Interest:
The authors declare no conflict of interest.