Intratumoral microbiome is associated with gastric cancer prognosis and therapy efficacy

ABSTRACT

The role of the intratumoral microbiome in gastric cancer (GC) has not been comprehensively assessed. Here, we explored the relationship between the microbial community and GC prognosis and therapy efficacy. Several cancer-associated microbial characteristics were identified, including increased α-diversity, differential β-diversity, and decreased Helicobacter pylori abundance. After adjusting for clinical features, prognostic analysis revealed 2 phyla, 14 genera, and 5 species associated with the overall survival of patients with GC. Additionally, 2 phyla, 14 genera, and 6 species were associated with adjuvant chemotherapy (ACT) efficacy in patients with stage II – III GC. Furthermore, we classified GC microbiome structures into three microbial subtypes (MS1, MS2 and MS3) with distinguishing features. The MS1 subtype exhibited high immune activity and enrichment of microbiota related to immunotherapy and butyric acid-producing, as well as potential benefits in immunotherapy. MS2 featured the highest α-diversity and activation of the TFF pathway, MS3 was characterized by epithelial-mesenchymal transition and was associated with poor prognosis and reduced ACT efficacy. Collectively, the results of this study provide valuable insights into the microbial characteristics associated with GC prognosis and therapy efficacy.

GRAPHICAL ABSTRACT

KEYWORDS:

• Prognosis
• adjuvant chemotherapy
• immunotherapy
• microbial subtype

Introduction

As detection technologies have advanced and our understanding of the tumor microenvironment (TME) has deepened, increasing evidence has confirmed the presence of microbiota within, and complex interactions with, tumors.^1,2 Intratumoral microbiota may affect cancer development and progression through myriad mechanisms, such as reshaping the TME, promoting tumor metastasis, and influencing therapy efficacy.³ However, the precise role that intratumoral microbiota have in cancer pathogenesis, remains unclear. Hence, elucidating the complex interactions between the microbiota, TME, and cancer cells can provide valuable insights into potential cancer therapies and aid the evaluation of therapy efficacy for individual patients.

Gastric cancer (GC) is a complex disease with multiple risk factors. In particular, the intratumoral microbiota has emerged as a potential risk factor. Historically, the human stomach was thought to serve as the exclusive habitat for Helicobacter pylori due to its acidity and other antimicrobial factors, making it unsuitable for other microorganisms. However, this view has recently been challenged as non-Helicobacter pylori-microbiota were found to impact GC progression by altering immune and metabolic homeostasis.^4–8 Moreover, several studies have extensively characterized the microbiota profiles of the stomach and GC, and those associated with GC development.^9–12 However, there remains a lack of research into the relationship between intratumoral microbiota and the clinical features, prognosis, and chemotherapy and immunotherapy efficacy in GC. Hence, the clinical translation of gut microbiota in GC research has been limited.

The rapid advancement of omics technologies has led to the molecular and genomic characterization of GC, representing a significant area of research. Indeed, new GC subtypes have been established based on genomic, transcriptomic, and proteomic analyses, improving the stratification of patient prognosis and advancing the development of targeted therapies.^13–15 However, our understanding of microbial subtypes (Msubtypes) remains limited. Current research has focused on individual pathogens or is constrained by small sample sizes and a lack of clinical outcome information. Hence, studies seeking to characterize the Msubtypes associated with GC and elucidate their clinical significance and molecular features, are greatly needed.

In this study, we conducted a comprehensive analysis combining microbiome (16S rRNA gene sequencing) data from 251 GC samples with paired noncancerous adjacent tissues (NATs), transcriptome data from 94 GC samples, and detailed clinical characteristics. Our study aimed to elucidate the different microbial characteristics among clinical subgroups and identify microbiota with independent prognostic significance for overall survival (OS) and adjuvant chemotherapy (ACT). Additionally, we performed microbiome-centered cluster analysis to identify subtypes with distinct molecular characteristics associated with GC prognosis and chemotherapy and immunotherapy efficacy. Taken together, our findings reveal the association between intratumoral microbiota and clinical outcomes in GC, providing valuable insights for precision medicine.

Results

GC intratumoral microbiome landscape

