DCBLD1 Overexpression Is Associated With a Poor Prognosis in Head and Neck Squamous Cell Carcinoma

Background DCBLD1 is highly expressed in several kinds of cancer and plays a potential prognostic factor. However, the prognostic value and immune infiltration in head and neck squamous cell carcinoma remain unclear and need further research. Materials and Methods DCBLD1 expression and clinical information were obtained from the Cancer Genome Atlas (TCGA) database. The mRNA level in cell lines (SCC25 and CAL27) and gingival fibroblasts were detected using quantitative PCR. Cox regression analysis was used to evaluate the prognostic values of DCBLD1 and clinical data in HNSCC. A nomogram was also established to predict the impact of DCBLD1 on prognosis based on Cox multivariate results. The methylation level of DCBLD1 in HNSC and its prognosis were analyzed in UALACN and MethSurv. Finally, the potential biological functions of DCBLD1 were investigated using gene set enrichment analysis (GSEA) and single-sample GSEA (ssGSEA). Results The mRNA and protein expression levels of DCBLD1 were highly expressed in HNSCC tissue and cell lines. The Cox analyses demonstrate that highly expressed DCBLD1 is an independent prognosis marker (p < 0.05). ROC curve analysis showed the performance of DCBLD1 (area under the ROC curve: 0.948, sensitivity: 93.2%, specificity: 84.7%). The methylation was increased in HNSCC patients compared with normal subjects (p < 0.05) and was associated with poor prognosis at sites cg27642470 and cg21104965. Additionally, DCBLD1 expression is poorly associated with immune cell infiltration and immunological checkpoints PD-L1 and TIM-3. Conclusion In head and neck squamous cell carcinoma, DCBLD1 is overexpressed, associated with poor patient prognosis. The detailed underlying mechanism merits further research.


INTRODUCTION
Head and neck cancer is a broad term that encompasses epithelial malignancies in the oral cavity, hypopharynx, nasopharynx, and larynx, mainly consisting of squamous cell carcinomas [head and neck squamous cell carcinoma (HNSCC)]. HNSCC is the sixth most frequent cancer worldwide, with an incidence of over 830,000 new cases annually (1). Although surgical and comprehensive treatment techniques over the last few decades have improved the quality of life and longevity, the prognosis of HNSCC patients remains poor due to local recurrence and distant metastasis.
The clinical significance of the mRNA expression has been previously described in several human cancer, including HNSCC (2). In addition, they play an important role in various physiological and pathological processes. For example, overexpression of SEC61G and HMMR has been associated with poor prognosis in HNSCC (3,4). DCBLD1, a transmembrane protein with extracellular CUB (an extracellular domain of approximately 110 residues that is found in functionally diverse, mostly developmentally regulated proteins), LCCL (including 94-97 amino acids, non-cytoplasmic but the function is undefined), and F5/8 type C domains 8, is a type I transmembrane protein (A-D) that binds to semaphorins (5,6) belonging to the DCBLD family, and is rarely studied and poorly characterized (5). DCBLD2 is another member of the DCBLD family and has been identified as a prognostic marker in HNSCC (7). It was reported that DCBLD2 was associated with lung cancer metastasis and involved cell growth and proliferation. DCBLD1 and DCBLD2 have similar amino acids, meaning that they may have related functions (5,8). The domain structure of DCBLD also closely resembles neuropilins harboring two CUB and discoidin domains (9). Additionally, both act as co-receptors for class 3 semaphorins and growth factors in axon guidance and angiogenesis (10). Neuropilin proteins are a class of transmembrane receptors that express on the surface of a large number of tumor cells and are widely involved in regulating tumor angiogenesis. Previous research showed that neuropilins could regulate multiple biological behaviors of tumor cells by acting as an oncogene, such as proliferation, migration, invasion, and apoptosis (9)(10)(11). In non-smoking women, a susceptible locus for lung cancer in the DCBLD1 gene was noted (12), suggesting the tumorigenic potential of DCBLD1 for the first time. Wang et al. (13) found that the knockout of DCBLD1 gene significantly affected the tumor formation ability of non-small cell lung cancer cells (A549) and also suppressed cellular proliferation (A549 and NCI-H1299). Cardin et al. (14) reported a potential association of DCBLD1 with the prognosis of HNSCC by studying the expression of rs6942067, a singlenucleotide polymorphism (SNP) located upstream of the DCBLD1 gene. However, there are no relevant studies about the role of DCBLD1 in the tumor microenvironment and prognosis in HNSCC.
In the present study, the RNA-sequencing data of DCBLD1 were downloaded from The Cancer Genome Atlas (TCGA) database. We determined the differential expression of DCBLD1 in pan-cancer and HNSCC from 546 patients. Culturing two HNSCC lines, SCC25 and CAL27, and human gingival fibroblasts (HGFs), the mRNA levels for DCBLD1 were quantified using real-time PCR among different groups. Subsequently, we detected the relationship between DCBLD1 expression with clinical characteristics and prognostication and established a predictive nomogram to investigate the functional role of DCBLD1 in HNSCC using functional enrichment analysis, protein interaction network, immune infiltration, mutation, and methylation.

