An integrated bioinformatic investigation of focal adhesion-related genes in glioma followed by preliminary validation of COL1A2 in tumorigenesis

Molecular classification based on FARGs in TCGA cohort

First, we show the workflow of this study in Figure 1. According to the flow chart, 192 FARGs were obtained from the intersection of TCGA, CGGA, and GSE16011 cohorts (Figure 2A). A novel molecular classification was identified by using an unsupervised clustering algorithm based on the expression profiles of the 192 FARGs extracted from the TCGA cohort. According to the area under the CDF curves, probably approximately correct algorithm and correlations between clusters, we determined the optimal number of clusters is two (k=2) (Figure 2B–2E). Principal components analysis (PCA) can effectively distinguish glioma patients with the FARG1 subtype and FARG2 subtype (Figure 2F).

Figure 1. Flow diagram of this study.

Basic analysis between different subtypes in the TCGA cohort

The OS and progression-free survival (PFS) of glioma patients, low-grade glioma (LGG) patients, and GBM patients with the FARG1 subtype were significantly longer than those with the FARG2 subtype (Figure 2G, 2H). Except for gender, there are dramatic differences in the distribution of the other four clinicopathological characteristics between FARG1 and FARG2 subtypes. In addition, most FARGs are differentially expressed between the two subtypes (Figure 2I). We simultaneously performed GSVA in TCGA, CGGA, and GSE16011 cohorts and identified 167 functional pathways with different relative activity between the two subtypes (Figure 2J and Supplementary Table 3), of which 65 functional pathways were more active in FARG1 subtype, such inositol phosphate metabolism, glycerol lipid metabolism and lysine degradation, and the other 102 functional pathways were more active in FARG2 subtypes, such as cell lung cancer, pancreatic cancer, and primary immunodeficiency.

Figure 2. Novel molecular classification based on FARGs in TCGA cohort. (A) The intersection of 199 FARGs between TCGA, CGGA, and GSE16011 cohorts. (B–D) Consensus clustering matrix based on 192 FARGs for k=2, k=3, and k=4. (E) The optimal number of consensus clustering matrices is determined by the cumulative distribution function of the unsupervised consensus clustering algorithm. (F) Obvious differences in the transcriptomes between FARG1 and FARG2 subtypes were analyzed by principal component analysis (PCA). (G, H) K-M survival curves of glioma patients. (I) The distributions of five clinicopathological characteristics and 199 FARGs between FARG1 and FARG2 subtypes. (J) GSVA between FARG1 and FARG2 subtypes. Red and blue represent the relative activation and inhibition of the pathways, respectively. *P<0.05, **P<0.01, ***P<0.001.

Immune-related analysis between different subtypes in TCGA cohort

Given our above finding, we further investigated whether different subtypes have different immunological characteristics. The results revealed that there were remarkable differences in most immune signatures (except for dendritic cells and mast cells) and all four TME-related scores among the two subtypes (Figure 3A, 3C–3F). Compared with the FARG1 subtype, the glioma tissue of the FARG2 subtype contains more immune cells and stromal cells. We also estimated the content of immune cells and found that the content of almost all immune cells was different between the two subtypes (Figure 3B and Supplementary Table 4). Moreover, we also found that except for TNFSF9, other immune checkpoints (ICPs) were differentially expressed between the two subtypes of FARG1 and FARG2 (Figure 3G). Considering the effect of gene mutations on tumorigenesis and progression, we found that the TMB was dramatically different between the two subtypes (Figure 3H).

Figure 3. Immune-related analysis between FARG1 and FARG2 subtypes in the TCGA cohort. (A) Difference analysis of 29 immune signatures and four types of TME-related scores between FARG1 and FARG2 subtypes. (B) The difference in the content of immune cells was calculated by CIBERSORT, MCPCOUNTER, and QUANTISEQ algorithms between the two subtypes. (C–F) Difference analysis of TME-related scores between the two subtypes. (G) Difference analysis of ICPs between two subtypes. (H) Difference analysis of TMB between the two subtypes. *P<0.05, **P<0.01, ***P<0.001, ****p < 0.0001.

Construction of a nine-gene scoring system

