Construction of lactate-myeloid-related prognostic risk model and prediction of immune-excluded microenvironment in gastric cancer via integrated scRNA-Seq and bulk RNA-Seq

Introduction

Gastric cancer (GC) remains a major global health issue, ranking as the fifth most common malignancy and the fourth leading cause of cancer-related death worldwide, with approximately 1.1 million new cases and 770,000 deaths reported in 2020 (1). Although overall incidence and mortality rates are declining, attributed primarily to reduced Helicobacter pylori infection and improved food preservation practices, recent epidemiological trends show increasing incidence among younger populations, suggesting novel underlying biological drivers (1,2). GC progression involves multiple precancerous stages, including chronic atrophic gastritis, intestinal metaplasia, early gastric cancer (EGC), and advanced gastric cancer (AGC), each characterized by distinct cellular transformations and molecular events (3,4).

The tumor microenvironment (TME) critically influences cancer initiation, progression, immune evasion, and metastasis (5,6). A range of immunosuppressive mechanisms—including regulatory T-cell activity, tumor-associated macrophage (TAM) polarization, and immune checkpoint upregulation—collectively fosters immune escape and disease progression (7-9). Understanding how epithelial cell states interact with stromal and immune components across the precancer–cancer continuum is therefore essential for identifying robust biomarkers and therapeutic vulnerabilities (10,11).

Single-cell RNA sequencing (scRNA-seq) has transformed our view of GC by resolving epithelial, immune, and stromal heterogeneity at cellular resolution (12-14). However, most studies analyze isolated disease stages or focus on either epithelial or immune compartments in isolation, making it difficult to connect epithelial programs to myeloid remodeling and clinical readouts. One potentially unifying axis is lactate metabolism, particularly the SLC16A3 (MCT4)-mediated lactate export that acidifies the local milieu (15). Lactate is an active metabolite capable of shaping stromal behavior, matrix remodeling, angiogenesis, and myeloid polarization. In tumors, macrophages frequently adopt immunoregulatory/repair-like states (e.g., SPP1⁺ and APOE⁺ TAMs) associated with T-cell exclusion and poorer outcomes (16-18); yet in GC, the stage-resolved emergence of an epithelial lactate-export program and its linkage to specific myeloid states remain incompletely defined.

Another unmet need is the translation of single-cell mechanisms into parsimonious bulk signatures that generalize across platforms, stratify survival, and illuminate drug response. Such integrative models could explain why apparently “genomically quiet” tumors exhibit strong immune evasion and might nominate therapies that exploit metabolic-myeloid-extracellular matrix (ECM) dependencies.

In this study, we combined scRNA-seq with complementary bulk analyses to chart the cellular and molecular landscape from adjacent normal tissue through chronic atrophic gastritis and intestinal metaplasia to EGC and AGC. We assembled an integrated single-cell reference spanning six datasets and 116 samples, resolved stage-specific niche occupancy, and delineated epithelial states—including a late-pseudotime SLC16A3⁺ (MCT4) program—together with a macrophage-centric myeloid compartment. Building on these single-cell insights, we developed a mechanism-grounded, 19-gene bulk prognostic model anchored on SLC16A3 co-expression and macrophage immunoregulatory modules, and evaluated its prognostic performance, external validity, and independence from clinical covariates. We further interrogated risk-stratified genomics [mutation count and tumor mutational burden (TMB)], tumor immune dysfunction and exclusion (TIDE), and predicted drug response to nominate actionable pathways. By linking epithelial metabolism to macrophage-dominated immune exclusion and by delivering a clinically applicable signature, this work aimed to provide a mechanistic and translational framework for risk stratification and combination therapy design in GC. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0633/rc).

Methods

Immunohistochemistry (IHC)

Formalin-fixed, paraffin-embedded (FFPE) GC tissue sections (4 µm thick) were deparaffinized in xylene and rehydrated through a graded series of ethanol. Endogenous peroxidase activity was blocked with 3% hydrogen peroxide for 10 min. Antigen retrieval was performed using ethylenediaminetetraacetic acid (EDTA) buffer (pH 9.0) in a microwave for 15 min. After blocking with 5% bovine serum albumin for 1 hour at room temperature, the sections were incubated overnight at 4 ℃ with primary antibodies against SLC16A3 (MCT4) (ab308528, 1:100). Following TBST washes, the sections were incubated with an HRP-conjugated secondary antibody for 1 hour at room temperature. The signals were visualized using 3,3'-diaminobenzidine (DAB) substrate, and cell nuclei were counterstained with hematoxylin. Finally, the slides were dehydrated, cleared, and mounted. Whole-slide images (WSIs) were digitally scanned using a (KFBIO KF-PRO-120 digital pathology slide scanner. The cohort was listed in Table S1. Tissue samples were obtained from patients with GC from The Second Affiliated Hospital, Zhejiang University School of Medicine, between 2015 and 2017. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Institutional Review Board of The Second Affiliated Hospital, Zhejiang University School of Medicine (No. 2023.1107). Written informed consent was obtained from all patients prior to enrollment and sample collection.

Data source and preprocessing

IHC analysis

To rigorously evaluate the spatial interplay between the epithelial metabolic state and the immune microenvironment, we employed a highly quantitative digital pathology workflow. Single-cell detection and morphological feature extraction [including nuclear area, eccentricity, and optical density (OD) of hematoxylin and DAB] were performed on the WSIs using QuPath software (version 0.5.0-x64). Subsequently, empirically determined DAB OD thresholds were applied to define biomarker positivity.

Bulk RNA-seq

The Cancer Genome Atlas stomach adenocarcinoma (TCGA-STAD) was used as the development framework for model training and internal testing, and GEO GSE84437 served as the external validation cohort (19,20). Baseline demographic and available clinicopathologic characteristics of the development and external validation cohorts was listed in Table S2. Bulk transcriptomic and survival data were harmonized at the patient level before prognostic modeling. In TCGA-STAD, sample-level expression profiles were aggregated to the patient level by averaging samples sharing the same 12-character TCGA patient barcode, and survival annotations were collapsed accordingly. Samples were eligible for prognostic-model analyses if they had available bulk expression data, survival time, and survival status. Samples lacking survival time, survival status, or complete expression values for genes required for score calculation were excluded from downstream prognostic analyses. The study size was determined by all eligible samples available in the public cohorts after patient-level harmonization and exclusion of samples with missing survival annotations or required predictor values. For both cohorts, gene symbols were harmonized to Human Genome Organisation Gene Nomenclature Committee (HGNC), and low-variance features and duplicated gene symbols were removed. To support cross-platform transportability of the locked model, validation data were winsorized, quantile-mapped to the training distribution, and standardized using the mean and standard deviation estimated from the training data before risk-score calculation.

Single-cell RNA-seq

Eligible studies were required to include human GC or gastric premalignant lesion samples, provide accessible raw or processed expression matrices suitable for reanalysis, contain clear sample-level annotations such as tissue source or disease stage, and include sufficient metadata to support cell-type annotation and cross-dataset integration. Single cell data were retrieved from GSE134520 (21), GSE163558 (22), GSE167297 (23), GSE206785 (24), HRA000704 (12), dbGaP identifier phs001818.v1.p1 (25). Cells with <200 detected genes or >20% mitochondrial transcripts were removed, and genes detected in <3 cells were excluded. Post-QC matrices were concatenated into a single AnnData object. Highly variable genes were selected with sc.pp.highly_variable_genes (min_mean =0.0125, max_mean =3, min_disp =0.5), followed by regression of cell-cycle (S, G2/M), total counts, and mitochondrial content (sc.pp.regress_out). Data were scaled and clipped to a maximum of 10 (sc.pp.scale), then reduced by principal component analysis (PCA).

Batch-aware integration and clustering

To integrate datasets across platforms and batches, we used scVI (single-cell Variational Inference, python version 1.2.2) (26). Models employed two hidden layers, dropout =0.2, and layer normalization on both encoder and decoder (use_layer_norm = "both"), with batch covariates explicitly encoded. The learned latent space (X_scVI) was used to build the k-nearest-neighbor graph (sc.pp.neighbors, use_reP="X_scVI"), compute the Uniform Manifold Approximation and Projection (UMAP) embedding, and perform Leiden clustering. Clusters were annotated as distinct cell populations by inspecting canonical lineage-specific markers from the published literature.

Compositional analysis

Immune cell-type composition across sample categories was quantified and compared using density UMAP plots, stacked bar plots of cell proportions, and heatmaps of observed-to-expected counts (Ro/e).

Differential gene expression and pathway analysis

Differential expression between specified subsets across tissue compartments was assessed with the Wilcoxon rank-sum test in Scanpy (Python, version 1.9.3) (27), (sc.tl.rank_genes_groups). Signature scores were calculated with sc.tl.score_genes and the genes were listed in Table S3. Pseudotime trajectories were inferred to characterize lineage structure and differentiation dynamics using Slingshot (R, version 2.12.0) (28).

Sample-level aggregation from single cells

For sample-level analyses, we generated pseudobulk matrices by averaging log-normalized expression within a given compartment/state per sample. Signature scores (e.g., Epi_SLC16A3; TAM/monocyte/neutrophil modules) were computed per sample. Associations were evaluated by Spearman correlation and partial correlation after residualizing for batch/study and stage covariates.

