Bioinformatics analysis of the clinical value and potential mechanisms of AHNAK2 in papillary thyroid carcinoma

Background: AHNAK2 has been recently reported as a biomarker in many cancers. However, a systematic investigation of AHNAK2 in papillary thyroid carcinoma (PTC) has not been conducted. Results: AHNAK2 is overexpressed in PTC tissues and could be an independent prognostic factor. AHNAK2 expression was significantly high in patients with advanced stage, advanced T classification, lymph node metastasis, increased BRAF mutations and decreased RAS mutations. Cell adhesion-, cell junction-, and immune-related pathways were the most frequently noted in gene set enrichment analysis. AHNAK2 expression in PTC was positively correlated with immune infiltration and negatively correlated with AHNAK2 methylation. AHNAK2 expression was significantly positively correlated with tumor progression and poor overall survival (OS) in pan-cancer patients. Conclusions: AHNAK2 is a good biomarker for the diagnosis and prognosis of PTC. AHNAK2 may promote thyroid cancer progression through cell adhesion-, cell junction-, and immune-related pathways. Methylation may act as an upstream regulator to inhibit the expression and biological function of AHNAK2. Additionally, AHNAK2 has broad prognostic value in pan-cancer. Methods: Based on The Cancer Genome Atlas (TCGA) data, we screened AHNAK2-related genes through weighted gene coexpression network analysis and explored the clinical value and the potential mechanism of AHNAK2 in PTC by multiomics analysis.


INTRODUCTION
Thyroid cancer (TC) is the fifth most common cancer among women in the United States, with an estimated 52,890 new TC cases nationwide in 2020 [1], and its incidence is increasing [2,3]. Papillary thyroid carcinoma (PTC) is the most common subtype of TC, comprising approximately 80-85% of all TCs [4,5]. Although most well-differentiated PTC patients have an excellent overall prognosis (the 5-year survival rate of such patients is greater than 97%) [4], PTC recurrence and metastasis still hinder clinical management and prognosis in certain patients. A study with a median follow-up of 27 years also reported that 28% of PTC patients had recurrent disease, and 9% of patients died from PTC [6]. Additionally, overdiagnosis and overtreatment are common problems associated with indolent diseases. The screening and identification of indolent TC and the treatment of these overdiagnosed cancers can increase the risk of injury to patients [7,8]. Therefore, the use of novel and sensitive biomarkers to effectively identify specific PTC patients and provide personalized treatment has become an urgent need.
AHNAK nucleoprotein 2 (AHNAK2) is a giant protein (600 kDa) with a PDZ domain that was first detected in 2004 [9,10]. AHNAK2 was recently reported as a biomarker for the diagnosis and prognosis of pancreatic AGING ductal adenocarcinoma (PDAC) [11][12][13], clear cell renal cell carcinoma (ccRCC) [14], thymic carcinoma [15], bladder cancer [16], gastric cancer [17], and uveal melanoma (UM) [18]. Although the biological function of AHNAK2 in cancer remains unclear, some progress has been made regarding its mechanism in tumors. AHNAK2 is a critical element of the stress-induced FGF1 export pathway [19]. Li et al. [18] reported that AHNAK2 may play a role in promoting the proliferation and migration of UM cells by regulating the PI3K signaling pathway. Wang et al. [14] also found that AHNAK2 is a target gene of HIF1α, which mediates epithelial-mesenchymal transition (EMT) and stem cell characteristics driven by the hypoxia pathway, thereby promoting the progression of ccRCC.
The purpose of this study was to explore the clinical value and potential mechanism of AHNAK2 in PTC. Based on the bioinformatics method, we conducted a systematic multiomics study of AHNAK2 in PTC. Multiomics single-gene studies on the candidate oncogene AHNAK2 have enriched our understanding of the molecular function of AHNAK2 and have provided a new piece of the puzzle to the mechanism of cancer development and a basis for subsequent research on AHNAK2. The significant diagnostic and prognostic value of AHNAK2 may also be used in clinical practice to assist doctors in the diagnosis and treatment of PTC patients.