To establish a scoring system, 175 FARGs were screened from 192 FARGs using uni-Cox regression analysis (Supplementary Table 5). These genes were further screened using the LASSO algorithm, and then 17 FARGs were obtained (Figure 4A, 4B). After that, multi-Cox regression analysis in both directions was performed and yielded nine FARGs (CDC42, COL1A2, ELK1, ERBB2, FLNC, PDGFA, PTEN, SPP1, and THBS4) (Figure 4C). The nine genes were used to construct a scoring system according to the following equation: risk score = 0.298 × (CDC42 expression) + 0.103 × (COL1A2 expression) + 0.408 × (ELK1 expression) + 0.237 × (ERBB2 expression) + 0.111 × (FLNC expression) + 0.179 × (FDGFA expression) - 0.577 × (PTEN expression) + 0.259 × (SPP1 expression) + 0.138 × (THBS4 expression).

Figure 4. Construction of scoring system and analysis of its prognostic value in the TCGA cohort. (A) LASSO coefficient profiles of 192 FARGs. (B) Selection of the penalty parameter (λ) in the LASSO model via 1,000 cross-validations. The vertical dotted line passing through the red dot is drawn at the optimal value. (C) Nine FARGs screened by multi-Cox regression analysis were used to construct the scoring system. (D, E) The K-M survival curves showed that the OS and PFS of glioma patients in the low-risk group were longer than those in the high-risk group. (F, G) ROC curves show that risk score can effectively predict the 1,3,5-year OS and PFS for glioma patients. (H) The distributions of five clinicopathological characteristics and nine FARGs between high and low-risk groups. (I, J) Uni- and multi-Cox regression analysis of risk score and five clinicopathological characteristics. *P<0.05, **P<0.01, ***P<0.001.

Prognostic analysis of scoring system

Since there are significant differences in the prognosis outcomes between GBM and LGG patients, IDH mutant and wild-type glioma patients, we should not only explore the prognostic value of the scoring system in pan-gliomas but also explore its prognostic value among different grades and different IDH mutation types. If the risk score of glioma patients was higher than the median value, they were classified as a high-risk group, otherwise, they were classified as a low-risk group. K–M survival curves showed that in the three cohorts, the OS of patients with pan-gliomas (Figure 4D and Supplementary Figure 1A, 1B), LGG (Supplementary Figure 2A, 2E, 2I), GBM (Supplementary Figure 2B, 2F, 2J), IDH mutation (Supplementary Figures 3A, 4A, 4G) and IDH wildtype (Supplementary Figures 3B, 4B, 4H) in the low-risk group tended to be longer than that of patients in the high-risk group. Moreover, the risk score can effectively predict the 1, 3, and 5-year OS rate of patients with pan-gliomas (Figure 4F and Supplementary Figure 1C, 1D), LGG (Supplementary Figure 2C, 2G, 2K), GBM (Supplementary Figure 2D, 2H, 2L), IDH mutation (Supplementary Figures 3C, 4C, 4I) and IDH wildtype (Supplementary Figures 3D, 4D, 4J) in the three cohorts (Supplementary Table 6). In addition, in the TCGA cohort, patients with pan-gliomas, IDH mutant and IDH wildtype in the low-risk group also tended to have longer PFS than patients in the high-risk group (Figure 4E and Supplementary Figure 3G, 3H), and the risk score can also effectively predict the 1, 3, 5-year PFS rate of patients with pan-gliomas, IDH mutant and IDH wildtype (Figure 4G, and Supplementary Figure 3I, 3J and Supplementary Table 6). Next, we explored the differential distribution of five clinicopathological characteristics and the differential expression of nine FARGs between the two risk groups in the three cohorts. The results showed that except for gender, there were significant differences in the distribution of the other four clinicopathological characteristics and the expression of nine FARGs between different risk groups (Figure 4H and Supplementary Figure 1E, 1G). Finally, we explored whether the scoring system can be used as an independent prognostic factor for patients with pan-gliomas, LGG, GBM, IDH mutant, or wildtype. The results showed that the independent prognostic value of risk score was superior to that of age, gender, grade, IDH status, and 1p19q status in TCGA, CGGA, and GSE16011 cohorts (Figure 4I, 4J, and Supplementary Figures 1F, 1H, 3E, 3F, 3K, 3L, 4E, 4F, 4K, 4L).

Genomic variation analysis between high and low-risk groups

Because genomic variation can impact tumor immunity and immune cell infiltration patterns [28], we explored the association between genomic variation and risk score in this study. The top-16 genes with the highest mutation frequency have significant differences between the two risk groups (Figure 5A, 5B). There are also obvious differences in the distribution of mutation and wildtype of six well-known genes (TP53, PTEN, IDH1, EGFR, ATRX, TTN) between high and low-risk groups (Figure 5C). Besides, we found that glioma patients with low-risk scores inclined to have a lower TMB (Figure 5D). Moreover, patients in the low-risk group had a lower frequency of CNV, either amplification or deletion (Figure 5E, 5F).

