Effect of cancer‐associated fibroblasts on prognosis and immune infiltration in head and neck squamous cell carcinoma

Cancer‐associated fibroblasts (CAFs) are the center of cross‐communication between various cells in the tumor stroma. However, how CAFs‐associated genes play an important role in Head and neck squamous cell carcinoma (HNSCC) prognosis has not been reported. Transcriptome data were downloaded from TCGA and GEO databases. Devtools, DPIC, xCell, MCPcounter, and Estimate packages were used to calculate CAFs scores and immune infiltration. Prognosis and weighted gene coexpression network analysis (WGCNA) analysis were performed between high or low risk populations based on CAF scores. Hub genes were identified, intersected, and enriched between TCGA and GEO databases. CAFs related genes were used to construct a prognostic model and the tumor immune dysfunction and exclusion database was used to evaluate the immune infiltration. Drug sensitivity, difference analysis and the HPA database were used to identify sensitive drugs and verify their expression. TCGA and GEO data suggested that CAFs scores had a role in HNSCC prognosis prediction. Based on CAFs scores, WGCNA and core gene enrichment analysis were performed to construct a CAFs‐related prognostic model. The prognostic model composed of a total of 12 CAFs genes could predict the prognosis well and was validated in the validation dataset, demonstrating its applicability to external data. According to the model, although there was no statistical difference in immune escape between the high and low risk groups, the proportion of patients who responded to immunotherapy was different. Drug sensitivity also differed between the two groups. This study suggests that CAFs associated genetic signatures may help to optimize risk stratification and provide new insights into individualized cancer treatment.

recovery.HNSCC is sensitive to radiotherapy, but it can lead to more severe mucosal damage, which is difficult for patients to tolerate and interrupt the treatment.Patients with HNSCC may exhibit different clinicopathological features, suggesting that traditional clinicopathological staging may not be completely accurate.Therefore, the search for new prognostic biomarkers and the identification of treatmentsensitive and high-risk populations have certain clinical reference significance for the treatment of HNSCC patients.
Cancer-associated fibroblasts (CAFs) are stromal cell populations with heterogeneity in cell origin, phenotype and function, and are a major component of the tumor microenvironment (TME).Their role in the progression of HNSCC is gradually being appreciated.CAFs promote tumor proliferation, angiogenesis, invasion, immune escape, drug resistance, and metastasis. 1 CAFs contribute to the growth of HNSCC through the secretion of proteins, miRNAs, and metabolites, and also secrete factors that act on the tumor cells, leading to a more aggressive phenotype. 2These are the theoretical basis for CAFs affecting tumors and microenvironment.However, the association between the genetic signature of CAFs and their prognostic and therapeutic predictive roles in HNSCC remains under-explored.In this study, we explored the correlation between CAFs and prognostic based on public databases and constructed relevant prognostic features, then validated them in a validation cohort.The aim of this study was to explore the clinical significance and prognostic impact of CAFs-related genes in HNSCC, and to construct a prognostic prediction model.The relationship between target genes and immune cells and drug sensitivity was investigated through a database.These results provide theoretical guidance for the formulation of therapeutic strategies for HNSCC.

