Transcriptome sequencing data of STAD patients and matching clinical features were sourced from the TCGA database. Validation cohorts were acquired from GSE84437.

We used the R package “ConsensusClusterPlus” to further investigate the associations between histone acetylation-related genes and STAD subtypes and determine number and stability of clusters. Based on the cumulative distribution functions, the optimal cluster number of k-value was obtained (20). The KM survival analysis was conducted using the R package “survival.” Differential genes between the two clusters and histone acetylation-related genes are shown in Supplementary Tables 1, 2 .

The univariate Cox regression algorithm was employed to filter out histone acetylation-related genes independently related to prognosis, and LASSO Cox regression was carried out to further narrow the range of candidate genes using the R package “glmnet.” The risk scores were calculated for STAD patients as R i s k   S c o r e = ∑ i n X i * Y i (X: coefficients, Y: gene expression level), and based on the middle value, the STAD patients were equally assigned into low- and high-risk groups. Then, Kaplan-Meier analysis was conducted, and PCA, t-SNE and ROC curves were plotted. The R package “pheatmap” was applied to draw the heatmap and applied the GSE84437 cohort to validate the risk signature. The prognostic significance of the signature was analyzed with univariate and multivariable Cox regression. The Eight-Gene Prognostic Model is shown in Supplementary Table 3 .

A nomogram based on independent prognostic parameters was developed using R package “rms” and “regplot” for STAD patients to predict the survival probability. In addition, calibration curve and AUC were applied to evaluate the nomogram regarding accuracy, and decision curves analysis was also performed to confirm the model with respect to its clinical effectiveness.

The TCGA database was used to retrieve sample mutation data, and the R package “maftools” was used to estimate the mutation status between the two groups. Immune, Stromal and ESTIMATE scores were applied to evaluate association between the risk score and tumor microenvironment (TME) with “ESTIMATE” R package. Immune cell subtype infiltration values of TCGA-STAD dataset samples were based on seven algorithms: XCELL, TIMER, QUANTISEQ, MCPcounter, EPIC, CIBERSORT-ABS and CIBERSORT (21–30). The correlation between infiltration immune cell subtypes and risk scores was established through the application of Spearman correlation analysis with R packages “ggpubr.” In addition, different immune-checkpoint genes in our two risk groups were confirmed with respect to their expression by using R packages “limma.” MSI, TMB, TIDE scores and the response to CTLA4 inhibitor were performed to assess and predict immunotherapy effectiveness of our prognostic signature. IC50 (Half-maximal inhibitory concentration) was used to determine the relationship between our prognostic signature and various drug susceptibility.

GO and KEGG analyses using “clusterProfiler” R package were performed to examine the mechanisms of the signature, and GSEA analysis was used to elucidate signaling pathways.

The human immortalized normal gastric epithelial cell (GES1) and gastric adenocarcinoma cell lines (MKN-28 and HGC-27) were provided by the Chinese Academy of Sciences (Shanghai). These cells were maintained at 37°C and in 5% CO₂ in RPMI-1640 (Gibco, USA) medium with 10% fetal bovine serum (Gibco, USA) and 1% penicillin-streptomycin. Over one year, Oxaliplatin-resistant GC cell line HGC-27R was established stepwise through constant exposure to increasing quantities of Oxaliplatin and culture in RPMI-1640 (Gibco, USA) medium with oxaliplatin (1 uM).

The small interfering RNAs (siRNAs) targeting DCLK1 and scrambled siRNA nontarget control (Si-NC) were supplied by TRANSHEEPBIO (China), and their sequences are as follows: Si-DCLK1-1: sense UCUGUCGGAUAACGUGAAUUU and antisense AAAUUCACGUUAUCCGACAGA; Si-DCLK1-2: sense ACCCGAACUCUGUCGGAUAAC and antisense GUUAUCCGACAGAGUUCGGGU ( Supplementary Table 4 ). Six-well plates were seeded with GC cell at a density of 2.0 × 10⁵/well and then Namipo™ was used to assist transfection in accordance with the instructions from manufacturer.

