Lysyl oxidase-like 1 predicts the prognosis of patients with primary glioblastoma and promotes tumor invasion via EMT pathway

Prognostic variants identified via the genomic sequencing of LOXs in GBM

cBioPortal integrates gene expression, alteration information, and clinical data from multiple databases (Gao et al., 2013). We integrated 429 GBM samples with complete mutation and copy number alteration (CNA) data (Table S1: Brain Tumor PDXs, Mayo Clinic (Vaubel et al., 2020); Glioblastoma, Columbia (Zhao et al., 2019); GBM, TCGA, PanCancer Atlas) to analyze the ratio and sites of each LOXs genetic alteration. The potential prognostic value of LOXs genetic alterations was also investigated with Kaplan‒Meier survival analysis.

TCGA and CGGA databases

The GBM cohorts included both isocitrate dehydrogenase (IDH)1-wild-type (WT) and IDH1-mutant patients from GlioVis (TCGA) (Bowman et al., 2017) and the CGGA website (Zhao et al., 2021b). We explored the potential diagnostic and prognostic value of LOXs via differential expression and survival analyses. After constructing the nomogram and PRSM, we further included another TCGA GBM cohort to test the PRSM.

Biological function analysis of LOXs

A total of 250 genes related to LOXs (50 of each member) in GBM based on the Pearson correlation coefficient were obtained from GEPIA (Tang et al., 2019). Based on the integrated groups of LOXs-related genes, a protein‒protein interaction network was constructed (Table S2) via STRING (von Mering et al., 2005). Further GO and KEGG biological function analyses were subsequently conducted, and the results were visualized with Cytoscape (version 3.8.2).

Survival analysis, model construction, and evaluation

The prognostic role of each LOXs family member in GBM was analyzed by Kaplan–Meier analysis in GBM cohorts (TCGA training cohort, n = 330; Table S3; CGGA validation cohort, n = 193; Table S4; and external TCGA validation cohort, n = 266; Table S5). Subsequently, the prognostic value of LOXs was explored through Cox survival analyses, which revealed clinicopathological variables such as age, sex, tumor classifications (including IDH mutation status and GBM subtypes), and treatment modalities (including chemotherapy and radiotherapy). Based on these discerned prognostic factors, a nomogram and the PRSM were systematically constructed, providing a robust framework for appraising the 1- and 2-year prognostic outcomes for GBM patients. The predictive risk score formula was as follows: Y = 0.115*Loxl1 + 0.362*Age-0.739*IDH1-1.010*Chemo-1.824*Radio. After excluding GBM patients with IDH1 mutations according to the 2021 WHO CNS tumor classification (Louis et al., 2021), the training TCGA (IDH1-WT, n = 309), CGGA (IDH1-WT, n = 161), and external TCGA (IDH1-WT, n = 246) GBM datasets were used to further assess the model. Time-dependent receiver operating characteristic (t-ROC) curves were used to analyze the model prediction accuracy. The net clinical benefits were judged by calibration and decision curve analyses. The entire spectrum of the data was subjected to rigorous analysis and subsequent visualization, a process meticulously executed through the use of R (version 4.0.3).

Single-cell RNA sequencing data processing

The GSE182109 datasheet (Abdelfattah et al., 2022) from the GEO database contained low-grade glioma (LGG) and newly diagnosed GBM (ndGBM, primary GBM) groups. The quality control metrics are listed below, and the following exclusion criteria were as follows: genes detected in <3 cells, <300 genes detected in cells, a proportion of mitochondria >10%, and a percentage of hemoglobin genes >3%. A total of 3,000 highly variable feature genes were marked by the function “FindVariableFeatures” from the R package “Seurat”. The quality-controlled data were scaled by the function “ScaleData”, and the data dimensionality was reduced by the function “RunPCA”, according to highly variable genes. Finally, the cells were clustered by the functions “FindNeighbors” and “FindClusters” (Dim = 20, Solution = 0.4). T-distributed stochastic neighbor embedding (t-SNE) was used to further reduce the data dimensions and visualize the distribution on the two-dimensional plot. As described in a previous study of GSE182109, marker genes for cell type annotation were obtained and are listed in Table S6 (Abdelfattah et al., 2022). We visualized the expression pattern of each marker gene by the function “DotPlot” and annotated the cell types.

