The expression profile and clinical application value of hsa_circ_0016148 in head and neck squamous cell carcinoma

Abstract Background Dysregulated circular RNAs (circRNAs) are involved in human cancers and may be used as biomarkers with the potential of clinical application. However, little is known regarding the mechanism of circRNAs and their clinical application value in head and neck squamous cell carcinoma (HNSCC). Methods In the current study, we established the profile of circRNAs in HNSCC using microarray and then measured the expression of hsa_circ_0016148 in 137 paired HNSCC tissues by qRT‐PCR technique, analyzed the relationship between hsa_circ_0016148 and clinicopathological data, and investigated its diagnostic and prognostic value. The hsa_circ_0016148‐miRNA‐mRNA interaction network was predicted and constructed by Cytoscape. Results Our study showed a circRNA expression profile and confirmed downregulated hsa_circ_0016148 in HNSCC tissues (p = 2.64E‐35). The hsa_circ_0016148 expression is remarkably correlated with lymph node metastasis (p = 0.001) and clinical stage (p = 0.029). Then, the area under the receiver characteristic curve (AUC) was 0.912 with 92% of sensitivity and 86.9% specificity, respectively. Besides, our study demonstrated that lower‐expressed hsa_circ_0016148 could independently predict poorer overall survival of HNSCC patients (hazard ratio [HR] = 0.456; 95% confidence interval [CI] = 0.265–0.784; p = 0.005). The hsa_circ_0016148‐miRNA‐mRNA interaction network was constructed, which included a total of nine targeted miRNAs. Conclusion Taken together, our results revealed that hsa_circ_0016148 might play a critical role in HNSCC tumorigenesis and may serve as an indicator with the potential of diagnosis and prognosis for HNSCC.

half cases were dead annually. 4 The high-risk factors for HNSCC include upper aerodigestive tract mucosa exposure to carcinogens such as tobacco, alcohol, and human papillomavirus (HPV), with a multiplicative effect. 5,6 While the use of tobacco is decreasing around the world, there is a recent upsurge in the incidence of HPV-related HNSCC. 7 Despite improvements in radiation technology and techniques, better surgical options for organ preservation, and multidisciplinary treatment modalities, HNSCC patients continue to have poor survival rates due to local recurrence and distant metastasis, particularly in advanced patients. 8 Recently, immune-checkpoint inhibitors nivolumab and pembrolizumab have been approved for treatment of platinum chemotherapyfailed HNSCC patients with recurrent or metastatic. This treatment has provided significantly better outcomes, but only a small proportion of patients benefit from it. The 5-year survival rate is still only 50% and has slightly improved in the last two decades. 2 The current plight is partly due to the lack of accurate and reliable tools for early diagnosis and prognosis prediction of HNSCC. 9 Therefore, identifying an effective and specific biomarker is urgent for severity assessment and improvement of personalized medical treatment for patients with HNSCC.
In the beginning, circular RNAs (circRNAs) are recognized as a class of noncoding RNAs that are covalently bonded to form a special closed-loop structure without 5′ caps and 3′ tails. 10 However, recent studies have shown that some circRNAs could also be translated, which might be driven by a single N6-methyladenosine residue. 11,12 CircRNAs are considered as abundant, stable, and conserved across evolutions, and many play important biological functions by acting as microRNA (miRNA) sponges that regulate protein function or by being translated themselves. 13 In addition, circRNAs can modulate several important physiological processes through interacting with RNA-binding protein, such as apoptosis. 14 Evidence is increasing that circRNAs are involved in several diseases including cancers. 13,15 One recent study by Zhang and his colleagues showed that circular RNA circ_0001287 inhibited the proliferation, metastasis, and radiosensitivity of non-small cell lung cancer cells by sponging miR-21. 16 Hao et al. reported that circDCUN1D4 suppressed tumor metastasis and glycolysis in lung adenocarcinoma by stabilizing TXNIP expression.
Additionally, due to its characteristics of abundance, diversity, stability, and disease-specific exhibition, circRNA was considered as ideal biomarker for auxiliary diagnosis and prognosis estimation of disease. 17,18 Nevertheless, to the best of our knowledge, the mechanism of circRNAs in HNSCC is still fairly elusive.
In the current study, a circRNA expression profile was constructed from eighteen paired HNSCC samples and corresponding normal tissues using a microarray platform. circRNA-hsa_circ_0016148 was found downregulated in the microarray and selected for further study. The gene of hsa_circ_0016148 is located at chromosome chr1:204410598-204411761. Its spliced sequence length is 201 nt; the associated gene symbol is phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 beta (PIK3C2B). We believe that this is the first time hsa_circ_0016148 expression in HNSCC tissues and normal controls has been investigated or that the relationship with clinicopathological factors of patients with HNSCC has been evaluated. Moreover, we explored its potential diagnostic and prognostic value and regulatory mechanisms for HNSCC. Potential biological functions were further predicted and annotated by bioinformatics analysis. Lihuili Hospital is obtained before the study begins. Informed consents were signed by all included patients.