Primary GC samples and paired NATs were collected from 251 patients with GC. The median age of the patients was 62 years (range: 36–84 years). Of the GC samples, 11.6% were Epstein-Barr virus (EBV)-positive (n = 29), 30.7% were mismatch repair deficient (dMMR; n = 77), and 57% were proficient mismatch repair (pMMR; n = 150). Moreover, 12.4% of the GC samples were stage I (n = 31), 36.7% were stage II (n = 92), 50.2% were stage III (n = 126), and 0.8% were stage IV (n = 2). Additionally, 66.5% (n = 167) of the patients received ACT after surgery. All clinical information is summarized in Table S1 and presented in Figure 1a. Information related to quality control is summarized in Table S2. The rarefaction curve indicated that the microecological characteristics of GC and NAT reached saturation with increased sequencing depth (Figure 1b).

Figure 1. Intratumoral microbiome landscape of gastric cancer (GC). (a) Summary of the clinical characteristics of 251 patients with GC. (b) Rarefaction curve of GC and noncancerous adjacent tissues (NAT). (c) α-diversity difference between GC and NAT. (d) Principal coordinate analysis (PCoA) of the microbiome data for GC and NAT. (e) Differences in microbiota abundances between GC and NAT. (f) α-Diversity differences and PCoA among dMMR, EBV and pMMR subtypes. (g) Microbiota abundance differences among dMMR, EBV and pMMR subtypes. (h) Correlations between genera (top 20) and EBV load. The size of the circle represents the magnitude of the Spearman correlation coefficient.

We evaluated the differences in the overall intratumoral microbial community structures between GC samples and NAT (251 pairs). Analyses of the GC and NAT genera profiles revealed a weak correlation (ρ = 0.08; Figure S1A). Moreover, GC samples exhibited increased α-diversity compared with those in NATs at both ASV and genus levels (Figure 1c). Furthermore, β-diversity analysis revealed a significant separation between GC and NAT samples (Figure 1d; PERMANOVA: p < 0.001). We also identified 16 phyla, 75 genera, and 34 species that were significantly enriched in GC compared with NAT (Adj p < 0.05; Figure 1e, Table S3–S5), including several well-established GC-associated bacteria such as Streptococcus, Prevotella, Fusobacterium, Parvimonas, and Lactobacillus.^16,17 In addition, the Helicobacter genus and Helicobacter pylori species were predominant in NAT with markedly decreased abundance in the GC samples (Figure S1B). These findings highlight the distinct microbial ecologies between GC and NAT.

Considering the lack of consensus on the relationship between microbiota and clinical features, we compared the intratumoral microbial community structures of different clinical subgroups based on age, sex, pathology, tumor grade, location, vascular invasion, TNM stage, Lauren subtype, and molecular subtype. The results revealed that only the molecular subtype exhibited significant differences in microbial community structures (Figure 1f and Figure S1C – S1D). Moreover, compared to pMMR subtype, both the EBV and dMMR subtypes showed a large number of different genera (Figure 1g). These findings suggest that dMMR and EBV subtypes have relatively specific microbial characteristics.

Currently, little is known regarding the association between EBV subtypes and microbiota. Hence, we conducted a transcriptome analysis of 94 tumor samples and aligned the data with the EBV genome to assess the viral load of each sample. As expected, the EBV subtypes exhibited higher EBV load. (Figure S1E). Further correlation analyses identified the top 20 genera positively correlated with the EBV load; Butyricicoccus was among these genera (Figure 1h). Interestingly, previous studies have reported a positive association between Butyricicoccus and EBV in the blood and plasma.¹⁸ Furthermore, we observed a significant increase in the abundance of Butyricicoccus in the EBV subtype (Figure S1F). These results suggest a potential regulatory relationship between EBV and other microbiota in the EBV subtype.

Exploring the microbial taxa associated with GC OS and ACT

Increasing evidence has shown that intratumoral bacteria are associated with the prognosis of cancer patients.^19–21 However, the relationship between GC prognosis and the microbial profile remains unclear. Hence, we performed univariate Cox analysis to identify prognostic microbiota associated with patient OS at the phylum, genus, and species levels. Results identified 2 phyla, 23 genera, and 9 species associated with GC OS (Figure S2A). Furthermore, multivariable Cox regression was used to adjust for multiple clinical features, including age, sex, pathology, grade, location, vascular invasion, TNM stage, Lauren subtype, and molecular subtype, and to identify independent prognostic microbiota. Ultimately, 2 phyla, 14 genera, and 5 species were identified as independent prognostic microbiota for the OS of GC (Figure 2a).