Figure 5. Comparisons of genomic variations and functional annotations between high and low-risk groups in the TCGA cohort. (A, B) Waterfall plots of the top-16 mutated genes. The genetic alteration types are listed at the bottom of the waterfall plot. The mutation frequencies of genes are listed on the right side of the waterfall Plot. (C) Differential distribution of mutation and wildtype of six well-known genes between the two risk groups. (D) Comparison of TMB between high and low-risk groups. (E) Circular diagram of chromosome amplification and deletion between high and low-risk groups. (F) CNV frequencies (amplification or deletion) were lower in the low-risk group than that in the high-risk score group. (G) GSVA between FARG1 and FARG2 subtypes. Red and blue represent the relative activation and inhibition of the pathways, respectively.

Immune-related analysis between high and low-risk groups

By performing GSVA in TCGA, CGGA, and GSE16011 cohorts, we found that the relative activities of 156 functional pathways were different between the two risk groups: the relative activity of 53 functional pathways increased in the low-risk group, mainly focusing on metabolic-related pathways, such as “butanoate metabolism”, “propanoate metabolism” and “glyoxylate and dicarboxylate metabolism”; the relative activity of other 103 functional pathways increased in the high-risk group, including immune-related pathways (Figure 5G and Supplementary Table 7). Next, we began an in-depth analysis of the relationship between immunologic characteristics and risk scores. We found that the TME-related scores were not only highly correlated with the risk score (Supplementary Figure 5C–5E) but also significantly different between the high and low-risk groups in the three cohorts (Figure 6A and Supplementary Figure 5A, 5B). Then, we also found that the enrichment scores of most immune signatures calculated by ssGSEA were different between the two risk groups in the three cohorts (Figure 6A and Supplementary Figure 5A, 5B). Considering that immune cells play a key role in tumorigenesis, tumor progression, and tumor immunity, we further explored the association between the scoring system and immune cells. The content of most immune cells was not only correlated with the risk score (Supplementary Figure 6A–6C) but also significantly different among different risk groups (Figure 6B and Supplementary Figure 6D, 6E). The results of the analysis between ICPs and risk scores (Figure 6C and Supplementary Figure 7A–7E) prompted us to continue to explore whether the scoring system can effectively predict the immunotherapy response of glioma patients. There was a significant difference of the risk score between the responder and non-responder groups (Figure 6D and Supplementary Figure 8A–8D). The proportion of responders in the high-risk group is much higher than those in the low-risk group, which indicates that glioma patients in the high-risk group are more likely to benefit from immunotherapy (Figure 6E and Supplementary Figure 8B–8E). The subclass mapping analysis indicated that glioma patients in the high-risk group are more sensitive to anti-PD-1 therapy (Bonferroni P = 0.016,0.032, 0.039 in GSE16011, TCGA, and CGGA cohorts, respectively) (Figure 6F and Supplementary Figure 8C–8F). If the Bonferroni corrected P > 0.05, the P-value would not be marked in the subgroup map, but from the color comparison of the corresponding patches in the subgroup map, it can be seen that the P-value of anti-CTLA4 therapy in the high-risk group is lower than that in the low-risk group, indicating that glioma patients in the high-risk group are more likely to respond to anti-CTLA4 therapy (Figure 6F and Supplementary Figure 8C–8F).

Figure 6. Immune-related analysis in the TCGA and IMvigor210 cohorts. (A) Difference analysis of immune signatures and four types of TME-related scores between high and low-risk groups in the TCGA cohort. (B) The difference in infiltration level of immune cells between high and low-risk groups. (C) Differential expression of six well-known ICPs between high and low-risk groups. (D) The difference in risk score between high and low-risk groups was stratified according to response to immunotherapy. (E) The difference in the distribution of responders and non-responders between the two risk groups. (F) The subgroup map predicted the response to ICB therapy between the two risk groups. (G) ROC curves verified the accuracy of the risk score in predicting the OS for patients in the IMvigor210 cohort. (H) The K-M survival curves of OS in patients between the high and low-risk groups. (I) Differential expression of CD274(PDL1) between high and low-risk groups. (J, K) Differences in the distribution of patients with and without response to ICI immunotherapy between high and low-risk groups. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease. CR/PR was identified as responders, and SD/PD was identified as non-responders. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001.

