Pan-Cancer Analyses Reveal Oncogenic and Immunological Role of Dickkopf-1 (DKK1)

WNT signaling pathway inhibitor Dickkopf-1 (DKK1) is related to cancer progression; however, its diagnostic and prognostic potential have not been investigated in a pan-cancer perspective. In this study, multiple bioinformatic analyses were conducted to evaluate therapeutic value of DKK1 in human cancers. The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) project served as data resources. The Wilcoxon rank test was performed to evaluate the expression difference of DKK1 between cancer tissues and normal tissues. A Kaplan-Meier curve and Cox regression were used for prognosis evaluation. Single-sample gene set enrichment analysis (ssGSEA) was used to evaluate the association of DKK1 expression with the immune cell infiltration. The potential function of DKK1 was explored by STRING and clusterProfiler. We found that the expression level of DKK1 is significantly different in different cancer types. Importantly, we demonstrated that DKK1 is an independent risk factor in ESCA, LUAD, MESO, and STAD. Further analysis revealed that DKK1 had a large effect on the immune cell infiltration and markers of certain immune cells, such as Th1 and Th2 cells. PPI network analysis and further pathway enrichment analysis indicated that DKK1 was mainly involved in the WNT signaling pathway. Our findings suggested that DKK1 might serve as a marker of prognosis for certain cancers by affecting the WNT signaling pathway and tumor immune microenvironment.

DKK1 could affect the function of immune cells, such as T lymphocytes and bone marrow-derived suppressor cells (Katoh and Katoh, 2017). By activating CD4 + and CD8 + T lymphocytes, DKK1 could eliminate myeloma cells in mouse models (Qian et al., 2012). DKK1 also inhibited the secretion of IFN-γ in Th1 cells and induced the production of interleukin (IL)-4, IL-5, IL-10, and IL-13 in Th2 cells (Bais et al., 2005). The inflammation caused by tumor-specific Th1 cells could prevent cancer, but Th2 cells have the opposite function ( Kennedy and Celis, 2008;Lefrancois et al., 2020). By inhibiting β-catenin to prevent clearance by natural killer (NK) cells, DKK1 helps to sustain the stem cell-like  properties of cancer cells (Malladi et al., 2016). Based on these findings, it is necessary to evaluate the role of DKK1 in the cancer immune microenvironment. For several years, numerous studies have been conducted to explore the role of DKK1 in various cancers and revealed the different roles of DKK1 in different cancer types. In this study, we evaluated the pan-cancer expression of DKK1 using The Cancer Genome Atlas (TCGA) dataset. Subsequently, we investigated the association of DKK1 expression with the survival time of patients with different cancers. Finally, we analyzed the effect of DKK1 expression on immune cell infiltration and immune cell markers. The overall process of this research is shown in Figure 1. Our findings deepened our understanding of the roles of DKK1 in cancer progression and prognosis.