Figure 2. Taxa associated with overall survival (OS) and adjuvant chemotherapy (ACT) after adjusting for clinical features. (a) Phyla, genera, and species associated with OS, after adjusting for clinical features. (b) Phyla, genera, and species associated with ACT, after adjusting for clinical features.

Considering that 5-fluorouracil-based ACT is the standard treatment for postoperative stage II and III GC, there is still a lack of awareness regarding the prognostic factors that influence the outcomes in this particular group.²² Therefore, 2 phyla, 22 genera, and 12 species were identified to be associated with the OS of stage II – III GC receiving ACT (Figure S2B). After adjusting for clinical features, 2 phyla, 14 genera, and 6 species were considered the key microbiota influencing ACT efficacy (Figure 2b).

Identification of intratumoral microbial subtypes

The coexistence of microorganisms with anti-tumor and pro-tumor effects within the TME emphasizes the need to comprehensively evaluate the microbial characteristics rather than focusing on a single species. To gain a comprehensive understanding of the relationship between the intratumoral microbiome and GC outcomes, we identified three Msubtypes in the intratumoral microbial profiles (Figure S3). The MS1 subtype contained the most subtype-specific genera, including Parabacteroides, Ruminococcus, Butyricicoccus, and Lactobacillus. The MS2 subtype contained the most tumor-enriched genera, including Rhodoplanes, Nocardioides, Kaistobacter, and Pseudonocardia. Meanwhile, the MS3 subtype was characterized by Sphingobium, Delftia, Comamonas, and Stenotrophomonas. Delftia was associated with poor OS and poor ACT efficacy (Figure 3a). These three subtype categories were further supported by co-occurrence analysis of genera exhibiting strong positive correlations (Figure S4A).

Figure 3. Msubtypes and their correlation with clinical gastric cancer (GC) outcomes. (a) Heatmap of differentially abundant genera (Adj p < 0.05, Log2 FC > 2, and AUC > 0.7) among the Msubtypes. (b) Kaplan-Meier curve for overall survival (OS) among the Msubtypes. (c) α-Diversity differences and PCoA among the Msubtypes. (d) Co-abundance patterns among the genera of the Msubtypes.

In the joint analysis with clinical features, EBV subtype was enriched within the MS1 subtype, dMMR subtype within the MS2 subtype, and pMMR subtype within the MS3 subtype. Moreover, no significant correlation was observed between the Msubtype and other clinical features (Table S6). Furthermore, survival analysis revealed that patients assigned to the MS3 subtype had the worst OS (Figure 3b). Additionally, the MS2 subtype had the highest α-diversity, and the MS3 subtype had the lowest. Meanwhile, the β-diversity analysis revealed that the microbial community exhibited phylogenetic differences within each subtype (PERMANOVA: p = 0.001; Figure 3c and Figure S4B).

To assess the metabolic differences between Msubtypes, we constructed a virtual metagenome for each sample’s microbiome using PICRUSt2 software. Several putative metacommunity-related functional shifts were identified between Msubtypes. More specifically, MS1 was characterized by the potential enrichment of carbohydrate metabolism pathways, MS2 exhibited increased potential for oxidative phosphorylation, TCA cycle, and amino acid metabolism, and MS3 might be associated with bacterial secretion system and glutathione metabolism (Figure 3d). To further investigate the variation patterns among the Msubtypes, we characterized the co-abundance patterns of each subtype. The MS1 subtype was characterized by the dominance of Halomonas and the hub genus Blautia. In contrast, MS2 was predominated by Burkholderia and the hub genus Chelativorans. The MS3 subtype was characterized by Acinetobacter and the hub genus Propionibacterium (Figure 3e). In total, 225 (MS1), 155 (MS2), and 145 (MS3) associations were identified (Figure S4C). Interestingly, the Msubtypes shared few associations, indicating significant changes in the symbiotic networks among the subtypes.

Multiple cohort analysis reveals that MS3 is associated with poor OS, RFS and poor chemotherapy efficacy

To further explore the biological characteristics of the three subtypes, we evaluated the transcriptome data for 94 samples. Gene Set Enrichment Analysis (GSEA) indicated that the MS1 subtype exhibited enrichment of numerous immune-related pathways, including PD1 signaling, suggesting its potential suitability for immunotherapy. The MS2 subtype was characterized by the TFF pathway, which plays an important role in maintaining the integrity of the gastrointestinal tract and promoting tissue repair. Meanwhile, the MS3 subtype was characterized by digestion and absorption, WNT signaling, and epithelial-mesenchymal transition (EMT; Figure 4a). According to previous reports, Helicobacter pylori potentiates EMT in GC.^23,24 Interestingly, we found that the MS3 subtype had the highest Helicobacter pylori abundance (Figure 4b).

