The study was performed under the approval of the Stockholm North Ethical Committee on Animal Experiments (permit number N397/13 and N487/11). The animals were euthanized using cervical dislocation.

C57BL/6 mice were purchased from Janvier labs. C3HeB/FeJ mice were received from Igor Kramnik (BU, Boston, MA). All mice were housed and handled at the Astrid Fagreus Laboratory, Karolinska Institutet, Stockholm, under specific pathogen-free conditions and according to directives and guidelines of the Swedish Board of Agriculture, the Swedish Animal Protection Agency and the Karolinska Institutet. For the infection assay, mice were infected with 150–200 Mycobacterium tuberculosis (Mtb) Harlingen strain by aerosol using a nose-only exposure unit (In-tox Products). At the indicated time after infection, mice were euthanized and lungs were extracted and fixed in 4% buffered paraformaldehyde for 24 h. Fixed left lungs were paraffin-embedded. From each lung sample, 8-μm sections (along the long axis of the lobe) were obtained and stored at −80°C. Slides were paraffin-removed and dehydrated directly before further processing.

ISS was performed as described in Carow et al. (13, 19). After deparaffinization and digestion with pepsin 100 μg/ml in 0.1 M HCl for 20 min at 37°C, sections were washed with dH₂O and PBS and then dehydrated through ethanol series. Hybridization chambers were mounted on top of the samples (SecureSeal Hybridization chambers, Grace Bio-Labs). Reverse transcription of mRNA was performed overnight at 37°C with the use of random decamer primers and specific primers that partially overlap with the sequence of the padlock probes in a mix containing 5 μM random decamers, 300 μM specific primers (IDT, Leuven, Belgium), 0.2 μg/μl BSA (NEB), 1 mM dNTPs (Blirt), 20 U/μl TranscriptMe reverse transcriptase (Blirt), and 1 U/μl RiboLock RNase Inhibitor (Blirt) in the TranscriptMe reaction buffer. The newly synthesized cDNA was crosslinked with 4% paraformaldehyde at room temperature for 45 min followed by degradation of mRNA and padlock probe hybridization and ligation performed in a mix containing 1× Tth buffer (20 mM Tris-HCl, pH 8.3, 25 mM KCl, 10 mM MgCl₂, 0.5 mM NAD, and 0.01% Triton X-100), 10 nM of each padlock probe, 0.2 μg/μl BSA, 0.5 U/μl Tth ligase (Blirt), 0.4 U/μl RNase H (Blirt), 50 mM KCl, and 20% formamide (Sigma). Samples were incubated at 37°C for 30 min and then 45°C for 45 min. Multiple padlock probes were designed for each of the 36 transcripts of interest. Each padlock probe set for a target gene contains a unique four-base barcode, which is used for identification of target genes. Rolling circle amplification mix included 1 U/μl phi29 (Monserate Biotechnology), 0.4 U/μl ExoI (Thermo Scientific), 0.25 mM dNTPs, and 5% glycerol in phi29 buffer [50 mM Tris-HCl, pH 7.5, 10 mM MgCl₂, and 10 mM (NH4)₂SO₄]. The samples were incubated at 37°C for 5 h and then at 30°C overnight, generating the rolling cycle product (RCP) for detection.

Sections were then washed and 0.5 μM Alexa750 anchor probe was hybridized to the amplification products for 45 min at room temperature in the dark in a mix containing 2×SSC and 20% formamide. For sequencing by ligation, 0.1 μM of each of the four interrogation probes conjugated with corresponding fluorophores (FITC, Cy3, Cy5, or TexasRed) were incubated for 1 h at room temperature in the dark with a mix containing 0.1 U/μl T4 ligase (Blirt), 1× T4 ligase buffer, 0.2 mg/ml BSA, 1 mM ATP (Thermo Scientific), and 100 ng/ml DAPI. After PBS washes, sections were dehydrated through ethanol series and air dried. Cover slips were mounted with SlowFade™ Gold anti-fade reagent (Invitrogen). Imaging was performed in a Axio Imager Z2 epifluorescence microscope (Zeiss) with 20× objective by acquiring Z-stacks of overlapping tiles that together cover the tissue section (10% overlap).

After imaging, the samples were prepared for the next sequencing cycle. For the removal of the interrogation probes, the samples were treated with 0.02 U/μl UNG (Thermo scientific) and 0.2 μg/ml BSA in UNG treating buffer for 15 min and washed twice with DEPC-PBS-T. Three washing steps with 100% formamide for 3 min followed. After washing with DEPC-PBS-Tw, the anchor hybridization, ligation, and imaging processes were repeated for the next base. A full list of padlock probes, reverse transcription primers, and sequencing-by-ligation oligonucleotides is available in the work of Carow et al. (13).