Validation of the clinical application value of the scoring system in the IMvigor210 (mUC) cohort

We selected the IMvigor210 (mUC) cohort to further validate the clinical application value of the scoring system. K-M survival curves indicated that the OS of patients with a low-risk score was longer than those with a high-risk score (Figure 6H). Moreover, the risk score can effectively predict the 1, 3, and 5-year OS rates of patients in the IMvigor210 cohort (Figure 6G). Then, the differential expression of PDL1 (CD274) between high and low-risk groups was also explored. The results showed that the expression level of PDL1 in the high-risk group was higher than that in the low-risk group (Figure 6I). Furthermore, the proportion of CR/PR and PD/SD in the high-risk group is higher than that in the low-risk group (Figure 6J, 6K).

COL1A2 is up-regulated in GBM

First, compared with the other eight FARGs, the mRNA expression of COL1A2 in GBM tissues was significantly higher than that in normal brain tissues (NBT), which was verified on the GEPIA website and three independent GBM cohorts (Figure 7A, 7B). Then, immunohistochemical staining of eight pairs of GBM tissues and adjacent tissues showed that COL1A2 was significantly up-regulated in GBM tissues (Figure 7C, 7D). In addition, at the tissue protein level, the expression COL1A2 in GBM tissues was also significantly higher than that in corresponding adjacent tissues (Figure 7E, 7F). Finally, the K-M survival curves of COL1A2 in four independent GBM cohorts (Rembrandt, GSE16011, TCGA, and CGGA cohorts) showed that the prognosis outcomes of GBM patients in the high expression group were significantly shorter than that of GBM patients in the low expression group (Figure 7G).

Figure 7. Basic analysis of COL1A2. (A) Summary of the differential expression analysis of nine FARGs in three independent GBM cohorts and the GEPIA website. Red represents up-regulated genes in GBM tissue, while blue represents a down-regulated gene in GBM tissue. Grey represents no statistical difference. (B) Differential expression analysis of COL1A2 between GBM tissues and NBTs in TCGA (HG-UG133A), Kamoun and Rembrandt cohorts, and GEPIA website. (C) Immunohistochemical analysis of COL1A2 protein in eight pairs of GBM tissues and adjacent nontumor tissues. (D) The IHC scores of COL1A2 in GBM tissues and matched adjacent nontumor tissues. (E, F) The protein expression level of COL1A2 in GBM tissues was significantly higher than that in corresponding adjacent tissues. (G) K-M survival curves of OS in GBM patients between the high and low COL1A2 expression groups. *P < 0.05, **P < 0.01.

COL1A2 promotes the malignant progression of GBM cells in vitro

By comparing the differential expression of COL1A2 between different cell lines at the transcription level and protein level, we found that the expression level of COL1A2 was the highest in U87 cells and the lowest in U251 cells (Figure 8A–8C). Then, we constructed knockdown and overexpression plasmids of COL1A2 to further explore the effect of COL1A2 on the biological behavior of GBM in vitro. The result of RT-qPCR and western blotting showed that the constructed knockdown and overexpression plasmids could significantly reduce or increase the expression of COL1A2 (Figure 8D–8I). The results of colony formation assay (Figure 8J–8M), CCK-8 assay (Figure 8N, 8O), and EdU staining (Figure 8P, 8Q) suggest that COL1A2 can significantly enhance the proliferation of GBM cells. In addition, the outcomes of wound healing and transwell assay support COL1A2 promoting the migration and invasion of GBM cells (Figure 8R–8U). These results suggest that COL1A2 plays an important role in the progression of GBM cells.

Figure 8. COL1A2 promotes the proliferation, migration, and invasion of GBM cells in vitro. (A–C) Differential expression of COL1A2 at transcriptional (A) and protein levels (B, C) between different GBM cell lines. (D–F) The efficiency of COL1A2 knock-down plasmid at transcriptional (D) and protein levels (E, F) was verified in U87 cells. (G–I) The efficiency of the COL1A2 overexpression plasmid at transcriptional (G) and protein levels (H, I) was verified in U251 cells. In (A–I), data are presented as the mean ± SD. X-axis: four different cell lines. Y-axis: expression level of genes. (J–Q) The proliferation activity of U251 and U87 cells were detected by colony formation assay (J–M), CCK-8 assay (N, O), and EdU staining (P, Q). (R–U) The migration and invasion ability of U87 and U251 cells were detected by wound healing (R, S) and Transwell assays (T, U). The representative photographs were photographed under a microscope. *P<0.05, **P<0.01, ***P<0.001.