Figure 4. MS3 is associated with poor overall survival (OS), regression-free survival (RFS), and poor chemotherapy efficacy. (a) GSEA of enriched pathways at the transcriptome level among Msubtypes. (b) Differences in Helicobacter pylori abundances among Msubtypes. (c) Kaplan-Meier curves for OS among the Msubtypes based on patients from the five transcriptome cohorts. (d) Kaplan-Meier curves for OS among the Msubtypes based on patients with stage II III gastric cancer (GC) who received adjuvant chemotherapy (ACT), from the three transcriptome cohorts. (e) Differences in response to neoadjuvant chemotherapy (NAC) among the Msubtypes, based on the published PUCH cohort data.

Immune infiltration analysis indicated that the MS1 subtype was dominated by CD8⁺ T cells (Figure S5A, Table S7), which may explain the better prognosis for patients with the MS1 subtype. We further determined whether the relationship between our subtypes and prognosis was universal by applying the Msubtypes to six published GC transcriptomic datasets (see Supplementary methods and Figure S5B). Survival analysis of multiple cohorts revealed that patients with the MS3 subtype had the worst OS and recurrence-free survival (RFS; Figure 4c and S5C). For patients with stage II – III GC who received ACT, the MS3 subtype was also associated with the worst OS and RFS (Figure 4d and S5D). Moreover, the MS3 subtype exhibited resistance to neoadjuvant chemotherapy (NAC) in a cohort of patients at our center (Figure 4e). These results demonstrate the clinical relevance of the Msubtypes.

MS1 is associated with favorable efficacy of immunotherapy

The MS1 subtype is characterized by increased immune-related pathways and PD1 signaling, indicating its potential suitability for anti-PD1/PDL1 therapy. To further validate this result, we performed an integrative analysis of a curated list of 166 genes, including immune checkpoint genes and other known and emerging viable immunomodulatory targets (Table S8). The results showed that the expression of genes from all six immune-related categories was significantly upregulated in the MS1 subtype, including CD274 and IDO1 (Figure 5a and Figure S6A).

Figure 5. MS1 is associated with favorable immunotherapy efficacy. (a) Differences in the expression of immune genes among the Msubtypes. (b) Differences in immunotherapy-associated mRNA signatures among the Msubtypes. (c) Differences in immunotherapy-associated microbiota signatures among Msubtypes. (d) Differences in CD8⁺ T cell-associated microbiota signatures among the Msubtypes. (e) Heatmap of butyric acid-producing genera among the Msubtypes. (f) Submap analysis of the three Msubtypes between our cohort and the response (R) and non-response (NR) subgroups of two gastric cancer (GC) cohorts who received immunotherapy.

Considering that transcriptome and microbiome studies have aided the prediction of immunotherapy efficacy in various cancer types, we collected six mRNA signatures and nine microbiota signatures associated with immunotherapy response and CD8⁺ T cell infiltration (Table S9). The scores for these signatures were significantly higher in the MS1 subtype than the other subtypes in our cohort (Figure 5b-d). Further analysis showed that the MS1 subtype could effectively distinguish potential immunotherapy-responsive patients from non-responsive patients within the dMMR, EBV, and pMMR subtypes based on the microbiota signatures (Figure S6B). Intriguingly, several butyric acid-producing genera were significantly enriched in the MS1 subtype (Figure 5e); notably, butyric acid has been demonstrated to activate cytotoxic CD8+ T cells and enhance the effectiveness of anti-PD1 therapy.^25,26

To further analyze the relationship between MS1 and immunotherapy, we performed a submap analysis between the three subtypes in our cohorts and the response (R) and non-response (NR) subgroups of the two GC cohorts receiving immunotherapy. The MS1 subtype was significantly associated with the responsive subgroups (Figure 5f), suggesting that the MS1 subtype may benefit from immunotherapy.

