Extraction of qualified microarray data and meta-analysis
According to the inclusion criteria, we screened 18 eligible GSE Microarray Data from the GEO database (Figure 2), as follows: GSE69002, GSE70289, GSE34496, GSE45238, GSE58911, GSE37472, GSE11163, GSE28100, GSE73171, GSE82064, GSE98463, GSE103931, GSE107591, GSE113956, GSE116994, GSE124566, GSE168227, GSE207030. There were total of 865 HNSCC samples and 292 normal samples from the GEO and TCGA datasets, and the relevant information was shown in Table 1. All the above data was used to evaluate miR-15b-5p expression in HNSCC by meta-analysis. As shown in Figure 3, I square(I2) was 82.3%, P< 0.05, indicating that there was high heterogeneity in the study data, so we chose the random effects model. The meta-analysis forest plot showed that the combined SMD was 0.49 (95% CI: 0.10-0.88), and the horizontal line of SMD 95 % CI fell to the right of the null line, indicating that miR-15b-5p was highly expressed in HNSCC. In the GSE34496, GSE45238, GSE37472, GSE28100, GSE103931, GSE124566 and TCGA datasets, miR-15b-5p expression was significantly upregulated in HNSCC than in normal group (P<0.05) (Figure 4).
Potential diagnostic role of miR-15b-5p in HNSCC
As mentioned before, with high heterogeneity among each dataset, random-effects model was used to comprehensively evaluate the potential diagnostic role of miR-15b-5p in HNSCC by sensitivity, specificity, PLR, NLR, and DOR. As shown in Figure 5A-S, the ROC curve showed the diagnostic accuracy of miR-15b-5p in HNSCC for each dataset. In addition, Diagnostic test demonstrated that the pooled sensitivity, specificity, PLR, NLR, and DOR values were as follows: 0.70 (95% CI:0.67-0.73), 0.72 (95% CI: 0.67-0.77), 2.15 (95% CI: 1.57-2.94), 0.48 (95% CI: 0.35-0.67), and 7.33 (95% CI: 4.95-10.86), respectively (Figure 6A-E). Based on the SROC curve analysis, the area under curve (AUC) was 0.81 (95% CI: 0.77-0.84) (Figure 6F), indicating that miR-15b-5p played a potential diagnostic role in identifying HNSCC tissues.
MiR-15b-5p was overexpressed in HNSCC cells
RT-qPCR data indicated that miR-15b-5p was overexpressed in HNSCC cells compared to normal cells (P < 0.05) (Figure 7A-B), which was consistent with the results of meta-analysis based on GEO and TCGA databases. This result indicated that miR-15b-5p expression was linked to HNSCC and therefore, we developed additional experiments to examine its mechanism of action in HNSCC.
Clinical features of miR-15b-5p in patients with HNSCC
We downloaded 527 HNSCC clinical samples from the TCGA database and performed statistical analyses. The result demonstrated that miR-15b-5p expression was higher in HNSCC than in normal tissues (P < 0.05) and statistically different with respect to gender, age, N stage and T stage (P < 0.05). Specifically, miR-15b-5p expression was higher in male than female patients, in those aged > 60 years than in those aged < 60 years, in N1- 4 stage than in N0 stage, and in T3- 4 stage than in T1-2 stage. However, miR-15b-5p expression level was not significantly correlated with tumor pathological grade, clinical stage, or M stage (P > 0.05) (Table 2). These findings indicated miR-15b-5p expression was closely associated with the clinical characteristics of patients and may be a potential marker to evaluate HNSCC progression.
Identifying MiR-15b-5p key genes in HNSCC
We searched 3336 downstream genes of miR-15b-5p in the miRWalk database, identified 848 differentially expressed genes in HNSCC using the GEPIA2 database, and finally obtained 54 intersecting genes with VENNY graph analysis that are likely miR-15b-5p key genes in HNSCC (Figure 8A; Table 3). we established a PPI network map of key genes, and hid proteins that did not interact. The network shown that the close connections among miRNA-15b-5p key targets (Figure 8B).
Identifying miR-15b-5p target genes in HNSCC
The 54 intersection genes were performed GO and KEGG enrichment analysis to farther research the role and signaling pathways of miR-15b-5p in HNSCC. GO analysis results demonstrated that key genes were highly enriched in extracellular matrix organization, extracellular structural organization, regulation of cellular component movement, and flavone metabolic process in BP. In term of CC, the most enriched terms were extracellular matrix component, proteinaceous extracellular matrix, basement membrane, anchoring junction, and apical part of cell. Those in the MF were closely related to filamin binding, extracellular matrix structural constituent, and protein kinase C binding (Figure 9A-D; Table S1). The KEGG enrichment analysis revealed that miR-15b-5p affected HNSCC through multiple pathways, including retinol metabolism, drug metabolism-cytochrome P450 and metabolism of xenobiotics by cytochrome P450 (Figure 9E; Table 4). The retinol metabolism pathway target genes included RDH12, UGT1A10, and UGT1A7, which ranked second in the KEGG pathway with a P < 0.05. Based on the enrichment analyses, miR-15b-5p possibly played an importance part in HNSCC by regulating retinol metabolism.
Functional analysis of miR-15b-5p target genes
Compared with normal tissues, two target genes (RDH12 and UGT1A7) associated with miR-15b-5p and retinol metabolism signaling pathway were downregulated in HNSCC (P < 0.05), whereas UGT1A10 expression was not statistically different (P > 0.05) (Figure 10 A-C). Furthermore, the methylation levels of target genes UGT1A10 and UGT1A7 were lower in HNSCC than in normal tissues (P < 0.05) (Figure 11A-B). Spearman correlation analyses suggested that RDH12 was inversely proportional to miR-15b-5p expression in HNSCC, whereas UGT1A7 was directly proportional to miR-15b-5p expression (P < 0.05). There was no statistically correlation between UGT1A10 and miR-15b-5p levels (P > 0.05) (Figure 12A-C). Together, the above results demonstrated that miR-15b-5p was involved in regulating the retinol metabolism pathway by targeting RDH12 and UGT1A7, thereby affecting the progression of HNSCC.
Evaluate the prognosis of target genes in HNSCC
Kaplan-Meier curve analysis was performed on 518 patients with HNSCC to evaluate the effect of target gene expression levels on HNSCC prognosis. Patients with high UGT1A7 expression had noticeably longer DFS than those with low expression (P < 0.05), and the effect of UGT1A7 expression on OS had no statistical significance (P > 0.05) (Figure 13A-B). However, there were no statistical differences between the effects of UGT1A10 and RDH12 on both OS and DFS in patients with HNSCC (P > 0.05) (Figure 13C-F). These findings suggested that high UGT1A7 expression is a favorable factor for good DFS in patients with HNSCC.