The RNA-seq data (workflow type, HTSeq-FPKM) and relevant clinical data for the ‘TCGA-HNSC’ cohort, including 502 tumor samples and 44 normal samples, were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/tcga/) (26).

Expression levels of CTHRC1 were analyzed in multiple cancer types using the ONCOMINE database (https://www.oncomine.org/resource/main.html) (27), and the cut-off P-value used was 0.05, while the log fold change cut-off was equal to 1.

The TIMER database (https://cistrome.shinyapps.io/timer/) (28) has incorporated 39 types of cancer in the TCGA database. Gene expression levels were presented as log2 transcripts per million. The TIMER database was used to evaluate the transcriptional level of CTHRC1 in multiple cancer types and its correlation with immune cell infiltration.

KM plotter contains data and clinical information from the Gene Expression Omnibus, TCGA and European Genome-Phenome Archive databases. The prognostic value of the mRNA expression levels of CTHRC1 in 21 cancer types was evaluated using the KM plotter (http://kmplot.com/analysis/). P<0.05 was considered to indicate a statistically significant difference.

UALCAN (http://ualcan.path.uab.edu) (29) contains level 3 RNA-seq data and clinical information from 31 cancer types from the TCGA database and is an interactive and comprehensive web resource for analyzing cancer omics data. In the present study, UALCAN was used to investigate the relationship between the levels of CTHRC1 and clinicopathologic parameters, including age, TP53 mutation, nodal metastasis (N0, no regional lymph node metastasis; N1, metastases in 1–3 axillary lymph nodes; N2, metastases in 4–9 axillary lymph nodes; N3, metastases in ≥10 axillary lymph nodes), individual cancer stages, tumor grades (grade 1, well differentiated; grade 2, moderately grade; grade 3, poorly differentiated; grade 4, undifferentiated) (30,31) and HPV status. Unpaired student's t-test was used to assess transcriptional expression and P<0.05 was considered to indicate a statistically significant difference.

Based on the validated leukocyte gene signature matrix (LM22), the CIBERSORT (http://cibersort.stanford.edu/) computational method was used to analyze the infiltration ratio of 22 TICs in each HNSCC sample, and the ratios of all were equal to 1.

R 64.4.0.0 software (https://www.r-project.org) was used to plot associations with TICs. The Barplot and corrplot were generated using the corrplot package (version 0.84, http://cran.r-project.org/src/contrib/Archive/corrplot/), while vioplot was produced using the BiocManager (version 3.12, http://www.bioconductor.org/install) and vioplot packages (version 0.3.4, http://cran.r-project.org/src/contrib/Archive/vioplot/). The Venn plot was generated using the Venn Diagram package (http://www.rdocumentation.org/packages/VennDiagram/versions/1.6.18). The scatter plot was generated using the ggplot2 (version 3.2.1, http://cran.r-project.org/src/contrib/Archive/ggplot2/), ggpubr (version 0.2.4, http://cran.r-project.org/src/contrib/Archive/ggpubr/) and ggExtra packages (version 0.9, http://cran.r-project.org/src/contrib/Archive/ggExtra/). All packages were used with standard settings.

GSEA were downloaded from the Broad Institute Website (version 4.0.2, http://software.broadinstitute.org/gsea/index.jsp). HALLMARK, Kyoto Encyclopedia of Genes and Genomes (KEGG) and IMMUNE SIGNATURE gene sets were obtained from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). The transcriptome of patients with HNSCC were assessed using GSEA software (version 4.0.2, http://software.broadinstitute.org/gsea/downloads.jsp) and all the samples with nominal (NOM) P<0.01 and false discovery rate (FDR) Q<0.06 were considered to indicate a statistically significant difference.

From January 2020 to December 2020, 70 patients with primary HNSCC who underwent radical resection of tumors, did not receive preoperative radiotherapy or chemotherapy in the People's Hospital of Tongxu County (Kaifeng, China) and provided written informed consent were enrolled in the present study, and their tumor tissue samples were collected. A total of 11 patients with benign head and neck lesions were included, informed consent was signed and normal tissue samples adjacent to benign lesions were collected. Patients with HNSCC were between 32 and 79 years old, with a mean age of 59.9 years. Human tissue specimens were obtained from patients who had provided written informed consent and the present study and tissue collection were approved by the Ethics Supervision Committee of The People's Hospital of Tongxu County (Kaifeng, China; approval no. TX20NP003) and the National Human Genetic Resources Sharing Service Platform (2005DKA21300). The study protocol was approved by the Ethics Committee of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology (Wuhan, China; approval no. 2020IEC-J050).

The tissue samples were fixed using 10% formalin solution at room temperature for 24 h, dehydrated, embedded in paraffin and sectioned into 4 µm sections. The expression of CTHRC1 protein was assessed using the immunohistochemistry ultrasensitive TMS-P method. After the sections were deparaffinized, hydrated and washed, they were placed in EDTA repair solution (pH 9.0), and boiled using a microwave oven for 2 min, and then cooled to room temperature naturally. Then sections were blocked using 3% H₂O₂ solution for 15 min at room temperature in the dark room. The sections were washed and incubated overnight at 4°C with CTHRC1 antibodies (1:200; cat no. 16534-1-AP; Wuhan Sanying Biotechnology). After washing three times with phosphate buffered saline, the sections were stained with secondary antibodies using the Ready-to-use Ultrasensitive^™ SP kit (cat. no. Kit-9720; Fuzhou Maxin Biotechnology Development Co., Ltd.) and incubated for 30 min in a wet box at 37°C. DAB horseradish peroxidase color development kit (cat no. P0202; Beyotime Institute of Biotechnology) was used for color development. The sections were re-stained using hematoxylin staining solution (cat no. C0107; Beyotime Institute of Biotechnology) for 2 min at room temperature, dehydrated using ethanol, sealed and imaged using a Pannoramic MIDI scanner (3DHISTECH, Ltd.) and ImageJ software (version 1.8.0, National Institutes of Health) was used for image analysis. The cells on the microarrays were scored according to the staining intensity of the marker, as follows: 0, no coloring; 1, light yellow; 2, brown-yellow; and 3, tan. The proportion of cells under each score was then recorded, and the H-score value of each specimen was calculated. H-scores were defined as: (1× percentage of cells staining at 1) + (2× percentage of cells staining at 2) + (3× percentage of cells staining at 3).

SPSS 26.0 software (IBM Corp.) was used for statistical analyses. Fisher's exact test was used to analyze the relationship between clinical characteristics and CTHRC1 mRNA expression levels. The Mann-Whitney U test was used for calculating statistical differences in the H-scores in TMA and the ratio differentiation level of 21 types of immune cells in high and low CTHRC1 expression groups in HNSCC. Pearson correlation coefficient analysis was used to assess the correlation value between two kinds of cells. P<0.05 was considered to indicate a statistically significant difference.

To evaluate the distinct prognostic and potential therapeutic value of CTHRC1 in HNSCC, the mRNA expression level of CTHRC1 was studied using the TCGA, ONCOMINE and TIMER databases. Data from the TCGA database showed that CTHRC1 was highly expressed in 16 cancer types, including HNSCC (Fig. 1A and B). Results from the ONCOMINE database also showed that the mRNA expression levels of CTHRC1 were significantly upregulated in multiple tumor tissues compared with normal tissues (Fig. 1C). Furthermore, the same result was demonstrated by HNSCC clinical tissues samples when compared with normal tissues (Fig. 1D and E). Taken together, these data indicated that CTHRC1 was enriched in multiple human cancer types.

Since clinical pathology can determine the progression and prognosis of diseases, the association of transcriptional levels of CTHRC1 with clinicopathological parameters in patients with HNSCC was investigated using UALCAN. As presented in Fig. 2 and Table I, the transcriptional level of CTHRC1 was significantly associated with age (Fig. 2A), TP53 mutation (Fig. 2B), HPV status (Fig. 2C), tumor grade (Fig. 2D), nodal metastasis status (Fig. 2E) and individual cancer stages (Fig. 2F). The mRNA expression level of CTHRC1 increased markedly with age. Patients with TP53 mutations and those who were HPV-negative also had higher CTHRC1 mRNA expression levels. CTHRC1 mRNA expression levels increased with tumor progression and, in general, the more lymph node metastases in patients, the higher the mRNA expression level of CTHRC1. Among the patients, the transcriptional level of CTHRC1 in the N2 group (4–9 axillary lymph node metastasis) was lower than that of the N3 group (≥10 axillary lymph node metastasis), which may be due to the limited sample size. Additionally, the mRNA expression level of CTHRC1 in pathological stage IV patients was significantly higher than those in stage III or II, but markedly lower than that in patients with stage I, which may also be due to the difference in sample size. In summary, these results suggested that the mRNA expression level of CTHRC1 was closely related to the clinicopathological parameters of patients with HNSCC.

Since CTHRC1 was differentially expressed in diverse cancer types, the relationship between CTHRC1 expression levels and the survival time of patients with HNSCC was next evaluated using the KM plotter. Patients were grouped, into CTHRC1 high and CTHRC1 low expression groups, using the auto select best cutoff function on the KM plotter website (32). The data indicated that the CTHRC1 high expression group had a worse OS (Fig. 3A) and RFS time (Fig. 3B). These results demonstrated that CTHRC1 had the potential to be a prognostic biomarker in HNSCC.

In order to further investigate the potential function of CTHRC1, GSEA of the data from patients with HNSCC was conducted. A total of eight HALLMARK gene sets (Fig. 4A), 14 KEGG gene sets (Fig. 4B) and 1,231 immune signature gene sets (Fig. 4C) were significantly enriched in the high CTHRC1 expression groups (NOM P<0.01; FDR Q<0.06). These gene sets included angiogenesis, apical junction, coagulation, epithelial mesenchymal transformation (EMT), the KRAS signaling pathway, the notch signaling pathway, glycosaminoglycan biosynthesis, chondroitin sulfate and the TGF-β signaling pathway. Most of these gene sets serve pivotal roles in tumorigenesis. Therefore, CTHRC1 may mediate immune cell infiltration during cancer development, and it may be regarded as an important indicator for cancer progression.

The immune system serves pivotal roles in the tumorigenesis of HNSCC (3). The aforementioned data showed that CTHRC1 may be involved in immune responses and tumor associated pathways. Therefore, the association between CTHRC1 expression and immune cell invasion in patients with HNSCC from the TIMER database was evaluated (Fig. 5A). The results indicated that CTHRC1 expression was significantly correlated with the degree of infiltration of B cells, CD4⁺ T cells, neutrophils and dendritic cells (DCs) in HPV negative patients with HNSCC (Fig. 5B). However, no significant correlation was demonstrated for these groups in HPV positive patients with HNSCC. The expression of CTHRC1 was significantly associated with the infiltration of macrophages in both HPV negative and positive patients with HNSCC (Fig. 5C). This suggested that the relationship between CTHRC1 expression and macrophage infiltration was not related to HPV infection. Collectively, these results demonstrated that CTHRC1 expression was significantly associated with immune cell infiltration in HNSCC.

The infiltration of 22 TICs in each HNSCC sample and the CIBERSORT algorithm was used to further evaluate the correlation between the transcriptional levels of CTHRC1 and the immune microenvironment (Fig. 6A and B). Correlation analysis indicated that the transcriptional levels of CTHRC1 were significantly positively correlated with M0 macrophages (P=2.1×10⁻¹³; Fig. 7A), M2 macrophages (P=6.7×10⁻⁵; Fig. 7B) and resting CD4⁺ memory T cells (P=5.2×10⁻⁴; Fig. 7C), and were significantly negatively associated with the levels of activated CD4⁺ memory T cells (P=3.8×10⁻¹⁰; Fig. 7D), M1 macrophages (P=2.6×10⁻⁴; Fig. 7E), activated natural killer (NK) cells (P=2.1×10⁻⁴; Fig. 7F), monocytes (P=2.1×10⁻²; Fig. 7G), CD8⁺ T cells (P=5.9×10⁻⁹; Fig. 7H), activated DCs (P=4.2×10⁻⁵; Fig. 7I) and follicular helper T (Tfh) cells (P=3.0×10⁻³; Fig. 7J) (Table II). Analysis demonstrated that the transcriptional levels of CTHRC1 were correlated with the infiltration of 10 types of TICs, including resting CD4⁺ memory T cells, memory B cells, activated CD4⁺ memory T cells, CD8⁺ T cells, Tfh cells, M0 macrophages, activated NK cells, M1 macrophages, M2 macrophages and activated DCs (Fig. 7K and L). The high CTHRC1 expression group tended to have a larger proportion of resting CD4⁺ memory T cells, M0 macrophages and M2 macrophages compared with the low CTHRC1 expression group, whereas the low CTHRC1 expression group had a larger proportion of resting CD4⁺ memory T cells, memory B cells, activated CD4⁺ memory T cells, CD8⁺ T cells, Tfh cells, activated NK cells, M1 macrophages and activated DCs compared with the high CTHRC1 expression group.