Introduction

Biofilms, contributing to approximately 80% of chronic and recurrent microbial infections in the human body¹, are complex microbial ecosystems characterized by a diverse array of bacterial cells existing in various physiological states^1–3. This heterogeneity within biofilms is influenced by multiple factors, including bacterial interactions (such as competition, symbiosis, and parasitism) that contribute to stable community structures^4,5, environmental factors (like host environment, aquatic conditions, and nutrient concentrations), and spatial organization. Bacteria occupying different positions within the biofilm perform distinct roles⁶: some mediate material exchange between cells and with the extracellular matrix, others facilitate complex communication systems between bacteria and with the host through signaling molecules, while certain bacteria participate in various energy conversion processes. This intricate division of labor not only contributes to bacterial heterogeneity but also enhances the biofilm’s overall resistance to various stresses⁷. The resulting cellular and functional diversity reflects the complex nature of the biofilm ecosystem, allowing it to adapt and persist in challenging environments. However, the study of biofilms faces significant limitations, primarily stemming from challenges in investigating heterogeneity within a bacterial population^8,9. Single-cell RNA-seq emerges as a promising avenue for addressing this heterogeneity^10–18.Expending on established protocols for cell fixation and permeabilization which facilitate in-cell barcoding while avoiding cell lysis, combinatorial barcoding-based bacterial scRNA-seq techniques, such as Prokaryotic Expression profiling by Tagging RNA in situ and sequencing (PETRI-seq)¹² and Microbial Split-Pool Ligation Transcriptomics (microSPLiT)¹⁴, have been developed. Nevertheless, these methods encounter challenges in terms of low transcript recovery rates due to overwhelmingly abundant rRNA, restricting the comprehensive analysis of within-population heterogeneity. In comparison to mammalian cells¹⁹, the absence of mRNA polyadenylated tails in bacteria necessitates an alternative approach for isolating mRNA (~5%) from the significantly more abundant rRNA (~95%). Here, by integrating a Ribosomal RNA derived cDNA Depletion protocol (RiboD) into a PETRI-seq, we developed RiboD-PETRI-seq that efficiently eliminates rRNA reads, thereby significantly improving mRNA detection rates and enabling exploration of within-population heterogeneity.

Results

In the RiboD protocol, we designed a set of probe primers that spans all regions of the bacterial rRNA sequence (Table S1). The core principle behind our probe design is twofold: the 3′-end of the probes are reverse complementary to the r-cDNA sequences, allowing for specific recognition of r-cDNA, while the 5′-end complements a biotin-labeled universal primer. This design enables the probes to be bound to magnetic beads, facilitating the separation of r-cDNA-probe-bead complexes from the rest of the library. Following template switching and RNaseH treatment on the barcoded cDNA from lysed cells to eliminate hybridized RNA, the library of probe primers and biotin-labeled universal primers is introduced to facilitate adequate hybridization. Pre-treated Streptavidin magnetic beads are then added to the hybridized rRNA-derived cDNA. The mRNA-derived cDNA remains in the supernatant and is collected for subsequent library construction and sequencing (Fig. 1A). To assess the efficiency of single-cell capture in RiboD-PETRI, we calculated the multiplet frequency¹² using a Poisson distribution based on our sequencing results (see details in Methods). The multiplet frequency for RiboD-PETRI ranges from 1.16% to 3.35% (Table S2), indicating the technique’s capability to effectively capture transcriptomes at the single-cell level.

Figure 1.

To assess the performance of RiboD-PETRI, we designed a comprehensive assessment of rRNA depletion efficiency under diverse physiological conditions, specifically contrasting exponential and stationary phases. This approach allows us to understand how these different growth states impact rRNA depletion efficacy. Additionally, we included a variety of bacterial species, encompassing both gram-negative and gram-positive organisms, to ensure that our findings are broadly applicable across different types of bacteria. By incorporating these variables, we aim to provide insights into the robustness and reliability of the RiboD-PETRI method in various biological contexts. The results highlight a substantial enhancement in rRNA-derived cDNA depletion, with mRNA ratio increases from 8.2% (Ctrl, the PETRI-seq we performed) to 81% (☒☒, RiboD-PETRI) for E. coli from exponential phase, from 10% (Ctrl) to 92% (☒☒) for S. aureus Trom stationary phase, and from 3.9% (Ctrl) to 54% (☒☒) for C. crescentus from exponential phase (Fig. 1B; Table S3). Additionally, we compared our findings with other reported methods (Fig. 1B; Table S4). The original PETRI-seq¹² protocol, which does not include an rRNA depletion step, exhibited an mRNA detection rate of approximately 5%. The MicroSPLiT-seq¹⁴ method, which utilizes Poly A Polymerase for mRNA enrichment, achieved a detection rate of 7%. Similarly, M3-seq¹⁶ and BacDrop-seq¹⁵, which employ RNase H to digest rRNA post-DNA probe hybridization in cells, reported mRNA detection rates of 65% and 61%, respectively. MATQ-DASH²⁰, which utilizes Cas9-mediated targeted rRNA depletion, yielded a detection rate of 30%. smRandom-seq utilizes a CRISPR-based rRNA depletion technique, reduced the rRNA proportion from 83% to 32%, increasing the mRNA proportion from 16% to 63%²¹. BaSSSh-seq’s employs a rational probe design for efficient rRNA depletion, though specific efficiency was not reported²². Among these, RiboD-PETRI demonstrated superior performance in mRNA detection while requiring the least sequencing depth. With equivalent sequencing depth, RiboD-PETRI demonstrates a significantly enhanced Unique Molecular Identifier (UMI) counts detection rate compared to PETRI-seq alone (Fig. 1C). This method recovered approximately 20,175 cells (92.6% recovery rate) with ≥ 15 UMIs per cell with a median UMI count of 42 per cell, which was significantly higher than PETRI-seq’s recovery rate of 17.9% with a median UMI count of 20 per cell (Figure S1A, B), indicating the number of detected mRNA per cell increased prominently. Notably, this enhancement was achieved while maintaining mRNA profiles consistent with non depleted samples (r = 0.93; Fig. 1D) and show a significant correlation with profiles from the traditional bulk RNA-seq method (r = 0.84; Fig. 1E).

