To elucidate the molecular mechanisms involved in the grazing resistance of V. cholerae to T. pyriformis, we applied MicroSPLiT (Kuchina et al., 2021), a single-cell transcriptomic approach for bacteria. We co-incubated V. cholerae and T. pyriformis in artificial seawater (ASW) with no addition of nutrients. A total of 5,344 bacteria were individually barcoded, with 222 cells identified under control conditions [ungrazed, starved cells incubated in ASW at room temperature (RT) for 6 h], and 5,122 cells identified in the experimental condition (grazed, cells incubated in ASW in the presence of T. pyriformis at RT for 6 h; Figure 1a).

Using the standard workflow of the R package seurat on the single-cell expression matrices, we resolved unique transcriptional states by clustering cells with a similar gene expression pattern, obtaining a total of 11 (0–10) clusters of cells with a distinctive gene expression profile (Figure 1b). The median number of Unique Molecular Identifiers (UMI) that are directly related to the number of mRNA transcripts per cell was low (Figure 1c) compared to other reports of bacterial single cell transcriptomics (average median for UMI/cell is >200; Thurman et al., 2010). Low UMI may be due to poor cell permeabilization. Therefore, permeabilization of the V. cholerae cells used to prepare the single-cell RNA library was confirmed using Ovalbumin Alexa Fluor 488 conjugate (Thermofisher, Waltham, Massachusetts, USA), resulting in the accumulation of green-fluorescent signal inside the cells. To further confirm that the low UMI/cell observed is a limitation of the experiment and not a barcoding inefficiency of cellular RNA, we compared the number of UMIs observed among the different clusters. The median UMI between clusters was significantly different (Figure 1d), indicating that even under low-nutrient conditions, where the overall number of UMI counts is low, the observed variation reflects genuine differences in metabolic activity across the cell clusters.

As observed in previous single-cell RNA-seq reports in bacteria, clusters of gene expression can be related to specific biological processes. For example, cluster 0 represents 33.54% (1718/5122) of all grazed cells in this experiment and reveals a slight upregulation of ompU (log₂FC = 0.58, adjustedPvalue = <1×10⁻⁶), a gene that we have previously found to be important for EFV production (Espinoza-Vergara et al., 2019). Cluster 7 reveals upregulation of flagellar genes (flhA, flhF, flaD, and flaC, log₂FC > 5.5, adjusted p-value = <1×10⁻⁶) and represents 3.68% (189/5122) of the total grazed population. The presence of specific subpopulations of cells overexpressing flagellar genes has been observed in previous barcoding-based single-cell RNA-seq reports (McNulty et al., 2023).

Genes with significant differential expression were identified among the different clusters (Figure 2; Supplementary File S1). Clusters 1, 5, 6, and 9 show upregulation of genes related to tRNAs and ribosomal proteins. In addition to the upregulation of ompU, cluster 0 shows upregulation of genes involved in the synthesis of the F₁F₀ -ATPase (VC2764, VC2766, and VC2770, avg. log₂FC > 1.4, adjusted p-value = <1×10⁻⁶), which is related to ATP synthesis (von Ballmoos et al., 2009) and for survival of acidic conditions by pathogenic bacteria (Sun, 2016).

Cluster 8 cells displayed upregulation of the transcriptional factor icsR (VC0747, log₂FC = 5.5, adjusted p-value = <1×10⁻⁶), a global transcriptional regulator which is involved in iron homeostasis and stress response (Lim and Choi, 2014), and the cysteine desulfurase gene (VC0748, log₂FC = 5.7, adjusted p-value = <1×10⁻⁶), which is important for mobilization of sulfur for iron–sulfur clusters involved in iron homeostasis (Zheng and Dean, 1994). VC2251 and VC2252, ompH, which is involved in the chaperoning of outer membrane proteins, and the outer membrane protein assembly factor BamA (Mathieu-Denoncourt and Duperthuy, 2025), respectively, are also upregulated (log₂FC > 5.9, adjusted p-value = <1×10⁻⁶). Cluster 10 shows strong upregulation of craA (VCA0849, log₂FC = 7.2, adjusted p-value = <1×10⁻⁶), a gene positively regulated by c-di-GMP, suggesting that this cluster identifies cells producing c-di-GMP, a key intracellular signal in V. cholerae involved in the positive regulation of biofilm genes (Beck et al., 2021). Cyclic AMP receptor protein (CRP, VC2614) is also shown to be upregulated in this cluster (log₂FC = 2.8, adjusted p-value = 2.9×10⁻²), and it has previously been demonstrated that HapR and CRP can bind to shared DNA sites, while also maintaining the ability to directly interact with each other, primarily in the blocking of the RNA polymerase binding site of CRP by HapR (Walker et al., 2023).