Detection of the viability of GBM cells in vitro

In eight cases of GBM tissues, we found that the tissues with high expression of COL1A2 were accompanied by high expression of CD8. On the contrary, the tissues with low expression of COL1A2 were accompanied with low expression of CD8 (Figure 9A), which to some extent suggested that the expression of COL1A2 might induce the infiltration of CD8 T cells or the infiltration COL1A2. Then, we studied whether COL1A2 played a of CD8 T cells might induce the expression of role in the process of Jurkat cells acting on GBM cells in vitro. Both methods, PMA, and ionomycin or anti- CD3 and anti-CD28 antibodies, can effectively activate Jurkat cells into CD8 + Jurkat cells (Figure 9B). CD8 T cells transiently overexpress CD69 in the early stage of activation, which is considered to be a costimulatory signal of T cell proliferation. However, we did not detect the expression of CD69 protein in Jurkat cells stimulated by anti-CD3 and anti-CD28 antibodies (Figure 9C), so we utilized PMA and ionomycin to stimulate Jurkat cells. The results of the cell viability assay showed that activated Jurkat cells could effectively inhibit the viability of U87 (Figure 9D) and U251 (Figure 9E) cells. With the increase of COL1A2 expression, the lethality of activated Jurkat cells to GBM cells decreases, which may be related to COL1A2 promoting the proliferation of GBM cells.

Figure 9. COL1A2 inhibits the effect of activated Jurkat cells on the viability of GBM cells in vitro. (A) Immunofluorescence staining of the nucleus, COL1A2, and CD8 in eight cases of GBM tissues. (B, C) Jurkat cells were stimulated by PMA and ionomycin, or by anti-CD3 and anti-CD28 antibodies. After 24 hours of stimulation, the expression of CD8 and CD69 proteins were detected. (D, E) The viability of U87 (D) and U251 (E) cells were analyzed by measuring luciferase activity. ***P<0.001.

Data sources

The flow diagram of our study was shown in Figure 1. RNA-seq data, somatic mutation data, and corresponding clinical annotation of glioma patients were downloaded from the TCGA database, CGGA database, GEO database, UCSC Xena portal, GlioVis website, and NIH website. Before comparing and analyzing the gene expression data from different platforms, we performed transcripts per kilobase million (TPM) values transformation and robust multichip averaging (RMA) values transformation on RNA-seq data and microarray data, respectively. Patients with incomplete clinical information (including OS, OS status, gender, age, and WHO grade) will not be included in the clinical correlation analysis. Demographics and clinical information of 1813 glioma patients are shown in Supplementary Table 1. Somatic mutation data were processed with the R package “maftools” [57]. The genomic identification algorithm of important targets in tumors is used to process CNV data [58]. We searched 199 FARGs from the Molecular Signatures Database and the details of these FARGs can be found in Supplementary Table 2. In addition, the IMvigor210 cohort treated with the PD-L1 inhibitor was also included in this study to confirm our findings.

Acquisition of GBM tissue samples

The current research was authorized by the Ethics Committee of The First People’s Hospital of Fuzhou City and informed consent was obtained from the patients. The surgical specimens of GBM patients were removed to liquid nitrogen immediately after surgical excision.

Unsupervised consensus clustering based on FARGs

Different glioma cohorts have different sequencing methods and gene annotation information. To facilitate follow-up analysis, we took the intersection of 199 FARGs in three independent glioma cohorts and finally obtained 192 FARGs. A novel molecular classification was identified based on the 192 FARGs extracted from the TCGA cohort by unsupervised consensus clustering with 1000 iterations using the “ConsensusClusterPlus” package [59]. The optimal value of molecular classification should consider not only the rate of increase of the area under the cumulative distribution function (CDF) curves but also consider the correlation between different molecular classifications.

Gene set variation analysis (GSVA) between different subtypes

After downloading “c2.cp.kegg.v7.2.symbols” from the Molecular Signatures Database, GSVA between different subtypes was performed using the R package “GSVA” [60]. Differences in the relative activity of functional pathways between different subtypes were explored by using the R package “limma” [61].

Immune-related analysis between different subtypes