We subsequently investigated the transcriptome coverage of RiboD-PETRI across different physiological states and bacterial species. For exponential phase E. coli cells, we sequenced a library with 60,000 cells, recovering approximately 30,004 cells (50% recovery), each with ≥15 UMIs (Fig. 2A and S1C). This analysis revealed 99.86% transcriptome-wide gene coverage across the cell population. The method achieved an average of 128.8 UMIs per single cell, with a median UMI count of 102 per cell. Further examination of high-quality cells showed varying levels of detection: the top 1,000, 5,000, and 10,000 cells exhibited median UMI counts of 462, 259, and 193, respectively (Fig. 2B), and median gene detection of 362, 236, and 188, respectively (Fig. 2C). These high-performing cells demonstrate the upper limits of the method’s capabilities. For stationary phase S. aureus cells, we sequenced a library with 30,000 cells, recovering approximately 9,982 cells (33.3% recovery), each with ≥15 UMIs (Fig. S2A). Analysis showed 99.96% transcriptome-wide gene coverage across the cell population. At the single-cell level, we observed an average of 153.8 UMIs and a median of 142 UMIs. Top high-quality cells exhibited the following median UMI counts: 378 (top 1,000 cells), 207 (top 5,000 cells), and 167 (top 8,000 cells) (Fig. S2B). These cells also demonstrated median gene detection of 308, 194, and 158 genes, respectively (Fig. S2C). For exponential phase C. crescentus cells, we sequenced a library with 30,000 cells, recovering approximately 13,897 cells (46.3% recovery), each with ≥15 UMIs (Fig. S2G). Analysis showed 99.64% transcriptome-wide gene coverage across the cell population. At the single-cell level, we observed an average of 439.7 UMIs and a median of 182 UMIs. Top high-quality cells demonstrated the following median UMI counts: 2,190 (top 1,000 cells), 662 (top 5,000 cells), and 225 (top 10,000 cells) (Fig. S2H). These cells also exhibited median gene detection of 1,262, 529, and 219 genes, respectively (Fig. S2I). These results underscore RiboD-PETRI’s ability to capture a wide range of transcripts across varying cell qualities and species, providing a comprehensive view of gene expression at the single-cell level.

Figure 2.

Our results affirm RiboD-PETRI’s reliability in capturing the bacterial single-cell transcriptome, providing ample coverage and sensitivity for various species. To provide a thorough evaluation of our sequencing depth and library quality, we performed sequencing saturation analysis on our sequencing samples. The findings reveal that our sequencing saturation is greater than 90% (Fig. S1D-F), indicating that our sequencing depth is sufficient to capture the diversity of most transcripts.

We further investigated its ability to consistently identify within-population heterogeneity across different bacterial species and growth conditions. In the exponential phase of E. coli, we recovered 1,464 cells and identified three major subpopulations (Fig. 2D), with 17 cells (1.2%) in a unique subpopulation characterized by pentose and glucuronate interconversions (Fig. 2E, F) and the marker genes of cluster 2 included yfgM, glpG, scpA, elyC, ptsA, carB, actP and pgpB (Fig. 2G). For the expression levels of marker gene shown in Fig. 2E, violin plots have been created to offer a more comprehensive view of the distribution across different cell populations (Fig. S4). In stationary phase S. aureus cells, we recovered 9,386 cells and found six major subpopulations (Fig. S2D), with 437 cells (4.7%) in a distinct subpopulation named cluster 4. The marker genes of cluster 4 included KQ76-13335, KQ76-00740 and KQ76-11725 (Fig. S5). In the stationary phase of C. crescentus cells, we recovered 5728 cells and identified four major subpopulations (Fig. S2J), with 603 cells (10.5%) in a unique subpopulation named cluster 3. The marker genes of cluster 3 included CCNA-00259, CCNA-03402, CCNA-02361 and CCNA-03119 (Fig. S6). These findings highlight RiboD-PETRI’s consistent ability to unveil within-population heterogeneity across different cell physiology and bacterial species (Fig. 2H, I and S2E, F, K, L), crucial for understanding bacterial population complexity. While RiboDPETRI consistently detects potential heterogeneity, further experimental validation would be required to confirm the biological significance of the observations.

We next focused on exploring biological heterogeneity of a biofilm at the early stage of development by utilizing the static biofilm system²³. E. coli cells were cultured in microtiter dishes overnight, adhered cells were fixed for RiboD-PETRI processing in duplicate experiments. For these two replicates, we sequenced libraries containing 20,000 and 40,000 cells, recovering 5,244 and 11,344 cells, which corresponded to recovery rates of 26% and 28%, respectively. The correlations between detected reads and unique molecular identifiers (UMIs) were found to be 0.87 and 0.90 for the two replicates, respectively (Fig. S3A). The median UMI counts for the recovered cells were 34 and 52 (Fig. 3A). After screening, the final datasets comprised 1,621 and 3,999 cells for the two replicates, respectively. While no significant batch effects were observed, we applied batch correction as a precautionary measure (Fig. S3B-D). In replicate 1, each cell was sequenced with an average of 1,563 reads, while in replicate 2, the average was 2,034 reads (Table S5), yielding median UMI counts of 283.5 and 239 per cell, respectively (Fig. 3B). For gene detection, the median counts were 219 and 193 genes per cell for the respective replicates (Fig. 3C). Additionally, UMAP visualization was employed to illustrate the distribution of cellular UMI numbers, revealing heterogeneity among populations that was independent of UMI counts (Fig. S3E, F). Unsupervised clustering analysis identified four major subpopulations in each replicate, with a consistently identified rare subpopulation (2.6%/2.1%) as cluster 2, driven by cell envelope genes (Fig. 3D-F). Marker genes for this cluster included yffO, lptE, rdgB, pdeI, sstT, fixA, yjjG, rlmI, accC and yaiA (Fig. 3G, H and S7).