Interestingly, one of the most differentially expressed genes in the various clusters is a hypothetical 65-amino acid protein, VC0713. We observed that the expression of this peptide positively correlates with the expression of ompU in the different clusters, suggesting a potential connection in the expression of these two genes.

The expression of virulence factors has been demonstrated to be heterogeneous in bacterial populations (Davis, 2020; Gonzalez and Mavridou, 2019). In these experiments, we observed that there is differential expression of toxins that are regulated by the master regulator of quorum sensing in V. cholerae, HapR. Upregulation of the motility-associated killing (mak) factor genes makA, makB, and makC (log₂FC > 1.2, adjusted p-value < 2×10⁻¹⁸), along with the hemagglutinin (HA)/protease gene (VCA0865), hapA (log₂FC = 1.47, adjusted p-value = 1.1×10⁻¹³), was observed under grazing pressure (Figure 3). Interestingly, clusters 2 and 3 that show the strongest upregulation of the toxins, also show upregulation of VC1590 (acetolactate synthase) and VC1591 (oxidoreductase) (log₂FC > 1.5, adjusted p-values < 2×10⁻¹⁵), which are involved in the production of 2,3 butanediol. It has recently been shown in Clostridium perfringens that acetate is involved in the regulation of the secreted pore-forming toxin NetB (McNulty et al., 2023). This could suggest that in V. cholerae, the production of Mak toxins could also be regulated by a metabolite-toxin-regulated system, as was shown in C. perfringens.

To investigate and validate the role of Mak toxins in the grazing defense of V. cholerae against T. pyriformis, we constructed an in-frame ΔmakA, ΔmakB, and ΔmakE deletion mutant. After 24-h incubation, the makA and makE mutants show a survival defect (16 and 31% survival compared to WT, respectively) in the EFVs. However, the survival is restored to WT levels after complementing the makA gene in trans in the ΔmakA strain (Figure 4a). In the makE mutant, complementation resulted in enhanced EFV survival when compared to the WT (Figure 4c). Additionally, it was observed that deletion and complementation of the makB mutant had no effect on bacterial survival (Figure 4b). Furthermore, after approximately 5 h of co-incubation, we observed that the WT strain kills a fraction of the T. pyriformis population (5.07%), and this killing effect is not observed in the ΔmakA, ΔmakB, or ΔmakE deletions. Complementation of these three genes in trans fully restores the killing effect in the ΔmakA, ΔmakB, and ΔmakE strains (4.87, 3.60, and 3.73%, respectively) (Figure 5; Supplementary Figures 1–2; Supplementary Videos 1–3). Additionally, ΔhapR and ΔflaA mutants were tested, and the mak mutants were unable to kill the protozoa. HapR is the regulator of the mak operon (Dongre et al., 2018), while FlaA is required for flagellum formation, the site through which the Mak tripartite toxin is secreted (Nadeem et al., 2021a). Together, these data suggest that makA and makE in V. cholerae are involved in intracellular survival, and all three parts of the mak toxin are required for the killing of T. pyriformis.

To investigate the potential importance of other genes identified in the single-cell transcriptomic data that are involved in bacterial survival, mutants were generated and incubated for 24 h. After 24 h, the only mutant that was shown to have a differential survival rate compared to the WT was ΔnlpD (Supplementary Figure S3).

Organisms used in this study are listed in Supplementary Table S1. Bacterial strains were routinely grown in lysogeny broth (LB) and on LB agar plates. V. cholerae mutants were constructed by splicing by overlap extension PCR (Horton et al., 1989), and natural transformation was performed on chitin flakes (Dalia et al., 2014). Complementation was done using the expression vector pBAD24. Bacteria carrying the expression vector were grown in LB broth at 37 °C containing ampicillin 100 μg mL⁻¹ and 0.2% arabinose for gene expression.

