Gut bacteria-derived sphingolipids alter innate immune responses to oral cholera vaccine antigens

Domain: neuroscience
Published: 2025-12-12
Editorial outcome: major_revision

Abstract
The degree of protection conferred after receiving an oral cholera vaccine (OCV) varies based on age, prior exposure to Vibrio cholerae, and unknown factors. Recent evidence suggests that the microbiota may mediate some of the unexplained differences in oral vaccine responses. Here, we use metagenomic sequencing of the fecal microbiota at the time of vaccination and relate microbial features to immune responses after OCV using a reference-independent gene-level method. We find that the presence of sphingolipid-producing bacteria is associated with the development of protective immune responses after OCV. We test these associations by stimulating human macrophages with Bacteroides xylanisolvens metabolites and find that sphingolipid-containing extracts increase innate immune responses to OCV antigens. Our findings demonstrate a new analytic method for translating metagenomic sequencing data into strain-specific results associated with a biological outcome, and in validating this tool, we identify that microbe-derived sphingolipids impact immune responses to OCV antigens.

Introduction
Cholera is a severe diarrheal disease caused by the bacterium Vibrio cholerae, and millions of cases occur annually1. Despite increased surveillance and ongoing efforts to provide safe water to vulnerable populations, there has been a recent re-emergence or new outbreaks of cholera in Haiti, Lebanon, and Syria, in addition to ongoing endemic spread in over 50 countries in Asia and Africa2–4. Oral cholera vaccines (OCV) are important tools in cholera control, but provide a short duration of protection and are less effective in children under five years of age5.

Long-term immunologic protection from cholera is incompletely understood. The vibriocidal titer is a commonly used measure of the humoral response to V. cholerae infection, and although higher titers are associated with protection from disease, titers wane prior to loss of protective immunity6,7. Antibody responses to natural infection are largely directed toward cholera toxin (CT) and the V. cholerae lipopolysaccharide (LPS)8. In contrast with other toxin-mediated infections, CT-specific antibody responses are associated with only short-term protection from cholera5,9,10. The V. cholerae O1 serogroup is responsible for the majority of epidemic cholera, and serogroup specificity is defined by the composition of the O-specific polysaccharide (OSP) of the LPS. The O1 serogroup is subdivided into serotypes, and Inaba and Ogawa serotypes are the most dominant, and differ only by a single methyl group of the terminal sugar on the O-antigen11. After infection, OSP-specific memory B cells (MBCs) are detectable for up to one year12, and upon exposure to V. cholerae, these cells are presumed to differentiate rapidly into antibody-secreting plasma cells to produce antibodies that act at the mucosal surface to protect from infection9,13–16. This is supported by longitudinal studies of household contacts of cholera patients, in which persons with elevated OSP-specific MBCs at the time of exposure to V. cholerae were protected from infection9,17. Recent in vivo studies established that OSP-specific antibodies reduce V. cholerae colonization in a motility-dependent manner, and thus, are thought to be a mediator of protection in human infection18. A collection of genetic, genomic, and laboratory-based studies has established that innate immune responses occur after acute infection and have an important role in shaping long-term memory B cell responses13,16,19–22. Innate immune responses contribute to the development and maintenance of OSP-specific MBCs that later reactivate by promoting upregulation of effector proteins and stimulating follicular helper T cell development for the provision of help to support B cell development, proliferation, and differentiation in germinal centers22–24.

On a population level, vaccination with a whole-cell killed OCV stimulates an OSP-specific MBC response in half of OCV recipients12,25–27. While demographic factors such as age correlate with the likelihood of developing a protective vaccine response, biological mechanisms of these response differences have not been identified26,28,29. The concept that gut microbiota differences may mediate differences in host response to OCV is supported by associations between the microbiota and responses to other oral vaccines30–33, the observation that immune responses to OCV were reduced in recipients with bacterial intestinal overgrowth34, and our prior study that detected associations between the gut microbiota and OSP-specific MBC responses to OCV in a separate cohort25. To understand how gut microbes present at the time of vaccination could impact the host response to OCV, we applied a novel gene-level association method to metagenomic sequencing data from the microbiota at the time of OCV receipt. This analytic technique was selected and adapted to our analysis to identify specific microbial isolates associated with our outcome of interest, and this enabled us to experimentally test our sequencing study findings. We found that the presence of Bacteroides species capable of producing sphingolipids was correlated with the development of protective immune responses after OCV. To validate these findings, we exposed human macrophages to Bacteroides-derived sphingolipids in vitro and found a greater innate immune response to V. cholerae antigens, which suggests that gut microbial metabolites may impact immune responses to OCV.

Results
Immune responses to OCV are not correlated with microbiome composition
Eighty-nine participants received a whole-cell killed V. cholerae OCV with recombinant cholera toxin B and underwent measurements of post-vaccine immune response in Dhaka, Bangladesh (Fig. 1a)27. Stool collected on the day of OCV administration underwent DNA extraction and metagenomic sequencing. Twenty-two participants were aged 2-5 years (11/22, 50% female), 22 participants >5-17 years (12/22, 54% female), and 45 participants ≥18-46 years (32/45, 71% female). Immune responses after OCV were consistent with prior studies of whole-cell killed vaccine responses from populations living in cholera-endemic areas12,25,29,35. As described in the parent study27, vibriocidal responses to V. cholerae O1 Ogawa and Inaba serotypes increased after vaccination as expected, and only some participants had an increase in MBC response to vaccination, as observed in prior studies26,36. Among participants with MBC increases after vaccination, the mean increases in CT- and OSP-specific IgG-producing MBCs before and after vaccination were significantly increased, while CT- or OSP-specific MBCs producing IgA were not, as previously observed27. We found that OCV recipients were more likely to develop IgG-directed OSP-specific MBC responses than IgA OSP-specific responses (Fig. 1b). Among the immunologic parameters measured in this study, we focused primarily on V. cholerae OSP-specific MBC results because the presence of detectable OSP-specific MBCs is associated with long-term protection from infection, and a potential mechanism for OSP-specific antibody protective immunity has recently been established9,17,18. To assess microbiome differences among vaccine recipients, participants were grouped according to the development of an MBC response between pre- and post-vaccination time points for each (CT or OSP) IgA and IgG response. Individuals with increased MBC counts following vaccination were defined as vaccine responders (R), and participants who had unchanged or decreased MBC measures after vaccination were defined as vaccine nonresponders (NR).

We compared microbial taxonomic diversity measures by age and sex and found no significant differences by age (Supplementary Fig 1a-b). Shannon diversity (a measure of both species diversity and abundance) demonstrated a difference in microbial diversity by sex (Supplementary Fig 1c-e), driven by fewer Lactobacillales, Veillonellaceae, and Eubacteriaceae, and more Erysipelotrichaceae, Ruminococcaceae, and Lachnospiraceae in women, consistent with microbiome-specific sex differences noted in other populations37. No significant differences were found by Shannon Index or by Bray-Curtis community structure between vaccine responders and nonresponders for the four antigen-specific MBC measures (IgG- and IgA-producing OSP-specific or CT-specific MBCs, Supplementary Fig 2). A rarefaction curve, which shows the number of species observed in the study based on the number of sequences sampled, demonstrated that our sampling method captured the majority of the bacterial diversity in our study population (Supplementary Fig 3).

Co-abundant gene (CAG) groups differentiate the microbiota of OCV responders from nonresponders
To address the high dimensionality of the metagenomic data and examine detailed features of the microbiota, we used gene groups as the fundamental unit of analysis. CAGs are likely to be encoded on a shared genetic element (e.g., chromosome or plasmid) across the organisms present in the cohort, and grouping genes together in this way has been successfully used to infer the relationship between specific genes and biological outcomes38–40. We identified CAGs that were associated with MBC responses and predetermined a significance threshold of a False-Discovery Rate (FDR)-adjusted q-value of ≤0.1, based on prior studies associating metagenomic data with biological responses (Fig. 1c)41. We calculated significant CAGs for each vaccine response measure separately (IgA OSP, IgA CT, IgG OSP, and IgG CT), and identified CAGs associated with an MBC response or lack of MBC response in each group. A total of 4147 CAGs were identified, and 323 CAGs containing ≥2 genes met our significance threshold across the four vaccine response measures (Supplementary Table 1). The greatest number of significant CAGs that mapped to bacterial genomes were associated with IgA OSP-specific MBC responses, and most of these CAGs were associated with a higher immune response (Supplementary Table 1). CAGs mapped to assembled bacterial genomes demonstrated that many CAGs aligned to multiple strains (Supplementary Fig 4). Many genes composing significant CAGs have unknown functions or hypothetical proteins (based on the results of the eggNOG functional annotation and NCBI nr/blastp database (Supplementary Data 1)). Because all human sequences were removed from metagenomic data prior to CAG construction, this represents an enormous underexplored and unidentified microbial diversity of possible clinical significance that may be bacterial, fungal, or viral in origin. Among annotated genes in significant CAGs, most encoded nonspecific cellular functions, were related to amino acid metabolism, or were phage-associated.

