All of the data analyzed in this article were public available. The RNA sequencing data of health tissues was downloaded from Genotype-Tissue Expression (GTEx) portal (https://www.gtexportal.org/) which collects organs from autopsies. The RNA sequencing data, simple nucleotide variations, copy number variations, as well as clinical data of TCGA projects were downloaded from GDC portal (https://portal.gdc.cancer.gov/), which contains datasets of adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal papillary cell carcinoma (KIRP), kidney renal clear cell carcinoma (KIRC), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), sarcoma (SARC), stomach adenocarcinoma (STAD), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma (UCS), uveal melanoma (UVM).

The immunohistochemistry (IHC) results of AKR1B10 in cancers and normal tissues were downloaded from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). STAD enrolled samples from stomach (T-6300), stomach upper (T-62350) and stomach lower (T-63700) Adenocarcinoma; COAD enrolled samples from colon (T-67000) adenocarcinoma; LIHC enrolled samples from liver (T-56000) hepatocellular carcinoma; PAAD enrolled pancreas (T-59000) adenocarcinoma; UCEC enrolled samples from endometrium (T-84000) adenocarcinoma; CESC enrolled samples from cervix (T-83000) squamous cell carcinoma.

The validation datasets (GSE107850, GSE7696, GSE41272, GSE8894 and GSE57303) were available in gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