Anchor-guided gene selection centered on SLC16A3

To capture the lactate-export axis, we anchored candidate-gene discovery on SLC16A3 in epithelial pseudobulk profiles. For each gene, we computed the Spearman partial correlation with SLC16A3 across samples while controlling for study/batch effects and metabolic background (hypoxia and glycolysis scores). P values were adjusted by the Benjamini-Hochberg method. We retained genes showing significant positive association with SLC16A3, bootstrap sign stability, and stage-aware differential correlation favoring AGC over EGC. This positive module was then overlapped with a curated macrophage immunoregulatory gene set, yielding the 19-gene candidate panel used for downstream multivariable prognostic modeling.

Risk model construction and validation

A gene-level penalized Cox regression pipeline was used for prognostic model development (lifelines, Python, version 0.27.8). Candidate genes were entered into a cross-validated penalized Cox framework, and model complexity was controlled through 5-fold cross-validation to favor a parsimonious multivariable solution. After gene selection, the final set of retained predictors was re-estimated in a lightly ridge-regularized Cox model (penalizer =0.05) to improve coefficient stability. The patient-level risk score was defined as the linear predictor of this final Cox model. Model performance was evaluated by Kaplan-Meier analysis with log-rank testing, time-dependent receiver operating characteristic analysis [sksurv, Python, version 0.22.2; Uno’s area under the curve (AUC) at 1, 3, and 5 years], and univariable and multivariable Cox regression adjusting for stage, age, grade, and gender. External validation was performed in GSE84437 by applying the training-derived scaling parameters and fixed model coefficients without coefficient re-estimation.

The final prognostic model comprised 19 genes: RNF144A, NMNAT1, OSBPL1A, ITGAV, IQCC, SPRY1, PNMA1, GLT8D1, ENTPD2, DNAAF3, COL1A1, SBNO2, DPH7, PPM1J, ARG1, PPP1R10, C2CD4D, ABCA7, and CUL4A. The regression coefficients of the final Cox model are provided in Table S4. For each patient, the risk score was calculated as the weighted linear predictor of the model, with higher values indicating a higher hazard of death and poorer expected overall survival. For external validation, the same locked coefficients were directly applied to the aligned GSE84437 expression matrix without re-estimation. For descriptive stratification, patients were dichotomized into low-risk and high-risk groups using the median risk score threshold defined in the development framework, and the same prespecified threshold was applied to the external validation cohort.

Somatic mutation calls and TMB

Mutation annotation format (MAF) files were parsed; non-synonymous variants (missense, nonsense, frameshift indels, splice-site, in-frame indels, translation start, nonstop) were retained. Sample-level mutation count and TMB were computed as non-synonymous count per 38 Mb exome (unless cohort-specific covered Mb was available). Group differences used Mann-Whitney tests; gene-wise high-low frequency differences used Fisher’s exact test with Benjamini-Hochberg-false discovery rate (BH-FDR) correction.

Immune-evasion metrics (TIDE)

TIDE, dysfunction, exclusion (29) were computed from normalized bulk expression using the TIDE framework with default settings. Group differences used Mann-Whitney tests; associations with the risk score used Spearman correlation.

Drug-response prediction

Expression-based predicted half maximal inhibitory concentration (IC50) values were obtained using oncoPredict (R, version 0.3.1) (30) and pRRophetic-style (R, version 0.5) (31), trained on public pharmacogenomic resources [Genomics of Drug Sensitivity in Cancer (GDSC)]. Validation expression was quantile-mapped and standardized with the training scaler before inference. For each drug, we computed Δ median IC50 (high − low) and tested high vs. low by Mann-Whitney with BH-FDR correction. Volcano plots display Δ median vs. −log10(FDR); representative drugs were visualized by box-and-jitter plots.

Statistical analysis

Unless noted, tests were two-sided. Multiple testing was controlled by Benjamini-Hochberg FDR. Effect sizes are reported as Spearman’s ρ, log-rank P, hazard ratios with 95% confidence intervals (CIs) (Cox), Δ median, and AUC as appropriate. Analyses were performed in Python and R with fixed random seeds for resampling steps.

Results

Integrated single-cell reference and stage-specific niche occupancy in GC

The full landscape of our integrated gastric single-cell atlas and stage-specific cellular niche remodeling is summarized in Figure 1. We assembled a gastric single-cell reference by integrating six datasets (116 samples) spanning paratumor (PT), chronic atrophic gastritis (CAG), moderate intestinal metaplasia (MIM), severe intestinal metaplasia (SIM), EGC, and AGC (Figure 1A). After stringent quality control and batch-aware integration, this atlas served both as a biological map and as the framework used later to derive and interpret the lactate-myeloid-related prognostic program in bulk cohorts. To assess integration quality, we compared the unintegrated PCA space with the scVI-integrated latent space. UMAP visualization colored by dataset and cell-type-stratified sample-level silhouette analysis showed reduced dataset/sample-specific segregation after integration, supporting effective cross-dataset correction while preserving biological lineage structure (Figure S1).

Figure 1 Integrated single-cell reference and stage-specific niche occupancy in gastric cancer. (A) Schematic of the analysis pipeline: tissue sampling → single-cell profiling → batch-aware integration and annotation → linkage to bulk cohorts for model training. (B) UMAP of 355,091 cells from 116 samples showing 18 lineages (legend at right) spanning epithelial, lymphoid, myeloid and stromal compartments. (C) Stage-wise density projection of PT, CAG, MIM, SIM, EGC and AGC onto the unified UMAP; red = high density. PT concentrates in normal epithelium and T cells; CAG/MIM show immune contraction with epithelial replacement; SIM displays an expansion of aberrant epithelium; EGC/AGC show increasing occupancy of abnormal epithelium together with myeloid/lymphoid territories. (D) Dot plot of canonical markers validating lineage annotations (dot size = fraction of cells expressing, color = mean expression within lineage). (E) Stacked bar chart of lineage proportions by clinical group, illustrating the rise of epithelium from CAG → SIM and the re-emergence of myeloid/lymphoid in EGC/AGC. (F) Heatmap summarizing stage-wise shifts in lineage abundance (color scale reflects relative enrichment), highlighting the late-stage increase of macrophages/monocytes/neutrophils and aberrant epithelium. AGC, advanced gastric cancer; CAG, chronic atrophic gastritis; DC, dendritic cell; EGC, early gastric cancer; GEO, Gene Expression Omnibus; MIM, moderate intestinal metaplasia; NK, natural killer; PT, paratumor; SIM, severe intestinal metaplasia; TCGA, The Cancer Genome Atlas; UMAP, Uniform Manifold Approximation and Projection.

A total of 355,091 cells resolved on a common UMAP into 18 lineages that span the four major compartments—epithelial, lymphoid, myeloid, and stromal. Lineages included gastric epithelial subtypes (chief, parietal, mucous neck, pit/aberrant epithelium), immune populations [T, B, natural killer (NK), plasma, dendritic cell, monocyte/macrophage, neutrophil, mast cell], and stromal/vascular cells (fibroblast, smooth muscle, endothelial) (Figure 1B). Annotation was supported by canonical markers summarized in the dot plot (Figure 1D), where circle size reflects the fraction of cells expressing a marker and color encodes mean expression, providing orthogonal validation of lineage identity across datasets.

To visualize disease-stage niches without re-embedding, we projected each clinical group’s cells onto the integrated UMAP and computed normalized 2D density (Figure 1C). The maps reveal a stepwise reallocation of niche occupancy along progression: PT concentrates within homeostatic epithelium and T-cell territories, consistent with a largely immune-surveilled, non-dysplastic environment. CAG/MIM show a contraction of lymphoid niches with compensatory epithelial expansion, indicating immune attenuation alongside metaplastic remodeling. SIM exhibits a pronounced bloom of aberrant epithelium, marking the inflection toward dysplasia. EGC/AGC progressively populate aberrant epithelial neighborhoods and myeloid-rich territories, with AGC showing the broadest epithelial coverage and substantial myeloid infiltration, consistent with an inflamed yet immune-dysregulated carcinoma microenvironment.

Stacked-bar compositions (Figure 1E) quantify these shifts: epithelial fractions increase from CAG to SIM, while myeloid and, to a lesser extent, lymphoid compartments re-emerge in EGC/AGC. An effect heatmap using the observed/expected ratio (Ro/e) by lineage and stage (Figure 1F) highlights the same pattern from a statistical perspective—macrophage/monocyte/neutrophil lineages and aberrant epithelium are enriched toward late disease (Ro/e >1), whereas several homeostatic epithelial subsets contract (Ro/e <1) as dysplasia and invasion take hold.

Together, this integrated atlas (I) harmonizes cellular identities across platforms, (II) captures stage-specific niche migration from normal epithelium and T-cell zones toward aberrant epithelium with myeloid predominance, and (III) establishes a reference scaffold for linking bulk-level risk signals to precise cellular contexts. In subsequent analyses, this map underpins the lactate-myeloid/ECM interpretation of the High-risk state and provides lineage-resolved anchors for cross-cohort validation.

Epithelial states along progression and the signaling profile of the SLC16A3⁺ program

