Development and validation of glycosyltransferase related-gene for the diagnosis and prognosis of head and neck squamous cell carcinoma

Public datasets collection

RNA-seq data and relevant clinical characteristics of patients with HNSCC from The Cancer Genome Atlas (TCGA) database were downloaded as the training cohort. Simultaneously, we used data from GSE41613 as the test cohort for the prognostic model and the matrix files of GSE127165 as validation data for the diagnostic model. The GSE41613 and GSE127165 datasets contained clinical information and gene expression profiles of patients with HNSCC obtained from the Gene Expression Omnibus (GEO) database. We identified 169 glycosyltransferase-related genes in a study by Mohamed Abd-El-Halim et al. [19] and downloaded them for further analysis (Supplementary Table 1). All the data used in our study are publicly available.

Identification and validation of prognostic gene signature

In the TCGA cohort, the R package “limma” was used to recognize the differentially expressed genes (DEGs) by comparing the different gene expression levels of tumor and normal cells. The thresholds were set as |log₂FC| > 0.5 along with a false discovery rate < 0.05. Univariate Cox analysis was performed on these DEGs, and valuable genes (P<0.05) were further analyzed using multivariate analysis. The risk score was calculated using the following formula: Risk score = expression level of gene-1 × coefficient of gene-1 + expression level of gene-2 × coefficient of gene-2 +. . .+ expression level of gene-n × coefficient of gene-n [20]. Glycosyltransferase-related genes screened using the aforementioned methods were used to build a prognostic signature of patients with HNSCC. The expression levels of each gene were multiplied by the corresponding Cox coefficient and added to calculate the risk score for all clinical cases. HNSCC samples were then segmented into high- and low-risk subgroups using the median risk score. The prognostic power of the model was evaluated using Kaplan-Meier (K-M) survival curves with time-dependent receiver operating characteristic (ROC) curves. The same process was completed for the GSE41613 dataset to validate the signature stability. Finally, to estimate the independent prognostic value of the glycosyltransferase-related gene signature, univariate and multivariate Cox regression analyses were performed.

Functional enrichment analysis and mutation gene analysis

To better comprehend the biological functions of the target genes, the “clusterProfiler” package was used to perform gene set enrichment analysis (GSEA) of the two subgroups based on the Kyoto Encyclopedia of Genes and Genomes. In addition, somatic mutation information was collected from the cBioPortal database and gene mutation analysis was performed on both subgroups using the “Maftools” package.

Evaluation of immune cell infiltration and immunotherapy