CAG group association with OCV responses was strain-specific
Because CAGs mapped to multiple reference strains, a priority score was assigned to each bacterial strain detected from participants in this study (see Methods). This process resulted in the identification of top bacterial strains associated with each MBC response, and many were from diverse taxonomic groups (Fig. 2a-b, Supplementary Fig 5a-b, Supplementary Table 2, list of strains with priority scores, Supplementary Data 2-5). Typically, only a small fraction of strains within a given taxonomic group had a high priority score, indicating that the association between MBC responses and gut microbes may be strain specific (Fig. 2c, Supplementary Fig. 5c). Among 412 total bacterial genera detected in the stool of the study population, only 34 contained five or more strains that mapped to a significant CAG (Fig. 2c, Supplementary Fig. 5c). To demonstrate the relationship between high-priority and unscored strains within one species, we identified Bacteroides xylanisolvens strains that were detected in the stool of our study population. To put these strains in the context of the observed diversity of this species, we next created a phylogenetic tree using NCBI genomes of all existing B. xylanisolvens strains (Supplementary Fig. 6). The four most highly ranked B. xylanisolvens strains identified by our analysis were not consistently clustered within the overall population of known B. xylanisolvens (Fig. 2d). This suggests that even closely related strains differ enough in gene content to be associated with different biological effects, and potential biological effects are not restricted to a single monophyletic lineage or clade. One exception to this finding is that most strains of Clostridioides difficile detected in our study population were highly positively associated with this MBC response (Fig. 2c). Among the 1,080 strains identified in study participant stool that positively correlated with OSP IgA MBC responses, Bacteroides species (n=31) and Prevotella copri (n=53) were disproportionately represented among the highest priority scored strains (Fig. 2b, Supplementary Data 2). Strains positively associated with OSP IgG responses were taxonomically varied and included members of Bifidobacterium and Ruminococcus/Blautia genera and several Enterobacteriaceae (including Escherichia coli, Enterobacter, and Klebsiella species) (Fig. 2a, Supplementary Data 3). Our results indicate that associations with biological outcomes were strain-specific, and thus, strain-specific methods of identifying these relationships are needed to select appropriate strains for experiments to test and validate a correlative host-microbe interaction. These findings support the concept that taxonomic groupings alone are unreliable markers for specific genes or associations with biological effects.

The microbiota of vaccine responders was characterized by bacterial species with sphingolipid-producing potential
The identification of bacterial strains associated with OSP-specific IgA-producing MBC responses featured genera with the ability to produce sphingolipids, such as Bacteroides and Prevotella (Fig. 2b-d). We were intrigued by this because animal studies have demonstrated that Bacteroides-derived sphingolipids reduce intestinal mucosal inflammation42,43. This was also a notable finding because bacteria rarely produce their own sphingolipids, and these bacterial-derived lipids are unique due to their odd-chain structures, in contrast to even-chain sphingolipids made by eukaryotes44–46. Because the metabolic capacity for sphingolipid production can vary widely across members of Bacteroides47, we used the strain-level information identified by the CAG analysis to measure the relative abundance of organisms encoding the serine palmitoyltransferase (spt) gene, the bacterial enzyme required for de novo assembly of microbe-derived sphingolipids43. We found that vaccine responders stool had a higher abundance of Bacteroides spt compared to nonresponder stool (Fig. 3a). For these reasons, we investigated the relationship between vaccine response and Bacteroides-derived lipids using a targeted, quantitative mass spectrometry-based lipidomic analysis of study participant stool. Several sphingolipid derivatives, including ceramides, were more abundant among vaccine responders (Fig. 3b and Supplementary Table 3). Ceramides are composed of a sphingoid base and a fatty acid, and bacterial sphingolipids are derived from this ceramide backbone48. Thus, ceramides are a structural backbone of all sphingolipids, and ceramide detections approximate sphingolipid content. There were also specific species of non-sphingolipid lipids that were more abundant in vaccine responders, including diacylglycerol, cholesterol esters, and free fatty acids, while triacylglycerols were found to be more abundant in vaccine nonresponders (Supplementary Table 3). While this stool lipidomics analysis does not fully distinguish between host- or bacterial-derived lipids, we were interested in differences in lipid content between OCV responders and nonresponders because in our prior study, fecal extracts from these two groups generated different innate immune responses in our model of OCV response25, and our metagenomic analysis above suggested that bacterial-derived sphingolipids could be contributing to these findings.

Bacteroides-produced sphingolipids blunted inflammatory cytokine production in human macrophages
Since innate immune responses impact the development of subsequent adaptive immune responses, including the differentiation and proliferation of MBCs, we next prepared to test whether Bacteroides sphingolipids impacted mucosal innate immune responses in our established human macrophage model of OCV response25,49. We first selected a spt-positive Bacteroides xylanisolvens (Bx) strain from our catalog of isolates derived from the stool of the study population and used this for additional experiments. We next generated lysates from Bx culture grown with or without the presence of myriocin, an spt inhibitor42. We then fractionated the Bx lysates to evaluate if lipid or non-lipid components of the lysate impacted immune responses in our model. Fractionation using the Bligh–Dyer method yielded culture metabolites, non-soluble cell membrane, and lipid components50. To ensure that myriocin had reduced sphingolipid production in the myriocin-treated Bx lipid extracts, we next measured lipid composition and abundance in the Bx lipid extracts using hydrophilic interaction liquid chromatography (HILIC) followed by ion mobility-mass spectrometry51. Bx lipid extracts exposed to myriocin had reduced ceramide abundance when compared to Bx lipid extracts not exposed to myriocin (Supplementary Fig. 7). We also sought to identify specific ceramides made by Bx and compared these to lipids found in vaccine responders' and nonresponders' stool. Lipid identification from stool has significant limitations, especially for full resolution of branching and structural patterns, and this is further complicated by incomplete knowledge about sphingolipids made by bacteria. Thus, full resolution of lipid structures from measurements in stool was not possible. However, we identified two bacterial-derived sphingolipids detected with high confidence, both in stool and in the Bx sphingolipid extract, CER(16:0) and CER(18:1). Both of these ceramides were more abundant in vaccine responders and more abundant in B. xylanisolvens lipid extracts untreated with myriocin (Supplementary Table 3 and 4).

After confirming that Bx lipid fractions treated with myriocin had reduced sphingolipids (Supplementary Fig. 7), we next examined the effect of sphingolipid-reduced (SL-reduced) Bx lysate and found significantly less interleukin-6 (IL-6) production compared to the untreated, sphingolipid-containing (SL-containing) lysate when applied to THP-1 macrophages (Supplementary Fig. 8). We then tested the fractionated Bx. Cytokine responses in our human macrophage model stimulated with metabolite and non-soluble fractions of Bx culture did not differ between SL-reduced and SL-containing fractions (Fig. 4a-b). We then tested the SL-reduced and SL-containing lipid fractions, and we identified that IL-6, tumor necrosis factor-α (TNF-α), and IL-8 were decreased in the presence of Bx sphingolipids (Fig. 4c). The production of IL-6 by cells exposed to Bx SL-reduced lipids was nearly 6-fold compared to cells exposed to the SL-containing lipid fractions. We also evaluated cell viability in the presence of lipid treatment and found that lipid extracts did not alter cell viability (Supplementary Fig. 9).

Bacteroides-produced sphingolipids increased innate immune responses to V. cholerae antigens
We observed an increase in innate cytokine responses when Bx SL-reduced lipid extracts were applied to human macrophages, and this was consistent with our prior findings that fecal extracts from vaccine nonresponders increased IL-6 production in this model25. Because this response was triggered by vaccine nonresponders only, we hypothesized that an increase in resting innate immune activation could dampen the immune response to vaccine antigens. To test this, we pretreated THP-1-derived human macrophages with SL-reduced and SL-containing lipid fractions and then exposed them to heat-killed V. cholerae strain JBK70, an A-B toxin negative vaccine strain, as a proxy for OCV, as in prior studies49. Macrophages that were pretreated with SL-containing lipid fractions and then stimulated with JBK70 were found to produce increased IL-6 and TNF-α (Fig. 5a). In contrast, macrophages pretreated with SL-reduced lipid fractions had lower cytokine responses to JBK70. To minimize the possibility that this effect was caused by non-SL lipid species, we next depleted phospholipids from the lipid extracts using a mild alkaline hydrolysis treatment42. This change potentiated the macrophage response to sphingolipids, with further increases in the magnitude of IL-1β, IL-6, IL-12p40, and TNF-α production in response to JBK70 stimulation (Supplementary Fig. 10). To test for this effect in another spt-positive Bacteroides isolate from our study population, we selected a strain of Bacteroides koreensis, a species identified as highly associated with OSP-specific MBC OCV responses (Fig. 2b). Our findings were replicated using this strain (Supplementary Fig. 11).