Data Acquisition
Level 3 HTSeq-FPKM (fragments per kilobase of transcript per million mapped reads) data of HNSCC patients, including 44 normal and 502 tumor cases, were downloaded from the TCGA data portal (http://tcga-data.nci.nih.gov/tcga/), then were transformed into transcripts per million reads (TPM) and log2-transformed for subsequent analyses. Gene expression data were divided into high and low groups according to the median expression levels of DCBLD1. As TCGA was an open public database, obtaining relevant information did not require additional ethics approval.

Analysis of Differentially Expressed Genes Between High and Low DCBLD1 Expression Groups in HNSCC Patients
Expression profiles (HTSeq-Counts) were compared between the high DCBLD1 expression group and the low DCBLD1 expression group to identify differentially expressed genes (DEGs) using Wilcoxon rank-sum test (15) in the R languagerelated software, DESeq2 (version 1.26.0). Differences with a | log2 fold change| > 2 and adjusted p-value < 0.05 were considered threshold values for identifying DEGs.

Enrichment Analysis
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the top 200 DCBLD1-related genes were performed by the "ClusterProfiler" package (16) and visualized by the "ggplot2" package. In addition, the protein-protein interaction network of DCBLD1 co-expressed genes was visualized by STRING (http://string-db. org; version 11.5) with a minimum level of confidence >0.4 to analyze the functional interactions among proteins (17).

Cell Culture
Two HNSCC cell lines (CAL27 and SCC25) were purchased from the American Type Culture Collection (ATCC). Moreover, gingival fibroblast cells were cultured from normal tissues.

RNA Isolation and qPCR
SCC25, CAL27, and gingival fibroblast cells were harvested using a cell scraper after washing 3 times with phosphate-buffered saline (PBS). RNA was extracted through the following steps: lysing cells in Trizol reagent (Cat No. 15596-026, Life Technologies, Carlsbad, CA, USA), followed by extracting RNA in trichloromethane, and then precipitating it in isopropanol, and finally resuspending it in RNase-free water. The RNA concentration and purity levels were determined using a Nanodrop2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Total RNA (2.5 µg) was subjected to cDNA synthesis using a qScript cDNA SuperMix (Quanta Biosciences, Beverly, MA, USA) through the following consequent cycles: firstly at 25°C for 5 min, followed by 42°C for 30 min and finally at 85°C for 5 min. A real-time PCR was performed to determine the mRNA levels of DCBLD1 and GAPDH using SYBR Green Master MIX (ABI, Vernon, CA, USA). Real-time PCR results were calculated using the 2 −DDCq method.

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) (18) is used to access the concordance between our data and a priori-defined gene set. This study performed GSEA with the R package ClusterProfiler (3.14.3) (16) to elucidate the significant function and pathway differences between the high and low DCBLD1 expression groups. Each analysis procedure was repeated 5,000 times. A function or pathway term with adjusted p-value < 0.05, false discovery rate (FDR) < 0.25, and normalized enrichment score |NES| > 1 was considered to be statistically significant enrichment.

Immune Infiltration Analysis by Single-Sample GSEA
Immune infiltration analysis of HNSCC samples was performed by the single-sample GSEA (ssGSEA) (21) method using the GSVA package in R for 24 types of immune cells, namely, neutrophils, mast cells, eosinophils, macrophages, natural killer (NK) cells, CD56dim NK cells, CD56bright NK cells, dendritic cells (DCs), immature DCs (iDCs), activated DCs (aDCs), plasmacytoid DCs (pDCs), T cells, CD8 + T cells, T helper (Th) cells, Th1 cells, Th2 cells, Th17 cells, T follicular helper cells, regulatory T cells (Treg), central memory T cells (Tcm), effector memory T cells (Tem), gamma delta T cells (Tgd), cytotoxic cells, and B cells. The relative enrichment score of each immunocyte was calculated using the reported signature genes for the 24 types of immunocytes and the gene expression profile for each tumor sample. Furthermore, the relationship between DCBLD1 and immune cell markers was analyzed by TIMER (http://timer.comp-genomics.org/) online tools (22).

Statistical Analysis
All data analyses and plots were performed using R language (Version 3.6.3) and visualized by R package ggplot2. Wilcoxon rank-sum test and Wilcoxon signed-rank test were used to analyze the expression of DCBLD1 in unpaired and paired samples, respectively. In addition, ROC analysis and the frequently used method for binary assessment were conducted using the pROC package (1.17.0.1) to assess the diagnostic capability of DCBLD1 in head and neck cancer. The computed AUC value from 0.5 to 1 indicates the discriminative potential from 50% to 100%. The prognostic data were obtained from Cell (23), while the survival package (3.2-10) for Cox regression analyses and the Kaplan-Meier method were used to evaluate prognostic factors. In all tests, p-value < 0.05 was considered statistically significant.

DCBLD1 mRNA Is Overexpressed in HNSCC
The TCGA database was utilized to identify the DCBLD1 mRNA expression in different cancers, which was highly expressed in 12 of the 33 tumor types, as shown in Figure 1A. DCBLD1 mRNA expression was significantly higher in 502 head and neck cancer patients than those in 44 non-cancer subjects (p < 0.001, Figure 1B) and was also higher in 43 tumor tissue paired with pericarcinous tissue (p < 0.001, Figure 1C). The immunohistochemistry was performed to validate the expression of DCBLD1 in normal and HNSCC tissue at the protein level, which was downloaded from the HPA database ( Figure 1E). To further analyze the significance of DCBLD1 expression, qPCR results were shown in HNSCC cell lines (SCC25 and CAL27) and gingival fibroblasts. ( Figure 1F). A relatively high AUC value showed an excellent predictive power of DCBLD1 expression to identify tumors from normal tissue (AUC = 0.948, Figure 1D).

Identification of Differentially Expressed Genes Between the High and Low DCBLD1 Expression Groups
TCGA data were analyzed using the DSEeq2 package in R (|log FC| > 2, adjusted p-value < 0.05), and a total of 379 DEGs were found between the groups with high and low DCBLD1 expressions, including 349 downregulated genes in blue color and 30 upregulated genes in red color ( Figure 2A). In addition, the top 20 DEGs are shown as a heatmap ( Figure 2B).   Figures 3A-D).
The nomogram plot showed the prognosis of DCBLD1 and the relative clinical situation. Based on the multivariate Cox analysis, a nomogram was assigned to the clinical characteristic of a point, and the sum of the points was awarded as the total points. The probability of DCBLD1 survival at 1, 3, and 5 years can be determined using the absolute point axis down to the outcome axis ( Figure 3E). The calibration curve evaluated the nomogram's performance of DCBLD1. The C-index of nomogram was 0.720 with 1,000 bootstrap replicates (95% CI: 0.693-0.747). The bias-corrected line in the calibration plot was close to the 45°line, which is the ideal curve ( Figure 3F).

Function and Pathway Enrichment Analysis
The top ten related proteins to DCBLD1 were established by STRING (Figures 4A, B). A total of 197 correlated genes with DCBLD1 in the TCGA-HNSCC database were evaluated (|cor Pearson| > 0.5, p-value < 0.05), including 193 positively correlated genes and 4 negatively correlated genes. The PPI network of those 197 co-expressed genes was established in STRING and Cytoscape ( Figure 4C). GO and KEGG enrichment analyses investigated the functional significance of DCBLD1 co-expressed genes in HNSCC. The top 20 enriched items, including cellular components, molecular function, biological process, and KEGG, are shown in Figure 4D.

Gene Set Enrichment Analysis
Based on significant differences (adjusted p-value < 0.05, FDR < 0.25), GSEA was used to identify signaling pathways associated with HNSCC between the high and low DCBLD1 expression groups, including 102 positive regulation pathways and 74 negative regulation pathways. The most significantly enriched pathways were the integrins 3 pathway, met promotes cell motility, syndecan interactions, met activates PTK2 signaling, matrix metalloproteinases, MIR5093P alteration of YAP1ECM axis, UPA UPAR pathway, and laminin interactions ( Figure 5).

Bioinformatical Analysis
The correlation between DCBLD1 expression and its mutation in pan-cancer was analyzed using cBioPortal, as revealed in Figure 6A.
In 523 HNSCC cases, the genetic alteration was found in eight patients, and the mutation rate was 1.53%. Moreover, the results revealed that genetic alteration in HNSCC was associated with superior OS of HNSCC, implying that the genetic mutation of DCBLD1 could also affect HNSCC patients' prognosis. The methylation level of DCBLD1 in HNSCC was examined by UALCAN based on TCGA. The methylation was higher in the high DCBLD1 expression group ( Figure 6B). MethSurv analysis showed that patients with high DCBLD1 methylation had a worse OS than patients with low DCBLD1 methylation ( Figure 6C). The 2 CpG sites are shown in Figures 6D, E.

The Correlation Between DCBLD1 and Immune Infiltration in HNSCC
The association between DCBLD1 expression and immune cell infiltration was visualized using ssGSEA and the TIMER    C). The correlation of DCBLD1 and immune cell markers analyzed by TIMER is shown in Table 3, which demonstrates that high DCBLD1 expression was associated with most immune cell markers but also with a poor correlation. Previous studies showed that higher B-cell and CD8 + T  infiltration were associated with poor prognosis. The SCNA module compares tumor infiltration relationships among tumors with different somatic copy number variations for DCBLD1 ( Figure 7B). Solid tumors, especially melanoma and renal cancer, have significantly benefited from immunotherapy, especially from immune checkpoint inhibitors (ICIs), which also offer potential predictive and therapeutic strategies for R/M HNSCC. In 2014, FDA approved using nivolumab and pembrolizumab, targeting PD-1/PD-L1 interaction. Both demonstrated better survival and response rates in clinical trials CheckMate-141 and KEYNOTE-012. TIM-3 and PD-L1 play an important role in tumor immune escape, which were also recognized as immune checkpoints in HNSCC. The expression of DCBLD1 was positively related to those immune checkpoints based on the HNSCC-TCGA database ( Figures 7D, E).