Figure 3.

PdeI, identified among marker genes, was predicted as a phosphodiesterase enzyme hydrolyzing c-di-GMP, a vital bacterial second messenger^24,25 (Fig. 4A, B). However, our comprehensive structural analysis revealed a more complex and novel role for PdeI. While PdeI contains an intact EAL domain typically associated with c-di-GMP degradation, it also possesses a divergent GGDEF domain, generally linked to c-di-GMP synthesis (Fig. S8). This dual-domain architecture suggested potential complex regulatory roles. To validate PdeI’s function, we created a PdeI-BFP fusion construct under the native pdeI promoter, integrated with a ratiometric c-di-GMP sensing system²⁶ in E. coli. Confocal microscopy revealed PdeI as a membrane protein (Fig. 4C). Single-cell level monitoring showed cell-to-cell variability in c-di-GMP levels and PdeI expression, with a positive correlation observed (Fig. 4D), indicating PdeI upregulated c-di-GMP synthesis rather than degradation. This finding was confirmed by HPLC LC-MS/MS, which showed an approximately elevenfold increase in c-di-GMP concentration in the PdeI overexpression strain compared to the control strain (Fig. 4E). These results align with previous studies showing that a point mutation (G412S) in PdeI’s divergent GGDEF domain in a strain lacking PdeH, the major phosphodiesterase in E. coli, resulted in decreased c-di-GMP levels²⁷. Our additional experiments with a PdeI(G412S)-BFP mutation strain showed constant c-di-GMP levels despite increasing BFP fluorescence, serving as a proxy for PdeI(G412S) expression levels (Fig. 4D). These results, combined with the presence of a CHASE (cyclases/histidine kinase-associated sensory) domain in PdeI, strongly suggest that PdeI functions as a membrane-associated sensor that integrates environmental signals with c-di-GMP production under complex regulatory mechanisms. This discovery challenges the initial prediction of PdeI as solely a phosphodiesterase and highlights its novel role as a c-di-GMP synthetase, contributing significantly to our understanding of bacterial signaling pathways. It’s worth noting that while the other marker genes in this cluster are co-expressed, our analysis indicates that they do not have a significant impact on biofilm formation or a direct relationship with c-di-GMP or PdeI.

Figure 4.

Confocal laser scanning microscopy provided further insights into the spatial distribution of PdeI-positive cells within the biofilm structure. In the PdeI-BFP fusion strain, PdeI-BFP positive cells, characterized by elevated c-di-GMP levels, were predominantly located at the bottom of the static biofilm (Fig. 4F). This localization corresponds to the region of cell-surface attachment, aligning with our hypothesis that PdeI functions as a membrane-associated sensor integrating environmental signals with c-di-GMP production through complex regulatory mechanisms²⁸. In contrast, in the control strain where BFP was expressed alone under arabinose-induced promoter, BFP-positive cells were observed to be distributed throughout the entire biofilm community (Fig. 4G). This distinct spatial distribution pattern between the PdeI-BFP fusion and the BFP-only control strains provides compelling evidence for PdeI’s specific role in biofilm formation, particularly at the biofilm-surface interface. The concentration of PdeI-positive cells at the bottom of the biofilm suggests that PdeI may be especially crucial in the initial stages of biofilm formation, potentially responding to surface-associated cues to modulate c-di-GMP levels and promote attachment. The uniform distribution of BFP in the control strains suggests that the localization of PdeI-BFP is not affected by BFP labelling. These observations further underscore the complex and nuanced role of PdeI in bacterial signaling and biofilm development, highlighting the importance of considering cell-to-cell heterogeneity in understanding the function of regulatory proteins in microbial communities.

To investigate the association of the PdeI-high cluster with bacterial drug tolerance in the early stages of biofilm development, we isolated PdeI-high cells using flow cytometry (Fig. 4H, I) and subjected them to an ampicillin antibiotic killing assay to determine their persister frequency. Our results revealed that the PdeI-high population produced a significantly higher ratio of persister cells (~7.3%) compared to the whole biofilm population (~0.6%). Notably, cells expressing high levels of BFP alone or PdeI(G412S)-BFP showed no increase in persister ratios (Fig. 4J). This finding suggests that the increased persistence is specifically linked to PdeI activity. Time-lapse imaging during the antibiotic killing process consistently demonstrated that persisters primarily originated from PdeI-GFP positive cells (Fig. 4K and Movie S1). These PdeI-GFP positive cells, displaying characteristics of dormancy, survived ampicillin treatment for 6 hours without visible growth or division. Upon antibiotic removal and replacement with fresh growth medium, the PdeI-GFP positive persister cells resumed activity, elongating, dividing, and forming new microcolonies (Fig. 4K and Movie S1). This dynamic behavior provides direct visual evidence of the persister phenotype associated with PdeI-high cells. These findings strongly suggest that c-di-GMP, a molecule whose intracellular levels are upregulated by PdeI, plays a significant role in generating a persister subpopulation during the early stages of biofilm development. The mechanism by which elevated c-di-GMP levels contribute to antibiotic tolerance may involve modulation of cellular metabolism or activation of stress response pathways, leading to a state of dormancy that enables survival under antibiotic stress. This discovery not only enhances our understanding of the link between biofilm formation and antibiotic tolerance but also identifies PdeI as a potential target for strategies aimed at combating persistent bacterial infections.

Discussion