The regulation and shaping of tumor immunity are influenced by microbial microenvironments and the metabolites produced by microbes^27,28 Currently, studies are establishing links between microbes and immune function.²⁹ The use of omics has accelerated the analyses of these relationships. Thus, to elucidate the association between microbes and immunity in these three Msubtypes, we applied the O2PLS model to perform a joint analysis of immune-related genes and microbial genera. The 15 most relevant immune genes and genera were identified for each Msubtype, which were used to construct networks. In the MS1 subtype, IDO1 and CD274 had the most associations with the microbiome, suggesting that they may be key targets for regulating the microbial microenvironment in this subtype. In the MS2 subtype, four genera exhibited extensive associations with immune genes, indicating that these genera might affect immunity in the MS2 subtype. However, the MS3 subtype lacked interactions between the immune system and microbiome (Figure S6C).

Discussion

In this study, we conducted intratumoral microbiome profiling of patients with GC. Our study sheds light on the relationship between the intratumoral microbiome and clinical features and outcomes of GC. We also identified three microbial subtypes with distinct microbial and genetic characteristics, offering potential guidance for personalized patient therapy. Moreover, the results of this study provide a valuable resource for further research into the clinical applications of the intratumoral microbiota.

Although Helicobacter pylori has been implicated in the malignant transformation, invasion, metastasis, and immune suppression of GC,³⁰ we found that the abundance of Helicobacter pylori was significantly reduced in GC compared to NAT and was not directly correlated with clinical outcomes. Meanwhile, the MS3 subtype, associated with the worst survival, had the highest Helicobacter pylori abundance. This suggests that Helicobacter pylori infection may play a complex role in GC initiation and progression.

While various omics studies have established a significant correlation between gut microbiota and cancer prognosis and treatment response,^31–33 research on intratumoral microbiota is limited. Our study addressed this issue and identified 2 phyla, 21 genera, and 8 species associated with OS or ACT efficacy in GC. These findings suggest that the microbiota composition can serve as a potential prognostic marker and may influence the efficacy of chemotherapy in GC. Interestingly, no beneficial microbiota associated with GC prognosis were found. Hence, from a translational perspective, the use of antibiotics rather than probiotics may improve GC prognosis. Indeed, antibiotics can inhibit the growth of GC cells and enhance the efficacy of chemotherapy for advanced GC.^34,35 Our findings suggest the possibility that selective bacterial eradication or manipulation might improve prognosis. Although these findings require further validation, they provide insights into the potential clinical translation of the intratumoral microbiome in GC.

We employed a microbiota-based classification method to identify three GC subtypes, each with distinct molecular features that connect clinical and molecular characteristics. Specifically, we found that the MS1 subtype exhibited high immune activity and upregulated expression of multiple immune checkpoints, such as CD274 and IDO1. This suggests that immunotherapy may have potential benefits for patients with this subtype. This was supported further by the immune-related microbial features. Although these features were derived from the studies of the gut microbiome in patients undergoing immunotherapy, including two other GC-based studies,^36,37 suggested that the intratumoral microbiome may originate from the gut microbiome.^38,39 Interestingly, we observed the enrichment of butyric acid-producing genera in the MS1 subtype, indicating the potential enrichment of butyric acid. Butyric acid, an important short-chain fatty acid, plays a crucial role in immune regulation and gastrointestinal disorders.^40,41 Recent studies have also indicated its positive impact on promoting CD8+ T cells activation and enhancing the efficacy of immunotherapy.^25,26 The MS2 subtype featured the highest α-diversity and enrichment of the TFF pathway. The TFF pathway is associated with preserved gastric function, which may account for the relatively favorable prognosis of patients with this subtype. In contrast, the MS3 subtype was characterized by EMT features at the transcriptome level, which have been implicated in various aspects of tumor development, including tumor invasion, metastasis, cancer stemness, and therapy resistance.^42,43 Additionally, microbiome analysis revealed the increased potential for glutathione metabolism that have been linked to platinum resistance.^44,45 These features may contribute to the poor clinical prognosis and response to ACT observed in patients with the MS3 subtype. Our results emphasize the potential of the Msubtype as a potential application with a future opportunity for identification of correlated microbiotas to facilitate clinical decision-making. In addition, our findings suggest the possibility that selective bacterial eradication or manipulation might improve therapy efficacy. For instance, the oral administration of Bifidobacterium has the potential to enhance immunotherapy in melanoma.⁴⁶ In addition, transplanting the gut microbiota of chemotherapy-responsive mice can reverse chemotherapy resistance.⁴⁷