Diagnostic value of AHNAK2
The mRNA expression level of AHNAK2 was significantly higher in PTC tissues than in normal tissues (P<0.0001) ( Figure 1A). By performing receiver operating characteristic (ROC) curve analysis to discriminate PTC from normal thyroid tissues, we found that the area under the curve (AUC) of the AHNAK2 expression level was 0.88, suggesting that it could be a good diagnostic biomarker ( Figure 1B).
Additionally, immunohistochemistry performed with four different AHNAK2 antibodies showed deeper staining levels in PTC tissues than in normal thyroid tissues, suggesting higher protein expression in PTC ( Figure 1C).

Prognostic value of AHNAK2
Kaplan-Meier (K-M) survival analysis showed that the high AHNAK2 (H-AHNAK2) group was associated with worse progression-free survival (PFS) than the low AHNAK2 (L-AHNAK2) group (P=0.0005) ( Figure 1D). Multivariate Cox regression analysis suggested that AHNAK2 could become an independent predictor after adjusting for other parameters, including age, gender, American Joint Committee on Cancer (AJCC) stage, BRAF mutations, RAS mutations and tumor mutation burden (TMB) ( Table 1) (PFS as the ending indicator). In summary, AHNAK2 could be an important tool for distinguishing PTC patient clinical outcomes.