In this study, we introduce RiboD-PETRI, an enhanced bacterial scRNA-seq method that offers a cost-effective (Table S6), equipment-free, and high-throughput solution. By incorporating a probe hybridization-based rRNA-derived cDNA depletion protocol, our approach efficiently removes rRNA reads and significantly improves mRNA detection rates, enabling a more comprehensive exploration of within-population heterogeneity. At $0.0049 per cell, RiboD-PETRI is substantially more economical than the original PETRI-seq ($0.056 per cell), making it an attractive option for budget-conscious researchers. The method demonstrates improved mRNA detection, versatility across various bacterial species and growth conditions, preservation of transcriptome profiles consistent with non-depleted samples (r = 0.93) and traditional bulk RNA-seq methods (r = 0.84), high transcriptome coverage (>99%), and robust single-cell resolution with median UMI counts ranging from 102 to 182 per cell across different species and conditions.

The application of RiboD-PETRI to investigate biofilm heterogeneity exemplifies its potential for exploring complex biological systems. Our analysis of early-stage biofilm development uncovered a rare subpopulation (2.1-2.6%) characterized by cell envelope genes, including the previously uncharacterized gene pdeI. Further investigation revealed PdeI’s novel role as a c-di-GMP synthetase rather than a phosphodiesterase, challenging initial predictions and contributing significantly to our understanding of bacterial signaling pathways. Moreover, we demonstrated that PdeI-high cells exhibit increased antibiotic tolerance, with a significantly higher proportion of persister cells compared to the general biofilm population. This finding establishes a link between elevated c-di-GMP levels, regulated by PdeI, and the generation of antibiotic-tolerant subpopulations during early biofilm development. The spatial distribution of PdeI-positive cells at the bottom of the static biofilm, corresponding to the cell-surface attachment region, supports our hypothesis that PdeI functions as a membrane-associated sensor integrating environmental signals with c-di-GMP production. This localization pattern suggests PdeI’s crucial role in the initial stages of biofilm formation, potentially responding to surface-associated cues to modulate c-di-GMP levels and promote attachment. While other marker genes in this cluster are co-expressed, our analysis indicates they do not significantly impact biofilm formation or directly relate to c-di-GMP or PdeI.

In conclusion, RiboD-PETRI represents a significant advancement in bacterial scRNA-seq methodology. Its ability to uncover hidden variations within bacterial populations, as demonstrated in our biofilm analysis, underscores its potential impact on advancing our understanding of microbial behavior and population dynamics. By providing a cost-effective and efficient tool for exploring bacterial heterogeneity, RiboD-PETRI opens new avenues for research in microbiology, potentially leading to novel insights into antibiotic resistance, biofilm formation, and other critical areas of bacterial biology.

Materials and Methods

Resource availability

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yingying Pu (yingyingpu@whu.edu.cn).

Materials availability

Plasmids generated in this study are available from the lead contact upon request.

Data and code availability

Sequencing data, the processing data and original code have been deposited to GEO repository (GSE260458).

Bacterial strains and growth conditions

The bacterial strains used in this study included Escherichia coli strains MG1655, Caulobacter crescentus NA1000 and Staphylococcus aureus strain ATCC. E. coli cultures were grown in Luria Broth (LB) medium. For the biofilm setup, bacterial cultures were grown overnight. The next day, we diluted the culture 1:100 in a petri dish. We added 2ml of LB medium to the dish. If the bacteria contain a plasmid, the appropriate antibiotic needs to be added to LB. The petri dish was then incubated statically in a growth chamber for 24 hours. After incubation, we performed imaging directly under the microscope. The petri dishes used were glass-bottom dishes from Biosharp (catalog number BS-20-GJM), allowing for direct microscopic imaging without the need for cover slips or slides. This setup allowed us to grow and image the biofilms in situ, providing a more accurate representation of their natural structure and composition. Caulobacter crescentus strain NA1000 was grown in peptone yeast extract (PYE) medium. And Staphylococcus aureus strain ATCC 25923 was grown in Mueller-Hinton Broth (MHB) medium. All bacterial strains were routinely grown at 37 °C and 220 rpm. To maintain plasmids, when necessary, media were supplemented with chloramphenicol (25 pg/ml) or kanamycin sulfate (50 μg/ml). For arabinose-induction system expression experiments, 0.002% or 0.00002% arabinose was supplemented in the medium.

Strains construction