The association of AKR1B10 expression with stages or grades in each cancer from TCGA database was analyzed by TISIDB (http://cis.hku.hk/TISIDB/), which is an integrated repository portal for tumor-immune system interactions.

The prognosis information including overall survival (OS), progression-free interval (PFI) and disease-free interval (DFI) from clinical follow-up data were integrated with AKR1B10 expression from RNA sequencing data. The cutoff value of AKR1B10 expression was determined by the optimal calculation result. The prognostic value of AKR1B10 was determined by Kaplan-Meier analysis.

The frequency and site of mutations (missense mutations, nonsense mutations, frame shift inserts, frame shift deletions and splice sites) in each cancer based on the simple nucleotide variation dataset processed by MuTect2 software (38) were integrated with AKR1B10 protein domain information acquired from maftools (version 2.2.10) in R software. The tumor mutation burden (TMB)scores were calculated by tmb function of maftool.

The copy number variation (CNV) datasets processed by GISTIC software (39) were integrated with log2(x+0.001) transformed AKR1B10 expression from RNA-seq datasets. A total of 26 types of cancers were analyzed with R software (version 3.6.4) when cancers with sample number less than 3 were excluded.

The correlation between methylation of AKR1B10 promoter and FPKM in each cancer was analyzed by DNMIVD (40) which is an interactive visualization database for DNA methylation profiles of diverse human cancers including TCGA databases.

The expression of AKR1B10 among different organs and cell types from single cell RNA sequencing were interrogated and visualized by Human Cell Landscape (http://bjs.zju.edu.cn/HCL/) (41), which contains 60 human tissue samples from 8 systems.

The stromal, immune and ESTIMATE scores of each cancer sample were calculated by ESTIMATE package (version 1.0.13) in R software based on the gene expression matrix after data normalization and transformation mentioned above. While the purity and homologous recombination deficiency (HRD) data of each cancer were downloaded from previous published literature (42). The expression of AKR1B10 was integrated with purity, stromal scores, immune scores, ESTIMATE scores, and HRD data, and then the correlation was analyzed with corr.test function of psych package (version 2.1.6) in R software.

Normalized pan-cancer data were annotated with gene symbols, and the expression was transformed with log2(x+0.001). The infiltration scores of B cells, T cells CD4, T cells CD8, neutrophils, macrophages and DCs in cancers were calculated with TIMER (43) algorithm from IOBR package (version 0.99.9) (44) in R software. Moreover, the infiltration scores of B_cell_naive, B_cell_memory, plasma_cells, T_cells_CD8, T_cells_CD4_naive, T_cells_CD4_memory_resting, T_cells_CD4_memory_activated, T_cells_follicular_helper, T_cells_regulatory_(Tregs), T_cells_gamma_delta, NK_cells_activated, Monocytes, Macrophages_M0, Macrophages_M1, Macrophages_M2, Dendritic_cells_resting, Dendritic_cells_activated, Mast_cells_resting, Mast_cells_activated, Eosinophils and Neutrophils in cancers were determined by deconvo_CIBERSORT (45) algorithm from IOBR package in R software. The association between AKR1B10 expression and immune cell infiltration scores in each cancer was analyzed with corr.test function of psych package (version 2.1.6).

The expression of 41 chemokine genes, 18 chemokine receptors, 21 MHC genes, 24 immunoinhibitors, 46 immunostimulators and 24 immune checkpoint inhibitors (42) was extracted from normalized pan-cancer RNA-seq data, the correlation of which with AKR1B10 in each cancer was analyzed with corr.test function of psych package (version 2.1.6) after log2(x+0.001) transformed.

The function of AKR1B10 in pan-cancers was analyzed with gene set enrichment analysis (GSEA) software (version 4.2.3). RNA-seq data of 33 cancers from TCGA project were loaded into GSEA software. C2.cp.kegg.v7.5.1.symbols was chosen as a gene set database, and the number of permutations was set to 1000. The no_collapse was used to gene symbols, and phenotype was set for permutation type. Enriched result lists were first filtered with threshold of |NES|≥1.5, NOM p-value<0.05, and FDR q-value<0.25. Filtered results were further presented with a heatmap of NES.

Gastric cancer tissue microarray (A962) was purchased from WeiAoBio (Shanghai, China). Multiplex immunofluorescence staining procedure on paraffin-embedded tissue section was followed by the instruction of supplier (Absin, Cat No. abs50013). In briefly, the tissue section was deparaffinize, rehydration, antigen retrieval by microwave, serum blocking, primary antibody incubation, secondary antibody incubation, fluorescence staining signal amplification (multiple rounds), staining nuclear with DAPI and mounting. The primary antibodies were AKR1B10 (Abnova, Cat No. H00057016-M01, 1:200), CD68 (Absin, Cat No. 171440, 1:200), CD19 (Abcam, Cat No. 134114, 1:100), PanCK (Absin, Cat No. 123684, 1:400). The staining with the horseradish peroxidase conjugated secondary antibody and the TSA tyramide signal amplification (520-TSA, 620-TSA, 570-TAS, 690-TSA) were supplied within the staining kit. Multispectral images were analyzed. The percentage of positive staining cell was calculated as number of positive cells divide total cell number (P). For the intensity (I), extensive staining was recorded as 3, moderate staining was recorded as 2, weak staining was recorded as 1. IHC   score = ∑ ( P × I ) = ( w e a k   p e r c e n t a g e × 1 ) + ( m o d e r a t e   p e r c e n t a g e × 2 ) + ( s t r o n g   p e r c e t a g e × 3 ) .

The difference of AKR1B10 expression between cancers and corresponding normal tissues was analyzed by unpaired t-test. The significance between two groups was analyzed by unpaired Wilcoxon Rank Sum and Signed Rank tests. The significant difference among multiple groups was analyzed by the Kruskal test. All of the correlations were presented as spearman’s correlation coefficient.

To demonstrate the expression of AKR1B10 in human normal tissues, we draw an anatomy map based on AKR1B10 expression from GTEx databases. The result showed that the AKR1B10 highly expressed in stomach and bladder, followed by esophagus, colon, liver, skin and adrenal gland, and AKR1B10 is barely expressed in brain and lymph node ( Figure 1A ). Overall, the gastrointestinal system showed a high-medium AKR1B10 expression, while development gastrointestinal system, genitourinary system, and gynecological system showed a medium-low AKR1B10 expression. The single-cell RNA sequencing (ScRNA Seq) data from Human Cell Landscape (HCL) also showed that the AKR1B10 is highly expressed in the jejunum, stomach, esophagus and transverse-colon ( Figures 1B, C ).

AKR1B10 expression is significantly altered in various types of cancers. In general, AKR1B10 mRNA down-regulated in most of gastrointestinal cancer, genitourinary cancer and PCPG, but up-regulated in most of developmental gastrointestinal cancer, some of gynecological cancer, glioma and lung cancer when compared with normal tissues from TCGA databases ( Figure 2A ). Consistently, AKR1B10 protein decreased in STAD and COAD, but increased in LIHC, PAAD, UCEC and CESC when compared with their individual normal tissues ( Figure 2B ). Moreover, we analyzed the association of AKR1B10 mRNA expression with clinical features. The results showed AKR1B10 mRNA expression was significantly correlated with grades of HNSC and KIRC ( Figure 2C ). AKR1B10 was also significantly correlated with stages of KIRC, KIRP, PAAD and ACC ( Figure 2D ).

Both genetic and epigenetic alteration can result in gene expression change. Thus, we analyzed the relationship of DNA mutation, copy number variation (CNV) and DNA methylation with AKR1B10 expression. Mutation analysis indicated that scatter missense mutations (41), nonsense mutations (3), frame shift inserts (1), frame shift deletions (3) and splice sites (1) were observed within the coding sequence of AKR1B10 in 17/33 cancers, of which UCEC had the highest incidence rate of 3.4% ( Figure 3A ). CNV analysis showed that there are no distinct change in AKR1B10 expression among different types of CNV in most cancers ( Figure 3B ). While the methylation analysis showed that AKR1B10 expression was significantly negatively associated with DNA methylation in 12 types of cancer. Specifically, LUSC, CHOL and BLCA showed a moderate correlation (0.4< r <0.6), while LUAD, HNSC, ESCA, LIHC, PAAD, BRCA and CESC showed a weak correlation (0.2< r <0.4) ( Figure 4 ). Taken together, epigenetic might contribute more to the alteration of AKR1B10 expression in pan-cancer.

Tumor microenvironment (TME) is the non-cancerous cells surround a tumor including fibroblasts and immune cells, which are reported to play a role in tumor growth and drug resistance Purity and ESTIMATE (including stromal score, immune score and ESTIMATE score) algorithm are novel computational methods to infer tumor purity from bulk cell genomic analysis. Purity patterns differ by cancer type ( Supplementary Figure 1A ). High purity is correlated with shorter OS of PCPG, ACC, ESCA, LUAD, READ,CHOL, CESC and SARC. It is also correlated with longer OS of UVM, MESO, STAD, PAAD, KIRC, BLCA and TGCT ( Supplementary Figure 1B ). High stromal score is correlated with shorter OS in LGG,GBM, LAML, UVM, MESO, LUSC, STAD, COAD, KIRC, KIRP and BLCA. It is also correlated with longer OS in THCA, ACC, SKCM, LUAD, CHOL, LIHC and SARC ( Supplementary Figure 2 ). High immune score is correlated with shorter OS in LGG, LAML, UVM, ESCA and KIRC. It is also correlated with longer OS in PCPG, THCA, ACC, SKCM, MESO, LUAD, HNSC, COAD, CHOL, LIHC, KICH, BRCA and CESC ( Supplementary Figure 3 ). High ESTIMATE score is correlated with shorter OS in LGG, GBM, LAML, UVM, LUSC, ESCA, STAD, PAAD, KIRC, KIRP and TGCT. It is also correlated with longer OS in PCPG, THCA, ACC, SKCM, LUAD, CHOL, LIHC, CESC and SARC ( Supplementary Figure 4 ).

From the HCL, we knew the AKR1B10 expressed in the epithelial cells, goblet cells, basal cells, enterocytes and enterocyte progenitors, gastric chief cells and hepatocytes, as well as some T cells, B cells, fibroblasts and endothelial cells in normal tissues ( Figures 1D–F ). The function of AKR1B10 in cancer have attracted the attention of scientists, the role of AKR1B10 in TME remains elusive. Herein, we analyzed the association of AKR1B10 expression with tumor purity ( Figure 5A ). There were 7 cancers stood out, of which the MESO showed positive correlation (r>0.2), but the GBM, THYM, THCA, READ, SKCM, and KIRP showed negative correlation (r<-0.2). Further analysis with ESTIMATE algorithm showed the LGG, GBM, THCA, ACC and UVM were top 5 cancers positively associated with ESTIMATE scores, immune scores and stromal scores, but the CHOL and LUSC were top 2 cancers negatively associated with ESTIMATE scores, immune scores and stromal scores ( Figure 5B ). The correlations in LGG, GBM, CHOL and non-small cell lung cancer (NSCLC) were validated with GEO datasets ( Supplementary Figure 5 ). Taken LGG, GBM, LUSC and CHOL as examples, we find the AKR1B10 expression is elevated in purity low group, stromal score high group, immune score high group and ESTIMATE score high group in LGG and GBM ( Figure 5C ); AKR1B10 expression is also elevated in the purity high group, stromal score low group, immune score low group and ESTIMATE score low group in LUSC and CHOL ( Figure 5C ). Together these data may imply the potential role of AKR1B10 in TME.

We analyzed AKR1B10 expression and immune cells infiltrations with both TIMER and CIBERSORT algorithms. TIMER analysis showed the AKR1B10 expression was associated with multiple immune cells, especially CD4+ T cell, neutrophils and dendritic cells, and these associations were positive in GBM, THCA, THYM and PRAD, but negative in UVM, MESO, LUSC ( Figure 6A ). CIBERSORT analysis showed that the association was mainly focused on T cells CD4+ memory resting, T cells regulatory (Tregs), and dendritic cell resting ( Figure 6B ). AKR1B10 exhibited a broad correlation with immunomodulators involved in 5 immune pathways, among which the correlation with MHC pathway was highest ( Figure 7A ). Our results also supported an association of AKR1B10 expression with dendritic cells and T cells, indicating that AKR1B10 may be involved in antigen presentation.

The immune therapy based on immune checkpoint inhibitors (ICI) is a novel and promising anti-tumor therapy mode. Our analysis showed the correlation between AKR1B10 expression and various ICIs in LGG, GBM, PCPG, THCA, ACC, UVM, KICH, and TGCT, but no correlations in THYM, MESO, LUSC, ESCA, and CHOL ( Figure 7B ).

Tumor mutation burden (TMB) was defined as the number of gene mutations (e.g., base substitutions, insertions and deletions in coding area) per million bases in somatic cells. It is an indicator of mutation frequencies. The higher TMB was usually associated with the favored ICI response regardless of PD-L1 status in various types of tumors. Our analysis showed AKR1B10 expression was significantly associated with TMB in 12 types of cancers, with HNSC correlation coefficient larger than 0.2, THYM, SKCM, ESCA and READ correlation coefficient smaller than -0.2 ( Figure 7C ). Together, these results indicate that AKR1B10 may be a potential predictor and target of immune therapies.

Homologous recombination deficiency (HRD) refers to HRR dysfunction caused by HRR related gene mutations or epigenetic inactivation. HRD often presents in various tumor and produces specific, quantifiable and stable genomic alterations. Clinical research showed HRD status was related to the sensitivity of PARP inhibitors and platinum chemotherapy. We found that AKR1B10 expression significantly correlated with HRD in 9 types of cancer, with ESCA, KIRP and PCPG correlation coefficient lager than 0.2 ( Figure 7D ).

We further explored AKR1B10 function in pan-cancer cells through GSEA. We noticed that AKR1B10 was enriched with cancer cell features and pathways, as well as apoptosis, cell cycle, autophagy, and ECM interaction processes in some cancers. Further studies revealed that AKR1B10 was most remarkably and positively enriched with multiple metabolism processes, including drugs, hormones, saccharides, and amino acids in pan-cancer ( Figure 8 ). These results indicate a potential role of AKR1B10 in cancer development and progression.

Subsequently, we evaluated the prognostic value of AKR1B10 in multiple types of cancers through Kaplan-Meier analysis. As shown in Figure 9 , the increased AKR1B10 expression was significantly related to poor overall survival (OS) in GBM, THCA, UVM, SKCM, MESO, LIHC, PAAD, KICH, KIRC, KIRP, and PRAD, but better OS in ACC, COAD, and READ. In addition, elevated AKR1B10 expression was associated with shorter PFI in GBM, SKCM, STAD, PAAD, KIRC, KIRP, TGCT, and BRCA, as well as decreased DFI in ACC, STAD, LIHC, BLCA, and TGCT. On the contrary, increased AKR1B10 expression was associated with longer PFI in ACC, COAD, and CESC, and also the prolonged DFI in THCA, BRCA, and CESC.

Our previous analysis showed stomach was the organ with highest AKR1B10 expression. Some studies indicated that decreased AKR1B10 expression promote M2 macrophage marker expression in THP1 macrophage cells (46), promote NK cell kill activity on liver cancer cell Huh7 (47), and knock down of AKR1B8 transgenic mice (homologous of AKR1B10 in mice) show a immunity disorder in colon (48), however, the correlation between AKR1B10 expression and immune infiltration in STAD remains elusive. We evaluated AKR1B10 expression in paired normal and gastric cancer tissues, and the correlation of AKR1B10 expression with immune cells in gastric cancer ( Figure 10 ). The results showed that AKR1B10 expression was dramatically decreased in gastric cancer compared to normal tissues ( Figures 10A, C ). Meantime, the CD19 expression (B cell marker) was also decreased in gastric cancer ( Figures 10A, D ), and CD68 (Macrophage marker) expression elevated in gastric cancer ( Figures 10A, E ). The PanCK (Epithelial marker) expression difference is not significant between gastric cancer and normal tissue ( Figures 10A, F ).

As AKR1B10 expression almost diminished in gastric cancer, we evaluated the relationship between AKR1B10 and immunity by dividing samples with AKR1B10 expression into AKR1B10 high group with H score > 10, AKR1B10 medium group with 10 > H score > 0.5 and AKR1B10 low group with H score < 0.5. We found CD68 show an elevated trend as AKR1B10 expression decreased ( Figures 10G, I ), which is consistent with the TIMER analysis ( Figure 6A ). However, either expression difference of CD19 or PanCK is clear ( Figures 10H, J ).