Authors: Ari Sarfatis (1Program in Cellular and Molecular Medicine, Boston Children’s Hospital; Boston, MA 02115 USA.; 2Department of Microbiology, Blavatnik Institute, Harvard Medical School; Boston, MA 02115 USA.), Yuanyou Wang (1Program in Cellular and Molecular Medicine, Boston Children’s Hospital; Boston, MA 02115 USA.; 2Department of Microbiology, Blavatnik Institute, Harvard Medical School; Boston, MA 02115 USA.), Nana Twumasi-Ankrah (1Program in Cellular and Molecular Medicine, Boston Children’s Hospital; Boston, MA 02115 USA.; 2Department of Microbiology, Blavatnik Institute, Harvard Medical School; Boston, MA 02115 USA.), Jeffrey R. Moffitt (1Program in Cellular and Molecular Medicine, Boston Children’s Hospital; Boston, MA 02115 USA.; 2Department of Microbiology, Blavatnik Institute, Harvard Medical School; Boston, MA 02115 USA.; 3Broad Institute of Harvard and MIT; Cambridge, MA 02142 USA.)
Categories: Article
Source: Science (New York, N.Y.)
Authors: Ari Sarfatis, Yuanyou Wang, Nana Twumasi-Ankrah, Jeffrey R. Moffitt
Single-cell decisions made in complex environments underlie many bacterial phenomena. Image-based transcriptomics approaches offer an avenue to study such behaviors, yet these approaches have been hindered by the massive density of bacterial mRNA. To overcome this challenge, we combine 1000-fold volumetric expansion with multiplexed error robust fluorescence in situ hybridization (MERFISH) to create bacterial-MERFISH. This method enables high-throughput, spatially resolved profiling of thousands of operons within individual bacteria. Using bacterial-MERFISH, we dissect the response of E. coli to carbon starvation, systematically map subcellular RNA organization, and chart the adaptation of a gut commensal B. thetaiotaomicron to micron-scale niches in the mammalian colon. We envision bacterial-MERFISH will be broadly applicable to the study of bacterial single-cell heterogeneity in diverse, spatially structured, and native environments.
Population-level bacterial dynamics often emerge from the heterogeneous behaviors of single cells. Notable examples include entry into and exit from antibiotic persistent states (1), bet hedging (2), virulence factor expression (3), and cellular specialization within biofilms (4). Recent advances in bacterial single-cell RNA-sequencing (scRNA-seq) offer an exciting avenue to study such phenomena by providing transcriptome-wide expression profiles for thousands of cells (5–14). Indeed, such methods have provided insights into antibiotic response (11, 13, 14), prophage activation (9, 11, 14), toxin expression (9, 10), sporulation (10), competence (9, 10), mobile genetic elements (11, 14), cell-cycle-dependent gene regulation (15), and functional heterogeneity within the rumen microbiome (16).
Missing from these studies is the natural spatial context in which many behaviors occur. Yet, spatial organization is an essential modulator of bacterial dynamics across a range of length scales. On the tens-of-micron-scale, spatial gradients in small molecule concentrations tune bacterial responses, define niches for commensal growth, or shape interactions in multi-species communities (17, 18). On the micron-scale, direct cell-to-cell contact controls effector protein delivery which mediates predation (19), self- versus non-self recognition (20), contact-dependent inhibition (21), and virulence (22). Finally, even sub-micron length scales are relevant, as growing evidence indicates that the bacterial transcriptome is internally organized with functional consequences (23). Unfortunately, such spatial information is lost during cell dissociation and RNA extraction in scRNA-seq; thus, current methods are not well suited for the study of such processes.
By contrast, image-based approaches to single-cell transcriptomics provide this spatial context by directly imaging and identifying RNAs within fixed cells in their native spatial environment (24–30). Moreover, by leveraging combinatorial optical barcodes to distinguish RNAs, these measurements can be massively multiplexed, producing spatially resolved, transcriptome-scale expression profiles that span intracellular to tissue-scale lengths (24–26). In eukaryotic systems such methods have mapped the intracellular RNA organization, explored regulatory networks, and defined, discovered, and charted cell types and states across a range of tissues (24–26). Unfortunately, current methods are not compatible with bacteria, as the massive density of bacterial RNA challenges combinatorial barcode detection. Non-combinatorial barcoding approaches can bypass this challenge, as recently illustrated with par-seqFISH (31). This method profiled the expression of 105 genes in sessile and planktonic Pseudomonas aeruginosa, revealing, among many features, a diversity of distinct cellular states in cell culture and patches of coordinated gene expression associated with distinct anaerobic processes in biofilms (31). This work highlights the potential of image-based approaches for single-cell transcriptomics in the study of bacteria. Yet, there are biological questions that would benefit from higher multiplexing, which would be challenging to accomplish with the non-combinatorial barcoding approach of this method.
Here we overcome this RNA density challenge and introduce a transcriptome-scale, image-based approach for bacterial, single-cell transcriptomics. This approach combines an expansion microscopy toolbox optimized for bacteria with multiplexed error-robust fluorescent in situ hybridization (MERFISH) (30) and allows single-cell profiling of up to 80% of the transcriptome. We demonstrate that this technique—bacterial-MERFISH—accurately profiles 97, 1,057, or 1,930 operons with large detection efficiency, accuracy, and throughput in log-phase Escherichia coli (E. coli) cells. To highlight the discovery potential of this technique, we first profile E. coli response to a carbon-source switch, revealing a heterogeneous sequential nutrient exploration program. Next, we chart the intracellular organization of the E. coli transcriptome, uncovering a previously unappreciated diversity in spatial patterning and a cooperative role for genome and proteome organization in shaping transcriptome organization. Finally, we map the adaptation of a human gut commensal—Bacteroides thetaiotaomicron (B. theta)—to the mouse colon, revealing micron-scale fine-tuning of gene expression based on local polysaccharide availability. More broadly, these measurements illustrate the potential for bacterial-MERFISH to reveal single-cell heterogeneity in a wide range of biological and spatial contexts.
MERFISH enables the identification of thousands of different mRNA molecules by using combinatorial, error-robust, fluorescent optical barcodes built from repetitive rounds of single-molecule FISH (smFISH) (30). However, to decipher barcodes, the fluorescent signal from different molecules must be optically resolvable. For conventional high-resolution optical microscopy, only a few molecules per μm^3^ can be distinguished (32). For eukaryotic systems, large transcriptome fractions can be targeted while satisfying this limit (30). By contrast, a log-phase E. coli cell contains ~8,000 mRNA molecules in a cell volume of ~3 μm^3^ (33), producing a total mRNA density nearly three orders of magnitude greater than that resolvable with diffraction-limited imaging (Fig. 1A). Thus, mRNA density restricts the imaging of more than a small number of bacterial mRNAs (34) and is a substantial challenge to transcriptome-scale imaging.
Density reduction is a natural avenue to address this challenge. Indeed, image-based approaches to single-cell transcriptomics have achieved modest degrees of RNA density reduction by spreading the barcode signal over more imaging rounds (35, 36) or by leveraging expansion microscopy (37) to physically swell the sample (35, 38). While the approximately 10-fold reduction in RNA density achieved by these approaches was sufficient to extend image-based transcriptomes to whole-transcriptome-scale in eukaryotes (35, 36), such density reduction is still nearly two orders of magnitude insufficient for similar profiling in bacteria.
To overcome bacterial mRNA density, we leveraged recent advances in expansion microscopy (39, 40) to develop a bacterial-FISH-optimized expansion toolbox capable of up to 1000-fold volumetric expansion (Fig. 1B), complementing recent bacterial expansion methods developed for non-RNA targets and with modest degrees of expansion (41–46). Briefly, we grew E. coli to mid-log phase, fixed them with paraformaldehyde (PFA), digested the cell wall, expanded them in a Ten-fold Robust Expansion (TREx) gel (39), and then re-embedded the sample in a non-expanding, stabilizing gel (Fig. 1B; Materials and Methods). Using custom expansion-optimized staining protocols, we labeled samples with a 16S ribosomal RNA (rRNA) probe and a MERFISH probe set targeting 97 E. coli operons (tables S1 and S2; Materials and Methods). We targeted operons rather than individual genes as the signals from different barcodes from genes on the same polycistronic mRNA would not be identifiable due to the overlapping signals corrupting these barcodes.
Expansion increased the width of E. coli 3.7±0.8-fold (standard deviation [STD], n=4; fig. S1A) with a modest cell-to-cell variation in width (fig. S1B). As this linear expansion is consistent with a ~50-fold volumetric expansion (fig. S1C), we term this approach the 50X expansion protocol. In expanded cells stained with the 97-operon library, individual fluorescent puncta were visible, and these molecules were identified with MERFISH (Fig. 1, C to H), indicating sufficient expansion to resolve the mRNA density of this targeted library. We noted an ample number of molecules even in the absence of specific RNA-gel anchoring chemistries in PFA-fixed but not methanol-fixed expanded samples, suggesting a PFA-dependent, RNA-anchoring mechanism (fig. S1, D to J).
To further explore the multiplexing possible with the 50X protocol, we stained 50X-expanded E. coli with a MERFISH probe set targeting 1,057 operons, roughly 40% of the E. coli transcriptome (tables S1 and S2; Materials and Methods). While the mRNA density was higher than that observed for the 97-operon measurement, individual puncta were still resolved and many molecules were identifiable with MERFISH (fig. S1, K to M). Nonetheless, we noticed an increased frequency of overlapping RNA signals (fig. S1L). To address this overlap, we developed an iterative expansion protocol that combined TREx (39) with a previous iterative strategy (40). Briefly, cells expanded and stabilized with the 50X protocol were expanded in a second TREx gel and then embedded in a second stabilizing gel (Fig. 1B; Materials and Methods). This protocol expanded the width of cells 11.1±2.2-fold (STD, *n=*4; fig. S1A) with modest variation in the expanded width from cell to cell (fig. S1B). As this linear expansion is consistent with a ~1,400-fold volumetric expansion (fig. S1C), we term this protocol the 1000X protocol. When cells were expanded with this protocol, the signal from individual RNAs and their identity were clearly distinguished when 1,057 operons were stained (Fig. 1, I to K). Notably, the 1000X protocol starts with 50X-expanded samples (Fig. 1B), facilitating exploration of the necessary expansion for a given sample.
Inspired by the ability to expand E. coli volumetrically by three orders of magnitude, we designed a MERFISH library that covers 80% of the transcriptome, corresponding to 1,930 operons (tables S1 and S2; Materials and Methods). 1000X-expanded samples stained with this probe set showed clear single-molecule signals with limited overlap, and these RNAs could be identified with MERFISH (fig. S1, N to P), suggesting that 1000X expansion is sufficient to allow MERFISH profiling of a substantial fraction of the E. coli transcriptome.
Notably, during the early development of bacterial-MERFISH, we observed a low correlation between the abundance determined via bacterial-MERFISH and that of bulk RNA-sequencing for lowly expressed operons (Fig. 1L). This loss of correlation occurred at a level much greater than the internal measurements of false-positive rates provided by the blank controls, which are barcodes for which no probes were assigned (Fig. 1L; Materials and Methods). This observation suggested a bacterial-expansion-dependent source of false positives, which we reasoned might be due to probe binding to genomic DNA melted during expansion. Supporting this hypothesis, DNase treatment reduced this apparent false-positive rate, bringing it into agreement with that measured with internal controls (Fig. 1M). Thus, we conclude that off-target binding to DNA is not a dominant source of false positives after DNase treatment.
To benchmark the performance of bacterial-MERFISH, we performed two replicate measurements of 97, 1,057, or 1,930 operons in combination with 50X or 1000X expansion in log-phase E. coli, segmented cells from these images, and partitioned RNAs into those cells (Fig. 1, F to K, and fig. S1, K to P). We observed strong correlation between these measurements and bulk RNA-sequencing across all conditions and for mRNAs expressed at multi-copy- to much less than single-copy-per-cell levels (Fig. 1, M to O, and fig. S2, A to C), indicating that bacterial-MERFISH can accurately profile RNA expression across four orders of magnitude in abundance with a lower limit set by the measured false-positive rate. Further supporting its accuracy, measurements with bacterial-MERFISH correlated strongly between biological replicates, expansion protocols, and multiplexing levels (fig. S2, D to F).
To determine the sensitivity of bacterial-MERFISH, we next calculated the detection efficiency—the fraction of targeted molecules actually detected (Fig. 1P and fig. S2G; Materials and Methods). Values ranged from 7% for near-whole transcriptome measurements (1,930 operons) to as high as 50% for more targeted measurements (97 operons), with additional expansion providing an increase from 10% to 25% for 1,057 operons measurements (Fig. 1P). As expected, the number of mRNA counts and unique operons observed per cell varied based on the multiplexing and detection efficiency (Fig. 1, Q and R). Finally, bacterial-MERFISH can image large numbers of cells, comparable or greater than those characterized previously with scRNA-seq (5–16), despite the decreased imaging throughput due to expansion (fig. S2H). With the ability to vary multiplexing and expansion, it is possible to balance different aspects of performance, e.g., the number of imaged cells versus the degree of expansion, to best suit the question.
While bacterial-MERFISH measurements correlate strongly with reported scRNA-seq measurements in E. coli in comparable media, protocol differences complicate a quantitative comparison of the sensitivity of the methods (fig. S3). Nonetheless, the detection efficiency of bacterial-MERFISH (Fig. 1P) compares favorably to reported values for single-cell RNA sequencing (8–10). Similarly, the total number of detected molecules per cell (Fig. 1R) can be comparable to that detected via sequencing methods, even though these methods have the capacity to detect all transcripts while bacterial-MERFISH is targeted. Thus, these measurements collectively indicate that bacterial-MERFISH is a high-performance, versatile, image-based approach to bacterial single-cell transcriptomics that complements scRNA-seq methods with its low false-positive rates and large dynamic range, detection efficiencies, and throughput.
One promise of single-cell bacterial transcriptomics is the ability to identify heterogeneous bacterial behaviors and to computationally resynchronize the desynchronized response of individual cells to environmental stimuli, unmasking dynamics obscured by population averages. One illustrative dynamic response is the diauxic shift caused by a switch in carbon source. Bacteria often consume different carbon sources in a preferential order regulated, in part, by carbon-catabolite repression (CCR) (47). A bacterial population grown in a mixture of two carbon sources will consume the preferred source, pause growth in a diauxic shift, and then resume growth on the less preferred substrate (48). The diauxic shift is classically interpreted as the time required to express utilization machinery for the second sugar; however, recent single-cell studies suggest that this response is instead shaped by differential dynamics of sub-populations (49). Nonetheless, the diversity and transcriptional profiles of such sub-populations remain poorly defined.
Thus, we revisited the classic diauxic shift with bacterial-MERFISH. We grew E. coli in a defined minimal medium with a mixture of glucose and xylose (Materials and Methods). As expected, the culture grew logarithmically until glucose was exhausted, paused growth in a diauxic shift, and then resumed logarithmic growth on xylose before entering stationary phase once both sugars were exhausted (Fig. 2A). We harvested E. coli cells throughout this process and profiled the expression of 1,057 operons (Fig. 2, A to C) in 296,666 cells with 50X expansion—a multiplexing and expansion level set to balance transcriptome-wide profiling with the number of profiled cells. For experimental efficiency, multiple time points were combined and profiled in a single MERFISH measurement by labeling cells with barcoded 16S rRNA probes (31) optimized for expansion protocols (Fig. 2B; Materials and Methods). Putative cells were segmented, RNA molecules assigned to these cells, and conditions demultiplexed with the 16S rRNA barcode (Fig. 2C; Materials and Methods).
Supporting these measurements, we found that the average RNA expression determined via MERFISH was strongly correlated between biological replicates, with bulk RNA-sequencing at different growth stages, and with the average abundance determined by scRNA-seq in comparable conditions (fig. S4, A to E). Moreover, cell populations displayed expected expression patterns (Fig. 2D) with cells harvested during growth expressing operons associated with amino acid synthesis (e.g., argG and ilvC), translation (e.g., tff-rpsB-tsf and cmk-rpsA-ihfB), and aerobic respiration (e.g., atpIBEFHAGDC), and cells harvested during the diauxic shift and stationary phase expressing operons associated with stress response (e.g., nlpD-rpoS) and gluconeogenesis (e.g., glpFKX and glpD). Finally, despite a difference in the effective detection efficiency between replicates, both replicates showed the same trend in the counts per cell with greater counts per cell during growth versus non-growth conditions (fig. S4, F and G). Similar trends have been seen with other single-cell methods during bacterial growth in other media (9, 31).
To explore the heterogeneity in cellular response to this carbon shift, we integrated the measurements from two biological replicates, visualized single-cell heterogeneity with Uniform Manifold Approximation and Projection (UMAP), and performed Leiden clustering to distinguish sub-populations (Fig. 2, E to I; Materials and Methods). Supporting our analysis, cells harvested at similar growth phases largely co-integrated (fig. S4, H and I). Nonetheless, even modest transcriptional differences, such as those between cells collected from log-phase growth in glucose or xylose, were resolved (Fig. 2E and fig. S4J). This analysis revealed a rich diversity in behavior at the single-cell level. Cells were organized into two major groups, corresponding to cells taken from conditions of growth or non-growth (Fig. 2E), and were collectively sub-divided into 14 different clusters (Fig. 2F). These clusters had unique gene expression profiles (Fig. 2, G and H) and abundances across different conditions (Fig. 2I). Supporting cluster validity, individual clusters were observed across both replicates (Fig. 2I), clusters were each marked by multiple operons associated with related biological processes (Fig. 2H), marker expression was conserved between replicates (fig. S5A), the operon covariances that drive the emergence of distinct clusters were conserved between replicates (fig. S5, B to E), and clusters co-expressing similar operon sets to those seen when all cells were jointly analyzed were identified when each condition was analyzed alone (fig. S6). Finally, single-molecule FISH (smFISH) in unexpanded E. coli further supported this analysis by confirming the co-expression of cluster markers (fig. S7).
During log-phase growth, we identified a diversity of sub-populations (Fig. 2, E to I), including clusters associated with translation and amino acid transport (G0 marked by tff-rpsB-tsf and lysP), nucleobase synthesis (G6 marked by carAB, codBA, and purHD), arginine synthesis (G7 marked by argG, argD, and argCBH), serine biogenesis (G4 marked by serA and serC-aroA), sulfate utilization (G11 marked by cysDNC and cysJIH), protein folding (G12 marked by dnaKJ, groSL, hslVU, clpB, and htpG), and motility and chemotaxis (G3 marked by motRAB-cheAW, fliDST, and tar-tap-cheRBYZ) (Fig. 2H). Many of these biochemical processes are required for growth in minimal media, yet our analysis revealed that only subsets of cells expressed high levels of these essential operons, suggesting a model in which the homeostatic, population-level expression of these pathways is produced not by uniform expression across all cells but rather by transient, coherent bursts of expression that time-average to required levels. This observation is consistent with transient bursts in promoter activity revealed with live-cell imaging (50, 51). Moreover, this observation is also supported by recent scRNA-seq measurements performed in E. coli in comparable defined medium (10), which also highlighted clusters of co-expressing genes associated with similar cellular processes (e.g., cell motility [G3 marked by motRAB-cheAW, flgBCDEFGHIJ, and tar-tap-cheRBYZ], and carbamoyl metabolism [carAB] feeding into arginine synthesis [G7 marked by argG, argD, and argCBH] and nucleobase synthesis [G6 marked by purHD and lapAB-pyrF-yciH]). Collectively, these measurements now suggest that such bursts are widespread and occur not at the individual operon level but rather upstream in common regulatory factors.
We also observed a similar degree of heterogeneity during the diauxic shift (Fig. 2, E to I). Only a subset of clusters expressed operons associated with xylose utilization (N2 and N10 both marked by xylAB, xylE, xylFGHR, and yagGH) while most clusters in the diauxic shift expressed utilization operons associated with carbon sources not included in our medium (Fig. 2, E to I). These sources comprised mannose and glycerol (N1 marked by manXYZ, glpABC, and glpTQ); maltose (N9 marked by malEFGH, malKM-lamB, and malPQ); arabinose (N10 marked by araE-ygeA, araFGH, and araBAD in addition to xylose-associated operons); mannitol and glycolate (N8 marked by mtlADR and glcDEFGBA); and acetate and putrescine (N5 marked by acs-yjcH-actP, puuAP-ymjE, and puuDRCBE) (Fig. 2, H and I). N13, like G3, was marked by chemotaxis and motility operons. However, these clusters were distinguished by operons associated with chemotaxis towards peptides (G3; e.g., tar-tap-cheRBYZ and tsr) or sugars (N13; e.g., trg and aer), reflecting the differential nutrient needs between the conditions associated with these clusters (Fig. 2H).
We also noted an apparent temporal order in which specific carbon sources were explored. To investigate this ordering, we computationally resynchronized cells using a pseudotime analysis from glucose-log-phase growth through the diauxic shift (Materials and Methods). Not only did this analysis reproduce the order of sampling during the diauxic shift (Fig. 2, J and K), supporting the pseudotime ordering, it also revealed a temporal cascade of carbon-utilization operon expression (Fig. 2L). This cascade started with operons associated with glucose (e.g., ptsG) and then proceeded through operons associated with glycogen (e.g., segA-pgm), peptidoglycan recycling (e.g., folP-glmM), and TCA cycle entry (e.g., pdhR-aceEF-lpd), consistent with the mobilization of carbon source reserves. The cascade then transitioned from mostly non-glucose phosphotransferase system (PTS) carbohydrate utilization (e.g., manXYZ, nagE, gatYZABCD, and mglBAC for mannose, N-acetylglucosamine, galactitol, and galactose, respectively) to mostly non-PTS carbon sources (e.g., malEFGH, ugpBAECQ, and glcDEFGBA for maltose, glycerol 3-phosphate, and glycolate)—consistent with the view that, after glucose, non-glucose PTS carbohydrates are preferred over non-PTS carbohydrates (52). This progression culminated with the upregulation of xylose (e.g., xylFGHR and xylAB) and, finally, arabinose (e.g., araFGH and araBAD) utilization operons. The transient, co-expression of xylose- and arabinose-utilization operons seen at late pseudotime values (Fig. 2L) may represent an overshoot in carbon-source progression due to the lag between transcription of xylose-utilization operons and the transition to the functional utilization of xylose.
Carbon utilization operons were not the only pathways differentially expressed throughout this pseudotime progression. Extending this analysis to any expressed operon revealed a rich remodeling of the transcriptome throughout the diauxic shift (fig. S8 and table S3). Notably, this analysis further supported the pseudotime ordering as stationary phase markers (e.g., katE, aldB, and puuDRCBE), presumably indicative of more advanced stress and starvation, were most enriched at late pseudotime values. More broadly, this analysis revealed that the transcriptional progression observed in carbon metabolism should be understood as part of a broader transcriptional response.
The range of pseudotime values observed for each shift condition overlapped with those of other conditions (Fig. 2M), consistent with a stochastic response model in which cells are desynchronized in their progression along this carbon-source hierarchy. The simplest model that explains this hierarchy is one in which each cell will eventually explore this full range of carbon sources; however, it is worth noting that, as a static snapshot in time, our data do not rule out a model in which a subset of carbon sources is only explored by a subset of cells.
Together these results reveal that the homeostatic levels of required pathways can be maintained by transient, coherent bursts of expression within entire regulons, and that, when faced with carbon starvation, E. coli adopts a responsive diversification strategy (53) in which the lack of glucose triggers a stochastic, hierarchical progression along carbon utilization operons. We suggest that, with the ability to characterize sub-populations and computationally resynchronize cells undergoing dynamic responses, bacterial-MERFISH could prove useful in dissecting the diverse regulatory mechanisms that encode such cellular heterogeneity.
Another advantage of image-based, single-cell transcriptomic methods is the ability to measure the intracellular transcriptome organization. Specific bacterial mRNAs are localized to the cytoplasm (54–56), membrane (54–56), poles (54, 56–59), chromosomal loci of origin (60), or the surface of intracellular organelles (61) as determined via low throughput mRNA imaging (54–61) or biochemical fractionation (54, 56). This organization has proposed functional roles in protein sorting, complex assembly, and mRNA turnover (23). However, as RNA localization has yet to be imaged at the transcriptome-scale, the extent and diversity of spatial organization remain unclear.
To explore the spatial organization of the E. coli transcriptome, we leveraged our 50X-expanded, 1,057-operon measurements of E. coli grown in LB. To define intracellular RNA localization, we mapped each mRNA molecule to its axial and radial position within each cell and normalized these coordinates to the cell length and width (Fig. 3A; Materials and Methods). We then computed the average axial and radial distributions for each mRNA across all measured cells (Fig. 3, B to D; Materials and Methods). Consistent with previous reports in E. coli (23), we identified mRNAs enriched in the cytosol (e.g., dnaKJ), at the membrane (e.g., ptsG), and toward the poles (e.g., dnaB) (Fig. 3, C and D). However, we also noticed a diversity of variations on these major patterns, with mRNAs enriched at different locations along the membrane, at multiple cytoplasmic locations, or in a single central focus (Fig. 3E). These patterns were reproduced between replicate measurements (fig. S9, A to C) and with 1000X-expanded MERFISH measurements (fig. S9D). Moreover, as sample expansion might introduce some degree of distortion (39, 40), we additionally validated a small number of operons with unexpanded smFISH and observed highly similar localization patterns (Fig. 3F).
To further categorize this spatial diversity, we leveraged a measure of pattern similarity to visualize and cluster spatial patterns (Materials and Methods). This analysis produced five major clusters of mRNA distributions (Fig. 3G) which showed some degree of continuous spatial variation within the clusters. The Cytoplasmic cluster comprised patterns including uniform filling of the cytoplasm (e.g., dnaKJ and glnS), multiple foci throughout the cell (e.g., uxuAB and groSL), or increased central density (e.g., metG). mRNAs with a strong central focus (e.g., adhE and rnb) defined the Midcell-cytosolic cluster. The Membrane cluster was defined by diffuse membrane-enriched patterns (e.g., ptsG and cydAB) or membrane enrichment at or adjacent to the poles (e.g., yifK, kgtP or yojI) while the Midcell-membrane cluster was defined by strong membrane enrichment only in the middle of the cell (e.g., dtpA and pntAB). Finally, the Polarized cluster was defined by diffuse (e.g., mreBCD) or sharp cytoplasmic foci (e.g., metBL) near the poles. Some mRNAs at cluster boundaries shared features similar to nearby clusters (e.g., proS and proP), underscoring a continuous variation in spatial patterns.
To investigate possible patterning mechanisms, we explored the correlation between clusters and mRNA features. Covariation with spatial distribution was modest for transcript length, GC content, abundance, and half-life (fig. S10, A to H) but strong for encoded protein location (Fig. 3, H to J, and fig. S10, I to N). Specifically, operons containing mRNAs that encode at least one inner-membrane protein were enriched in membrane-associated RNA localization clusters, whereas mRNAs that encode cytoplasmic proteins were enriched in clusters found within the cytosol (Fig. 3J). These observations are consistent with previous measurements for individual mRNAs (54) as well as transcriptome-wide mRNA groups (55); with the co-translational insertion mechanism of inner-membrane proteins, which would concentrate mRNAs at the membrane during translation (62); and with the report of membrane-associated sequence features enriched in inner-membrane-protein-coding mRNAs (63, 64). To test the extent to which co-translational insertion is responsible for membrane localization, we used the presence of an annotated transmembrane domain in inner membrane proteins as a proxy for co-translational insertion. Of the 213 RNAs in the Membrane or Midcell-membrane clusters, we found 143 inner-membrane-protein-encoding RNAs with at least one transmembrane annotation (fig. S10O), as opposed to only 33 of all 297 cytosol-enriched RNAs, further supporting the role of co-translational insertion in membrane localization.
In parallel, the genomic locus from which the mRNA was transcribed also correlated strongly with spatial patterning. mRNAs within the Midcell-cytosolic or Midcell-membrane clusters were preferentially encoded from genomic loci near the terminus of replication (terC) while the genes for mRNAs in other clusters were depleted in this chromosomal region (Fig. 3K and fig. S10P). During conditions of fast growth, the E. coli chromosome is organized with terC in the cell center and oriC and the left and right chromosomal arms replicating near the poles. This structure is maintained for the majority of the cell cycle, with the polar terC macrodomain formed at the pole of a newly divided cell moving rapidly to the cell center (65). Indeed, the axial distribution of mRNAs averaged across all (Fig. 3, L and M) or portions (fig. S10Q) of the cell cycle were consistent with this organization.
Our measurements now unify two previous models suggested for bacterial transcriptome organization—that either genomic (60) or proteomic (54–56) features dictate organization. Specifically, we show that both features play a role in the global organization of mRNAs, with protein location shaping cytoplasmic versus membrane enrichment and genomic feature shaping axial enrichment both within the cytoplasm or on the membrane. Nonetheless, we identified multiple mRNAs with spatial patterns that deviate from these global rules, including membrane enriched mRNAs that do not encode known inner-membrane proteins (fig. S11, A and B), cytoplasmic-enriched mRNAs that encode inner-membrane proteins with annotated transmembrane domains (fig. S11, C and D), and many mRNAs with spatial patterns inconsistent with those predicted for their genomic loci (fig. S11, E to G). These exceptions raise the possibility that there are mRNA-specific localization mechanisms that remain to be discovered and that these exceptions, in addition to global patterns, might have functional significance. Importantly, bacterial-MERFISH offers a direct approach to measuring such patterns, which should greatly enable mechanistic and functional studies of the intracellular organization of the bacterial transcriptome.
Image-based approaches to single-cell transcriptomics also promise the ability to explore gene expression within complex, spatially structured environments. To explore this capability of bacterial-MERFISH, we leveraged germ-free mice monocolonized with the human gut commensal B. theta (Fig. 4A). As Bacteroides have the ability to harvest a remarkable diversity of polysaccharides—both from the rich dietary pool as well as those deposited onto the mucus layer by the host (66)—we reasoned that B. theta might modulate its polysaccharide utilization based on local polysaccharide availability, as suggested previously (67). To explore this possibility, we designed a MERFISH panel that covers the diverse importers (i.e., SusC proteins) within polysaccharide utilization loci (PUL), hybrid two-component systems (HTCS) often associated with the regulation of PUL expression, and a handful of genes associated with central metabolism and other bacterial functions, targeting 159 operons in total (Fig. 4B, and tables S1, S2, and S4).
We then resected the colon from a monocolonized mouse, and the sample was methacarn-fixed, paraffin-embedded, and sectioned (Fig. 4A; Materials and Methods). In an unexpanded section, 16S rRNA staining revealed depletion of B. theta from the sterile inner mucus layer next to the host epithelial layer, enrichment near the outer mucus layer (68), and variable density throughout the colonic lumen (Fig. 4, A and C). We then expanded a section with a modified 50X-protocol and stained it for MERFISH. The 16S rRNA signal revealed that the distribution of B. theta was largely preserved despite a notable loss in dietary debris (Fig. 4, C and D; Materials and Methods). Individual mRNAs could be identified by MERFISH in single B. theta cells (Fig. 4, E and F). Across this colonic cross section, we observed that the expression of many B. theta operons was spatially variable. Some mRNAs were expressed throughout the lumen (e.g., BT_0364 and BT_4671) while others were more prominently expressed near the mucus layer (e.g., BT_3958 and BT_2894) (Fig. 4G). Multi-color smFISH in unexpanded samples supported these expression patterns (fig. S12).
To quantify this spatial variation, we performed these measurements in six slices sampled from the colon of two mice. Each colon was fixed with one of two different methods to control for potential fixation-dependent artifacts (Materials and Methods), which were minor as evidenced by strong correlation between these measurements (fig. S13, A to C). As PUL expression can be low (fig. S12), we averaged gene expression over small spatial patches containing ~5 B. theta cells and visualized transcriptional variation across patches using UMAP and diffusion maps (Fig. 4H and fig. S13, D to F). Individual genes expressed more prominently near the mucus layer or throughout the lumen defined different regions of this UMAP (fig. S13, D and E), and the diffusion analysis suggested that an important axis of gene expression variation, captured by the first diffusion component (DC1), is related to mucus proximity. Patches of low DC1 values were found enriched near the mucus whereas patches of high DC1 values were distributed more uniformly throughout the lumen (Fig. 4, H to K, and fig. S13F), revealing mucus proximity as a covariate of gene expression. More broadly, our data suggest that other spatial covariates remain to be described, as some samples showed luminal regions with unique patterns of gene expression (fig. S13G) or low-DC1 expression signatures throughout the entire lumen (fig. S13, H to M). Nonetheless, while these observations reveal that mucus proximity is a statistically significant covariate of B. theta expression variation for most samples (fig. S13, H to M), the presence of low-DC1 patches deeper into the lumen and many high-DC1 patches near the mucus suggest that local spatial heterogeneity may produce similar scale variations (fig. S13N) and that host proximity is not the only determinant of B. theta expression.
To determine the operons that underlie the mucus-proximal adaptation of B. theta identified with DC1, we separated spatial patches into low-DC1 (mucus-associated) or high-DC1 (lumen-associated) compartments. We found 38 operons with statistically significant enrichment between these two compartments (Fig. 4L and table S4), with substantial overlap in the enriched operons when samples from the two fixation methods were analyzed separately (fig. S13, O to R). There was a clear difference in the basic gene categories differentially expressed between the mucus- or lumen-associated compartments. 9 of the 15 lumen-associated operons were connected with central metabolism while none of the 23 mucus-associated operons had this functional annotation (Fig. 4L and table S4), suggesting that luminal B. theta may upregulate elements of central metabolism relative to B. theta closer to the mucus layer. By contrast, all 23 mucus-associated operons contain SusC or PUL-associated HTCS genes, as opposed to only 5 of the 15 lumen-associated operons, suggesting mucus-associated niches support harvesting of a greater polysaccharide diversity.
We next examined the known substrates for the PUL enriched in each compartment. Of the mucus-associated PUL with known substrates, 10 target host-mucus polysaccharides while only 1 targets dietary polysaccharides (Fig. 4L and table S4). By contrast, in the luminal compartment, 2 of the 3 PUL with known substrates target dietary polysaccharides (Fig. 4L and table S4). These observations are consistent with a greater availability of host-derived polysaccharides near the mucus, as expected, supporting our spatial analysis. Moreover, of the 23 PUL enriched in mucus-associated patches, 10 have unknown substrates (Fig. 4L and table S4). Given that we observe these operons enriched in mucus-associated patches, we predict that these PUL target host-derived polysaccharides or other carbohydrates enriched in this local environment.
Our measurements complement a previous microdissection study (67) of the spatial variation of B. theta gene expression by providing a direct micron-scale measure of this variation and by extending the list of mucus-enriched PUL due, perhaps, to increased spatial resolution, improved sensitivity, or biological variability between the studies. Further underscoring the importance of the improved spatial resolution offered by bacterial-MERFISH, we found that a simple categorization of patches based solely on host proximity was underpowered for the detection of differentially expressed operons with the sample number we profiled despite revealing luminal and mucus-proximal enrichment trends largely consistent with that observed above (fig. S14).
In total, these measurements reveal that bacteria fine tune gene expression to adapt to micron-scale niches in the gut. With bacterial-MERFISH it should now be possible to profile such micron-scale adaptation to a wide variety of complex environments.
Here we introduced bacterial-MERFISH, an image-based approach to single-cell transcriptomics that overcomes the massive mRNA density within bacterial cells by combining an optimized expansion microscopy toolbox with MERFISH. Bacterial-MERFISH offers complementary benefits to the growing suite of bacterial scRNA-seq methods (5–14). As a targeted method, it can sidestep abundant RNA (e.g., rRNA) challenges while providing an opportunity for high detection efficiency, which may prove essential for the many bacterial mRNAs that are very lowly expressed. In parallel, bacterial-MERFISH can be highly multiplexed, providing the ability to screen transcriptional changes with minimal prior knowledge of relevant targets. Bacterial-MERFISH can also image large numbers of cells, which may prove useful in the characterization of rare phenotypes. As an image-based technique, it naturally links gene expression to cell morphology or to intracellular molecular organization. Finally, with recent advances in all-optical readouts of pooled genetic screens (69, 70), it may now be possible to combine transcriptome-wide profiling with genome-wide perturbations in bacteria.
However, we anticipate that one of the most substantial advantages of bacterial-MERFISH will be the ability to profile bacterial behaviors in situ. Whether it is interactions between specialized cellular states within single-species biofilms or between different species in mixed communities, bacterial-MERFISH provides a means of linking spatial proximity, cellular micro-environment, and global architecture to gene expression. Studying bacteria in their native environment also bypasses the need for culture; thus, bacterial-MERFISH may offer an avenue for the in situ characterization of the diverse range of unculturable bacteria. Finally, as MERFISH can now target both eukaryotic and prokaryotic mRNAs, bacterial-MERFISH may allow the simultaneous profiling of host and bacterial gene expression, which may deepen our understanding of host-microbe interactions such as commensal colonization or pathogenic infection. More broadly, the ability to directly profile single-cell bacterial transcriptomes in their native, complex environments may offer a new window into the substantial range of bacterial behaviors not well captured in a culture flask.
To design a probe set that targets 101 operons in E. coli, we randomly selected operons to cover a range of expected abundance using published bulk RNA-sequencing data (55). MERFISH encoding probes are comprised of a target region which is complementary to the RNA of interest and flanking readout sequences. There is one readout sequence associated with each bit in the barcode assigned to targeted RNAs (30). To select the target regions for this library, we used a previously described pipeline (71) with the following 30-nt length; a melting temperature between 65 °C and 75 °C; a GC content between 40% and 74%; a gene specificity between 0.75 and 1; no predicted homology with rRNA; and no allowed overlap between target regions. For this initial library we chose to design target regions that would bind to a single gene within polycistronic operons. However, in a few instances, genes were too short, and we extended target region selection to the entire predicted operon. No final target had less than 54 target regions. The target region design pipeline is available at https://github.com/ZhuangLab/MERFISH_analysis.
To encode target identity, we assigned to each target a unique 16-bit binary barcode drawn from a set with a minimum Hamming distance of 4 and a constant Hamming weight of 4. Out of the 140 possible unique 16-bit barcodes satisfying our Hamming parameters, we used the 39 barcodes that were not assigned to any target mRNA as false-positive (‘blank’) controls. Each bit in the barcodes was then associated with a unique readout sequence (71, 72). As our barcodes use a Hamming weight of 4, each mRNA target is, thus, defined by four unique readout sequences. To create the sequences of the final encoding probes for a given mRNA, we concatenated two of the four readout sequences associated with that RNA to each of its target regions. To create template molecules for the construction of these encoding probes, a T7 promoter and two unique PCR primers were concatenated. Table S1 lists all targeted genes, operons, and their associated barcodes while table S2 contains the sequences of the probe template molecules and the associated fluorescently labeled probes complementary to these readout sequences, which we term readout probes. After the template probes were created for this library, we noticed two instances in which the genes we targeted are likely transcribed on the same operon (infB and pnp; and bamA and dnaE). As different barcodes assigned to the same transcript will likely be frequently corrupted, we chose to discard these four barcodes from all subsequent analysis, reducing this library effectively to a 97-operon library, which is how we refer to this library throughout. Nonetheless, these barcodes and their associated probes are maintained in tables S1 and S2 and the molecules decoded from these barcodes are present in the deposited MERFISH data for completeness.
For the design of the 1,057- and 1,930-operon libraries, we developed a more direct method for considering polycistronic operon structure. Specifically, we leveraged the operon structure annotations provided in Ecocyc v27.1 (73). Operons can have alternative transcription start and stop sites, leading to differences in the genes contained on polycistronic messages. As such variation raises the possibility that two different barcodes could be assigned to genes on the same message, we adopted a conservative approach to defining polycistronic messages. Specifically, we merged overlapping operon annotations to generate consensus sequences that did not share genes with any other consensus operon. We named the consensus operons based on the set of all the genes encompassed by the merged operon variants. Next, we used the same pipeline as described above to design target regions to these consensus operons with two notable exceptions. First, we excluded a handful of the most abundant mRNAs (e.g., some ribosomal protein mRNAs), which would have disproportionately affected target mRNA density. Second, we allowed some overlap in the potential target regions as some operons were too short to support 50 non-overlapping target regions. For the 1,057-operon library we allowed probe homology regions to overlap by 10 nt while for the 1,930-operon library we allowed target regions to overlap by up to 20 nt. This degree of overlap might contribute to the detection efficiency drop observed for the 1,930-operon library relative to other libraries. In both cases, these libraries targeted all consensus operons which could accommodate at least 50 probes. However, we generated probes for more available target regions if an operon could support it. To encode the targets for the 1,057- and 1,930-operon libraries, we employed a 31-bit or a 40-bit barcoding scheme, respectively, with the same Hamming weight and distance constraints as described above. The encoding probes and their template molecules were also created as described above. The barcode sets and template probe sequences are listed in tables S1 and S2. The consensus operon sequences are available at the Gene Expression Omnibus (GEO) at accession GSE268480.
To design a 159-operon MERFISH library to profile the adaptation of B. theta to different niches in the mouse colon, we leveraged a published list (74) of SusC polysaccharide importers and hybrid-two-component systems (HTCS) associated with polysaccharide utilization loci (PUL) as well as published bulk RNA-sequencing data (75). In addition, we included a random selection of operons associated with core elements of central metabolism, including glycolysis, gluconeogenesis, and fermentation (tables S1 and S4). The same approach to target region design described above was used with a slight adjustment in the constraints on GC content (40% to 74%) and the overlap permitted for target regions (no overlap). We targeted individual genes within most operons but extended probe design to nearby genes in the same operon if the target gene could not accommodate enough probes. We designed at least 49 encoding probes for each gene with two readout sequences per probe. An 18-bit binary barcode set with the same Hamming properties as described above was used to define these 159 operons with 14 barcodes used as blank controls. Encoding probe templates were built for these target regions as described above. The targeted operons, utilized barcodes, and template sequences for the encoding probes are included in tables S1 and S2.
Encoding probe template libraries were purchased as complex Oligo Pools (Twist Biosciences) and amplified to create large quantities of single-stranded DNA encoding probes using protocols described previously (30, 76). Briefly, template oligo pools were amplified with limited-cycle qPCR to create in vitro template molecules, which were used to create large quantities of single-stranded RNA via a high yield in vitro transcription with T7 RNA polymerase. Encoding probes were then created via reverse transcription of these RNA templates. The RNA templates and the reverse transcription primer (which contained a terminal ribonucleoside) were then removed via alkaline hydrolysis. Probes were purified by solid-phase reversible immobilization (SPRI beads; assembled as described previously (76)), and encoding probes were concentrated to a stock concentration of 1 mM with ethanol precipitation. Probes were stored as single-use aliquots at −20 °C.
All smFISH probes used for validation experiments in E. coli and B. theta were designed using the approach described for the 1,057-operon library for E. coli or the 159-operon library for B. theta. We designed at least 55 target regions per targeted operon and concatenated a single readout sequence to all target regions associated with a specific operon. The sequences of these probes are provided in table S2. smFISH probes were synthesized as oPools by Integrated DNA Technologies (IDT) and used for staining directly without amplification.
We designed FISH probes to target and barcode the 16S ribosomal RNA (rRNA) of E. coli to provide dense staining of cells or, in the case of the diauxic shift measurements, to also distinguish cells harvested from multiple time points which were combined and measured in a single MERFISH measurement. These rRNA probes were based upon the common Eub388 probe sequence, which targets a region of the 16S rRNA largely conserved across all bacteria (77). We modified this sequence by extending the targeted sequence to 30 nt of homology to the E. coli 16S rRNA to allow hybridization in the same conditions as used for MERFISH encoding probes. This target sequence was concatenated to different readout sequences, and these probes were synthesized with a 5’-acrydite moiety to enable copolymerization in expansion gels (38). Additionally, these probes were synthesized with ribonucleotides to prevent their digestion during DNase treatments (RA83-RA48; table S2). We also designed a DNA version of the extended Eub388 probe without the acrydite for experiments that did not involve expansion microscopy (DA48; table S2). The rRNA probes were synthesized by IDT.
For samples cultured in lysogeny broth (LB; VWR, J106–500G), a single colony of E. coli K-12 MG1655 (CGSG #7740, Yale Coli Genetic Stock Center) grown on an LB-agar plate was inoculated in ~5 mL LB Lennox medium (Fisher, AC612725000) to create an overnight culture. The next day this overnight culture was diluted 100 in 30 mL of fresh LB, grown to an optical density (OD600) of 0.32–0.34, and 7–14 mL of culture was harvested. Unless otherwise specified, all bacterial liquid cultures were grown at 37 °C with shaking at 220 rpm.
For samples associated with the glucose-xylose diauxic shift, an overnight culture of E. coli was created by inoculating a single colony grown on an LB-agar plate into 5 mL MOPS defined minimal medium (DMM; Teknova, M2106) supplemented with 0.2% w/v glucose (Teknova, G0520). The following day this overnight culture was pelleted, the supernatant was removed, the pellet was resuspended in 2 mL of fresh DMM, and the OD600 was measured. This resuspended overnight culture was then diluted into 400 mL of DMM supplemented with 0.025% w/v glucose and 0.025% w/v xylose (Sigma, X3877) to a final OD600 of 0.004. This culture was grown, and small volumes (~2 mL) were drawn to monitor growth by their OD600. To harvest cells for MERFISH, a variable quantity of culture was collected for fixation with the volume adjusted (between 7 and 21 mL) such that samples harvested at different stages of the growth curve had comparable numbers of cells for the filtration and rRNA probe hybridization steps below.
For paraformaldehyde-fixed (PFA) samples, fresh 32% w/v PFA (Electron Microscopy Sciences, 15714) was added directly into the LB or DMM cultures to a final concentration of 4% w/v PFA, and the sample was fixed for 30 minutes at room temperature. The sample was then transferred into a 30 mL syringe and drawn through a 0.22 μm mixed-cellulose-ester syringe filter (Fisher, 09–720-004) to separate the cells from the fixation solution. The syringe was replaced with a new syringe, and reverse flow of 10 mL of 1×Phosphate Buffered Saline (PBS; Thermo, AM9625) across the filter was used to resuspend the cells. This solution was then pushed back through the filter to re-immobilize the cells, and this gentle wash was repeated for a total of three times. Upon the final elution from the filter, cells were resuspended in 3 mL nuclease-free water in a fresh syringe and transferred to a 50 mL tube containing 7 mL 100% ethanol (KOPTEC 200 Proof Ethanol; VWR 71001–866) for permeabilization. After 1 hour of incubation at room temperature, cells were pelleted, and the supernatant was removed. The pellet was then resuspended in 1 mL 1×PBS before proceeding to subsequent staining. For the diauxic shift experiments, samples collected at early time points were stored in 70% ethanol on ice until all other samples could be collected, fixed, and permeabilized as described above.
For methanol-fixed samples, the unfixed culture described above was pelleted, the supernatant was discarded, and the pellet resuspended in 100 μL water. The sample was fixed by adding 25 mL of 100% methanol (Sigma, MX0480) and incubating for 2 hours at room temperature. The sample was pelleted, the methanol removed, and the sample was resuspended in 1 mL of 1×PBS. This process was repeated for a total of three washes before the pellet was resuspended in a final 200 μL of 1×PBS before proceeding with subsequent staining.
Nuclease-free water for all protocols was generated by reverse osmosis (MilliporeSigma, Synergy UV) with a Biopack polisher (MilliporeSigma, CDUFBI001). Unless otherwise specified here and below, all bacterial pelleting steps were performed at room temperature via a 600×g spin for 10 minutes on a tabletop centrifuge.
40 mm-diameter coverslips (Bioptechs, 40–1313-03193) were used for holding samples throughout the expansion protocol or during imaging. When necessary, these coverslips were silanized as described previously (78). Briefly, coverslips were first cleaned in 1 37% HCl (Sigma, 258148) and methanol and then coated with an allyl-silane layer. The allyl moiety was selected as it incorporates into polyacrylamide gels, creating a covalent bond between the stabilization gel and the coverslip. Silanized coverslips were stored at room temperature in a desiccated environment prior to use.
In some cases, silanized or non-silanized coverslips were treated with poly-L-lysine (PLL) to enhance the adherence of samples to the coverslips. PLL-coated coverslips were prepared by washing coverslips in 70% ethanol, covering them with a solution of 0.1 mg/mL PLL (Santa Cruz Biotechnology, sc-286689), incubating at room temperature for 10 minutes, and then washing away the excess PLL twice with 1×PBS.
All mice were used in accordance with animal care guidelines from the Harvard Medical School Standing Committee on Animals and the National Institutes of Health under a protocol (IS00003215) approved by the Harvard Institutional Animal Care and Use Committee (IACUC).
We obtained mice monocolonized with B. theta from the Massachusetts Host-Microbiome Center. Germ-free C57BL/6 female mice were maintained in gnotobiotic isolators under a strict 12-hour light cycle and at a constant temperature (21 ± 1 °C) and humidity (55% - 65%). 10-week-old female mice were colonized with B. theta VPI 5482 by oral gavage (4.8 × 10^8^ colony-forming units [CFU]/mL) and maintained on a standard chow (LabDiet, 5021). Mice were euthanized with CO2 asphyxiation followed by cervical dislocation.
The entire colon was rapidly dissected and fixed either with methacarn or periodate-lysine-paraformaldehyde (PLP). We pursued two different methods for fixation to confirm that the observed spatial patterns in B. theta gene expression were not fixation dependent, selecting these particular methods as they can better preserve the delicate mucus layer.
For methacarn fixation, the colon was immersed in 60% methanol, 30% chloroform (Sigma, C2432), and 10% acetic acid (Sigma, 695092) and incubated at 4 °C for 48 hours. The samples were then washed with 100% methanol for 35 minutes at 4 °C for a total of two times, then washed with 100% ethanol for 30 minutes at 4 °C for another two washes.
For PLP fixation, the colon was rapidly immersed in 2% w/v paraformaldehyde, 2 mg/mL sodium (meta)periodate (Sigma, S1878), and 0.075 M L-lysine (Sigma, L8662) in 1×PBS and incubated for 3 hours at 4 °C. These samples were then washed twice with 1×PBS, and then dehydrated with successive ethanol washes increasing in concentration from 50%, 70%, 95%, to 100%. Each wash was performed at 4 °C for 30 minutes, and the final 100% ethanol wash was performed twice. Samples prepared with either fixation method were stored at −20 °C pending further processing.
Samples were then embedded in paraffin and sectioned by the Rodent Histopathology Core at the Dana-Farber/Harvard Cancer Center. Briefly, the samples were incubated in xylene (VWR, 89370–088) for 1.5 hours at room temperature twice and infiltrated with paraffin (Leica, 3801340) for 4–5 hours at 60 °C, followed by embedding in a paraffin (Leica, 3801320) block. Paraffin blocks were then cross-sectioned to 5 μm slices, which were mounted on PLL-coated coverslips. The six measured samples were taken from multiple locations in multiple fecal pellets from two mice, each fixed with either methacarn or PLP. Coverslips were dried at room temperature and stored at 4 °C prior to expansion or smFISH staining.
To hybridize the rRNA probes, permeabilized cells were washed twice in 1 mL RNA FISH wash buffer (30% v/v formamide [Thermo, AM9342] in 2×SSC [Thermo, AM9765]), and then resuspended in 50 μL RNA FISH hybridization buffer (30% v/v formamide, 10% w/v dextran sulfate [Millipore, S4030], and 1 mg/mL yeast tRNA [Thermo, 15401029] in 2×SSC) supplemented with the appropriate rRNA probe at a final concentration of 8 μM. Samples were incubated at 37 °C overnight to hybridize these probes. The next day, samples were pelleted, the supernatant removed, and the pellet was resuspended in 400 μL RNA FISH wash buffer. Samples were incubated 30 minutes at 47 °C to promote the melting of off-target probes, pelleted, and the supernatant was removed. This wash step was repeated a second time. For multiplexed samples, one additional wash with RNA FISH buffer was performed, and a similar number of cells from each time point, as judged by pellet size, were mixed. To remove the formamide carried over from the RNA FISH wash buffer, samples were pelleted, the supernatant removed, and the pellet was resuspended in 1×PBS. This wash was repeated a total of three times. Pelleting of samples in RNA FISH hybridization buffer was performed with a 600×g spin for 15 minutes given the increased viscosity of this buffer.
As the bacterial cell wall can resist expansion, the cell was digested prior to embedding in the expansion gel (41). Cells were pelleted and resuspended in cell wall digestion buffer (400 U/mL mutanolysin [Sigma, M9901] and 0.1 U/μL RNasin [Promega, N2615] in 1×PBS). The volume of the digestion buffer was adjusted to bring the concentration of the cells to an OD600 of ~1. 40 μL of cells were placed on a PLL-coated coverslip, which was placed in a Petri dish, and the solution was spread over a 1–2 cm^2^ area with the side of a P1000 pipette tip while avoiding direct contact with the coverslip surface. The coverslip was placed in a humidified oven for 2 hours at 37 °C to perform the digestion, and then washed delicately with 1×PBS for a total of three times. Critically, aspirating too strongly or dropping liquid directly onto the fragile cells at this step can severely damage the samples.
Our expansion gel recipe follows that introduced with the TREx method (39) with notable modifications in the handling of these gels for MERFISH. Briefly, single-use, 980-μL aliquots of expansion gel solution (1.1 M sodium acrylate, 2 M acrylamide [Bio-Rad, 1610140], and 60 ppm bis-acrylamide [Bio-Rad, 161–0142] in 1×PBS) were prepared and frozen at −20 °C. We found that the quality of commercially available sodium acrylate varied from lot to lot, and, as reported by other groups (79), we avoided sodium acrylate lots with a strong yellow or orange color or a cloudy appearance when dissolved in water as they produced inconsistent expansion. Due to these quality issues, we sourced sodium acrylate from multiple vendors (Sigma, 408220; or Santa Cruz Biotechnology, sc-236893). Like other groups (39), we also made sodium acrylate in house by neutralizing acrylic acid (Sigma, 147230) with sodium hydroxide (Sigma, 72068) to a pH range of 8–8.5 and adding water to adjust the final concentration of sodium acrylate to 38% w/v. We observed no obvious links between the quality of bacterial-MERFISH and the source of sodium acrylate; however, variations in the sodium acrylate source may have contributed to some variation in the degree of expansion (fig. S1, A and C).
To embed samples in expansion gels, a gel monomer solution aliquot was thawed and mixed with N’-tetramethylethylenediamine (TEMED; Sigma, T7024) and ammonium persulfate (APS; Sigma, 215589) to a final concentration of 0.10% v/v and 0.10% w/v, respectively, to initiate polymerization. Glass slides (Ted Pella, 260439) were coated with GelSlick (Lonza, 50640) according to the manufacturer’s instructions. 200 μL of this initiated gel solution was then spotted onto the coated glass slides. Excess 1×PBS solution on the digested sample coverslips was removed by gently tapping the side of the coverslips on a Kimwipe (Kimtech, Kimwipes), and each coverslip was inverted onto an initiated gel solution droplet to embed cells within a thin film of gel between the coverslip and the glass slide. The samples were placed in a nitrogen chamber (Embrient, MIC-101), a damp Kimwipe was added to maintain humidity within the chamber, and the chamber was then sealed and purged with nitrogen. The chamber was then placed at 37 °C for 2 hours to allow the gels to polymerize in an oxygen-free environment. Samples were removed from the chamber, and the coverslips were gently detached from the glass slides. The gels were then cut with a rigid razor blade (WB Mason, ATSP591915) to an asymmetric shape to distinguish the orientation of the gels after expansion.
To further digest cellular contents that would restrict expansion, the gels—still attached to the coverslip—were transferred into 6 cm diameter Petri dishes (VWR, 25384–092) and covered with 4 mL protease digestion buffer (8 U/mL proteinase K [New England Biolabs, P8107S] in 50 mM Tris [Thermo, AM9856], 1 mM ethylenediaminetetraacetic acid [EDTA; Thermo, AM9849], 0.5% v/v Triton X-100 [Sigma, T8787], and 0.8 M guanidine hydrochloride [Sigma, G7294]). The digestion was performed in a humidified oven at 37 °C overnight. The next day, the digestion buffer was removed, and the gels were quickly rinsed twice with nuclease-free water. The gels either detached from the coverslips during the digestion or, if still partially attached after digestion, were gently dislodged with a wash in nuclease-free water. The gels were then transferred into fresh 6 cm diameter Petri dishes for expansion. To further expand the samples, residual buffer was removed, gels were covered in ~10 mL nuclease-free water, incubated for 20–60 minutes, the water was removed, and this wash process was repeated (typically 3–5 times) until the gel visibly stopped increasing in size.
Because expansion gels change size in different buffers, which would be incompatible with the buffer exchanges required for MERFISH, we embedded samples within an acrylamide stabilization gel that prevented this buffer-dependent size change (38). Water was removed from the Petri dishes, and gels were covered in ~5 mL of stabilization gel solution (4% v/v 1 acrylamide/bis-acrylamide [Bio-Rad, 1610144] in water) for 30 minutes at 4 °C. The gel solution was then replaced with ~5 mL of initiated stabilization gel solution (4% v/v 1 acrylamide/bis-acrylamide, 0.05% w/v APS, and 0.05% v/v TEMED in water) spiked with a 10,000 dilution of 0.1 μm-diameter carboxylate-modified orange-fluorescent beads (Thermo, F8800), which were added to serve as fiducials during MERFISH imaging. The samples were incubated for another 30 minutes at 4 °C to allow this activated gel solution to fully penetrate. Next, excess stabilization gel solution was removed, and the gels—infused with initiated stabilization gel solution—were placed onto a nylon mesh sheet (McMaster-Carr, 9318T46), which was itself placed on top of an acrylic cutting board covered with parafilm (VWR, 13–374-12). The nylon and parafilm provided a convenient support for the manipulation of fragile expanded gels. Each gel was then cut into ~4 mm × 4 mm squares with a rigid razor blade, and these squares were transferred, sample-side-up, onto a second parafilm-coated acrylic board. A silanized coverslip was then placed on top of each sample, and the sample and coverslip were placed in a humidified, nitrogen-purged chamber (as described above) and incubated at 37 °C for two hours to allow the gel to polymerize and crosslink to the coverslip. We noted that gel stabilization caused a 30–40% linear shrinkage of the expansion gels relative to the pre-stabilization expanded size likely due to the addition of charged polymerization initiators. For this reason, our 50X expansion protocol does not reach the same final degree of expansion reported for the TREx protocol (39). The stabilized samples were either immediately further processed with the following steps of the 50X expansion protocol, used as input for the following steps of the 1000X protocol (described below), or stored in 2×SSC at 4 °C for later use.
MERFISH requires the penetration of encoding and readout probes into the gel, and we found that thick gels substantially increase both the time required for this penetration and background fluorescence. Thus, we trimmed gels to a uniform, thin thickness to enhance the rate of probe penetration. Briefly, the silanized coverslips carrying the samples prepared above were placed on a flat RNase-free surface with the stabilized gel facing up and briefly washed with nuclease-free water for a total of three times. A 2 mil (~50 μm) thick slotted shim (McMaster-Carr, 9722K25) was placed around the gel to serve as a height spacer, and the gel was gently thinned with a thin flexible razor blade (Electron Microscopy Sciences, 72003–01) using the shim to set the height of the razor blade and the thickness of the gel.
During the development of bacterial-MERFISH, we found that very lowly expressed genes were measured more frequently with MERFISH than expected from bulk RNA-sequencing, producing an apparent false-positive rate higher than the false-positive rate we estimated from our blank control barcodes. In our experience, the blank controls are an excellent predictor of false-positive rates for MERFISH in eukaryotic systems. Thus, we reasoned that there must be an additional source of false positives in these bacterial samples. RNA FISH probe binding to DNA is not typically considered to be a possible source of false positives in unexpanded samples as the DNA duplex prevents the binding of probes; however, we surmised that during expansion some portions of the bacterial genome might be effectively melted due to topological constraints, allowing the binding of MERFISH encoding probes to the genome. This effect may not have been noticed in previous eukaryotic expansion-MERFISH (35, 38) due, perhaps, to the modest degrees of expansion in those works, a differential propensity for topological-stretch-induced genomic melting in bacteria, or the greater fraction of low-abundance transcripts in bacteria relative to mammalian cells. To eliminate this source of false positives, we digested the genome in our samples. Briefly, excess water from the gel trimming step above was removed, and a hydrophobic barrier was drawn around gels with an Aqua-Hold PAP pen (Electron Microscopy Sciences, 71311). Gels were then briefly rinsed with water, the excess water was removed, and the gels were covered with 50 μL of DNA digestion buffer (2.4 U/μL murine RNase inhibitor [New England Biolabs, M0314L] and 0.1 U/μL DNase I in 1×DNase reaction buffer [Thermo, EN0521]). Samples were digested for 40 minutes at room temperature.
To hybridize MERFISH encoding probes, the DNA digestion buffer was aspirated, and the samples were washed with a high-salt hybridization-and-wash (HSHW) buffer (30% v/v formamide in 10×SSC) for 10 minutes at room temperature for a total of three times. To hybridize the probes, excess wash buffer was removed, and each gel was covered with 50 μL of HSHW buffer supplemented with 10 μM (97-operon library) or 100 μM (1,057- or 1,930-operon libraries) of the MERFISH encoding probe libraries constructed above. The samples were then hybridized in a humidified oven at 37 °C for 48–72 hours. After hybridization, the samples were washed in HSHW buffer at 47 °C for 30 minutes for a total of two times. Samples were then washed in 2×SSC for 2 minutes for a total of three washes before proceeding with imaging. We found that the increased salt concentration in these buffers, relative to the standard hybridization buffers used for MERFISH and smFISH (71, 78), improved binding of the encoding probes. This improvement may be caused by the negatively charged expansion gel inhibiting the penetration of the negatively charged probes, and high salt may shield these charges.
To perform a second round of expansion to produce 1000X-expanded samples, 50X-expanded samples were taken after the formation of the stabilization gel and before thinning. These samples were first thinned as described above; however, as we found it easier to handle thicker gels during the second round of expansion, a thicker, 3 mil (~75 μm) shim was used. The thinned gels were then incubated in 300–400 μL of the same expansion-gel monomer solution described above initiated with 0.10% v/v TEMED and 0.10% w/v APS. The gels were gently shaken on an orbital shaker for 30 minutes at 4 °C or room temperature to allow the full penetration of the gel solution. Excess gel solution was removed, and this step was repeated a second time with fresh initiated gel solution. The samples on coverslips were then inverted onto a fresh 200 μL droplet of initiated gel solution spotted on a Gel-Slick-coated glass plate. These samples were transferred to a humidified chamber, which was then purged with nitrogen, and the gels were allowed to polymerize for 2 hours at 37 °C. The coverslips with the gel attached were then gently removed from the glass plate, and any newly formed gel outside of the original 50X expansion gel was removed with a razor blade. Gels were rinsed with nuclease-free water, expanded (without additional cell wall or protease digestion), and embedded in a second stabilization gel as described for the 50X expansion protocol. As fiducial beads were already incorporated into the first stabilization gel, such beads were not added to the second stabilization gel solution. Samples were either stored in 2×SSC at 4 °C before proceeding or were immediately thinned to 50 μm, DNase-treated, and stained with MERFISH encoding probes using the same protocol described above for the 50X samples.
Unlike the original iterative expansion microscopy protocols (80), this approach does not degrade the previous gels but rather expands them together with the second expansion gel. This form of iterative expansion was inspired by the expansion revealing protocol (40), which showed that degradation of previous expansion and stabilization gels was not necessary for a second round of expansion.
E. coli in LB or glucose-xylose DMM were fixed with PFA, permeabilized, and stained with FISH probes using the same protocols described above for ribosomal probe staining but with a final concentration of 2 μM of the smFISH probes and 1 μM of the rRNA probe (DA48; table S2). Cells were washed, pelleted, and resuspended in 1×PBS. Cells in 1×PBS were deposited on a PLL-coated coverslip prepared as described above, incubated for 1 hour at room temperature, and then washed with 1×PBS.
Colon slices prepared above were deparaffinized by heating the samples at 60 °C for 20 minutes followed by four washes of xylene (Sigma, 534056), each performed at room temperature for 2.5 minutes. The slices were then washed in 100% ethanol for 3 minutes at room temperature for a total of two times, then washed in 95% ethanol for 1 minute and then 70% ethanol for 1 minute. Sections were post-fixed on PLL-coated coverslips by incubating in 4% v/v PFA in 1×PBS at room temperature for 10 minutes, and the PFA was removed with two 1×PBS washes.
We made two notable changes to the 50X expansion protocol for colon slices. First, in early pilot work, we found that a PFA-treatment alone did not seem to efficiently anchor RNAs into the gel for paraffin-embedded samples; thus, we adopted a modified form of an RNA anchoring method developed previously (81, 82). Second, we found it difficult to orient expanded samples without the clear landmark provided by the nuclei of the host gut; thus, we did not perform DNase treatment in these samples.
To create an RNA-reactive gel-crosslinking agent, we were heavily inspired by the MelphaX compound which combines melphalan (a bi-functional nucleic acid alkylating agent) with acryloyl-X (a gel-reactive acrylamide) to create a crosslinking agent that decorates nucleic acids with gel-reactive acrylamide groups (82). As Acryloyl-X is expensive, we reasoned that we could replace it with a much less expensive NHS-ester containing a methacrylic group (MA-NHS; Sigma, 730300), creating a MelphaX-like crosslinker we term MelphaMA. To synthesize MelphaMA, an 8 mM melphalan (Cayman Chemicals, 148–82-3) stock in DMSO (Invitrogen, D12345) was mixed with a 100 mM stock of MA-NHS in DMSO in a 1 ratio, and the reaction was incubated at room temperature overnight with gentle shaking. As NHS reactions are efficient in DMSO and the MA-NHS was in excess of the melphalan, we assumed that the reaction proceeded to completion and treated the final solution as 6 mM MelphaMA. Single-use aliquots of MelphaMA were stored at −20 °C in a desiccated container.
To treat samples with MelphaMA, the slices were washed once in 20 mM MOPS pH 7.7 (Sigma, M9381) for 30 minutes at room temperature and then incubated with 50 μL of a 1 mix of this buffer with 6 mM of MelphaMA in a humidified oven at 37 °C overnight. The samples were then washed twice in 1×PBS for 5 minutes at room temperature. To digest the cell wall, the samples were covered with ~50 μL cell wall digestion solution (800 U/mL mutanolysin and 0.1 U/μL RNasin in 1×PBS buffer) and incubated for 2 hours at 37 °C. After digestion, samples were washed with 1×PBS for a total of three times.
To expand the samples, sections were washed with 200 μL tissue expansion gel monomer solution which contained the gel monomer solution described above supplemented with 0.01% w/v 4-Hydroxy-TEMPO (4-HT; Sigma, 176141), 0.20 % v/v TEMED, and 0.20% v/v APS for 20 minutes at 4 °C for a total of two times. 4-HT was added to this solution to further slow the polymerization and allow better penetration of the solution into the tissue slice (37). The gel was then polymerized, the sample digested with protease, the gel expanded and stabilized, and MERFISH and rRNA probes were hybridized using the protocol described above for liquid cultures.
Colon samples were fixed, paraffinized, sliced, deparaffinized, and rehydrated as described above. smFISH staining was performed as described for E. coli samples with hybridization solutions that contained a total concentration of 1 μM per smFISH probe set and 1 μM of the DA48 rRNA probe (table S2). As this probe was based on the eubacterial Eub388 probe, we found that it also labeled B. theta 16S rRNA despite a few nucleotides of mismatch with the B. theta 16S rRNA sequence. Unlike the cultured E. coli samples, we found higher degrees of background in smFISH imaging of colon slices; thus, after the hybridization and wash of smFISH probes, we used a previously described embedding and clearing approach to reduce background (78). Briefly, the samples were washed in a non-expanding gel solution (4% v/v 1 acrylamide/bis-acrylamide, 50 mM Tris-HCl, 300 mM NaCl [Thermo, AM9759], 0.03% w/v APS and 0.15% v/v TEMED) for 2 minutes, then inverted onto a GelSlick-coated glass plate with a 65 μL droplet of the same gel solution. This thin film of gel was polymerized for 2 hours at room temperature, and then the coverslip with the gel was gently removed from the glass plate, washed with 2×SSC, and digested in a clearing buffer (2% v/v sodium dodecyl sulfate [SDS; Thermo, AM9822], 8 U/mL proteinase K, and 0.25% v/v Triton X-100 in 2×SSC) at 37 °C overnight. The next day samples were washed with 2×SSC for 30 minutes at room temperature for each wash with a total of five washes and then stored in 2×SSC at 4 °C prior to imaging.
All samples were imaged on an automated home-built microscope system as described previously (76). Briefly, this system is comprised of an epifluorescence microscope body (Nikon Ti-2) with a Celesta laser light engine (Lumencor, 90–10521), a 60× CFI PlanApo oil objective (Nikon), a 10× CFI PlanApo air objective (Nikon), and two CMOS cameras (Hamamatsu, ORCA Flash 4.0) with a twin-camera color splitter (Cairn, TwinCam). The sample was illuminated at 750 and 635 nm for MERFISH and smFISH signals; 750 nm, 635 nm, and 488 nm for rRNA probes; 555 nm for fiducial beads; and 405 nm for 4’,6-diamidino-2-phenylindole (DAPI). Samples were mounted in a closed-flow chamber (Bioptechs, FCS2) with a 1 mm-thick flow gasket (Bioptechs, 1907–1422-1000). Buffers were programmatically delivered via a home-built flow system comprised of a peristaltic pump (Gilson, MINIPULS 3) and a daisy-chained set of valves (Hamilton, MVP).
The first set of readout probes complementary to the readout sequences associated with the first two bits were stained prior to loading the sample on the microscope. The storage buffer was aspirated from samples, and they were covered with 5 mL of readout hybridization buffer (10% v/v ethylene carbonate [Thermo, A15735–36] and 0.125% v/v Triton X-100 in 2× SSC) containing 3 nM of each of the fluorescently labeled readout probes. Samples were incubated for 30 minutes at room temperature in the dark with gentle orbital shaking, and then washed with readout hybridization buffer that did not contain readout probes for 20 minutes at room temperature for a total of two washes. For samples that were not DNase treated, the wash solutions were supplemented with 8 μg/mL DAPI (Thermo, D1306). Samples were then washed once with 2×SSC before imaging.
For experiments that require multiple rounds of readout staining, imaging, and removal, a flow cartridge was assembled containing readout hybridization buffers containing 6 nM (for MERFISH samples) or 3 nM (for smFISH samples) of the appropriate fluorescently labeled readout probes for each of the hybridization and imaging rounds. Once loaded onto the microscope, the sample was covered with an imaging buffer (4 μM Trolox-quinone (83), 0.5 mg/ml Trolox [Abcam, AB120747], 500 recombinant protocatechuate 3,4-dioxygenase [rPCO; OYC Americas, 46852004], and 5 mM protocatechuic acid [Sigma, 37580] in 2×SSC) designed to reduce photobleaching and enhance fluorophore brightness, and the desired fields-of-view (FOV) were imaged. 10 to 15 z-planes, separated by 1 μm, were imaged per FOV with the number of planes set by the degree of expansion. Once imaging was complete, the buffer was exchanged with a readout cleavage buffer (50 mM Tris(2-carboxyethyl)phosphine [GoldBio, TCEP25] in 2×SSC) and incubated, with gentle flow, for 15 minutes. This buffer was then exchanged with 2×SSC to remove residual cleavage buffer, and the next readout hybridization solution was added to the sample. The sample was incubated in the readout hybridization buffer for a total of 30 minutes, with gentle flow during the incubation, and then this buffer was exchanged with a wash buffer identical in composition to the readout hybridization buffer but lacking readout probes, in which the sample was incubated for 11.5 minutes. This buffer was then exchanged with the imaging buffer, and the sample was incubated in this buffer for 10 minutes prior to the start of the next imaging round to allow the rPCO to scavenge any oxygen introduced into buffers during flow. For E. coli samples, we performed 8 (97-operon library), 16 (1,057-operon library), or 20 (1,930-operon library) rounds of staining and imaging to identify all readout sequences. Where necessary, additional rounds of staining and imaging were performed to identify all utilized rRNA probes. For B. theta samples, we performed 9 rounds of staining and imaging for the 159-operon library. Readout probes, conjugated to Alexa488, Cy5, or Alexa750 via disulfide bonds, were synthetized by Bio-Synthesis, and the sequences are provided in table S2.
To identify molecules from the raw MERFISH images collected above, we utilized a published analysis pipeline (https://github.com/ZhuangLab/MERFISH_analysis) (71, 72). Briefly, this pipeline registered and aligned images taken of the same FOV from different imaging rounds with affine corrections built from the location of the fiducial beads to account for imperfect stage movements and affine corrections determined by multi-color samples to correct for small chromatic aberration in the imaging optics. Next, background noise was removed with a high-pass filter, and RNA signals tightened with Lucy-Richardson deconvolution. The intensity profile was then normalized across all imaging rounds, and the normalized intensity profile across all imaging rounds for each pixel was then compared via Euclidean distance to the profiles expected from all barcodes, and a pixel was assigned to a barcode if it was within a distance threshold set by a single bit flip. Adjacent pixels assigned to the same barcode were merged to create an identified RNA. Background signal was rejected by removing putative molecules with signal spread over too few pixels (the area) or with an average intensity too low (the brightness). In addition, molecules identified in the bottom-most z-plane were excluded from some datasets due to an increased false detection rate caused by nonspecific probe binding to the coverslip surface. For E. coli measurements, where we wanted to provide a conservative estimate of our detection efficiency, we identified likely instances in which the same RNA molecule was detected in different z-planes using DBSCAN to identify RNAs of the same type within 2 μm in the axial dimension and 0.2 μm in the lateral dimension and kept only the brightest molecule in each group. The RNAs identified with this pipeline are provided for all measured datasets at https://doi.org/10.5061/dryad.n5tb2rc4d.
To provide a validation of the mRNA abundance determined by bacterial-MERFISH, we collected bulk RNA-sequencing data. E. coli were cultured in LB or glucose-xylose DMM as described above. Samples were harvested at an OD600 of 0.33 (LB), 0.18 (glucose growth phase), or 0.47 (xylose growth phase). Total RNA was extracted using the RNAsnap protocol (84). Briefly, 4.5 mL of bacterial culture was mixed with 0.5 mL ice-cold stop solution (10% v/v phenol, pH 6.6 [Fisher, BP1750I-100] in ethanol), cells were pelleted with a 600×g spin for 8 minutes, and the supernatant was discarded. The pellet was then resuspended in 0.5 mL RNAsnap extraction buffer (9.5 mL 100% formamide, 360 μL 0.5 M EDTA, 25 μL 10% SDS, and 100 μL 100% 2-mercaptoethanol [Fisher, O3446I-100]), and the sample was incubated at 95 °C for 1 minute. Samples were stored at −20 °C as needed before proceeding. RNA was purified using the Zymo Direct-zol RNA kit (Zymo, R2051) by mixing 100 μL of the total RNA in the RNAsnap extraction buffer with 200 μL of the Zymo RNA binding buffer and then following the manufacturer’s instructions.
Library preparation and sequencing for the DMM samples were performed by the Harvard Medical School Biopolymers Core using the Ribo-Erase kit for rRNA depletion, the KAPA HyperPrep kit for library preparation, and 150-bp paired-end sequencing on an Illumina NextSeq 500. Library preparation and sequencing for the LB sample were performed by Azenta using TURBO DNase for DNA depletion, the QIAGEN FastSelect rRNA HMR Kit for rRNA depletion, and the NEBNext Ultra II RNA Library Preparation Kit for Illumina according to the manufacturer’s instructions. The library was 150-bp paired-end sequenced on an Illumina NovaSeq 6000.
Sequencing data were mapped to transcript abundance using salmon v1.8.0 (85). First, to facilitate comparison to the consensus operon structures used for the design of the MERFISH target regions, a salmon index was built on the sequences of these consensus operons. Salmon was then used, via the quant function, to measure the mapped count, effective length, and transcripts per million reads (TPM) for each operon. For consensus operons with disjoint probe design regions (e.g., when a noncontiguous subset of genes on the operon were targeted), we recomputed the TPM of the operon using the average mapped counts of contiguous regions weighted by their effective lengths. Both sequencing data and the sequences of the consensus operons were deposited in GEO under accession GSE268480.
To segment cells, an interactive segmentation toolkit, Ilastik 1.4.0 (86), was used to perform 3D pixel classification on the rRNA probe images. For computational efficiency, these images were downsized from 2048×2048 pixels to 512×512 pixels for 50X-expanded samples or to 256×256 pixels for 1000X-expanded samples. We used a human-in-the-loop Ilastik workflow to create a pixel classifier. For the LB measurements, this classifier classified pixels as in or out of a cell. For the diauxic shift measurements, pixels were classified as in a cell for a specific rRNA channel or out of a cell. The pixel classifiers for all data types were built using color, edge, and texture namely, we used the GaussianSmoothing, LaplacianOfGaussian, GaussianGradientMagnitude, DifferenceOfGaussians, StructureTensorEigenvalues, and HessianOfGaussianEigenvalues feature options provided by Ilastik with a series of 2D (1–40 pixels) and 3D sigma values (0.3–3.5 pixels). The range of the sigma values was tuned on a dataset-by-dataset basis based on the density of the rRNA stain, the image downsampling, and the degree of expansion for individual samples within the ranges provided. Ilastik was then used to predict class probabilities for the pixels of these downsampled images, which were then used to generate segmentation masks.
For 50X-expanded samples, preliminary masks were produced by thresholding on the pixel probabilities of the in-cell classification (>0.5), followed by morphological opening (with a radius of 1–2 pixels depending on the dataset) to separate contacting cells. These preliminary masks often had small internal regions of pixels that were below the in-class probability threshold. To fill these holes, a morphological closure operation (with a radius of 1–2 pixels depending on the dataset) was applied. This closure operation was applied to the mask of each cell separately to prevent the merger of adjacent cells. Contiguous regions within these masks were then identified and labeled as individual cells using scikit-image (87).
For 1000X-expanded samples the density of rRNA staining was not always sufficient to provide high-probability in-cell classification throughout the entire cell volume; thus, we adopted a seeded-watershed approach to define cell masks. To define the ‘seeds’, i.e. the locations in the image that most likely contain a cell, we used a dataset-specific, stringent in-cell probability threshold (0.4–0.7). To then define a foreground mask, we used a dataset-specific, less-stringent in-cell probability threshold (0.35–0.4). We then inverted the in-cell probabilities to create an image where high values correspond to regions of low cell probability. We next performed seeded-watershed on these inverted in-cell probability images in combination with the seeds and foreground masks using the watershed routine provided by scikit-image to create masks for individual cells.
Masks created for the 50X-expanded or 1000X-expanded samples were then filtered by volume (>60 pixels for 50X-expanded samples; >300–500 pixels for 1000X-expanded samples) and the integrated rRNA signal intensity (>1000 arbitrary units [AU] for 50X-expanded samples; >100 AU for 1000X-expanded samples) to remove spurious masks. This approach often produced mask boundaries that had a few-pixel roughness to the edge, as opposed to the smooth boundaries expected. Thus, a binary closure operation (1 pixel-radius for 50X-expanded samples or 5-pixel-radius for 1000X-expanded samples) was performed on the binary masks of each cell. To further smooth the edges of masks, the binarized mask of each cell was individually blurred with a Gaussian filter (sigma of 2), and a new binarized mask was generated by thresholding the values of this filtered image with a threshold of 0.2–0.3 for 50X-expanded samples, depending on the dataset, or 0.33–0.4 for 1000X-expanded samples.
For the diauxic-shift measurements with multiple rRNA channels, the Ilastik in-cell probabilities for each rRNA channel were processed separately as described above. This produced a set of segmentation masks for each channel, effectively associating each cell with a specific rRNA channel. However, in some instances the same pixels were assigned to multiple rRNA channel masks. To reconcile these conflicts and create consensus cell masks, these disagreements were settled by randomly selecting a priority for the rRNA channels and discarding overlapping masks from all but one channel based on this priority. To avoid preferentially favoring one channel over the others, this rRNA channel priority was randomized for each FOV.
RNAs were assigned to individual cells if they fell within the boundaries of the 3D masks as derived above. A small fraction of RNAs were found close to specific masks but just outside the boundary of these masks. Thus, we also assigned any mRNA within 1.2 μm (after expansion) of a mask to its nearest cell. When measuring the internal organization of the transcriptome, we noticed that outlier RNAs could have a disproportionate effect on the apparent coordinate system within each cell; thus, we adopted a more conservative distance threshold of 1.0 μm for these measurements.
To determine the degree of expansion in expanded E. coli cells, we started with 3D segmented masks generated as described above. For each mask, we measured the cell solidity and used it to eliminate segmentation artifacts or cells that were highly curved (which would challenge width measurements). Solidity was defined as the ratio of the volume of the segmentation mask to the volume of the convex hull that contains all the pixels within the mask. Segmentation artifacts or curved cells will have solidity values much less than 1, so, for this analysis, all cells with a solidity value less than 0.8 were removed. To define the width of cells, the mRNAs associated with each cell were projected onto a single 2D plane and the convex hull containing them was constructed. Given that the cells were oriented in various angles in space, cells were oriented by calculating a minimum bounding box around each convex hull. The dimensions of this box were then used to define a long axis which corresponded to the maximum length of the cell and a short axis which corresponded to the maximum width. As the poles can be narrower, we removed both poles by dividing the hull into four evenly spaced parts along the long axis and discarding the first and last quarters. The width of the cell was then estimated from the length of the intersection of lines perpendicular to the long axis with the boundaries of the convex hull. To estimate the width of each cell, multiple, perpendicular lines were measured at random positions along the central line, and the median of these values was used for the final width. To obtain an estimate of the linear expansion factor of each cell, the measured width was compared to a previous report of the width of E. coli MG1655 cell grown in similar conditions (1 μm) (88) (fig. S1A). To estimate the volumetric expansion factor, this linear expansion factor was cubed (fig. S1C). The expansion factor of each replicate was computed as the mean of the expansion factor of each cell. To compute the normalized widths, we divided each cell’s width by its sample expansion factor (fig. S1B). Variation in the normalized widths was compared against previous measurements of the distribution of the width of E. coli grown in the same conditions (89).
To compare the average abundance per mRNA determined with bacterial-MERFISH to that determined via bulk RNA-sequencing, the Pearson correlation coefficient was computed between the log10 expression of each mRNA measured with both techniques. For datasets that were not segmented, MERFISH abundance was measured in average counts per FOV. For datasets that were segmented, MERFISH was measured in average counts per segmented cell.
As MERFISH is an image-based method, some cells may only be partially imaged. For example, cells that were aligned vertically or at an angle with respect to the coverslip would not be fully captured within the imaged z-stack. Similarly, cells that fell across the boundary of a FOV would also not be fully imaged. In parallel, errors in our segmentation masks can group two or more cells, effectively creating multiplet cells. For comparisons to bulk RNA-sequencing, all cells, including these partially imaged cells or potential multiplets, were retained in the analysis. However, to provide a measure of the detection efficiency, the number of mRNA copies, and the number of unique expressed operons per cell, we excluded potential multiplet cells and cells that were only partially imaged. For the 50X-expanded samples, cell multiplets were excluded by setting a minimum solidity threshold of 0.6, while partially imaged cells were identified as cells that were too small (less than 20 μm^3^ mask volume), too close to a FOV boundary (within 120 pixels), or whose fraction of imaged cell volume within the top-most imaged z-plane was too large (more than 10% of mask pixels). For 1000X-expanded samples, we used a minimum mask solidity of 0.5, a minimum volume of 2000 μm^3^, and a maximum of 7% pixels in the top z-plane. In both approaches, cells that were vertically aligned on the coverslip were identified and removed by requiring a minimum length to width ratio of 2, and occasional spurious masks were removed by eliminating all putative cells with less than 10 measured mRNAs. Finally, focus errors and large auto fluorescent debris can occasionally corrupt MERFISH measurements in some FOVs, leading to a dramatic reduction in the number of properly identified RNAs. Corrupted FOVs were identified from abnormally low numbers of detected mRNAs using the LocalOutlierFactor function of scikit-learn with the default parameters, and all cells within those FOVs were excluded. Cells that passed all these cuts were deemed ‘Whole’ cells. Finally, the blank controls were not included in the number of detected mRNAs or operons per cell.
To estimate the detection efficiency of bacterial-MERFISH, published and calibrated bulk RNA-sequencing data (i.e., where abundance was measured in units of counts per cell) for E. coli grown in the same conditions was used to calibrate our own bulk RNA-sequencing data. First, the measured mRNA abundances in Bartholomäus et al. (33) (measured in reads per kilobase per million mapped reads) was converted to mRNA copy numbers by cell by scaling these reads by the ratio of the reported total copy number of mRNAs per cell (7,800) to the sum of the RNA abundances for all mRNAs. A conversion factor between the measured RNA abundances in our own bulk RNA-sequencing data (measured in TPM) to copy numbers per cell was calculated as the ratio of the total copy number per cell of all monocistronic operons in the Bartholomäus et al. measurements to the total TPM measured in our own data for the same operons. Bartholomäus et al. report abundances for genes instead of operons; thus, to avoid the challenge of properly estimating polycistronic mRNA abundances from those of the multiple genes they contain, the calculation of the conversion factor was restricted to monocistronic operons. The measured TPM of all the operons in our sequencing data were then multiplied by this conversion factor to obtain counts per cell. To determine the detection efficiency of MERFISH, we fit a line of slope 1 to the copy numbers per cell in each MERFISH measurement to the calibrated bulk RNA-sequencing abundances for the operons included in the measurement, using log10 expression to weight lowly and highly expressed genes more comparably. The detection efficiency was determined from the intercept of this line fit.
As a further validation of bacterial-MERFISH measurements, we compared the average abundance of bacterial-MERFISH measurements to published measurements for E. coli grown in similar conditions from multiple scRNA-seq methods. To compare measurements made for E. coli grown in LB, we averaged the measured counts per cell across both replicates of the 1,057-operon bacterial-MERFISH in 1000X-expanded samples for all ‘Whole cells’. We then downloaded the published E. coli measurements for PETRI-seq, microSPLiT, and M3-seq, filtered cells following the filters described by the authors of those works, and then, where needed, converted counts per gene to counts per operon.
For PETRI-seq (8) we used the counts provided in table S6 of that publication for exponential E. coli in experiment 2.01, where E. coli MG1655 cells were grown in MOPS EZ rich defined medium at 37°C and collected at an OD600 of 0.4 (cell names with the prefix “SB442” with barcodes from plate columns 1–6). We then excluded cells with fewer than 128 total counts (including rRNA). PETRI-seq reports their measurements in terms of total counts per operon; however, because PETRI-seq and bacterial-MERFISH occasionally use different operon definitions, we restricted our analysis to the 1,006 operons that include the same genes in both datasets.
For microSPLiT (9), we used the count matrix files provided by the authors on GEO (GSM4594094) and extracted the E. coli cells grown in LB at 37 °C and collected at an OD600 of 0.5 (cells from wells >25 whose name contains “K12”). We then excluded cells with fewer than 200 total reads (including rRNA and tRNA) or 100 mRNA reads (using Ecocyc v28.1 “All-Genes Class: tRNA” [“ECOCYC-CLASS&object=BC-2.2.5”] to classify tRNAs). microSPLiT reports counts based on genes; thus, we leveraged the same operon structure we used for MERFISH probe design to aggregate gene counts into operon counts. Specifically, we added all counts associated with all genes found in each of our operons. We did note some differences in the gene name conventions; thus, out of an abundance of caution, we discarded all operons for which the set of gene names was not identical.
For M3-seq (14), we used the count matrix provided on GEO (GSM7306272) for rRNA-depleted E. coli grown in LB at 37°C and harvested at an OD600 of 0.6. The count matrix was pre-filtered by the authors of M3-seq to include all cells with more than 10 UMIs; thus, we performed no additional filtering. Since M3-seq, like microSPLiT, reports counts per gene, we used the same aggregation strategy to determine the counts per operon and adopted the same rejection criteria for the small number of genes with different names.
As seen in fig. S3, we noticed that these three techniques produce counts per operon that are larger for longer operons than for shorter operons relative to bacterial-MERFISH counts, with mRNA length determined from our operon definitions. As all three methods can generate a sequencing read from any region of an mRNA, this length-dependent capture efficiency is expected. To confirm that this effect was not dependent on bacterial-MERFISH counts, we compared all three techniques to the bulk RNA-sequencing data for E. coli in LB we collected as described above. Indeed, this enriched frequency of detection of longer mRNAs is still apparent. The standard adjustment for this effect in bulk RNA-sequencing methods where longer mRNAs also have a greater probability of producing a sequencing read than shorter RNAs is to divide by mRNA length. We applied this correction to the total counts per operon for each of these techniques and found that this length-dependent increased detection of longer mRNAs was removed. This length-dependent detection probability complicates a direct comparison between the counts per cell determined with these sequencing techniques and that determined via bacterial-MERFISH; nonetheless, our counts per cell remain comparable to those detected with these methods, consistent with the detection efficiency reported by these methods.
We also compared bacterial-MERFISH measurements collected during the diauxic shift to previously published scRNA-seq measurements collected in similar conditions for E. coli. As above, we averaged both replicates of the bacterial-MERFISH measurements of 1,057 operons in 50X-expanded samples harvested during the glucose log phase. We then compared these to ProBac-seq (10) measurements made in E. coli grown in a similar medium (M9 minimal media supplemented with glucose) at 37°C. Specifically, we downloaded the count matrices from GEO (GSM6994840) and did not apply any cell filtering, as none was described by the authors. Unlike the random priming methods above, ProBac-seq uses targeted probes to profile gene expression, and most genes are targeted with the same number of probes. Thus, a length correction is not required. We, therefore, compiled operon abundance by averaging the counts for all genes within a given operon. The gene counts were obtained by selecting the maximum probe count for each gene, as recommended by the ProBac-seq authors.
Three 1,057-operon MERFISH measurements were collected from two biological replicates of E. coli grown in glucose-xylose DMM, cells were fixed, 50X-expanded, rRNA-barcoded, measured with MERFISH, and segmented as described above. For Replicate 1, two bacterial-MERFISH technical replicates of the same biological replicate were combined. Counts per cell and cell metadata from all measurements were combined into a single anndata object (90) for downstream analysis.
To exclude segmentation artifacts and cells for which only a small fraction of the cell was measured, cells were filtered by volume (>25 μm^3^), solidity (>0.6), MERFISH read counts (≥15 per cell), and number of unique operons detected per cell (≥7). Long cells tend to have lower solidity due to their increased curvature; thus, to better retain such cells, the minimum solidity requirement was relaxed to 0.5 for cells with a volume greater than 300 μm^3^. The dimensions and orientation of each cell were identified from the principal components (PC) and eigenvalues of a principal component analysis (PCA) on the coordinates of the pixels within the masks. As cells are longer than they are wide, we found that the first PC aligned with the length of cells while the second and third PCs provided orthogonal measurement directions for the diameter of the cells. This analysis allowed us to further exclude cells with odd volumes given their lengths or with disagreements between two different estimates of the diameter, which are likely the result of segmentation artifacts. Specifically, cells were removed if their length (estimated from the first eigenvalue) to volume ratio was less than 0.35 or the ratio of the two diameter estimates (defined as the second eigenvalue over the third) was above 1.9.
The filtered anndata was then analyzed with the single-cell analysis toolkit scanpy (90) with default parameters unless otherwise specified. First, the number of operon counts per cell was normalized to a fixed target of 100 counts per cell, a pseudocount of 1 was added and these expression values were log normalized (with the natural base) and z-scored. Batch effects between the replicate MERFISH measurements were corrected with harmony (91) using the first 20 principal components as a basis. Differentially expressed genes were identified using the Wilcoxon rank-sum method on the log2-transformed normalized counts (with a pseudocount) without batch correction. A nearest neighbor graph with a neighborhood size of 20 was computed using the cosine distance metric on the first 20 harmonized PCs and was then used to compute a Uniform Manifold Approximation and Projection (UMAP) embedding (computed with a minimum distance of 0.2), a diffusion map, and to perform Leiden clustering with a resolution parameter of 1.3. Leiden clustering produced a total of 17 clusters, including three small clusters (<1% of all cells) which were excluded from subsequent analysis. These clusters did not display strong differential expression of any operons and were enriched in blank barcodes, suggesting that they were contaminated by false positive counts. The reported marker genes for the Leiden clusters were selected from the top 15 operons per cluster with a minimum log fold change of log2(2.5) and the highest marker significance scores reported by scanpy using the Wilcoxon rank-sum method. To facilitate comparison with a prior study (10), we additionally manually appended lapAB-pyrF-yciH to the retained markers.
To compute the pseudotime for cells transitioning from glucose log-phase growth through the diauxic shift, a cell in the late shift cluster N10 was selected as a starting point for the pseudotime trajectory, and the standard scanpy algorithm for diffusion pseudotime was applied (tl.dpt). We chose to root the pseudotime analysis in cluster N10 as opposed to glucose log-phase growth due to the transcriptional heterogeneity seen in log-phase growth. For visualization purposes, pseudotime calculated in this fashion was then reversed to match the actual direction of the experiment from glucose into the diauxic shift. The second half of the pseudotime range was then monotonically rescaled to 0.9–1 with a piecewise-linear function to more equally weight the pseudotime ranges distributed across each diauxic condition and the adjusted pseudotime was divided into 10 equal bins. As we were interested in the response to glucose starvation, we did not include samples collected during the xylose-growth or stationary phase of the growth curve in this analysis. The reported operons for the pseudotime analysis were expressed in at least 3% of cells harvested in the remaining conditions and which contained at least one gene carrying the GO annotation of “carbohydrate metabolic process” or “carbohydrate transport” in Ecocyc v27.1. When extending the pseudotime analysis to the entire library, we used a less stringent expression threshold of 1 count per 100 cells to increase the number of retained operons. We note that this new threshold remains higher than our estimated false-positive rate for the diauxic shift experiment (fig. S4, C and D). The mean log2-transformed normalized expression was computed for each pseudotime bin and then z-scored operon-wise. Finally, the z-scored pseudotime expression profiles were clustered with spectral clustering using scikit-learn and the following n_clusters=8, assign_labels=“discretize”, and random_state=0.
To assess the transcriptional heterogeneity within each condition, we applied the preprocessing pipeline described above to each sample individually using a Leiden resolution of 0.6 for clustering. As expected, this analysis produced a diversity of clusters, and we kept only clusters which had at least one marker with a significance score greater than 10 as determined with the Wilcoxon rank-sum method of scanpy. The marker genes shown in fig. S6 represent the top 10 operons with the largest significance scores subject to the requirements that they have an absolute log2 fold change greater than 0.5 and were in the list of markers derived from any clusters when all cells were analyzed (Fig. 2H).
The main analysis of the E. coli transcriptome organization was performed on cells taken from both replicates of the 50X-expanded, 1,057-operon MERFISH measurement of E. coli grown in LB; however, supporting analysis was also performed on 1000X-expanded samples. To further simplify this analysis, ‘Whole’ cells (as defined above) were further cut to retain only straight cells with sufficient RNA counts using a solidity threshold of 0.7 and a count per cell threshold of 40 for 50X-expanded samples, or 0.6 and 80 for 1000X-expanded samples.
To enable the comparison of intracellular RNA localization between different cells, a normalized cylindrical coordinate system was built for each cell, where each measured mRNA was localized by its position along the axial (a) or radial (r) axis. Formally, a and r were respectively defined as the distance from the RNA to the midcell along the long and short axes of the cell, normalized by the estimated length and diameter of the cell, respectively. We note that this definition ignores the distinction between ‘old’ and ‘new’ poles, which our measurements were not designed to identify, and assumes a radial symmetry to mRNA distributions.
To define this coordinate system for cells measured with 50X expansion, the orientation of each cell was determined using PCA on the 3D coordinates of the mRNAs associated with that cell, where PC1 indicated the long cell axis, and PCs 2 and 3 identified the radial plane. mRNAs were mapped to our intracellular coordinate system by projecting their 3D coordinates onto these three PCs and then normalizing by the dimensions of the cell. The cell dimensions were defined for each PC as the maximum projected distance between any two mRNAs in the cell. For each RNA, a was then given by the absolute value of the projection along PC1 and r by the L2 norm of the normalized projected distances on PCs 2 and 3.
Because of the lower final RNA density and relaxed solidity requirements in the 1000X-expanded samples, we found that the approach used for the 50X-expanded samples did not perform well. Thus, for the 1000X-expanded samples, we instead used the segmentation masks to define the intracellular coordinates. First, to correct for the differential resolution in X and Y versus Z, the segmentation masks were upsampled in Z to match the pixel size in X and Y. Each mask was then skeletonized with scikit-image to identify the central line of the cell. As skeleton endpoints are associated with the cell poles, cells with more than two end points in the skeleton were discarded as segmentation artifacts. This skeletonization process slightly eroded the cell central line from each of the two poles. To recover the original position of the poles, a second order polynomial was fit on the set of 10 pixels at each terminal end of the skeleton, and the intersection of these fit curves with the original cell mask was the used as an estimate of the pole location. The two polar intersection points as well as the pixel coordinates of the skeleton central line were then interpolated with a smooth spline, which was then resampled at equally spaced points to produce the final central line. To determine the values of a and r for each mRNA within the cell, r was defined as the Euclidean distance from that mRNA to the nearest point on the central line and a was defined as the integrated distance along the central line from this nearest point to the center of the cell. The length of the cell was measured as the full length of the center line, and the average radius of the cell was determined as the radius of a cylinder of cell length which would produce the volume of the cell segmentation mask.
To estimate the average distribution of individual operons, we used kernel density estimation (KDE). The volume of a cylindrical section depends non-linearly on the radius; thus, to create a volume-corrected 2D density projection, we first squared r to r^2^. To address boundary effects introduced by mapping mRNAs to the nonnegative quadrant of a symmetric cell, mRNAs were mirrored along the a and r axes. These mirrored mRNAs were then used to compute a KDE with scikit-learn and a bandwidth of 0.07. The KDEs were finally evaluated on 300 linearly spaced points along a and 100 linearly spaced points along r^2^, extending the evaluation range beyond the expected cell dimensions on each side by ~10% for visualization purposes. As the quality of the density estimates decreases when fewer observations are available and these spatial estimates could be influenced by false positive counts, mRNAs with fewer than 500 measured molecules across all cells in both replicates were excluded. For the 50X expanded samples, this threshold left 510 operons for the analysis presented in Fig. 3.
To numerically compare different spatial patterns, the 2D densities described above were flattened into 1D vectors, these vectors were reduced in dimension to 150 components with PCA, and the pairwise cosine distances between different spatial patterns were calculated. These distances were then used to create a nearest neighbor graph using 15 nearest neighbors, which was then used, in turn, to create a UMAP embedding (with a minimum distance of 0.5) and to define Leiden clusters with a resolution parameter of 0.3. To determine the reproducibility of these spatial patterns, 2D KDEs for each mRNA were computed for each of the two replicates separately using the same approach. For this analysis, the minimum distance used to create the UMAP embedding was set to 0.25 and the Leiden clusters were defined with a resolution of 0.5.
To explore the correlation between the predicted location of the encoded proteins and the spatial patterns of mRNAs, PSORTdb 4.0 (92) was used to annotate E. coli operons with the predicted cellular location of the encoded protein. As polycistronic operons contain multiple genes which may have different protein locations, we defined the annotations of polycistronic operons as the union of the predicted locations of the individual genes they contain with one notable exception. In a co-translational insertion model of mRNA localization, co-translational insertion of inner-membrane-protein-encoding mRNAs should enrich such mRNAs at the membrane even if the mRNA contains genes that encode proteins that ultimately reside in other cellular compartments. Thus, in this model, a PSORTdb annotation of ‘Cytoplasmic Membrane’, which we redefine for clarity as ‘Inner membrane’, would supersede other annotations. Conversely, as proteins that are found within the periplasm or outer membrane are post-translationally transported, a co-translational-insertion model for membrane RNA localization would not predict membrane enrichment for mRNAs that encode such proteins. To test this model, we modified the annotations associated with each polycistronic any operon that had at least one gene predicted to encode an inner-membrane protein had all other predicted locations for other genes discarded, whereas polycistronic operons that did not include a gene predicted to encode an inner membrane protein kept the full union of predicted locations in their annotations. While this model was applied to test the significance and enrichment of predicted encoded protein locations for the measured RNA localization clusters, all annotations are shown in fig. S10 for completeness. The significance and enrichment calculations were performed with the GOATOOLS package (93) using the Fisher exact test and a Benjamini-Hochberg false discovery rate (FDR) correction with an adjusted p-value threshold of 0.05 to define significant enrichments.
To determine whether an operon encodes at least one transmembrane domain, we determined whether any of its gene components contains a transmembrane annotation in the Ecocyc database. More specifically, the protein-features.dat file in the offline version of Ecocyc 27.1 was parsed to associate each protein (as named by the “FEATURE-OF” field) to the types of its features (“TYPE” field). A protein was considered to have a transmembrane domain if it possessed at least one feature with a “transmembrane-regions” type. To link the gene names to their corresponding proteins, we converted the protein names to Ecocyc gene names using the associations described in the protein-links.dat file and the Ecocyc names to the gene symbols as described in the gene-links.dat file. Both files were downloaded from Ecocyc 27.1.
To compute the 1D axial distributions of the mRNAs, the radial component of the intracellular coordinates of the identified RNA molecules was discarded, the axial component was mirrored to enforce axial symmetry, and a 1D Gaussian KDE with a bandwidth of 0.2 was computed with scipy for each mRNA with more than 500 identified molecules across all cells. To explore how RNA localization may vary with cell cycle, we used cell length as a proxy for division stage and grouped cells by length. First, the pre-expansion length of the cells was estimated by dividing the length of their segmentation mask by their dataset-specific linear expansion factor computed above. Then, the cells were grouped by pre-expansion cell length in 0.6 μm increments and the axial KDEs were computed for each group separately. Because the number of cells in each group was sometimes too low to accurately estimate density for individual mRNAs, the mRNAs were binned by chromosome position in 10° increments (~130 kbp), and KDE was performed as described above on the chromosome bins instead of individual mRNAs.
To determine the spatial distribution of mRNAs measured with smFISH in unexpanded samples, where the spatial resolution was often insufficient to resolve individual RNA molecules, we developed an alternative approach to estimating average RNA localization. Specifically, 3D smFISH images were deconvolved, z-plane by z-plane, with a 2D Lucy-Richardson deconvolution with a Gaussian kernel of 7 pixels and a sigma of 2 pixels using the routines provided by scikit-image. In parallel, cells were segmented with Ilastik on the rRNA images as described above. Cells were then filtered based on the solidity (≥0.8), volume (≥50 voxels), planarity (the skeletonized mask must be contained within two z-planes), and length (a skeleton volume ≥ 7 pixels) of the segmentation masks. To eliminate cells that might span multiple FOVs, masks within 50 pixels of a FOV side were discarded. To align cells, PCA was performed on the 3D mask pixel coordinates to identify the orientation of the cell, and an affine transformation was applied to rotate both the mask and the deconvolved smFISH into a new coordinate system in which all cells are centered and have a uniform direction for the long axis of the cell. As cell cycle is correlated with cell length, cells were sorted into different length groups based on the length of the cell and resized with scipy such that all cells in the group had the same dimensions. The intensity profiles of all cells within a length group were then averaged to produce the final estimate of the spatial distribution of each mRNA measured with smFISH. The presented spatial distributions correspond to the average of cells post division, i.e., lengths between 3.3–3.8 μm, and were shown for a mid-cell z-plane.
Because expression levels for even abundant PULs can be low (fig. S12) and because the variation in B. theta density led to regional differences in the frequency of segmentation artifacts (when using similar Ilastik protocols as described above), we adopted a cell-free spatial analysis for the distribution of B. theta gene expression within the mouse colon. Moreover, as the density of mRNA expression varied substantially across regions within the colon and there was variability in expansion between samples, we adopted an approach that allowed for adaptive spatial binning. Specifically, we generated local patches of gene expression with fixed counts and dynamic size. A spatial patch was defined for each measured mRNA as the RNA itself plus its 49 nearest neighbors. To remove low mRNA density regions dominated by false positives, any patch that contained two mRNAs separated by more than 30 μm was discarded. By design, the patch construction process created partially redundant patches that could share many of the same mRNA molecules. Since making a non-redundant patch set is a form of the optimal packing problem, for which no computationally efficient exact solution exists, we applied a greedy algorithm to create a non-redundant set. In this approach, patches were selected at random and kept only if all mRNAs within that patch were unassigned to any other kept patch. The mRNAs within a kept patch were then marked as kept, and the process was repeated until all patches were either kept or discarded.
To eliminate patches disproportionately enriched in false positives, patches were further cut on the number of blank counts (greater than 3), the presence of unusually dim molecules (average mRNA brightness below the 2^nd^ percentile for all patches), or the presence of unusually bright molecules (average mRNA brightness above the 96^th^ percentile for all patches). The position of the host mucosa was determined via the intensity of the DAPI signal for all samples, and the distance between each spatial patch and the mucosal boundary was measured. The rare spatial patch that mapped to the host tissue was also removed. Based on visual inspection and the average count per cell determined from well segmented B. theta, individual spatial patches typically contained ~3–7 B. theta cells.
To identify reproducible spatial variation in the expression of B. theta in the colon, the spatial patches determined from four samples measured from a mouse prepared with the methacarn fixation and two samples measured from a mouse prepared with PLP fixation were combined in a single anndata structure, and a nearest neighbor graph of 15 nearest neighbors was created using a cosine distance metric. This graph was used to create a UMAP embedding with a minimum distance of 0.2 and to compute diffusion components (DCs) using the default parameters in scanpy.
To identify operons that were up or down regulated at low or high DC1 values, patches were categorized as low, high, or intermediate DC1 by dividing DC1 into 3 equally sized bins. The expression of each operon was averaged across all patches within each sample and DC bin, and the sample-averages were then used to calculate the log2-fold change in expression between the high- and low-DC1 bins. A t-test was performed with these sample-averages to determine the significance of this difference. The reported p-values were corrected for multiple hypothesis testing with the Benjamini-Hochberg FDR method as implemented in scipy. An FDR-corrected p-value threshold of 0.05 was used to define significance. The annotations of the target of each PUL were sourced from dbCAN-PUL (94). In the rare instances in which a PUL listed multiple target polysaccharides with both host and dietary origins, we annotated the PUL as host-mucus-associated.
To explore the degree to which a simple spatial categorization could recapitulate similar trends, we categorized each spatial patch based on its distance from the host tissue. For each dataset, we selected a distance threshold that reflected the local enrichment in bacterial density near the mucus and the overall size of the sample. From this distance threshold, we categorized patches greater than this distance threshold as high-distance patches and those less as low-distance patches. As expected, we noted that the fraction of high- or low-DC1 patches within these categories reflected the partial covariation of DC1 value with spatial proximity. The same approach described above was used to calculate the log2-fold changes and significances. The gene annotations were drawn from KEGG (95), dbCAN-PUL (94), and CAZy (96).
MERFISH probe design and decoding were performed with public code (github.com/ZhuangLab/MERFISH_analysis) and were run in MATLAB r2021a. All other analysis was done in python 3.11.8 using scanpy (v1.9.8), anndata (v0.10.5.post1), harmonpy (v0.0.9), leidenalg (v0.10.2), umap-learn (v0.5.5), numpy (v1.26.4), scipy (v1.12.0), opencv-python (v4.9.0.80), scikit-image (v0.22.0), scikit-learn (v1.4.1.post1), tables (v3.9.2), tifffile (v2023.2.28), pandas (v2.2.1), and h5py (v3.10.0).