The construction of recombinant plasmids was performed using the 2 x MultiF Seamless Assembly Mix (ABclonal, RK21020). For the detection of c-di-GMP levels using c-di-GMP sensor and the detection of persister. the PdeI gene, pdeI, along with its native promoter (250bp), was fused with either gfp or bfp and cloned into the pBAD backbone. The original promoter region of the pBAD vector was removed to avoid any potential interference. This construction allows the expression of the BFP, PdeI-BFP and PdeI(G412S)-BFP fusion proteins to be driven by pdeI’s native promoter, thus maintaining its physiological control mechanisms. And the BFP coding sequence was fused to the pdel gene to create the PdeI-BFP fusion construct. Besides, for membrane localisation and localisation in biofilm community, the pdel-gfp and pdel-bfp with native promoter of pdel were cloned into the pBAD backbone. For the contro group, bfp and gfp was cloned into the pBAD backbone under the control of Arabinose-induction system. For HPLC LC-MS/MS analysis, pdeI and empty vector were cloned into the pBAD backbone, induced by arabinose. GFP and BFP were used in different experiments. GFP was used for imaging and time-lapse imaging to observe persister cell growth. BFP was used for cell sorting and detecting the proportion of persister cells. For the c-di-GMP sensor (Addgene: #182291), the plasmid origin was replaced with the p15A ori.

RiboD-PETRI

Cell preparation

E. coli MG1655 cells were cultured overnight and subsequently diluted at a ratio of 1:100 into fresh LB medium and grown statically for 24h at 37 °C. For 3h exponential period E. coli sample, E. coli MG1655 cells were grown overnight and then diluted 1:100 into fresh LB medium and grown for 3 h at 37 °C and 220 rpm. Caulobacter crescentus strain NA1000 cells were grown overnight and then diluted 1:100 into fresh Mueller-Hinton Broth (MHB) medium and grown for 9 h at 37 °C and 220 rpm. And Staphylococcus aureus strain ATCC 25923 cells were grown overnight and then diluted 1:100 into fresh peptone yeast extract (PYE) medium and grown for 3 h at 37 °C and 220 rpm. All the culture was vigorously shaken using a vortex, and the cells were then centrifuged at 5,000g for 2 minutes at 4 °C. The pellet was resuspended in 2 ml of ice-cold 4% formaldehyde (F8775, Millipore Sigma, diluted into PBS). These suspensions were rotated at 4 °C for 16 hours.

Cell permeabilization

1 ml of fixed cells were centrifuged at 5,000g for 5 minutes at 4 °C, then resuspended in 1 ml washing buffer (100 mM Tris-HCl pH 7.0, 0.02 U/μl SUPERase-In RNase Inhibitor, AM2696, Invitrogen). After another centrifugation at 5,000g for 5 minutes at 4 °C, the supernatant was removed. The pellet was then resuspended in 250 μl permeabilization buffer (0.04% Tween-20 in PBS-RI, PBS with 0.01 U/μl SUPERase-In RNase Inhibitor) and incubated on ice for 3 minutes. 1 ml cold PBS-RI was added, and the cells were centrifuged at 5,000g for 5 minutes at 4 °C. The pellet was resuspended in 250 μl Lysozyme Mix (250 μg/ml Lysozyme or 5 pg/ml Lysostaphin for S. aureus in TEL-RI buffer, comprising 100 mM Tris pH 8.0 (AM9856, Invitrogen), 50 mM EDTA (AM9261, Invitrogen), and 0.1 U/μl SUPERase In RNase Inhibitor). The samples were incubated at 37 °C and mixed gently every minute. Then 1 ml cold PBS-RI was added immediately, and cells were centrifuged at 5,000g for 5 min at 4 °C. The cells underwent another wash with 1 ml cold PBS-RI. Subsequently, cells were resuspended in 40 μl DNaseI-RI buffer (4.4 μl 10× reaction buffer, 0.2 μl SUPERase In RNase inhibitor, 35.4 μl H₂O), followed by addition of 4 μl DNaseI (AMPD1, Millipore Sigma). The samples were incubated for 30 min at room temperature and mixed gently every 5 minutes. 4 μl Stop Solution was added, and the samples were incubated for 10 minutes at 50 °C with gentle mixing every minute. Following centrifugation at 5,000g for 10 minutes at 4 °C, cells were washed twice with 0.5 ml cold PBS-RI. Finally, cells were resuspended in 200 μl cold PBS-RI, and their count and integrity were assessed using the ACEA NovoCyte flow cytometer with a 100× oil immersion lens.

Primer preparation

For round 1 RT, round 2 and round 3 ligation reactions, all primers design and preparation as previously described¹². All primers were purchased from Sangon Biotech (Table S1). For ligation primers preparation, mixtures were prepared as follows: 31.1 μl each R2 primer (100 μM), 28.5 μl SB83 (100 μM) and 21.4 μl H₂O, were splitted to 2.24 μl for one sample. Mixtures containing 63.2 μl each R3 primer (70 μM) and 58 μl SB8 (70 μM), were splitted to 3.49 μl for one sample. Before use, ligation primers were incubated as follows: 95 °C for 3 min, then decreasing the temperature to 20 °C at a ramp speed of −0.1 °C s−1, 37 °C for 30 min. For blocking mix preparation, 50 μl primer SB84 (400 μM) and 80 μl primer SB81 (400 μM) were incubated as follows: 94 °C for 3 min, then decreasing the temperature to 25 °C at a ramp speed of −0.1 °C s−1, 4 °C for keeping. Round 2 blocking primers were mixed as follows: 37.5 μl 400 μM SB84, 37.5 μl 400 μM SB85, 25 μl 10× T4 ligase buffer, 150 μl H₂O. Round 3 blocking primers were mixed as follows: 72 μl 400 μM SB81, 72 μl 400 μM SB82, 120 μl 10× T4 ligase buffer, 336 μl H₂O, 600 μL 0.5 M EDTA.

Round 1 RT reaction

About 3×10⁷ cells were introduced into an RT reaction mix composed of 240 μl 5× RT buffer, 24 μl dNTPs (N0447L, NEB), 12 μl SUPERase In RNase Inhibitor and 24 μl Maxima H Minus Reverse Transcriptase (EP0753, Thermo Fisher Scientific), 132 μL PEG8000 (50%). Nuclease-free water was added to achieve a total reaction volume of 960 μl, and the mixture was thoroughly mixed by vortexing. Subsequently, 8 μl of the reaction mixture was dispensed into each well of a 96-well plate, where 2 μl of each RT primer had been added previously. The sealed 96-well plate was inverted repeatedly for thorough mixing, followed by a brief spin. The plate was then incubated as follows: 50 °C for 10 min, 8 °C for 12 s, 15 °C for 45 s, 20 °C for 45 s, 30 °C for 30 s, 42 °C for 6 min, 50 °C for 16 min, and finally held at 4 °C. After the RT process, all 96 reactions were pooled into one tube. 75 μl of 0.5% Tween-20 was added, and the reactions were incubated on ice for 3 minutes. Cells were centrifuged at 7,000g for 10 minutes at 4 °C and then resuspended in 0.4 ml PBS-RI. Thirty-two microliters of 0.5% Tween-20 was added, and the cells underwent centrifugation at 7,000g for 10 minutes at 4 °C.

Round 2 ligation reaction

Cells were resuspended in 500 μl 1× T4 ligase buffer, followed by the addition of 107.5 μl PEG8000, 37.5 μl 10× T4 ligase buffer, 16.7 μl SUPERase In RNase Inhibitor, 5.6 μl BSA, and 27.9 μl T4 ligase (M0202L, NEB). The reaction solution was thoroughly mixed by vortexing. Subsequently, 5.76 μl of the reaction mixture was dispensed into each well of a 96-well plate, where 2.24 μl of each round 2 ligation primer had been added previously. The sealed 96-well plate was inverted repeatedly for thorough mixing and then subjected to a short spin. The plate was incubated at 37 °C for 45 minutes. Following this, 2 μl of round 2 blocking mix was added to each well and incubated at 37 °C for an additional 45 minutes. All 96 reactions were pooled into one tube after incubation.

Round 3 ligation reaction

A mixture comprising 89 μl H₂O, 26 μl PEG8000, 46 μl 10× T4 ligase buffer, and 12.65 μl T4 ligase was prepared and thoroughly mixed by vortexing. Subsequently, 8.51 μl of the reaction mixture was dispensed into each well of a 96-well plate, where 3.49 μl of each round 3 ligation primer had been added previously. The sealed 96-well plate was inverted repeatedly for thorough mixing and then subjected to a brief spin. The plate was incubated at 37 °C for 45 minutes. Following this, 10 μl of round 3 blocking mix was added to each well and incubated at 37 °C for an additional 45 minutes. All 96 reactions were combined into one tube after incubation.

Cells lysis

42 μl of 0.5% Tween-20 was added, and cells were centrifuged at 7,000g for 10 minutes at 4 °C. The cells underwent two washes using 200 μl TEL-RI containing 0.01% Tween-20, each time centrifuged at 7,000g for 10 minutes at 4 °C. Subsequently, cells were resuspended in 30 μl TEL-RI buffer. Cell counting and integrity checks were performed using the ACEA NovoCyte flow cytometer with a 100× oil immersion lens. A moderate amount of cells was then added to the lysis buffer (50 mM Tris pH 8.0, 25 mM EDTA, 200 mM NaCl, 0.5% Triton X-100), and 5 μl of proteinase K (AM2548, Invitrogen) was introduced. Samples were incubated at 55 °C for 60 minutes and gently mixed every minute.

Library construction

To facilitate template switching, lysates were purified with VAHTS DNA Clean Beads (N411, Vazyme) at a ratio of 2.0×, and cDNA was eluted in 12 μl of water. The purified cDNA was then combined with 4 μl of 5× RT buffer, 1 μl of dNTPs (N0447L, NEB), 0.5 μl of SUPERase In RNase Inhibitor, 0.5 μl of Maxima H Minus Reverse Transcriptase, and 2 μl of the TSO²⁹ primer (100 mM, Table S1). This reaction solution underwent incubation as follows: 25 °C for 30 min, 42 °C for 90 min, 85 °C for 5 min, and then held at 4 °C. Subsequently, 1 μl of RNaseH was added, and the reaction solution was incubated at 37 °C for 30 min. The cDNA was purified once again with VAHTS DNA Clean Beads at a ratio of 2.0× and eluted in 13 μl of H₂O. The integrity of the cDNA was assessed using primers TSO-2 and R1 or R2 or R3 by qPCR (Table S1).

Ribosomal RNA derived cDNA depletion (RiboD)

We developed a set of cDNA probe primers to selectively deplete r-cDNA (Table S1). These probe primers possess the ability to specifically hybridize with r-cDNA and also hybridize with a biotin-labeled universal primer. In the reaction, 5 μl of r-cDNA probe primers (10 μM), 2.5 μl of 10× hybridization buffer (Tris-HCl pH 8.0 100 mM, NaCl 500 mM, EDTA pH 8.0 10 mM), and 5 μl of biotin primer (10 μM) were added to 12.5 μl of purified cDNA. The reaction solution underwent incubation as follows: 95 °C for 2 min, followed by a temperature decrease to 20 °C at a ramp speed of −0.1 °C s−1, and then held at 37 °C for 30 min. Subsequently, 20 μl of Streptavidin magnetic beads (BEAVER, 22307) was washed twice using 1 ml of 1× B&W buffer (Tris-HCl pH 7.5 10 mM, EDTA 1 mM, NaCl 1M, Tween-20 0.05%) and resuspended in 25 μl of 2× B&W buffer. Twenty-five microliters of washed Streptavidin magnetic beads were added to 25 μl of annealed cDNA. The reaction solution was incubated at room temperature for 30 min with gentle mixing per minute. Following this, the reaction solution tube was placed into a magnetic stand to collect the supernatant. The cDNA depleted of r-cDNA was purified using VAHTS DNA Clean Beads at a ratio of 2.0× and eluted in 12.5 μl of H₂O. The depletion of r-cDNA could be repeated using the above protocol, and ultimately, the cDNA was eluted in 20 μl of H₂O. We designed separate probe sets for E. coli, C. crescentus, and S. aureus. Each set was specifically constructed to be reverse complementary to the r-cDNA sequences of its respective bacterial species. This species-specific approach ensures high efficiency and specificity in rRNA depletion for each organism.

Library amplification and sequencing

To the 20 μl cDNA solution, the following components were added: 2.4 μl R3 primer (10 mM, Table S1), 2.4 μl TSO-2 primer (10 mM, Table S1), 40 μl 2× KAPA HIFI mix (KAPA, 2602), 1.6 μl Sybr green (25×), 0.8 μl MgCl₂ (0.1 M), and 12.8 μl H₂O. This PCR reaction solution was placed in a thermocycler and incubated with the following parameters: 98 °C for 45 s, followed by cycling of 98 °C for 15 s, 60 °C for 30 s, and 72 °C for 60 s. Cycling continued on a qPCR machine until the reaction approached saturation. PCR products were then purified using VAHTS DNA Clean Beads at a ratio of 0.9× and eluted in 25 μl of H₂O. Finally, the purified PCR products underwent end repair and adaptor ligation using the VAHTS Universal DNA Library Prep Kit for Illumina V3 (Vazyme, ND607).

Bulk RNA-seq library construction

Total RNA of the samples was extracted utilizing the Bacteria RNA Extraction Kit (R403-01, Vazyme). Subsequently, the RNA underwent mRNA enrichment (N407, Vazyme), fragmentation, cDNA synthesis, and library preparation using the VAHTSTM Total RNA-seq (H/M/R) Library Prep Kit for Illumina® (NR603, Vazyme).

Bioinformatics analysis Methods

Single-Cell Analysis

The sequencing data underwent processing into matrices using scripts and a pipeline as previously described¹² in Python 2.7.15, with some modifications (the detailed original code and all the data were deposited in the GEO repository). After the count tables were made, subsequent analysis of single-cell data was conducted using Seurat³⁰ package (version 4.3.0; http://satijalab.org/seurat/) in R (https://www.r-project.org/). Since there were two replicates of static E. coli biofilm, these two datasets were merged into one SeuratObject and batch effects were removed. However, the samples for exponential period E. coli, S. aureus and C. crescentus only had one sample, so they did not need this process. At the beginning of doing the scRNA-seq analysis, we screened the data of all samples. For preprocessing of static E. coli biofilm data, cells were filtered with UMI per cell more than 100 and less than 2000 for replicate 1 and replicate 2 to obtain 1621 and 3999 cells, respectively. For data of exponential period E. coli, the data was screened for cells with UMIs greater than 200 and less than 5000 to obtain 1464 cells. The screening criteria of S. aureus were cells with UMIs greater than 15 and less than 1000 and genes greater than 30 (1000>UMIs>15, gene counts>30). The screening criteria of C. crescentus were cells with UMIs greater than 200 and less than 5000 and gene counts greater than 30 (5000>UMIs>200, gene counts>30). After screening, all the data were normalized using a scale factor of 10000 through a global-scaling normalization method called “LogNormalize”. Highly variable features were then identified, returning 500 features per dataset. Then we combine the data of the two replicates of static E. coli biofilm into a single SeuratObject by FindIntegrationAnchors and IntegrateData functions. Then all the data underwent scaling using the ScaleData function, followed by dimension reduction through principal component analysis. To avoid subtle batch effects influencing downstream analyses, we removed batch effects using RunHarmony³¹ for the two replicates of static E. coli biofilm data. Then a graph-based clustering approach was employed in all data to identify clusters of gene expression programs using the Louvain algorithm (Seurat 4.3.0). The dims we choosed were 6. And the resolution was 0.3 for C. crescentus and S. aureus or 0.4 for E. coli data. Marker genes for each cluster were computed using the Wilcoxon Rank-sum test. Specifically, marker genes for each cluster were initially obtained using the FindMarkers function of Seurat. Then we performed pathway enrichment analysis of marker gene by clusterProfiler function³² within R. For transcriptome-wide gene coverage across the cell population, we counted the number of genes expressing at least one UMI. Then we calculated the percentage of these genes out of all the genes in each bacterium.

Comparison of scRNA-seq with Bulk RNA-Seq

The bulk RNA-seq clean data reads were mapped to the E. coli MG1655 k12 genome (EnsemblBacteria Taxonomy ID: 511145) using the BWA aligner software (v0.7.17-r1188, https://github.com/lh3/bwa.git). Sam files were converted to bam files using samtools (v1.9). The mapping results were counted by featureCounts (v2.0.1, https://github.com/topics/featurecounts) to generate expression results. Single-cell and bulk transcriptomes of E. coli were compared by computing the Pearson correlation of log₂ reads per gene of bulk RNA-seq and log₂ UMI per gene of scRNA-seq.

Sequencing saturation of the libraries

To assess sequencing saturation, we generated five subsamples from the single-cell sequencing data, representing 20%, 40%, 60%, 80%, and 100% of the total data. Each subsample was analyzed independently following the previously described single-cell sequencing data analysis process. We then created gene expression matrices for each subsample based on the analysis results. The number of UMIs or genes was counted for each cell detected in each subsample. Next, we sorted the cells in descending order based on their UMI counts or gene counts and selected different numbers of cells from this sorted list, starting from those with the highest UMI or gene counts. For each selection of cells, we calculated the median number of UMIs or genes. Finally, we created a line graph representing the median UMI counts and genes using GraphPad Prism 9 software, allowing us to visualize the sequencing saturation across the different subsamples. In addition, we used the saturation calculation method of 10× Genomics to further detect the saturation of the data. The formula for calculating this metric is as follows: Sequencing Saturation = 1 - (n_deduped_reads / n_reads). Given the differences between RiboD-PETRI and 10× Genomics datasets, we have adapted the calculation as follows: n_deduped_reads: The number of UMIs as a measure of unique reads. n_reads: The total number of confidently mapped reads.

Multiplet frequency determination

Determination of the multiplet frequency was essential in assessing the efficiency of single-cell capture in RiboD-PETRI. This frequency is defined as the probability that a non-empty barcode corresponds to more than one cell. To calculate it, we used a Poisson distribution-based approach involving several key steps. Initially, we calculated the proportion of barcodes corresponding to zero cells using the formula . Then, we determined the proportion for one cell, , and derived the proportions for more than zero cells P(≥1)=1−P(0) and more than one cell P(≥2)=1−P(1)−P(0). These values allow for the calculation of the multiplet frequency as . The parameter λ plays a vital role in this model, representing the ratio of the number of cells to the total number of possible barcode combinations.

Flow cytometry sorting of bacteria and analysis

All samples were measured using a Beckman CytoFLEX SRT flow cytometer with a 70 pm nozzle, using normal saline as sheath fluid. During the 24-hour static biofilm growth phase, strains labeled with BFP, PdeI-BFP, PdeI(G412S)-BFP or c-di-GMP sensors were washed and resuspended in sterile PBS. Microorganisms were identified based on forward scatter (FSC) and side scatter (SSC) parameters. Cells were sorted into distinct groups according to their fluorescence intensity, with V450 used for BFP, FITC for mVenusNB, and ECD for mScarlet-I. The resulting data were subsequently analyzed using FlowJo V10 software (Tree Star, Inc.).

Antibiotic killing and persister counting assay

Cells sorted by flow cytometry were resuspended in fresh LB broth supplemented with 150 μg/ml ampicillin. The suspension was then incubated at 37 °C for 3 hours with continuous shaking at 220 rpm. To determine the initial cell count, an aliquot of the cell suspension was taken before the ampicillin challenge, serially diluted, and plated on LB agar plates for colony-forming unit (CFU) enumeration. These plates were incubated overnight at 37 °C. Following the ampicillin challenge, cells were harvested by centrifugation, washed once with sterile PBS to remove residual antibiotic, and resuspended in fresh PBS. This suspension was then serially diluted and plated on LB agar plates for post-challenge CFU counts. The persister ratio was calculated as the number of CFUs after ampicillin challenge divided by the number of CFUs before challenge. All experiments were performed in triplicate, with the results presented as means ± standard deviations of three independent biological replicates.

Microscopy

Bright-field and fluorescence imaging

Inverted microscopes, specifically the Nikon Eclipse Ti2 and Leica Stellaris 5 WLL, were employed for imaging, utilizing different lasers for illumination: a 405 nm laser for BFP and a 488 nm laser for GFP. Fluorescence emission signals were captured using an sCMOS camera (pco.edge 4.2 bi). Dedicated filter sets corresponding to the spectral characteristics of each fluorophore were utilized. Image analysis was performed with ImageJ software (Fiji). For the analysis of the c-di-GMP sensor, the ratio of mVenusNB to mScarlet-I (R) displayed a negative correlation with c-di-GMP concentration. Consequently, the value of R⁻¹ demonstrated a positive correlation with c-di-GMP concentration.

Time-lapse imaging

To investigate the processes of antibiotic killing and bacterial resuscitation, cells labeled with PdeI-GFP during the 24-hour static growth phase were collected and washed twice with PBS. These cells were then imaged on a gel pad composed of 3% low melting temperature agarose in PBS, which was prepared as a gel island in the center of the FCS3 chamber. The cells were observed under either bright-field or epifluorescence illumination. Following the imaging, the gel pad was surrounded by LB broth containing 150 pg/ml ampicillin, and the cells were incubated for 6 hours at 35 °C. Fresh LB broth was subsequently flushed in, and the growth medium was refreshed every 3 hours, ensuring sufficient recovery time for the cells.

Determination of c-di-GMP concentration by HPLC-MS/MS

The determination of c-di-GMP concentration by HPLC-MS/MS involved a series of steps. Initially, MG1655 Aara ΔBAD::pdel and MG1655 Δara pBAD::empty-vector strains were grown to mid-exponential growth phase, followed by induction with 0.002% arabinose. After a 2-hour incubation period, cells were harvested and washed with PBS. The washed cells were then rapidly frozen using liquid nitrogen.

Simultaneously, another portion of washed cells was stained with SYTO^TM24 and quantified using flow cytometry. The determination of c-di-GMP concentration was conducted by Wuhan Lixinheng Technology Co. Ltd. through high-pressure liquid chromatography-tandem mass spectrometry (HPLC-MS/MS). In the experiment, first, for cell samples, addition of 500 μL of extract solvent (precooled at -20 °C, acetonitrile-methanolwater, 2:2:1), the samples were vortexed for 30 s, homogenized at 38 Hz for 4 min, and sonicated for 5 min in ice-water bath. The homogenate and sonicate circle were repeated for 3 times, followed by incubation at -20 °C for one hour and centrifugation at 12000 rpm and 4 °C for 15 min. An 80 μL aliquot of the clear supernatant was transferred to an auto-sampler vial for LC-MS/MS analysis. The UHPLC separation was carried out using an Waters ACQUITY H-class plus UPLC System, equipped with Agilent ZORBAX Eclipse Plus C18 column (2.1 mm × 150 mm, 1.8 μm). An Waters Xevo TQ-XS triple quadrupole mass spectrometer, equipped with an electrospray ionization (ESI) interface, was applied for assay development. All strains were assayed in biological triplicates, and the measured values were converted into intracellular c-di-GMP concentrations (pg) per cell.

Quantification and statistical analysis

Statistical analysis was conducted using GraphPad Prism 9 software for Windows. The significance of results was determined using a two-tailed Student’s t-test. Error bars in the data representation indicate the standard deviations of the mean from a minimum of three independent experiments. A significance threshold of P < 0.05 was applied throughout the analysis. To denote significant differences in the results, a system of asterisks was used: * for P < 0.05, ** for P < 0.01, *** for P < 0.001, and **** for P < 0.0001. This comprehensive approach ensured a thorough and statistically sound analysis of the c-di-GMP concentration in the studied bacterial strains.

Additional Declarations

The authors declare no competing interests.

Author Contributions

Conceptualization: Y.Z., H.L., Y.P.; Methodology: H.L., X.Y., C.W., C.H.; Investigation: H.L., X.Y., C.W., C.H.; Visualization: H.L., C.W.; Bioinformatics analysis: X.Y.; Clinical stain isolation: C.G.; Supervision: Y.P.; Writing - original draft: H.L., Y.P.; Writing - review & editing: H.L., Y.P.

Competing Interest Statement

Authors declare that they have no competing interests.

Classification

Biological Sciences

Acknowledgements

We thank Prof. Fan Bai (Peking University) for valuable discussions. We thank Drs. Jidong Xing and Ziyang Liu for support in bioinformatics. We also thank the members of our laboratory for helpful discussion. This work is supported by the grants to Y.P. from the National Key R&D Program of China (2021YFC2701602), the National Natural Science Foundation of China (31970089), Science Fund for Distinguished Young Scholars of Hubei Province (2022CFA077), Major Project of Guangzhou National Laboratory (GZNL2024A01023), and the Fundamental Research Funds for the Central Universities (2042022dx0003). This work is also supported by the grants to Y.Z. for the National Science Fund for Distinguished Young Scholars (82025011) and C.G. from the Natural Science Foundation of Yunnan Province of China (202001BB050005). We also thank all the staff in the Core Facilities of Medical Research Institute at Wuhan University and the Core Facilities at School of Life Sciences at Peking University for their technical support.

Supplementary Figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.