Pan-Cancer DKK1 Expression Profile Analysis
The Genotype-Tissue Expression (GTEx) project and RNA-seq datasets from The Cancer Genome Atlas (TCGA) were downloaded from UCSC Xena (https://xena.ucsc.edu/) and used for pan-cancer analysis of DKK1. TOIL was used to reprocess the raw RNA-seq data from GTEx and TCGA databases to correct batch effects and allow for data merging across GTEx and TCGA datasets (Vivian et al., 2017). Expression differences of DKK1 were examined using the Wilcoxon rank test with the threshold of |log 2 FC| >1 and p-value < 0.05.

Survival Analysis
Patients with different types of cancer were segregated into high and low expression groups by the median of the expression level of DKK1. Kaplan-Meier (KM) survival analysis was conducted by R survival and survMiner packages. Cox regression analysis was used to evaluate the relationship between DKK1 expression and overall survival (OS) based on TCGA data. Univariate Cox analysis was performed to select relevant variables, and a multivariate Cox model was used to evaluate the independent prognostic factors. Differences were considered significant when p-values were less than 0.05. All analyses were carried out using R language (version 3.6.3).

Correlations Between DKK1 Expression and Infiltration Immune Cells
Single-sample gene set enrichment analysis (ssGSEA) was used to assess the immune cell infiltration signatures of each individual with LUAD according to the expression level of DKK1 by using the R GSVA package (Hänzelmann et al., 2013). The gene set for immune cell markers was retrieved from the Laboratory of Integrative Cancer Immunology (LICI) (Bindea et al., 2013).

Correlation Between DKK1 Expression and Immune Cell Markers
The correlation module in TIMER (http://timer.cistrome.org/) was used to analyze the correlation between the expression of DKK1 and diverse immune cell markers. The list of different immune factors was obtained from the tumor immune system interaction database (Ru et al., 2019).

PPI Network Construction and Functional Enrichment Analysis
The protein-protein interaction network (PPI network) of DKK1 was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, http://string-db.org) (von Mering et al., 2003) with a minimum required interaction score > 0.7. To identify the hub genes of DKK1 PPI, the maximal clique centrality (MCC) algorithm was analyzed by the Cytohubba (Chin et al., 2014) plugin based on Cytoscape (Shannon et al., 2003). To evaluate the biological functions that DKK1 is involved in, Gene Ontology (GO, http://geneontology.org/) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG, http:// www.kegg.jp/) pathway analyses were performed using the R package clusterProfiler (Yu et al., 2012).

DKK1 Expression Difference in Pan-Cancer
In this study, expression difference analyses of DKK1 were performed between cancer tissues and adjacent normal tissues. DKK1 mRNA expression in cancer tissues from the TCGA database was inconsistent with that in GTEx and TCGA normal tissues (Figures 2A,B). DKK1 expression was significantly higher in cancer tissues with cholangiocarcinoma (CHOL), ESCA, HNSC, LIHC, lung squamous cell carcinoma (LUSC), and STAD than that in their respective adjacent normal tissues. However, DKK1 expression was significantly decreased in bladder urothelial carcinoma (BLCA), kidney chromophobe (KICH), kidney renal papillary cell carcinoma (KIRP), and prostate adenocarcinoma (PRAD).

DKK1 Expression as an Independent Prognostic Factor
Differential DKK1 expression was associated with poor overall survival (OS) in several types of cancer (p < 0.05; Figure 3).  Table 4). Taken together, these results demonstrated that differential DKK1 expression had a non-

Association Between DKK1 Expression and Immune Responses in Cancer
To further validate the role of DKK1 as a potential immune influencer, the relationship between DKK1 expression and immune cell infiltration was estimated. It turned out that DKK1 was strongly correlated with the immune cell infiltration in many types of cancer (|r| > 0.4, p < 0.05). The infiltration level of Th2 cells was positively correlated with the expression of DKK1 in ACC, KICH, MESO, and PAAD ( Figures  5A-D). DKK1 level was also correlated with immune cell infiltration of macrophages (GBM), neutrophils (PRAD), Th1 Figures S1A-F).

Recognizable immune cell markers included B cells, T cells (general), CD8 + T cells, T cells with different functions, M1
and M2 macrophages, TAMs, monocytes, NK cells, neutrophils, and dendritic cells. More directly, our findings provided evidence that the expression of DKK1 was correlated with the level of Th1 markers (STAT1 and IFNG) and Th2 markers (GATA3 and STAT6) in various cancers, such as ACC, KICH, MESO, and PAAD ( Figure 6; Supplementary Table S1). In addition, DKK1 was also correlated with the level of other immune cell markers, such as GBM (macrophage markers), PRAD (neutrophils markers), and TGCT (NK cell markers) ( Supplementary  Table S2). These results suggested that DKK1 might have a large impact on the tumor immune microenvironment.

Protein-Protein Interaction Network, Gene Ontology, and KEGG Pathway Analyses
A PPI network of DKK1 included 46 nodes and 506 edges. The hub genes were screened out using the Cytoscape app cytoHubba plugin. The top five hub genes were DKK1, WNT1, WNT2, WNT3A, and WNT5A ( Figure 7A), which indicated that DKK1 might be involved in the functional regulation of WNT family genes. To further clarify the biological functions of DKK1, KEGG pathway and GO terms analyses were performed. KEGG pathway analysis showed that genes in the DKK1 PPI were mainly enriched in the WNT signaling pathway, basal cell carcinoma, breast cancer, gastric cancer, and hepatocellular carcinoma pathways ( Figure 7B). GO terms were mainly concentrated in regulation of the WNT signaling pathway (biological process, BP), WNT signalosome (cellular component, CC), and frizzled binding (molecular function, MF) ( Figure 7C). These results implied that genes in the DKK1 PPI might work together to participate in cancer progression by the WNT signaling pathway.