| Total RNA extraction and quantitative realtime PCR assay
Following the manufacturer's instructions, we extracted total RNA from tissues using TRIzol reagent (Invitrogen). Then, cDNA was reverse-transcribed using the GoScript Reverse Transcription (RT) System (Promega). The primers for the amplified sequences of hsa_ circ_0016148 were: forward, 5′-TCTTCAAGGGAATCCTCCGC-3′, and reverse, 5′-CCTTAACCAGCAGACTGGGG-3′. GAPDH was simultaneously amplified as internal control. 20 RT-PCR was performed using the RocheLightCycler 480II System (Roche). The PCR was con- and GAPDH were recorded. The relative quantification level of hsa_ circ_0016148 was calculated by the ΔCq method using GAPDH as a reference gene. All data were presented as means ± SD (standard deviation) through at least three independent experiments.

| Diagnostic value analysis
We generated a receiver operator curve (ROC) and calculated the area under curve (AUC) to assess the diagnostic efficacy of hsa_ circ_0016148. 21 The maximum Youden index was established to identify the optimal cut-off value of ROC.

| Survival analysis
According to the optimal cut-off value of ROC, HNSCC patients were divided into high and low hsa_circ_0016148 groups. Then, Kaplan-Meier analysis was performed between these two groups.
Univariate Cox regression analysis was used to identify potential prognostic factors, which were then involved in the multivariate Cox regression analysis to evaluate the independent risk factors associated with overall survival of HNSCC.