Figure 2 characterizes the transcriptional heterogeneity of aberrant epithelium and its functional signaling programs that drive microenvironmental reprogramming across gastric lesion progression. To dissect epithelial drivers of these microenvironmental shifts, we sub-clustered aberrant epithelium from SIM/EGC/AGC and identified five transcriptional states (Epi_ANPEP, Epi_LLGL2, Epi_OLFM4, Epi_SLC16A3, Epi_SPINK4) (Figure 2A). Marker-by-cluster dot plots (Figure 2B) confirmed state identity and purity: each cluster shows a characteristic marker peak with limited spillover, supporting discrete epithelial phenotypes rather than a single continuous gradient. Across disease stages, Epi_SLC16A3 (harboring SLC16A3/MCT4, a canonical lactate-exporter) expanded progressively from SIM to EGC to AGC, whereas homeostatic states (Epi_ANPEP, Epi_OLFM4, Epi_SPINK4) contracted (Figure 2C). A lineage-by-stage enrichment summary (Figure 2D) recapitulates these trends statistically: Epi_SLC16A3 shows positive enrichment toward late disease, while several homeostatic compartments display negative enrichment. These dynamics place Epi_SLC16A3 on the main trajectory of malignant progression and suggest a metabolic shift toward glycolytic/lactate-export programs as histologic severity increases. A curated ligand–receptor screen revealed that Epi_SLC16A3 is selectively enriched for chemokine outputs spanning CXCL families (Figure 2F), together with signal-readout capacity via CXCR-class receptors (Figure 2E), suggesting autocrine tuning and paracrine recruitment of CXCR-bearing myeloid lineages. This chemokine-rich signature provides a mechanistic link between SLC16A3⁺ epithelium and the myeloid-dominant, immune-excluded niches observed in high-risk tumors. Epi_SLC16A3 further displays elevated growth-factor programs (VEGF/EGF/FGF families) and cognate receptor tyrosine kinases (e.g., MET, EGFR/ERBBs, FGFRs), consistent with proliferative and pro-angiogenic signaling. The co-occurrence of VEGF ligands and VEGF receptor pathways is congruent with vascular remodeling and supports microenvironmental conditions permissive to myeloid influx and stromal activation. Checkpoint-focused panels indicate that Epi_SLC16A3 co-expresses a mixed immune-modulatory repertoire comprising co-stimulatory (TNFSF/TNFRSF) and co-inhibitory ligands/receptors (e.g., PD-L1/PD-L2 classes and related axes). Finally, cell-adhesion modules (integrin/selectin/cadherin families) and antigen-presentation genes (HLA class I and inducible class II components) are comparatively upshifted in Epi_SLC16A3, indicative of (I) enhanced epithelial-stromal and epithelial-immune contacts and (II) a cytokine-responsive epithelial state capable of presenting antigens under inflammatory cues. Together, these features suggest that the SLC16A3⁺ program is not merely metabolic but signaling-competent, able to shape immune cell trafficking, retention, and tone. Collectively, Epi_SLC16A3 marks a late-stage, lactate-exporting epithelial program that expands with disease and is molecularly wired to attract and condition myeloid cells, to stimulate angiogenic/RTK pathways, and to modulate checkpoint tone while reinforcing tissue contacts and inducible antigen presentation. This state provides a cellular and signaling scaffold for the myeloid/ECM-dominated, T-cell-excluded microenvironment associated with high-risk tumors, mechanistically connecting epithelial metabolism to immune evasion and progression.

Figure 2 Epithelial states along progression and the signaling profile of the SLC16A3⁺ program. (A) UMAP of integrated single cells showing five epithelial states (Epi_ANPEP, Epi_LLGL2, Epi_OLFM4, Epi_SLC16A3, Epi_SPINK4). (B) Marker dot plot for the five states (dot size: fraction of cells in the state expressing the marker; color: mean expression). (C) Stacked bar showing the proportion of epithelial states across SIM, EGC, and AGC. (D) Heatmap summarizing stage-wise shifts in epithelial states abundance (color scale reflects relative enrichment). (E) Module-level heatmaps (row = epithelial state; column = gene family members) summarizing chemokine receptors. (F) Module-level heatmaps (row = epithelial state; column = gene family members) summarizing chemokines (G) Module-level heatmaps (row = epithelial state; column = gene family members) summarizing growth factors and growth-factor receptors. (H) Module-level heatmaps (row = epithelial state; column = gene family members) summarizing co-stimulatory/co-inhibitory ligands and receptors. (I) Module-level heatmaps (row = epithelial state; column = gene family members) summarizing cell-adhesion molecules, and antigen-presentation genes. AGC, advanced gastric cancer; CAG, chronic atrophic gastritis; EGC, early gastric cancer; MIM, moderate intestinal metaplasia; PT, paratumor; Ro/e, observed/expected ratio; SIM, severe intestinal metaplasia; UMAP, Uniform Manifold Approximation and Projection.

Epithelial trajectory places the SLC16A3⁺ state at late pseudotime and ties it to ECM/angiogenic and myeloid-recruiting programs

Using pseudotime inference on aberrant epithelium from SIM/EGC/AGC, five epithelial states—Epi_ANPEP, Epi_LLGL2, Epi_OLFM4, Epi_SLC16A3, and Epi_SPINK4—were arranged along a continuous manifold. The ordering places Epi_ANPEP/LLG2/OLFM4 at early to intermediate positions and Epi_SLC16A3 at late pseudotime, indicating that the SLC16A3⁺ (MCT4/lactate-export) program emerges toward the terminal, progression-linked end of the trajectory (Figure 3A,3B). To validate the clinical relevance and spatial context of this terminal metabolic state, we evaluated epithelial SLC16A3 expression in situ using IHC. Strikingly, quantitative assessment across a patient cohort (n=9) revealed that SLC16A3 positivity was markedly enriched in Lauren intestinal and mixed subtypes compared to the diffused subtype (Figure 3C). This subtype-specific enrichment strongly suggests that the dense, cohesive glandular architecture typical of intestinal and mixed tumors creates localized hypoxic bottlenecks. This physical constraint exerts severe selective pressure, driving the aberrant epithelium toward an MCT4-dependent metabolic climax.

Figure 3 Trajectory of aberrant epithelium and program dynamics highlight the late SLC16A3⁺ state. (A) UMAP of aberrant epithelial cells from SIM/EGC/AGC showing five states (Epi_ANPEP, Epi_LLGL2, Epi_OLFM4, Epi_SLC16A3, Epi_SPINK4) with inferred trajectories (black curves). (B) Pseudotime map (color scale) placing Epi_SLC16A3 at late pseudotime, whereas Epi_ANPEP/LLGL2/OLFM4 occupy earlier positions. (C) Left: representative IHC staining images illustrating the spatial distribution of epithelial SLC16A3 (MCT4) in gastric cancer tissues. Scale bar: 250 µm; right: quantitative comparison of SLC16A3 positivity rates among gastric cancer patients, stratified by Lauren classification subtypes (Lauren intestinal/mixed vs. Lauren diffused, n=9). (D) Heatmap of pathway scores across stages (SIM, EGC, AGC). (E) Representative violin plots of selected pathways showing stage-dependent shifts (x-axis: SIM → EGC → AGC). (F) Pseudotime trends (loess) for key programs demonstrate rising ECM/angiogenesis and chemokine activity along the trajectory, with the SLC16A3⁺ branch displaying the strongest increase; intestinal/secretory modules decrease or flatten. AGC, advanced gastric cancer; ECM, extracellular matrix; EGC, early gastric cancer; IHC, immunohistochemistry; SIM, severe intestinal metaplasia; UMAP, Uniform Manifold Approximation and Projection.

Consistent with this phenotypic climax, gene-set scoring revealed coordinated functional shifts across stages. In AGC, intestinal/secretory and cell-polarity/epithelial differentiation programs were attenuated (32-43), whereas pathways related to ECM remodeling/cell adhesion (collagen/integrin/MMP modules), growth-factor signaling (TGF-β/VEGF/FGF axes), and myeloid-recruiting chemokines (CCL/CXCL families) were preferentially elevated (7,8,44-59) (Figure 3D). Stage-wise violin plots (Figure 3E) show the same pattern at single-cell resolution, with progressive up-shift of matrix, angiogenic, and chemokine signatures from SIM to EGC to AGC, accompanied by depletion of homeostatic epithelial features.

Smoothing signature scores along pseudotime demonstrates monotonic or biphasic increases in ECM/adhesion, VEGF/FGF/TGF-β signaling, and chemokine programs as cells transition toward late states (Figure 3F). The Epi_SLC16A3 branch shows the steepest late-pseudotime up-trend across these modules, while differentiation and secretory signatures decline early and remain low thereafter. These dynamics position SLC16A3⁺ epithelium as a metabolically rewired, pro-angiogenic, matrix-remodeling, chemokine-competent state.

Together, the trajectory and program analyses indicate that the SLC16A3⁺ (MCT4) epithelial state marks a late, progression-linked phenotype transcriptionally wired to (I) remodel ECM, (II) activate VEGF/FGF/TGF-β growth signaling, and (III) recruit/retain myeloid cells via chemokine outputs. This provides a mechanistic bridge from epithelial metabolic rewiring to the myeloid-rich, T-cell-excluded microenvironment that characterizes the high-risk tumor state in bulk analyses.

