This article is mentioned in:

Introduction

In total, ~80% of malignant brain tumors are gliomas (1). The histological subtypes of glioma include astrocytoma, oligodendroglioma, ependymoma and mixed oligoastrocytoma. The World Health Organization (WHO) has also classified glioma into four grades according to the degree of malignancy (2). Glioblastoma (GBM), which has the highest frequency of occurrence, is the most aggressive cerebral malignancy (3). The prognosis of patients with GBM is highly unfavorable, with a median survival time of 1 year (4). Based on specific genetic mutation profiles, GBM can be further subdivided into four molecular subtypes: Proneural, mesenchymal, neural and classical (5). Furthermore, other molecular subpopulations of GBM have been proposed by Lee et al (6). Numerous studies have indicated that various signaling pathways act in the progression of glioma, including purinergic, receptor tyrosine kinase, TGF-β and STAT signaling pathways (7).

Ubiquitin-specific protease 18 (USP18) (enzyme commission: 3.4.19.12), originally termed UBP43, was primitively discovered as a differentially expressed gene in mice with acute myeloid leukemia (AML) carrying the AML t(8;21) gene translocation (8). Subsequently, Zhang et al (9) and Kang et al (10) verified that USP18 exists in virus-infected porcine alveolar macrophages and human melanoma cells, respectively. USP18 is expressed in multiple tissues; for example, USP18 is upregulated in the spleen and liver, while it is downregulated in the lungs and bone marrow (8). At present, the known primary function of USP18 is the role it has in the de-ISGylation mediated by interferon (IFN)-stimulated gene 15 (ISG15) (11). Furthermore, USP18 inhibits ISGylation-independent IFN-I signaling through the IFNAR2/JAK/STAT pathway and it mediates the inflammatory response accordingly (12,13). The mechanism by which USP18 inhibits the IFN-I response to promote cancer progression has been elucidated in detail (14). Furthermore, USP18 is involved in several mechanisms involved in carcinogenesis, such as in promoting breast cancer growth and epithelial-mesenchymal transition (EMT) in glioma (15).

In the present study, a detailed analysis of the role of USP18 in glioma was performed. USP18 was upregulated in glioma and was associated with clinical traits such as tumor stage and a worse prognosis. USP18 was also associated with the IFN-I pathway in glioma T cells. In glioma, USP18 expression was also associated with the proportion of a number of immune cells. Genes co-expressed with USP18 were mainly enriched in the immune pathways such as IFN-I and complement activation. Cell Counting Kit-8 (CCK-8) experiments demonstrated that USP18 knockdown reduced glioma cell viability. Furthermore, colony formation assays demonstrated that knockdown of USP18 significantly reduced the proliferation of glioma cells. Specifically, the results of the present study suggested that USP18 may influence glioma initiation and progression, with potential regulatory effects on IFN-I signaling pathways.

Materials and methods

Data preparation