This study used a large independent cohort of microbiome samples as well as detailed clinical and therapy data to identify and elucidate the molecular characteristics and clinical significance of microbial subtypes. These findings were validated in multiple independent cohorts. Despite the strengths of this study, there were certain limitations. First, the sequencing resolution of the 16s rRNA was limited and did not allow for precise microbiota identification. Second, we observed a significant enrichment of butyric acid-producing genera in the MS1 subtype, suggesting a potential increase in butyric acid. However, this finding needs to be further validated by metabolite detection. Third, owing to the lack of research on intratumoral microbiota in GC and small sample sizes, as well as the lack of prognostic information, we were unable to identify an appropriate cohort of microbiome samples for external validation. Further research on the microbial function and interaction, together with improved sequencing techniques, are required to validate and expand upon these findings.

In conclusion, our study provides a comprehensive microbial perspective that enhances our understanding of GC. First, we identified the microbiota associated with GC prognosis and ACT efficacy. Notably, we discovered that the intratumoral microbiome pattern, referred to the Msubtype, is associated with OS as well as chemotherapy and immunotherapy efficacy in GC. Importantly, our study provides an integrated and extensive microbial resource for the further exploration of the pathogenesis of GC.

Material and methods

Clinical specimens of GC

GC samples were collected from the Peking University Cancer Hospital in China. None of the enrolled patients had received therapy before surgery, and all were followed up for five years after resection. After radical gastrectomy, resected tumor samples and NATs (>5 cm from the tumor edge) were obtained and snap-frozen in liquid nitrogen in less than 30 min after resection. Then the tissues were stored at − 80°C. To ensure the quality of tissues, routine histological evaluation was performed for each sample. All of the samples were obtained, stored, and transported in a sterile environment. Clinical information was collected for the patients, including age, sex, pathology, grade, location, vascular invasion, TNM stage, Lauren subtype, molecular subtype, and ACT information (Table S1). Patients who received ACT were defined as those who received at least one cycle of 5-fluorouracil-based ACT.⁴⁸ OS was defined as the time interval between surgery and death. This study was approved by the Research Ethics Committee of Peking University Cancer Hospital (2019KT111), and written informed consent was obtained from all patients.

dMMR status determination and EBV in situ hybridization

Immunohistochemistry (IHC) for MLH1, PMS2, MSH2, and MSH6 was performed as previously described to determine the deficient mismatch repair (dMMR) status.⁴⁹ The monoclonal antibodies used specific for MLH1 (Clone ES05, DAKO), MSH2 (Clone FE11, DAKO), PMS2 (Clone EP51, DAKO), and MSH6 (EP49, DAKO). Tumors with a negative stain for any of these markers were considered dMMR, and tumors with positive staining for all four were defined as proficient mismatch repair (pMMR).

Chromogenic in situ hybridization with EBV-encoded small RNA (EBER) was performed to detect EBV infection using fluorescein-labeled oligonucleotide probes (INFORMEBER Probe; Ventana). Specimens displaying EBER nuclear expression in > 20% of the tumor cells were considered EBER-positive.

Primary GC dataset collection

We collected five transcriptome datasets, TCGA-STAD (T = 375), GSE84437 (T = 433),⁵⁰ GSE62254/ACRG (T = 300),¹⁴ GSE26942 (T = 217),⁵¹ and GSE13861 (T = 90),⁵² from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases; 1415 samples with clinical data were included in the study. The three cohorts ACRG, GSE26942, and GSE13861 contained information on ACT. In addition, we collected a transcriptome dataset of patients who responded to neoadjuvant chemotherapy (NAC) from our previous study.⁵³ Gene expression data was log2 transformed if necessary.

Patients with GC who received immunotherapy dataset collection

Two transcriptomic datasets with immunotherapy response information were included in our study: patients with GC (T = 45) treated with anti-PD1 therapy⁵⁴ and patients with GC (T = 24) treated with anti-PD1/PDL1/CTLA4 therapy.⁵⁵ Gene expression data in all datasets was log2- and z-score-transformed, if necessary.

16S rRNA gene sequencing