Macrophage diversity resolves into anti-inflammatory/angiogenic versus antigen-presenting programs that shift with progression

Figure 4 delineates the myeloid compartment architecture, subset-specific molecular signatures, and their dynamic association with the SLC16A3⁺ epithelial program along gastric lesion progression. To resolve the myeloid compartment, we delineated the myeloid compartment into monocyte, macrophage, and neutrophil subsets using canonical markers. Classical/intermediate monocytes expressed CD14, CCR2, VCAN, LST1, FCN1 with variable HLA-II and S100A8/A9 (Figure 4A). Neutrophils were identified by CXCL8, S100A8/A9 and granule-associated genes, alongside subclusters demarcated by PTMA and TMSB4X (Figure 4C). In contrast, macrophages segregated into four reproducible states—APOE⁺, FOLR2⁺, HLA-DRA⁺, and SPP1⁺—each supported by state-specific marker expression (Figure 4B). Sample-level correlation of cell-type abundances revealed a tight co-clustering of macrophage states with SLC16A3⁺ epithelium (red box), whereas neutrophils and classical monocytes showed weaker coupling (Figure 4D). This positioning implicates lactate-exporting epithelial neighborhoods as preferential niches for macrophage accumulation.

Figure 4 Myeloid decomposition reveals a macrophage-centric program linked to epithelial lactate export and disease progression. (A) Dot plot of monocyte subclusters (rows) showing selected marker genes (columns). Dot size indicates the fraction of cells in the subcluster expressing the gene; dot color indicates mean expression (scaled within gene). Subclusters include CD14⁺, CD14⁺CD16⁺, CD16⁺, CXCL8⁺, and LGALS3⁺ monocytes; markers include antigen-presentation (e.g., HLA-DRA), inflammatory (S100A8/A9, CXCL8), scavenger/adhesion (LGALS3, VCAN), and trafficking/chemokine genes (e.g., CCR2, CCL2). (B) Dot plot of macrophage subclusters: APOE⁺, FOLR2⁺, HLA-DRA⁺, and SPP1⁺ macrophages, with canonical markers (e.g., APOE, TREM2, FOLR2, SPP1, LILRB1, HLA-DRA, CD68). (C) Dot plot of neutrophil subclusters (e.g., CXCL8⁺, PTMA⁺, TMSB4X⁺, S100A8⁺ neutrophils) highlighting lineage-specific programs. (D) Sample-level co-variation heat map (Spearman correlations) of cell-type abundances across the atlas. Hierarchical clustering groups related states; the red box highlights a positive block linking SPP1⁺/APOE⁺ macrophages with the SLC16A3⁺ (MCT4) epithelial state, consistent with a macrophage-dominated niche emerging alongside epithelial lactate-export programs. (E) Stage-resolved enrichment of macrophage subclusters across PT → CAG → MIM → SIM → EGC → AGC. Color intensity reflects effect size (scaled change in fraction relative to PT). (F) Heat maps of chemokine receptors (left) and chemokines (right, two panels) across macrophage subclusters (rows). SPP1⁺/APOE⁺ tumor-associated macrophages preferentially express CCR/CXCR receptor and are associated with ligands such as CCL and CXCL family members, consistent with recruitment/retention and ECM-angiogenic coupling. Values are mean expression scaled within gene. (G) Violin plots showing the pathway score including anti-inflammatory score and pro-inflammatory score of SPP1⁺ and APOE⁺ macrophages across disease stages (PT, CAG, MIM, SIM, EGC, AGC). AGC, advanced gastric cancer; CAG, chronic atrophic gastritis; ECM, extracellular matrix; EGC, early gastric cancer; MIM, moderate intestinal metaplasia; PT, paratumor; Ro/e, observed/expected ratio; SIM, severe intestinal metaplasia.

Across the clinical trajectory (PT → CAG → MIM → SIM → EGC → AGC), macrophages increased with the strongest enrichment in AGC, driven predominantly by SPP1⁺ and APOE⁺ states (Figure 4E). FOLR2⁺ macrophages were present but less expanded late, while HLA-DRA⁺ macrophages—consistent with antigen-presenting/IFN-γ-responsive programs—did not show the same late-stage amplification. These trends refine “myeloid infiltration” in advanced disease to a macrophage-dominated phenomenon.

Ligand-receptor/module screens showed that SPP1⁺ and APOE⁺ macrophages up-regulate CCR/CXCR-class receptors together with chemokines from the CCL/CXCL families, consistent with recruitment and retention loops (e.g., CCL2-CCR2, CCL5-CCR5, CXCL12-CXCR4 axes) (Figure 4F). This paracrine architecture aligns with epithelial outputs mapped in SLC16A3⁺ neighborhoods and provides a plausible mechanism for macrophage crowding in those regions.

State-specific scoring indicated that SPP1⁺ and APOE⁺ macrophages preferentially express anti-inflammatory/angiogenic and tissue-remodeling programs—including lipid-handling and efferocytosis modules (e.g., APOE/MERTK/AXL/GAS6) and pro-angiogenic mediators (e.g., SPP1, VEGF-axis components)—with scores rising from SIM to EGC/AGC. Inflammatory scores were comparatively lower in these states, whereas HLA-DRA⁺ macrophages retained the highest HLA-II/co-stimulatory features, consistent with antigen presentation rather than pro-angiogenic remodeling (Figure 4G).

Together, these data pinpoint macrophages—especially SPP1⁺ and APOE⁺ TAM states—as the principal myeloid partners of the SLC16A3⁺ (MCT4) epithelial program, linking epithelial lactate-export programs with a chemokine-rich, immuno-modulatory, and angiogenic niche. This macrophage bias, rather than neutrophil or classical monocyte dominance, provides a cellular explanation for the myeloid/ECM-rich, T-cell-excluded microenvironment that characterizes high-risk tumors in bulk analyses.

Development and external validation of a 19-gene prognostic model reflecting an epithelial lactate-export-associated myeloid anti-inflammatory program

After preprocessing and patient-level harmonization, 355 The Cancer Genome Atlas stomach adenocarcinoma (TCGA-STAD) patients were retained within the development framework, including 147 outcome events (41.4%) and 208 censored observations. The median observed time-to-event or censoring in the development framework was 476 days. For external validation, 343 GSE84437 patients had matched expression and clinical data and were included in downstream analyses; among them, 207 patients (60.3%) experienced the event of interest, and the median overall survival time was 2,100.2 days. Baseline demographic and available clinicopathologic characteristics of the development and external validation cohorts were listed in Table S2.

To capture the lactate-export axis, we anchored on SLC16A3 (MCT4) and computed Spearman partial correlations on epithelial pseudobulk profiles while controlling for study/batch and metabolic background (hypoxia and glycolysis). This yielded a broad, statistically significant set of positively co-expressed genes and a smaller inverse module. We intersected the positive set with a curated macrophage anti-inflammatory/angiogenic program (e.g., MRC1, CD163, TREM2, FCN1, SBNO2, F13A1), then applied bootstrap stability selection and stage-aware differential correlation (AGC vs. EGC) to prioritize robust, progression-linked features. The resulting 19-gene core—RNF144A, NMNAT1, OSBPL1A, ITGAV, IQCC, SPRY1, PNMA1, GLT8D1, ENTPD2, DNAAF3, COL1A1, SBNO2, DPH7, PPM1J, ARG1, PPP1R10, C2CD4D, ABCA7, and CUL4A—served as candidate predictors for multivariable prognostic modeling.

We next translated these candidate predictors into a gene-level prognostic model using a penalized Cox regression pipeline with 5-fold cross-validation. During model development, penalization was used to shrink unstable coefficients and favor a parsimonious multivariable solution, after which the retained genes were re-estimated in a lightly ridge-regularized Cox model to improve coefficient stability. The coefficient path (Figure 5A) illustrates shrinkage of unstable predictors as penalization increased, while the cross-validation curve (C-index versus log penalty; Figure 5B) showed a broad near-optimal region supporting robust internal performance. The final 19-gene model was selected from this cross-validated procedure and carried forward as a locked risk-score formula for subsequent analyses.

Figure 5 Construction and validation of the 19-gene prognostic model in gastric cancer. (A) Penalized Cox coefficient paths during model selection across log(λ), showing progressive shrinkage of unstable predictors. (B) Five-fold cross-validated performance versus λ; the selected penalty balanced predictive performance and model parsimony. (C) Internal TCGA test-set risk-score distribution (top) and survival scatter (bottom) sorted by risk; the vertical line indicates the prespecified median risk-score threshold. (D) TCGA test set Kaplan-Meier curves with numbers at risk; high-risk patients show significantly worse survival (log-rank P=8.86×10⁻⁸). (E) TCGA test set expression heatmap of the 19 genes ordered by risk score, illustrating a consistent expression gradient. (F) TCGA test set time-dependent ROC at 1/3/5 years (AUC 0.695/0.744/0.757). (G) Univariable Cox forest plot (TCGA test set): risk_z strongly associated with poor survival (P~1.62×10⁻¹⁷). (H) Multivariable Cox forest plot (TCGA test set): the standardized risk score remains an independent prognostic factor after adjustment for stage, age, gender and grade (risk_z P=4.79×10⁻¹⁰; stage and age significant; grade and gender not). (I) GEO validation (GSE84437) risk score distribution and survival scatter, showing the same monotonic pattern. (J) GEO validation (GSE84437) Kaplan-Meier curves; high-risk patients have inferior outcomes (log-rank P=6.68×10⁻³). AUC, area under the curve; GEO, Gene Expression Omnibus; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Applied within the TCGA-STAD development framework, the model-derived risk score displayed a monotonic gradient across patients and co-segregated with deaths, with event barcodes enriched toward the high-risk end (Figure 5C). Kaplan-Meier analysis showed a marked separation between high- and low-risk groups with a significant log-rank P value (Figure 5D). A heatmap ordered by risk score demonstrated a coherent, patient-level expression gradient across the 19 core genes (Figure 5E), consistent with coordinated pathway activation rather than isolated markers. Time-dependent ROC curves yielded AUCs of 0.695/0.744/0.757 at 1/3/5 years, respectively (Figure 5F), indicating durable discrimination over clinically relevant horizons.