Tetrahymena pyriformis was routinely passaged in 10 mL of growth medium containing peptone-yeast-glucose (PYG) (20 g l⁻¹ proteose peptone, 1 g l⁻¹ yeast extract) supplemented with 0.1 × M9 minimal medium (6 g l⁻¹ NaH₂PO₄, 3 g l⁻¹ K₂PO₄, 0.5 g l⁻¹ NaCl, 1 g l⁻¹ NH₄Cl) and 0.1-M sterile-filtered glucose in 25 cm² tissue culture flasks with ventilated caps (Sarstedt Inc., Nümbrecht, Germany) and incubated statically at room temperature (RT).

Prior to experiments, 250 μL of T. pyriformis were passaged in 10 mL of ASW medium (5.6 g l⁻¹ NaCl, 0.470 g l⁻¹ Na₂SO₄, 0.026 g l⁻¹ NaHCO₃, 0.08 g l⁻¹ KCl, 0.013 g l⁻¹ KBr, 0.600 g l⁻¹ MgCl₂.6H₂O, 0.130 g l⁻¹ CaCl₂.2H₂O, 0.002 g l⁻¹ SrCl₂.6H₂O and 0.002 g l⁻¹ H₃BO₃) supplemented with 1% heat-killed Pseudomonas aeruginosa PAO1 (HKB) in a 25 cm² tissue culture flask, and further incubated at RT statically for 24 h before enumeration and use. This process is necessary to remove the nutrient media and to acclimatize the ciliate to phagotrophic feeding (Noorian et al., 2017; Espinoza-Vergara et al., 2019).

To prepare heat-killed bacteria (HKB), P. aeruginosa was grown overnight in LB at 37 °C with shaking at 200 rpm and adjusted to (OD₆₀₀ = 1.0; 10⁹ cells mL⁻¹) in ASW. The tubes were then transferred to a water bath at 65 °C for 2 h and then tested for viability by plating on LB agar plates at 37 °C for 2 days. HKB stocks were stored at −20 °C.

V. cholerae A1552 was co-incubated with T. pyriformis in ASW. Briefly, T. pyriformis were enumerated by microscopy and adjusted to 10³ cells mL⁻¹ and added to co-cultures of V. cholerae A1552 adjusted to 10⁸ cells mL⁻¹ in ASW using a spectrophotometer (OD₆₀₀ nm). Co-incubation was performed for 5 h at RT, 60 rpm. After incubation, TritonX-100 was added at a final concentration of 1% to release intracellular and EFV-encased bacteria. The cells were centrifuged for 5 min (3,220 × g) at 4 °C, and the cell pellet was resuspended in a freshly made 4% ice-cold paraformaldehyde solution. Fixed cells were incubated overnight at 4 °C.

Fixed cells were centrifuged (5,500 × g) for 10 min. All centrifugation steps were performed at 5500 × g, 4 °C, 10 min. After centrifugation, cells were washed with 1 mL of 0.1 M Tris–HCl pH7 (Sigma) and centrifuged again. Cells were then resuspended in 250 μL of cold-PBS (pH 7.4, Sigma) supplemented with RNAase inhibitors (Enzymatics RI and Superase-IN) as indicated by the MicroSPLiT protocol (Kuchina et al., 2021). Following resuspension, 250 μL of 100% cold ethanol was added and incubated for 1 min at RT. Cells were again centrifuged and resuspended in a solution containing 0.04% Tween20 and 50 mM EDTA. The suspension was incubated for 3 min on ice. Following incubation, 750 μL of cold PBS supplemented with RNAase inhibitors was added. After centrifugation, cells were resuspended in 100 μL of a lysozyme buffer (0.1 M Tris pH7, 50 mM EDTA, 1 mg mL⁻¹ lysozyme) and incubated at RT for 15 min. Permeabilization was assessed by mixing 50 μL of permeabilized cells with Ovalbumin Alexa Fluor 488 conjugate (Thermofisher, Waltham, Massachusetts, USA), following the manufacturer’s recommendations.

In situ polyA tailing, reverse transcription, ligation barcoding, and library preparation were performed following the steps published in the MicroSPLiT protocol (Kuchina et al., 2021).