Total RNA from cell lines was obtained with RNA Isolater Total RNA Extraction Reagent (Vazyme, Nanjing, China) in accordance with the instructions from manufacturer. Subsequently, cDNA was generated utilizing the Evo M-MLV RT Premix kit (Accurate Biotechnology). Then, qRT-PCR assay was performed using SYBR Green Master Mix (Vazyme), with beta-actin as the normalization control and the 2−ΔΔCt method for the quantification of RNA expression levels. The primer sequences (5′-3′) were as follows: DCLK1 forward, ACTTCGACGAGCGGGATAAG and reverse, GGGCCTCAAAAGATCGGAACC; β-Actin forward, CACCATTGGCAATGAGCGGTTC and reverse, AGGTCTTTGCGGATGTCCACGT ( Supplementary Table 4 ).

Cell viability and proliferation were assayed by Cell Counting Kit-8 kit (MeilunBio, Guangzhou, China) following the instructions from manufacturer. On the first day, 2 × 10³ cells per well were seeded into 96-well plates for culture. Then on day 1, 2, 3, 4 and 5 follow cell seeding, the plates were added with 10 ul of CCK-8 reagent into each of wells and incubated for 1 h at 37°C to detect cell viability. The 450-nm optical absorption values were recorded using a microplate reader (Thermo, U.S.A.) on days 1, 2, 3, 4, and 5.

Each well of a 96 well plate was seeded with 10,000 cells and treated with different concentrations of oxaliplatin (Aladdin, Shanghai) after 24 h. The Cell Counting Kit-8 (CCK8) (MeilunBio, Guangzhou, China) was employed to assess cell viability 48 h following the addition of oxaliplatin. The growth-inhibitory curves of oxaliplatin-resistant HGC-27 cells versus the parental cells were drawn, and the IC50 of oxaliplatin of each cell line was calculated.

The indicated cells (1000 cells/well) were plated into six-well plates and the plates were then cultured at 37°C for 10 to 14 days. After 15-min fixation with formaldehyde (10%), cells were stained with 1.0% crystal violet for 30 s before counting the colonies. Finally, cell clones (> 50 cells) were counted and analyzed.

For wound healing assay, cells were inoculated into 6-well plates for culture to form a confluent layer and wounds were created by scraping using sterile tips. The cells were then rinsed thoroughly with PBS and further cultured in serum-free medium. Images were acquired at 0 and 24 h with microscope and the rate of cell migration was obtained by Image J.

For apoptosis analysis, the cells were seeded into plates of six wells and then treated with oxaliplatin (Aladdin, Shanghai) for 24 h; afterward, they were digested by trypsin without addition of EDTA and washed for two times using PBS; after addition of AnnexinV- FITC and PI staining solution, they were incubated at ambient temperature in dark place for 15 min. Then, the apoptotic rate was determined with use of the flow cytometry.

We collected primary gastric cancer tissues and adjacent non-tumor tissues from 14 gastric cancer patients who underwent radical gastrectomy at The First Affiliated Hospital of Sun Yat-sen University between 2010 and 2017, and from another 14 patients who received the same procedure at The Seventh Affiliated Hospital of Sun Yat-sen University from 2021 to 2024. Clinicopathological features were gathered for patients undergoing radical gastrectomy, all of whom had not received any radiotherapy or chemotherapy before surgery. Tumor stages were classified based on the American Joint Committee on Cancer (AJCC) TNM staging system. This study received approval from the Ethics Committee of The First Affiliated Hospital of Sun Yat-sen University (No.KY-2020–024-01) and The Seventh Affiliated Hospital of Sun Yat-sen University (No.KY-2021–072-01), and all participants provided written informed consent.