To further query the effect of bacterial sphingolipids on human macrophages, we next performed a gene expression analysis on RNA extracted from THP-1 macrophages preconditioned with Bx SL-containing or SL-reduced lipid fractions with and without exposure to heat-killed JBK70 (Fig. 5b). Using a Gene Set Enrichment Analysis52, we found that preconditioning macrophages with SL-reduced Bx lipids prior to JBK70 exposure increased activation of several inflammatory pathways compared to cells preconditioned with SL-containing lipids (Supplementary Fig. 12a). Compared to the SL-containing fraction, SL-reduced fraction preconditioning triggered pathways associated with interferon-α, interferon-γ, and TNF-α cytokine responses, and IL-6—Janus kinase (JAK)—signal transducer and activator of transcription (STAT3) pathways. After preconditioning and subsequent exposure to heat-killed JBK70, TNF-α signaling via the nuclear factor kappa-light-chain-enhancer of B cells (NF-κB) pathways was upregulated independent of SL exposure (Supplementary Fig. 12b, c, individual genes listed in Supplementary Table 5).

To identify immune activation pathways that distinguished SL-reduced from SL-containing lipid stimulation of human macrophages after JBK70 stimulation, we repeated our analysis to exclude gene expression specific to heat-killed JBK70 that was independent of bacterial sphingolipid exposure. Here we found that human macrophage gene expression after SL-containing lipid treatment and JBK70 stimulation was enriched in TNF-α signaling via NF-κB, IL-6—JAK—STAT3 signaling, inflammatory response, epithelial/mesenchymal transition, and complement pathways, compared to the macrophage response after treatment with SL-reduced lipids (Fig. 5c). The 19 genes differentially associated with enriched pathways are shown in Fig. 5d and Supplementary Fig. 13. Several of the genes were associated with multiple innate inflammatory pathways, including inhibin A (INHBA), C-X-C motif chemokine ligand 7 (CXCL7), and CXCL1.

Lastly, to validate our findings in THP-1 macrophages using additional immune cell types, we tested the effect of OCV antigens with and without sphingolipid exposure in human peripheral blood mononuclear cells (PBMCs). We preconditioned PBMCs with B. xylanisolvens sphingolipid fractions and replicated our findings of increased IL-6 after JBK70 stimulation, as well as TNF-α in some donors (Supplementary Fig. 14). These data in combination support our hypothesis that the correlations we observed between sphingolipid-producing bacteria and vaccine response in humans may be partly attributable to sphingolipid-driven modulation of the innate immune responses to OCV antigens.

Methods
Study design, subject enrollment, and sample collection
Adult and child participants over 2 years of age and <60 years of age living in urban Dhaka, Bangladesh, were enrolled through the International Centre for Diarrheal Disease Research, Bangladesh (icddr,b) between Jan 1, 2017, and Dec 30, 2017, if they had not previously been vaccinated with an OCV and had no major comorbid conditions89. After obtaining written informed consent from patients or guardians, study participants received two doses of WC-rCTB (Dukoral), per the package insert, on study Day 0 and Day 14. All time points are reported in reference to Day 0, the day of administration of the first dose. Blood specimens were collected from the participants on Day 14 (two weeks prior to the first dose of vaccination), Day 2, Day 7, and Day 44. Stool samples were collected on Day 0 on site at the icddr,b in cryovials and immediately frozen at −80°C. Vibriocidal titer was assessed at each of these time points; OSP-specific and CT-specific MBC were measured at Day 14 and Day 44. Demographic data were collected at enrollment, and stool samples were collected at baseline prior to the first dose of WC-rCTB. Additional details and methods for vibriocidal titers, serum antibodies, and MBC measures are described in the parent study27.

Stool processing and DNA preparation
Stool samples were selected for sequencing if the participant had MBC ELISPOT results for both pre- and post-vaccination timepoints for at least two out of four MBC responses. Application of this criterion resulted in sequencing the stool from 90% of study participants. Stool was collected in cryovials and immediately frozen at −80°C and shipped to US investigators. At the University of Washington, stool microbial DNA was extracted using the PowerSoil DNA isolation kit (Qiagen) with a modified protocol as previously described in ref. 78. Briefly, the stool was thawed on ice, and approximately 100 mg of stool was added to bead-beating tubes. The stool was treated with C1 solution, heated at 65°C for 10 min, 95°C for 3 min, and vortexed for 10 min. Samples were treated using C2–C5 solutions provided by the Powersoil kit. The resulting DNA was eluted in DNase- and RNase-free water; quality and quantity were measured using a NanoDrop ND-1000 (Thermo Scientific).

Library preparation, sequencing, read processing, and CAG grouping
Library preparation was performed using uniquely dual-indexed libraries using the Nextera DNA Flex kit. Sequencing was performed on the NovaSeq 6000 PE150 platform using standard Illumina adapters targeting a minimum of 10 M reads per sample with 150 bp pair-end reads in one sequencing run. Fastq files were processed using geneshot, a pipeline developed to streamline the analysis of gene-level metagenomic data and its association with experimental metadata38. Using default settings for this pipeline (v0.9), raw WGS data from each individual specimen was quality-filtered and de novo assembled using Megahit90; human reads were excluded using human reference genome GRCh37 hg19, and protein-coding sequences were predicted using Prodigal91 and deduplicated across all specimens using linclust (at 90% amino acid similarity); coding sequences were annotated with taxonomic labels using DIAMOND (by alignment against NCBI RefSeq 01/22/2020) and with functional labels using eggNOG-mapper v2 (accessed 05/22/2024)92 and NCBI blastp (accessed 05/22/2024) with database nr_clustered (experimental)); WGS reads from each specimen were aligned against the de novo generated gene catalog with DIAMOND; coding sequences were grouped by co-abundance into CAGs39. The centroid was selected for each deduplicated gene using the mmseqs2 algorithm93. The depth of sequencing used for CAG generation was an average of the depth of sequencing across all genes in each CAG. The relationship between CAGs containing ≥2 genes was correlated with participant vaccine response measures using the corncob algorithm implementation of beta-binomial regression and accounting for multiple comparisons using the Benjamini-Hochberg false-discovery rate correction procedure94,95. CAGs containing only one gene were excluded to reduce the total number of tested hypotheses, considering the relative lack of biological interpretability for singleton coding sequences in the context of gene-level metagenomic analysis. Top CAGs were defined prior to the analysis as those with a q-value ≤0.1.

Assessment of gut microbial communities
Reads from all participant samples were used to estimate the relative abundance of taxonomic groups using metaphlan296, which uses the alignment of raw reads against a set of previously identified marker genes that are conserved across taxonomic groupings of genomes. The output of metaphlan2 was merged across all samples, and non-bacterial taxa were removed. Alpha diversity measures by Shannon Index and beta diversity measures by Bray-Curtis Dissimilarity Index were calculated using divnet (v3.6) in R (v4.0.2)97. Principal coordinate of analysis (PCoA) plots to visualize these results were generated using the R package vegan. Rarefaction on species-level identification was calculated using mothur (v.1.44.1) command rarefaction.single() with 10 iterations/randomizations with sequence sampling frequency set at 10098.

