Introduction

Biofilms, comprising approximately 80% of chronic and recurrent microbial infections in the human body, contain bacterial cells existing in a diverse array of physiological states 13. This diversity in states reflects a division of labor within the biofilm, contributing to increased resistance to various stresses4. However, the study of biofilms faces significant limitations, primarily stemming from challenges in investigating heterogeneity within a bacterial population. Single-cell RNA-seq emerges as a promising avenue for addressing this 512. 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) 7 and Microbial Split-Pool Ligation Transcriptomics (microSPLiT)9, 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. The 3’-end of the probes is specifically complementary to the rRNA-derived cDNA, while the 5’- end complements a biotin-labeled universal primer. 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 eliminate the hybridized rRNA-derived cDNA. The mRNA-derived cDNA remains in the supernatant and is collected for subsequent library construction and sequencing (Fig. 1A). The Multiplet frequency of this method is 1.16% – 3.35%, demonstrating the ability of RiboD-PETRI to capture transcriptomes of single cells (SI)7.

Development of RiboD-PETRI and validation of its technical performance in studying population heterogeneity.

(A) Graphic summary of the RiboD-PETRI method illustrating the incorporation of RiboD after cell pooling and lysis in PETRI-seq. The RiboD protocol is represented by the dashed-line box. (B) Comparison of non-rRNA (tRNA, mRNA and other non-rRNA) and rRNA UMI counts ratio among different bacterial scRNA-seq methods. Data from PETRI-seq (E col/)8, MicroSPLIT-seq (E coli)9, M3-seq (E coli)11 cited from from previous studies. Error bars represent standard deviations of biological replicates. (C) Comparison of UMI counts per cell between RiboD-PETRI and PETRI at the same unsaturated sequencing depth. (D) Assessment of the effect of rRNA depletion on transcriptional profiles. The Pearson correlation coefficient (r) of UMI counts per gene (log2 UMIs) between RiboD-PETRI and PETRI was calculated for 3790 out of 4141 total genes, excluding those with zero counts in either library. Each point represents a gene. (E) Evaluation of the correlation between RiboD-PETRI data and bulk RNA-seq results. The Pearson correlation coefficient (r) of UMI counts per gene (log2 UMIs) among RiboD-PETRI data and the reads per gene (log2 reads) of bulk RNA-seq data was calculated for 3814 out of 4141 total genes, excluding those with zero counts in either library. Each point represents a gene. (F-I) RiboD-PETRI data from E coli in exponential phase (E coli 3h). (F) Distributions of mRNA UMIs captured per cell in RiboD-PETRI data of E coli in exponential phase. (G) Uniform Manifold Approximation and Projection (UMAP) visualization of E coli bacteria during the exponential phase. Data were filtered for cells with UMIs between 200 and 5,000, resulting in 1,464 cells. Each dot represents a cell. (H) Heatmap illustrating the normalized gene expression levels of marker genes in different clusters of E. coli during the exponential phase. Marker genes with relatively high expression levels are depicted in yellow, while lower expression levels are shown in purple. Each row represents a gene, and each column represents a cell. (I) Functional enrichment analysis of marker genes in cluster 2. Marker genes were selected based on screening criteria of p-value < 0.001 and log2 fold change (FC) > 0.2.

To assess the performance of RiboD-PETRI, we conducted a comprehensive evaluation of rRNA depletion efficiency across gram-negative and gram-positive bacterial species. The results highlight a substantial enhancement in rRNA-derived cDNA depletion, with mRNA ratio increases from 8.2% to 81% for E. coli from exponential phase, from 10% to 92% for S. aureus from stationary phase, and from 3.9% to 54% for C. crescentus from exponential phase (Fig. 1B). With equivalent sequencing depth, RiboD-PETRI demonstrates a significantly enhanced Unique Molecular Identifier (UMI) counts detection rate compared PETRI-seq alone (Fig. 1C), 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 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. Analysis revealed 99.86% gene expression collectively and an average of 128.8 UMIs per single cell. Top 1,000, 5,000, and 10,000 high-quality cells showed averages of 462, 259, and 193 detected UMIs, respectively (Fig. 1F). 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. Analysis showed 99.96% gene expression collectively, and at the single-cell level, an average of 153.8 UMIs. Top 1,000, 5,000, and 8,000 high-quality cells exhibited averages of 273, 122, and 94 detected UMIs, respectively (Fig. S1A). In the exponential phase of C. crescentus cells, a library with 30,000 cells was sequenced, recovering approximately13,897 cells (46.3% recovery), each with >15 UMIs. Analysis showed 99.64% gene expression collectively, and at the singlecell level, an average of 439.7 UMIs. Top 1,000, 5,000, and 10,000 high-quality cells demonstrated averages of 2190, 662, and 225 detected UMIs, respectively (Fig. S1F). This performance surpassed other reported bacterial scRNA-seq methods7,911,13. Considering that this large library was sequenced with 1,508 reads per cell for E. coli, 2,372 for S. aureus, and 8,790 for C. crescentus, it is anticipated that a higher number of UMIs per cell will be detected with increased sequencing depth.