Single-cell gene expression matrices were obtained by mapping the cDNA reads (Read 1) to the reference genome of V. cholerae O1 El Tor strain N16961 using the STAR aligner (v2.7.10b) with the addition of several STARsolo parameters (−-genomeDir, −-readFilesIn, −-soloType, −-soloCBposition, −-soloUMIposition, −-soloCBwhitelist, −-soloCBmatchWLtype, −-soloUMIdedup, −-soloFeatures, −-soloMultiMappers) (Dobin et al., 2013). The resulting matrices were analyzed using the standard workflow of the R package seurat (v4.3.0) with default parameters except otherwise indicated (Butler et al., 2018). The data were log-transformed using the ‘NormalizeData(normalization.method = ‘LogNormalize’, scale.factor = 10,000)’ function. Five hundred most variable genes were selected using ‘FindVariableFeatures(selection.method = ‘vst’, nfeatures = 500)’ and z-scored using ‘ScaleData’. Linear dimensional reduction was performed using principal component analysis (PCA). Louvain clustering algorithm was used to generate clusters using the function ‘FindNeighbors(dims = 1:10)’ and ‘FindClusters(resolution = 0.5)’. Uniform manifold approximation and projection (UMAP) non-linear dimensional reduction techniques were used to visualize and explore the datasets using ‘RunUMAP(dims = 1:10)’ function. Differential markers in each cluster were identified using the ‘FindAllMarkers(logfc.threshold = 0.2)’ function. The average expression of the top 10 differentially expressed markers in each cluster was used to generate a heatmap using the ‘DoHeatMap’ function. The FeaturePlot function was used to visualize specific marker expression in the clustered population.

To assess the number of dead T. pyriformis, co-incubations were performed as described before. Briefly, V. cholerae and T. pyriformis were co-incubated in ASW for 5 h at RT. After approximately 5 h, dead cells accumulated at the bottom of the well and were counted using an inverted epifluorescence microscope (Nikon Eclipse Ti inverted microscope).

To produce EFVs, V. cholerae strains, ΔlacZ (WT) and mutants, were co-incubated with T. pyriformis independently in ASW. Briefly, T. pyriformis were enumerated by microscopy and adjusted to 10³ cells mL⁻¹ and added to co-cultures of V. cholerae, adjusted to 10⁸ cells mL⁻¹ in ASW using a spectrophotometer (OD₆₀₀ nm). After overnight incubation at RT, samples were analyzed using an inverted epifluorescence microscope (Nikon Eclipse Ti inverted microscope) to detect the presence of EFVs in the supernatant. To purify V. cholerae-EFVs, supernatants were filtered (by gravity) several times through 8-μm filters (Millipore, Darmstadt, Germany), and the filters containing EFVs were suspended in 1 mL ASW. The EFVs were incubated for 1 h with gentamicin 300 μg mL⁻¹ at RT to kill any remaining extracellular bacteria. After gentamicin treatment, V. cholerae-EFVs pellets were collected by centrifugation (3,220 × g for 20 min), washed three times in ASW, and suspended in 1 mL of ASW. Finally, Triton X-100 was added to a final concentration of 1% to lyse the EFVs and release bacterial cells. The lysed EFV aliquots (WT and mutant) were mixed 1:1 and plated on LB supplemented with X-Gal (60 μg mL⁻¹) to differentiate between the V. cholerae ΔlacZ (white colonies) and the mutant strains (blue colonies).

Vibrio cholerae WT, ΔnlpD, ΔVC1590, ΔVCA1097, and ΔhapA strains were co-incubated with T. pyriformis independently in ASW. T. pyriformis were enumerated by microscopy and adjusted to 10³ cells mL⁻¹ and added to co-cultures of V. cholerae adjusted to 10⁸ cells mL⁻¹ in ASW using a spectrophotometer (OD₆₀₀ nm). After 24 h, samples were treated with Triton X-100 to a final concentration of 1% to lyse the EFVs and release bacterial cells from the protozoa. These aliquots were plated onto LB and enumerated for survival.

Statistical analysis was performed using GraphPad Prism version 8.4.3 for Windows, GraphPad Software, La Jolla, California, United States.¹ Data that did not follow a Gaussian distribution was determined by analyzing the frequency distribution through the Shapiro–Wilk normality test. Two-tailed Student’s t-tests were used to compare means between experimental samples and controls using Welch’s t-test correction. For experiments including multiple samples, one-way ANOVA was used for the analysis, the Kruskal-Wallis test was used to compare the medians of non-normally distributed data, and Tukey’s test was used to compare the averages of normally distributed data.