After ISS, the cover slips were removed from the glass slides, the tissue sections were washed in PBS, and TB Auramine–Rhodamine T (AR) staining was performed using the TB Fluorescent staining kit T (BD). Briefly, the tissue sections were incubated at 37°C for 15 min in pre-warmed TB AR (BD), washed with tap water, and decolorized for 4 min at room temperature with Decolorizer. After being washed with tap water, the tissue sections were stained with DAPI for 15 min in room temperature, washed with tap water, and air dried, and cover slips were mounted with SlowFadeTM Gold anti-fade reagent (Invitrogen). Imaging was performed in an Axio Imager Z2 epifluorescence microscope (Zeiss) with 20× objective by acquiring Z-stacks of overlapping tiles that together cover the tissue section (10% overlap).

The z-stacks of acquired images were merged to a single image using the maximum-intensity projection (MIP) in the Zeiss ZEN software. The resulting images were stitched into a single image containing the entire scanned area. The images from the four sequencing rounds were aligned based on the DAPI staining. For the analysis of the aligned images, a CellProfiler (v2.1.1) automated pipeline was used (Github repository, Blobidentification_NucleiCounting), which applied ImageJ plugins for image registration. A manual threshold for identification of primary objects was set between 0.0015 and 0.0020. The general stain from RCPs (anchor probe) was used for saving x and y coordinates and the fluorescence intensities for each RCP for each of the four barcode bases were saved in a.csv file and decoded using a MATLAB script (Github repository, InSituSequencing_1). For each RCP and sequencing round, the base with the highest intensity was assigned to the corresponding RCP and a fixed quality threshold of 0.35 that was defined as the maximum signal divided by the sum of all signals for that base was applied.

For signal visualization the InSituSequencing_1 MATLAB script was used to plot the transcripts on DAPI, hematoxylin–eosin (HE), or AR staining images.

AR staining images were used for the single bacteria identification. Fluorescence from FITC channel was subtracted from the Cy3 channel to reduce background fluorescence (ImageSubstraction MATLAB script) and the images were tiled with InSituSequencing_V3 MATLAB script (both in Github repository). The tiled images were used as input for a Cell profiler pipeline (IdentifyImagedBacteria) (v 4.0.6) for the identification of all imaged objects (Github repository). The intensity range for the input images was set between 0 and 0.5 except for sections in which we observed oversegmentation of background and then the threshold was set to 0.008/0.015–0.5. The selected thresholding method was Otsu 3-classes, which calculates the threshold separating the three classes of pixels (foreground, background, and middle intensity) by minimizing the variance within each class. Threshold correction factor was set to 0.4 and the size of the smoothing filter was set to 8. The generated .properties file and training set were used as input for the Cell profiler-Analyst (v 2.2) to classify single bacteria, bacterial clumps (when few bacteria are very close to each other and create one identifiable object), and background. A.txt file was generated containing 20 rules for optimal object classification (rules_Mixed_20, Github repository).

The obtained object classification rules were used in a final Cell profiler (v 4.0.6) pipeline with tiled AR images as input. The parameters used were similar to the Cell Profiler pipeline for the identification of all the objects (Github repository, FinalSingleBacteria).

The generated .csv files with the x and y coordinates of identified single bacteria and bacterial clumps were merged and used as input for the RelateBacteriaAndGene MATLAB script together with the (details) .csv file obtained from the MATLAB decoding step that contains the x and y coordinates of the ISS-identified transcripts from the same tissue (Github repository). From this analysis, a QT_details_ParentBact.csv file was obtained, which contained the distance of each transcript to its closest bacteria. Transcripts within given bacteria distances were extracted, and the frequency of each immune transcript among all immune transcripts of the given distance was used for analysis and normalized to the frequency of the transcript in the total slide.

AR staining images were re-sized to 20% of the original image and used for the bacterial cluster identification. The re-sizing allowed the process of one (smaller size) image per tissue section, avoiding the tiling that is needed for processing original files and thus not interrupting the recognition of fluorescence signal from cluster areas. Subtracted images created as described above were processed with an automated MATLAB script (Github repository, BacterialClusters_wGradientOutput_Areas) for the cluster identification and gradient border assignment. For the bacterial cluster identification, a manual threshold between 78% and 95% was set for optimal identification depending on how widespread the staining for each tissue was. A threshold of 90% means the top 10% pixels are considered as foreground. The identified clusters were numbered and clusters with area smaller than 7,000 pixels² were excluded from further analysis. A gradient of selected size 300 pixels (97,5 μm) was drawn around the selected clusters. We were able to identify the transcripts reads falling in each cluster area or surrounding gradient by using the x and y coordinates of each identified transcript, which was produced during decoding with the InSituSequencing_1 MATLAB script. Finally, the density of each transcript as the number of transcript reads divided by cluster area or gradient step area was used for further analysis. For comparisons between sections, the transcript density was normalized by a performance factor as the ratio of total number of identified immune transcripts and total section area.