Clinical value of AHNAK2
We analyzed the relationships between AHNAK2 and clinical parameters, including age, gender, stage, metastasis, N classification, T classification, pathologic type, and mutations in BRAF and RAS ( Table 2, Supplementary Figure 1). AHNAK2 expression was significantly high in patients with an advanced stage, an advanced T classification, lymph node metastasis, increased BRAF mutations, decreased RAS mutations, and lower follicular PTC and higher tall-cell PTC proportions.
Interestingly, in the study of the relationship between AHNAK2 expression and gene mutations in PTC (Supplementary Figure 1A), we found that AHNAK2 expression was positively correlated with BRAF mutations (Supplementary Figure 1B

Screening modules and genes related to AHNAK2
To construct a weighted coexpression network and screen for modules and genes related to AHNAK2, 2696 differentially expressed genes (DEGs) between PTC and normal tissues from The Cancer Genome Atlas (TCGA) were selected and subjected to weighted gene coexpression network analysis (WGCNA) (475 patients with complete clinical information were selected) (Figure 2A and 2B). After a series of adjustments for WGCNA parameters, the DEGs were divided into 11 modules by average linkage hierarchical clustering ( Figure 2C-2F). The blue module, which contains 407 genes, exhibited the highest correlation with AHNAK2 expression (Pearson's correlation coefficient =0.65, P<0.0001) ( Figure 2G). Fifty-seven genes in the blue module were selected as hub genes (absolute module membership [MM] > 0.5 and absolute gene significance [GS] > 0.5) ( Figure 2H).

Analysis of the potential mechanism of AHNAK2
Cell adhesion-and cell junction-related pathways were the most frequently noted pathways in the functional enrichment analysis of the 407 blue module genes. "Cell-cell adhesion via plasma membrane adhesion molecules", "cell-cell junction", and "cell adhesion   AGING mediator activity" were the most significant Gene Ontology (GO) terms for cellular components, biological processes and molecular functions, respectively ( Figure 3A). "Cytokine−cytokine receptor interaction" was the most significant pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis ( Figure 3B). In addition, we constructed a proteinprotein interaction (PPI) network with the 57 hub genes and found that ICAM1 and FN1 play key roles ( Figure  3C).
Additionally, gene set enrichment analysis (GSEA) was used to identify the mechanism and functional differences between the L-AHNAK2 group and H-AHNAK2 group. The H-AHNAK2 group was enriched in the following terms: 1. Cell adhesion-and cell junction-related pathways: cell adhesion molecules, focal adhesion, adherens junction, and tight junction; 2. Immune-related pathways: cytokine-cytokine receptor interaction, JAK-STAT signaling pathway, natural killer cell-mediated cytotoxicity, chemokine signaling pathway, and leukocyte transendothelial migration; and 3. Cancer-related pathways: p53 signaling pathway, small cell lung cancer, and pathways in cancer ( Figure  3D). Conversely, the L-AHNAK2 group was enriched in metabolism-related pathways: glycine serine and threonine metabolism, fatty acid metabolism, beta alanine metabolism, peroxisome, valine leucine and isoleucine degradation, oxidative phosphorylation, and tryptophan metabolism ( Figure 3E).

Further analysis of the relationship between AHNAK2 and immune infiltration
The results of the ESTIMATE analysis suggested that the H-AHNAK2 group had a higher immune score and stromal score than the L-AHNAK2 group, while the tumor purity score was lower ( Figure 4A-4D). These results indicate that the proportions of immune cells and stromal cells are higher, and that the proportion of tumor cells is lower in the H-AHNAK2 group than in the L-AHNAK2 group. The TIMER analysis showed that AHNAK2 has significant positive correlations with B cells, CD4 + T cells, macrophages, neutrophils and dendritic cells (DCs) in PTC ( Figure 4E). The results of the TISIDB analysis suggested that AHNAK2 has positive correlations with 28 tumor-infiltrating lymphocyte (TIL) types ( Figure 4F) and human leukocyte antigens (HLAs) ( Figure 4G) across human cancers, which is particularly significant in THCA (papillary thyroid carcinoma in TCGA).

Analysis of AHNAK2 methylation in PTC
Our data show that the methylation level of AHNAK2 in PTC tissues is significantly lower than that in normal tissues ( Figure 5A). A significant negative correlation was identified between the methylation level of AHNAK2 and the expression level of AHNAK2 (Spearman r=-0.678, P<0.0001) ( Figure 5B). Among the DNA methylation sites of AHNAK2, cg06799735, cg01513078, cg11138227 and cg23385208 had the most significant negative correlations with the AHNAK2 expression level (Table 3).
Moreover, we found significant connections between AHNAK2 methylation and clinical parameters such as tumor purity, stage, T classification, N classification and histological type ( Figure 5C-5G).
To understand the biological significance of AHNAK2 methylation in THCA, a functional module of LinkedOmics was used to examine AHNAK2 coexpression pattern in the THCA cohort. Based on RNAseq, we screened 9,944 genes related to AHNAK2 methylation (false discovery rate (FDR) <0.01) ( Figure  5H). The top 50 significant genes that were positively and negatively correlated with AHNAK2 methylation are presented as two heatmaps in Figure 5I. The GO (biological process) analysis results derived by GSEA were significant; the results indicated that AHNAK2 methylation coexpressed genes participate primarily in metabolism-related pathways, while genes related to activities such as the adaptive immune response, granulocyte activation, leukocyte migration, regulation of leukocyte activation, and positive regulation of the defense response were inhibited. KEGG pathway analysis also showed that genes related to the complement and coagulation cascades, cytokinecytokine receptor interaction, cell adhesion molecules, and Epstein-Barr virus infection, among other pathways, were inhibited. These results illustrate the extensive effects of AHNAK2 methylation on transcriptome and immune-related pathways.

Generalization value of AHNAK2 in pan-cancer
To investigate whether AHNAK2 has broad value, we performed a series of studies on AHNAK2 across all cancers. Gene Expression Profiling Interactive Analysis (GEPIA) showed that AHNAK2 expression status varies in different cancers ( Figure 7A). TISIDB showed that the high expression of AHNAK2 in pan-cancer tended to be accompanied by an advanced tumor stage ( Figure  7B) and short overall survival (OS) ( Figure 7C). K-M survival analysis showed that the high AHNAK2 group had significant associations with short OS in bladder urothelial carcinoma (BLCA), glioblastoma multiforme (GBM), lung adenocarcinoma (LUAD), mesothelioma (MESO), pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM) and uveal melanoma (UVM) (Figure 7D-7J). AHNAK2 is an extremely significant gene in the survival statistics of a large number of pan-cancer samples (N = 9497, HR = 1.4, P<0.0001) ( Figure 7K).

DISCUSSION
Our research suggests that AHNAK2 can be a biomarker for the diagnosis and prognosis of PTC. AHNAK2 functions mainly in cell adhesion-, cell junction-, and immune-related pathways. AHNAK2 also has extensive prognostic value in pan-cancer.
We previously found that AHNAK2 has a significant negative correlation with the prognosis of PTC and used AHNAK2 as one of the model genes for incorporation into prognostic models. Using this preliminary finding, we conducted a systematic study on the role of AHNAK2 in PTC. In recent years, many studies have reported that AHNAK2 could be a biomarker for the diagnosis and prognosis of PDAC, ccRCC, thymic carcinogenesis, bladder cancer, gastric cancer, and UM [11][12][13][14][15][16][17][18]. Kim et al. [20] also noted that high CXCL16 expression combined with the high expression levels of related genes (e.g., AHNAK2 and THBS2) was associated with poor recurrence-free survival in PTC patients. However, a systematic investigation of AHNAK2 in PTC has not been conducted. Our research fills the knowledge gap in the study of AHNAK2 in PTC and reveals the potential mechanism of AHNAK2 in PTC. This study also provides strong support for AHNAK2 as a prognostic biomarker for pan-cancer.
AHNAK2 has been described as a diagnostic biomarker in PDAC [11,12], thymic carcinoma [15], and bladder cancer [16]. We found that overexpression of AHNAK2 in PTC tissues can be used to distinguish between PTC and normal tissues. We further verified the outstanding clinical diagnostic efficacy of AHNAK2 in 8 GEO verification cohorts.
AHNAK2 has been described as a prognostic biomarker in ccRCC [14] and PDAC [11,13]. Our results show AGING that AHNAK2 can be used as an independent prognostic predictor in PTC and to classify the PTC population and provide personalized treatment to reduce harm to patients caused by overtreatment. AHNAK2 may play a role in promoting the proliferation and migration of UM cells by regulating the PI3K signaling pathway [18]. A recent study also found that AHNAK2 is a target gene of HIF1α, which mediates EMT and stem cell characteristics driven by the hypoxia pathway, thereby promoting the progression of ccRCC [14]. AHNAK2 expression is significantly high in patients with an advanced cancer stage, an advanced T classification, AGING lymph node metastasis, and increased BRAF mutations, showing a positive correlation with PTC progression. These findings imply that AHNAK2 may exert a tumorpromoting effect in PTC.
Interestingly, the two genes with the highest mutation rates in PTC [21] exhibited a notable trend: BRAF mutations were positively associated with AHNAK2 expression, while RAS mutations were negatively correlated with AHNAK2 expression.
Through functional enrichment analysis of the WGCNA-selected module eigengenes, AHNAK2 was determined to be most closely related to cell adhesion and cell junction pathways and inflammatory pathways in PTC. The GSEA results also revealed that AHNAK2 may promote progression in PTC through cell adhesion-, cell junction-, immune-, and cancer-related pathways.
Previous studies have noted that AHNAK2, as a scaffolding protein, is associated with voltage-gated calcium channels [9], and calcium plays a key role in the regulation of cell adhesion and cell junctions [22,23]. Thus, these results can help explain the enrichment of cell adhesion-and cell junction-related pathways in our study.
Related reports of AHNAK2 in systemic lupus erythematosus (SLE) indicate its potential link with immunity [24,25]. AHNAK2 is a critical element of the stress-induced FGF1 export pathway [19]. FGF1 is a nonclassically released growth factor that regulates carcinogenesis, angiogenesis and inflammation in tumors. Based on the results of our previous research and functional analysis, we further explored the relationship between AHNAK2 and immunity in PTC. Through ESTIMATE analysis, we found that high AHNAK2 expression indicates high degrees of immune and stromal cell infiltration and low tumor purity in the PTC microenvironment. The TIMER analysis showed that AHNAK2 is positively correlated with immune infiltration. The correlation analysis between AHNAK2 and 28 immune cell types in TISIDB showed that the entire immune system and AHNAK2 are positively correlated in PTC. The above results suggest that AHNAK2 may promote the occurrence and development of PTC through immune-related pathways.
AHNAK2 is uniquely monomethylated by the protein lysine monomethyltransferase (SMYD2) and is involved in human esophageal squamous carcinoma (ESCC) [26]. Ohmura et al. [27] reported that the methylation of AHNAK2 is associated with chemotherapy resistance in Epstein-Barr virus (EBV)-associated gastric cancer. Our research on PTC epigenetics has also shown that methylation may act as an upstream regulatory mechanism of AHNAK2 expression to inhibit its biological function.
Our study found that although the expression levels of AHNAK2 differ in various tumors, they are significantly positively correlated with tumor progression and short OS in pan-cancer patients, indicating that AHNAK2 also has universal prognostic value in pan-cancer.
Our study found that AHNAK2 can be used as a diagnostic and prognostic biomarker in PTC. AHNAK2 is believed to be able to identify critical cases of PTC and avoid overdiagnosis. The prognostic value of AHNAK2 should assist doctors in providing personalized treatment to patients and avoiding harm to patients caused by overtreatment. Studies in pancarcinomas have shown the broad clinical application value of AHNAK2. Multiomics single-gene studies on the candidate oncogene AHNAK2 have enriched our understanding of the molecular function of AHNAK2 and have provided a new piece of the puzzle to the mechanism of cancer development and a basis for subsequent research on AHNAK2.
However, this study has some limitations. First, because TCGA is the only public database with sufficient data on both PTC gene expression and prognosis, we did not use other independent cohorts to verify the prognostic value of AHNAK2 in PTC. Second, the molecular functions of AHNAK2 in PTC were analyzed through bioinformatics, and the conclusions have not been further confirmed experimentally. However, related experimental research is ongoing.
AHNAK2 can be used as a diagnostic and prognostic biomarker in PTC. AHNAK2 may promote the progression of PTC through cell adhesion-, cell junction-, and immune-related pathways, and methylation may serve as an upstream regulator to inhibit the expression and biological functions of AHNAK2. AHNAK2 also has extensive prognostic value in pancancer. Research on the cancer candidate gene AHNAK2 will help us understand the development and prognosis of cancer and provide a piece of the puzzle to overcome cancer.

Comparison of AHNAK2 expression between THCA and normal thyroid tissues
The mRNA expression levels of AHNAK2 in PTC and normal thyroid tissues were analyzed using TCGA data.
We defined the top quarter as the H-AHNAK2 group (based on mRNA expression rank in the THCA dataset) and the remainder as the L-AHNAK2 group. Additionally, a direct comparison of AHNAK2 protein expression between PTC and normal thyroid tissues was performed by analyzing immunohistochemistry images from The Human Protein Atlas (https://www. proteinatlas.org) [28].

WGCNA reveals key modules and hub genes related to AHNAK2
First, DEGs between PTC and normal tissues in the TCGA dataset were identified with the "limma" package in R (adjusted p-value <0.05 and |log 2 FC| >1).
WGCNA aims to identify coexpressed gene modules and to explore associations between gene networks and phenotypes of interest, as well as core genes in the network [29]. We divided the DEGs into different gene modules using the "WGCNA" package in R (minModuleSize = 30).
We defined the module with the highest absolute module significance as the key module. We defined the hub genes as those that satisfied the following criteria: 1. In the key module; 2. An absolute value of MM >0.5, which represents Pearson's correlation between the gene and the module; and 3. An absolute value of GS >0.5, which represents Pearson's correlation between the gene and the clinical parameter.

GO and KEGG enrichment analyses
The genes in the key module were subjected to GO and KEGG enrichment analyses with DAVID 6.7 (https://david-d.ncifcrf.gov/).

PPI network analysis
The hub genes were input into the STRING database (version 11.0, https://string-db.org/) to analyze the interaction network of the hub gene-encoded proteins. The minimum required interaction score was set to 0.4 (medium confidence), and the results were visualized using Cytoscape (v3.7.2).

ESTIMATE
ESTIMATE [30], a method that uses gene expression signatures to infer the fractions of stromal and immune cells in tumor samples, was used to evaluate the levels of immune cell infiltration (immune score), the stromal content (stromal score), the stromal-immune comprehensive score (ESTIMATE score) and tumor purity for each THCA sample.

TIMER
The TIMER online database [31], a web server for comprehensive analysis of tumor-infiltrating immune cells, was used to analyze and visualize associations between AHNAK2 expression and the abundance of 6 tumor-infiltrating immune cell subtypes (B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and DCs).

TISIDB
We used TISIDB, a comprehensive repository portal for tumor-immune system interactions, to determine the Spearman correlations between AHNAK2 expression and 28 TIL types across human cancers [32].

Analysis of AHNAK2 methylation in PTC
Based on DNA methylation data downloaded from TCGA, we compared the methylation levels of AGING AHNAK2 between PTC and normal tissues. We analyzed the correlation between AHNAK2 methylation and AHNAK2 expression in PTC. We also performed a correlation analysis between each methylation site of AHNAK2 and AHNAK2 expression.
LinkedOmics is an online platform for analyzing multiomics data across 32 cancer types (http://www.linkedomics.org/) [33]. We downloaded and visualized the correlation results between AHNAK2 methylation and clinical parameters from LinkedOmics. Pearson correlations were used to analyze the AHNAK2 methylation-induced coexpression of genes, and the results are displayed with a volcano map and a heatmap. The functional modules of LinkedOmics are based on AHNAK2 methylation data, and GSEA was performed with data from GO (biological process) and KEGG pathway databases. We set the minimum number of genes (Size) =10 and simulations=1000 and selected a weighted set cover to reduce redundancy.

Generalization of conclusions in pan-cancer
GEPIA (http://gepia.cancer-pku.cn/) is used as a tool to generalize conclusions in pan-cancer [34]. A comparison of AHNAK2 gene transcripts between PTC and normal thyroid tissues in pan-cancer was performed using TCGA and Genotype-Tissue Expression (GTEx) data (p-value <0.05 and |log 2 FC| >1). K-M curves in pan-cancer were generated using the cohort of TCGA data. We used TISIDB to analyze relationships between AHNAK2 expression and tumor stage in pan-cancer, as well as associations between AHNAK2 expression and OS.

Statistical analysis
K-M survival analysis was performed with the "survival" package (PFS was used as the endpoint). The chi-square test was used to test differences in clinical parameters between the L-AHNAK2 and H-AHNAK2 groups. The Spearman method was used to test for correlations. The Mann-Whitney test was used to compare two groups. The log-rank method was used to calculate significant P values for survival. R (v3.6.0) and SPSS version 25.0 software programs were used for statistical processing. Data were visualized with GraphPad Prism V.8.0 software and R.

AUTHOR CONTRIBUTIONS
ZX designed the analytical strategies, performed data analyses and wrote the manuscript. YL, XL, YZH, SW, SYW, JS and YCH performed data analyses, JZ conceived the research and wrote the manuscript.