Interferon Regulatory Factor Family Genes: At the Crossroads between Immunity and Head and Neck Squamous Carcinoma

Objective This study is aimed at investigating the regulating mechanisms of the interferon regulatory factor (IRF) family genes in head and neck squamous cell carcinoma. Methods Based on the HNSC data in the ‘The Cancer Genome Atlas (TCGA)' database, the expression pattern of IRF family genes was investigated. The association of IRFs family genes and survival outcomes were analyzed by Kaplan–Meier plotter web portal. The relation of IRF genes and tumor stages was evaluated by using stage plots and based on GEPIA portal. 50 genes interacting with IRFs were identified using the NetworkAnalyst's protein-protein interaction (PPI) network construction tool. The top 200 correlated genes with similar expression patterns in HNSC were obtained by the similar gene detection module of GEPIA. Furthermore, functional enrichment analysis was performed to determine the biological functions enriched by the interacting and correlated genes. The potential implication of IRFs in tumor immunity was investigated in terms of tumor-infiltrating immune cells, a pair of immune checkpoint genes (CD274 and PDCD1), and ESTIMATE-Stromal-Immune score. Results The unpaired sample analysis shows that all of the IRF family genes were highly expressed in HNSC tumor samples compared to control samples. The survival analysis results showed that the overexpression of IRF1, IRF4, IRF5, IRF6, IRF8, and IRF9 was associated with better overall survival in HNSC, while the other IRFs genes (IRF2, IRF3. and IRF7) did not show prognostic values for overall survival outcome of HNSC. Four genes (STAT1, STAT2, FOXP3, and SPI1) were overlapping among 50 interacted genes in the PPI network and top 200 correlated genes identified by GEPIA. The 50 interacting genes in the PPI network and top 200 correlated genes were integrated into 246 genes. These 246 genes were found to be overrepresented in multiple KEGG pathways, e.g., Th17 cell differentiation, T cell receptor signaling pathway, cytokine-cytokine receptor interaction, natural killer (NK) cell-mediated cytotoxicity, FOXO signaling, PI3K-Akt signaling, and ErbB signaling. Most correlations between IRF gene members and TIICs were positive. The strongest positive correlation was identified between IRF8 and T cells (r = 0.849, p < 0.001). The majority of correlation between IRF family genes and ESTIMATE-Stromal-Immune score was found to be positive. The highest positive correlation was found to be between IRF8 and Immune score (r = 0.874, p = 1.09E − 158). Most correlations between IRFs and two immunoinhibitor genes (CD274 and PDCD1) were positive. IRF1 and PDCD1 were found to show the highest positive correlation (r = 0.764, p < 2.2e − 16). Conclusions The current analysis showed IRFs were differentially expressed in HNSC, indicated significant prognostic values, were involved in tumor immunity-related signaling pathways, and significantly regulated tumor-infiltrating immune cells. IRF family genes could be potential therapeutic biomarkers in targeting tumor immunity of head and neck cancer.