Our results affirm RiboD-PETRI’s reliability in capturing the bacterial single-cell transcriptome, providing ample coverage and sensitivity for various species. We 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. 1G), with 17 cells (1.2%) in a unique subpopulation characterized by pentose and glucuronate interconversions (Fig. 1H, I). For stationary phase S. aureus cells, we recovered 9,386 cells and found six major subpopulations (Fig. S1D), with 437 cells (4.7%) in a distinct subpopulation named cluster 4. In the stationary phase of C. crescentus cells, we recovered 5728 cells and identified four major subpopulations (Fig. S1I), with 603 cells (10.5%) in a unique subpopulation named cluster 3. These findings highlight RiboD-PETRI’s consistent ability to unveil within-population heterogeneity, crucial for understanding bacterial population complexity.

We next focused on exploring biological heterogeneity of a biofilm at the early stage of development by utilizing the static biofilm system14. E. coli cells were cultured in microtiter dishes overnight, adhered cells were fixed for RiboD-PETRI processing in duplicate experiments, each consisting of 1621 or 3999 cells. No significant batch effect was observed (Fig. S2C, H). Replicate 1 was sequenced 1,563 reads per cell, while replicate 2 with 2,034 reads per cell, yielding an average of 283/239 UMIs per cell (Fig. 2A). 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. 2B, D). Marker genes for this cluster included yffO, lptE, rdgB, pdel, sstT, fixA, yjjG, rlml , yncD, accCand yaiA (Fig. 2C, E, F). PdeI, identified among marker genes, was predicted as a phosphodiesterase enzyme hydrolyzing c-di-GMP, a vital bacterial second messenger (Fig. 2G, H). To validate this, we created a Pdel-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. 2I). Single-cell level monitoring showed cell-to-cell variability in c-di-GMP levels and PdeI expression, with a positive correlation observed (Fig. 2J). PdeI upregulated c-di-GMP synthesis rather than degradation, confirmed by both microscopy and 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. 2K).

RiboD-PETRI Resolves Biofilm Heterogeneity and Identifies New Marker Genes.

(A-D) RiboD-PETRI data from 24-hour static biofilms of E coli (E coli 24h). (A) UMI counts per cell in RiboD-PETRI data of 24-hour static biofilms of E coli. Two replicates were screened for cells with UMIs between 100 and 2000, resulting in 1621 and 3999 cells, with median UMIs per cell of 283.5 and 239, respectively. (B) UMAP visualization of 24-hour static biofilms of E. coli, revealing two small populations of heterogeneous cells in clusters 2 and 3. (C) Inferred expression levels of marker genes from 24-hour static biofilms of E. coli across different clusters. (D) Enrichment pathways for marker genes in cluster 2, selected based on screening criteria of p-value < 0.001 and log2 fold change (FC) > 0.2. (E & F) Dot plot displaying scaled expression levels of marker genes in different clusters of E. coli in exponential phase (E) and E. coli in 24- hour static biofilms (F). These genes were markers of static E. coli biofilms in cluster 2, identified with screening criteria of p-value < 0.001 and logFC > 3. Dot size represents the percentage expression of the gene in the cluster, while color indicates the average expression level normalized from 0 to 1 across all clusters for each gene. (G & H) UMAP plots showing the distribution of pdel in single-cell data of E. coli in exponential phase (G) and E. coli in 24-hour static biofilms (H). Each dot represents a cell colored by normalized expression levels of genes. (I) Subcellular localization of Pdel-GFP. (J) c-di-GMP levels (R-1 score) in E. colicells with different PdeI-GFP expression levels (low or high) in 24-hour static biofilms. c-di-GMP levels are measured using the c-di-GMP sensor system15 integrated into E. coli cells. Scale bar, 2 μm. R-1 score was determined using the fluorescent intensity of mVenusNB and mScarlet-I in the system. The fluorescent intensity is measured by sorted by fluorescence-activated cell sorting (FACS). (K) Determination of cellular concentrations of c-di-GMP by HPLC-MS/MS in cells overexpressing Pdel under the control of arabinose promoter, with 0.002% arabinose induction for 2 h (n=3). (L) Localization of Pdel-high cells in the biofilm matrix. Cells expressing Pdel-BFP under the control of the pdel native promoter were grown in a glass-bottom cell culture dish and stained with SYTO™ 24 for bacterial DNA. (M & N) Heterogeneous expression of Pdel in single-cell data of E. coli in exponential phase (M) and E. coli in 24- hour static biofilms (N). Biofilm cells with high or low expression levels of Pdel-BFP were sorted by FACS. (O) Persister counting assay using 150 pg/ml ampicillin on cells with high or low expression levels of Pdel-BFP from 24h static biofilm E. coli, sorted by FACS. (P) Time-lapse images of the persister assay observed under a microscope. Static biofilm cells of the Pdel-GFP strain were spotted on a gel pad and treated with 150 μg/ml ampicillin in LB broth. lmages were captured over 6 hours at 37°C, followed by the replacement of fresh LB broth to allow persister cell resuscitation. Scale bar, 2 μm. Error bars represent standard deviations of biological replicates. Significance was ascertained by two-tailed Student’s t test. Statistical significance is denoted as *P <0.05, **P < 0.01, ***, P<0.001, ****, P<0.0001.