DISCUSSION
DKK1 is a well-established crucial participant in the WNT signal pathway and also acts as a major target in drug design. DKN-01, an anti-DKK1 mAb, could block the immunosuppressive effects of DKK1 in the tumor microenvironment (TME) (Haas et al., 2021) and also perform potential antiangiogenic and immunomodulatory activity in combination therapy with gemcitabine/cisplatin in advanced biliary tract cancer (Goyal et al., 2020). Other clinical trials using anti-DKK1 mAb, such as BHQ880 (Fulciniti et al., 2009) and PF-04840082 (Betts et al., 2010) were also carried out to treat certain cancers. The differential expression of DKK1 has been reported in many different cancers; however, there is a lack of comprehensive pan-cancer analysis of DKK1. In the current study, after analyzing the expression level of DKK1 in cancer and normal tissues of 33 cancer types, we found that DKK1 was upregulated in CHOL, ESCA, HNSC, LIHC, LUSC, and STAD, but downregulated in BLCA, KICH, KIRP, and PRAD. Similarly, several studies demonstrated that DKK1 was differentially expressed in a variety of cancers and affected cancer progression by changing cancer proliferation and invasion capabilities (Shi et al., 2013;Zhuang et al., 2017;Fezza et al., 2019). Our prognostic analysis showed that in most cancers (ACC, HNSC, LAML, LUAD, MESO, PAAD, and STAD), the upregulated expression of DKK1 was associated with poor prognosis. This further supported the results from previous studies which also confirmed the relationship between the overexpression of DKK1 and the lower overall survival in HNSC (Gao et al., 2018), NSCLC (Yamabuki et al., 2007), andPAAD (Han et al., 2015). However, in ESCA and KIRC, our data showed that a lower DKK1 level contributed to poor prognosis. Other pancancer analysis also provided evidence that the same gene could result in a controversial consequence in different kinds of cancer. For example, a pan-cancer study showed that there was a significant different expression of NLRP3 in 15 different cancers, and it was used as an independent posterior factor of SKCM (Ju et al., 2021). Similarly, when comparing the expression level of Fam20C in cancer tissues with that in neighboring normal tissues, we found a large variation across different kinds of cancers which indicated the different roles in different cancers . The effect of DKK1 on various cancers may be the result of abnormal activation of WNT signaling (Niida et al., 2004). Our PPI analysis and the pathway enrichment analysis showed that the genes interacting with DKK1 were mainly involved in the WNT signaling pathway, basal cell carcinoma, breast cancer, gastric cancer, and liver cancer. DKK1 has been reported to inhibit the interaction of LRP 5/6 with a WNT signal and the formation of the Fzd-WNT-LRP5/6 complex (Wirths et al., 2003;Ahn et al., 2011). By interfering with DKK1, the activated WNT3a/β-catenin signal had a great influence on cell proliferation, cell cycle acceleration, invasion, and migration (Ren et al., 2021). These results suggested that DKK1 may participate in cancer progression through the WNT signaling pathway.
To unravel the potential mechanism of the predictive value of DKK1 alterations for the tumor immune microenvironment, types of infiltrating immune cells were surveyed. Chronic inflammation is a well-acknowledged risk factor of cancers, we hypothesized that DKK1 influenced cancer prognosis through immune cell infiltration. After conducting immune infiltration analysis, we found that the expression DKK1 was correlated with certain immune cell markers on Th1 and Th2 cells. Wang's study showed that Th2 cells were identified as prognostic immune cells in gastric cancer (Wang et al., 2020). Th2 cells also demonstrated therapeutic potential for adoptive cell therapy (ACT) (Lorvik et al., 2016). Recent studies demonstrated that Th2 responses played a critical role in the pathogenesis of cancers, such as luminal breast cancer (Zhang et al., 2015), prostate and advanced melanoma cancer (Dulos et al., 2012), and myeloma (Tian et al., 2019). The deep understanding of the relationship between DKK1 and Th2 responses will help us comprehend the mechanism of local antitumor response. In recent years, with the development of immune checkpoint inhibitors, infiltrating immune cell markers can not only be used as prognostic markers, but also have received extensive attention as a new type of treatment (Ladanyi, 2015).
Taken together, our study annotated DKK1 expression in a pan-cancer manner and identified that DKK1 could be used as an independent prognosis factor. DKK1 was significantly expressed in various cancers, and it might also be a biomarker for tumor immunity or even targeted therapy. This study also provided evidence of the effect of DKK1 on immune cell infiltration. However, this study has its limitations. Since all analyses were based on online datasets, experimental confirmation from a laboratory is still needed.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://xena.ucsc.edu/