Transcripts for 33 types of cancer were obtained from the University of California Santa Cruz (UCSC) database (https://xenabrowser.net/datapages/), all of which were uniformly processed using the TOIL workflow (16). All transcriptome profiles [in high throughput sequencing (seq) transcripts per million (TPM) formats] and clinical information were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) lower grade glioma (LGG)/GBM projects. The immunophenoscore (IPS) of TCGA glioma samples was obtained from The Cancer Immunome Atlas (https://tcia.at/) (17). The transcript and clinical data of glioma samples from TCGA with complete survival information and clinical features were extracted using R software (https://www.R-project.org/) (v.4.2.1). The RNAseq and clinical information of 693 and 325 glioma samples were acquired from the Chinese Glioma Genome Atlas [CGGA; Dataset IDs, mRNAseq_693 (CGGA_693) and mRNAseq_325 (CGGA_325); http://www.cgga.org.cn/download.jsp] as additional test sets. Visualization of USP18 expression in the pan-cancer analysis was conducted using the R package, ‘ggplot2’ (https://ggplot2.tidyverse.org/).

Differential expression and clinical correlation analyses

A pan-cancer analysis of USP18 expression in 33 types of cancer and the corresponding normal tissues was conducted using the R package, ‘ggplot2’. To perform clinical correlation analyses, all patients with glioma were separated into subgroups based on clinical features and were analyzed using R software (18).

Single-cell sequencing analysis

The single cell RNAseq glioma tissue dataset, GSE131928_10X and the glioma T cell dataset, GSE163108_10X, were obtained from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo) (19). The single-cell expression matrix, Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction and cell type metadata were downloaded from the Tumor Immune Single-cell Hub 2 (TISCH2) database (http://tisch.comp-genomics.org/) (20). The dimensionality reduction and cell annotation of single-cell matrices were conducted using the R package, ‘Seurat’ (21). The UMAP visualization, violin plots, bubble plots and differential expression analyses were also conducted using the R package, ‘Seurat’. Differential expression analyses were performed in cell groups with specific USP18 expression patterns. The threshold for differential genes was set to |log fold change (FC)|>0.5 and adjusted P<0.05. Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) analysis was conducted using the R package, ‘clusterProfiler’ (22). Histogram visualization was completed using the R package, ‘ggplot2’. The list of USP18 co-expressed genes in CD4 and CD8 T cells was obtained from TISCH2 and the correlations between genes with an average logTPM of >0.5 or a maximum logTPM of >2 were calculated.

Prognosis analysis

Kaplan-Meier (KM) curves and Cox regression analyses were performed using the R package, ‘survival’ (23). The prognostic data used for KM survival curve analysis were obtained from a previous study (23) and samples were separated into corresponding cohorts based on the median expression value of USP18. The time-dependent receiver operating characteristic (ROC) curve was constructed using the R package ‘timeROC’ (https://CRAN.R-project.org/package=timeROC). The two aforementioned CGGA datasets served as validation sets in the KM survival analysis. The nomogram was constructed with the R package, ‘regplot’ (https://CRAN.R-project.org/package=regplot). The bootstrapping method was used and repeated 1,000 times to validate the nomogram for 1, 3 and 5 years. All validation was conducted using the R package, ‘rms’ (https://CRAN.R-project.org/package=rms).

Tumor immune profile analysis

The checkpoint inhibitors and major histocompatibility complex (MHC) molecule genes were acquired from the study by Charoentong et al (17). The relationship between the USP18 expression level and the score of 24 immune cells was assessed using the R package, ‘GSVA’ (24). The immune correlation analysis used data from 24 types of immune cell markers (25). The immune microenvironment analysis of the stromal score, estimate score, immune score, tumor purity and USP18 expression was conducted with the R package, ‘estimate’ and the score was calculated as previously described by Yoshihara et al (26). The USP18 expression stratified gene mutation files were differentiated using the Perl language and visualized using the R package, ‘maftools’ (https://github.com/PoisonAlien/maftools).

Functional enrichment analyses

Gene Ontology (GO), KEGG and Gene Set Enrichment Analysis (GSEA) functional analyses for the differentially expressed genes (DEGs) were performed using the R packages, ‘enrichplot’ (https://github.com/YuLab-SMU/enrichplot), ‘org.Hs.eg.db’ (https://bioconductor.org/packages/org.Hs.eg.db/) and ‘clusterProfiler’ (https://guangchuangyu.github.io/software/clusterProfiler/) (27-29). The DEGs were selected with the R package, ‘DESeq2’ (https://github.com/thelovelab/DESeq2) and the thresholds were |logFC|>1.5 and adjusted P<0.05. All DEGs were included in the GSEA, using the USP18 low-expression group as a reference.

Cell culture

U87MG ATCC (glioblastoma of unknown origin) and U251 cell lines were purchased from the Chinese National Infrastructure of Cell Line Resource and cultured in DMEM supplemented with 10% fetal bovine serum, 2 mM glutamine and 1% penicillin-streptomycin (all from Gibco; Thermo Fisher Scientific, Inc.). The cell lines were cultured in a 37˚C incubator containing 5% CO2. U87MG ATCC was authenticated using STR profiling.

USP18 knockdown

Lipofectamine RNAiMAX (Invitrogen; Thermo Fisher Scientific, Inc.) was used for transfecting USP18 targeting and control (Ctrl; 40 nM) small interfering RNAs (siRNAs; 40 nM) into U87MG ATCC and U251 cells, according to the manufacturer's protocols and as previously described (30). The siRNA sequences used were as follows: siUSP18#1, 5'-GGAAUUCACAGACGAGAAATT-3'; siUSP18#2, 5'-GGAAGAAGACAGCAACAUGTT-3'; and siCtrl, 5'-UUCUCCGAACGUGUCACGU-3'. The USP18 knockdown efficiency was tested by reverse transcription-quantitative PCR (RT-qPCR) analysis. The subsequent experiments were conducted 48 h after transfection.

RNA isolation and RT-qPCR

Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Inc.). The RNA was reverse transcribed using the reverse transcription kit (cat. no. A3500; Promega Corp.) according to the manufacturer's instructions, and qPCR assays were performed with SYBR Green Master Mix (Takara Bio, Inc.) on a 7500 fast (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The thermocycling conditions were as follows: Pre-denaturation at 95˚C for 5 min and 40 cycles of denaturation at 95˚C for 15 sec, annealing at 60˚C for 20 sec and elongation at 72˚C for 30 sec. The following primers were used for qPCR: USP18 forward, 5'-CCTGAGGCAAATCTGTCAGTC-3' and reverse, 5'-CGAACACCTGAATCAAGGAGTTA-3'; and β-actin forward, 5'-CATGTACGTTGCTATCCAGGC-3' and reverse, 5'-CTCCTTAATGTCACGCACGAT-3'. The 2-ΔΔCq method was used to quantify the fold changes in gene expression (31).

CCK-8 assay

Cell viability was assessed using CCK-8 (Sigma-Aldrich; Merck KGaA), following the manufacturer's protocols as previously described (32).

Colony formation

The colony formation assay was used to examine the cell colony formation ability of U87MG ATCC and U251 cells following USP18 expression knockdown. The colony formation assay was performed as previously described (32). Clusters of >50 cells were regarded as colonies. The colonies were calculated using ImageJ software (v.1.47; National Institutes of Health).

Apoptosis assay

The transfected U87MG ATCC and U251 cells were cultured for 48 h at 37˚C. Subsequently, the cells were thoroughly mixed with propidium iodide (5 µl) and binding buffer (500 µl), and then incubated with fluorescein isothiocyanate-conjugated anti-annexin V antibody (5 µl) at room temperature for 10 min, according to the manufacturer's protocols. Finally, the cells were evaluated using a CytoFLEX LX (Beckman Coulter, Inc.) and the data were analyzed using CytTExpert software (v.2.4; Beckman Coulter, Inc.).

Statistical analysis

The gene expression data were presented as the mean ± standard deviation and the survival hazard ratios (HRs) were presented as HR (95% confidence interval). The counting and statistical analyses were repeated three times. The cell line experimental data were analyzed using GraphPad Prism 8.0 software (Dotmatics). The age, isocitrate dehydrogenase (IDH) mutation and 1p19q co-deletion status of different groups were compared using the Wilcoxon rank-sum test. The WHO grade and histological type were compared using Kruskal-Wallis test with Dunn post hoc tests. Survival analysis was conducted using the log-rank test. Gene correlation analysis was conducted using Spearman's correlation analysis. For multiple comparisons vs. the same control group, including mRNA levels, cell viability, colony numbers and apoptosis percentage, one-way ANOVA followed by Dunnett's test was used. P<0.05 was considered to indicate a statistically significant difference.

Results

Differential expression and clinical correlation analyses of USP18 in glioma

As indicated in Fig. 1, A pan-cancer analysis of USP18 expression using data from TCGA was performed and it was found that USP18 was upregulated in most cancer types (Fig. 1A). To ensure TCGA samples were sufficient and comparable, the analysis was repeated using data from the UCSC. USP18 was significantly differentially expressed in various types of cancers in the UCSC dataset, with high expression observed in LGG and GBM and low expression observed in adrenocortical carcinoma and kidney chromophobe (Fig. 1C). Furthermore, the expression level of USP18 in the GBM dataset from TCGA was examined and an upregulation of USP18 was observed (Fig. 1B). The same result was observed in the LGG and GBM datasets from the UCSC (Fig. 1D). Patients with LGG or GBM, aged >60 years, with wild-type IDH and without the 1p/19q co-deletion had a higher expression of USP18 (Fig. 1E-G). However, the expression level of USP18 revealed no sex-specific differences (Fig. 1H). In more advanced glioma WHO grades, the USP18 expression level was significantly higher (Fig. 1I). In addition, patients with GBM expressed notably higher USP18 levels than patients with other histological types and the expression of USP18 was not significantly different between patients with oligoastrocytoma, astrocytoma or oligodendroglioma (Fig. 1J).

Figure 1 Differential expression and clinical correlation analyses of USP18 in glioma. (A) USP18 expression in pan-cancer using The Cancer Genome Atlas data. Red indicated tumor samples and blue indicated normal samples. (B) USP18 expression in GBM using TCGA data. (C) USP18 expression in pan-cancer using University of California Santa Cruz data. (D) USP18 expression in LGG and GBM using UCSC data. Clinical correlation analyses comparing (E) age in years, (F) IDH status, (G) 1p/19q co-deletion, (H) sex, (I) WHO grade and (J) histological type with USP18 expression in patients with LGG and GBM. Red indicates patients aged >60 years, IDH mutation, 1p/19q non-co-deletion, male, WHO grade 3 and oligodendroglioma. Blue represents patients aged <60 years, IDH wild-type, 1p/19q co-deletion, female, WHO grade 2 and oligoastrocytoma. Green represents groups with WHO grade 4 and astrocytoma, while purple indicates the GBM cohort. *P<0.05, **P<0.01 and ***P<0.001; ns, no significance. The abbreviation for cancers in A and C may be found at https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/tcga-study-abbreviations. USP18, ubiquitin-specific protease 18; N, normal tissue; T, tumor tissue; WT, wild-type; Mut, mutant; codel, codeletion; GBM, glioblastoma; LGG, lower grade glioma; IDH, isocitrate dehydrogenase; WHO, World Health Organization.

USP18 expression patterns in the glioma single-cell landscape

The glioma tissue single-cell dataset was selected for dimensionality reduction (Fig. 2A). The overall expression level of USP18 was low in this dataset and no characteristic cell expression cluster was observed (Fig. 2B). The USP18 expression level in different cell clusters was further compared (Fig. 2C) and it was determined that USP18 was significantly upregulated in terminally exhausted CD8 T cells (CD8 Tex). Since the USP18 negative regulation mechanism of the IFN-I pathway in cancer has been demonstrated, it was hypothesized that USP18 might play a specific role in glioma T cells. Therefore, the glioma T cell dataset was selected for further analysis of the USP18 expression pattern (Fig. 2D), and it was demonstrated that USP18 was highly expressed in a small number of CD8 and CD4 T cells (Fig. 2E). However, there was no significant difference in the USP18 expression level in each of the cell clusters shown in the violin plot (Fig. 2F). The pathways enriched in the CD8 Tex subcluster of the glioma tissue dataset were the T cell subtype differentiation and T cell receptor-related pathways, antigen presentation and programmed cell death protein 1 in cancer (Fig. 2G). The genes co-expressed with USP18 in the CD4 (Fig. 2H) and CD8 (Fig. 2I) T cell subsets of the glioma T cell dataset included the IFN-induced protein with tetratricopeptide repeats (IFIT) family, the 2'-5'-oligoadenylate synthetase family and the MX dynamin-like GTPase family.

Figure 2 USP18 expression patterns in glioma tissues and glioma T cells. UMAP visualization of cell types in (A) glioma tissue and (D) glioma T cells. UMAP visualization of USP18 expression in (B) glioma tissue and (E) glioma T cells. Violin plots of USP18 expression in (C) different cell types in glioma tissue and (F) glioma T cells. (G) Kyoto Encyclopedia of Genes and Genomes enrichment analysis of the CD8 Tex cluster in the high USP18 expression cohort. USP18 co-expressed genes in the (H) CD4 T cell and (I) CD8 T cell clusters. USP18, ubiquitin-specific protease 18; UMAP, Uniform Manifold Approximation and Projection; CD8 Tex, exhausted CD8 T cells.

Upregulation of USP18 is related to a worse survival outcome in patients with glioma

A KM survival curve analysis based on the USP18 expression levels in TCGA cohort was conducted and the results demonstrated that patients with glioma and upregulated USP18 expression had a worse overall survival (OS) (Fig. 3A), disease-specific survival (DSS) (Fig. 3B) and progression-free interval (PFI) (Fig. 3C). The area under the curve values at 1, 3 and 5 years for OS (Fig. 3D), DSS (Fig. 3E) and PFI (Fig. 3F) were all >0.6. A validation analysis of the KM survival curves using the CGGA_325 (Fig. 3G) and CGGA_693 (Fig. 3H) datasets was then performed. Patients with high USP18 expression in both validation sets had a worse prognosis. The 5-year prognostic predictive power was the best for OS, DSS and PFI. To build a prognostic model to predict survival risk, a prognostic nomogram was constructed at 1, 3 and 5 years based on USP18 expression and the patient clinicopathological parameters (Fig. 3I). Consequently, the constructed calibration curves had a satisfactory consistency at 1, 3 and 5 years (Fig. 3J).

Figure 3 Survival and ROC analyses for different USP18 expression levels. (A) OS, (B) DSS and (C) PFI analyses of patients with GBM or LGG with high or low USP18 expression. Red indicates the high USP18 expression cohorts and blue represents the low USP18 expression cohorts. Time-dependent ROC curve analysis at 1, 3 and 5 years for (D) OS, (E) DSS and (F) PFI. Blue indicates the ROC curve for USP18 expression at 1 year, red indicates 3 years and green indicates 5 years. Kaplan-Meir survival analysis according to the USP18 expression level in the (G) CGGA_325 and (H) CGGA_639 test cohorts. Red indicates the high USP18 expression cohorts and blue indicates the low USP18 expression cohorts. (I) Nomograms at 1, 3 and 5 years based on USP18 expression and clinicopathological parameters. (J) Calibration curve analysis at 1, 3 and 5 years. Blue indicates the calibration curve of USP18 at 1 year, red indicates 3 years and green indicates 5 years. The grey line indicates the ideal curve. ROC, receiver operating characteristic; USP18, ubiquitin-specific protease 18; OS, overall survival; DSS, disease-specific survival; PFI, progression-free interval; GBM, glioblastoma; LGG, lower grade glioma; CGGA, Chinese Glioma Genome Atlas; TCGA, The Cancer Genome Atlas; WHO, World Health Organization; IDH, isocitrate dehydrogenase.

Crosstalk of USP18 expression on the glioma immune microenvironment and immunotherapy

The MHC molecule (Fig. 4A), Immunostimulator (Fig. 4B) immunoinhibitor (Fig. 4C), chemokine (Fig. 4D) and chemokine receptor (Fig. 4E) gene levels were compared between the two USP18 expression cohorts and it was revealed that most were upregulated in the high USP18 expression cohort. However, the immunostimulatory gene, TNF superfamily member 9 and the chemokine receptor genes, C-C motif chemokine ligand (CCL)3 and CCL19, were downregulated in patients with high USP18 expression. Overall, the proportions of different immune cell types, such as activated dendritic cells (DCs), macrophages, T cells, neutrophils and natural killer (NK) CD56dim cells, were significantly positively correlated with the USP18 expression level, while the other immune cell types were negatively correlated, such as plasmacytoid DCs and NK CD56-bright cells (Fig. 4F). Furthermore, the high USP18 expression cohort was characterized by a high stromal score, immune score and low tumor purity (Fig. 4G-J). The frequency of PTEN and EGFR mutations in the high USP18 expression cohort was also higher (Fig. 5A). However, the low USP18 expression cohort had a higher IDH1 and ATRX mutation frequency (Fig. 5B). A higher tumor mutation burden (TMB) was observed in the high USP18 expression group, which indicated tumor neoantigen potential (Fig. 5C). It is also worth noting that TMB had a prognostic effect on glioma (Fig. 5D) and patients with a high TMB and high USP18 expression had a worse prognosis (Fig. 5E). Furthermore, patients with glioma and upregulated USP18 expression and pos_PD1 had a higher IPS, suggesting a higher immunotherapeutic potential (Fig. 5F-I).

Figure 4 Correlation analysis between glioma immune cell infiltration and USP18 expression. Differential analysis of (A) major histocompatibility complex molecules, (B) immunostimulatory, (C) immune checkpoints, (D) chemokines and (E) chemokine receptors in the low and high USP18 expression groups. Red indicates the high USP18 expression group and blue indicates the low USP18 expression group. (F) Lollipop plot of the correlation between the proportion 24 immune cell types and USP18 expression. The size of the circle represents the correlation coefficient, and the color indicates the P-value. The P-value decreases as the color becomes redder. (G-J) Immune microenvironment analysis based on USP18 expression, including the (G) estimate score, (H) stromal score, (I) immune score and (J) tumor purity. *P<0.05, **P<0.01 and ***P<0.001. USP18, ubiquitin-specific protease 18.

Figure 5 Effect of USP18 expression on immunotherapy of glioma. Visualization of the gene mutation frequency in the USP18 (A) high and (B) low expression groups. (C) Differences in TMB between samples with high and low USP18 expression levels. The prognostic difference between (D) TMB and (E) TMB combined with USP18 expression in glioma samples. (F-I) Differences in immunotherapy scores between the high and low USP18 expression groups in patients with (F) neg-CTLA4 and pos-PD1, (G) pos-CTLA4 and pos-PD1, (H) neg-CTLA4 and neg-PD1 and (I) pos-CTLA4 and neg-PD1. *P<0.05, **P<0.01 and ***P<0.001. USP18, ubiquitin-specific protease 18; TMB, tumor mutation burden; neg, negative; pos, positive; CTLA-4, cytotoxic T-lymphocyte associated protein 4; PD1, programmed cell death protein 1.

USP18 is associated with multiple signaling pathways in glioma

Statistically significant DEGs between the USP18 low and high expression cohorts in the LGG and GBM datasets were identified, which included 462 upregulated and 175 downregulated genes. GO analysis chord diagrams incorporated the 10 most significant GO terms (P<0.001), while the KEGG pathway analysis chord diagram incorporated the 8 most significant pathway terms (P<0.001). The term, ‘biological process’, included genes enriched in the ‘pattern specification process’ and the ‘development of the embryonic skeletal system’ (Fig. 6A). The term, ‘cellular component’, included genes enriched in the ‘transmembrane transporter complex’, an integral component of the postsynaptic membrane (Fig. 6B). While the ‘molecular function’ term included genes enriched in ‘ligand-gated ion channel activity’ and ‘ion-gated channel activity’ (Fig. 6C). The KEGG analysis revealed an association with ‘neuroactive ligand-receptor interaction’ and ‘taste transduction’ (Fig. 6D). Next, all DEGs were incorporated into GSEA. The main enriched pathways in the high USP18 expression group were ‘IFN alpha beta signaling’, ‘transcriptional regulation of granulopoiesis’ and ‘initial triggering of complement’ (Fig. 6E). The pathways enriched in the low USP18 expression group were the ‘neuronal system’, ‘neurotransmitter release cycle’ and ‘neurotransmitter receptors and postsynaptic signal transmission’ (Fig. 6F). Genes co-expressed with USP18 included STAT1 and OSAS3 (Fig. 6G), and the genes negatively correlated with USP18 were mainly mitochondrially encoded genes (Fig. 6H).

Figure 6 Functional enrichment analyses of DEGs between two USP18 cohorts. (A-D) Functional enrichment analysis chord diagrams for GO enrichment analyses terms, (A) BP, (B) CC and (C) MF and (D) Kyoto Encyclopedia of Genes and Genomes analysis of the DEGs. (E and F) Gene Set Enrichment Analysis for the visualization of pathways enriched in the USP18 (E) high and (F) low expression groups. (G and H) Heatmaps of genes significantly (G) positively and (H) negatively correlated with USP18 expression. DEGs, differentially expressed genes; USP18, ubiquitin-specific protease 18; GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

Knocking down USP18 expression suppresses proliferation and colony formation and induces the apoptosis of glioma cells

To investigate the biological role of USP18 in the development of glioma, two siRNAs targeting USP18 expression (siUSP18#1 and siUSP18#2) were transfected into glioma cells to knockdown USP18 expression. The RT-qPCR validation assay confirmed that USP18 expression was significantly reduced in U87MG ATCC and U251 cells transfected with siUSP18 (Fig. 7A). Next, the CCK-8 assay was used to evaluate the impact of USP18 knockdown on glioma cell proliferation. The viability of U87MG ATCC and U251 cells following USP18 knockdown was significantly reduced on days 3 and 4 (Fig. 7B). Analogous results were obtained in the colony formation assay. Compared with the siCtrl group, the number of colonies in the siUSP18#1 and siUSP18#2 groups decreased significantly (Fig. 7C). Moreover, it was also found that USP18 knockdown significantly induced the apoptosis of U251 and U87MG ATCC cells (Fig. 7D). Collectively, these results indicated that USP18 knockdown inhibited the proliferation and induced the apoptosis of glioma cells.

Figure 7 Knockdown of USP18 expression impairs the proliferation of glioma cells. (A) USP18 was significantly downregulated in U87MG ATCC and U251 cells following siUSP18 transfection, as demonstrated by reverse transcription-quantitative PCR. β-actin was used as the internal control. (B) Cell viability of U87MG ATCC and U251 cells transfected with siCtrl, siUSP18#1 or siUSP18#2, as determined by Cell Counting Kit-8 assay. (C) Colony formation assay results of U87MG ATCC and U251 cells transfected with siCtrl, siUSP18#1 or siUSP18#2. (D) Apoptosis assay results of U87MG ATCC and U251 cells transfected with siCtrl, siUSP18#1 or siUSP18#2. *P<0.05, **P<0.01 and ***P<0.001. USP18, ubiquitin-specific protease 18; si, small interfering; Ctrl, control.

Discussion

The role of USP18 in viral infectious diseases, particularly those related to the central nervous system (CNS), has been the focus of previous research (33,34). The role of USP18 in the IFN pathway has therefore encouraged research in the field of malignant tumors (35). Since USP18 was first identified in human melanoma cell lines in 2001, research on the effects of USP18 on cancer has increased (8). There is also growing evidence that USP18 is upregulated in melanoma and liver cancer (36-39). A pan-cancer analysis of USP18 expression was also conducted in the present study. Although the function of USP18 in the modulation of signaling pathways in glioma has been previously described (7), its role has not been exhaustively elucidated. Therefore, in the present study, the high expression of USP18 in glioma was verified using TCGA and UCSC database and these findings were consistent with a previous study (15). The upregulation of USP18 expression may be induced by the IFN-I response elicited by glioma cells, which is similar to other IFN-induced genes (15). The mechanism may involve free ISG15 blocking S-phase kinase-associated protein-Cullin-F-box protein-mediated ubiquitination and subsequent proteasomal degradation of USP18(40). In the present study, USP18 expression in glioma was further verified using single-cell data and in the glioma microenvironment, significantly high expression of USP18 in CD8 T cells was observed. The CD8 T cell cluster is associated with IFN-I signaling and T cell differentiation. Furthermore, USP18 has been revealed to regulate T cell activation and T helper 17 cell differentiation through deubiquitination of the TAK1-TGF-β activated kinase 1 complex (41). This suggests that USP18 may play a similar role in regulating differentiation and inhibiting IFN-I signaling in glioma T cells, thereby promoting glioma progression. In the present study, USP18 co-expressed genes were further identified in the glioma T cell subgroups and most of these genes were IFN-inducible. For example, the IFIT genes are important mediators of innate immunity (42), and the IFN-induced MX genes enhance endoplasmic reticulum stress-related cell death (43). USP18 localizes to the nucleus, is recruited to IFNAR2 by STAT2 and interferes with the assembly of the ternary IFN-IFNAR1-IFNAR2 complex, thereby acting as a negative regulator of the IFN-I signaling pathway (44).

Glioma is characterized by heterogeneity and differences in glioma classification, prognosis and treatment based on clinical features have been widely discussed (45,46). To verify whether USP18 expression is associated with clinical characteristics, a clinical correlation analysis was performed in the present study and it was found that USP18 was highly expressed in patients >60 years old, with wild-type IDH, non-co-deletion of 1p19q, WHO grade 3 or 4 and GBM. A survival analysis of all three dataset cohorts was also performed and it demonstrated that USP18 upregulation was related to poor survival outcomes.

Since the tumor immune microenvironment has been deemed a key element involved in tumor development, immune infiltration and differences in immune-related gene expression were evaluated in the present study (47). It was found that the proportion of various immune cells was related to the USP18 expression level, and the immuno-genes were mostly upregulated in the high USP18 expression subgroup. Since functional lymphangion was described in the meninges in 2015, the tumor microenvironment (TME) in GBM has been considered an immunologically distinct organ (48,49). The TME of GBM consists of microglia (permanent macrophages of the CNS), macrophages (peripheral bone marrow-derived macrophages), tumor-infiltrating lymphocytes and tumor-infiltrating DCs (47,50,51). Microglia and macrophages in glioma are tumor-associated macrophages (TAMs) and represent distinctly different cell groups (52). Furthermore, TAMs have been demonstrated to promote GBM infiltration and GBM proliferation (47). Neutrophils in GBM, or tumor-associated neutrophils, excrete elastase to boost tumor proliferation and angiogenesis (53,54). The results of the present study demonstrated that the USP18 expression level was positively correlated with macrophage and neutrophil infiltration, indicating that USP18 may promote GBM progression through this infiltration. T cells engaged in the adaptive immune response can restrict glioma growth (55). The nuclear factor-κ light chain enhancer of activated B cells (NF-κB) pathway, a T cell activating pathway, has also been reported (56). Regulatory T cells, which express the Forkhead Box P3 transcription factor, downregulate the NF-κB pathway to promote immunosuppression in the adaptive immune response (57,58). In addition, USP18 has been revealed to suppress the NF-κB pathway in innate immune cells (59). Thus, further research is necessary to explore the USP18 and NF-κB pathways in the adaptive immunological response. However, these results suggested that USP18 may crosstalk with immunotherapy in glioma. In the present study, the potential impact of USP18 expression on immunotherapy was explored. The results suggested that patients with glioma with high USP18 expression may exhibit improved response to immunotherapy.

GO, KEGG and GSEA pathway analyses of DEGs were conducted in the present study and it was found that the GO terms, CC and MF, were both enriched in intercellular signal transduction processes, such as those involving the ion channel complex and ligand-gated ion channel activity. KEGG analysis revealed a similar result, with neuroactive ligand-receptor interactions being the most significant term. The GSEA indicated that the DEGs were mainly enriched in gene sets associated with the IFN, complement activation signaling pathways. Further co-expression analysis indicated that USP18 may play a role in the regulation of IFN signaling and mitochondrial oxidative respiratory chain transmission. Mitochondrial-encoded genes involved in the regulation of apoptosis in glioma cells have also been reported (60). The respiratory chain transmission and oxidative phosphorylation metabolism associated with these genes may provide the feasibility of an in-depth study of the mechanism of USP18 in glioma.

IFN-I are the largest IFN class of the three types (61). IFN-I is regularly used as a biotherapeutic agent in clinical malignancy as it can directly restrict tumor growth and progression and improves the immune surveillance against cancer (62-64). However, GBM is strongly resistant to IFN-I-induced apoptosis and therefore IFN-I achieves little efficacy in GBM therapy (65-68). TNF-related apoptosis-inducing ligand (TRAIL), which exerts specific pro-apoptotic activity against transformed cells, belongs to the TNF superfamily (69,70). Downregulation of USP18 was initially reported to enhance apoptosis induced by IFN-α or TRAIL in different cell lines, activating the extrinsic apoptosis pathway and the associated upregulation of TRAIL (71). In addition, USP18 can also maintain resistance against IFN-I treatment of GBM, which is TRAIL-dependent (72). These findings provided a unique perspective on the resistance of GBM to IFN-I, and USP18 may play a key role. GBM currently lacks effective therapies (73), and the first-line drugs for the treatment of GBM are temozolomide and bevacizumab. Research to identify targeted agents for future immunotherapy, gene therapy and antiangiogenic strategies are required to combat GBM tumor heterogeneity (74).

In the present study, to determine whether USP18 exerted an effect on glioma proliferation, RT-qPCR, CCK-8, colony formation and apoptosis experiments were performed. Knockdown of USP18 expression suppressed the proliferation capacity and colony formation of glioma cells. USP18-mediated glioma proliferation may be related to EMT. EMT maintains epithelial cells in an unstable quasi-mesenchymal cell state (75,76). Although the detailed mechanisms underlying EMT in GBM remain unknown, evidence indicates that EMT is engaged in the invasion and growth of GBM (77). Furthermore, USP18 has been discovered to stabilize TWIST and its expression, which may contribute to promoting EMT in GBM (15). A recent study has also indicated that USP18 affects EMT and promotes glioma stem cell growth (78).

This study still has limitations. For instance, the specific mechanism by which USP18 affects glioma cell proliferation was not elucidated and the experiments were only conducted at the cellular level, requiring further verification by animal experiments, such as in vivo tumorigenicity assay in athymic nude mice.

In conclusion, the present study indicated that USP18 may be involved in glioma invasion, growth, metastasis, immune response and tolerance. USP18 may be valuable as a latent outcome indicator and a latent therapy target in glioma. Furthermore, the results of the present study provided additional information regarding the molecular mechanisms involving USP18 in glioma. Therefore, future studies should focus on the role of USP18 in glioma.

Acknowledgements

Not applicable.

Funding

Funding: The present study was supported by the National Natural Science Foundation of China (grant no. 81201991), the Basic Research Program of Shanxi Province (grant no. 202103021224400), the Key Medical Scientific Research Projects of Shanxi (grant no. 2021XM35) and the Four ‘Batches’ Innovation Project of Invigorating Medical through Science and Technology of Shanxi (grant no. 2022RC07).

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in TCGA (https://portal.gdc.cancer.gov/), UCSC Xena database (https://xenabrowser.net/datapages/), CGGA (http://www.cgga.org.cn/) and TISCH2 (http://tisch.comp-genomics.org/) repositories.

Authors' contributions

YC, RL and ZL conceived and designed the study and conducted the experiments. ZL and PL acquired the data. JL, ZZ and JH analyzed the results. YC, RL and BY operated software and visualized the results. YW, YZ and GG supervised the experiments and validated the results. YC and RL confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References