In univariable Cox analyses, the standardized continuous risk score derived from the locked 19-gene model (risk_z) was strongly associated with overall survival (Figure 5G). In multivariable Cox models adjusting for stage, age, grade, and gender, risk_z remained an independent predictor of poor outcome, alongside stage and age (grade and gender were not significant; Figure 5H). These results indicate that the model contributes prognostic information beyond routine clinicopathologic variables.

Cross-platform evaluation reproduced the expected monotonic risk gradient and event enrichment (Figure 5I), and yielded a significant KM separation between risk strata (Figure 5J), supporting the robustness and transferability of the model after cross-study alignment.

As an unbiased benchmark, we additionally constructed a common-gene whole-transcriptome agnostic LASSO-Cox model without predefined SLC16A3, lactate-related, or macrophage-related anchors. The agnostic model retained 17 genes and significantly stratified survival in both TCGA-STAD and GSE84437 (Figure S2), indicating that the prognostic signal was not solely dependent on the SLC16A3-centered feature-selection strategy. Meanwhile, the lactate-myeloid 19-gene signature retained a clear mechanism-guided interpretation of the epithelial lactate-export and macrophage-rich immune-exclusion axis.

The 19-gene panel is anchored in epithelial SLC16A3/MCT4 and enriched for macrophage-tolerogenic/angiogenic programs (e.g., ARG1, COL1A1, ABCA7), coherently linking lactate export, matrix/RTK signaling, and myeloid conditioning. Its prognostic power in TCGA and GEO, and independence from Stage/Age, indicate that it captures a microenvironment-coupled risk axis that is not reducible to conventional clinicopathologic staging. In conjunction with our single-cell maps, these results position the signature as a readout of a lactate-myeloid/ECM-dominated, T-cell-excluded niche, offering mechanistic traction for risk stratification and therapeutic hypothesis generation.

Risk-stratified genomic and immune-evasion landscape of GC

We next compared somatic alterations between risk strata. In the low-risk cohort, 168/174 (96.6%) tumors harbored at least one nonsynonymous mutation, versus 148/173 (85.5%) in the high-risk cohort (Figure 6A,6B). Canonical gastric-cancer genes—including TTN, TP53, MUC16, LRP1B, ARID1A, SYNE1, CSMD3, FAT4, FLG, PCLO, OBSCN, HMCN1, SPTA1, DNAH5, ZFHX4, RYR2, ACVR2A, PCDH15, and PIK3CA—showed higher mutation frequencies in low-risk than high-risk tumors (Figure 6C). A gene-wise volcano plot of Δfrequency (high − low) displayed a global left shift (most genes Δ<0), consistent with lower alteration rates in the high-risk group; after multiple-testing correction, few genes reached FDR significance (Figure 6D).

Figure 6 Somatic mutation patterns and mutational burden by risk group. (A) Co-oncoplots for the low-risk (168/174 altered, 96.6%) cohort. Rows: top mutated genes ordered by overall frequency; columns: patients within each risk group. Cell colors indicate mutation types (Missense, Nonsense, Frame_Shift_Del/Ins, In_Frame_Del/Ins, Splice_Site; black = Multi_Hit). Upper bars show per-sample mutation counts (proxy TMB). Right bars display per-gene mutation frequency. Bottom tracks indicate risk group. (B) Co-oncoplots for the high-risk (148/173 altered, 85.5%) cohort. (C) Bar chart of mutation frequencies for the top genes in each group (red = high, blue = low), showing consistently higher frequencies in low-risk for multiple genes (TTN, TP53, MUC16, LRP1B, ARID1A, SYNE1, CSMD3, FAT4, FLG, PCLO, HMCN1, DNAH5, OBSCN, RYR2, ZFHX4, ACVR2A, PCDH15, PIK3CA). (D) Volcano plot of Δ mutation frequency (high − low) versus −log10 FDR (Fisher’s exact test). The dashed vertical line marks Δ=0; the dashed horizontal line marks FDR=0.05. The predominant left-shift (Δ<0) indicates lower mutation prevalence in HIGH-risk; few genes cross the FDR threshold. (E) Groupwise comparisons of non-synonymous mutation count between low- and high-risk tumors. Box-and-jitter plots (dashed boxes, no outliers shown) with Mann-Whitney tests (both P<0.001). (F) Groupwise comparisons of TMB (mut/Mb) between low- and high-risk tumors. Box-and-jitter plots with Mann-Whitney tests (both P<0.001). FDR, false discovery rate; RS, risk score; TMB, tumor mutational burden.

Consistent with these frequency patterns, the low-risk group harbored significantly higher nonsynonymous mutation counts and TMB (mutations per mega base) than the high-risk group (Mann-Whitney U, P<0.001 for both; Figure 6E,6F). Thus, low-risk tumors are more mutationally active, whereas high-risk tumors are comparatively genomically quieter—aligning with our lactate-myeloid/ECM axis, in which the high-risk state reflects immune exclusion rather than hyper-mutated immunogenicity.

To quantify immune evasion, we applied TIDE (29), a transcriptome-based surrogate that integrates two mechanisms: (I) T-cell dysfunction in T-cell-inflamed tumors and (II) T-cell exclusion in non-inflamed tumors; higher TIDE scores indicate transcriptome-inferred immune-evasion features and a lower predicted likelihood of response to immune checkpoint blockade (Figure 7). Relative to low-risk, high-risk tumors exhibited higher TIDE (Figure 7A), greater exclusion (Figure 7C) and greater dysfunction (Figure 7E), but a lower microsatellite instability (MSI) expression signature (a transcriptomic surrogate of MSI-H status) (Figure 7G). These patterns argue against an MSI-H, hyper-mutated phenotype in the High-risk state and instead support a myeloid/ECM-dominated, T-cell-excluded microenvironment.

Figure 7 TIDE-defined immune evasion by risk group and its association with the risk score. (A) Box-and-jitter plots comparing low RS (blue) and high RS (red) for TIDE. Dashed boxes mark distribution cores; two-sided Mann-Whitney tests are reported above each panel (P<0.001). (B) Patient-level correlations between the risk score and TIDE components with least-squares trend lines: TIDE (ρ=0.231, P=1.12×10⁻⁵). (C) Box-and-jitter plots comparing low RS (blue) and high RS (red) for T-cell exclusion. Dashed boxes mark distribution cores; two-sided Mann-Whitney tests are reported above each panel (P<0.001). (D) Patient-level correlations between the risk score and T-cell exclusion (ρ=0.467, P=1.27×10⁻²⁰). (E) Box-and-jitter plots comparing low RS (blue) and high RS (red) for T-cell dysfunction. Dashed boxes mark distribution cores; two-sided Mann-Whitney tests are reported above each panel (P<0.001). (F) Patient-level correlations between the risk score and dysfunction (ρ=0.24, P=4.76×10⁻⁶). (G) Box-and-jitter plots comparing low RS (blue) and high RS (red) for MSI expression signature. Dashed boxes mark distribution cores; two-sided Mann-Whitney tests are reported above each panel (P=0.0181). (H) Patient-level correlations between the risk score and MSI expression signature components with least-squares trend lines (ρ=−0.184, P=4.86×10⁻⁴). TIDE is a transcriptomic metric of tumor immune evasion that models two mechanisms—T-cell dysfunction in inflamed tumors and T-cell exclusion in non-inflamed tumors; higher values denote stronger immune escape and lower probability of response to immune checkpoint inhibitors. MSI, microsatellite instability; RS, risk score; TIDE, tumor immune dysfunction and exclusion.

At the patient level, the risk score correlated positively with TIDE (Figure 7B) and dysfunction (Figure 7F), with the strongest association for exclusion (Figure 7D), correlated negatively with MSI-signature (Figure 7H). These orthogonal metrics support a coherent interpretation that the adverse prognosis captured by our signature may be associated with a T-cell-exclusion phenotype, consistent with a lactate-myeloid/ECM axis in high-risk tumors rather than heightened tumor immunogenicity.

Expression-based drug-response prediction across risk strata