Confocal laser scanning microscopy confirmed that PdeI-positive cells, exhibiting elevated c-di-GMP levels, were predominantly located at the bottom of the static biofilm (Fig. 2L), corresponding to the region of cell-surface attachment. To investigate the association of the PdeI-high cluster with bacterial persistence in the early stage of biofilm development, PdeI-high cells were isolated by fluorescence- activated cell sorting (FACS) (Fig. 2M, N) and subjected to an ampicillin antibiotic killing assay to determine their persister frequency. Results show that the PdeI-high population produced a significant higher ratio of persister cells (~7.3%) compared to that of the whole biofilm population (~0.6%) (Fig. 2O). Time-lapse imaging during the antibiotic killing process consistently revealed that persisters mainly originated from PdeI-GFP positive cells (Fig. 2P and Movie S1). PdeI-GFP positive cells, displaying dormancy, survived ampicillin treatment for 6 hours without 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. 2P and Movie S1). These findings 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.

Discussion

In conclusion, we report an improved approach, called RiboD-PETRI, which offers a cost-effective ($0.0049 per cell in RiboD-PETRI, and $0.056 per cell in PETRI-seq), equipment-free, and high-throughput solution for bacterial scRNA-seq. By integrating a probe hybridization-based rRNA derived cDNA depletion protocol, our method efficiently removes rRNA reads and significantly enhances mRNA detection rates. This improvement enables us to explore and analyze within-population heterogeneity. One important advantage of RiboD is the depletion process takes place at the cDNA level, following cell pooling and lysis. Consequently, RiboD only requires customization in the library preparation step. By implementing RiboD, we effectively reduce the necessary sequencing depth and overall sequencing costs, making it an attractive option for cost-conscious researchers. Beyond serving as a proof of principle, our study exemplifies the application of RiboD-PETRI to investigate the heterogeneity within biofilms, specifically delving into the early stages of biofilm development. This showcases the versatility and efficacy of RiboD-PETRI in exploring complex biological systems. The ability to uncover hidden variations within bacterial populations, as demonstrated in our biofilm analysis, underlines the potential impact of RiboDPETRI on advancing our understanding of microbial behavior and population dynamics.

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 (Uploading).

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 25923 were included. Cultures of E. coli were grown in Luria Broth (LB) medium. 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 μg/ml). For arabinose-induction system expression experiments, 0.002% arabinose was supplemented in the medium.

Strains construction