Normalized transcript reads or frequencies were uploaded to ClustVis, a web tool for visualizing clustering of multivariate data using heatmaps and principal component analysis (https://biit.cs.ut.ee/clustvis/). ClustVis calculates principal components using one of the methods in the pcaMethods R package and plots heatmaps using the heatmap R package (version 0.7.7).

Student’s t-test and ANOVA were used as indicated with p < 0.05 marked as significant (*), p < 0.01 (**), p < 0.001 (***), and p < 0.0001 (****).

All relevant scripts for ISS image analysis, the single bacteria identification, and the bacterial cluster identification are available at https://github.com/Moldia/Tuberculosis-bacterial-analysis.

In this study, we extended the analysis of the distribution of 33 immune and 3 control transcripts in lung sections of Mtb-infected mice (13) with the focus of describing the immune environment of Mtb bacteria in different mouse infection models and at different time points post infection. We have generated new datasets in order to compare two different Mtb mouse models, C57BL/6 and C3HeB/FeJ, since they show different pathologies after Mtb infection: the former does not develop necrotic lesions after Mtb infection, while the latter is susceptible to Mtb and shows necrotizing granulomas (20).

By aligning the generated transcript signals to the HE staining of the analyzed lung tissue section, we found the predicted patterns of immune transcripts and that our control Cc10 transcript (coding for Clara cell 10-kDa protein) localized as expected in airway epithelium ( Figure 1A and Supplementary Figure 1A ). For base identification in the barcode reading of the sequencing rounds, we chose a quality score of 0.35, defined as the maximum signal divided by the sum of all signals for that base. We further identified comparable expression patterns and a strong linear correlation (r ² = 0.97) when using different ISS quality thresholds ( Supplementary Figures 1B, C ).

By exploring the distribution of transcripts to identify major immune cell subsets, we found that transcripts of macrophage-indicating Cd68 and of inducible nitric oxide synthase-2 (Inos) that are expressed in activated macrophages increased from 3 weeks post Mtb infection (wpi) to 8 and 12 wpi. These transcripts were absent in encapsulated granulomas from C3HeB/FeJ mice at 10 wpi ( Figure 1A ). Cd11c and Cd86 transcripts expressed by antigen-presenting cells showed a scattered distribution ( Figure 1A ). However, the frequency of both Cd11c and Cd86 significantly increased at 8 and 12 weeks as compared to 3 wpi ( Supplementary Figure 1D ). As expected, frequency and occurrence of T cell-related transcripts like Cd8a were elevated after 3 wpi. In C57BL/6, the distribution was even across the lung tissue, whereas in C3HeB/FeJ mice, they localized mainly to organized granuloma structures. Cd19 transcripts localized to lymphoid-rich histology structures ( Figure 1A ). These potential B-cell follicles formed at the later time points in C57BL/6 lungs, but were reduced in the lungs of C3HeB/FeJ mice, confirming our previously published results (13) in this new C57BL/6 ISS dataset.

To determine immune transcripts that localized in the proximity of Mtb, we performed AR staining for the analyzed sections ( Figure 1B ). Bacteria accumulated in high-density regions of Cd68 and Inos transcripts, as reported previously (13) ( Figure 1C ). As we were able to locate single bacteria and immune transcripts generated by ISS ( Figures 1C, D ), we aimed to describe the immune molecules in the proximal to single bacteria and bacteria accumulated in lesions by developing two different automated pipelines.

As we were able to locate single bacteria and immune transcripts generated by ISS ( Figures 1C, D ), we aimed to describe the immune molecules in the proximal to single bacteria and bacteria accumulated in lesions by developing two different automated pipelines.

To objectively and automatically detect the Mtb bacteria in the tissue, we applied a machine learning approach to identify all bacteria in lung sections, to exclude background noise, and to connect the bacteria coordinates to ISS data ( Figure 2A ). As input, we used the Cy3 channel of AR images that had been subtracted by the FITC channel of the same image to reduce the impact of autofluorescence. In a first Cell Profiler pipeline, objects in tiled images were detected. These objects were then loaded in Cell Profiler analyst and manually sorted for a correct classification into single bacteria, clumps (2–4 bacteria that are not distinguishable), and background. Using the fast gentle boosting algorithm in Cell Profiler analyst, we retrieved 20 rules based on object shape and intensity that were of more than 88% accuracy to correctly classify unknown objects ( Supplementary Figure 2A ). Implementing the identified rules in the initial Cell Profiler pipeline (Github repository, FinalSingleBacteria) allowed the detection of objects, retrieving the identification and coordinates for single bacteria and clumps of bacteria ( Figure 2B ).

Next, we used a MATLAB script (Github repository, RelateBacteriaAndGene) to measure the distance of each transcript to the nearest bacteria, and found distinct frequencies of immune transcripts at different distances from bacteria ( Figures 2C, D ). Results for single bacteria and bacteria clumps were similar, and therefore pooled data were analyzed ( Supplementary Figure 2B ). In a C57BL/6 lung section, 3 wpi showed the most distinct transcript frequencies of innate immune transcripts (e.g., Cd11b, Cd11c, Cd68, Inos, Il12p35, and Il6) at subcellular distances (0–3 µm and 3–10 µm) from Mtb, which we interpret as high expression within the infected cells or in a cell neighboring an infected cell ( Figures 2C, D ). By computing the ratio of a transcript frequency at a given distance to the transcript frequency across the whole slide, we estimate over/underrepresentation of a given transcript in the selected distance. As transcript frequencies varied between the time points and models, only comparisons between the distances within one time point/model were performed ( Supplementary Figure 1D ).

A significant increase of transcripts indicating activated macrophages was found at subcellular distances compared to longer distances in C57BL/6 lungs, indicated by Inos expression at all time points and for Cd68, Cd11b, and Tnf at some time points ( Figure 3A ), analyzing three sections per condition. In contrast, applying this analysis of single bacteria surroundings in C3HeB/FeJ lungs at 10 wpi did not show any preferential expression pattern of macrophage and activation markers. The elevated Cd11c transcript ratio at subcellular distances at 3 wpi could indicate the infection of Cd11c-expressing cells as monocytes, alveolar macrophages, or dendritic cells (DCs) by Mtb ( Figure 3B ). However, at 8 wpi, no differences between subcellular and longer distances were observed. At 12 wpi, an opposite pattern was found with a preferential Cd11c expression at longer distances. Transcripts for co-stimulatory molecules Cd80 and Cd86 involved in antigen presentation and T-cell activation as well as Th1 differentiation cytokine Il12p40 similarly displayed an increase with distance at 12 wpi. C3HeB/FeJ lungs showed a similar pattern with increased DC-related transcript ratios for longer distances to identified bacteria ( Figure 3B ). As expected, T cell-indicating transcripts were preferentially found for distances >10 µm ( Figure 3C ). However, the transcript ratio for Cd8a was increased at subcellular distances to Mtb in lungs of C57BL/6 mice at 3 wpi and of C3HeB/FeJ mice.

The single bacteria identification analysis gave a robust overview of the Mtb immune environment in the analyzed tissue; however, it might be a more relevant tool for studying the disease development in the early stages upon infection. That is, because in early stages of the disease, the bacteria are more spread throughout the tissue, compared to later time points when Mtb aggregates in lesions creating distinct microenvironments. Our above pipeline described an average immune environment around single bacteria or small clumps without the ability to display the environment’s heterogeneity. Thus, we aimed to automatically identify individual bacteria clusters corresponding to tuberculous lung lesions to analyze their associated immunogenic heterogeneity.

For this purpose, we developed a MATLAB script that used a 20% image of the Auramine staining (Cy3 channel subtracted by FITC channel as above) as this enabled us to analyze the whole tissue section without dividing images by image tiling. The bacteria cluster identification was based on the intensity of the AR staining applying a top-hat filter, smoothening and thresholding of the image using parameters that consistently gave a satisfactory separation (Github repository, BacterialClusters_wGradientOutput_with_areas) ( Figure 4A ; Supplementary Figures 3A, B ). The identified clusters were connected to the ISS coordinates resulting in extracted reads from each assigned cluster. To study areas within distance to the cluster, we additionally extracted ISS reads from areas in 100-µm steps starting from the cluster borders. As the ISS assay performance varies between sections, we normalized the transcript density per cluster to the density of immune transcripts per slide, which allowed the comparison between distances and time points.

The percentage of tissue area that was occupied by bacteria clusters increased from 10% at 3 wpi to around 30% at 8 and 12 wpi, showing the dissemination of bacteria in the lung tissue over time ( Figure 4B ). We compared the ratios of all expressed transcripts at 3 and 12 wpi in C57BL/6 lungs in the bacteria clusters containing the highest transcript numbers to ensure data robustness (5 highest at 3 wpi, 10 highest at 10 wpi). This analysis showed that any cluster from 3 wpi was distinct from 12 wpi, although some transcript ratios could overlap at the individual cluster level ( Figures 4C , 5 ). The range of ratios within one section could reflect the immune heterogeneity of bacteria clusters. Only 3 transcripts (Inos, Cd68, and Socs3) increased over time within the analyzed bacteria clusters ( Figures 5A, B ). We expected that the most relevant immune responses to influence the pathogen–host interaction would occur within the identified cluster and within the surrounding 100-µm radius ( Supplementary Figure 4A ). In line with our single bacteria identification, innate immune response-associated transcripts (Inos and Tnf) showed a strong decrease with an increased distance to the bacteria cluster ( Figure 5A ). The increase of Socs3 from 3 to 12 wpi within the bacteria cluster showed a preference for the bacteria cluster and lower levels in the cluster-surrounding 100 μm ( Figure 5B ). SOCS3 is a negative feedback regulator of cytokine and growth hormone signaling and has been shown to regulate the immune response against Mtb within several immune cell types (21, 22). However, transcripts of possible SOCS3-inducing cytokines—Il6 and Il10—did not follow the same pattern as Socs3 ( Figure 5B ). Most T cell-related transcript density ratios decreased from 3 to 12 wpi, which was unexpected ( Figure 5C ). Importantly, no differences between transcript ratios within the cluster and its 100-μm gradient could be detected at 12 wpi, indicating an even distribution of T cells in the tissue.

Next, we applied our bacteria cluster identification script to find and describe bacteria clusters in Mtb-infected lung sections of C3HeB/FeJ mice. After automated cluster identification, we manually grouped the clusters based on their underlying histopathology into three different lesion types that we called BOG (big, organized granulomas), IOG (intermediate, organized granulomas), and SC (small clusters) ( Figure 6A and Supplementary Figure 5A ). BOG and IOG bacteria clusters showed a round structure with a center of necrotic cells surrounded by a lymphocytic rim, whereas SC bacteria clusters were of uneven structure, mainly consisting of epithelioid cells ( Figure 6A and Supplementary Figure 5A ). Comparing the transcript ratios of all immune transcripts contained in the cluster types, we observed a clear distinction between BOG and SC, with IOG showing overlap within both groups ( Figure 6B and Supplementary Figure 5B ).

The grouping of bacterial clusters allowed a refinement of the single bacteria analysis that had shown lower transcript levels associated with innate-mediated protection compared to C57BL/6 mice. We found reduced levels of macrophage and dendritic cell-associated transcripts (Cd68, Inos, Tnf, and Cd11c) in BOG compared to SC ( Figures 6C, D ). T cell-associated transcripts showed higher expression levels of Tcrb and Cd8b1 in SC and increased Cd8a, Ifng, and Foxp3 levels in BOG ( Figures 6D, E ). Cd4 transcript levels were comparable between BOG, IOC, and SC ( Figure 6E ). Most of the differences in transcript levels were very localized as they were only present within the bacteria clusters and transcript patterns of BOG and SC overlapped when applying the 100-μm radius around the identified clusters ( Figures 6C–E ; Supplementary Figures 4B , 5B ).

Finally, to illustrate differences in localized immune environments of Mtb in lungs of C57BL/6 and C3HeB/FeJ mice considering cluster types and time points after infection, we compared expression patterns in the bacteria cluster of C57BL/6 mice at 3 and 12 wpi with BOG and SC in C3HeB/FeJ mice at 10 wpi ( Supplementary Figure 6 ). A principal component analysis of normalized transcript densities placed the SC of C3HeB/FeJ mice between 3 and 12 wpi bacteria clusters of C57BL/6 mice ( Figure 7A ). BOG clustered distinct from SC and from C57BL/6 clusters. The expression levels of Cd68, Inos, and Tnf were reduced in BOG clusters not only compared to SC but also in comparison to clusters in early and late time points in Mtb-infected C57BL/6 lungs ( Figure 7B ). In contrast, the BOG characteristic elevated Cd8a and Foxp3 transcripts were found at increased levels compared to SC but also in comparison to 3 and 12 wpi bacteria clusters of C57BL/6 mice ( Figure 7C ). Altogether, this “cross-comparison” further emphasizes the compartmentalization in two different bacteria environments in C3HeB/FeJ mice.

Recent developments in the field of spatial transcriptomics and metagenomics have been focusing either on the identification of multiple microbial species in situ (23, 24) or on the identification and transcriptional analysis of selected microbial species (25). Other approaches have associated the local expression profile with the manually annotated histopathological features of the tissue or the expression of specific disease markers (26).

Immunological studies of the TB granuloma have applied several spatially resolving techniques like laser caption microscopy in combination with proteomics as well as multiplexed ion beam imaging by tune of flight (MIBI-TOF) on human tissue sections containing active TB granulomas to find that granuloma structure and immune cell function are connected (27, 28). McCaffrey et al. imaged 37 proteins in TB granulomas constructing an atlas identifying 19 cell subsets and 8 spatial microenvironments. Importantly, one identified microenvironment indicated immune dysfunction characterized by IFN-γ depletion and an enrichment for TGF-β, regulatory T cells, and IDO1- and PD-L1-expressing myeloid cells. Their analysis increases our understanding of granuloma immunology and suggests key events leading to active TB, which opens for new therapies. However, all the aforementioned approaches lack the association of the microbial presence with the local immune profile of the tissue.

Mtb is unique with regard to several aspects of pathogen–host interactions, the formation of granulomas, a non-sterilizing T cell-mediated immunity, and a bacterial persistence program as some characteristics. As one step in understanding local immune responses to Mtb, we present automated analysis pipelines to combine the localization of Mtb bacteria and immune transcripts, which allows the definition of differences in the environment that is proximal to bacteria. The single bacteria detection pipeline used allowed us to study the immune response locally in selected distances from Mtb in infected mouse lungs. This approach proved valuable in the early stages of infection and in sections of tissue where Mtb is spread or scarce, and where only frequencies of immune transcripts in given distances could be calculated (no normalization to area possible). In later stages of Mtb infection, the bacteria form clusters corresponding to histopathological lesions or tuberculous granulomas depending on the murine genetic background. With the cluster analysis, we were able to identify the bacteria clusters according to the intensity of the AR Mtb staining and the adjacent areas allowed to calculate the immune transcript density. Combining results gained from the single bacteria and the bacteria cluster analysis gave insights into different aspects of the occurring pathogen–host interactions. Regarding our results, the presentation of one sample per time point (in triplicate) limits the conclusions that can be drawn, but the robustness and multiplexity of the generated data still allow us to generate suggestions that may be evaluated in extended studies.

Mtb has developed mechanisms to escape the microbicidal mechanisms of macrophages and can instead survive and replicate within them (29). We found a high expression of Cd68 and Inos at subcellular distances to bacteria as well as within bacteria clusters in C57BL/6 mice lungs. Inos and Tnf levels surrounding infected cells and clusters were also elevated. Similarly, we found elevated Cd68 and Inos levels in small bacteria clusters (SC) of C3HeB/FeJ mice but not in their organized granulomas. No subcellular increase for those transcripts was found in the single bacteria analysis for C3HeB/FeJ mice, probably since the average Mtb surrounding represents a mixture of SC and BOG. The expression of Inos at subcellular levels is indicative of classically activated M1 macrophages upon infection of the cell. However, as no M2 markers were included in the panel, we cannot compare the M1/M2 differentiation for infected cells and the Mtb immune environment, which remains to be studied in the future.

Inos encodes for the enzyme producing bactericidal nitric oxide (NO) and is central in the innate defense against Mtb in mice, and its absence leads to mice that are highly susceptible to Mtb (30, 31). The frequency of Inos mRNA at subcellular distances to Mtb increased from 3 to 8 and 12 wpi as well as the lung area occupied by bacteria clusters. We found the highest expression of Inos mRNA within the infected cells but not in neighboring cells, excluding a direct transcriptional repression by Mtb. Despite its correct cellular location and high expression, there are several possibilities that may account for the insufficient Mtb control: Toxic NO levels may not have been reached due to limitation of the iNOS substrate L-Arginine or insufficient enzymatic activity of iNOS. Another possibility is that Mtb’s evasion mechanisms prevented the bacterial clearance. Mtb has been shown to produce proteins that repair NO-damaged bacterial proteins at the proteasome level (32–34) and inhibits the recruitment of the enzyme iNOS to the phagosome (35).

Another transcript density that increased over time in C57BL/6 bacteria clusters besides Cd68 and Inos was Socs3. SOCS3 is a negative feedback regulator of cytokines of the IL-6 family, and several growth factors and hormones and SOCS3 in myeloid cells contribute to defense against Mtb (21, 22). As transcripts of its main inducing cytokines, IL-6 and IL-10 did not increase within bacterial cluster over time, and we speculate that other stimuli like bacterial products may contribute to the increase in Socs3 expression (36). We have previously shown that Socs3 expression induced by Mtb in macrophages is dependent on NF-κB and thereby at least partially Il-6 and Il-10 independent (21). It is worth mentioning that the local expression of immune transcripts in the bacterial cluster ( Figure 5 ), as measured here by ISS, can differ from the bulk measurement of cytokines, e.g., in lung homogenates. Therefore, comparison of ISS data with a measurement of cytokine levels in the lungs will not automatically match our findings but would be an interesting follow-up.

CD11c is commonly used to identify DCs but can also indicate the presence of alveolar macrophages, monocytes, and immature macrophages. DCs can take up Mtb antigens but also get infected by bacteria. Alveolar macrophages reconstitute a significant proportion of macrophages within the Mtb-infected mouse lung (37). They are the first cell contact of Mtb with the host after inhalation and contribute to Mtb dissemination, which matches our finding of an enriched Cd11c expression at subcellular distances at 3 wpi but not later (1). At 12 wpi, a transcript pattern indicating activated and antigen-presenting cells (Cd80, Cd86, and Il12p40) was preferentially found in uninfected cells. This is in line with previous observations showing a negative spatial association of CD11c+ cells with Mtb bacteria in mouse lungs in contrast to a positive correlation of CD11b+ cells indicating macrophages in neighborhood analysis (14, 38). It remains to be investigated whether these non-infected CD11c+ cells contribute to bacterial control. Interestingly, elevated levels of Cd11c and transcripts indicating dendritic cell activation (Cd86 and Il12p40) were found in SC in contrast to BOG of C3HeB/FeJ mice confirming their functional diversity.

It has been long questioned if a granuloma structure enables the direct contact of primed T cells with infected macrophages as required for their activation. Recent studies redirected the focus of T-cell functionality in, e.g., lung homogenates to the importance of the localization of T cells close to Mtb-infected cells to execute their function (16, 17, 39). In the single bacteria pipeline, we expected a preferential accumulation of T-cell transcripts especially indicating a differentiation of macrophage activating CD4 Th1 cells (e.g., Cxcr3, Ifng, and Tbx1) at a 10- to 30-µm distance compared to further distances. This would have indicated a functional immune interaction between an antigen-specific T cell and an Mtb-infected macrophage. However, we did not confirm this hypothesis. Instead, for 8 and 12 wpi, an increase in T-cell transcript frequencies for 10–30 µm and 30–100 µm correlated with the onset of adaptive immunity without a localized accumulation. Contradicting our expectations was the reduced T-cell transcript density at 12 wpi vs. 3 wpi within the bacteria cluster. The BOG formed in C3HeB/FeJ mice resembles human necrotic granulomas despite their susceptible genetic phenotype and will be discussed in detail below.

C3HeB/FeJ mice are very susceptible to Mtb infection despite having a functional immune system with unimpaired CD4 Th1 cell differentiation and macrophages that respond to IFN-γ with NO production (40). Their survival is shorter than iNOS-deficient mice and susceptibility has been mapped to the sst1 locus and the interferon-inducible gene Ipr1 that switches the mechanism of cell death in macrophages (41, 42). C3HeB/FeJ mice are used as a TB infection model as they develop necrotic granulomas not seen in C57BL/6 mice, and this severe pathology within a generally immunocompetent mouse resembles in that aspect human TB (20). The formation of necrotic granulomas in C3HeB/FeJ mice is not a consequence of high bacteria loads, as it has been shown to precede the strong increase in bacteria amplification. The comparison of bacteria containing parts of necrotic granulomas (BOG) that were identified by the bacteria cluster pipeline to the non-necrotic small bacteria cluster (SC) showed a different immune profile within the cluster types. It also confirmed the ability of macrophages to induce Inos and Tnf. Our analysis clearly demonstrates the compartmentalization of immune parameters looking at SC and BOG, respectively ( Supplementary Figure 6 ). However, it is important to bear in mind that the mechanisms leading to necrotic granulomas without control of bacteria growth in C3HeB/FeJ mice may differ from local immune events leading to the progression from latent to active TB in humans. BOG clusters were characterized by an increased expression of Foxp3, Cd8a, and Ifng transcripts previously identified as characteristic for necrotic centers in these mice (13). It is uncertain if Cd8a indicates the presence of T cells as Tcrb and Cd8b levels are low and Cd8a may also be expressed by dendritic cell subsets in mice during infections (43). As previously indicated, the low signal density did not allow a cell phenotyping with several markers detected on one cell, and the absence of additional T-cell transcripts emphasized this uncertainty.

Here, we used AR staining, which is an acid-fast fluorescent staining, a diagnostic method used as an alternative to classical Ziehl-Neelsen staining. Problems with AR stainability of bacteria in the necrotic core of C3HeB/FeJ mice at 10 wpi have been reported previously, as well as a reduction of general acid-fast staining of Mtb during dormancy in humans; a change in phenotype/cell wall composition was suggested as a reason for this phenomenon (44, 45). In our study, Mtb bacteria are stainable in all C57BL/6 sections and in the necrotic center at 10 wpi in C3HeB/FeJ mice, and an explanation could be that the pepsin treatment of the tissue before the ISS assay made the bacteria accessible for the staining. As a long-term goal, we would like to combine the detection of immune transcripts with mycobacteria transcripts in our ISS assay, to detect transcripts indicating viability, stress responses, and replication of Mtb in connection to their immune environment. The technical challenge in that approach is to digest the mycobacteria cell wall to make mycobacterial RNA accessible for detection probes without damaging the surrounding tissue and degrading host RNA.

Our study is based on a padlock probe panel against 33 immune transcripts of selected immune cell markers, activation molecules, and cytokines, which were all identified in parallel within one section. Although cell typing would likely reveal even more details of the immune environments, our focus was to study the expression levels of immune transcripts, in relation to the Mtb bacteria localization in the infected mouse lung tissue and not doing cell typing per se. Cell typing would require a higher density of transcript signals as achieved in our samples in combination with a much larger padlock probe immune panel, in order to base cell identity in multiple marker reads per cell. The nature of our samples, which are formalin-fixed and paraffin-embedded (FFPE), a biosafety requirement for the collection of Mtb-infected tissues, shows lower signal densities compared to fresh frozen samples. Methodological advancements have been achieved since the beginning of our study, e.g., the direct RNA detection for ISS (46), that would allow for increased transcript detection and subsequent cell typing in the future. However, we want to emphasize that ISS has been shown to be highly reproducible and representative for transcripts even at low signal densities (13, 26, 47) and that several transcript probes have been confirmed by immunohistochemistry in our previous publication (13).

We demonstrate here the advancement of combining ISS data (identified transcripts as coordinates) and bacteria localization. ISS was conducted first and the same slide was used for bacteria (AR) and HE staining. The fact that all the experimental procedures are conducted on the same slide allows us to acquire the exact coordinates of transcripts or bacteria localization and connect them with the tissue histopathology. However, if coordinates for identified cells by immunohistochemistry are available, these can be used as input in the pipelines as well. High multiplexed immunohistochemistry could not be performed on the same tissues in our case, since the processing with multiple experimental techniques (ISS, AR, and HE staining) could have possibly affected the structure of several antibody epitopes.

Altogether, the here-presented tools constitute a toolbox for analysis of biologically complex tissues and have the potential to increase our understanding of the localization of immune components with respect to Mtb. The advantage of using formaldehyde-fixed paraffin-embedded tissues allows the use of archived human biopsies. As padlock probes are custom-designed based on target sequences, this analysis is open to all species with known target sequences.

Aiming for the translational goal to not only understand the local immune control of Mtb in mice but also predict disease progression and status in human TB patients by investigating the peripheral blood, we suggest the following next steps. Using the here-presented approach, Mtb immune environments correlated to Mtb control or progression could be defined by studying non-human primate (NHP) samples that show a spectrum of latent and active granulomas and stainable Mtb bacteria. NHPs mirror the morphology and physiology observed in human TB disease, and immune correlates of progression in lungs of NHP could be correlated to immune correlates in the blood of TB progressors (48, 49). In addition, studies comparing BCG-vaccinated and non-vaccinated animal groups would allow the definition of differences in the Mtb environment induced by vaccination and correlated to bacterial control. Examples of human granuloma samples could be applied to confirm findings of identified expression patterns in Mtb immune environments in the different conditions. In vitro generated granulomas would give the opportunity to manipulate and challenge the effect of identified protective expression patterns. Although important progress in the in vitro granuloma field has been achieved, generating, e.g., three-dimensional spheroid human granulomas forming single organized structures consisting of human lung-derived alveolar macrophages surrounded by layer of autologous T cells to study innate and adaptive stages of the TB granuloma, the reduced complexity of in vivo granulomas must be considered (50, 51). Lastly, a transcriptomic meta-analysis of peripheral blood from patients with active and latent TB could correlate identified patterns of Mtb’s immune environments during disease progression to systemic measurable parameters.

Overall, our here-presented tools allowed for spatial analysis of transcript expression according to the tissue histopathology and the causative agent (Mtb) localization, which enabled us to compartmentalize the disease lesions and their distinct immune environments. As the identification of bacteria cluster is based on fluorescence intensity, the analysis pipeline is not limited to Mtb infections. It could be applied in other infectious diseases in which pathogens can be visualized in tissues either by direct staining or by antibody- or hybridization-based detection, enhancing our understanding of tissue and disease heterogeneity. Another application could be in the cancer field, in which a combination of transcripts identified by ISS with areas overexpressing tumor antigens like c-myc and p53 detected by immunohistochemistry would lead to a better definition and understanding of tumor microenvironments.