Different types of cell components in tumor tissues can be estimated by calculating TME-related scores using the R package “ESTIMATE” [62]. The content of different types of immune cells was reckoned using MCPCOUNTER, CIBERSORT, and QUANTISEQ algorithms [63, 64]. Enrichment scores of 29 immune signatures were determined by using a single sample gene set enrichment analysis (ssGSEA), which to some extent represents the immune activity within the tumor. [65]. The analyzed immune checkpoint proteins (ICPs) were selected from a previous study [66]. The above immune-related characteristics were analyzed in different subtypes.

Construction of a scoring system and verification of its prognostic value

We screened 175 FARGs with prognostic significance from 192 FARGs using uni-Cox regression analysis. 175 FARGs were further screened by lasso algorithm and multi-Cox regression analysis, and nine FARGs were finally obtained to construct a scoring system: risk score = sum (gene expression × coefficient). High and low-risk groups were classified according to the median value of the risk score. The prognostic value of the scoring system was evaluated by using Kaplan–Meier (K–M) survival curves (R package “survminer”), receiver operating characteristic (ROC) curves (R package “survivalROC”), uni- and multi-Cox regression analysis (R package “survival” and “forestplot”). The different distribution of different clinicopathological characteristics between high and low-risk groups was analyzed by R package “limma”.

Genomic variation analysis between high and low-risk groups

The mutation type and frequency of genes between high and low-risk groups were analyzed by the R packages “maftool” [57] and “GenVisR” [67]. The relationship between TMB, as a marker for the therapeutic efficacy of anti-PD-1 antibodies in other cancers [32, 68], and risk score was also analyzed by R package “ggpubr”. As genomic alterations can impact tumor immunity and immune infiltration patterns [69, 70], we compared amplifications and deletions between high and low-risk groups and visualized the results as circle graphs using the R package “RCircos” [71].

Immune-related analysis of risk score

At present, only a sub-fraction of patients achieved long-lasting clinical benefits from ICBs treatment. To effectively predict the response of tumor patients to ICBs, the TIDE algorithm was developed mainly to model two primary mechanisms of tumor immune evasion: inducing the dysfunction of cytotoxic T lymphocytes (CTLs) in tumors and preventing the infiltration of CTLs in tumor tissues [72]. T cell dysfunction was identified by measuring the interaction between each gene and the infiltration level of CTL to influence patient survival. The TIDE algorithm explored the association between gene expression data and markers of T cell dysfunction. Tumor samples that were highly positively correlated with markers of T cell dysfunction were identified as non-responders and otherwise as responders. Finally, the TIDE algorithm and an unsupervised subclass mapping method [73] were used together to forecast the response of glioma patients to anti–PD-1 and anti–CTLA-4 immunotherapy.

Cell culture, transfection, and activation

Normal human astrocytes (NHA) and GBM cell lines (U251, LN229, U118, and U87) were obtained from the Shanghai Institute of Biosciences and Cell Resources Center. Jurkat cells were purchased from The Global Bioresource Center. Except that Jurkat cells were cultured with RPMI-1640 medium (Gibco, USA), other cell lines were cultured with MEM (Gibco, USA) or DMEM (Gibco, USA). 50 IU/ml penicillin, 50 μg/ml streptomycin, and fetal bovine serum (FBS, Gibco, USA) were added to the culture medium. The COL1A2 knockdown and overexpression plasmids were constructed by the Sheweisi Biotechnology Company (Tianjin, China). U6-MCS-SV40 polyA-CMV-EGFP-SV40-NeoR was used to construct COL1A2 overexpression plasmid, and its three target sequences are shCOL1A2-1, 5’-GGTGTAAGCGGTGGTGGTTAT-3’, shCOL1A2-2, 5’-GCACTATGGATGCTATCAAAG-3’ and shCOL1A2-3, 5’-GCAACAGCAGGTTCACTTACA-3’ respectively. CMV-MCS-EF1a-ZsGreen1-SV40-Neomycin was used to construct COL1A2 overexpression plasmid, and its upstream and downstream primers are “AGCTGGCTAGCGTTCTCGAGGCCACCATGCTCAGCTTTGTGGATACGC” and “GTCTTTTTATTGCCGGGTACCTTATTTGAAACAGACTGGGCCAATG” respectively. Lipofectamine 3000 (Thermo Fisher, L3000075, USA) was used to transfect plasmids when the cell adhesion density reached about 80%. Two methods were used to activate Jurkat cells. Method 1: Jurkat cells (1 × 10⁶/ml) were stimulated by adding 20 ng/ml phorbol 12-myristate 13-acetate (PMA, Acmec, P33390-1mg) and 500 ng/ml ionomycin (70-CS0002, MultiSciences, China). Method 2: Jurkat cells were stimulated by adding 5μg/ml of anti-CD3 antibody (17617-1-AP, Proteintech, China), which was coated in a 6-well plate with the 5μg/ml of anti-CD28 antibody (65099-1-Ig, Proteintech). The expression of CD8 protein was used to detect whether Jurkat cells were activated, and the expression of CD69 protein was used to detect whether Jurkat cells proliferated. Jurkat and GBM cells were co-cultured (Jurkat: GBM cells = 1:1) in 96-well plates with a total volume of 200 μl.