The construction of recombinant plasmids was performed using the 2 × MultiF Seamless Assembly Mix (ABclonal, RK21020). To label the Pdel gene, pdel, along with its native promoter, was fused with either gfp or bfp and cloned into the pBAD backbone. For gene overexpression, the target genes of interest were cloned into either a p15A ori plasmid or a pUC ori plasmid, under the control of the Arabinose-induction system. For the c-di-GMP sensor (Snapgene: #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. coliMG1655 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 9h 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). This suspension was rotated at 4°C for 16 hours.

Cell permeabilization

1 ml of fixed cells were centrifuged at 5,000g for 5 minuts at 4°C, followed by resuspension in 1 ml washing buffer (100 mM Tris-HCl pH7.0, 0.02 U/μl SUPERase-In RNase Inhibitor, AM2696, Invitrogen). After another centrifugation at 5,000g for 5 minuts 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 minuts. 1 ml cold PBS-RI was added, and the cells were centrifuged at 5,000g for 5 minuts at 4°C, followed by resuspension in 250 μl Lysozyme Mix (Dissolving 250 μg/ml Lysozyme or 5 μg/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 cold PBS-RI (1ml) 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 H2O), 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 minuts. 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 7. 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 H2O, 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 10x T4 ligase buffer, 150 μl H2O. 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 H2O, 600 μL 0.5 M EDTA.

Round 1 RT reaction

About 3×107 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). 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 H2O, 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 0.5 μ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 H2O. 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 H2O. The depletion of r-cDNA could be repeated using the above protocol, and ultimately, the cDNA was eluted in 20 μl of H2O.

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 MgCl2 (0.1 M), and 12.8 μl H2O. 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 H2O. 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 described7 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 are two replicates of E coli static biofilm, these two datas need to be merged into one SeuratObject and removed batch effects. 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 datas of all samples. For preprocessing of E coli static 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 exponential period E coli data. 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 datas 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 E coli static biofilm into a single SeuratObject by FindIntegrationAnchors and IntegrateData functions. Then all the datas underwent scaling using the ScaleData function, followed by dimension reduction through principal component analysis. After principal component analysis, we removed batch effect by RunHarmony of the two replicates of E coli static biofilm datas. Then a graph-based clustering approach was employed in all datas 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 datas. Marker genes for each cluster were computed using the Wilcoxon Ranksum 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 function16 within R.

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). Converting sam files to bam files using samtool (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 log2 reads per gene of bulk RNA-seq and log2 UMI per gene of scRNA-seq.

Flow cytometry and FACS analysis

All samples were measured on a Beckman CytoFLEX SRT flow cytometer with a 70 μm nozzle in which normal saline was used as sheath fluid. PdeI-BFP or c-di-GMP sensor labeled strain in the 24 h static biofilm growth phase was washed and resuspended in sterile PBS. Microorganisms were identified by FSC (forward scatter) and SSC (side scatter) parameters. Cells were sorted into different groups based on their fluorescence intensity (PB450 for BFP, FITC for mVenusNB, ECD for mScarlet-I). The results were analyzed by FlowJo V10 software (Treestar, Inc.).

Antibiotic killing and persister counting assay

Cells sorted from FACS were suspended in fresh LB broth containing 150 μg/ml ampicillin and incubated for an additional 3 hours at 37°C with continuous shaking at 220 rpm. Following this, the cells were plated on LB plates for colony-forming unit (CFU) counts before the ampicillin challenge, and the plates were incubated overnight at 37°C. After the ampicillin challenge, cells were washed with PBS buffer and plated again for CFU counts. The persister ratio was defined as the ratio of CFU after ampicillin challenge to CFU before ampicillin challenge. Averages and standard deviations presented are representative of three biological replicates.

Microscopy

Bright-field and fluorescence imaging

Inverted microscopes (Nikon Eclipse Ti2 and Leica Stellaris 5 WLL) were utilized, with illumination provided by different lasers: 405 nm for BFP, 488 nm for GFP, respectively. The fluorescence emission signal was captured by an sCMOS camera (pco.edge 4.2 bi). Filter sets tailored to each fluorophores spectrum were employed. Image analysis was conducted using ImageJ software (Fiji). For c-di-GMP sensor analysis, the ratio of mVenusNB to mScarlet-I (R) exhibited a negative correlation with the concentration of c-di-GMP. Therefore, R-1 displayed a positive correlation with the concentration of c-di-GMP.

Time-lapse imaging

To observe antibiotic killing and bacteria resuscitation processes, cells labeled with Pdel-GFP in the 24-hour static growth phase were collected, washed twice with PBS buffer, and imaged on a gel pad containing 3% low melting temperature agarose in PBS [REF]. The gel pad was prepared in the center of the FCS3 chamber as a gel island. The cells were then observed under bright field or epifluorescence illumination. Subsequently, the gel pad was surrounded by LB broth containing 150 μg/ml ampicillin for 6 hours at 35°C. Fresh LB broth was flushed in, and the growth medium was refreshed every 3 hours, allowing cells to recover sufficiently.

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

Determination of c-di-GMP concentration by HPLC-MS/MS involved growing MG1655 Δara pBAD::pdeI and MG1655 Δara to mid-exponential growth phase, followed by induction with 0.002% arabinose. After 2 hours of incubation, cells were harvested and washed with PBS. The washed cells were rapidly frozen using liquid nitrogen. Simultaneously, another portion of washed cells was stained with SYTOTM24 and quantified using flow cytometry. c-di-GMP concentration detection was conducted by Wuhan Lixinheng Technology Co. Ltd. through high-pressure liquid chromatography-tandem mass spectrometry (HPLC-MS/MS). All strains were assayed in biological triplicates, and 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. Significance was determined by a two-tailed Student’s t-test. Error bars indicate the standard deviations of the mean from a minimum of three independent experiments. A significance threshold of P < 0.05 was applied. Asterisks were used to denote significant differences (*, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001).

Declarations

Acknowledgments

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), and Science Fund for Distinguished Young Scholars of Hubei Province (2022CFA077). 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.

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.