To explore the correlation between the glycosyltransferase-related gene model and immune cell levels, the CIBERSORT (https://cibersort.stanford.edu/) algorithm was used to handle the expression data of each case to assess the proportion of different immune cells in the two risk groups of the TCGA cohort [21]. In addition, we analyzed the association of the prognostic risk score with immune checkpoint (IC) expression levels and immunotherapy scores.

Construction and evaluation of diagnostic model

In the TCGA cohort, two machine learning algorithms, the least absolute shrinkage and selection operator (LASSO) and support vector machine (SVM), were used to select genes with diagnostic value. SVM is a machine-learning technique widely used for tumor gene expression profile analysis and tumor marker detection in different types of cancers [22]. LASSO regression can reduce the number of estimated parameters while maintaining a high prediction accuracy, thereby reducing data overfitting and making the model easy to visualize and interpret [23]. The intersection genes of the two algorithms were used to construct predictive diagnostic models. ROC curve was used to examine the effectiveness of the model. As with the prognostic model, we used matrix data from GEO for model validation.

Cell culture

FaDu cells were obtained from the American Type Culture Collection (USA), and TU177 cells were procured from Otwo Biotech, Inc. (Shenzhen, China). These cells were cultured in a 37° C incubator with a 5% CO₂ atmosphere using DMEM or RPMI-1640 medium supplemented with 1% penicillin-streptomycin and 10% fetal bovine serum.

Lentiviral infection

TU177 cells were stably transduced with lentiviral vectors (GV248, GenePharma, Shanghai, China) carrying specific short hairpin RNAs (shRNAs) targeting B4GALT3. The shRNA sequences used were as follows: shB4GALT3-1, 5’-GGGATGAACTCACTGACATAC-3’, shB4GALT3-2, 5’-GGACGCAAGATGGGATGAACT-3’, and control scrambled shRNA (shSc), 5’-TTCTCCGAACGTGTCACGT-3’. Cell lines overexpressing B4GALT3 were established using the lentiviral vector GV492 containing the complete B4GALT3 sequence (GenePharma).

Western blot analysis

Immunoblotting was performed according to standard procedures [24]. Briefly, cell lysates were prepared using NuPAGE 4-12% Bis-Tris gels and subsequently transferred onto PVDF membranes (Millipore, USA). Before the application of primary antibodies, the membranes were blocked with 5% nonfat milk solution to prevent nonspecific binding of antibodies. Following overnight incubation at 4° C with primary antibodies and subsequent one-hour incubation with secondary antibodies, chemiluminescence was used to visualize the results. The commercial antibody anti-B4GALT3 (11041-1-AP, Proteintech, China) was used in these experiments.

In vitro functional assays

In vitro functional assessments, including CCK-8 assay (TargetMol; Shanghai, China), colony formation, wound healing, and Transwell assays, were employed to assess the functional roles of B4GALT3.

Statistical analysis

R software and various packages were used for the statistical analyses. Univariate and multivariate Cox regression analyses were used to screen the prognostic genes and independent prognostic parameters. Survival analysis was performed using the K-M analysis with the log-rank test. Spearman’s correlation analysis was performed to analyze the relationships. The Wilcoxon test was used for comparison between subgroups.

Data availability statement

Publicly available datasets were analyzed in this study. The raw data of this study are derived from the TCGA database (https://portal.gdc.cancer.gov/) and the GEO data portal (https://www.ncbi.nlm.nih.gov/geo/), which are publicly available databases.

Glycosyltransferase-related DEGs were identified

In the TCGA cohort, we analyzed the differential expression of 169 glycosyltransferase-related genes in 502 tumors and 44 normal tissues and obtained a total of 58 DEGs (Figure 1A), of which 42 were upregulated and 16 were downregulated (Figure 1B).

Figure 1. Differentially expressed glycosyltransferase-related genes. (A) Volcano map of differential genes, the red nodes represent upregulated genes while the blue nodes represent downregulated genes. (B) Heatmaps showing differentially expressed genes in the TCGA dataset. The color indicates the level of expression of the gene (red represents upregulation, blue represents downregulation).

Establishment and validation of the glycosyltransferase-related gene prognostic signature

Univariate Cox analysis was performed on the candidate genes obtained in the above process, and a total of 12 genes exhibited statistically significant differences (P<0.05) (Supplementary Figure 1). To further identify independent prognostic genes in patients with HNSCC, multivariate Cox analysis was performed on 12 glycosyltransferase-related genes. As shown in Table 1, five significant genes (B4GALT3, PYGL, GALNT14, FUT2, and GALNT16) were screened to construct a prognostic signature. The risk score of all tumor samples was counted as follows: Risk score = (0.974* B4GALT3) + (0.418* PYGL) + (0.195* GALNT14) + (−0.413* FUT2) + (−0.886* GALNT16).

Using the median risk score as the cut-off value, all patients with HNSCC were classified into high- and low-risk subgroups (Figure 2A), and high-risk samples exhibited poorer overall survival (OS) in comparison to low-risk samples (log-rank P < 0.001, Figure 2B). ROC curves were used to assess the capability of the signature and the AUC at one, three, and five years in Figure 2C shows moderate sensitivity and specificity of the prognostic signature. Additionally, validation was performed using the GSE41613 dataset. The same formula was used to compute the risk score, and the 97 patients were classified into two subgroups based on the median risk score (Figure 2D). The K-M curve and AUC of the ROC curves showed significance similar to that of the TCGA cohort (Figure 2E, 2F). These results demonstrate the reliability and stability of the prognostic signature.

Figure 2. Construction of prognostic signature of five glycosyltransferase-related genes. (A–C) The prognostic signature was constructed in the TCGA dataset. (A) The distribution of risk score, OS status, and the heatmap of the expression profiles of signature genes, (B) K-M survival curves based on the prognostic signature, and (C) the AUCs of glycosyltransfer related gene signature. (D–F) Evaluate the performance of the prognostic signature in the GSE41613 dataset. (D) The distribution of risk score, OS status, and the heatmap of the expression profiles of signature genes, (E) K-M survival curves based on the prognostic signature, and (F) the AUCs of glycosyltransfer related gene signature.

Independent analysis of prognostic signature

Risk scores and other clinicopathological characteristics (age, sex, grade, and stage) were combined for independent analyses. Univariate and multivariate Cox analysis suggested that risk score was a significantly independent prognostic factor for OS (TCGA cohort univariate: hazard ratio (HR) = 1.990, 95% confidence interval (CI) = 1.529−2.589, P < 0.001, Figure 3A; TCGA cohort multivariate: HR = 2.184, 95% CI = 1.631−2.926, P < 0.001, Figure 3B). This was confirmed again in the GSE41613 validation cohort (Figure 3C, 3D).

Figure 3. The independence identification of the risk model. Univariate and multivariate Cox regression analysis has been performed in the TCGA cohort (A, B) and GSE41613 cohort (C, D).

GSEA analysis and mutant landscape

The results of GSEA revealed that the high-risk subgroup samples were usually enriched in transcription and protein regulation-related pathways, such as glycosaminoglycan biosynthesis of chondroitin sulfate, proteasomes, and ribosomes (Figure 4A). Groups with low risk tended to activate arachidonic acid metabolism, cell adhesion molecules (CAMs), chemokine signaling pathways, and primary immunodeficiency pathways, among others (Figure 4B).

Figure 4. Molecular characteristics of high- and low-risk groups. (A) Gene sets enriched in the high-risk group. (B) Gene sets enriched in the low-risk group. The mutation profile of the top 20 mutation genes in high-risk patients (C) and low-risk patients (D).

Next, we performed gene mutation analysis according to the glycosyltransferase-related gene signature and observed that the mutation count in the high subgroup was higher than that in the low subgroup. Moreover, the top 20 mutated genes in the high- (Figure 4C), and low-risk (Figure 4D) samples were contradistinguished. Additionally, we observed that missense variations were the most common mutation type, and all samples were dominated by mutations in TP53, titin (TTN), and FAT atypical cadherin 1 (FAT1) in both groups.

The associations with immune microenvironment and ICI therapy

To further explore the differences in immune cells between the two subgroups, the CIBERSORT method was used to quantify the immune infiltration scores, and the Wilcoxon test was performed for comparison. We found that low-risk samples were enriched in naive B cells, plasma cells, CD8 T cells, regulatory T cells (Tregs), and resting mast cells, whereas resting memory CD4 T cells, resting NK cells, M0 macrophages, and activated dendritic cells were more abundant in high-risk samples (Figure 5A). We then calculated the Spearman correlation coefficient between the risk score and common immune checkpoint, the outcomes suggested that the risk score was positively correlated to fibroblast activation protein (FAP) and lysyl oxidase-like 2 (LOXL2), and negatively correlated with cytotoxic T lymphocyte-associated protein 4 (CTLA4), ICOS, and programmed cell death protein 1 (PDCD1) (Figure 5B).

Figure 5. Correlation of tumor microenvironment (TME) landscape and ICI therapy with risk score. (A) The proportions of TME cells in high- and low-risk groups. Significant statistical differences between the two subgroups were assessed using the Wilcoxon test. (B) Linear regression among immune checkpoints (CTLA4, FAP, ICOS, LOXL2, and PDCD1) and risk scores, and the numbers placed on the right of the plot represent coefficients. (C) TIDE, MSI, and T cell exclusion and dysfunction score in two groups. (*p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001).

In addition, we used the tumor immune dysfunction and exclusion (TIDE) and microsatellite instability (MSI) scores to evaluate the potential value of immunotherapy. In our study, the high-risk group had lower TIDE and MSI scores, which suggested that this group may benefit significantly from immunotherapy than the low-risk group (Figure 5C). Finally, the low-risk group had lower T cell exclusion scores and higher T cell dysfunction.

Diagnostic model construction

After evaluating the prognostic value, we performed gene selection for diagnosis. 37 glycosyltransferase-related genes were extracted by SVM and 18 genes were screened by LASSO analysis, and after intersecting the two algorithms, we obtained 15 genes with diagnostic value (Supplementary Table 2). A diagnostic model was constructed based on the corresponding coefficients. In the TCGA cohort, the AUC value of the ROC curve was 0.997 (95CI%: 0.948–0.999), which suggested excellent diagnostic efficiency of the model (Figure 6). In addition, the above diagnostic model was also verified in GSE127165 (AUC:0.978, 95CI%: 0.947–0.991).

Figure 6. Establishment of multigene diagnostic signature. A total of 15 genes have been selected through SVM and LASSO analysis, and the AUC of the ROC curve is 0.997 (TCGA) and 0.978 (GSE127165).

B4GALT3 promotes tumor progression in vitro

To elucidate the biological functions of B4GALT3 in HNSCC, we created a stable TU177 cell line with a B4GALT3 knockdown. In contrast, we ectopically overexpressed the B4GALT3 gene in the FaDu cells. Western blotting was performed to assess the effectiveness of B4GALT3 knockdown and overexpression in these cell lines (Figure 7A). CCK-8 and colony formation assays revealed a notable decrease in the cell viability and colony-forming ability of TU177 cells upon knockdown of B4GALT3 (Figure 7B, 7D, 7E). Conversely, the overexpression of B4GALT3 in FaDu cells stimulated cell proliferation and enhanced colony formation (Figure 7C, 7F, 7G). Subsequently, we investigated whether B4GALT3 plays a role in HNSCC cell motility. Wound healing and Transwell assays revealed a significant enhancement in cell migration and invasion following B4GALT3 overexpression, whereas these migratory and invasive properties decreased upon B4GALT3 knockdown (Figure 7H–7O).

Figure 7. B4GALT3 promotes the proliferation and migration of HNSCC cells. (A) The expression of B4GALT3 has been assessed by western blotting. (B, C) CCK-8 assay has been performed to evaluate the growth rates of the indicated cells. The proliferative and migratory abilities of the indicated cells have been measured by colony formation (D–G) and wound healing assays (H–K), respectively. (L–O) Cell migration and invasion abilities of the indicated cells have been measured using Transwell assays.