A total of 251 GC tissue pairs were selected for microbiological analysis. Total genomic DNA samples were extracted using the OMEGA Soil DNA Kit (M5635–02) (OmegaBio-Tek, Norcross, GA, USA), following the manufacturer’s instructions, and stored at −20°C for further analysis. The quantity and quality of extracted DNA were measured using a NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively. PCR amplification of the bacterial 16S rRNA gene V3-V4 region was performed using the forward primer 338F (5’-ACTCCTACGGGAGGCAGCA-3’) and the reverse primer 806 R (5’-GGACTACHVGGGTWTCTAAT-3’). Sample-specific 7-bp barcodes were incorporated into the primers for multiplex sequencing. The PCR components contained 5 μl of buffer (5×), 0.25 μl of Fast pfu DNA Polymerase (5 U/μl), 2 μl (2.5 mM) of dNTPs, 1 μl (10 uM) of the forward and reverse primers, 1 μl of DNA Template, and 14.75 μl of ddH2O. Thermal cycling consisted of initial denaturation at 98°C for 5 min, followed by 25 cycles of denaturation at 98°C for 30 s, annealing at 52°C for 30 s, and extension at 72°C for 45 s, with a final extension cycle for 5 min at 72°C. The PCR amplicons were purified using Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and quantified using a Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). After the individual quantification step, the amplicons were pooled in equal amounts, and pair-end 2 × 250bp sequencing was performed using the Illumina NovaSeq platform with the NovaSeq 6000 SP Reagent Kit (500 cycles) at Shanghai Personal Biotechnology Co., Ltd (Shanghai, China).

16S rRNA data analysis

Sequence data analyses were performed using QIIME2.⁵⁶ Taxonomic information for each ASV is obtained by aligning it with reference sequences in the Greengenes database. Alpha diversity indices were calculated using the ASV table in QIIME2. Beta diversity analysis was performed to assess the structural variation of microbial communities across samples using Bray-Curtis metrics and was analyzed via principal coordinate analysis (PCoA). The significance of differential microbiota among groups was assessed using PERMANOVA. SparCC⁵⁷ was used to compute the co-occurrence between microbiota. The network was visualized and analyzed using the iGraph package. To measure and identify the hub microbiota in the network, Kleinberg’s hub centrality scores were calculated for each network and scaled across all networks. O2PLS modeling was then employed to delineate the relationships between microbiota and immune-related genes; the main immune-related genes and core microbial genera with strong correlations were selected. Microbial functions were predicted using PICRUSt2⁵⁸ in the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/) database.

RNA extraction

A total of 94 tumors were analyzed using RNA-seq. Total RNA was isolated using Trizol reagent (Invitrogen Life Technologies), and the concentration, quality, and integrity of RNA were determined using a NanoDrop spectrophotometer (Thermo Scientific). Three micrograms of RNA were used as the input material for RNA sample preparation. Sequencing libraries were generated by first purifying mRNA from the total RNA using poly T oligo-attached magnetic beads. Fragmentation was performed using divalent cations at elevated temperatures in the Illumina proprietary fragmentation buffer. Then, first-strand cDNA was synthesized using random oligonucleotides and SuperScript II reverse transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities, and the enzymes were removed. After adenylation of the 3′ ends of the DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare them for hybridization. To select 400–500 bp cDNA fragments, library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, CA, USA). DNA fragments with ligated adaptor molecules at both ends were selectively enriched using the Illumina PCR Primer Cocktail in a 15 cycle PCR reaction. The products were purified (AMPure XP system) and quantified using an Agilent High-Sensitivity DNA Assay on a Bioanalyzer 2100 system (Agilent Technologies). The library was sequenced using a NovaSeq 6000 platform (Illumina).

RNA-Seq data analysis

Quality control and raw data processing were performed using FastQC (v0.11.9) and Cutadapt (v1.15).⁵⁹ The reads were aligned to the reference genome (hg38) using HISAT2 (v2.0.5).⁶⁰ HTSeq (0.9.1)⁶¹ was applied to obtain FPKM. The FPKM values were transformed into TPM values. Gene expression was then log2-transformed. The EBV reference genome was downloaded from https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_027943465.1/and the EBV genome was mapped using the STAR software (v2.7.9a),⁶² generating the read counts for the EBV genome. The EBV load of each sample was defined as the total read counts of the EBV genome.

Deconvolution of GC transcriptome data using the MCP-counter

To characterize the immune characteristics of the GC transcriptome, the abundance of eight immune cells (T cells, CD8⁺ T cells, cytotoxic lymphocytes, NKs, B cells, monocytes, dendritic cells, and neutrophils) and two stromal cells (endothelial cells and CAFs) were estimated using the MCP-counter.⁶³

Consensus clustering analysis