RT-qPCR and Western blotting

The primer sequences of COL1A2 and GAPDH (internal reference gene) purchased from Ribobio (Guangzhou, China) are as follows: COL1A2 F primer: 5’-GGCCCTCAAGGTTTCCAAGG-3’, R primer: 5’-CACCCTGTGGTCCAACAACTC-3’; GAPDH F primer: 5’-TGTGGGCATCAATGGATTTGG-3’, R primer: 5’-ACACCATGTATTCCGGGTCAAT-3’. The primary antibodies include anti-COL1A2 (Servicebio, GB13022-2, 1:3000 dilution), anti-CD8 (Proteintech, 66868-1-Ig, 1:4000 dilution), anti-CD69 (Proteintech, 10803-1-AP, 1:1500 dilution) and anti-β tubulin (Proteintech, 11224-1-AP, 1:8000 dilution). The secondary antibodies included goat anti-rabbit IgG (Abcam, ab6721 1:5000 dilution) and goat anti-mouse IgG (Abcam, ab6789 1:10000 dilution). The extraction and detection methods of protein and total RNA in different cell lines were completely consistent with those described in our previous study [74]. Here, we will describe the extraction steps of tissue protein in detail. First, the tissues were weighed and sheared into 2ml EP tubes. RIPA and PMSF were mixed in a ratio of 100:1 to prepare tissue lysates. 1ml of lysate was added to every 50mg of sheared tissue and homogenized with a homogenizer, and then lysed on ice for 30 minutes. Finally, the mixture of lysate and tissue was centrifuged with a low-temperature high-speed centrifuge at 4° C and 12000 rpm for 10 minutes. The subsequent steps are the same as the extraction of cellular proteins. Both RT-qPCR and western blotting assays were repeated three times.

Immunohistochemical and immunofluorescence staining

First, eight pairs of GBM and adjacent tissues were fixed with 10% formalin for seven days, and then these tissues were embedded with paraffin and sectioned. After dewaxing and dehydration, the tissue sections were treated with 3% hydrogen peroxide for 10 minutes. Subsequently, the tissue sections were blocked with 5% BSA and then incubated overnight with a primary antibody against COL1A2 (Servicebio, GB13022-2, 1:1000 dilution) at 4° C. Then the tissue sections were treated with the corresponding secondary antibodies at room temperature for one hour. Finally, DAB staining, target molecule detection, and hematoxylin re-staining were carried out in turn. For immunofluorescence staining, tissue sections were immune-stained overnight with primary antibodies against COL1A2 (Servicebio, GB13022-2, 1:1000 dilution) and CD8 (Proteintech, 1:400) at 4° C, and then incubated with fluorochrome-conjugated antibodies. DAPI was added as a nuclear counterstain.

Cell proliferation assay

Cells were seeded into 6-well plates at 1000 cells/well and cultured for about 15 days. It was predicted that the cells would proliferate for about 5-7 generations. The culture medium was changed every three days, and the colony formation was closely observed. Cell culture was stopped when the number of cells in a single colony approached 50. The colonies were fixed with 4% ice-precooled paraformaldehyde for about 20 minutes, then stained with 0.1% crystal violet for about 20 minutes.

The Cell Counting Kit-8 (CCK-8, Beyotime, Shanghai, China) was employed to assay cell proliferation activity. Cells were seeded into 96-well plates at 2000 cells/well. Adding 10ul CCK-8 reagent to each well and continuing to culture in the incubator for two hours, the absorbance of cells was measured at 450nm every 24 hours for 5 consecutive days.

