The dataset for this study consisted of the mRNA transcriptome data and the clinical information of 169 glioblastoma tissues, as well as the mRNA transcriptomic data and the clinical information of 207 normal brain tissues. The former was downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov), while the latter from the Genotype-Tissue Expression (GTEx) database (https://commonfund.nih.gov/gtex).

The 207 normal brain tissues were set as control group, and the 169 glioblastoma tissues were set as experimental group. The limma package in R was used to screen differential genes under the threshold of log2FC absolute value > 1 and a P value < 0.05 (15).

The Gene Expression Profiling Interactive Analysis (GEPIA) database is a web tool based on the TCGA and GTEx data for cancer and normal gene expression profiles as well as interactive analysis (16). We used the Expression module in the GEPIA database to analyze the differential expressions of β3GNT9 in normal tissues and in glioblastoma tissues. The median of β3GNT9 expression was taken as the reference point to divide the tissues into high and low expression groups, and the “survival” package in R was employed to draw the Kaplan-Meier survival curves to analyze the impact of β3GNT9 on the prognosis of glioblastoma patients. Later, based on the survival data of glioblastoma patients in the TCGA-GBM database combined with their β3GNT9 expression, univariate and multivariate COX regression analyses were conducted to figure out the effects of clinical factors and β3GNT9 on patient survival.

Gene Set Enrichment Analysis (GSEA) is a bioinformatics analysis method used to determine whether there is such a gene set that consists of genes with biological function, chromosomal location, or regulation (17). GSEA can be applied to predict the potential working mechanism of the studied genes. In this study, gene enrichment analysis was performed using the GSEA 4.1.0 software. We included gene sets with a P value < 0.05 and false discovery rate (FDR) < 0.05 as significant enriched gene sets.

CIBERSORT (http://cibersort.stanford.edu/), a method for assessing the composition of tissue cells based on gene expression profiles, is a gene expression-based deconvolution algorithm that evaluates the expression changes of a set of genes relative to other genes in a sample (18). We assessed the immune response of 22 tumor immune-infiltrating cells in glioblastomas using CBIERSORT to identify the relationship between survival and molecular composition. We uploaded the data-standardized glioblastoma mRNA dataset to CIBERSORT, and obtained 22 immune tumor immune-infiltrating cells data when setting permutation displacement 1000 times. Glioblastoma patients were divided into screening gene high expression groups and low expression groups, and then the differences of immune cells between β3GNT9 high expression group and low expression group were identified and then visualized.

TIMER (https://cistrome.shinyapps.io/timer/) database, a web server for comprehensive analysis of tumor-infiltrating immune cells, covers 32 tumors with a total of 10897 samples (19). The TIMER database includes six modules: gene module, survival module, mutation module, somatic copy number module, and gene expression module. We used the gene module to analyze the correlation between β3GNT9 and the infiltration levels of six immune cells (i.e., B cells, CD8+T cells, CD4+T cells, macrophages, neutrophils and dendritic cells) in glioblastomas.

The human astrocyte cell line NHA and the human glioblastoma cell lines A172, U87MG, U251 were purchased from Procell Life Science & Technology Company (Wuhan, China) and cultured in the Dulbecco’s Modified Eagle Medium (DMEM; Gibco, Grand Island, NY, USA) containing 10% fetal bovine serum (FBS; Gibco, Grand Island, NY, USA) and 1% penicillin/streptomycin (Cyagen Biosciences, Guangzhou, China) at 37°C in a cell incubator with an atmosphere of 5% CO₂ and humidity saturation (20). The culture medium was replaced once every 2 days. Cells in the log growth phase were taken for the experiment.

When the cells reached 70% confluence, β3GNT9-shRNA was transfected into U251 cells informed by the procedure instruction of the transfection reagent Lipofectamine TM2000 (Invitrogen, Carlsbad, CA, USA) (21). The β3GNT9 knockdown shRNA plasmid was purchased from Shanghai Genechem Co., Ltd. (Shanghai, China).

Four logarithmically growing cell lines were taken and RNA was extracted according to the Trizol instructions (Toyobo (Shanghai) Biotech, Shanghai, China). Reverse transcription reactions were performed according to the instructions by steps on the test kit (Toyobo (Shanghai) Biotech, Shanghai, China) and the quantitative real-time PCR according to the manufacturer’s instructions (Roche, Basel Switzerland). The primer sequences of β3GNT9 (forward: 5′-GCTCAAGGAGATCCACTTTCT-3′, reverse: 5′-CAGGAGATTTCCCACGTTCA-3′) and the primer sequences of the internal reference GAPDH (forward: 5′-CTGGGCTACACTGAGCACC-3′, reverse: 5′-AAGTGGTCGTTGAGGGCAATG-3′) were used. The expression levels of β3GNT9 mRNA in the cells were analyzed through the 2-ΔΔCt assay (22).

At 6 hours post-transfection, the cells were digested, centrifuged, and re-suspended in FBS-free cultures. The cell concentration was adjusted to 2×10⁵ cells/ml. After cell mounting, the cell layer was then scratched with a 200 μl of sterile pipette tip. After removing cell culture medium and suspended cells and cell debris, each well was added into serum-free medium and cultured in the incubator for 24 and 48 hours when cell migration areas were photographed. The scratch assays were evaluated for differences in terms of the healing ability.

Matrigel (Corning, NY, USA) was melted at 4°C overnight. Then 100 μl of diluted matrix gel was added to the chamber. After that, 200 μl of serum-free medium was added to the upper chamber, and 500 μl of DMEM containing 10% FBS was added to the lower chamber. The collected cells with a density of 2×10⁵ cells/ml were planted in the upper chamber of culture plate and cultured in the incubator for an additional 48 hours. Subsequently, the invasion chamber was removed, and the cells on the polycarbonate membranes were fixed with 4% paraformaldehyde (Beyotime, Shanghai, China), and then stained with 0.1% crystal violet (Beyotime, Shanghai, China). Three fields were selected at random, and the invading cells were counted under the microscope.

We utilized R programming language (version 3.6.3; https://www.r-project.org) and GraphPad Prism 7 (GraphPad Software, San Diego, CA, USA) for statistical analysis. The research data are presented as the means ± standard error of the means (SEM) or standard deviation (SD). A P value less than 0.05 was considered significant. The comparison of the experimental data was evaluated using the Student’s t test for two groups and one-way analysis of variance (ANOVA) for multiple groups. Differences in β3GNT9 expression between normal tissues and tumor tissues were analyzed using Wilcoxon signed-rank tests and one-way analysis of variance. The association between β3GNT9 expression and clinicopathological features was analyzed using the Kruskal-Wallis test, Dunn’s test, and Wilcoxon rank-sum test. Moreover, the evaluation of the prognostic significance of β3GNT9 expression was based on Kaplan-Meier method and Cox regression analysis. Independent prognostic factors for glioblastoma were identified via the R package. Univariate and multivariate Cox proportional hazard regression analyses were conducted to exclude β3GNT9 and four clinicodemographic characteristics (age, sex, radiation and pharmaceutical) from overall survival (OS) analysis.

Differential genes were screened using the limma package under the threshold of log2FC absolute value > 1 and a P value < 0.05. 18689 differential genes were obtained with their heatmap shown in Figure 1. These differential genes displayed different levels of expression only in tumor and normal tissues. Prognostic genes were obtained when further combining the clinical survival data with differential genes. We searched for those genes that have high clinical correlation but had not been studied previously, and finally we identified β3GNT9.

Difference analysis was conducted through combining TCGA-GBM and GTEx data. The results showed that the expression of β3GNT9 in glioblastoma tissues was significantly higher than that in normal tissues (P<0.05, Figure 2A). Further exploration of the relationship between β3GNT9 and clinically relevant indicators of glioblastoma suggested that patients older than 65 years of age had higher β3GNT9 expression than patients of or under 65 years of age (P<0.05, Figure 2B).

To explore the effect of β3GNT9 on glioblastoma patients’ survival, we evaluated the prognostic value of β3GNT9 using the “survival” package in R (Figure 3A). In addition, we also evaluated the prognostic value of β3GNT9 using the online database GEPIA. These Kaplan–Meier survival curves show high expression of β3GNT9 was associated with shorter OS (P<0.05, Figures 3B, C). In total, higher expression of β3GNT9 tends to have a worse prognosis compared to those with low expression.

According to the univariate COX regression analysis (Figure 4A), age, radiotherapy, chemotherapy, and β3GNT9 were significantly associated with OS in glioblastoma patients (P<0.05). Incorporating the above variables and gender, the multivariate COX regression analysis showed that radiotherapy and β3GNT9 were independent prognostic indicators for patients. Radiotherapy is a protective prognostic indicator for glioblastoma (P<0.05), whereas high expression of β3GNT9 is an independent prognostic indicator for glioblastoma (Figure 4B).

Genes were subjected to GSEA analysis, and gene sets with P<0.05 and false discovery rate FDR<0.05 were considered as significantly enriched gene sets. The results showed ten significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in the β3GNT9 high expression phenotype, including cardiac muscle contraction, extracellular matrix receptor interaction, glycosaminoglycan biosynthesis chondroitin sulfate, glycosaminoglycan biosynthesis heparan sulfate, glycosaminoglycan biosynthesis keratan sulfate, glycosaminoglycan degradation, Huntington’s disease, oxidative phosphorylation, Parkinson’s disease, riboflavin metabolism (Figure 5A). Five significant GO pathways in the β3GNT9 high expression phenotype included sulfuric ester hydrolase activity, acylglycerol-O-acyltransferase activity, phosphatidylinositol 3-phosphate bindin, endoplasmic reticulum lumen, Golgi cisterna membrane (Figure 5B). Five significant GO pathways in the β3GNT9 low expression phenotype included regulation of ubiquitin protein ligase activity, ubiquitin-like protein conjugating enzyme activity, exit from mitosis, neural nucleus development, substantia nigra development (Figure 5B).

We investigated whether the expression of β3GNT9 was associated with immune infiltration in glioblastoma using the online database TIMER. The results revealed that β3GNT9 expression was significantly and positively associated with the infiltration levels of CD4+T cells, macrophages, dendritic cells (P<0.05), but negatively correlated with tumor purity, B cells, CD8+T cells, and neutrophils (P<0.05, Figure 6).

Glioblastomas were divided into β3GNT9 high expression half and low expression half, and the differences between the 22 immune cell expression levels were analyzed using the CBIERSORT algorithm. It was found that T cells follicular helper, activated NK cells, monocytes, M2 macrophages and eosinophils were the main immune cells affected by β3GNT9 expression. Among them and compared with the low expression group, helper follicular T cells (P=0.033), activated NK cells (P=0.031), mononuclear leukocytes (P=0.023), M2 macrophages (P=0.048), and eosinophils (P=0.014) in the high β3GNT9 expression group were significantly reduced (Figure 7). Afterwards, the heat maps of 22 immune cells were drawn as shown in Figure 8. Based on the immune infiltration analysis, we speculated that the high expression of β3GNT9 may lead to the change of activated NK cells and macrophages in tumor microenvironment, which in turn leads to poor prognosis of patients.

The analysis of β3GNT9 expression in various tumors using the online database GEPIA revealed significant high β3GNT9 expression in esophageal cancer, glioblastoma, head and neck squamous cell carcinoma, low-grade glioma, pancreatic cancer, and thymic carcinoma (P<0.05), and significant low expression in chromophobe renal cell carcinoma, ovarian serous cystadenocarcinoma, prostate cancer, and endometrial cancer (P<0.05, Figure 9).

We evaluated the expression of β3GNT9 in a variety of glioblastoma cell lines and normal astrocytes and found that β3GNT9 expression was higher in glioblastoma cell lines than in NHA cells (*P<0.05, **P<0.01, Figure 10A). The mRNA levels declined after targeted β3GNT9 silencing (**P<0.01, Figure 10B).

We performed cell scratch and invasion assays to examine whether the downregulation of β3GNT9 affected the invasion or migration ability of glioma cells. Cell migration was assessed using a scratch wound healing assay and the extent of cell migration to the scratch area was measured. The wounds in the wells of the control transfected cells healed rapidly, with almost no gaps left after 48 hours. However, the wells of β3GNT9-shRNA transfected cells showed a much lower healing capacity than the control wells (P<0.05, Figures 11A, B). Furthermore, the number of cells transfected by invading β3GNT9-shRNA was significantly lower than that of the control transfected cells (P<0.05, Figures 11C, D). These suggested that β3GNT9 silencing inhibited the migration and invasion capacity of glioblastoma cells.