Microbial genera that were expressed in more than 60% of the patient samples were selected for clustering analysis. Consensus clustering was conducted using the ConsensusClusterPlus package with the Pearson correlation as the distance measure.⁶⁴ The consensus matrix for k = 3 showed clear separation between the subtypes. In addition, the consensus CDF and delta plots showed that the relative change in the area under the cumulative distribution function curve increased the most from two to three subtypes, whereas the others exhibited no significant increase. Furthermore, the three subtypes were significantly associated with patient survival (Figure S3). Taken together, the Msubtypes were defined using consensus clustering with k = 3. Subtype-specific genera were defined as those that were highly enriched in a specific subtype (Adj p-value <0.05, Log2 FC > 2, and AUC > 0.7).

Functional enrichment analysis

To further analyze biological characteristics of different subtypes, we performed gene set enrichment analysis (GSEA) to identify enriched pathways from the KEGG, Biocarta, Reactome, PID, Wikipathways, and Hallmarker databases associated with the Msubtypes with at least 10 overlapping genes. GSEA was performed on the log2 (FC) values from the differential abundance analysis between Msubtypes and others in the transcriptome. Pathways with an FDR < 0.1 and a normalized enrichment score (NES) > 0 were considered significantly enriched.

GC subtype classification and prediction evaluation

We used the support vector machine (SVM)⁶⁵ algorithm to divide samples from the published transcriptome cohorts into subtypes. The SVM algorithm was performed using e1071 package. A transcriptome dataset was used as the training set to train the model, and 10-fold cross-validation was performed to evaluate the accuracy of the model.

Subgroup correlation analysis of independent cohorts

An unsupervised subclass mapping method (SubMap) (GenePattern module “SubMap”) (https://cloud.genepattern.org) was used to identify similar subgroups between independent cohorts despite their clinical differences.⁶⁶

Quantification and statistical analysis

Statistical analyses were performed using R (v. 4.2.1) unless otherwise stated. Standard statistical tests were used to analyze data from different groups. For categorical variables, Fisher’s exact test was applied, while continuous variables were evaluated via Wilcoxon rank-sum test to compare differences between two groups. The Spearman’s rank correlation coefficient was used to evaluate the correlations between two continuous variables. For survival analysis, the Kaplan-Meier method and log-rank test were used to compare survival distribution. Additionally, Cox analysis was applied to evaluate the correlation between survival and different variables and to calculate hazard ratio (HR) values. All statistical tests were two-sided, and p < 0.05 was considered statistically significant.

Author contributions

J.F.J., H.J.W. and X.F.X. conceived the study. X.F.X. contributed to biospecimen processing. G.J.W., H.J.W., X.J. and W.J.J. analyzed the data. G.J.W., X.F.X., X.J. and L.M. contributed to collection of clinical samples and data and critical review of data. G.J.W., H.J.W., and H.J.W. drafted the manuscript. J.F.J., H.J.W., H.J.W., X.F.X., X.J. and G.J.W. critically revised the manuscript. J.F.J., H.J.W., H.J.W., X.J., X.F.X., T.W., Y.Z., W.J.J., L.M. and G.J.W. supervised the study. J.F.J., T.W., H.J.W., and Y.Z., obtained funding. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The sequencing data generated by this study is available through the Genome Sequence Archive for Human (GSA-Human): HRA005570 (16s rRNA, https://ngdc.cncb.ac.cn/gsa-human/submit/hra/subHRA007973/finishedOverview.) and HRA005478 (RNAseq, https://ngdc.cncb.ac.cn/gsa-human/submit/hra/subHRA007857/finishedOverview).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19490976.2024.2369336

Additional information

Funding

This work was supported by grants from the joint fund for key projects of National Natural Science Foundation of China [U20A20371], National High Technology Research and Development Program of China [863 Program, No. 2014AA020603], National Natural Science Foundation of China [No. 82272889, 82073312, 82103528, 82203579], Beijing Municipal Administration of Hospitals Incubating Program [PX2019040, PX2019039], Clinical Medicine Plus X-Young Scholars Project, Peking University [the Fundamental Research Funds for the Central Universities, PKU2020LCXQ001], 2021 Taihu Talent Program Top Medical Expert Team [2021-THRC-DJ-PWK], Sanya Yazhou Bay Science and Technology City Doctoral Research Innovation Fund Project [HSPHDSRH-2023-11-001].