| Finding the potential target miRNAs of hsa_circ_0016148 and ceRNA network analysis with bioinformatics
First, we used the online circinteractome database (https://circi ntera ctome.irp.nia.nih.gov/index.html) to predict miRNA targets of hsa_circ_0016148. Then, the target genes of these miRNAs were predicted from the TARBASE v.8 and miRTARBASE databases.
Cytoscape software was used to construct the ternary circRNA-miRNA-mRNA interaction network. Finally, GO enrichment analysis was conducted for all target genes using org. Hs.eg.db and the clus-terProfiler package in R 3.6.0 software.

| Statistical analyses
SPSS version 18.0 (SPSS Inc.) was used for all statistical analysis. The comparison of hsa_circ_0016148 expression was made using paired sample Student's t test between HNSCC samples and controls. The association between hsa_circ_0016148 and clinical parameters was explored by independent sample Student's t test. We considered estimates to be statistically significant if the p value was less than 0.05 with a two-tailed test.

| Comparison of expression of hsa_ circ_0016148 in HNSCC
In this study, we established a circRNA expression profile from eighteen paired HNSCC specimens and adjacent normal tissues using a microarray platform. The volcano plot is shown in Figure 1.
Hsa_circ_0016148 was found downregulated in the microarray and selected for further study. We confirmed hsa_circ_0016148 expression levels in 137 paired HNSCC tissues and corresponding normal tissues by quantitative real-time polymerase chain reaction. The sequence of qRT-PCR product was identical to that from circBase (http://circb ase.org/). We found that hsa_circ_0016148 expression in HNSCC tissues was remarkably decreased than that in adjacent normal tissues (p = 2.64E-35, Figure 2A). The expression of hsa_ circ_0016148 was downregulated in 90.5% (124/137) of HNSCC tissues compared with the corresponding non-tumorous tissues ( Figure 2B).

| Association between hsa_circ_0016148 expression and clinicopathological parameters in HNSCC patients
Subsequently, we analyzed the association between hsa_ circ_0016148 expression level and clinicopathological characteristics, including age, gender, smoking behavior, alcohol behavior, tumor site, differentiation grade, tumor stage, lymph node metastasis, and clinical stage. As shown in Table 1, low hsa_circ_0016148 expression in HNSCC tissues was significantly associated with lymphatic metastasis (p = 0.001) and clinical stage (p = 0.029). The rest of clinicopathological characteristics showed slight association between hsa_circ_0016148 expression levels in HNSCC patients.

| The receiver operating characteristic (ROC) curve of hsa_circ_0016148 in HNSCC
The potential diagnostic value of hsa_circ_0016148 was identified by ROC curve in HNSCC patients. The value of Youden index was calculated by the following formula: sensitivity + specificity-1, and the maximum Youden index was defined as a cut-off point. In our study, the area under the ROC curve (AUC) was 0.912 at a cut-off F I G U R E 1 Volcano plot of circRNA microarray in head and neck squamous cell carcinoma (HNSCC). Log FC >3, p < 0.05 value of 11.715 with 92% sensitivity and 86.9% specificity, respectively ( Figure 3).

| DISCUSS ION
Despite tremendous research efforts and huge improvements in a combination of surgery, chemotherapy, radiotherapy, immunotherapy, and molecular-targeted agents for HNSCC, most patients are diagnosed with a locally advanced stage and have poor outcomes. 22 Identifying reliable biomarkers for early diagnosis is a key step to improve prognosis and quality of life for HNSCC patients. 23 Therefore, there is an urgent need to develop and verify an effective and specific biomarker to improve the prevention, diagnosis, and prognosis prediction for HNSCC. Recent functional studies have revealed many tissue-specific and cell-specific circRNAs and functionally characterized them as efficient biomarkers in human diseases, especially cancers. 24 For example, hsa_circ_0065149 has been demonstrated to be an indicator for early gastric cancer screening and prognosis prediction. 25 However, the mechanism and the clinical application value of circRNAs in HNSCC are still largely unknown.
Utilization of high-throughput sequencing and bioinformatic analysis can be efficient strategies to identify aberrantly expressed circRNAs in cancers. In the current study, we focused our atten- Patients in whom the primary neoplastic site is the oral cavity, upper glottis, or nasopharynx often present with cervical lymph node metastasis as their first presenting sign. 28 Additionally, conventional tumor-related blood biomarkers, such as carcinoembryonic antigen (CEA) and carbohydrate antigens 19-9 (CA19-9), do not have satisfactory sensitivity for HNSCC detection, especially in the early stages. 29 Therefore, due to non-specific symptoms in the early stages and lack of effective diagnostic biomarkers for HNSCC, the majority of HNSCC patients are diagnosed at an advanced stage (III and IV), at which point survival rates are low relative to those for early stages (I and II). 30 Due to the stability, tissue specificity, and

F I G U R E 3 Receiver operating characteristic (ROC) curve.
The cut-off point was defined as the maximum Youden index, highlighted by the arrow in the figure

F I G U R E 4 Association between hsa_circ_0016148 expression
and overall survival in head and neck squamous cell carcinoma (HNSCC). Low hsa_circ_0016148 expression is associated with poor overall survival (OS) in HNSCC patients the spatiotemporal specificity of circRNA, it could be a novel biomarker in early diagnosis of cancer patients. 31,32 In this study, we also calculated ROC curves to assess the value of hsa_circ_0016148 in HNSCC diagnosis. The higher AUC value indicates the better diagnostic capability. 33  These findings give us a hint that hsa_circ_0016148 may work through hsa_circ_0016148-miRNA-mRNA to regulate numerous pathophysiological processes associated with HNSCC tumorigenesis and progression, a conjecture which should be validated by rigorous in vitro and in vivo experiments in future.
In conclusion, our research revealed that hsa_circ_0016148 is remarkably downregulated in HNSCC and may play an important role in HNSCC tumorigenesis via the hsa_circ_0016148-miRNA-mRNA axis. We also discovered that hsa_circ_0016148 is an effective biomarker for the diagnosis of HNSCC with a high degree of accuracy, specificity, and sensitivity. Additionally, our findings supported hsa_circ_0005986 as an independent unfavorable factor for overall survival in HNSCC patients.

CO N FLI C T S O F I NTE R E S T
The authors declare that there is no conflict of interest regarding the publication of this article.

AUTH O R CO NTR I B UTI O N S
Liuqian Wang and Dong Ye wrote the manuscript. Zhisen Shen reviewed and edited the manuscript. Both authors have read and approved the final manuscript and take responsibility for its integrity.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used to support the findings of this study are available from the corresponding author upon request. F I G U R E 6 The top 15 biological processes and pathways with the most significant p-values. The x-axis represents the number of target genes involved in the pathway