Functional gene enrichment analyses based on scRNA-seq

To identify the enriched functions between the LGG and GBM subgroups, gene set enrichment analysis (GSEA) was conducted with the Molecular Signatures Database (MSigDB) (c2.cp.kegg.v2023.1.Hs.symbols.gmt) via the R package (Liberzon et al., 2015). Gene set variation analysis (GSVA) was also performed to calculate the pathway scores of each cell using the R package “GSVA”. The correlation between Loxl1 and enriched biological processes in GBM could be further explored.

Analysis of tumor-infiltrating immune cells

The correlations between Loxl1 and tumor-infiltrating immune cells (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells) in GBM were explored and visualized with TIMER (https://cistrome.shinyapps.io/timer/) (Li et al., 2017).

Cell culture and construction of the Loxl1 knockdown cell line

The U87 MG human glioblastoma cell line was acquired from the American Type Culture Collection (Virginia, USA) and cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, Billings, MT, USA) supplemented with 10% fetal bovine serum (FBS, F8318-500ML) in a 37 °C and 5% CO2 humidified incubator.

Loxl1 knockdown lentiviral vectors containing Loxl1-specific small hairpin RNA (shLoxl1-1: CCTGGGAACTACATCCTCA, shLoxl1-2: GCATTAAAGCAGCGTATC) and nontargeting shRNA (shNT: CTCGCTTGGGCGAGAGTAA) were purchased from GeneChem (Shanghai, China). In accordance with the manufacturer’s instructions, U87 MG cells were cultured at a density of 10 × 10⁴ cells/well in a six-well plate before transfection. The medium was changed to serum-free DMEM. Then, shNT or shLoxl1-1&2 lentivirus was separately added to the wells at an MOI of 1.5. After 6 h, the lentivirus medium was discarded, and the medium was supplemented with puromycin (1.5 μg/ml). The drug selection medium was replaced with complete medium 24 h later, and Loxl1-knockdown cell lines were constructed for further experiments.

Cell viability and proliferation assay

A cell counting kit-8 (CCK-8) assay was used to detect the effect of Loxl1 knockdown on the viability and proliferation of U87 MG cells. A total of 2,500 stable knockdown cells per well were seeded in a 96-well plate. We measured the absorbance of each well every 24 h for 96 h. Ten microliters of CCK-8 reagent (521942; Biosharp, Heifei, China) was added to each well, and the plate was placed in an incubator at 37 °C for 1 h. The absorbance of the different groups at 450 and 670 nm was read with a microplate reader. Cell viability and proliferation curves were generated and analyzed with GraphPad Prism and SPSS. Multiple comparisons were conducted between the shNT group and the sh-Loxl1-1 & 2 group at each time point.

Western blot analysis

