FNDC3B mRNA expression in various cancers
We first assessed the expression of FNDC3B in different tumors and normal tissues of multiple cancer types using the Oncomine database. The results showed that abnormal FNDC3B expression was retrieved from a total of 368 datasets. Among them, FNDC3B expression levels were significantly upregulated in tumor tissues in 43 datasets, including brain and central nervous system (CNS), head and neck, esophageal, kidney, cervical, bladder, and colorectal cancers, etc. (Figure 2A, P < 0.0001, Fold Change > 2, and Gene Rank < Top 10%). In addition, its expression in leukemia, lymphoma, breast cancer, and sarcoma lymphoma was shown to be downregulated in multiple datasets. In summary, FNDC3B is generally upregulated in several tumors. In the brain and CNS cancers dataset, there were seven studies on the upregulation of FNDC3B and no studies on its downregulation. Comparison of FNDC3B across the seven studies showed a Median Rank = 302, suggesting that FNDC3B was highly expressed in glioma tissues and concentrated both in LGG and GBM (Figure 2B, P <0.0001). As for the outlier analysis of FNDC3B, 822 different studies on FNDC3B have been included in Oncomine database. Among them, 11 studies indicated the upregulation of FNDC3B and four studies the downregulation in brain and CNS cancers (Figure 2A). Then, we evaluated the differentially expressed level of FNDC3B in TCGA pan-cancer data using GEPIA2. Our findings revealed that an elevated expression of FNDC3B was involved in a variety of tumors (Figure 2C). In TCGA LGG and GBM cohorts, FNDC3B expression was significantly higher in tumors compared to matched normal tissues (Figure 2D). Specifically, FNDC3B expression was 2.71-fold in LGG and 7.16-fold in GBM vs. normal brain tissue. Moreover, the expression of FNDC3B was negatively correlated with FNDC3B DNA methylation; the methylation levels were reduced in the FNDC3B high group based on the TCGA-LGG dataset (Figure 2E-2F). These results suggested that FNDC3B may be negatively regulated by methylation, which may lead to its high expression in glioma samples.
FNDC3B protein expression was explored using HPA database. The immunohistochemistry result revealed upregulated FNDC3B in glioma samples (Figure 2G). In human cancer tissues, FNDC3B protein expression in glioma was ranked as the top 12 out of 20 distinct cancer types (Supplementary Figure 1A). Furthermore, we detected the FNDC3B mRNA level in human normal tissues using the GTEx database, and FNDC3B expression was mainly found in lung, adipose tissue, thyroid gland, endometrium, and ovary. It is worth noting that the brain displayed the lowest expression levels of global FNDC3B transcript across all normal tissues (Supplementary Figure 1B). In addition, we systematically elucidated the expression levels of the FNDC3B in different cancer cells by querying the CCLE database and found that FNDC3B was highly expressed in glioma cell lines (Supplementary Figure 1C).
High expression of FNDC3B predicts poor prognosis of glioma
OncoLnc and GEPIA2 are free online resources and databases for the analysis and visualization of datasets from the TCGA and GTEx projects. To examine the function of FNDC3B on OS in various cancers, we used OncoLnc online tool to perform Cox regression analysis. We found that FNDC3B expression in LGG was ranked first among 21 different cancer types based on the FDR correction (Table 1). Moreover, we analyzed the relationships between FNDC3B expression and prognostic values in 33 types of cancer using GEPIA2 databases. As shown in Supplementary Figure 2, high FNDC3B expression levels were associated with poorer prognosis of OS and disease-free survival (DFS) in adrenocortical carcinoma (ACC), GBM, kidney Chromophobe (KICH), LGG, liver hepatocellular carcinoma (LIHC), mesothelioma (MESO), and pancreatic adenocarcinoma (PAAD); OS in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) and lung adenocarcinoma (LUAD); DFS in colon adenocarcinoma (COAD) and uveal Melanoma (UVM). A high FNDC3B expression was correlated with a better prognosis of OS and DFS in skin cutaneous melanoma (SKCM) as well as OS in acute myeloid leukemia (LAML).
To assess the prognostic significance of FNDC3B in glioma patients from TCGA and CGGA, samples were first split into two groups according to the median expression of FNDC3B for each dataset. Among pan-glioma, LGG, and GBM in the TCGA datasets, patients with higher FNDC3B levels presented shorter OS and DFS (Figure 2H) compared to patients expressing low levels of FNDC3B. Similarly, high FNDC3B expression was significantly associated with poor prognosis in all the three CGGA datasets (seq_325, seq_693, array_301) (Supplementary Figure 3). In particular, highly expressed FNDC3B was significantly related to reduced DFS in GBM (p<0.001), however, just marginally correlated with worse OS (P<0.05). We postulate that the less obvious but significant results of GBM may be due to the insufficient statistical power of a small sample size. Furthermore, we investigated the associations between FNDC3B expression and clinical characteristics, such as age, gender, grade, and isocitrate dehydrogenase (IDH) status. FNDC3B expression was higher in high-grade and IDH wildtype patients; there was no significant difference between age and gender based on the TCGA datasets (Figure 2I). Our findings revealed that higher FNDC3B expression is closely correlated with the malignant clinical characters of gliomas.
Subsequently, univariate and multivariate Cox regression analyses were performed to identify whether FNDC3B expression represented an independent prognostic factor. Univariate Cox analysis showed that FNDC3B (HR = 1.64; 95% CI = 1.50-1.80; P < 0.001), grade (HR = 3.37; 95% CI = 2.28-4.98; P < 0.001) and age (HR = 1.06; 95% CI = 1.04-1.07; P < 0.001) were high-risk factors, and IDH mutation (HR = 0.18; 95% CI = 0.07-0.48; P < 0.001) was a low-risk factor (Table 2). In multivariate Cox regression analysis, FNDC3B was independently associated with overall survival, suggesting it could be an independent prognostic biomarker for glioma (HR = 1.72; 95% CI = 1.07-2.80; P < 0.05). In addition, age may also be an independent prognostic factor (Table 2 and Figure 2J).
Construction and validation of a prognostic nomogram
A quantitative prognostic nomogram model to predict individual survival chances was established based on the TCGA-LGG dataset using Cox regression (Figure 3A). According to the stepwise Cox multivariate regression analysis, age, grade, IDH status, and FNDC3B expression were features that were included in the nomogram, and the risk scores were calculated based on that model. The concordance index (C-index) for OS prediction was 0.775, indicating high predictive performance of the model. A calibration curve was implemented to reflect the degree of consistency between the predicted risk and actual occurrence risk, and it can be used to estimate the accuracy of the model in predicting the probability of an individual outcome in the future . In our study, the calibration curve showed acceptable agreement between nomogram-predicted and observed 2-, 3- and 5-year OS in the CGGA_325, CGGA_693, and CGGA_301 validation cohorts (Figure 3B).
The patients in the training dataset were divided into high-risk and low-risk groups according to the positive and negative values of the risk score, respectively. The K-M survival curve showed a good discriminating ability of the nomogram (P = 1.4e-4) (Figure 3C). Moreover, the area under the ROC curve for OS was 0.755, indicating a reliable predictive ability in the TCGA LGG dataset (Figure 3D). The CGGA_325, CGGA_693, and CGGA_301 datasets were used to validate the performance of the nomogram. A risk score for each patient was generated by the same method. Consistently, the patients in the high-risk group had a notably poorer prognosis in all the three validation datasets (p < 0.0001) (Figure 3C). The area under the curve (AUC) for CGGA_325, CGGA_693, and CGGA_301 OS was 0.8, 0.816, and 0.807, respectively (Figure 3D). In conclusion, these results indicated that the nomogram had adequate performance in predicting the OS of glioma patients.
DEGs identification and weighted co-expression network construction based on FNDC3B expression
To further elucidate the role of FNDC3B expression in the glioma microenvironment, the median expression value was used to create a categorical variable for the TCGA-LGG cohort. The DEGs between the two groups were detected using Wilcoxon test. After setting FDR < 0.05 and fold change ≥ 2 in either direction, a total of 2,099 DEGs, including 1,631 upregulated and 468 downregulated genes were screened out in FNDC3B highly expressed group compared to the low expressed group (Figure 4A). In order to obtain more reliable results, 670 DEGs were filtered after setting the mean expression value to 1 for the raw DEGs (Figure 4B). GSEA was conducted to assess the potential functions of these DEGs. Our results suggested that various immune-related gene signatures were enriched in LGG samples, such as response to cytokine, cytokine secretion, immune system process, and inflammatory response (Figure 4C). According to the KEGG analysis, DEGs were mainly concentrated in PI3K-Akt, p53, and MAPK signaling pathways (Figure 4D). We used TIMER2.0 to explore the relationship between FNDC3B with the top 10 most significantly up or downregulated DEGs in pan-cancer. The corresponding heatmap showed that the correlation was consistent with the gene expression direction of the 20 DEGs only for LGG (Figure 4E).
WGCNA was applied to build a co-expression network based on the 2,099 DEGs. Before constructing the co-expression network, we screened the DEGs by setting expression value larger than 1 in at least 10% of all the samples, and the average level at least 0.5. Finally, 673 DEGs were filtered and selected for subsequent analysis. Power eight was chosen as the appropriate soft threshold because it was the first value to make the degree of independence reach 0.90 and the corresponding average connectivity was close to zero (Figure 5A-5B). A total of five gene modules were excavated (genes in the grey module that were not co-expressed; Figure 5C). As shown in Figure 5D, the turquoise module was the most relevant in FNDC3B expression level (R = 0.51, p = 7e-35). Therefore, the 282 hub genes included in the turquoise module were extracted for further analysis. GO enrichment analyses of these genes indicated that the response to interferon-gamma, neutrophil activation and degranulation, regulation of immune effector process, antigen processing and presentation, T cell activation, leukocyte migration, lymphocyte proliferation, mononuclear cell proliferation, interleukin-8 production and acute inflammatory response was related to FNDC3B-mediated immune events. The top 30 GO items, as ranked by their P-values, are shown in Figure 6A. A total of 210 nodes and 1,061 edges were mapped for the turquoise module genes in the PPI network (Supplementary Figure 4). Using CytoHubba in Cytoscape plug-in, we selected the top 10 genes ranked by the MCC method as hub genes, including TLR2, TLR7, PTPRC (CD45), CCR1, CCL5, TLR1, FN1, VCAM1 (CD106), CXCL10, and TLR6. They were immune-related genes and positively correlated with FNDC3B (R > 0.3 and P < 0.0001; Figure 6B).
A gene-gene interaction network for the 10 hub genes was built, and their functions were analyzed through the GeneMANIA database (Figure 6C). We found that toll-like receptor signaling pathways, ERK1 and ERK2 cascade and NIK/NF-kappaB signaling, were enriched in LGG. Then, OncoLnc online tool was applied to investigate the function of these hub genes on OS of LGG. The results showed that all hub genes were independent risk factors for evaluating OS (Table 3), and the K-M curves based on GEPIA2 displayed that higher expression of these genes predicted shorter OS in LGG (Supplementary Figure 5). Furthermore, the gene expression analysis of the LGG samples from TCGA showed that the combined expression of the 10 hub genes had a significant effect on overall survival (Figure 6D).
Relationship between FNDC3B expression and tumor immune infiltrates
Since previous studies have reported that TILs are independent predictors in cancers [19, 20], we used TISIDB database to infer the correlations between the expression of FNDC3B and the abundance of 27 types of TILs across TCGA pan-cancers. As shown in Figure 7A, FNDC3B expression was positively correlated with TILs in several human cancer types, especially in LGG. Moreover, we investigated the associations between FNDC3B expression and immune subtypes across human cancers, and the landscape of correlations between FNDC3B expression and immune subtypes in different types of cancer (Figure 7B). Among all cancer types, LGG showed the most significant results via Kruskal-Wallis test (p = 1.26e-23). In TISIDB, we further analyzed FNDC3B expression in different immune subtypes of LGG. We found FNDC3B was mainly expressed in three types, including C3 (inflammatory type), C4 (lymphocyte depleted type), and C5 (immunologically quiet type). FNDC3B expression was the highest in the C3 (inflammatory) type and the lowest in the C5 (immunologically quiet) type (Figure 7C). These results indicated that FNDC3B may play an important role in immune infiltration in glioma. Notably, FNDC3B expression was correlated with the abundance of central memory CD8 T cells (r = 0.497, p < 2.2e-16), effector memory CD8 T cells (r = 0.412, p < 2.2e-16), central memory CD4 T cells (r = 0.468, p < 2.2e-16), regulatory T cells (r = 0.521, p < 2.2e-16), natural killer (NK) cells (r = 0.532, p < 2.2e-16), natural killer T (NKT) cells (r = 0.64, p < 2.2e-16), memory B cells (r = 0.64, p < 2.2e-16) , and macrophages (r = 0.349, p < 2.2e-16) in LGG (Figure 7D). The positive correlations between FNDC3B expression and TILs were also observed in GBM.
We then assessed the relationships between FNDC3B and eight genes previously reported to be targets of immune checkpoint inhibitors, including CD274 (PD-L1), PDCD1 (PD-1), CD152 (CTLA-4), CD276 (B7-H3), HAVCR2 (TIM-3), CD223 (LAG-3), TNFRSF4 (OX40), and VTCN1. There were significant positive correlations between FNDC3B with B7-H3 (R = 0.69, p < 2.2e-16), PD-L1 (R = 0.58, p < 2.2e-16), TIM-3 (R = 0.43, p < 2.2e-16), PD-1 (R = 0.42, p < 2.2e-16), CTLA-4 (R = 0.34, p = 1.3e-15), and OX40 (R = 0.31, p = 1.3e-12) (Figure 8A). Thus, the six genes were significantly upregulated in the FNDC3B high group compared with the low group (Figure 8B). In summary, these results suggested that FNDC3B was correlated with clinically relevant immune checkpoint molecules in glioma.
Subsequently, to investigate whether FNDC3B expression was correlated with immune infiltration patterns in LGG, we compared the degree of immune cell infiltration between high and low expression groups using the ESTIMATE algorithm. The immune, stromal, and ESTIMATE scores were higher in the high-expression group than in the low-expression group (Figure 8C). Furthermore, we explored the proportions of 22 types of immune cells for LGG using CIBERSORTx to acquire a deeper understanding of the relationship between FNDC3B expression and tumor immune infiltrates. Among the 529 TCGA-LGG samples, 265 samples were in the high expression group and 264 samples in the low expression group. Figure 8D shows the differences in the proportions of the 22 subpopulations of immune cells in these two groups. Naive B cells, plasma cells, naive CD4 T cells, resting memory CD4 T cells, activated memory CD4 T cells, follicular helper T cells (Tfh), M1 and M2 macrophages, activated dendritic cells, and neutrophils were the main immune cells affected by FNDC3B expression. Among them, there were more proportions of resting memory CD4 T cells (p < 0.001), activated memory CD4 T cells (p < 0.01), M1 macrophages (p < 0.01), M2 macrophages (p < 0.01), dendritic cells activated (p < 0.01), and neutrophils (p < 0.0001) in the high expression group. In contrast, the proportions of naive B cells (p < 0.01), plasma cells (p < 0.0001), naive CD4 T cells (p < 0.05), and follicular helper T cells (p < 0.05) were lower in the high expression group compared with the low expression group. The correlation matrix heatmap of 22 immune infiltration cells in LGG samples was shown in Figure 8E. Moreover, we analyzed the association between FNDC3B expression and gene markers of various TILs, including B cell, plasma cells, T cell, CD4+ T cell, Tfh, M1 and M2 macrophages, dendritic cells, and neutrophils (Table 4). Overall, FNDC3B expression was strongly positively correlated with gene markers of B cells, T cells, M1 and M2 macrophages, dendritic cells, and neutrophils for TCGA-LGG and TCGA-GBM.