2.1 Identification of the common DEGs between the datasets used
We obtained three gene expression profiles (GSE33630, GSE58545, and GSE60542) from the GEO database. GSE33630 includes 45 normal thyroid samples and 60 THCA samples; GSE58545 includes 18 normal thyroid samples and 27 THCA samples; GSE60542 includes 30 normal thyroid samples and 33 THCA samples. According to the conventional criteria (adjusted P < 0.05 and |log2FoldChange| ≥ 2.0), 263 genes were identified as DEGs in GSE33630, of which 133 were downregulated and 130 were upregulated; moreover, GSE58545 included 270 DEGs, including 144 downregulated genes and 126 upregulated genes, and GSE60542 contained 228 DEGs, including 107 downregulated genes and 121 upregulated genes (Figure 1A-C). A Venn diagram showed that 98 genes were differentially expressed in the three data sets, of which 46 were downregulated and 52 were upregulated (Figure 1D, E).
2.2 Functional enrichment analysis of the common DEGs
In this study, we performed GO functional analysis of the common DEGs using the tool DAVID. Then, we filtered the results to improve the confidence according to the standard criteria (P < 0.05 and gene counts ≥ 5) (Table S1). GO analysis showed that the common DEGs were mainly enriched in the CC category, including plasma membrane, extracellular exosome, extracellular region, and extracellular space. DEGs annotated as BP were enriched in cell adhesion, signal transduction, extracellular matrix organization, positive regulation of gene expression, positive regulation of cell proliferation, blood coagulation, wound healing, platelet degranulation, and nervous system development, all of which are processes associated with the occurrence and development of tumors.
2.3 Identification and analysis of the hub genes
We selected 46 DEGs involved in specific BP (P < 0.05) to build the PPI network (Figure 2A, Table 1). The most important module was obtained using the plugin MCODE in Cytoscape (Figure 2B). The top eight genes, including FN1, TIMP1, SERPINA1, COMP, PROS1, MMRN1, KIT, and TNFRSF11B, were identified as potential hub genes according to the degree score generated by the plugin CytoHubba (Figure 2C, Table 2). This was consistent with their enrichment in the top module determined using MCODE.Among of them, FN1 had the highest degree of connectivity in the PPI network. Furthermore, logrank regression and multivariate Cox regression analysis were used to calculate the correlation between gene expression and PFS (Figure D-E), indicating FN1 appeared to be the most attractive drug target and prognostic marker., GSEA was used to perform kegg analysis for FN1. The results suggested that the most of the involved significant pathways included chemokine signaling pathway, cytokine cytokine receptor interaction, leishmania infection, natural killer cell mediated cytotoxicity, and T cell receptor signaling pathway, as it has been established that FN1 plays a crucial role in tumor architecture and controls metastasis (Figure 2F) (14).
2.4 Expression levels of FN1 in THCA and evaluation of its value as a prognostic marker
We analyzed the expression levels of FN1 in THCA using the Oncomine database. We found that FN1 expression was upregulated in almost all different subtypes of THCA, including PTC and ATC (Figure 3A). We validated these results using TCGA database, in which FN1 was also significantly more expressed in THCA samples than in normal thyroid samples (P < 0.05) (Figure 3B). FN1 expression levels were correlated with PFS (P < 0.05) (Figure 3C). The area under the curve (AUC) of FN1 from TCGA dataset was 0.8971 (Figure 3D), highlighting the value of FN1 as a diagnostic marker in THCA. We compared the clinical characteristics (including age, sex, location, clinical stage, progression state, lymph node metastasis, and distant metastasis) between the FN1-high expression and FN1-low expression groups, and observed statistically significant differences in lymph node metastasis, distant metastasis, and progression state, although no significant differences were identified for other clinical features (Table 3). Furthermore, univariate and multivariate Cox regression models revealed that clinical stage and FN1 expression were independent prognostic factors for PFS in patients with THCA (Table 4). In addition, analysis of the UALCAN database showed that FN1 expression levels were closely correlated to cancer stage (Figure 3E) and that the degree of methylation of the FN1 promoter was lower in THCA than that in normal tissues (Figure 3F), indicating that FN1 may be involved in the development of THCA.
2.5 Analysis of the potential genetic andepigenetic alterations underlying FN1 dysregulation
Next, we investigated the underlying mechanism of FN1 dysregulation in THCA. To determine whether copy number alterations (CNAs) are responsible for the abnormal expression of FN1 in THCA, we analyzed 397 cases from TCGA database for which CNAs data was available. No differences were observed for CNAs in the FN1-high and FN1-low groups (Figure 4A). Another mechanism that could underlie the altered FN1 expression profile observed in THCA samples is promoter hypermethylation, which plays an important role in the occurrence and development of several types of tumors [33, 34]. To investigate whether DNA methylation results in FN1 dysfunction, we examined the status of CpG sites in 507 THCA cases from TCGA database using the tool MEXPRESS. We found that 27 CpG sites had associated data; among which 20 CpG sites showed a negative correlation with FN1 expression (Figure 4C). The Pearson correlation coefficient was calculated for the five CpG sites with the highest correlation coefficient, including cg21494132, cg09040552, cg11309217, cg03228449, and cg03228449 (Figure 4D, all P < 0.05). Our results demonstrate the correlation between DNA methylation and abnormal expression levels of FN1 in THCA, highlighting the need for further investigation of the underlying mechanism.
2.6 Correlation between FN1 and inflammatory activities
Considering the strong association between FN1 expression levels and THCA prognosis, we hypothesized that FN1 may be associated with inflammatory responses, leading to enhanced survival rates. To identify the FN1 associated immune signature in THCA, we assessed the degree of immune infiltration with CIBERSORT. FN1 expression was closely correlated to the degree of M2 macrophages, resting memory CD4+ T cells, follicular helper T cells, and CD8+ T cell infiltration. The increase in FN1 expression was associated with an increase of the proportion of M2 macrophages and resting memory CD4+ T cells and a decrease of the proportion of follicular helper T cells and CD8+ T cells (Figure 5A). To further investigate the correlation between FN1 and inflammation, we analyzed the co-expression of FN1 and members of the B7-CD28 ligand-receptor family, including CD274, CD80, CD86, CD276, CD273, CD275, B7-H4, B7-H5, CD28, B7-H7, CD152, CD279, CD278, TLT-2, and NKp30, which are closely correlated to T cell function. The result demonstrated that FN1 exhibited a significant co-expression trend with CD273, CD274, CD275, CD276 and B7-H4 (Figure 5B). Furthermore, we investigated the correlation between FN1, CD273, CD274, CD275, CD276, and B7-H4 expression in the THCA cohort, and found a close positive correlation between FN1 and CD276 expression (R = 0.55, P = 0) (Figure 5C-D). This result was validated using TIMER and GEPIA, which provided R values of 0.682 and 0.79, respectively, for the correlation between FN1 and CD276 in THCA (P < 0.001) (Figure 5E-F).
2.7 Down-regulation of FN1 inhibited cell proliferation and invasion and decreased CD276 expression levels in THCA samples
To explore the biological significance of FN1 in THCA tumorigenesis, KTC-1 and B-CPAP cells were transfected with siRNA targeting FN1 (siFN1) or negative control siRNA (siNC). Efficient depletion of FN1 expression was confirmed via RT-qPCR (P < 0.05, Figure 6A). Moreover, we found that downregulation of FN1 significantly reduced the expression of CD276 in KTC-1 and B-CPAP cells compared with siNC transfection (Figure 6A). Then, we studied the effect of FN1 on THCA cell proliferation and invasion in vitro. The wound-healing assay, cell viability and cytotoxicity CCK-8 assay, and clone formation assay revealed that downregulation of FN1 in both cell types significantly inhibited cell proliferation and invasion compared to that in the control cells (P < 0.05, Figure 6B-D). These data suggest that downregulation of FN1 reduces the viability of THCA cell lines.
2.8 FN1 is positively correlated with CD276
In order to assess FN1 correlation with CD276, we analyzed FN1 and CD276 expression in tumor sites and the adjacent no-tumor samples (Figure 7A). We found that FN1 and CD276 showed significantly higher expression in tumor sites than in the adjacent no-tumor samples (P<0.05, Figure 7C-D). Furthermore, we divided the tumor sites into two groups according to FN1 expression and found that CD276 levels were positively correlated with the level of FN1(P<0.05, Figure 7B, E).