Transfected U87 MG cells were cultured in six-well plates and harvested when they reached >80% confluence. Each well was washed with PBS twice and lysed with RIPA buffer (P0013B; Beyotime, Shanghai, China) supplemented with PMSF for 15 min on ice. Then, we scratched off the cells with scrapers and centrifuged them for 20 min at 4 °C. The supernatant was carefully collected, quantified with a bicinchoninic acid (BCA) kit (P0011; Beyotime, Shanghai, China) and diluted to 1 mg/ml with RIPA buffer. Then, the loading buffer was mixed with the supernatant, and the mixture was denatured at 100 °C for 10 min. The denatured samples were separated by sodium dodecyl sulfate‒polyacrylamide gel electrophoresis (120 V, 60 min) and transferred to a nitrocellulose filter membrane (Millipore, Burlington, MA, USA), whose weight was based on the protein molecular weight. After blocking with 5% BSA at room temperature for 1 h, the membrane was incubated with primary antibodies (anti-Loxl1 (sc-166632, Santa Cruz, CA, USA), anti-N-cadherin (#13116; Cell Signaling Technology, Danvers, MA, USA), anti-E-cadherin (#3195; Cell Signaling Technology, Danvers, MA, USA), anti-Vimentin (#5741; Cell Signaling Technology, Danvers, MA, USA), anti-Snai1 (13099-1-AP; Proteintech, Rosemont, IL, USA), anti-GAPDH (60004-1-Ig; Proteintech, Rosemont, IL, USA), and anti-β-actin (20536-1-AP; Proteintech, Rosemont, IL, USA)) on a shaker at 4 °C overnight, followed by incubation with secondary antibodies (anti-mouse, SA00001-1 & anti-rabbit, SA00001-2; Proteintech, Rosemont, IL, USA) for 1 h at room temperature. Enhanced chemiluminescence reagent (P0018S; Beyotime, Shanghai, China) was added to visualize the bands with a ChemiDoc MP (Bio-Rad, Hercules, CA, USA). The gel bands were quantified by ImageJ (NIH, Bethesda, MD, USA), and the GAPDH/β-actin ratio was used as an internal control. The target protein/internal control in the shNT group was compared with the ratios in the Loxl1 knockdown groups.

Wound healing assay

U87 MG cells were subjected to different groups in six-well plates and reached >90% confluence. Then, a 10-μl pipette tip was used to scratch the bottom, and the medium was changed to serum-free DMEM. The wound gaps were photographed at 0 and 24 h at the same location. The wound gaps were measured by ImageJ software, and the ratio of the closure gap was analyzed between the shNT group and the sh-Loxl1-1&2 groups to detect differences in invasion ability.

Transwell invasion assay

U87 MG cells from different groups were counted and separately cultured in Transwell chambers (8.0 µm pore) at 1.5 × 10⁴ cells per well. The chamber was placed in a 24-well plate, and its membrane was coated with Matrigel (100 µg/cm²) (BD Bios, Franklin Lakes, NJ, USA). U87 MG cells were cultured in the upper compartment, and 10% FBS-containing medium was added to the lower compartment. After 24 h, the chambers were washed twice with PBS, fixed in 4% paraformaldehyde for 15 min and washed. Then, the cells were stained with 0.5% crystal violet for 15 min and washed twice. Photography was conducted with a microscope (ECLIPSE Ni-U, Nikon, Tokyo, Japan), and the number of cells per field was counted for the shNT group and the sh-Loxl1-1 & 2 groups.

Statistical methods

Statistical distinctions between the two cohorts were assessed through Student’s unpaired two-sided t test. Within the arms of multiple groups, one-way analysis of variance (ANOVA) was used, accompanied by subsequent post hoc least significant difference assessments for multiple pairwise comparisons. Survival evaluations, encompassing both Kaplan–Meier survival analyses and univariate/multivariate Cox regression analyses, were performed with SPSS (Chicago, IL, USA) and R software (R Core Team, 2020). A significance threshold of P < 0.05 was consistently upheld as the standard for statistical significance across all analytical appraisals.

Genomic sequencing of LOXs and prognostic variant identification

To delineate the correlation between genetic aberrations in LOXs and the prognostic landscape of GBM patients, we analyzed a cohort comprising 429 specimens obtained from three distinct genomic profiling studies of GBM patients. The LOXs family, which includes all five constituents, exhibited varying degrees of genetic perturbations within the GBM patient pool. The cumulative frequency of such genetic alterations reached 2.6%, accounting for 11 of the 429 samples under scrutiny. As shown in Figs. 1A and 1B, 0.4% of GBM patients presented missense mutations in Lox. This mutation often results in Lox dysfunction, with 269 and 274 amino acid changes (aa), and affects the phosphorylation and O-linked glycosylation of Lox. Missense mutations and inframe mutations were detected in Loxl1 with the same rate of genetic alteration as Lox. This missense mutation led to 421 amino acid changes in Loxl1. Loxl2 was the gene with the most genetic alterations (0.8%), and its mutation may contribute to changes in gene phosphorylation. A total of 0.4% of GBM samples presented missense mutations, which caused 545 amino acid changes in Loxl3. A total of 0.69% of GBM samples showed a missense mutation and amplification in Loxl4, of which 583 amino acid changes were detected in Loxl4.

To further determine the prognostic value of LOXs with genetic alterations in GBM, we performed a Kaplan‒Meier analysis of 11 GBM patients with LOXs genetic alterations from three GBM cohorts (Fig. 1C). GBM patients exhibiting genetic modifications pertaining to LOXs displayed notably truncated overall survival intervals compared to those of their counterparts devoid of such alterations (Fig. 1D). The prognostic genetic alteration of LOXs in GBM indicated the important clinical value of LOXs.

LOXs expression analysis and correlated biological functions

We performed differential expression analyses of LOXs in TCGA and CGGA at the transcriptional level, as shown in Figs. 2A and 2B. In contrast to their expression levels in normal tissue counterparts, the expression profiles of all LOXs (Lox, Loxl1, Loxl2, Loxl3, Loxl4) were notably elevated within the domain of GBM tissues, indicating the plausible oncogenic functions of LOXs in GBM. To confirm the functional involvement of LOXs in GBM, we subjected the 250 most relevant genes (50 for each respective member) to comprehensive scrutiny through GO analyses. The processes of extracellular structure organization (biological process, BP) and extracellular matrix structural constituent (molecular function, MF) in the endoplasmic reticulum lumen (cellular component, CC) were correlated with the expression of LOXs in GBM based on the GO and KEGG analyses. Moreover, KEGG analysis revealed that the biological roles of LOXs in GBM were associated with the pathways of protein digestion and absorption (Fig. 2C). Together, these results suggested that LOXs serve as important regulators of ECM protein expression.

LOXs and clinical GBM classifications

GBM is classified into the proneural, classical, and mesenchymal subtypes. The classifications were related to patient survival, of which mesenchymal subtypes displayed the shortest median survival (Verhaak et al., 2010). As shown in Figs 3A, 3B, we found that, except for Loxl1 in the CGGA, all LOXs in GBM were more highly expressed in the mesenchymal subtypes than in the classical and proneural subtypes, indicating that LOXs could serve as diagnostic indicators for predicting mesenchymal subtypes in GBM patients.

Another important indicator widely used in GBM diagnosis and treatment is IDH mutation (Yan et al., 2009). It is associated with better outcomes in IDH1-mutant patients than in IDH1-WT patients. We next explored the correlation between IDH mutation status and LOXs expression. Interestingly, we found that all LOXs in GBM were more highly expressed in the IDH1-WT subgroup than in the IDH1-mutant subgroup (Figs. 3C, 3D). The correlation between high LOX expression and IDH1-WT expression was confirmed both in the TCGA and CGGA cohorts, suggesting that LOXs could act as promising markers in GBM.

Independent prognostic variables and nomogram construction

To determine the effective prognostic factors of LOXs in GBM patients, two separate databases were used. Overall survival of GBM patients with high Lox, Loxl2, Loxl3, and Loxl4 expression was similar to that of patients with low Loxl2, Loxl3, and Loxl4 expression. Moreover, patients with higher Loxl1 levels might have a shorter survival time, which might be a prognostic indicator for clinical diagnosis and treatment (Fig. 4A).

The elucidation and validation of prognostic indicators in GBM encompassed an intricate interplay between patient-specific clinical attributes and Loxl1 expression patterns. Parameters spanning age, sex, tumor classifications encompassing IDH mutation status and distinct GBM subtypes, as well as the therapeutic modalities of chemotherapy and radiotherapy, were critically engaged in this multifaceted endeavor. The results of univariate Cox survival analyses revealed that younger age, IDH mutation status, chemotherapy, radiotherapy, and low Loxl1 expression were protective variables, which translated into an extended survival period for patients with GBM (Fig. 4B).

To ascertain the clinical predictive efficacy inherent to Loxl1 expression, we integrated the aforementioned clinical prognosticators into a multivariate Cox survival analysis framework. Through the deliberate exclusion of sex and GBM subtypes as confounding elements, we proceeded to gauge the adjusted hazard ratios for these prognostic determinants along with their corresponding 95% confidence intervals (Fig. 4C). To ascertain the clinical predictive efficacy inherent to Loxl1 expression, we integrated the aforementioned clinical prognosticators into a multivariate Cox survival analysis framework. Through the deliberate exclusion of sex and GBM subtypes as confounding elements, we proceeded to gauge the adjusted hazard ratios for these prognostic determinants along with their corresponding 95% confidence intervals (Fig. 4D). The nomogram included age, IDH mutation status, chemotherapy, radiotherapy, and Loxl1 expression as pivotal parameters contributing to its predictive power.

Verification and evaluation of the Loxl1-based nomogram in GBM cohorts

Within the training cohort of GBM patients from the TCGA, the nomogram and PRSM were constructed based on the combination of Loxl1 and clinicopathological parameters, including age, IDH mutation status, chemotherapy, and radiotherapy (Fig. 5A). The computation of individual patient risk scores facilitated an insightful determination, where elevated risk scores exhibited a noteworthy correlation with diminished survival intervals (HR = 1.82–2.97, P < 0.001) (Fig. 5D). To demonstrate the generalizability of the nomogram model, the CGGA GBM cohort and another different TCGA GBM cohort were utilized to confirm the model’s predictive efficiency. The risk scores of both validation sets were also calculated and are presented in Figs. 5B, 5C. Kaplan–Meier survival analyses of the CGGA and external TCGA validation cohorts revealed results similar to those of the training TCGA cohort (Figs. 5E, 5F), suggesting that the Loxl1-based PRSM classification functions stably in different populations.

To substantiate the prognostic precision of the model, we calculated the area under the curve (AUC) through a comprehensive analysis of t-ROC curves. Our findings indicate that all cohorts exhibited notably substantial AUC values. Specifically, the AUC values for the Loxl1-based nomogram for predicting 1- and 2-year overall survival durations were 0.760 and 0.792, respectively, in the TCGA training set, 0.684 and 0.730, respectively, in the CGGA validation set, 0.733 and 0.685, respectively, in the external TCGA validation set (Figs. 5G–5I), suggesting that the Loxl1-based nomogram and the PRSM were highly accurate and robust. We also generated a calibration curve of the nomogram in different cohorts. The diagrams show a congruence between the actual and projected outcomes at the 1- and 2-year time points concerning the survival duration of GBM patients, which was validated across both the training cohort (Fig. 5J) and the independent validation cohorts (Figs. 5K, 5L). The clinical applicability of the established Loxl1-based nomogram was evaluated by decision curve analysis (DCA). More net benefits could be achieved by the nomogram for almost all threshold probabilities when compared with other prognostic variables in the training (Fig. 5M) and validation cohorts (Figs. 5N, 5O). These results suggested that the Loxl1-based PRSM had a better ability to predict survival than other variables in primary GBM patients.

Cox regression analysis of the training cohort notably revealed a statistically significant association between the presence of IDH mutations and overall survival duration among GBM patients, including those with IDH1- mutations and IDH1-WT. However, the 2021 WHO classification further defined patients with IDH1-WT as GBM patients (Louis et al., 2021). We further excluded IDH-mutant patients from all cohorts and utilized the Loxl1-based nomogram without the IDH mutation parameter. The risk scores were calculated again (Figs. S1A–S1C), and the high-score groups had shorter survival times (Figs. S1D–S1F) in both the training and validation sets. Based on the new t-ROC (Figs. S1G–S1I), calibration curve (Figs. S1J–S1L), and DCA (Figs. S1M–S1O) analyses, we found that the prediction accuracy and net benefits of the model were still stable regardless of whether GBM was characterized by IDH1 mutation.

Loxl1-related enrichment pathways and tumor immune infiltrates in GBM

As shown in Fig. 6A, we processed the GSE182109 scRNA-seq data, and the distribution of all the raw cell samples is shown. The expression of Loxl1 in glioma cells of GBM tissues could help us understand its functional state. The stringent quality control metrics previously described were used to screen 29,402 high-quality cells. These cells were scaled and clustered into 13 clusters (Fig. 6B). Based on the annotated marker genes (Table S6), we classified the cells into myeloid cells (C0, C3), glioma cells (C1, C2, C5, C7, C8, C10, C11), T cells (C4), pericytes (C9), endothelial cells (C12), and oligodendrocytes (C6) (Fig. 6C). The glioma cells were separated and are shown in Fig. 6D according to tumor grade (LGG or GBM). Combined with visualization of the Loxl1 expression pattern (Fig. 6E), we found that Loxl1 was highly expressed in the GBM groups.

To explore the enriched biological functions in GBM, GSEA was employed, and the results revealed that ECM receptor interaction (P < 0.05) was positively associated with GBM (Fig. 6F). Furthermore, GSVA was performed to calculate the pathway scores of each cell. The expression of Loxl1 was also found to be positively correlated with the interaction of ECM receptors (Fig. 6G). Previous results of GO and KEGG analyses of GBM bulk sequencing data support the above results, suggesting that Loxl1 might participate in the progression of tumor invasion. We also explored the correlation between Loxl1 and various immune cells. Based on the CNA levels, we found that arm-level gain of Loxl1 was significantly associated with CD8+ T cells and dendritic cells infiltration, while its arm-level deletion was correlated to CD4+ T cells and neutrophils infiltration (Fig. 6H). At the mRNA level, Loxl1 expression was positively related to dendritic cells and negatively correlated with B cells and neutrophils (Fig. 6I), which implies that Loxl1 has promising potential for assessing the efficacy of GBM immunotherapy.

Decreased Loxl1 inhibited the invasion of GBM cells

Loxl1-knockdown U87 MG cells were generated to further confirm the invasive role of Loxl1 in GBM (Figs. 7A, 7B). CCK-8 assays revealed that the viability and proliferation of cells in the sh-Loxl1 group were significantly lower than those in the control group (Fig. 7C). The EMT pathway plays vital roles in tumor cell invasion and migration. As previously reported, the expression of E-cadherin, N-cadherin, Vimentin and Snai1 is associated with the EMT pathway. First, we found that N-cadherin was decreased (Figs. 7D, 7E) and that E-cadherin was increased (Figs. 7F, 7G) in both sh-Loxl1 groups. Decreases in Vimentin and Snai1 were also observed by western blot analyses following the downregulation of Loxl1 (Figs. 7H, 7I; Figs. 7J, 7K). To determine the effect of Loxl1 on the motility and invasion of U87 MG cells, wound healing and Transwell assays were performed. Wound closure was lower in the sh-Loxl1 group than in the control group (Figs. 7L, 7M). The number of invasive glioma cells in the Transwell assay also decreased significantly after shLoxl1 knockdown (Figs. 7N, 7O), confirming that Loxl1 contributes to the invasion of glioma cells. These data confirmed that Loxl1 plays an invasive role in tumor progression by regulating the EMT pathway, which might be a promising therapeutic strategy for GBM patients.