Immunohistochemical analysis of DCLK1 expression was conducted in paraffin-embedded gastric cancer specimens. In brief, paraffin-embedded slices were subjected to heating at 65°C for 2 h. Tissue sections were sequentially dewaxed and rehydrated through immersion in dimethylbenzene followed by a descending ethanol series. Following antigen retrieval in a microwave oven with an EDTA solution for 15 minutes, the slices were cooled at room temperature. After washing with PBS, the slices were treated with H₂O₂ to deactivate endogenous catalase. Subsequently, the sections were blocked in 5% goat serum for 30 min and incubated with primary antibody against DCLK1 (Proteintech, 21699-1-AP) overnight at a temperature of 4 °C. The next day, the slices were incubated with secondary antibody for 30 min at room temperature. Next, we performed color development using the DAB Detection Kit (Gene Tech, #GK600511) according to the manufacturer’s instructions. The slices were photographed and recorded by the Whole Slide Scanner VS200 (Olympus). Protein expression was assessed and verified by two independent pathologists based on staining intensity and the proportion of positively stained cells.

All statistical analyses were completed using Prism (Version 10.1) and R (Version 4.2.1) software.

One-way analysis of variance (ANOVA) and unpaired Student’s t-test were adopted to assess differences between groups. A two-tailed P < 0.05 was set to determine statistical significance.

To analyze the connections between the histone acetylation-related genes and STAD, consensus clustering analysis was performed in the TCGA-STAD cohort for STAD patients. The results of the CDF curves and tracking plot showed that the optimal cluster number was K = 2 ( Figures 1A, B ). Figure 1C demonstrated that the maximum intra-group correlations with low inter-group correlations were observed when K = 2 and STAD patients could be split into two clusters according to differential expression patterns. K-M analysis was performed for comparison of overall survival advantage between two clusters. As indicated in Figure 1D , the survival probability was higher in cluster 1 (C1) than in cluster 2 (C2) (P= 0.021). Then, the immune checkpoint-related genes between the two clusters were examined with respect to their expression, and there were 12 DEGs (CD244, CSF1R, TGFBR1, BTLA, CTLA4, IL10, CD160, PDCD1LG2, KDR, CD96, CD274 and HAVCR2) between C1 and C2 ( Figure 1E ). These 12 genes showed higher expression levels in C2 ( Figure 1E ). In addition, we investigated the difference of TMB and MSI between two clusters and results indicated that C1 had higher TMB and MSI than C2 ( Figures 1F, G ).

We also analyzed different immune cell infiltration level between C1 and C2. The results in Figures 2A–C showed that C2 have higher ImmuneScore, StromalScore and ESTIMATEScore ( Figures 2A–C ). Furthermore, we found that significant differences exist in some immune cell types including CD8 T cells, Cytotoxic lymphocytes, Endothelial cells, Monocytic lineage, Myeloid dendritic cells and NK cells with a p value <0.05 for all, specifically, 0.047, 0.019, 0.036, 0.0019, 6e-05 and 0.0002, respectively, and C2 had more these types of immune cell ( Figures 2D–L ).

We screened out the differentially expressed genes between C1 and C2 through the “limma” package. Then prognostic significance of histone acetylation-related DEGs was assessed by using univariate Cox regression analyses. The 10 genes (ASCL2, ADH1B, GPR87, F13A1, HDAC11, DCLK1, GCG, FABP4, AXIN2, OGN) were retained for subsequent analysis ( Figure 3A ). Using LASSO Cox regression, a signature of eight genes was generated based on optimum λ value ( Figures 3B, C ). These 8 genes (ASCL2, GPR87, F13A1, HDAC11, DCLK1, GCG, FABP4, AXIN2) were regarded as prognostic factors for STAD, and for which the calculation of risk score was: risk score = (−0.016*ASCL2 exp.) + (0.108*GPR87 exp.) + (0.052*F13A1 exp.) + (−0.085*HDAC11 exp.) + (0.111*DCLK1 exp.) + (0.088*GCG exp.) + (0.024*FABP4 exp.) + (−0.088*AXIN2 exp.). In addition, division of STAD patients was conducted as two groups of high- and low-risk on the basis of their median risk score ( Figure 3D ). Patients with STAD in the high-risk one possessed a greater risk of death and less survival time ( Figure 3E ). The expression patterns of these eight prognostic genes between the groups were demonstrated by heatmap ( Figure 3F ). Patients with STAD in the high-risk group presented a lower overall survival (OS) than those in the low-risk group ( Figure 3G ). Figure 3H indicated that the area under the ROC curve was 0.658, 0.699 and 0.632, respectively, at 1 year, 3 years and 5 years. The PCA and t-SNE analysis suggested that patients with varying risk levels in STAD can be effectively divided into two directions ( Figures 3I, J ).

The GSE84437 dataset was used for validation of robustness of our signature. Similar to the previously described results, STAD patients in GSE84437 dataset were also assigned into two groups of low- and high-risk in light of the median risk score ( Figure 3K ). The STAD patients with different risks could be assigned into two directions, which was confirmed by PCA and t-SNE analysis ( Figures 3P, Q ). The patients showed a greater risk of death, less survival time and decreased OS in the high-risk group than those in the low-risk group ( Figures 3L, N ). The expression patterns of eight prognostic genes between the two groups were shown by heatmap ( Figure 3M ). The area under the ROC curve reached 0.632 at 1 year, 0.624 at 3 years, and 0.680 at 5 years ( Figure 3O ).

To determine the feasibility of this novel gene signature being used as a prognostic predictor independent of other factors, the Cox regression of univariate and multivariate was carried out. For the univariate analysis, we discovered that the identification of risk score as an independent prognostic predictor for STAD patients was demonstrated in both the TCGA cohort (HR = 4.015, 95% CI: 2.038–7.911, Figure 4A ) and GEO cohort (HR = 1.998, 95% CI: 1.073–3.719, Figure 4C ). For the multivariate analysis similar results were found, indicating that the risk score is an independent prognostic predictor (TCGA cohort: HR = 4.686, 95% CI: 2.223–9.875, Figure 4B ; GEO cohort: HR= 2.238, 95% CI: 1.216–4.121, Figure 4D ). A nomogram was established to predict 1, 3, and 5-year OS in TCGA-STAD cohort ( Figure 4E ). The calibration curve confirmed that the nomogram operated with an extremely high consistency with an ideal model ( Figure 4F ). The results of the ROC curves showed that the nomogram had the largest AUC (AUC=0.694) than other clinicopathological features ( Figure 4G ). Decision curves of Figure 4H also showed the nomogram presented larger net benefits.

We investigated different somatic mutation landscapes between TCGA-STAD patients defined as high- and low-risk ( Figures 5A, B ). Patients in the low-risk group had higher mutation rates than the high-risk group (94.59 vs 81.36%). In terms of gene mutation frequency, TTN, TP53 and MUC16 were the most altered gene. The most frequent mutation was missense mutation in patients with STAD. According to Figure 5C , risk score and different immune subtypes were conspicuously correlated with each other. In addition, a positive correlation was also present between the risk score and lymph node metastasis, distant metastasis, clinical stage and tumor grade except for age and T stage ( Figures 5D–I ). The K-M curves implied that high-risk patients exhibited a less favorable prognosis than those classified as low-risk in age >65, age<= 65, Female, Male, G3, M0, N1-3, Stage III-IV and T3-4 ( Figure 6 ).

The relationships between risk score and TME were determined, using calculations of Immune, Stromal and ESTIMATE scores. As shown in Figures 7A–C , the Immune, Stromal and ESTIMATE scores were higher in the high-risk group. Based on seven known algorithms, the correlation of tumor-infiltrating immune cells to the risk score was analyzed ( Figure 7D ). Since tumor-associated macrophages (TAMs) are critical in TME, we explored the relationships between our risk score and TAMs. Figure 7E implied that the level of M2 macrophages was higher in the high-risk group than in the low-risk group. Furthermore, results of functional analysis suggested that our histone acetylation-related gene signature was notably associated with immune activity and the high-risk subgroup exhibited elevated infiltrating levels of most immune cells except for Th2_cells ( Figure 7F ). Figure 7G proved that the high-risk group showed more activity for most immune pathways excepting MHC_class_I, which showed less activity. The activity of APC_co_inhibition pathway revealed no significant difference between two subgroups ( Figure 7G ).

As shown in Figure 8A , 15 immune checkpoint-related genes (HAVCR2, CD274, ADORA2A, CD96, KDR, PDCD1LG2, TGFB1, CD160, IL10, TIGHT, BTLA, TGFBR1, CSF1R, CD244 and LAG3) were conspicuously overexpressed in the high-risk group, excepting one immune checkpoint-related gene (NECTIN2), which was obviously overexpressed in the low-risk group. Then, we investigated the different MSI status (MSS, MSI-H and MSI-L) between two subgroups. The data suggested that the MSI of high-risk group is lower ( Figure 8B ). We also observed that in the high-risk group, the TMB was lower, and there was a negative correlation between risk score and TMB (R=-0.29, P=1.9e-08) ( Figures 8C, D ). Figure 8E shows that the high-risk group had remarkably high TIDE scores than the low-risk group through the use of the TIDE algorithm, indicating that the low-risk group may have more immunotherapeutic effects than the high-risk group. Furthermore, TCIA immunotherapy prediction showed that the low-risk patients had a superior response to CTLA4 inhibitor ( Figure 8F ).

IC50 was used to estimate the difference between the two risk subgroups for 12 drugs. Our findings observed that patients from the high-risk group showed a higher estimated IC50 for five drugs including Doxorubicin, Etoposide, Gefitinib, GSK690693 and KIN001–266 but showed a lower estimated IC50 for the remaining seven drugs including AMG-706, CGP-60474, Cyclopamine, Dasatinib, JNK Inhibitor VIII, PIK-93 and PD-173074 ( Figure 9 ).

To analyze the differences between the two groups regarding gene functions and enrichment, we extracted DEGs between the two groups with use of R package (“limma”). With these DEGs, GO and KEGG analyses were performed, and as shown in Figures 10A–D . The results of GO analysis showed that the DEGs were related to muscle system process, etc. in BP, collagen-containing extracellular matrix, etc. in CC, receptor ligand activity, etc. in MF ( Figures 10A, C ). According to the data from KEGG, the DEGs were abundant in cytokine-cytokine receptor interaction, signaling pathways of cGMP-PKG and Wnt ( Figures 10B, D ). GSEA analysis was used to find which pathways were enriched in the two groups and the results are shown in Figures 10E , and F . The cell adhesion molecules cams, cytokine-cytokine receptor interaction and dilated cardiomyopathy were conspicuously enriched in the high-risk group. Figure 10F showed prominently enriched arginine and proline metabolism, cell cycle and DNA replication in the low-risk group.

Among the eight histone acetylation-related genes, DCLK1 had the largest correlation coefficient in the prognostic model. Therefore, we carried out in vitro analysis to determine the functions of DCLK1 in GC cells. The RT-qPCR experimental results suggested that the expression of DCLK1 was significantly higher in MKN-45, MKN-28 and HGC-27 of GC cell lines than in the normal human gastric epithelial cell line GES-1 ( Figure 11A ). DCLK1 exhibited the highest expression levels in HGC-27 and MKN-28. Therefore, these two GC cell lines were selected for subsequent studies. To determine whether DCLK1 affected progression in gastric cancer cells, we designed siRNA to silence DCLK1 expression in MKN-28 and HGC-27 cells. The efficiency of DCLK1 silencing in MKN-28 and HGC-27 was confirmed through RT-qPCR ( Figure 11B ). Following the successful gene knockdown demonstrated previously, the CCK-8 assay was performed to validate the effects of DCLK1 on the proliferation and viability of MKN-28 and HGC-27 cells. The results proved that the proliferation of these two gastric cancer cell lines were impaired after knocking down DCLK1 compared to the Si-NC group ( Figure 11C ). Furthermore, colony formation experiments showed that colony formation was lower in the DCLK1 down-regulation group ( Figure 11D ). Figure 11E also manifested that the migration of HGC-27 was weakened in si-DCLK1 groups. Since oxaliplatin is commonly used as a first-line chemotherapeutic agent for gastric cancer, we proceeded to investigate the impact of DCLK1 expression on the effectiveness of this chemotherapy drug. To check the DCLK1 expression in Oxaliplatin (Oxa)-resistant GC cells, RT-PCR was performed (HGC-27/Oxa), and the findings manifested that DCLK1 was more expressed in Oxa-resistant GC cells (HGC-27/Oxa) than in the parental ones ( Figure 11F ). When treated with increasing levels of Oxaliplatin, HGC-27/Oxa cells exhibited greater viability ( Figure 11G ). In addition, when comparing the DCLK1 knockdown transfection group to the Si-NC group, the 50% maximal inhibitory concentration (IC50) value of oxaliplatin was lower in si-DCLK1 group ( Figure 11H ). The apoptotic rate of HGC-27 and MKN-28 cells treated with Oxaliplatin increased through the suppression of DCLK1 ( Figure 11I ). In summary, the knockdown of DCLK1 could repress gastric cancer cell in terms of cell proliferation, migration and oxaliplatin resistance. To investigate the clinical relevance of DCLK1 in GC, IHC analysis was performed on paraffin-embedded gastric cancer tissue sections from 14 patients who underwent radical gastrectomy at each of two different hospitals (The First Affiliated Hospital of Sun Yat-sen University and The Seventh Affiliated Hospital of Sun Yat-sen University). IHC staining results of tissue sections from GC patients in these two hospitals revealed that DCLK1 expression in gastric cancer tissues was higher than in normal gastric tissues ( Figures 11J, K ). In addition, elevated DCLK1 expression was associated with Borrmann’s Type, depth of invasion, lymph node metastasis, and TNM stage in GC patients ( Supplementary Table 5 ).