2 × 104 cells were seeded into each well of the 24-well plate and cultured until the cell adhesion concentration reached about 75%. BeyoClick™ EdU Cell Proliferation Kit (Beyotime, Shanghai, China) was used to measure cell proliferation activity according to the manufacturer’s instructions. The nuclei of all cells with blue fluorescence and the positive cells with red fluorescence were photographed by fluorescence microscopy, and the results were analyzed by ImageJ software.

Cell migration and invasion assay

When the adhesion concentration of U251 and U87 cells in the six-well plates reached 85-90%, the tip of a 200μl sterile spear cuts through the cell layer at the bottom of the plate to form an artificial wound. After washing the non-adherent cells with phosphate-buffered saline (PBS), the adherent cells were further cultured in a serum-free medium. At 0 and 24 hours, wound closure was photographed using an inverted Leica microscope.

Transwell chambers (Corning, USA) were utilized to assay the invasion and migration of GBM cell lines. First, the Matrigel (Corning, 356234, USA) with a concentration of 500 μg/ml was used to cover the upper chambers. 200μl serum-free medium and 600μl medium containing 10% FBS were added to the upper chamber and the lower chamber, respectively, and 8×10⁴ transfected cells were seeded in the upper chamber. After 24 hours of incubation at 37° C and 5% CO₂ for 24 hours, the non-invasive cells on the upper chambers were removed with a microcell scraper. The following operations were carried out at room temperature. The cells at the bottom of the upper chamber were fixed with 4% ice-precooled formaldehyde (Solarbio, P1110, China) for 30 minutes and then washed twice with PBS, and then the chambers were placed in 0.1% crystal violet (Solarbio, G1075, China) for 20 minutes. Finally, the chambers were washed with PBS, dried, and photographed with a microscope (Leica Microsystems, D-35578). In addition, we also performed the migration assay in the same way as the invasion assay, except that the upper chamber was not covered with Matrigel. Both invasion and migration assays were repeated three times.

Tumor cell viability assay

In light of the instructions offered by the manufacturer, the lentivirus expressing luciferase designed by Sheweisi Biotechnology Company (Tianjin, China) was transfected into U251 and U87 cells and screened with neomycin. Then COL1A2 knockdown and overexpression plasmids were transfected into cells stably expressing luciferase, respectively. When these cells grew to about 90%, they have inoculated with Jurkat cells in a ratio of 1:1 to 96 well plates with a white and clear bottom (total 1.0×10⁴ cells/ml). Co-cultured cells were taken at different time points (8, 16, 24, 36, 48 h), the supernatant was discarded, washed with PBS, cell lysate was added, and then D-luciferin potassium salt (D1009, UElandy, China) was added. After the above operations are completed, the fluorescence intensity is detected using a multifunctional microplate reader. The experiment was repeated three times.

Statistical analyses

Spearman and Pearson analysis methods were used to determine correlations between different grouping variables. The unpaired Student’s t-test and Mann–Whitney U-test were used to analyzing normally and non-normally distributed data, respectively. The Wilcoxon rank-sum test was used to compare two groups or categories; for more than two groups or categories, the Kruskal–Wallis test was used. Survival was assessed by K–M analysis with the log-rank test. Cox regression analyses were carried out to evaluate the prognostic value and stability of the risk prediction model. ROC curves were used to assess the prognostic value of 1-, 3-, and 5-year OS. TIIC infiltration level was analyzed using three algorithms (CIBERSORT, QUANTISEQ, and MCPCOUNTER). All statistical tests were two-sided and P<0.05 was considered statistically significant.

Availability of data and materials

The Cancer Genome Atlas (TCGA) database: https://portal.gdc.cancer.gov/. Chinese Glioma Genome Atlas (CGGA) database: http://www.cgga.org.cn/. Gene Expression Omnibus (GEO) databse: https://www.ncbi.nlm.nih.gov/geo/. National Human Genome Research Institute Home (NIH): https://www.genome.gov/. Data Visualization Tools for Brain Tumor Datasets (GlioVis): http://gliovis.bioinfo.cnio.es/. Molecular Signatures Database: https://www.gsea-msigdb.org/gsea/index.jsp. Gene Expression Profiling Interactive Analysis (GEPIA) database: http://gepia.cancer-pku.cn/detail.php. Tumor Immune Dysfunction and Exclusion (TIDE): http://tide.dfci.harvard.edu/. IMvigor210 cohort: http://research-pub.gene.com/IMvigor210CoreBiologies).