| Transcriptome data download and CAF and immune scores
The HNSCC transcriptome data were downloaded from TCGA (https:// cancergenome.nih.gov/) and GEO GSE65858 (https://www.ncbi.nlm.nih.gov/) databases, respectively.The GEO contained the main clinical information and transcriptome data of 270 HNSCC patients (Table 1

| CAF score-based prognosis and weighted gene coexpression network analysis (WGCNA) analysis
Among the results of different immune infiltration scores, patients in TCGA and GEO database were divided into two groups of high and low risk according to the median score.The survival and survminer packages were analyzed for prognostic.
WGCNA identified gene modules with similar expression patterns. 3gulatory networks between genes in gene modules were mapped and the relationship between key regulatory genes and CAF score was identified.Analysis was performed with Limma, gplots, and WGCNA packages, and genes with high fluctuations were selected for analysis.

| Database intersection of core genes and enrichment analysis
The core genes obtained from the above TCGA and GEO were intersected to narrow down the screening range using VennDiagram package, and clusterProfiler, 4 enrichplot and ggplot2 were used for screening and enrichment analysis of the intersected core genes.

| Construction and evaluation of prognostic models
Circulation of genes was performed, and univariate Cox analysis was used to find prognosis-related genes.TCGA was used as the training

| Immune correlation analysis
Based on the CAF score and risk score of TIDE, the plyr, ggplot2, and

| Drug sensitivity and difference analysis
The oncoPredict and parallel packages were used to perform drug sensitivity analysis of gene expression in TCGA data.The two groups were set up according to the high and low risk scores, and the drugs were cycled and analyzed for differences, and box plots were drawn.

| Human protein atlas (HPA) database validation
The human proteome database HPA (https://www.proteinatlas.org) was used to validate the protein expression of model risk genes between tumors and normal tissues.The GEPIA database (http://gepia.cancer-pku.cn)was used to analyze the influence of single gene expression on prognosis.

| Immunoinfiltration score and prognosis analysis of CAFs
CAFs and immune infiltration scores were performed using TCGA and GSE65858, and survival differences were analyzed based on the

| WGCNA analysis and core gene enrichment analysis based on CAF scores
In order to further analyze the relationship between gene sets with similar expression patterns and CAF, WGCNA analysis was performed and key regulatory genes were screened out.In TCGA, the most

| Construction of CAF-related prognostic models
To explore the prognosis of these intersecting genes, we performed a univariate Cox regression analysis.The results showed that 19 of the 180 intersecting genes in TCGA had prognostic roles, of which 5 were protective genes and the remaining 14 were risk genes (Figure 3A).Lasso analysis was performed to construct the prognostic model (Figure 3B,C).In GO analysis, these genes were mainly enriched in cytoplasmic translation, ribosome, focal adhesion, cadherin binding, and cytokine receptor binding.In KEGG analysis, these genes were mainly enriched in COVID-19, ribosome, lysosome, reactive oxygen species, and phagosome (Figure 3B,C).After that, the importance of the selected prognostic genes was ranked by random forest tree.
The maximally selected rank statistics analyzed the risk scores of continuous variables and the relationship analysis after binary classification of them (Figure 3E) to calculate the optimal cutoff value.

| Prognostic and score correlation analysis of CAFs model
In

| CAF-related risk typing and immune function analysis
Next, we further investigated the role of CAFs-associated risk typing in immune function.The results of the TIDE database analysis did not reveal any significant difference in immunotherapy escape between the high-and low-risk groups, but the proportion of non-responder patients in the high-risk group was higher than in the low-risk group (Figure 5A,B).The AUC value of the CAFs-related risk model predicting immunotherapy efficacy was 0.571, which showed the predictive ability (Figure 5C).GSEA analysis was performed separately on the CAFs-related risk types in the high and low risk groups, and it was found that in the high-risk group, CYTOKINE_RECEPTOR_INTERAC-TION, FOCAL_ADHESION, HYPERTROPHIC_CARDIOMYOPATHY_HCM, PROTEASOME and RIBOSOME pathways were active, whereas in the low-risk group ARACHIDONIC_ACID_METABOLISM, DRUG_ METABOLISM_CYTOCHROME_P450, GLYCOSPHINGOLIPID_BIO-SYNTHESIS_LACTO, and IMMUNODEFICIENCY pathways were active (Figure 5D,E).

| Correlation analysis of tumor mutation load
In the correlation analysis of tumor mutation load, no correlation was found with the risk score (Figure 6A,B).The results of the correlation analysis between CAFs scores and tumor mutation load showed that all CAFs scores were correlated with TMB, except for the risk scores, which was consistent with previous studies (Figure 6C).

| Drug sensitivity analysis
To further investigate the response of these two groups of patients to the treatment, a drug sensitivity analysis was performed.Among the results analyzed with sensitivity differences, the following differential drugs are currently used in clinical or clinical trials (Figure 7A-O).
In addition to Trametinib, which was less sensitive in the high-risk group of patients, other differential clinical or clinical trial drugs were less sensitive in the low-risk group.

| Immunohistochemical expression analysis of risk genes in prognostic model
To The prognostic effect of these genes individually was verified in GEPIA, and it was found that OAF and ITGA5 had prognostic effects, demonstrating the prognostic role of these genes (Figure 8C,F,I).

| DISCUSSION
With the development of imaging, the popularity of health insurance policies, and the increased coverage area of medical institutions, the detection rate of head and neck tumors has gradually increased.However, the special location of head and neck tumors, complex anatomical The purpose of the WGCNA analysis method is to find coexpressed gene modules, explore the association between gene networks and CAF, find the core genes in the network, quickly locate target genes from a large number of genes, and take the intersection in different databases to obtain the genes with higher potency.7][18][19][20][21] In contrast, ENC1 has been classified as a tumor suppressor gene due to its loss of expression in many human cancers. 22merous studies have shown that CAF is intricately cross-linked with TME.As one of the most important stromal components of the TME, CAFs exhibit biological heterogeneity in many aspects, including cellular origin, phenotype, and function. 23CAFs also exert certain immunosuppressive effects by interacting with immune cells.For example, CAFs are able to promote the recruitment of inhibitory cells and limit the recruitment of effector cells in the immune microenvironment by secreting different cytokines. 24,25In this study, we demonstrated that HNSCC patients were divided into high-risk and lowrisk groups based on 12 CAF-related gene prognostic models, and found that although immune escape was not statistically significant, there were differences in response to immunotherapy between the two groups.Drug sensitivity analysis showed that the sensitivity of  This study suggests that CAF-associated genetic signatures may help to optimize risk stratification and provide new insights into individualized tumor therapy.However, out work exists some limitations.
For example, our experiment did not conduct biological functions and investigate the underline mechanisms, providing limited evidence for new means for clinical treatments.Moreover, this is a retrospective study, which needs prospective evaluations to enhance the robustness of our results.Finally, there are also many subtypes of head and neck tumors, which have not been classified and discussed due to insufficient cases.In the following days, we have to make our efforts on the research of a specific tumor to make our research more meaningful.

AUTHOR CONTRIBUTIONS
Yuanyuan Xu and Huanfeng Zhu designed the study.Dan Zong and Luxi Qian reviewed, analyzed and interpreted the data.Yi Cai and Nan Xiang developed the original concept, wrote and revised the manuscript.All authors approved the final manuscript.
ggpubr packages were used to calculate the difference in TIDE score and the difference in the proportion of response to immunotherapy in different risk groups.The enrichment pathways in both high-and lowrisk groups were analyzed by GESA.TCGA mutation data were downloaded, and the correlation between tumor mutation load and risk score was analyzed by the reshape2 package.The correlation between CAF score and tumor mutation load was analyzed by ggExtra package.

1
Immune infiltration scores and prognostic analysis.(A) TCGA database TIDE score prognosis; (B) TCGA database xCell score prognosis; (C) TCGA database stromal score prognosis; (D) GEO database EPIC score prognosis; (E) GEO database MCPcounter score prognosis; (F) GEO database Stromal score prognosis.scores.The results showed that TIDE, xCell, and stromal immune infiltration scores of three CAFs in the TCGA database could distinguish patients with different prognoses (Figure1A-C).In the GEO database, the immune infiltration scores of EPIC, MCPcounter, and stromal also predicted the prognosis of patients (Figure1D-F).This suggested CAFs played a role in the prognostic prediction of HNSCC.

F I G U R E 2
Weighted gene coexpression network analysis (WGCNA) analysis of head and neck tumors.(A,B) WGCNA analysis in TCGA database; (C,D) WGCNA analysis in GEO database; (E) core intersecting genes; (F) GO analysis in core intersecting gene; (G) KEGG analysis in core intersecting gene.significant difference was found in the MEblue module, while in GEO, the most significant difference was found in the MEpink module (Figure2A-D).By intersecting key genes of TCGA and GEO, a total of 180 intersected genes were obtained (Figure2E).According to the analysis by GO and KEGG, these genes were found to be mainly enriched in collagen-containing extracellular matrix, extracellular matrix structural constituent, protein digestion and absorption, PI3K-Akt signaling pathway, and focal adhesion (Figure2F,G).

F
I G U R E 2 (Continued) F I G U R E 3 Construction of prognostic models associated with cancer-associated fibroblasts.(A) Univariate Cox regression analysis of intersecting genes; (B,C) Lasso regression cross-validation model; (D) Random forest tree for importance ranking of features; (E) risk scores and distributions.
order to verify the prognostic role of CAFs-related model, patients were divided into two groups of high and low risk expression according to the cutoff values obtained above in the training set TCGA and the validation set GEO, and prognostic analysis were performed.There were 12 CAFs genes in the prognosis model, including HEY2, TSPAN9, OAF, ENC1, SRPX2, MME, EGFL6, COL8A2, EFEMP1, ITGA5, SRPX, and THBS1.We found that the model predicted prognosis well, and was validated in the validation set, demonstrating its applicability to external data (Figure4A,B).The correlation between the risk score of the prognostic model and CAF score was further analyzed.The results showed that the CAFs scores in xCell was negatively correlated with the risk scores, while the rest were positively correlated (Figure4C).We then analyzed the relationship between the CAFs-related genes in the risk model and 23 genes for CAFs reported in the literature and risk scores.The results showed that CAF-related genes were differentially expressed in both high and low risk groups, and there was a significant correlation between ASPN, F I G U R E 4 Survival analysis based on cancer-associated fibroblasts (CAFs) models.(A) CAFs model-related survival analysis in TCGA; (B) CAFs model-related survival analysis in GEO data; (C) CAFs score correlation analysis; (D) risk heat map; (E) the correlation between CAFs-related genes and risk scores.CAV1, FAP, FN1, FOXF1, MFAP5, MMP2, OGN, PDGFRA, PDPN, S100A4, SLC16A4, and TNC and risk scores (Figure 4D,E).

F I G U R E 5
Cancer-associated fibroblasts (CAF)-related risk typing and immune function analysis.(A,B) Immunotherapy analysis of the TIDE database; (C) accuracy of the ROC curve analysis model in predicting the effect of immunotherapy; (D,E) GSEA analysis of CAFs-related typing.
verify the prognostic model of gene expression in clinical samples, we searched the HPA database (https://www.proteinatlas.org/) for the expression of relevant proteins in normal nasopharyngeal and head and neck tissues.The results showed that within the retrievable range, MME was highly expressed in tumors (Figure 8A,B), OAF was lowly expressed in tumors (Figure 8D,E), and ITGA5 was moderately expressed in tumors (Figure 8G,H), as validated by clinical specimens.

F I G U R E 6 5 F
Construction and evaluation of prognostic models.(A,B) Correlation analysis of tumor mutation load; (C) Correlation analysis of cancer-associated fibroblasts score and tumor mutation load.structure, and highly lymph node metastatic characteristics pose certain challenges to treatment.Therefore, it is important to find effective biomarkers to diagnose and predict the prognosis of HNSCC patients.Epithelial-mesenchymal transition (EMT) is often reported to be associated with strong stromal infiltration, which further drives malignant tumor derivation and leads to poor prognosis.I G U R E 7 Analysis of differences in drug sensitivity.(A) Vorinostat; (B) Venetoclax; (C) Trametinib; (D) Temozolomide; (E) Tamoxifen; (F) Sorafenib; (G) Savolitinib; (H) Ribociclib; (I) Linsitinib; (J) Ipatasertib; (K) Crizotinib; (L) Axitinib; (M) Alpelisib; (N) AZD5991; (O) Afuresertib.Emerging evidence from tumor ecosystem studies suggests that CAFs are two major components of the tumor stroma and a key component of the TME, with multiple functions including stromal deposition and remodeling, extensive mutual signaling interactions with cancer cells,and crosstalk with infiltrating leukocytes.6It has been shown that CAFs can exert their immunosuppressive effects through the secretion of various effectors involved in the regulation of tumorinfiltrating immune cells as well as other immune components in the tumor immune microenvironment (TIME), thereby promoting tumorigenesis, progression, metastasis, and drug resistance.7 Studies on the prognosis and clinical stratification of CAFs in a variety of tumors have been reported in recent years, 8-10 however, whether CAFrelated genes or scores play an important role in the prognosis of HNSCC has not been reported.In this study, we identified and validated CAF score-related prognostic features and demonstrated that CAF treated HNSCC through its regulatory effect on immune cells for the assessment of prognostic and immunotherapeutic response.The reliability of the results was verified by survival analysis, multifactorial regression analysis, TCGA, GEO, and various public databases.The effect of the prognostic characteristics of this study on clinical decision-making in cancer patients is unclear, and further researches are needed to clarify it.
temozolomide, tamoxifen, and sorafenib, which were commonly used in clinical practice, was different in high-and low-risk patients.Undeniably, a deeper understanding of the multidimensional interactions between CAF and immune cells within the TME will help us to better define the mechanisms of CAF-induced immunosuppression, and further exploration of these interactions may identify more potential molecular targets for CAF-targeted therapy.

F I G U R E 8
The immunohistochemical expression of the prognostic model genes was verified by the human protein atlas database.(A) MME protein expression in nasopharyngeal normal tissue (expression level, low); (B) MME protein expression in Head and neck squamous cell carcinoma (HNSCC) tissue (expression level, high); (C) Survival curve based on MME from GEPIA database; (D) OAF protein expression in nasopharyngeal normal tissue (expression level, high); (E) OAR protein expression in HNSCC tissue (expression degree, low); (F) Survival curve based on OAR from GEPIA database; (G) ITGA5 protein expression in nasopharyngeal normal tissue (expression degree, low); (H) ITGA5 protein expression in HNSCC tissue (expression degree, medium); (I) survival curve based on ITGA5 from GEPIA database. ).
T A B L E 1 Clinical characteristics of TCGA and GEO datasets.
152][13][14] first reported as the causative gene of peroneal muscular atrophy in the population, but it is an oncogene in tumors, with a better prognosis for high expression in breast, bladder, and prostate cancers, and is a key regulator gene for esophageal cancer metastasis.11OAF,anothername for OA3, is an autosomal dominant ocular albinism gene, and its correlation with tumor prognosis and immunotherapy has not been reported.There are many studies on the effect of ITGA5 gene on tumor prognosis, including cervical cancer, glioma, and gastric cancer.[12][13][14]Fengetal.also verified the role of ITGA5 in head and neck tumors based on public databases.15The