Introduction
Interferon regulatory factors (IRFs) are well-known to be transcription factors involved in the regulation of interferon genes [1]. The interferon regulatory factor (IRF) family in mammals consists of nine members: IRF1, IRF2, IRF3, IRF4, IRF5, IRF6, IRF7, IRF8, and IRF9 [2]. Interferon regulatory factors (IRFs) are required for the activation of innate immune responses in response to pathogens and the subsequent induction of adaptive immunological responses [3]. It is believed that dysregulation of IRF signaling has a key role in the pathogenesis of autoimmune diseases, due to its inappropriate regulating role in immune cell activation and differentiation [3]. Although it was originally believed that IRFs function largely in the immune system by contributing to innate immune response, recent discoveries based on increasing evidence have revealed that IRFs play key roles in the governing oncogenesis [4][5][6]. Several outstanding and detailed reviews [4][5][6] have examined the role of IRFs in immunological disorders and malignancies in detail. They concluded that IRFs are placed at the nexus of immunity and cancer, connecting the processes governing both.
Increasing evidence has shown the role of IRF family genes in activating antitumor immunity and inhibiting immunosuppression. IRF1 was found to regulate C-X-C motif chemokine 10 (CXCL10)/chemokine receptor 3 (CXCR3) signaling axis, thereby further activating antitumor immunity in hepatocellular carcinoma (HCC). IRF4 was found to negatively regulate the development of myeloid-derived suppressor cells (MDSC) and its immunosuppressive function in tumors [7]. In addition, IRF family genes were found to be implicated in the antitumor immunity by controlling the differentiation and maturation process of tumor-associated immune cells including dendritic cells, NK cells, B lymphocytes, and T lymphocytes. For instance, IRF8 was shown to suppress tumor progression by inducing the maturation and differentiation of antigenpresenting cells (APCs) (e.g., MФ, DCs, and B lymphocytes) [8]. IRF4 was found to regulate the apoptosis of B lymphocytes by targeting the Fas apoptosis inhibition molecule [9]. Based on all this evidence, IRF family genes should be regarded as key mediators in regulating tumor immunity and thus might provide potential strategies in cancer immunotherapy.
Squamous cell carcinoma of the head and neck (HNSC) is the sixth most prevalent type of cancer worldwide, with a 5-year survival rate of less than 50% [10]. Previous studies demonstrated that IRF family genes were aberrantly expressed and also identified as diagnostic and prognostic markers in various hematological malignancies such as colorectal cancer [11], breast cancer [12], pancreatic cancer [13], and gastric cancer [14]. However, the effect of IRF family genes remains largely unclear in HNSC. The roles of IRFs in HNSC have not yet been investigated by using a systematic approach such as bioinformatics analysis. To the best of the authors' knowledge, this is the first study to examine the prognostic value and regulatory role of IRF4 especially focusing on immunity involvement.
The link between IRF family gene mRNA expression and clinical features of HNSC patients was studied using genomic and clinical data from The Cancer Genome Atlas (TCGA) database. Furthermore, the significant value of IRF family genes in HNSC prognosis was evaluated. In addition, functional enrichment analysis was used to determine IRF4's putative biological functions in HNSC. More importantly, the relationship between IRF4 and immunity was evaluated from several significant aspects including TIICs (tumor-infiltrating immune cells), a pair of classic immune checkpoint genes comprising of a receptor and its ligand (PD1 and PDL1), as well as the tumor microenvironment. The current research on IRFs may give a theoretical foundation for their mechanistic roles in tumor formation and immunology, as well as guidelines for medication therapy selection. Figure 1 used a flowchart to show the study design of the current research. Firstly, the genetic mutation, mRNA expression, protein expression, relation to tumor stages, and prognostic values of IRF family genes in HNSC were investigated. Afterwards, the interacted genes and correlated genes of the IRF family genes were, respectively, identified. After integrating the interacted and correlated genes, this group of genes were used for the functional enrichment analysis. More importantly, tumor immunity involvement of IRFs in HNSC was analyzed by investigating three aspects including tumorinfiltrating immune cells, two classic immunoinhibitory genes, and ESTIMATE-Stromal-Immune score.  3 Disease Markers assessed. The Kaplan-Meier curves were used to determine the link between target gene mRNA expression levels and relapse survival rate (RFS) and overall survival (OS) rates in the HNSC tumor group. The Kaplan-Meier survival graphs were used to illustrate the results. Using an online tool, the hazard ratio (HR) and 95% confidence interval (CI) were calculated automatically. The mean standard deviation is used to express the values for each group. When the Log-rank test was applied, a p value of 0.05 was considered statistically significant.