Additional Declarations

The authors declare no competing interests.

Supplementary figures

Technical Application of RiboD-PETRI in S. aureus and C. crescentus. (A-E) Single-cell RiboD-PETRI results for Staphylococcus aureus strain ATCC 25923 (S. aureus), cultured for 9 hours in MHB medium at 37°C. (F-J) Single-cell RiboD-PETRI results for Caulobacter crescentusstrain NA1000 (C. crescentus), incubated at 37°C for 3 hours. (A & F) Distribution of mRNA UMIs captured per cell in RiboD-PETRI data of(A) S. aureus and (F) C. crescentus, presented as violin plots showing the upper quartile, median, and lower quartile lines. “SA” denotes S. aureus, and “CC” denotes C. crescentus. For S. aureus, 1,000, 5,000, and 8,000 cells were selected in descending order of UMI counts per cell, with median UMIs per cell of 273,122, and 94, respectively. For C. crescentus, 1,000, 5,000, and 10,000 cells were selected with median UMIs per cell of 2190, 662, and 225, respectively. (B & G) Screening of single-cell data for (B) S. aureus and (G) C. crescentus, resulting in 9386 and 5728 cells, respectively. Screening criteria included UMIs between 15 and 1000 for S. aureus and UMIs between 200 and 5000 for C. crescentus, along with gene counts greater than 30. (C & H) Distribution of UMIs on the UMAP results for (C) S. aureus and (H) C. crescentus. UMAP results reveal heterogeneity among populations, with each point representing a cell and color shading indicating UMI counts. (D & I) UMAP visualization of (D) S. aureus and (I) C. crescentus, demonstrating the ability of RiboD-PETRI to distinguish population heterogeneity. (E & J) Expression of major marker genes for (E) S. aureus and (J) C. crescentus overlaid on the UMAP plot, highlighting cells with high expression levels in blue.

Technical performance of RiboD-PETRI. (A) Scatterplot illustrating the relationship between reads per cell and counts of UMIs per cell detected from exponential phase E. coli data. Each dot represents a cell. (B) Scatterplot demonstrating the relationship between reads per cell and counts of UMIs per cell detected from static biofilm E. coli data. Two replicates of the sample are included. (C) Calculation of the Pearson correlation coefficient (r) of UMI counts per gene between replicate 1 and replicate 2 of static biofilm E. coli. The analysis involved 4,062 out of 4,141 total genes, with a significant correlation (p-value < 0.0001, r = 0.96), indicating good replication between samples. Each dot represents a gene. (D & F) Normalized and Principal Component Analysis (PCA) performed on screened data of (D) exponential phase E. coli, resulting in 1,464 cells, and (F) two replicates of static biofilm E. coli, resulting in 1,621 and 3,999 cells, respectively. The resulting scatterplots show heterogeneity among the populations, with each point representing a cell. (E & G) UMAP plots generated based on the distribution of UMI counts of cells, illustrating the distribution of UMI counts for different cells in (E) and (G). The distribution of clusters in the UMAP results is independent of UMI counts. (H) UMAP plot based on the original identity of static biofilm E. co//samples (replicate 1 and replicate 2). Each dot represents a cell, with red indicating replicate 1 and green indicating replicate 2.

Marker Genes Identified in E.coli 3h Single-Cell Data by RiboD-PETRI. Expression levels of various marker genes during the 3-hour exponential period of E. coli overlaid on the UMAP plot. Cells with high expression levels are depicted in blue. Marker genes were selected based on a p-value greater than 0.001 and a log2FC greater than 3.

Marker Genes Identified in E. coli in 24-hour static biofilms by RiboD-PETRI. Expression levels of different marker genes during the 24-hour static biofilm stage of E. coli overlaid on the UMAP plot. Cells with high expression levels are highlighted in blue. Marker genes were selected based on a p-value greater than 0.001 and a log2FC greater than 3.