We contrasted predicted drug response between high- and low-risk tumors using a volcano-style analysis in which the x-axis is the difference in median drug response (Δ median = high − low, reported as AUC/viability area) and the y-axis is −log10(FDR) from a two-group ANOVA (Figure 8A). Vertical and horizontal dashed lines mark no effect (Δ=0) and the significance threshold (FDR =0.05), respectively. Points above the horizontal line (colored red) denote compounds with FDR-significant differential activity; gray points are non-significant. Because lower AUC implies predicted greater sensitivity, negative Δ (left of zero) indicates predicted greater sensitivity in the high-risk group, whereas positive Δ (right of zero) indicates predicted greater sensitivity in the low-risk group.

Figure 8 Risk-stratified pharmacogenomic landscape in gastric cancer. (A) Volcano plot contrasting predicted drug response between high- and low-risk tumors. The x-axis shows the difference in group medians (Δ median = high − low) of predicted response AUC/viability area (lower values = predicted greater sensitivity). The y-axis shows −log10(FDR) from Mann-Whitney U test. Vertical and horizontal dashed lines denote no effect (Δ=0) and the significance threshold (FDR =0.05), respectively. Red points are FDR-significant; grey points are non-significant. By definition, Δ<0 (left) indicates predicted greater sensitivity in high-risk, and δ>0 (right) indicates predicted greater sensitivity in low-risk. Most significant hits cluster at Δ>0, indicating broader predicted sensitivity in the low-risk stratum; selected top hits are annotated. A limited number of compounds fall at Δ<0, nominating candidate High-risk-associated signal. (B) Per-compound distributions (jittered dots with boxplots) for representative agents highlighted in (A). Blue = low-risk; red = high-risk. Boxes show median and interquartile range; whiskers denote 1.5× IQR. Consistent with ERBB/EGFR inhibitors—lapatinib, gefitinib, erlotinib, osimertinib, sapitinib—and OSI-027 (mTORC1/2) show lower median AUC in low-risk (greater predicted sensitivity), whereas AZD8055 (mTORC1/2) shows lower median AUC in high-risk (high-risk-selective). Additional examples include BMS-754807 (IGF-1R/IR), GSK1904529A (IGF-1R), ulixertinib (ERK), NU7441 (DNA-PK), ML323 (USP1), MK-1775 (WEE1), dihydroorotenone (mitochondrial complex II), AZD3759 (EGFR), LGK974 (porcupine/WNT), afatinib (ERBB), and BMS-345541 (IKK). All two-sided tests yielded P<0.001 (Mann-Whitney U test). AUC, area under the curve; FDR, false discovery rate; IQR, interquartile range.

Among significant hits, the distribution is right-skewed, with many red points clustering at Δ>0, indicating that a broad set of compounds shows computationally inferred greater activity in the low-risk stratum. These agents occupy the upper-right quadrant (large effect size with strong statistical support), and several top candidates are annotated in the plot. In contrast, fewer significant agents lie at Δ<0 (upper-left quadrant), indicating a smaller subset predicted to show lower AUC in the high-risk stratum; a handful of compounds display notably negative Δ values, consistent with pronounced high-risk selectivity.

For representative compounds highlighted in the volcano plot, we visualized patient-level predicted responses as jittered dot-and-box plots (blue, low-risk; red, high-risk) (Figure 8B). Most EGFR/ERBB inhibitors—including lapatinib, gefitinib, erlotinib, osimertinib, and sapitinib—show an upward shift of the High-risk group (higher median AUC), indicating greater predicted sensitivity in Low-risk tumors. The mTOR axis shows a similar pattern for OSI-027 (mTORC1/2), whereas AZD8055 exhibits the opposite direction (lower median AUC in high-risk), marking a candidate High-risk-associated signal. Additional agents [BMS-754807 (IGF-1R/IR), GSK1904529A (IGF-1R), ulixertinib (ERK), NU7441 (DNA-PK), ML323 (USP1/USP1-UAF1), MK-1775 (WEE1), dihydroorotenone (mitochondrial complex II)] also display Δ>0 with a lower central tendency in Low-risk, reinforcing Low-risk predominance among significant hits. In contrast, Acetalax shows higher AUC in High-risk (Δ>0 overall), again indicating relative sensitivity in Low-risk tumors.

Across panels, median separation rather than isolated outliers underpins most significant results, while the vertical spread within groups underscores inter-patient heterogeneity. Together with Panel A, these distributions show that—although the High-risk state is immune-excluded and myeloid/ECM-dominated—it showed only a limited number of candidate predicted drug-response signals, including AZD8055, whereas Low-risk tumors exhibited broader computationally predicted sensitivity patterns, including to ERBB and mTOR pathway inhibitors.

Discussion

In this study we combined single-cell and bulk transcriptomic analyses to chart the stepwise remodeling of the gastric tumor ecosystem and to translate these cell-resolved insights into a compact, generalizable prognostic model. Across six datasets (116 samples) spanning PT through AGC, we observed a reproducible reallocation of niche occupancy: immune contraction with epithelial expansion in premalignant stages, followed by concomitant enrichment of aberrant epithelium and myeloid lineages in EGC/AGC. Within the epithelial compartment, we resolved five states and identified a late-pseudotime SLC16A3⁺ program (MCT4/lactate export) that expands with progression and is transcriptionally coupled to ECM remodeling, pro-angiogenic signaling, and chemokine outputs.

Crucially, our in situ histological validation grounds this transcriptomic prediction in physical tumor architecture. We demonstrated that epithelial SLC16A3 positivity is markedly enriched in Lauren intestinal and mixed subtypes compared to diffused tumors. This spatial pattern implies that the dense, cohesive glandular nests typical of intestinal/mixed progression create severe, localized hypoxic bottlenecks. To survive these physical constraints, the aberrant epithelium is driven toward an MCT4-dependent metabolic climax, actively extruding lactate and chemokines into the immediate stroma.

This localized metabolic rewiring may contribute to immune-evasion-associated remodeling. Decomposition of the myeloid compartment indicated that this “myeloid” signal is macrophage-centric rather than neutrophil- or monocyte-driven, with preferential expansion of SPP1⁺ and APOE⁺ macrophages—states classically linked to immunoregulation, tissue remodeling, and angiogenesis. Together, these observations provide a mechanistic, cellular, and molecular basis for how physical tumor architecture shapes an immune-excluded, myeloid-rich microenvironment in late disease.

Mechanistically, our data support an association between epithelial SLC16A3/MCT4-linked lactate-export programs and a macrophage-rich immunoregulatory niche. The SLC16A3⁺ state was enriched for CXCL chemokines and growth-factor axes (VEGF/FGF/EGF) and showed ligand–receptor connectivity to macrophage receptors (CXCR), consistent with potential recruitment and retention of TAMs and their angiogenic/ECM functions. In contrast, neutrophil-directed axes (e.g., CXCL8-CXCR1/2) were not the dominant signal in this epithelial program. Functionally, macrophage states showed higher immunoregulatory/anti-inflammatory scores that rose from SIM to EGC/AGC, whereas pro-inflammatory signatures were comparatively attenuated. These layers suggest a putative link between epithelial metabolic remodeling and macrophage-rich immune exclusion; however, they should be interpreted as associative rather than directly causal (15,54,55,60-63).

An important alternative explanation is that hypoxia may independently contribute to the observed SLC16A3-associated lactate-myeloid phenotype. Hypoxia can induce glycolytic reprogramming and SLC16A3/MCT4 expression, while also promoting angiogenesis, ECM remodeling, and myeloid recruitment. Therefore, the association between SLC16A3-high epithelial states and macrophage-rich immune-exclusion features may partly reflect a hypoxia-adapted tumor niche rather than lactate export alone. Although we controlled for hypoxia and glycolysis scores during SLC16A3-associated gene selection, our transcriptomic data cannot fully disentangle the causal hierarchy among hypoxia, lactate export, extracellular acidosis, and macrophage conditioning. Future spatial and functional studies incorporating hypoxia markers, lactate measurements, and MCT4 perturbation will be required to clarify whether lactate export directly contributes to macrophage immunoregulatory remodeling in GC.

We distilled these single-cell insights into a 19-gene bulk prognostic model anchored on SLC16A3 co-expression and macrophage immunoregulatory modules, with candidate predictors prioritized by bootstrap stability and stage-aware differential correlation and finalized through a penalized Cox modeling pipeline. The resulting risk score stratified patients within the TCGA-STAD development framework (KM P=8.86×10⁻⁸; 1/3/5-year AUCs 0.695/0.744/0.757), validated in GEO (GSE84437) (KM P=6.68×10⁻³), and remained independent of clinicopathologic variables. These results suggest that a microenvironment-aware, mechanism-grounded prognostic model can retain predictive power across platforms and cohorts while preserving biological interpretability, thereby linking epithelial metabolism to macrophage-dominated immune exclusion.

Risk-stratified genomic and immunologic features further contextualize the model. Low-risk tumors exhibited higher non-synonymous mutation counts/TMB, whereas High-risk tumors were genomically quieter yet displayed higher TIDE, with the most pronounced difference in T-cell exclusion and positive correlations between the risk score and TIDE components. Thus, adverse outcome in the high-risk group appears not to be driven by hyper-mutation but by immune-exclusion biology, coherent with the SLC16A3⁺–macrophage axis and with the ECM/angiogenic programs observed in late epithelium and SPP1⁺/APOE⁺ TAMs.