The Relationship between IRF Family Genes and Tumor
Stages. The expression DIY module of the web platform Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn) was used to plot gene expression by tumor stage, based on the TCGA clinical annotation. The F test was used to determine the expression of IRFs in various stages of HSNC tumors. The pathological stage map was created by using GEPIA platform to compare the expression of a certain IRF family member in various stages of HSNC tissue. p < 0:05 was judged as statistically significant.
2.9. NetworkAnalyst Webtool Analysis. NetworkAnalyst (URL: https://www.networkanalyst.ca/) is a visual analytics platform for comprehensive gene expression profiling and meta-analysis. Firstly, the organism was specified to be H. sapiens (human), and Entrez ID of all IRF family genes were uploaded: 3659, 3660, 3661, 3662, 3663, 3664, 3665, 3394, and 10379. Secondly, generic protein-protein interactions (PPI) was selected to be constructed based on the STRING interactome database. The confidence score cutoff was set as 900, and all the interactions required experimental evidence. Afterwards, a network comprising of 59 nodes, 87 edges, and 9 seeds was built and viewed. 2.11. Venn Diagram Analysis. In the previous steps, the 50 genes interacting with IRFs were identified by building a PPI network by using NetworkAnalyst, and the top 200 genes highly correlated with the IRF gene signature were identified by using the GEPIA webtool. A Venn diagram webtool (URL: http://bioinformatics.psb.ugent.be/webtools/ Venn/) was used to identify the gene list which were overlapped between these two groups of genes. This tool is able to calculate the intersections of list of gene elements. The correlations between the overlapping genes and IRFs were identified and shown in a heatmap.
2.12. Metascape Analysis. Metascape database (URL: https:// metascape.org) is able to combine functional enrichment and gene annotation into a single integrated gateway, leveraging over 40 separate knowledge bases. The IRFsinteracted genes and IRFs-correlated genes were integrated and the union genes of these two groups were obtained. This gene list was uploaded in the Metascape webtool, and custom enrichment analysis was performed. Four functional sets (GO biological processes (3683), GO molecular functions (577), GO cellular components (391), and KEGG pathways (212)) were, respectively, selected.
2.13. Tumor-Infiltrating Immune Cell (TIIC) Analysis. Pearson correlation coefficient analysis was used to determine the relationship between IRF family genes and tumorinfiltrating immune cells in HNSC tumor samples. The statistical method used the ssGSEA algorithm from the GSVA package (version 1.34.0). The lollipop plot was used to illustrate the link between the expression of IRFs and 24 different types of TIICs in HNSC samples.
2.14. ESTIMATE-Stromal-Immune Score. An analysis was carried out to analyze the correlation between IRF family genes and tumor immune microenvironment (TIME).   IRF1  IRF2  IRF3  IRF4  IRF5  IRF6  IRF7  IRF8  IRF9 Paired samples Normal Tumor e expression levels log 2 (TPM + 1)  7 Disease Markers tumor-immune system interactions that incorporates a variety of data types. This web portal was used to investigate the relations between two typical immunoinhibitor genes (PDCD1 and CD274) and expression, copy number, methylation, or mutation of each specific IRF gene. The reason of selecting these two genes is based on the fact that programmed cell death-1 (PDCD1) and its ligand programmed cell death 1 ligand 1 (CD274/PDL1) have been established to implicated in the T cell-mediated suppression of antitumor immunity.     The results showed that the samples did not pass the normality test (p < 0:05); thus, the Mann-Whitney U test (also named as Wilcoxon rank sum test) was adopted for this investigation. Figure 3(a) shows that the mean expression levels of all IRF genes were significantly greater in HNSC tumor samples than in healthy control samples (p < 0:05).
The paired sample analysis was based on the 43 HNSC tumor samples and their adjacent 43 healthy control samples. The results revealed that the differences between the tumor and control groups did not meet the criteria of the normality test (p < 0:05); thus, the Wilcoxon signed rank test was used for the current study. Figure 3(b) shows that the expression level values of all IRF genes in the HSNC tumor sample group were greater than that in the healthy control sample group. There was no statistically significant difference in IRF4 and IRF8; however, a statistically significant difference was observed for the other genes.   Disease Markers indicated that there was a strong correlation among the IRF family genes. There was a significant positive correlation among most IRF genes. The negative correlation was observed in the following pairs: IRF5 and IRF7 and IRF6 and IRF8. Figure 4 shows that the protein level of IRF1, IRF2, and IRF8 were significantly downregulated in HNSC compared with healthy control samples. Except for these three genes, the protein level of the other IRF family genes (IRF3, IRF4, IRF5, IRF6, IRF7, and IRF9) were significantly upregulated in HNSC compared with healthy control samples.

The Expression of IRF Family Genes in Different Tumor
3.9. Identification of IRFs' Correlating Genes by GEPIA. By using the similar gene detection module of GEIPIA webtool, the top 200 genes with similar expression pattern with IRFs in HNSC were identified. Table S2 shows that the top 10 genes highly correlated with IRF family genes signature are as follows: IRF1, TIGIT, ICOS, IL2RB, SLA2, AKNA, FMNL1, CD2, FGD2, and ARHGAP30.

The Overlapping Genes between Interacted Genes and Correlated
Genes. The Venn diagram in Figure 9(a) shows that 4 genes were at the intersection between 50 interacted genes in the PPI network and 200 correlated genes identified by GEPIA. These four genes are as follows: STAT1, STAT2, FOXP3, and SPI1. The correlation analysis results in Figure 9(b) and Table 2 showed that the majority of gene pairings had a substantial positive correlation. Between IRF6 and SPI1, a negative association was detected (r = −0:148, p < 0:001). IRF8 and SPI1 had the strongest positive connection (r = 0:864, p < 0:001).

Identification of IRFs Involved in Biological Functions.
The 50 IRFs-interacting genes and 200 IRFs-correlated genes were integrated into 246 genes. Figure 10 shows that these 246 genes were implicated in immune cell-related biological processes (BPs) (e.g., leukocyte activation, regulation of lymphocyte activation, positive T cell selection, and negative regulation of leukocyte activation) and cytokine-related BPs (e.g., positive regulation of cytokine production, response to cytokine, and regulation of cytokine tumor necrosis factor (TNF) production). In addition, these 246 genes were significantly enriched in multiple molecular   Figure 7: Pathological stage plot demonstrating the expression of a specific member of the IRF family at various tumor stages in HNSC patients.

The Correlation between IRF Expression and Immune
Cells in HNSC. Figure 11(a) shows that most correlations between IRF gene members and TIICs were positive. The correlations between IRF6 and most TIICs (e.g., B cells, CD8 T cells, cytotoxic cells, dendritic cells, immature dendritic cells, macrophages, NK CD56dim cells, NK cells, plasmacytoid 3.13. The Correlation between IRF Family Genes and ESTIMATE-Stromal-Immune Score. Figure 11(b) and Table 3 show that majority of correlations were found to be positive. The highest positive correlation was found to be between IRF8 and Immune score (r = 0:874, p = 1:09E − 158 ). The highest negative correlation was found to be between IRF6 and ESTIMATE score (r = −0:236, p = 8:75E − 08).

Discussion
The main findings of the current research showed that IRF family genes were involved in several signaling pathways, e. g., FOXO signaling pathway, PI3K-Akt signaling pathway, and ErbB signaling pathway. IRF family genes were significantly implicated in tumor immunity by being mainly positively correlated with tumor-infiltrating immune cells and regulating tumor immunosuppression. Figure 9(b) shows that the majority of IRF family genes were positively correlated with their interacting and correlating four genes (STAT1, STAT2, FOXP3, and SPI1) in head and neck cancer. A significant role of the TAT1 and STAT2 proteins lies in an important role in interferon (IFN) signaling and cellular antiviral responses and adaptive immunity [9]. STAT1 and STAT2 proteins were found to associate with IRF9 to form a heterotrimeric transcription factor complex known as ISGF3 [15]. This should be the reason leading to the positive correlation between STAT1-2 and IRF9. The forkhead/winged-helix family transcriptional repressor FOXP3 is a surface marker specifically expressed in CD4 +CD25+ Treg cells and plays a pivotal role in regulating its development and differentiation [16]. A previous study found that IRF1 was found to negatively regulate the function of CD4+CD25+ Treg cells by directly and specifically repressing the expression of FOXP3 gene, suggesting the negative correlation between IRF1 and FOXP3 in tumor immunity [17]. However, the prediction of our computational biology analysis showed a positive correlation between IRF1 and FOXP3 in HNSC. The transcription factor SPI1, as a protein of ETS family, is a surface marker expressed in macrophages and neutrophils [18]. The present research showed that SPI1 was positively correlated with the majority of IRF family genes except for IRF6. IRF4 and IRF8 were found to cooperate with SPI1 and thereby upregulate the expression of pro-inflammatory genes (e.g., CD20, Ig light chain enhancers, IL-18, and IL-1β) [19][20][21]. SPI1-IRF coactivating complexes such as IRF4-SP1 and IRF8-SPI1 were found to bind to a SPI1-IRF composite element in the promoter of the proinflammatory cytokine IL-1β [22]. Based on the role of interleukin-1β (IL-1β) in modulating the tumor microenvironment and further promoting tumor growth, the SPI1-IRF coactivating complexes might play a role in contributing to the formation of an inflammatory tumor microenvironment.
The functional enrichment analysis showed that IRF family genes' 246 correlated and interacted genes were significantly involved in several KEGG signaling pathways including FOXO signaling pathway, PI3K-Akt signaling pathway, and ErbB signaling pathway. A previous conducted by Yuan et al. identified IRF7 to be an important target gene of FOXO3 by carrying out genome-wide location analysis and gene deletion experiments [23]. FOXO3 as a member of the forkhead family was shown to negatively regulate a subset of antiviral genes (e.g., Gbp2, Ccl5, Ifit1, Irf7, and Oasl1) [24]. A regulatory circuit formed by FOXO3, IRF7, and IFN-I was found to limit inflammatory consequences resulting from antiviral responses [25]. In addition, the relationship of IRF family genes and PI3K/Akt pathway has been demonstrated by previous evidence [24]. Dhamanage et al. discovered that activation of the PI3K/Akt signaling pathway is required for IRF-7 translocation from the cytoplasm to the nucleus in plasmacytoid dendritic cells [26,27]. The PI3K/Akt pathway was shown to diminish IRF-3- 13 Disease Markers dependent promoter activity and impede dimerization of IRF-3, thereby further resulting in the loss of host antiviral activity [26]. The ErbB family of proteins consists of four receptor tyrosine kinases (e.g., Her1 (EGFR, ErbB1), Her2 (Neu, ErbB2), Her3 (ErbB3), and Her4 (ErbB4)) that share structural homology with the epidermal growth factor receptor (EGFR) [27]. The upregulation of IRF1 was found to be induced by the activation of EGFR signaling was manifested in an EGF dose-dependent manner. The induced expression of tumor suppressor IRF1 was found to alert antitumor immunity by activating immune effector cells and inhibiting cell proliferation of human cancer cells [28].
In terms of IRF family genes' relationship to tumor immunity, the present research revealed that the IRF family    gene expression was associated with many types of immune cells, for example, NK cells, Treg cells, Th1 cells, macropahages, and Th17 cells. Our research showed that IRF3 was positively correlated with NK cells in HNSC (Figure 8(a)).
In accordance with such findings, a previous study conducted by Andersen et al. found that exosomes derived from head and neck cancer cells were able to increase the expression of interferon regulatory factor 3 (IRF3) and further promoted the function of NK cells with respect to proliferation, cytotoxicity, and the release of perforin and granzyme M [29,30]. Our research showed that IRF1 was significantly positively correlated with Treg cells and Th1 cells. In a related finding, a previous study conducted by Chen et al. showed that IRF1 exhibited immunomodulatory activities by controlling Treg depletion and governing Th1 polarization [31]. Additionally, it has been well established that a high infiltration of macrophages is associated with decreased overall survival and metastatic progression. Our research found that IRF8 was positively correlated with macrophages in HSNC. A study by Buccione et al. showed that IRF8 was able to transcriptionally modulate the response of macrophages by driving macropshages towards a more protumor and prometastasis phenotype [32]. Additionally, the present study found that IRF8 was positively correlated with Th17 in HNSC. However, a previous study using experimental colitis showed a negative correlation between IRF8 and Th17. This  Disease Markers study found that IRF8 is required for Th17 cell differentiation to be suppressed: IRF8 transduction decreased the expression of genes linked with Th17; and vice versa, the lack of IRF8 was able to enhance the immune response induced by Th17 [33]. The tumor microenvironment of HNSC was found to strongly induce Th17 cells through cytokines; and in another way, the presence of Th17 cells was able to promote the proliferation and angiogenesis of HNSCC [34]. Based on the promoting role of Th17 in the carcinogenesis of HNSCC as well as the promoting role of IRF8 in HNSCC, the correlation between Th17 and IRF8 should be speculated to be positive; however, it still needs to be verified in the future experimental investigations.
Regarding the involvement of IRF family genes in tumor immunosuppression, the present study found that IRF family genes were significantly correlated with the immune checkpoint genes (PDL1 and its receptor PDCD1) in HNSC. The correlation between PDL1 and interferon-γ (IFNγ) has been demonstrated by showing that IFN upregulated PDL1 in a variety of HNSC cell lines regardless of HPV infection. The high link between IFN and PDL1 is explained by the fact that an antitumor cellular immune response mediated by natural killer (NK) cells and CD8+ tumor-infiltrating lymphocytes (TILs) produces IFN, which in turn induces PDL1 expression on tumor cells [35]. However, the current research regarding the correlation between PDL1 and IRF family genes especially in HNSC is lacking. According to the authors' knowledge, there is a paucity of evidence showing the PDL1 and IRF family genes in other types of cancers. For instance, a previous study investigating lung cancer found that interferon-β (IFN-β) induced the upregulation of PDL1 in lung cancer cells via the activation of the IRF9 pathway [36]. Although a different tumor type, the current research also found a positive correlation between PDL1 and IRF9 in HNSC; however, such finding needs to be verified by future experimental research. Another research investigating hepatocellular carcinoma (HCC) also found that IRF-1 was able to upregulate the IFN-γ-induced PDL1 mRNA and protein expressions, indicating a positive correlation between IRF1 and PDL1 [37]. In accordance with these findings, the present study also showed that PDL1 was positively correlated with IRF1 in HNSC. Additionally, the current study showed that IRF3 was not associated with PDL1 in HNSC (p = 0:0928); however, a previous research found that IRF3 promoted the induction of PDL1 by forming a transcriptional complex with NF-κB/p65, and further resulting in the solar ultraviolet radiation-(UVR-) induced immune suppression [38].
In summary, the current research provided a comprehensive view regarding the implication of IRF family genes in HNSC from a variety of aspects, for example, expression pattern, prognostic values, relationship with clinical stages, interacting and correlating genes, enriched biological functions, involved signaling pathways, correlation with tumor cells, and classic immune checkpoint genes. The current research has identified directions towards a better understanding of mechanisms of IRF family genes in HNSC, and these findings based on computational prediction warrant experimental work for validation.

Conclusion
The upregulation of the IRF1, IRF4, IRF5, IRF6, IRF8, and IRF9 genes was identified as a potential prognostic indicator in HNSC. IRF family genes played a tumor-promoting role in HNSC by involving several pathways including Th17 cell differentiation, T cell receptor signaling pathway, cytokinecytokine receptor interaction, NK cell-mediated cytotoxicity, FOXO signaling, PI3K-Akt signaling, and ErbB signaling. Given such findings, the genes of the IRF family should be considered possible therapeutic targets in the treatment of head and neck cancer.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.