DISCUSSION
This study analyzed DCBLD1 expression in pan-cancer using the TCGA data and found upregulated expression in BLCA, COAD, DLBA, ESCA, GBM, HNSC, KIRC, LAML, LGG, LUSC, OV, PAAD, READ, STAD, and THYM as compared to healthy individuals. Subsequently, HNSCC was focused, and a higher expression of DCBLD1 was demonstrated in cancer tissues than in normal tissues (p < 0.001). Combining HPA data and qPCR results, it was found that DCBLD1 proteins were moderately expressed in OSCC tissues but were absent in normal oral tissues. SCC25 and CAL27 cancer cell lines are derived from human OSCC commonly used in HNSCC studies (24), while gingival fibroblasts originated from the gingival tissues resected in normal people after the removal of impacted teeth. The qPCR results indicated that the expression of DCBLD1 mRNA in OSCC cells  was higher than that in normal gingival fibroblasts, which further suggested that high DCBLD1 expression might associate with the development or progression of HNSCC. Previous studies about DCBLD1 were scarce, so the carcinogenic role of DCBLD1 remains unclear. Concurrently, Wang et al. (13) found that DCBLD1 was associated with tumor formation by knocking down this gene. Furthermore, the AUC score reached a value up to 0.948, demonstrating good diagnostic power of DCBLD1 expression toward HNSCC (25). Previous studies also showed the associations of the DCBLD1 expression with lung cancer and breast cancer. All these findings suggest the oncogenic potential of DCBLD1. The present study also revealed the significant prognostic value of DCBLD1 in the survival outcome of HNSCC. According to Kaplan-Meier curves, high DCBLD1 expression predicted a shorter OS time in patients with HNSCC. Consistent findings were also demonstrated in subgroup analysis. This also implied the potential of DCBLD1 as an oncogene. Additionally, DCBLD1 was identified as an independent prognostic factor for the OS and DSS of HNSCC using univariate and multivariate Cox regression analyses. Collectively, DCBLD1 might be a potential therapeutic target for treating HNSCC. The nomogram is a reliable, effective tool for cancer prognosis (26). Here, a nomogram combining DCBLD1 expression, primary therapy outcome, lymphovascular invasion, and radiation therapy was established and identified to be effective in the prognosis of HNSCC according to the calibration curve. We noted that high DCBLD1 expression accounted for most of the contribution to the nomogram in the prognosis of HNSCC, suggesting that DCBLD1 is a potent prognostic factor. In addition, the present study also explored the role of mutation and promoter methylation of the DCBLD1 gene in HNSCC prognosis. Mutations in the DCBLD1 gene were associated with good survival outcomes. It was also noted that DCBLD1 gene promoter methylation was much more prevalent in HNSCC tissue than in normal tissue. The current study indicated that methylation at cg27642470 and cg21104965 predicted a poor prognosis.
The tumor microenvironment is significantly associated with tumor growth, development, and prognosis (27). Inflammatory cells, essential components of the tumor microenvironment (27,28), can advance the degradation of the extracellular matrix (ECM), potentiate angiogenesis and tissue remodeling, and increase cellular migration. Koneva et al. (29) analyzed lymphocyte subtypes in HNSCC and reported many tumor-infiltrating lymphocytes. Another retrospective study about HNSCC showed that a high infiltration abundance of CD4 + and CD8 + T cells indicated a higher risk of death (30), demonstrating that upregulated infiltration of the two immune cell types was prognostic for the recurrence-free survival in HNSCC patients. Here, we found that the expression of DCBLD1 was positively associated with the infiltration abundance of eosinophils, iDC, macrophages, neutrophils, NK cells, Tcm, Tgd, Th1 cells, and Th2 cells in the tumor microenvironment and negatively associated with the infiltration abundance of B cells, CD8 T cells, cytotoxic cells, NK CD56 bright cells, pDCs, and Th17 cells. However, the correlation was poor in all correlations. This indicated that DCBLD1 is not associated with immune cell infiltration in head and neck cancer.  The copy number of the DCBLD1 gene was also observed to play a role in immune infiltration in HNSC. According to the data retrieved from the GEPIA database, DCBLD1 gene expression was positively correlated not only with the expression of immune checkpoints PD-L1 and TIM-3 but also with a poor association.
Immunotherapy is a step forward from traditional chemoradiotherapy based on interactions between tumor cells and the immune microenvironment (31,32).
In the tumor microenvironment, infiltrating immune cells experience functional impairment by expressing multiple inhibitory signals on the cell surface, such as PD-1/PD-L1, CTLA-4, TIGIT, and TIM-3, leading to tumor immunosuppression. These inhibitory signals are also known as immune checkpoints (33). PD-1/ PD-L1 is among the earliest used in treating advanced melanoma (34). The tumor cell-intrinsic PD-1 receptor is an important modulator of immune responses and primarily suppresses them, negatively regulating T cells and immune response. PD-1:PD-L1 binding inhibits T-cell activation, decreasing proliferation and effector cytokine production. The FDA has approved pembrolizumab and nivolumab for clinical treatment. In the clinical tests, pembrolizumab combined with chemoradiotherapy compared to chemoradiotherapy alone improved tumor response and progression-free survival in patients with recurrent/metastatic (R/M) HNSCC. Close to 20% of them survive more than 4 years. With nivolumab, patients had a longer OS than those with standard treatment (7.5 months versus 5.1 months). PD-L1 overexpression was demonstrated in 18%-96% of OSCC patients and associated with lymph node metastasis (35). Zhang et al. believed that PD-L1 overexpression predicted a higher risk of lymph node metastasis (36). Other reports also revealed that PD-L1 overexpression was correlated with the cytokine-induced immune silence of lymphocytes or regional recurrence. These studies suggest that PD-L1 is intimately associated with metastasis to peri-tumor lymph nodes (37). TIM-3 expresses on T cells, NK cells, and some antigen-presenting cells. It is believed to be positively associated with lymph node metastasis and tumor recurrence. Previous research found that high expression of TIM-3 led to effector T-cell depletion, which might cause an ineffective anti-tumor immune response and tumor clearance, resulting in metastasis and recurrence (38). The current study identified a poor association between DCBLD1 expression and PD-L1 and TIM-3.
Following the GSEA study, the integrin pathway was the most significant (39). It has been established that the integrin pathway will be triggered upon binding with ECM, accompanied by activation of downstream actin reorganization and MAPK signaling cascade (40). In HNSCC, the integrin pathway is involved in disease progression, tumor stem cells, and radioresistance (41,42). We also noted that GO annotations related to DCBLD1 gene included focal adhesion, cell-substrate junction, ECM organization, cell adhesion molecule binding, and other biological processes involved in ECM and cell adhesion. This suggests the important role of DCBLD1 in tumor progression and migration. However, the role of DCBLD1 in developing SCC of the oral cavity and other sites of the head and neck region has not yet been fully elucidated. The potential of DCBLD1 as an oncogene in HNSCC has not yet been demonstrated and is a matter of future studies. More in vivo and in vitro experiments are needed to explore this issue further. We plan to conduct cytological experiments in our future experiments.

CONCLUSION
For the first time, this study reported the effects of the DCBLD1 gene on the tumor microenvironment and prognosis in HNSC.
Our study suggests that DCBLD1 is a potential prognostic factor for survival outcomes but is not associated with immune cell infiltration. DCBLD1 may also play a role in HNSC progression and metastasis via the integrin pathway and cell adhesion. DCBLD1 is expected to be a novel diagnostic and prognostic factor for HNSCC.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: TCGA.

AUTHOR CONTRIBUTIONS
L-lF: conceived the study, supervised the experiments, and drafted the manuscript. MY: conceived the study, supervised the experiments, and drafted the manuscript. M-XM: data evaluation and manuscript preparation. YL and MS: data evaluation and manuscript preparation. MG, RF, TV, and RS: analyzed the data and revised the manuscript. H-cF: conceived the study, designed the data evaluation, and prepared the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported in part by grants from the Guiyang Municipal Health Bureau Fund for Science and Technology projects. The funder is in charge of the publication fees.