Our in silico pharmacogenomic screen reveals an asymmetric sensitivity landscape across risk strata. Most ERBB/EGFR inhibitors (lapatinib, gefitinib, erlotinib, osimertinib, sapitinib) and the mTORC1/2 inhibitor OSI-027 exhibit lower median AUC in the Low-risk group, indicating greater predicted sensitivity in this stratum. By contrast, only a limited subset of agents shows high-risk-selective activity; among them, AZD8055 (mTORC1/2) displays a consistent lower AUC in high-risk, nominating mTOR dependence as a potential vulnerability of the immune-excluded state. Agents such as Acetalax show higher AUC in high-risk (Δ>0), consistent with relative resistance in high-risk and preferential activity in low-risk.

Mechanistically, these patterns suggest that the low-risk transcriptomic program—which is relatively less myeloid/ECM-dominated—retains broader RTK-MAPK/mTOR signaling liabilities that are responsive to ERBB and certain mTOR inhibitors and extends to IGF-1R (BMS-754807, GSK19045figure29A), ERK (ulixertinib), DNA-PK (NU7441), WEE1 (MK-1775), and other targets. In contrast, the High-risk phenotype, characterized by lactate/ECM/myeloid remodeling and T-cell exclusion, appears selectively susceptible to mTOR blockade with AZD8055, but comparatively less to ERBB inhibitors and OSI-027. This divergence within the mTOR class (AZD8055 vs. OSI-027) may reflect differences in target engagement, network rewiring, or off-target profiles, and underscores the need to treat drug names—not just pathways—as testable, drug-specific hypotheses.

From a translational standpoint, these expression-based predictions should be interpreted as hypothesis-generating rather than direct evidence of therapeutic sensitivity. The predicted differences in ERBB/EGFR, mTOR, and other drug-response patterns nominate candidate vulnerabilities for future validation, but treatment-linked clinical cohorts, patient-derived organoids, and functional perturbation assays will be required before these findings can inform therapeutic decision-making.

Strengths of this work include (I) a stage-resolved, multi-dataset single-cell reference with batch-aware integration; (II) a trajectory-based dissection of epithelial states linking SLC16A3⁺ epithelium to matrix, angiogenesis, and chemokines; (III) myeloid decomposition that separates macrophages from neutrophils/monocytes and maps functional programs across stages; and (IV) translation into a compact bulk signature with independent validation and mechanistic interpretability, followed by risk-stratified genomics, immune-evasion metrics, and drug-response hypotheses.

We also acknowledge limitations. First, our analyses are retrospective and integrate public cohorts with heterogeneous sampling and clinical annotation. Although we applied predefined dataset-selection criteria, stringent quality control, batch-aware integration, and external validation, several potential sampling biases should be acknowledged. The included datasets were generated by different institutions using different sequencing platforms, tissue-processing protocols, and patient-enrollment strategies. Disease-stage groups were not uniformly balanced, and samples across the premalignant-to-malignant continuum were not longitudinally collected from the same individuals. In addition, public scRNA-seq datasets may be affected by cell-recovery and dissociation biases, which can preferentially retain or deplete specific epithelial, immune, or stromal populations. Therefore, prospective cohorts with standardized sampling, treatment annotation, and matched multi-omics or spatial validation will be required to further confirm the generalizability of our findings. Second, scRNA-seq lacks spatial context; although the ligand–receptor and program analyses are consistent with macrophage recruitment and immune exclusion, spatial transcriptomics/multiplex imaging are required to confirm co-localization of SLC16A3/CA9 hotspots with CD68/CD163 macrophages and to quantify neighborhood enrichment relative to neutrophils. Third, drug-response results derive from expression-based predictors and should be considered preliminary; functional perturbation in organoids and co-culture models (e.g., MCT4 inhibition or acidosis buffering combined with ERBB/mTOR blockade) will help determine causality. Finally, although the 19-gene signature is mechanistically anchored, further assay development (e.g., RT-qPCR/NanoString panels) and cut-point optimization are needed for clinical deployment.

Future directions include: (I) spatial multi-omics to map epithelial-myeloid-ECM neighborhoods and lactate/pH gradients; (II) metabolite-aware profiling (lactate flux, pH, MCT activity) and causal perturbations in patient-derived models to validate the epithelial-to-macrophage axis; (III) longitudinal sampling across premalignant to malignant stages to observe state transitions in vivo; and (IV) prospective studies embedding the 19-gene assay to test risk-adapted therapeutic strategies and to evaluate its value alongside MSI/TMB and standard clinicopathologic factors.

In sum, our data support an epithelial lactate-export–associated, macrophage-dominated model of the gastric TME. The single-cell atlas explains how a late SLC16A3⁺ epithelial state can wire chemokine, angiogenic and ECM programs that favor T-cell exclusion, and the derived 19-gene signature captures this biology in bulk to stratify prognosis, parse immune evasion, and nominate therapeutic hypotheses. This mechanism-to-model framework may aid risk stratification and combination-therapy design targeting metabolism-myeloid-ECM axes in GC.

Conclusions

We present a progression-resolved single-cell atlas of GC that links an epithelial SLC16A3-driven lactate-export program to expansion of SPP1⁺/APOE⁺ macrophage states and immune exclusion. Crucially, our in situ immunohistochemical validation anchors this metabolic climax to specific physical tumor architectures, demonstrating robust SLC16A3 enrichment specifically within the dense cohesive nests of Lauren intestinal and mixed subtypes. From this biology we developed a 19-gene prognostic model whose derived risk score stratifies survival independently of clinical covariates and validates across cohorts. Risk groups diverged in genomic burden and predicted candidate pharmacologic vulnerabilities requiring validation: low-risk tumors showed higher mutation counts/TMB and broader computationally predicted drug-response signals for EGFR/ERBB inhibitors and OSI-027, whereas high-risk tumors were genomically quieter yet TIDE-high, with AZD8055 emerging as a candidate predicted signal requiring validation. These findings provide a mechanistic bridge between epithelial metabolism and microenvironmental remodeling, yielding testable hypotheses for patient selection and for combination strategies targeting lactate transport and macrophage reprogramming; prospective, treatment-linked validation and orthogonal assays (IHC/spatial/proteomics) are warranted for clinical translation.

Acknowledgments

We are grateful to the contributors who contributed data to the public databases used in this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0633/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0633/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0633/prf

