Five raw PTC-associated gene expression datasets, including GSE3678, GSE6004, GSE33630, GSE53157, and GSE58545, were downloaded from NCBI-GEO (https://www.ncbi.nlm.nih.gov/gds/). Among them, the GSE3678, GSE6004, GSE33630 and GSE53157 datasets were based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2 Array platform. The GSE58545 dataset was analyzed on the GPL96 [HG-U133A] platform, which was used to test the correlation of screened hub genes. The raw data from these datasets were processed using the R language statistical software (version 3.6.1). The horizontal linear model, affy (13), and affyPLM packages (14) were used for the probe horizontal. Box plots of RNA degradation and Normalized Unscaled Standard Error (NUSE) were plotted to test the trend consistency of the data after the quality of the data was analyzed. RNA degradation was checked using the AffyRNAdeg function and microarray datasets with consistent trends and good RNA quality were selected for follow-up analyses. Data from The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/)) was used as an external validation dataset.

To ensure comparability and completeness of the dataset, gene expression profiling data were standardized and batch effects were removed. Normalization was performed using the RMA function in the Affy package (13). Differentially expressed genes (DEGs) between the tumor and peri-tumor tissues were screened using the LIMMA package (15). Threshold values of log₂FC>1 and adjusted P<0.05 were considered statistically significantly meaningful and the dataset (GSE3678, GSE6004, GSE33630, and GSE53157) DEGs were subjected to overlapping analyses using the website (http://bioinformatics.psb.ugent.be/beg/).

GO was used to make predictions based on the molecular function (MF), cellular components (CC), and biological processes (BP) of potential functions of target genes. KEGG enrichment analysis was used to determine the functional properties of DEGs. An online annotation, visualization, and integrated discovery database (16) (DAVID; http://david.abcc.ncifcrf.gov/) and Metascape (17) (http://metascape.org) were used to access GO and KEGG pathway enrichment analyses. P<0.05 was deemed statistically significant.

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database search tool is an online tool used to evaluate PPI information (18). In this study, DEGs with co-expression coefficients greater than 0.4 were extracted from the STRING database. The Cytoscape software, a publicly available bioinformatics software platform for visualizing and analyzing molecular interaction networks, was used to build a PPI network (19). Hub genes were selected from the intersection of the top 65 genes calculated using 12 topological analysis methods by using the Cytoscape plugin cytoHubba (20).

The cBioPortal for Cancer Genomics (http://www.cbioportal.org/) is an open source platform that supports visualization, analysis, and downloads of multiple cancer genomic datasets (21). The frequency of genetic alterations in 13 hub genes in patients with PTC and their relationship to survival outcomes were explored by using cBioPortal.

Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/) is a newly developed interactive web server for analyzing the RNA sequencing expression data of 9,736 tumors and 8,587 peri-tumor samples from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) databases using a standard processing pipeline (22). Disease-free survival (DFS) is the duration of time from the start of treatment to the date of progression, the start date of second-line treatment, or the date of death, whichever occurs first. OS is the length of time from the date of diagnosis or first treatment to the date of death or the last follow-up visit. We used GEPIA to compare gene expression and to evaluate the contribution of the screened candidate genes to OS and DFS of PTC patients.

Tumor and paired peri-tumor samples were collected from 10 histopathologically and clinically diagnosed PTC patients at Tianjin Medical University General Hospital from January 2017 to March 2017. Total RNA was extracted using TRIzol reagent and the quality of the extracted RNA was checked by Agilent 2100 Bioanalyzer (Santa Clara, CA, USA). Then, the cDNA libraries were constructed and sequenced on a HiSeq 2000 at Beijing Genomics institution. The study was planned according to the ethical guidelines following the Declaration of Helsinki. This project was approved by the Ethics Committee of Tianjin Medical University General Hospital. All participants gave informed consent prior to enrollment in the investigation.

UALCAN (http://ualcan.path.uab.edu) is a comprehensive, user-friendly, interactive web resource using TCGA RNA-seq data and the latest clinical data for 31 cancer types (23). It was designed to allow for the analysis of tumor and peri-tumor samples as well as different tumor types based on individual cancer stages, tumor grades, or other clinicopathological features and the relative gene expression of the subgroups. We selected key genes from the 13 hub genes to correlate to multiple clinical disease characteristics.

The Human Protein Atlas (https://www.proteinatlas.org/) is a website containing immunohistochemistry (IHC) data to help investigate protein expression in human tissues and cells (24). Patient information, staining intensity, staining location, and sample number are available for each type of cancer. We assessed the expression of representative proteins in PTC and peri-tumor tissues of six key genes using IHC images after screening.

To assess the prognostic value of the hub genes, we downloaded the THCA clinical dataset from TCGA database and tested the six gene signature via multifactorial Cox regression analysis using the Sangerbox (http://sangerbox.com/Tool) online tool.

DNA Methylation Interactive Visualization Database (DNMIVD); (http://119.3.41.228/dnmivd/index/) is a comprehensive annotation and interactive visualization database for DNA methylation profiles of diverse human cancers constructed with high throughput microarray data from TCGA and the GEO databases (25). To determine the best-fitting prognostic model, we established a multivariate proportional hazard regression model for the six hub genes. We utilized a Kaplan-Meier multivariate proportional hazards regression model to divide patients into high-risk and low-risk groups.

Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/) is a web resource for systematic evaluation of the clinical impact of various immune cell populations on different cancer types, and consists of 10,897 samples from 32 cancer types. In this study, we analyzed the expression of six hub genes in THCA involved in tumor purity and the abundance of their immune infiltration (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells (DCs). Furthermore, we explored the relationship between gene copy number variation and the abundance of immune infiltration.

GSEA is the process of sorting genes by the degree of differential expression between two types of samples using a predefined set of genes, usually from functional annotations or results of previous experiments, and then verifying whether the predefined set of genes is enriched at the top or bottom of this sorting table. We downloaded the THCA gene expression matrix from TCGA database and performed GSEA analysis using the Sangerbox website (http://sangerbox.com/Tool) to predict potential hallmarks. A permutation test with 1,000 permutations was used to identify the significantly changed pathways. Thresholds of an adjusted p-value <0.05 were considered statistically significant.

Data were expressed as mean ± standard deviation. Analyses were done using GraphPad Prism (version 8.0; San Diego, CA, USA). Comparisons between the tumor and peri-tumor groups were performed using Student’s t-test. P value less than 0.05 were considered statistically significant.

Workflows are presented for DEGs identification, function enrichment, and the validation analysis ( Figure 1 ). To ensure more reliable quality of the datasets, GSE3678, GSE6004, GSE33630, and GSE53157 dataset deviations with those with higher levels were removed, leaving 76 PTC samples and peri-tumor samples for comparison. A total of 1,547, 1,050, 1,025, and 699 DEGs were identified for the GSE3678, GSE6004, GSE33630, and GSE53157 datasets, respectively ( Figures 2A–D ). Among the DEGs, 82 genes showed the same expression trend across the four datasets, including 55 upregulated genes ( Figure 2E ) and 27 downregulated genes ( Figure 2F ). Using the data profile GSE3678 as a reference, the heat map showed the distribution of significant differences among the 82 DEGs ( Figure 2G ). After filtering out 332 overlapping DEGs, 164 upregulated and 168 downregulated genes were discovered at the intersection of at least three microarray datasets for PTC and peri-tumor for GO and KEGG analyses. The data were visualized with Sangerbox software. The top 10 for each GO term and the top 10 KEGG pathways are listed. For BPs, DEGs were mainly enriched in blood coagulation, the bone morphogenetic protein signaling pathway, angiogenesis, wound healing, and extracellular matrix organization ( Figure 3A ). The CCs analysis revealed that DEGs were mostly enriched in the extracellular space, proteinaceous extracellular matrix, extracellular region, extracellular exosome, and costamere ( Figure 3B ). In molecular function, DEGs were primarily enriched for heparin binding, calcium ion binding, insulin-like growth factor binding, growth factor activity, and small molecule binding ( Figure 3C ). The KEGG pathway analysis indicated that the DEGs were predominantly enriched in the complement and coagulation cascades, extracellular matrix-receptor interaction, tyrosine metabolism, mineral absorption, proteoglycans in cancer, tight junction, cell adhesion molecules, p53 signaling pathway, PPAR signaling pathway, and oxytocin signaling pathway ( Figure 3D ).

332 DEGs previously screened were used in the PPI network analysis and a total of 242 DEGs (118 upregulated genes and 124 downregulated genes) were filtered into a PPI network containing 242 nodes and 466 edges using the STRING online database and were visualized using Cytoscape software ( Figure 4 ). Thirteen hub genes were identified at the intersection of the top 65 DEGs calculated with 12 cytoHubba algorithms, including CTGF, CDH3, CHRDL1, FGF13, LPAR1, OGN, TIMP1, CYR61, CHI3L1, SERPINA1, PROM1, CCL21, and WFS1 ( Figure 5A ). The correlation between the 13 hub genes were calculated using the TCGA-THCA dataset ( Figure 5B ). A Pearson’s correlation analysis revealed significant correlations among them. To validate these results, another GEO dataset (GSE58545) was used to measure the mRNA expression of the 13 hub genes. Rank clustering showed that the hub genes can distinguish PTC samples from peri-tumor group samples ( Figure 5C ). Furthermore, Metascape software was used for GO and KEGG pathway analysis. The top 13 GO enrichment terms( Figures 6A, B ) and BPs included: cartilage development, positive regulation of protein kinase activity, regulation of MAPK cascade, multicellular organismal homeostasis, negative regulation of cell projection organization, positive regulation of MAPK cascade, response to 1-oleoyl-sn-glycer, negative regulation of ATF6-mediated unfolded protein response, acute-phase response. CCs include: extracellular matrix, and prominosome. Molecular function includes: growth factor activity, and extracellular matrix structural constituent. The top 4 KEGG pathways were Rap1 signaling, melanoma, complement and coagulation cascades, and the NF-kappa B signaling pathway ( Figure 6C ). Most of the pathways are involved in tumorigenesis and immune processes.

In order to analyze the genetic alterations of the 13 hub genes in THCA in the cBioPortal online databases, the THCA (TCGA; Firehose Legacy), PTC (TCGA; Cell 2014), and THCA (TCGA, PanCancer Atlas) databases were used. Results showed that the percentages of genetic alterations in each dataset were 3.17% (16/505), 3.02% (15/496), 2.6% (13/500), respectively ( Figure 7A ). Specific to PTC, the alteration frequency of the 13 hub genes was 4% for CDH3, 3% for CCL21, 2.6% for CHRDL1, 3% for FGF13, 5% for LPAR1, 4% for OGN, 2.1% for TIMP1, 6% for CYR61 (CCN1), 4% for CHI3L1, 4% for SERPINA1, 2.3% for PROM1, 6% for CTGF (CCN2), and 3% for WFS1. In addition, the predominant alterations were expression changes instead of mutations ( Figure 7B ).

We also performed an analysis of the correlation between cases with hub gene alterations and survival outcomes. Accordingly, results showed that cases with hub gene alterations were significantly associated with worse OS and DFS (p=1.146E-3; p=0.0420; Figures 7C, D ).

We used GEPIA to evaluate the contribution of the 13 hub genes to clinical outcomes. The results revealed that CDH3 (p=0.017) was significantly associated with OS ( Figure 8A ). CYR61 (p=0.0017), OGN (p=0.039), CHRDL1 (p=0.0028), and FGF13 (p=0.047) were significantly correlated with DFS ( Figures 8D–G ). Moreover, CTGF was involved in both OS (p=0.022) and DFS (p=0.0004; Figures 8B, C ). The remaining hub gene survival analyses did not show any statistical significance ( Supplementary Figure 1 ). In summary, we demonstrated that the 13 hub genes are involved in PTC, and six genes (CDH3, CYR61, OGN, CHRDL1, FGF13, and CTGF) were significantly correlated to patient outcomes.

With GEPIA, we detected the transcript levels of six hub genes in 59 peri-tumor tissues and 512 THCA tissues from TCGA database. CDH3 was found to be highly expressed, while CTGF, CYR61, CHRDL1, OGN and FGF13 had significantly lower expression in tumor tissues compared with peri-tumor tissues ( Figure 9A ). In order to validate the bioinformatics analysis results, 10 paired PTC tumor and peri-tumor samples were collected and performed for transcriptome sequencing studies. Comparing with peri-tumor controls, the expression of CTGF, CYR61, CHRDL1, OGN and FGF13 were significantly decreased, whereas CDH3 significantly increased in PTC tissues ( Figure 9B ) which were consistent with the bioinformatics results obtained by TCGA dataset. Additionally, the expression of the six hub genes were analyzed according to various clinicopathological characteristics, which included tumor stage, lymph node metastases, sex, and age by using UALCAN software. CDH3 transcript levels were significantly increased, while the expression of CHRDL1, CTGF, CYR61, FGF13, and OGN were significantly decreased in THCA patients compared with healthy subjects ( Figures 10A–F ). The expression levels of the six genes were not different when comparing tumor stages, lymph node metastasis, age, and sex ( Figures 10A–F ). To verify the transcriptome analysis results, immunohistochemistry data of the six genes were obtained from the Human Protein Atlas website. The protein expression levels of CDH3, CHRDL1, CTGF, CYR61, and OGN showed similar patterns of changes as the transcript levels ( Figure 11 ).

To validate the classification performance of the six gene signature model using different data platforms, we used TCGA data as an external dataset. We then calculated the risk score for each sample using the Sangerbox website and used the cut-off values of the training set to classify the samples into high and low risk groups. The low risk group had a significantly different prognosis than the high risk group ( Figure 12A ). An ROC analysis showed 1-, 3-, and 5-year AUCs of 0.81, 0.67, and 0.72 for the six gene signatures, respectively ( Figure 12B ). The data were presented for TCGA samples focusing on risk score, survival time, survival status, and expression of the six gene signature ( Figure 12C ). CDH3 was highly expressed in tumor samples, whereas high expression of CTGF, CYR61, OGN, FGF13, and CHRDL1 was associated with better prognosis and appeared to be a protective factor. Furthermore, we performed methylation survival analysis of the six gene signature using the DNMIVD tool. Patients were then classified into either high risk or low risk groups, using the median risk score as the cut-off point. Patients in the high risk group had a significantly shorter median OS than those in the low risk group (p<0.022; Figure 12D ). We also calculated the Z-score distribution of the prognostic classifier and patient survival status ( Figure 12E ). Overall, our data showed that all six genes were significantly associated with THCA and could be used as a signature to sensitively and accurately determine prognosis in tumor patients.

The tumor microenvironment (TME) is composed of tumor cells, stromal cells, and infiltrating immune cells. The TIMER database was utilized to analyze the association of specific gene expression levels with THCA immune populations. Interestingly, CDH3 (r=-0.101, p=2.59E-02) and OGN (r=-0.104, p=2.09E-02) was negatively correlated with tumor purity. Moreover, the expression of CDH3, CTGF, OGN, and CHRDL1 was significantly correlated with the infiltration of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and DCs. The expression of CYR61 was significantly correlated with the infiltration of CD8+ T cells, CD4⁺ T cells, and neutrophils. The expression of FGF13 was significantly correlated with the infiltration of B cells, CD4+ T cells, macrophages, neutrophils and DCs ( Figure 13A ). Interestingly, the copy number variation (CNV) of FGF13 and CHRDL1 were significantly correlated with infiltration levels of B cells and DCs. OGN was associated with B cells and CYR61 was associated with CD4+ T cells. CTGF and CDH3 were both associated with B cells, neutrophils, and DCs, and the CNV of CDH3 was also significantly correlated with infiltration levels of CD8+ T cells, CD4+ T cells, and macrophages ( Figure 13B ).

To further explore the potential function of the six key genes in THCA, we performed GSEA on TCGA-THCA RNA-seq data. A total of 100 significant genes were acquired by GSEA with both positive and negative correlations. GSEA was used to conduct a hallmark analysis for CDH3, CTGF, CYR61, FGF13, OGN, and CHRDL1. FGF13, CTGF, CHRDL1, OGN, and CYR61 were all enriched through hypoxia and UV response downregulation pathways. Meanwhile, CTGF, CHRDL1, and OGN were all enriched in UV response upregulation pathways. The five genes (FGF13, CTGF, CHRDL1, OGN, and CYR61) were also enriched through pathways such as TGF-β signaling, adipogenesis, hedgehog signaling, androgen response, apical surface, apoptosis, notch signaling, KRAS signaling upregulated and down regulated genes, and early and late estrogen response. The signaling pathways most involved with CDH3 include the P53 pathway, apical junction, apoptosis, interferon-alpha response, coagulation, and glycolysis ( Figures 14A–F ). Some of these pathways are associated with tumor development. Remarkably, these pathways were significantly enriched in high risk samples.