CAG alignment to bacterial genomes
Top CAGs found to be significantly associated with vaccine response measures were aligned to bacterial genomes to identify the strains containing the most top CAGs. A strain is defined here as a single reference genome present in the RefSeq database. A total of 61,918 bacterial reference genomes were selected from the National Center for Biotechnology Information prokaryote genome database using DIAMOND v0.9.10 via the AMGMA (Annotation of Microbial Genomes By Microbiome Association) workflow [https://www.ncbi.nlm.nih.gov/genome/browse#!/prokaryotes/] and [https://github.com/fredhutch/amgma]39,99,100. We used all bacterial sequences in the database apart from bacterial sequences with low assembly quality and completeness (chromosomes and scaffolds). We included all complete genomes and 85% of contig-level assemblies and excluded 15% of contig-level assemblies with the lowest levels of completeness (>200 contigs). We expect that genomes with lower assembly quality will be more difficult to detect from metagenomic datasets using this or any other analytical approach. Because this database contains multiple representatives from many species and genera, it is expected that homologous copies of many CAGs may be found in multiple genomes. For that reason, alignments of a single CAG to multiple genomes were all retained for the analysis. Using the alignment of coding sequences from the metagenomic analysis against this set of partial and full genomes, a relative abundance measure was generated for each genome in each sample by summing the proportion of reads from a given sample that aligned to those genes found in each genome. Using that strain-level relative abundance measure, we applied the same statistical approach (beta-binomial modeling using the corncob algorithm) to estimate the coefficient of association with vaccine response measures on a per-genome basis. Our gene-based analysis described here is reference-free due to our method of CAGs generation and quantification across metagenomes as a function of vaccine response; the CAG alignment to bacterial genomes is based on the genomes contained in the RefSeq database.

Development of a CAG-based strain identification method
Once CAGs were aligned to bacterial genomes38,39, we next quantified the association between these strains and vaccine response measures. For this purpose, we created a priority score designed to (1) account for the strength of the association between top CAGs aligning to a genome and the vaccine response measure of interest, (2) preserve the directionality of the association between the vaccine response measure and top CAGs, (3) be independent of CAG size, in order to avoid bias toward large CAGs containing more genes that could potentially align, and (4) discriminate between strains within the same species.
For a given strain with n number of significant CAGs aligning to strain genomes, this priority score was calculated as follows:
PS = [ (Sum of (NumGenes Align CAG1 * ± est coefficient CAG1^2 + wald CAG1^2)) / (Sum of NumGenes Align CAGn) ] * [ (Mean est coefficient StrainX^2 + mean wald StrainX^2) / (Mean est coefficient StrainX^2 + mean wald StrainX^2) ]
The estimated (est) coefficient reflects the abundance of a CAG in responders compared to nonresponders for each vaccine response measure, and the Wald score integrates the standard error, which reflects the uniformity of the association39. The ± in the first term corresponds to the sign of the CAG’s estimated coefficient. The number of genes aligning from a CAG to the strain genome is represented in the estimated coefficient and Wald score, creating a weighted measure of CAG-strain association, which is summed in the first term’s numerator for all top CAGs with genes aligning to the strain. By normalizing this weighted sum in the numerator of the first term by the unweighted number of genes aligning from top CAGs, we minimize bias based on CAG size. In the second term, the score is further weighted by integrating the number of genes aligned to reference strains without regard for CAG grouping, and thereby discriminates between strains that may have otherwise identical first terms due to the alignment of the same number of genes from the same top CAG(s). A + priority score indicates a positive association between the vaccine response measure and strain, and a - indicates an inverse relationship.

Log10 (P-value), which was obtained using multiple Mann–Whitney tests.

Bacterial strains, growth conditions, and identification
Bacteroides xylanisolvens was isolated from human fecal samples from the study participants27 using brain-heart infusion (BHI) agar (BD) supplemented with 5% sheep’s blood (Hemostat Labs) and Laked Brucella Blood agar with Kanamycin and Vancomycin (LKV) plates (Anaerobe Systems). Bacteroides koreensis was isolated using Gifu Anaerobic Media (HiMedia) supplemented with porcine hemin (10 mg/L; MP Biomedicals) and menadione (5 mg/L; MP Biomedicals) and LKV plates. Bacterial identity was confirmed using ~1000 bp Sanger sequencing of the 16S rRNA gene with >97% identity. Liquid cultures of Bacteroides strains were prepared in BHI broth (BD Diagnostic Systems) supplemented with 0.005 mg/mL porcine hemin, 0.5 g/L L-cysteine (Sigma), and 0.05 ng/mL Vitamin K1 (Sigma) and cultivated at 37°C in a vinyl anaerobic chamber (Coy) with <50 ppm O2. Glycerol stocks were prepared using BHI supplemented media, 20% glycerol, and stored in -80°C ultra-low freezers. Bacteroides strains were maintained on Tryptic Soy Agar with 5% sheep’s blood plates and grown in BHI supplemented media for 24–48 h. The V. cholerae strain JBK70 was gifted by Dr. Jim Kaper49. V. cholerae was maintained in lysogeny broth (LB) or LB agar at 37°C in aerobic conditions.

Mammalian tissue culture cell line and maintenance
THP-1 cells were purchased from and authenticated by ATCC (ATCC cat# TIB-202). PBMCs were purchased from Bloodworks Northwest. Four PBMC donors were acquired. All cells were maintained in RPMI-1640 media (containing L-glutamine and 25 mM HEPES) (Corning) supplemented with 10% fetal bovine serum (Fisher) and penicillin and streptomycin (Gibco) in 5% CO2 and 37°C conditions.

Bligh–Dyer lipid extraction and mild alkaline hydrolysis
To inhibit spt and de novo sphingolipid production, 5 µM myriocin (Cayman Chemicals) was added to bacterial culture. Forty-eight-hour liquid bacterial cultures in BHI supplemented with or without myriocin with OD > 1.00 were used for lipid extraction. Bacterial cultures were normalized to an optical density (OD 600nm) of 1–1.5 using sterile phosphate-buffered saline (PBS; Corning) with a total volume of 5 mL. Samples were then centrifuged at 4000×g for 20 min at 4°C. The supernatant was removed, and the bacterial pellets were resuspended in sterile PBS and transferred to 1.5 mL microtubes and washed twice using PBS. Lipids were then extracted using the Bligh–Dyer method50. Briefly, bacterial pellets were then resuspended in 100 µL sterile water and treated with 400 µL chilled 1:2 parts chloroform:methanol (Fisher Chemical). Samples were vortexed for 5 min at max speed using a vortex microtube adapter. Samples were then treated with 100 µL chilled chloroform and 100 µL sterile molecular-grade water (Corning) and vortexed for 1 min. Samples were then centrifuged at 2000×g for 10 min. The bottom layer containing chloroform and lipids was collected into a new tube and air-dried in a sterile environment. Samples were stored at −80°C until resuspended in 500 µL of cell culture media corresponding to the mammalian cells used for the assay. Reconstituted lipids were stored at −20°C. Mild alkaline hydrolysis for the removal of phospholipids from lipid extracts was performed as previously described in ref. 42. Briefly, the lipid extracts before drying were treated with 0.02 N sodium hydroxide (NaOH; Acros Organics) for 30 min at 37°C. Lipids were then dried, stored, and reconstituted as above.

Lipidomics analysis of Bacteroides xylanisolvens cultures
B. xylanisolvens cultures were grown in technical triplicates (n=3) for myriocin and no myriocin treatment. Bacterial cultures were pelleted by centrifugation, washed by resuspension and centrifugation in PBS, and total lipids were extracted using the Bligh–Dyer method50. Briefly, extracts were dried in a vacuum concentrator, then reconstituted in 2:1 acetonitrile-methanol. Lipids were analyzed by HILIC coupled with ion mobility-mass spectrometry (IM-MS). Chromatographic separations were carried out with a Phenomenex Kinetex HILIC column (50×2.1 mm, 1.7 µm) on a Waters Acquity FTN UPLC (Waters Corp., Milford, MA) with 95% acetonitrile/5% water/5 mM ammonium acetate as mobile phase A and 50% acetonitrile/50% water/5 mM ammonium acetate as mobile phase B51,102,103. Collisional cross-section calibration and IM-MS analysis were conducted on a Waters Synapt XS HDMS (Waters Corp.) in the positive ionization mode using a series of phosphatidylcholines, covering the range of m/z 454–980 and CCS 213–330 Å2, as the calibrants and negative ionization mode using a series of phosphatidylethanolamines, covering the range of m/z 410–786 and CCS 199–274 Å2, as the calibrants102–104. Data alignment and peak detection were performed in Progenesis QI (Nonlinear Dynamics; Waters Corp.) with normalization to all compounds. Retention time calibration and lipid identification were calculated with the Python package LiPydomics (version 1.6.8) with default settings [https://github.com/dylanhross/lipydomics]105.

THP-1 cell differentiation and treatment with bacterial lipids and JBK70
THP-1 cells were treated with 50 ng/mL PMA (phorbol 12-myristate 13-acetate) (Invitrogen) and plated into 24-well tissue-culture treated plates using 5×10^5 cells/well. Cells were incubated for 48 h prior to exposure to additional stimuli. Adherent cells were rinsed using sterile PBS, and the cell media was replaced. Bacterial lipids were applied to differentiated THP-1 cells in 1:100 concentrations for 24 h for pre-conditioning. Cells were rinsed with sterile PBS and then treated with heat-killed JBK70 (1:10) for a final dilution of 1:100. After 18 h, the supernatant of the cells was collected. For cell counts, adherent cells following lipid treatment were rinsed with sterile PBS and removed using Trypsin-EDTA (Corning). Cells were counted on a hemocytometer with Trypan-blue (Gibco).

Preparation of heat-killed JBK70
V. cholerae strain JBK70 was inoculated into LB broth from a glycerol stock and incubated at 37°C with agitation for 24 h. Opaque cultures were centrifuged and resuspended in an equal volume (4 mL) of cell culture media corresponding to the mammalian cells used for the assay. Samples were then diluted 1:10 in cell media and heated at 95°C for 30 min. Heat-treated samples were cooled to room temperature before application to mammalian cells. Aliquots of heat-killed JKB70 were plated on LB agar to ensure no live bacteria were present.

Enzyme-linked immunosorbent assays (ELISA) and cytokine measurements
Cell supernatants were analyzed using ELISA for cytokine production. ELISA kits were purchased from R&D Systems and used according to the manufacturer’s instructions, including IL-1β, IL-6, IL-10, TNF-α, MCP-1, and IL-8. Cytokine levels were determined using the 4-parameter logistic regression model calculated in GraphPad Prism (v10) according to the manufacturer’s instructions. For multiplex assays, samples were analyzed using the human cytokine proinflammatory focused 15-Plex Discovery Assay (HDF15) (Eve Technologies Corporation). Briefly, samples stored at −80°C were thawed and centrifuged at 3000×g for 10 min at 4°C. The supernatants of the samples were collected and shipped to Eve Technologies on dry ice for processing.

PBMC stimulation and treatment with bacterial lipids and JBK70
Frozen PBMCs were thawed and plated into 24-well tissue-culture treated plates using 1×10^6 cells/well in RPMI media containing 50 ng/mL PMA. Cells were incubated for 24 h. Adherent cells were rinsed with PBS and then used for lipid treatment and JBK70 stimulation as outlined above. Briefly, the adherent cells were rinsed, and bacterial lipids were applied for 24 h. Cells were then rinsed and treated with heat-killed JBK70 for a final dilution of 1:100. After 18 h, the supernatant was collected and tested for cytokines.

RNA extraction and sequencing
THP-1 cells were collected after exposure to V. cholerae and/or bacterial sphingolipids for RNA sequencing. After collecting the cell supernatant for cytokine analysis, adherent cells were treated with TRIzol (Invitrogen) and transferred into a 1.5 mL microtube. Chloroform was added, centrifuged at 4°C, and the aqueous layer was transferred into a new 1.5 mL tube and treated with 70% ethanol (Spectrum). Samples were then processed using the Qiagen RNeasy kit (Qiagen) according to the manufacturer’s instructions with on-column DNase treatment. RNA was eluted using 30 µL molecular-grade water and additionally treated with TURBO DNase (Ambion). RNA quality was assessed using TapeStation (Agilent). RNA libraries were prepared for sequencing using Illumina Stranded Total RNA Prep with Ligation with Ribo-Zero Plus according to the manufacturer’s instructions. Sequencing was performed using P3 Reagents (Illumina) and RNA UD Indexes Set A, Ligation (Illumina) on a NextSeq 2000 as 150-bp pair-ended reads.

RNAseq preprocessing and analysis
Fastq files were obtained, and RNA reads were curated and aligned using the RNA snakemake pipeline SEAsnake106. The raw data can be accessed through NCBI, BioProject PRJNA1170288. Briefly, RNA quality was assessed with FastQC, adapters were removed, and low-quality sequences were filtered using AdapterRemoval. Reads were aligned to the Homo sapiens genome GRCh38 (release 108) using STAR. Alignments were quality checked and filtered using samtools and Picard, and read counts were determined using Subread. RNA samples were then quality assessed and removed if libraries had <1,000,000 sequences, greater than 1 coefficient of variation (CV) coverage, and less than 75% alignment to the genome. No sample or library was removed using these criteria. Low-abundance genes with a minimum CPU of 0.1 and present in a minimum of 2 samples were filtered out. 4244 (21.22%) of 20001 genes were removed. Gene counts were then normalized by trimmed mean of M (TTM) and voom (log2 transformation of counts per million (CPM)). Analysis was then performed on the following comparisons: SL-containing with and without HKJBK70, SL-reduced with and without HKJBK70, SL-containing and SL-reduced without HKJBK70, and SL-containing and SL-reduced with HKJBK70. Linear modeling was performed using kmFit (Kimma software on R; [https://github.com/BIGslu/kimma]107 with weighted factors to obtain the estimate (log2 fold change) and FDR values. Genes were considered significant with FDR < 0.3 and log2 FC greater than 1 or less than -1. Enrichment pathways were determined using the R package BIGslu/SEARchways [https://github.com/BIGslu/SEARchways] with BIGprofiler that employs clusterProfiler using the gene set database and hallmark gene sets available in the Broad Molecular Signatures Database (MSigDB)108. Pathways were considered significant with FDR < 0.05.

Statistical analysis
Visualizations, figures, and statistical testing were generated using R (v4.0.2) and GraphPad Prism (v10) unless otherwise specified. Strain trees were generated using GToTree (v1.5.51)109 and visualized using the Interactive Tree of Life (v6.1.2)110. Shannon indices were compared between sexes using a two-tailed Mann–Whitney U-test. A Kruskal–Wallis test with Dunn’s adjustment for multiple comparisons was used to assess differences in Shannon index between age groups, as well as between responders and nonresponders for each vaccine response measure.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Weil and colleagues present an interesting premise that variations in response to the oral cholera vaccine is influenced by gut microbial sphingolipids. Authors present a set of in vitro experiments to suggest that sphingolipid-producing bacteria may directly contribute to vaccine response through the interaction of bacterial sphingolipids with immune cells. Major concerns with the manuscript are in the strength of the assertions that gut microbiome derived sphingolipids are involved in vaccine responsiveness through mediating innate immunity. In all the bacterial work, the attribution of response was based on lipid extracts derived from treating Bacteroides with a drug that inhibits de novo sphingolipid synthesis, yet no lipid profiling of these extracts was presented to determine the level of inhibition achieved. Moreover, if events are to be attributed to bacterial sphingolipids, then it is preferable that these sphingolipids are isolated, identified, and then used to assay the system. Other concerns include the strength of the connection of the responder phenotype to microbiome community composition and the premise for the focus on sphingolipids as a driver of the response. Overall, the responses to bacterial lipid extracts are interesting but more details on the lipid profiles and a stronger connection to the human cohort study would be needed to support a meaningful role of gut bacterial-derived sphingolipids on oral cholera vaccine response. Below are specific comments on the manuscript.\n1. How was it determined that the myriocin treated Bx cells were sphingolipid free? Would need to measure the lipids in the sphingolipid fractions in these experiments to make stronger statements.\n2. Differing levels of mainly host-derived sphingolipids in stool suggest that there may be other non-microbiome related intestinal track/lipid processing characteristics between the groups.\n3. Were the sphingolipid profiles measured in Bx present in the lipids in responder stool?\n4. What were the other lipids that were high in responders – maybe they are more responsible for phenotype\n5. Is there any cell death associated with lipid treatment in cell culture model?\nMinor:\n1. Line 490 Bacteroides serine palmitoyltransferase (spt) abundance was determined by querying BT_0780. Perhaps this is meant to be BT_0870. BT_0780 is an uncharacterized acetyltransferase.\n2. Figure 2C – color scale not showing up in provided version\n3. SPT gene abundance is almost a direct readout of presence of Bacterodies/Prevotella species – any insight into why this enrichment didn't come out in the community composition analysis?", "reviewer_2": "The manuscript entitled” Gut bacteria-derived sphingolipids alter innate immune responses to oral cholera vaccine antigens” is a follow up study from the parent mSphere article entitled “Cholera toxin and O-specific polysaccharide immune responses after oral cholera vaccination with Dukoral in different age groups of Bangladeshi participants”. Here the authors utilise the faecal materials collected at the time of vaccination and undertake metagenomic sequencing to infer microbiome associations to the cholera toxin vaccine responder and non-responders.\nThe authors subsequentially undertake a fundamental unit of analysis assessment of the metagenomic data. Using this approach, the authors find an association between Bacteroides, specifically xylanisolvens strains and Prevotella and memory B cells- which appears to be strain specific.\nNext the authors undertake a quantitative mass spectrometry-based lipidomics approach to the responder/non-responder stool samples. They find that several sphingolipid derivatives were more abundant is vaccine responders.\nSubsequently the authors undertake a series of experiments to validate these findings in human macrophages. They choose a particular Bacteroides xylanisolvens, isolated from one of the participants, for these studies.\nComments\nThe article is well written and easy to understand. It details the metagenomic sequencing of 98 participants in 3 defined cohorts between 2-5 years 5-17 years, and 18-46 years. The authors are clearly state the male/female ratios and the first 2 cohorts are equally balanced, with only the final cohort having a female bias.\nAs the parent article is published elsewhere, the authors are unfortunately unable to detail how the responder/non-responder ratios were defined. Instead the authors focus directly on the metagenomics data. Their findings indicate a possible difference in microbiomal diversity by sex; however, they do not find a significant difference by Shannon diversity or Bray-Curtis community structure between the vaccine responders and non-responders and hence.\nThe reviewer also has reservations regarding the THP-1 validation of these findings. This is the most important part of this manuscript, as it attempts to validate this important association. The authors take a series of experiments using a mammalian cell line, THP-1, when they should be using human cells derived from PMBCs. The initial assessment of innate responses includes ELISA validation of Il-1b, Il-6, Tnf-a and Il-8 (figure4). However, the authors do go on to expand this model (figure 5) and undertake a more thorough analyses using a 15-plex multi assessment of cytokines and RNA sequencing of resultant THP-1 cells. However, this assessment should be followed by macrophages derived from participant/control PMBCs and include responder material if possible and or in vivo validation in mouse models. Currently this association is not sufficiently validated\nMethods:\nAuthors need to detail fecal sample collection and timing information.\nAuthors do not detail the library preparation process at all.\nFigure 1: Inclusion of schematic figure detailing sample collection timepoints and participant cohorts would be a useful inclusion\nFigure 2:2c- Missing priority scoring colour information.\nFigure 3: Stats are not shown on the graph", "reviewer_3": "Chac et al. present a study investigating mechanisms by which gut microbiota may influence response to an oral cholera vaccine.\nThe first half of the study focusses on clustering de novo assembled bacterial genes into co-abundant gene groups (CAGs), which are subsequently correlated with several measures of humoral immune response to vaccination a cohort of 89 individuals.\nThe second half of the study focusses on exploring the impact of microbially-produced sphingolipids on innate immune responses through in vitro stimulation of THP-1 derived macrophages with heat-killed V. cholerae strain JBK7.\nWhile both of these halves present interesting and potentially important work, I don’t see strong justification for combining them into a single manuscript. In the first half of the manuscript, much emphasis is justifiably placed on de novo functional and strain-level characterization of the microbiome in relation to vaccine response. However, the only reason given for focussing on sphingolipids in the second half of the study is that genera identified as containing CAGs correlating with vaccine response (e.g. Bacteroides and Prevotella - lines 181-182) are known to produce sphingolipids.\nCritically, there seems to be no evidence that any of the CAGs identified as correlating with vaccine-induced IgA/IgG are associated with production of sphingolipids. There is also no effort to identify evidence of enrichment in sphingolipid production potential in strains identified as important using the humoral CAG-based strain identification method.\nWithout any bioinformatic evidence for sphingolipids relating to humoral vaccine response it is unclear a) why a metagenomic shotgun sequencing and complex CAG-based analysis approach is necessary in order to identify sphingolipid producing genera, and b) why the study subsequently chose to focus on sphingolipids rather than follow up on functions encoded in any of the 323 CAGs that were identified as correlating with vaccine response.\nUnfortunately, without clearer evidence to justify linking the two halves of this study I can find no rationale for publishing it as it’s currently presented. I recognise the potential importance of this study as well as the careful and innovative work that’s been done to date, and am including additional minor comments below in the hope that they are helpful towards a future publication.\nMINOR COMMENTS:\nLINE 396: Would be very useful to clarify if stool samples were collected on site, or whether they were collected at home and brought to the study centre? How long were they out of the freezer?\nLINE 90: Use of the word simple in section header is subjective.\nLINE 98: Serotypes (Ogawa and Inaba) are introduced here without prior explanation. A small amount of detail for unfamiliar readers would be very helpful.\nLINE 414: Given the emphasis on CAG-based strain identification (lines 148, 457), it would be useful to see more justification of the approach taken for taxonomic annotation. In particular that use of translated alignment and non-redundant gene clustering (90% sequence identity) are not limiting/skewing the accuracy of taxonomic identification. Why not just assign taxonomy to the original contigs generated by Megahit, thereby avoiding multiple unnecessary levels of nucleoitide-protein translation that are likely to reduce the accuracy of this important step?\nLINE 100: It's unclear whether the result reporting differences in IgA vs IgG MBC response are general, or referring to only those individuals where an increase in MBC response was observed (line 99).\nLINE 101: I can’t easily discern this result from the data presented (Fig 1A?) and more explanation would be helpful. i) why are the numbers of OSP-associated MBC responders and non-responders different between IgA vs IgG comparisons (62 vs 46, respetively)? ii) what evidence is there that the numbers developing an IgG directed OSP-specific MBC response (23) is greater than the number developing an IgA directed OSP-specific MBC response (19)?\nLINE 112: Given 22/89 study participants are aged 2-5, this is surprising and looks to be driven by a subset of individuals aged 20+ with low diversity scores. Could the authors not provide rarefaction/collectors curves to demonstrate all individuals were sequenced to an adequate depth to ensure robust diversity measures?\nLINE 450-453: A bit more detail would be helpful. How were read counts divided for non-redundant genes that aligned to multiple reference genomes? Were summed counts adjusted to account for variation in reference genome size?\nFIGURE S3: It would be very useful to see how, on these plots, the number of genomes aligned to relates to the number of genes in a CAG.\nLINE 153: Given the emphasis on strain-specificity, it would be useful to provide a little more detail/i) What is the working definition of a strain within the RefSeqq database used to map CAGs, ii) how does the strain identification method (Line 457) account for differences in variation in the number of reference genomes across taxonomic groups, iii) what evidence is there that differential alignment success of CAGs to genomes isn't simply an artefact of the variable quality of the underlying genome assemblies?\nLINE 229-233: I found this sentence difficult to understand as it appears to contrast SL-containing Bx supernatant to SL-containing lipid fractions, which doesn't accord with the referenced figure panel (Fig 5A).\nLINE 239: Incorrect figure reference?\nLINE 287: I disagree with the statement – here and elsewhere – that the approach is reference free. Unless I've misunderstood, strain annotation of CAGs is dependent on mapping them to genomes in the RefSeq database?"}, "author_response": "Overall Response (Reviewer 1):\nWe appreciate this critique provided to help us improve our manuscript. To address the reviewer concerns about the overall strength of assertions that bacterial-derived sphingolipids impact our model of vaccine responsiveness, we have performed additional experiments and added to our collaborative group Dr. Libin Xu (bacterial lipid expert from our Department of Medicinal Chemistry at the University of Washington) who has performed ion mobility-mass spectrometry to determine the lipids present in bacterial supernatant from our lipid producing bacteria. These additional experiments, in addition to other experimental additions and explanations below, have strengthened the assertions we present in this manuscript.\nResponses:\n1. How was it determined that the myriocin treated [Bacteroides xylanisolvens] Bx cells were sphingolipid free? Would need to measure the lipids in the sphingolipid fractions in these experiments to make stronger statements.\n● Response: We agree that this additional analysis will strengthen our manuscript. To address this, we measured the lipids in myriocin-treated B. xylanisolvens culture lipid extracts using hydrophilic interaction liquid chromatography (HILIC) followed by ion mobility-mass spectrometry (IM-MS). We partnered with Dr. Libin Xu, a medical chemist and expert in bacterial-derived lipids for this work and for the interpretation of the results. The methods used to extract lipids from Bx cells were identical to those used in our cell culture experiments. In HILIC-IM-MS, ceramides are the primary readout for estimating sphingolipid abundance. Ceramides are an essential structural backbone of all sphingolipids and are measured to approximate sphingolipid content. The use of myriocin inhibits spt (serine palmitoyltransferase) function, and this is the first step required for bacterial de novo sphingolipid synthesis from a ceramide.\nThe total ceramides we measured in myriocin-treated Bx lipids fractions were significantly lower than in non-myriocin treated lipid fractions. These results were added the manuscript as Supplemental Figure S7. In these results, we observed that the Bacteroides lipid fraction has reduced, but not zero, ceramides. Therefore, we have revised the text to refer to the myriocin-treated samples as sphingolipid-reduced or SL-reduced (rather than sphingolipid-free, as in our initial submission).\nFigure S7. Myriocin treatment of Bacteroides xylanisolvens culture lipid extract reduces ceramide levels.\nB. xylanisolvens (Bx) was grown in BHI supplemented media in anaerobic conditions with or without the spt inhibitor, myriocin. Lipids were then extracted using the Bligh-Dyer and measured using hydrophilic interaction liquid chromatography and ion mobility-mass spectrometry. Data is shown as the total intensity of detected positive (Cer+) and negative (Cer-) mode ceramides tested in triplicate. Intensity refers to the signal strength of an ion detected by the mass spectrometer and is proportional to the abundance of that ion. Ions are measured in either negative or positive mode based on their chemical structure and functional groups. Unpaired t tests with SEM are shown. ***, P≤0.001.\n2. Differing levels of mainly host-derived sphingolipids in stool suggest that there may be other non-microbiome related intestinal track/lipid processing characteristics between the groups.\n● Response: We agree and acknowledge that the sphingolipids found in human stool may be host derived, diet derived, or bacterial derived. Only some sphingolipid structures can be determined to be derived from bacteria based on branching structure. To our knowledge, it is currently not possible to determine the origin of some stool lipids, because Bacteroides species produce some similar sphingolipid structures to those derived from humans (this is described in An et al., 2010 - PMID: 20855611; Johnson et al., 2020 - PMID: 32424203) and because plant materials and other dietary sources may also contain non-human sphingolipids. We have expanded on our limitations to include that host-derived sphingolipids in stool, that are unrelated to the microbiota, may also be responsible for the findings we observed in human stool.\n3. Were the sphingolipid profiles measured in Bx present in the lipids in responder stool?\n● Response: To address this important comment, and further link our human data with our experimental results, we sought to compare the lipid profiles from study participant stool and the B. xylanisolvens-produced lipids. Due to limitations of branching and structural resolution using ion mobility-mass spectrometry and incomplete scientific knowledge about sphingolipid structures (especially those made by bacteria), our confidence in fully resolving structures of even some lipids contained in stool was extremely low. However, in our analyses of the stool and bacterial lipids, we identified two bacterial-derived sphingolipids detected with high confidence both in stool and in the bacterial lipid culture. These two identifiable ceramides were more abundant in vaccine responders and also found in B. xylanisolvens lipids – these are CER(16:0) and CER(18:1). We have revised the manuscript accordingly and added the supplementary data Table S9 that shows the abundances of these two lipids in extracts from Bacteroides xylanisolvens lipid fractions. Study participant stool measures of these two ceramides are also shown in the context of other lipids detected in human stool in Figure 3B, and quantities are listed in Table S8. To illustrate this point and the comparison, we have created a hybrid table here for this reviewer question.\nTable S9 and excerpt of Table S8, combined data: Ceramides detected in B. xylanisolvens lipid fractions with or without spt inhibitor myriocin and in study participant stool. The SL-reduced condition in bacterial lipidomics consist of B. xylanisolvens grown in culture with the spt inhibitor, myriocin. Lipids were the extracted from culture using the Bligh-Dyer method and measured by ion mobility-mass spectrometry. Fecal lipids from study participants were measured using the Sciex Lipidyzer platform described in the methods and as in PMID: 31832778.\n4. What were the other lipids that were high in responders – maybe they are more responsible for phenotype\n● Response: The other lipids that were differentially abundant in vaccine responders include diacylglycerol, cholesterol esters, and free fatty acids. The most abundant lipids and the association with responder and nonresponder status and abundances can be found in our new Table S8. We have also added comment on this to the manuscript and study limitations. These lipid species are primarily produced by the host and were not further investigated given our focus on the microbiome and microbiome-derived metabolites.\n5. Is there any cell death associated with lipid treatment in cell culture model?\n● Response: We assessed the cell numbers for 24 hours and did not observe any cell death associated with the lipid treatment in our THP-1 model. This data has been added to the manuscript as Figure S9.\nFigure S9. Lipid treatment did not alter THP-1 cell viability. The number of adherent cells were counted after THP-1 derived macrophages were incubated with Bacteroides xylanisolvens lipid extracts that was grown with or without myriocin, SL-reduced and SL-containing, respectively, for 24 hours. Statistical analysis was performed using one-way ANOVA. ns, P>0.05. Bars indicate mean with SEM.\nMinor:\n1. Line 490 Bacteroides serine palmitoyltransferase (spt) abundance was determined by querying BT_0780. Perhaps this is meant to be BT_0870. BT_0780 is an uncharacterized acetyltransferase.\n● Response: This is a typo, and the correct gene is BT_0870; this has been corrected. Thank you for this important correction.\n2. Figure 2C – color scale not showing up in provided version\n● Response: The color scale for Figure 2C has been revised.\n3. SPT gene abundance is almost a direct readout of presence of Bacteroides/Prevotella species – any insight into why this enrichment didn't come out in the community composition analysis?\n● Response: The primary reason that these genera did not emerge from our computational analysis is because not all the species in the Bacteroides and Prevotella genera have the spt gene. Spt is strain specific and is also not consistent across one species within these genera. To exemplify this, in Figure 2d we show a phylogenetic tree of Bacteroides xylanisolvens strains found in our study population, and only 7 of the 21 Bacteroides xylanisolvens strains contained spt. Therefore, due to the abundance of spt negative Bacteroides strains, our study was not powered enough to detect this difference at the genus or species level. In particular, the Prevotella genus is appreciated to encompass an extremely wide genetic and functional diversity (reviewed in PMID: 34050328). Therefore, while we did not examine this genus as closely in this manuscript, the Prevotella detected in our study population are also very likely to have a mix of spt+ and some spt- strains, similar to our findings in Bacteroides.\n\nOverall Response (Reviewer 2):\nWe appreciate this reviewer’s evaluation of our manuscript. First, we have clarified in the manuscript the immunologic measures we used and the definitions of vaccine recipients as responders or nonresponders. Vaccine recipients were grouped according to development of cholera toxin- and O-specific polysaccharide-specific memory B cell (MBC) response measures in peripheral blood, developed between the pre- and post-vaccination timepoints. Both IgA and IgG-specific responses for each of these antigens were measured in each participant. Individuals with increased MBC counts following vaccination were defined as vaccine responders (R) while those that had unchanged or decreased MBC response after vaccination were defined as vaccine nonresponders (NR). We have revised the manuscript with a clearer description of this information to address this reviewer comment.\nWe agree that validation of our in vitro human macrophage experiments is needed, and we wanted to address this using more diverse human immune cells. We used peripheral blood mononuclear cells (PBMCs) donated by healthy adults. We repeated our THP-1 experiment using PBMCs; specifically, we used pre-treatment of cells with B. xylanisolvens lipids and subsequent heat-killed JBK70 stimulation. We then measure innate cytokine responses and found increased IL-6 across PBMC donors between the SL-containing lipids condition compared to SL-reduced lipids. Although we did find donor-to-donor variability in the results, we overall found that PBMC experiments across four donors validated our results from human macrophages and we have added this data to the manuscript as a supplemental figure: Figure S14.\nFigure S14. Stimulated peripheral blood mononuclear cells (PBMCs) have increased inflammatory cytokine responses to JBK70 following SL-containing lipid pretreatment.\nFold change in cytokine response in supernatant from PBMCs after preconditioning with B. xylanisolvens lipid extracts and treated with heat-killed JBK70. Bars represent mean with SEM. Unpaired T-test comparing fold change: *, P<0.05; **, P<0.01; ***, P<0.001. Each number represents a different PBMC donor. #1 = 31 yr female; #2 = 25 yr male; #3 = 36 yr female; and #4 = 35 yr male. Three to four technical replicates were performed on each experimental condition depending on the PBMCs available after thawing.\nFor the validation experiment, we used human-derived models rather than mouse models because V. cholerae infection is human-restricted and there is not a suitable model for modeling oral cholera vaccine responsiveness. Studies have been done using mice, but these often use artificial routes of antigen delivery and poorly predict human immune responses. Mouse models for immune response to oral vaccines also do not model the human microbiota, and precise experimental manipulation of the gut environment and in situ analysis of host-microbe interactions is very challenging in animals. For these reasons, we pursued validation studies in human-derived models rather than in animals. Unfortunately, additional samples from study participants such as peripheral blood or responder material (fecal samples) were not available for further validation assays.\nResponse to questions and minor comments (Reviewer 2):\nMethods:\nAuthors need to detail fecal sample collection and timing information.\n● Response: The collections and timing of fecal samples have been revised regarding sample collection and timing, and this is listed in the “Study design, subject enrollment, and sample collection” section of the methods, and shown in a revised Figure 1a.\nAuthors do not detail the library preparation process at all.\n● Response: We have added details to the library preparation in the revised methods section of the manuscript.\nFigures:\nFigure 1: Inclusion of schematic figure detailing sample collection timepoints and participant cohorts would be a useful inclusion\n● Response: A schematic was added as Figure 1a to illustrate and clarify sample collection timepoints, the vaccination schedule.\nExcerpt of Figure 1a. Schematic of study design, OCV dosing, and sample collection timepoints, adapted from parent study. Responders were defined as an OCV recipient who developed an increase in MBC response between pre- and post-vaccination timepoints, and nonresponders were defined as participants with unchanged or decreased MBC responses between timepoints.\nFigure 2:2c- Missing priority scoring colour information.\n● Response: The color scale for Figure 2C has been corrected in the revised manuscript.\nFigure 3: Stats are not shown on the graph\n● Response: The figure has been corrected in the revised manuscript.\n\nOverall Response (Reviewer 3):\nThank you for this assessment of the structure of our study and description of a disconnect between two “halves” of work. Our overall goal in presenting this work as a whole and to illustrate a full computational to experimental testing pipeline, a process that starts with human samples and progresses to exploring mechanisms using in vitro models. Thus, we acknowledge this paper contains a range of methods and results, and we feel there is inherent value in presenting translational science such as this in a single manuscript. To achieve a smoother and more linear presentation of this process, we have reorganized and significantly revised the manuscript to clarify the sequential steps in this work that led to each next phase. With these changes, we feel this version more accurately reflects the series of steps we undertook as the study progressed, and this unites the two “halves” by demonstrating more clearly our rationale. We have added data to illustrate the steps that led to choosing specific in vitro experiments, based on the computational analysis, and we feel this has helped to unite the two “halves” of the manuscript.\nSpecifically, we have added:\n● Serine palmitoyltransferase gene(spt) counts from the stool of Responder and Nonresponder study participants. This is the enzyme required for bacteria to synthesize sphingolipids. Spt was more abundant in the stool of vaccine responders (Figure 3a), and this finding then resulted in our selection of a spt-positive B. xylanisolvens strain for experimental testing. We have reordered the description of results to demonstrate that the specific genera highlighted by our metagenomic results (sphingolipid-producing microbes) led us to next quantify spt in study participant stool.\n● We examined lipid content and abundance in study participant stool and compared this to the culture supernatant of lipids produced by Bacteroides xylanisolvens, and found specific ceramide groups (the head group required in bacterial sphingolipid synthesis) that were differentially abundant in these two groups, described in Table S8 and new Table S9 (link between these areas is further detailed in the revised manuscript and in the explanation to point #3 raised by Reviewer #1).\n● Revisions to the manuscript to reflect that the CAG analysis did not single out in the Bacteroides and Prevotella genera as associated with vaccine responsiveness because not all the species in the Bacteroides and Prevotella genera have the spt gene. Spt is strain specific and is thus not consistent across one species within these genera. To exemplify this, we show in Figure 2D a phylogenetic tree of Bacteroides xylanisolvens strains found in our study population, and only 7 of the 21 Bacteroides xylanisolvens strains contained spt. Therefore, due to the abundance of spt negative Bacteroides strains, our study was not powered enough to detect this difference at the genus or species level.\n● A clarification that CAGs are identified in this study for the purpose of mapping to genomes for the identification of specific strains for experimental testing. It would have been ideal if CAGs could be considered as functional units that we could query experimentally, but, unfortunately, a query of functions encoded by a CAG is not possible because all CAGs contain multiple genes, most of which are unannotated or nonspecific (Table S2).\nWe have revised the manuscript on all of the above points, and appreciate this reviewer’s comments about the disjointedness of our initial draft. This review has helped us to significantly improve the manuscript, in our opinion.\nResponse to questions and comments (Reviewer 3):\nMINOR COMMENTS:\nLINE 396: Would be very useful to clarify if stool samples were collected on site, or whether they were collected at home and brought to the study centre? How long were they out of the freezer?\n● Response: The collections and timing of fecal samples have been reworded for clarity and moved to the “Study design, subject enrollment, and sample collection” in the Methods. Stool samples were collected at the International Centre for Diarrheal Disease Research, Bangladesh, by study staff on Day 0 of the study (the day of first OCV dose) and immediately and frozen at -80C, and then shipped on dry ice to US investigators for storage, DNA extraction, and fecal lipidomics. We have added a schematic as Figure 1a to further demonstrate our study design in the revised manuscript.\nLINE 90: Use of the word simple in section header is subjective.\n● Response: This section header has been revised in the new draft of the manuscript.\nLINE 98: Serotypes (Ogawa and Inaba) are introduced here without prior explanation. A small amount of detail for unfamiliar readers would be very helpful.\n● Response: We have revised the manuscript introduction to include a description of the Ogawa and Inaba V. cholerae O1 serotypes.\nLINE 414: Given the emphasis on CAG-based strain identification (lines 148, 457), it would be useful to see more justification of the approach taken for taxonomic annotation. In particular that use of translated alignment and non-redundant gene clustering (90% sequence identity) are not limiting/skewing the accuracy of taxonomic identification. Why not just assign taxonomy to the original contigs generated by Megahit, thereby avoiding multiple unnecessary levels of nucleoitide-protein translation that are likely to reduce the accuracy of this important step?\n● Response: We have two primary motivations for using protein-level sequences instead of nucleotide-level contigs. (1) Using a comparison contig-level abundances across all samples would require some manner of co-assembly which can be prohibitively computationally intensive at this scale. (2) Functionally identical biological functions may be performed by genes with high amino acid sequence identity and secondary structure despite their nucleotide sequences being dissimilar enough to prevent assembly as a single contig.\nLINE 100: It's unclear whether the result reporting differences in IgA vs IgG MBC response are general, or referring to only those individuals where an increase in MBC response was observed (line 99).\n● We have revised the text to make clear that the summary data presented for vaccine responders refers only to study participants with increases in MBC responses after vaccination (vaccine responders).\nLINE 101: I can’t easily discern this result from the data presented (Fig 1A?) and more explanation would be helpful. i) why are the numbers of OSP-associated MBC responders and non-responders different between IgA vs IgG comparisons (62 vs 46, respetively)? ii) what evidence is there that the numbers developing an IgG directed OSP-specific MBC response (23) is greater than the number developing an IgA directed OSP-specific MBC response (19)?\n● Response: (i) The denominator of participants with OSP-associated MBC IgA responses differ from those with IgG responses due to the incomplete data (missingness) for some participants in the IgG MBC measurement. (ii) For each participant, IgG and IgA antibody secreting cells were assessed separately using ELISPOT assays (these results were published in PMID: 38391226), and thus, the two groups differ. There were more individuals who had an increase in OSP-specific memory B cells IgG responses (23/46 or 50%) compared to OSP-IgA MBC responses (19/62 or 30.6%).\nLINE 112: Given 22/89 study participants are aged 2-5, this is surprising and looks to be driven by a subset of individuals aged 20+ with low diversity scores. Could the authors not provide rarefaction/collectors curves to demonstrate all individuals were sequenced to an adequate depth to ensure robust diversity measures?\n● Response: We have added a rarefaction curve to demonstrate that all samples yielded at least 600,000 reads and appear to capture the majority of observed species. This information has been added to the manuscript as Supplemental Figure S3.\nFigure S3. Rarefaction curve of microbiome data. Each line represents an individual sample. Rarefaction performed on species-level identification of microbiome species.\nLINE 450-453: A bit more detail would be helpful. How were read counts divided for non-redundant genes that aligned to multiple reference genomes? Were summed counts adjusted to account for variation in reference genome size?\n● Response: We have added detail to clarify the justification for the approach used with alignment of non-redundant genes to multiple reference genomes, in which the biological expectation is that multiple representatives from a species or genus may contain homologous genes (see revised Methods section “CAG alignment to bacterial genomes”). While the summed counts were not adjusted for genome size, this would not have any impact on our statistical analysis which used the marginal association of each genome in isolation, and the genome size adjustment factor would have been applied equally to the abundance values across all samples for the same genome.\nFIGURE S3: It would be very useful to see how, on these plots, the number of genomes aligned to relates to the number of genes in a CAG.\n● Response: Figure S4 has been revised to illustrate the relationship between the number of genes in significant CAGs with the number of genomes that were aligned to that CAG.\nFigure S4. Number of NCBI reference bacterial genomes mapped to CAGs that were found to be associated with vaccine responsiveness. Only significant CAGs associated with vaccine response are shown, and each is represented as a circle or X. Significant CAGs are defined as a CAG containing ≥2 genes with a q value of ≤0.1 for an association with the listed vaccine response measure. Red x’s indicate CAGs that had no genomes aligned. The x-axis indicates the number of genes within each CAG while the y-axis represents the number of reference genomes to which genes in those CAGs align to with >90% sequence identity. Significant CAGs aligning to only one bacterial genome are not shown (20 CAGs for OSP IgA, 3 CAGs for OSP IgG, 12 CAGs for CT IgA, 5 CAGs for CT IgG).\nLINE 153: Given the emphasis on strain-specificity, it would be useful to provide a little more detail/i) What is the working definition of a strain within the RefSeqq database used to map CAGs, ii) how does the strain identification method (Line 457) account for differences in variation in the number of reference genomes across taxonomic groups, iii) what evidence is there that differential alignment success of CAGs to genomes isn't simply an artefact of the variable quality of the underlying genome assemblies?\n● Response: (1) The RefSeq database does not present a formal definition of a strain, and the term is being used here as a reference to a single reference genome which is present in the RefSeq database. We have added this definition in the methods for clarification in the revised manuscript. (2) The strain identification method is not expected to correct for differences in reference database composition, and is simply referring to the detection of genomic sequences which are homologous to the reference genome for a particular strain. (3) We have revised the limitations described in the manuscript to make clear that genomes with lower assembly quality will be more difficult to detect using this analysis approach.\nLINE 229-233: I found this sentence difficult to understand as it appears to contrast SL-containing Bx supernatant to SL-containing lipid fractions, which doesn't accord with the referenced figure panel (Fig 5A).\n● Response: The sentence has been revised for clarity in the manuscript.\nLINE 239: Incorrect figure reference?\n● Response: Thank you for this correction, and the figure references have been corrected in the revised manuscript.\nLINE 287: I disagree with the statement – here and elsewhere – that the approach is reference free. Unless I've misunderstood, strain annotation of CAGs is dependent on mapping them to genomes in the RefSeq database?\n● Response: The reference-free aspect of our study is the portion of the analysis in which CAGs are generated, quantification of the relative abundance of CAGs across metagenomes, and identifying CAGs which are differentially abundant across metagenomes as a function of vaccine response. The result of that portion of the analysis is a collection of organisms (defined by the presence of a particular collection of genes encoded in their genome) which are associated by relative abundance with vaccine response. The subsequent step of alignment to the RefSeq genome database is reference-dependent, but it is independent from and ancillary to the CAG analysis itself. We have revised this in the manuscript to clarify what we are referring to as “reference-free”."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have addressed my comments on those of the other reviewers in depth. I have no further comments to make."}, "author_response": ""}]

Decision letter