Funding: This work was supported by the National Natural Science Foundation of China (Nos. 92474106 and 82403738), the Noncommunicable Chronic Diseases—National Science and Technology Major Project (No. 2024ZD0520600), and the Medical Science and Technology Project of Zhejiang Province (No. WKJ-ZJ-2310).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0633/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Institutional Review Board of The Second Affiliated Hospital, Zhejiang University School of Medicine (No. 2023.1107). Written informed consent was obtained from all patients prior to enrollment and sample collection.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Aslani A, Soheili A, Mousavi SE, et al. Incidence trends of gastric cancer in the United States over 2000-2020: A population-based analysis. PLoS One 2024;19:e0310040. [Crossref] [PubMed]
3. Correa P. Human gastric carcinogenesis: a multistep and multifactorial process--First American Cancer Society Award Lecture on Cancer Epidemiology and Prevention. Cancer Res 1992;52:6735-40.
4. Correa P, Piazuelo MB. The gastric precancerous cascade. J Dig Dis 2012;13:2-9. [Crossref] [PubMed]
5. Li J, Jiang Y, Ma M, et al. Epithelial cell diversity and immune remodeling in bladder cancer progression: insights from single-cell transcriptomics. J Transl Med 2025;23:135. [Crossref] [PubMed]
6. Ma M, Sun J, Liu Z, et al. The Immune Microenvironment in Gastric Cancer: Prognostic Prediction. Front Oncol 2022;12:836389. [Crossref] [PubMed]
7. Sahai E, Astsaturov I, Cukierman E, et al. A framework for advancing our understanding of cancer-associated fibroblasts. Nat Rev Cancer 2020;20:174-86. [Crossref] [PubMed]
8. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
9. Estrella V, Chen T, Lloyd M, et al. Acidity generated by the tumor microenvironment drives local invasion. Cancer Res 2013;73:1524-35. [Crossref] [PubMed]
10. Binnewies M, Roberts EW, Kersten K, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med 2018;24:541-50. [Crossref] [PubMed]
11. Zhou Z, Xia G, Xiang Z, et al. A C-X-C Chemokine Receptor Type 2-Dominated Cross-talk between Tumor Cells and Macrophages Drives Gastric Cancer Metastasis. Clin Cancer Res 2019;25:3317-28. [Crossref] [PubMed]
12. Sun K, Xu R, Ma F, et al. scRNA-seq of gastric tumor shows complex intercellular interaction with an alternative T cell exhaustion trajectory. Nat Commun 2022;13:4943. [Crossref] [PubMed]
13. Kim J, Park C, Kim KH, et al. Single-cell analysis of gastric pre-cancerous and cancer lesions reveals cell lineage diversity and intratumoral heterogeneity. NPJ Precis Oncol 2022;6:9. [Crossref] [PubMed]
14. Huang KK, Ma H, Chong RHH, et al. Spatiotemporal genomic profiling of intestinal metaplasia reveals clonal dynamics of gastric cancer progression. Cancer Cell 2023;41:2019-2037.e8. [Crossref] [PubMed]
15. Manoharan I, Prasad PD, Thangaraju M, et al. Lactate-Dependent Regulation of Immune Responses by Dendritic Cells and Macrophages. Front Immunol 2021;12:691134. [Crossref] [PubMed]
16. Wei C, Ma Y, Wang M, et al. Tumor-associated macrophage clusters linked to immunotherapy in a pan-cancer census. NPJ Precis Oncol 2024;8:176. [Crossref] [PubMed]
17. Kim SJ, Kim M, Kim KH, et al. APOE-high myeloid cells are uniquely associated with metastatic intrathoracic lymph nodes obtained by EBUS-TBNA in primary lung cancer. NPJ Precis Oncol 2025;9:292. [Crossref] [PubMed]
18. Obradovic A, Chowdhury N, Haake SM, et al. Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages. Cell 2021;184:2988-3005.e16. [Crossref] [PubMed]
19. Cheong JH, Yang HK, Kim H, et al. Predictive test for chemotherapy response in resectable gastric cancer: a multi-cohort, retrospective analysis. Lancet Oncol 2018;19:629-38. [Crossref] [PubMed]
20. Yoon SJ, Park J, Shin Y, et al. Deconvolution of diffuse gastric cancer and the suppression of CD34 on the BALB/c nude mice model. BMC Cancer 2020;20:314. [Crossref] [PubMed]
21. Zhang P, Yang M, Zhang Y, et al. Dissecting the Single-Cell Transcriptome Network Underlying Gastric Premalignant Lesions and Early Gastric Cancer. Cell Rep 2019;27:1934-1947.e5. [Crossref] [PubMed]
22. Jiang H, Yu D, Yang P, et al. Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing. Clin Transl Med 2022;12:e730. [Crossref] [PubMed]
23. Jeong HY, Ham IH, Lee SH, et al. Spatially Distinct Reprogramming of the Tumor Microenvironment Based On Tumor Invasion in Diffuse-Type Gastric Cancers. Clin Cancer Res 2021;27:6529-42. [Crossref] [PubMed]
24. Kang B, Camps J, Fan B, et al. Parallel single-cell and bulk transcriptome analyses reveal key features of the gastric tumor microenvironment. Genome Biol 2022;23:265. [Crossref] [PubMed]
25. Sathe A, Grimes SM, Lau BT, et al. Single-Cell Genomic Characterization Reveals the Cellular Reprogramming of the Gastric Tumor Microenvironment. Clin Cancer Res 2020;26:2640-53. [Crossref] [PubMed]
26. Gayoso A, Lopez R, Xing G, et al. A Python library for probabilistic analysis of single-cell omics data. Nat Biotechnol 2022;40:163-6. [Crossref] [PubMed]
27. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol 2018;19:15. [Crossref] [PubMed]
28. Street K, Risso D, Fletcher RB, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics 2018;19:477. [Crossref] [PubMed]
29. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
30. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
31. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 2014;9:e107468. [Crossref] [PubMed]
32. Haber AL, Biton M, Rogel N, et al. A single-cell survey of the small intestinal epithelium. Nature 2017;551:333-9. [Crossref] [PubMed]
33. Adekola K, Rosen ST, Shanmugam M. Glucose transporters in cancer metabolism. Curr Opin Oncol 2012;24:650-4. [Crossref] [PubMed]
34. Wei X, Liu X, Tan C, et al. Expression and Function of Zinc-α2-Glycoprotein. Neurosci Bull 2019;35:540-50. [Crossref] [PubMed]
35. Austin RC. The unfolded protein response in health and disease. Antioxid Redox Signal 2009;11:2279-87. [Crossref] [PubMed]
36. Walter P, Ron D. The unfolded protein response: from stress pathway to homeostatic regulation. Science 2011;334:1081-6. [Crossref] [PubMed]
37. Citi S, Fromm M, Furuse M, et al. A short guide to the tight junction. J Cell Sci 2024;137:jcs261776. [Crossref] [PubMed]
38. Anderson JM, Van Itallie CM. Physiology and function of the tight junction. Cold Spring Harb Perspect Biol 2009;1:a002584. [Crossref] [PubMed]
39. Hansson GC. Mucus and mucins in diseases of the intestinal and respiratory tracts. J Intern Med 2019;285:479-90. [Crossref] [PubMed]
40. Keum YS. Regulation of Nrf2-Mediated Phase II Detoxification and Anti-oxidant Genes. Biomol Ther (Seoul) 2012;20:144-51. [Crossref] [PubMed]
41. Dasari S, Ganjayi MS, Yellanurkonda P, et al. Role of glutathione S-transferases in detoxification of a polycyclic aromatic hydrocarbon, methylcholanthrene. Chem Biol Interact 2018;294:81-90. [Crossref] [PubMed]
42. Grasberger H, Gao J, Nagao-Kitamoto H, et al. Increased Expression of DUOX2 Is an Epithelial Response to Mucosal Dysbiosis Required for Immune Homeostasis in Mouse Intestine. Gastroenterology 2015;149:1849-59. [Crossref] [PubMed]
43. Poirier Y, Antonenkov VD, Glumoff T, et al. Peroxisomal beta-oxidation--a metabolic pathway with multiple functions. Biochim Biophys Acta 2006;1763:1413-26. [Crossref] [PubMed]
44. Grivennikov S, Karin E, Terzic J, et al. IL-6 and Stat3 are required for survival of intestinal epithelial cells and development of colitis-associated cancer. Cancer Cell 2009;15:103-13. [Crossref] [PubMed]
45. Xiong X, Liao X, Qiu S, et al. CXCL8 in Tumor Biology and Its Implications for Clinical Translation. Front Mol Biosci 2022;9:723846. [Crossref] [PubMed]
46. Barker N, van Es JH, Kuipers J, et al. Identification of stem cells in small intestine and colon by marker gene Lgr5. Nature 2007;449:1003-7. [Crossref] [PubMed]
47. van der Flier LG, Haegebarth A, Stange DE, et al. OLFM4 is a robust marker for stem cells in human intestine and marks a subset of colorectal cancer cells. Gastroenterology 2009;137:15-7. [Crossref] [PubMed]
48. Yan KS, Kuo CJ. Ascl2 reinforces intestinal stem cell identity. Cell Stem Cell 2015;16:105-6. [Crossref] [PubMed]
49. Saxton RA, Sabatini DM. mTOR Signaling in Growth, Metabolism, and Disease. Cell 2017;168:960-76. [Crossref] [PubMed]
50. Vences-Catalán F, Levy S. Immune Targeting of Tetraspanins Involved in Cell Invasion and Metastasis. Front Immunol 2018;9:1277. [Crossref] [PubMed]
51. Yarden Y, Pines G. The ERBB network: at last, cancer therapy meets systems biology. Nat Rev Cancer 2012;12:553-63. [Crossref] [PubMed]
52. Lamouille S, Xu J, Derynck R. Molecular mechanisms of epithelial-mesenchymal transition. Nat Rev Mol Cell Biol 2014;15:178-96. [Crossref] [PubMed]
53. Faubert B, Li KY, Cai L, et al. Lactate Metabolism in Human Lung Tumors. Cell 2017;171:358-371.e9. [Crossref] [PubMed]
54. Zhang D, Tang Z, Huang H, et al. Metabolic regulation of gene expression by histone lactylation. Nature 2019;574:575-80. [Crossref] [PubMed]
55. Semenza GL. HIF-1 mediates metabolic responses to intratumoral hypoxia and oncogenic mutations. J Clin Invest 2013;123:3664-71. [Crossref] [PubMed]
56. Moon JW, Kong SK, Kim BS, et al. IFNγ induces PD-L1 overexpression by JAK2/STAT1/IRF-1 signaling in EBV-positive gastric carcinoma. Sci Rep 2017;7:17810. [Crossref] [PubMed]
57. Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity's roles in cancer suppression and promotion. Science 2011;331:1565-70. [Crossref] [PubMed]
58. Strelez C, Battaglin F, Perez R, et al. GABAergic signaling contributes to tumor cell invasion and poor overall survival in colorectal cancer. Oncogene 2025;44:3924-38. [Crossref] [PubMed]
59. Wang Y, Zhang Y, Li W, et al. GABA(A) receptor π forms channels that stimulate ERK through a G-protein-dependent pathway. Mol Cell 2025;85:166-176.e5. [Crossref] [PubMed]
60. Becker HM. Carbonic anhydrase IX and acid transport in cancer. Br J Cancer 2020;122:157-67. [Crossref] [PubMed]
61. Contreras-Baeza Y, Sandoval PY, Alarcón R, et al. Monocarboxylate transporter 4 (MCT4) is a high affinity transporter capable of exporting lactate in high-lactate microenvironments. J Biol Chem 2019;294:20135-47. [Crossref] [PubMed]
62. Choi SH, Kim MY, Yoon YS, et al. Hypoxia-induced RelA/p65 derepresses SLC16A3 (MCT4) by downregulating ZBTB7A. Biochim Biophys Acta Gene Regul Mech 2019;1862:771-85. [Crossref] [PubMed]
63. Lundø K, Trauelsen M, Pedersen SF, et al. Why Warburg Works: Lactate Controls Immune Evasion through GPR81. Cell Metab 2020;31:666-8. [Crossref] [PubMed]