Identification of gene expression and DNA methylation of SERPINA5 and TIMP1 as novel prognostic markers in lower-grade gliomas

Lower-grade glioma datasets

TCGA LGG dataset was downloaded from the University of California Santa Cruz cancer browser https://genome-cancer.ucsc.edu/ (version: 2015-02-24) as training dataset. In total, 473 samples (225 grade II, 248 grade III gliomas) having clinical data were profiled for class discovery and survival analysis. A total of 131 samples (97 grade II, 34 grade III gliomas) from CGGA repository (http://cgga.org.cn/) was included in our analysis as validation dataset, and all samples’ clinical data were downloaded for survival analysis. Overall survival (OS) was defined as the time interval from resection until the date of death. Relapse-free survival (RFS) is the period from resection to the radiological evidence of first tumor recurrence.

Gene expression data analysis

Gene expression data of TCGA LGG are from the Illumina HiSeq 2000 RNA Sequencing platform, and all counts data is then log2(count+1) transformed. The differential gene expression analysis and the adjusting for multiple testing was performed with edgeR package (Robinson, McCarthy & Smyth, 2010). The CGGA microarray dataset was generated by Agilent Whole Human Genome Array, and probe intensities were normalized using GeneSpring GX 11.0 (Yan et al., 2012). In the TCGA LGG dataset, 299 patients were classified as short-term survivors (<2 years) and the remaining 174 patients as long-term survivors (≥2 years). In order to exclude the influence of loss of follow-up or short follow-up time, we excluded the samples that did not reach the end time, and only samples that had died and had a clear survival time were included for screening of differential expressed genes (DEGs). Thus, a total of 44 short-term survivors (<2 years) and 48 long-term survivors (≥2 years) were included in the gene differential expression analysis. Genes were considered to have significant difference in expression if —log fold change (FC)— ≥ 1.0 and adjusted P < 0.05. The prognostic value of DEGs generated from TCGA lower-grade glioma dataset were then validated using Kaplan–Meier survival analysis with the CGGA LGG microarray dataset.

DNA methylation analysis

DNA methylation data of TCGA LGG are generated by the Illumina Infinium HumanMethylation450 platform. Mapping between probes on the RNA-seq and DNA methylation probes on the methylation array was performed. The β value was used to estimate the methylation level of probes. Probes with β ≥0.5 was considered as hyper-methylated sites, and β <0.5 was considered as hypo-methylated sites. The correlation between gene expression levels and DNA methylation levels were assessed using Spearman’s correlation analysis. —Spearman r— ≥0.6 was indicated a strong correlation and P < 0.05 was considered as statistically significant (Zeng et al., 2018).

To investigate the correlation between gene expression and DNA methylation, we performed a parallel DNA methylation analysis of the candidate genes. Mapping SERPINA5 and TIMP1 to DNA methylation probes identified 23 and 14 methylation sites, respectively. To obtain differentially methylated sites, patients were divided into two groups according to the median of gene expression. Of the 37 methylation sites, 4 differential methylation sites (SERPINA5 cg15509705; TIMP1 cg27151711; TIMP1 cg16523424; TIMP1 cg04791822) were identified (Table S2).

Functional enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) was used to identify the potential biological functions of co-differentially expressed genes (co-DEGs) in LGGs (https://david.ncifcrf.gov/home.jsp) (Huang da, Sherman & Lempicki, 2009; Yan et al., 2019). Gene ontology (GO) analysis, including biological processes (BP), molecular function (MF), and cellular composition (CC), was performed to annotate these genes and determine functional enrichment. P <0.05 was considered as statistically significant for pathway enrichment. Then, the protein and protein interaction network of SERPINA5 and TIMP1 were constructed with GENEMANIA online database (https://genemania.org/) (Warde-Farley et al., 2010; Zhang et al., 2017).

Statistical analysis

Statistical analyses were performed with SPSS 13.0 (SPSS Inc., Chicago, IL, USA) and GraphPad Prism 5.0 (Graphpad Inc., San Diego, CA, USA). Nonparametric test was performed to identify genes that were differentially expressed between two groups. Spearman’s correlation analysis was used to determine the correlation between gene expression levels and DNA methylation status. Kaplan–Meier survival analysis was carried out to assess the survival distribution and the log-rank test was performed to determine the significance of the differences between two groups (Wang et al., 2019). For multivariable analysis, a Cox proportional hazards model was constructed for OS with a limited backward-LR procedure and was adjusted by potential confounding covariates. Hazard ratio (HRs) and 95% confidence intervals (CIs) were used to describe the risk. P <0.05 was considered as statistically significant.

Identification of prognostic DEGs in LGGs (Fig. 1)

We first sought to identify DEGs between 44 short-term survivors (<2 years) and 48 long-term survivors (≥2 years) in the LGGs of TCGA microarray dataset. In total, 106 genes (78 upregulated genes and 28 down-regulated genes) were identified to be differentially expressed (Fig. 1A, Table 1). The results of two-dimensional hierarchical cluster indicated that the mRNA expression profiles of short-term survivors and long-term survivors distributed in separate clusters (Fig. 1B).

Functional enrichment analysis of DEGs in LGGs

To explore the potential biological functions of prognostic related DEGs in LGGs, GO categories and KEGG Pathway enrichment analysis were performed using the DAVID online database. GO analysis revealed that 106 DEGs were significantly enriched in cell components such as proteinaceous extracellular matrix, extracellular space, extracellular exosome and basement membrane, and was involved in the biological processes such as multicellular organism development, thyroid gland development and anterior/posterior pattern specification (Table S1 ). KEGG Pathway enrichment analysis showed that DEGs mainly enriched in phenylalanine metabolism and histidine metabolism, but the result was not significant (Table S1).

Survival value of top significant DEGs in LGGs

Then, TCGA-LGG and CGGA-LGG datasets were used to verify the prognostic value of the top eight significant changed genes, including EMP3, FBXO17, METTL7B, SERPINA5, SSTR5, TIMP1, TMEM61, VAV3. In the TCGA LGG dataset, SERPINA5 high expression indicated the worse OS and RFS of LGG patients (Figs. 2D, 2L). However, TIMP1 high expression can only predict the shorter OS of patients and has no significant correlation with the RFS of LGG patients (Figs. 2F, 2N). Moreover, the expression of SERPINA5 and TIMP1 in the CGGA dataset were coincident with those of TCGA dataset, with significant differences in OS (Fig. 3). Additionally, SERPINA5 and TIMP1 mRNA expression were significantly increased with the increase of glioma grades in both TCGA LGG dataset and CGGA dataset (Fig. 4).

The association between expression of the candidate genes and clinical characteristics in CGGA lower-grade glioma patients is presented in Table 2. Univariate analysis indicated that karnofsky performance score (KPS) ≥70 was significantly associated with better survival in patients (HR = 0.127, 95% CI [0.037–0.434], P = 0.001), and tumor grade (grade III vs grade II) was significantly correlated with poor survival in patients (HR = 8.883, 95% CI [3.670–21.501], P = 0.000001). Multivariate analysis, adjusted for KPS and tumor grade, also showed that the top two DEGs were significantly associated with survival. TIMP1 high expression group exhibited poor survival as compared to low expression group (HR = 8.656, 95% CI [2.578–29.060], P = 0.014). Unfortunately, SERPINA5 high expression was not an independent poor predictor for OS of LGGs patients (HR = 0.473, 95% CI [0.192–1.164], P = 0.103).

DNA promoter hypermethylation silences SERPINA5 and TIMP1 mRNA expression

To investigate the correlation between gene expression and DNA methylation, we performed a parallel DNA methylation analysis of the candidate genes. Mapping SERPINA5 and TIMP1 to DNA methylation probes identified 23 and 14 methylation sites, respectively. To obtain differentially methylated sites, patients were divided into two groups according to the median of gene expression. Of the 37 methylation sites, 4 differential methylation sites (SERPINA5 cg15509705; TIMP1 cg27151711; TIMP1 cg16523424; TIMP1 cg04791822) were identified (Table S2). As shown in Fig. 5, the methylation status of 4 methylation sites were remarkably lower in high gene expression group than low gene expression group. Furthermore, the correlation analysis revealed that 4 methylation sites were negatively correlated with gene expression levels (Spearman r <  − 0.4, P < 0.0001, Fig. 5).

SERPINA5 and TIMP1 methylation are potent prognostic markers for LGGs

In order to identify the effect of these methylation sites on prognosis, we assessed the association between 4 methylation sites and prognosis with the TCGA LGG DNA methylation dataset. The samples were divided into two groups with methylation β value of 0.5 as the cut-off value and the prognostic difference were compared. As shown in Figs. 6A–6D, hyper-methylation of 4 methylation sites indicated better OS (SERPINA5 cg15509705: HR = 1.66, 95% CI [1.101–2.501], P < 0.0001; TIMP1 cg27151711: HR = 5.375, 95% CI [2.435–11.87], P < 0.0001; TIMP1 cg16523424: HR = 3.978, 95% CI [2.245–7.049], P <0.0001; TIMP1 cg04791822: HR = 7.284, 95% CI [2.458–21.59], P <0.0001). In addition, except for SERPINA5 cg15509705 (HR = 1.411, 95% CI [0.967–2.058], P = 0.0762, Fig. 6E), hypo-methylation of 3 other methylation sites the high-risk group exhibited significantly worse RFS (TIMP1 cg27151711: HR = 4.700, 95% CI [2.127–10.38], P <0.0001; TIMP1 cg16523424: HR = 3.037, 95% CI [1.738–5.307], P < 0.0001; TIMP1 cg04791822: HR = 5.653, 95% CI [1.936–16.51], P < 0.0001, Figs. 6F–6H).

Functional enrichment analysis of SERPINA5- and TIMP1-associated co-expressed genes

Then, we used GENEMANIA online database to identify the proteins interacting with SERPINA5 and TIMP1. As shown in Fig. 7A, SERPINA5 mainly interacts with PROC, which has serine hydrolase activity and functions as a negative regulation of hemostasis and coagulation. While TIMP1 interacts with extracellular matrix proteins MMP1, MMP3 and MMP9, and participants in the processes of extracellular matrix disassembly and organization, collagen catabolism and metabolism (Fig. 7B).