CytoMAP: A Spatial Analysis Toolbox Reveals Features of Myeloid Cell Organization in Lymphoid Tissues

Abstract

Recently developed approaches for highly multiplexed imaging have revealed complex patterns of cellular positioning and cell-cell interactions with important roles in both cellular- and tissue-level physiology. However, tools to quantitatively study cellular patterning and tissue architecture are currently lacking. Here, we develop a spatial analysis toolbox, the histo-cytometric multidimensional analysis pipeline (CytoMAP), which incorporates data clustering, positional correlation, dimensionality reduction, and 2D/3D region reconstruction to identify localized cellular networks and reveal features of tissue organization. We apply CytoMAP to study the microanatomy of innate immune subsets in murine lymph nodes (LNs) and reveal mutually exclusive segregation of migratory dendritic cells (DCs), regionalized compartmentalization of SIRPα^- dermal DCs, and preferential association of resident DCs with select LN vasculature. The findings provide insights into the organization of myeloid cells in LNs and demonstrate that CytoMAP is a comprehensive analytics toolbox for revealing features of tissue organization in imaging datasets.

Keywords: cellular organization; dendritic cell positioning; lymphoid tissue anatomy; machine learning; quantitative image analysis; quantitative microscopy; spatial analysis; tissue microanatomy; tissue organization; tumor microenvironments.

Figure 1.. Workflow and Features of CytoMAP

(A) CytoMAP is designed to extract quantitative information on cellular localization and composition within tissue regions, revealing how local cell microenvironments form global tissue structure, as well as allowing comparison of intra- and inter-sample tissue heterogeneity. (B) The workflow starts with multi-parameter imaging of either thin sections or large 3D tissue volumes. Next, hierarchical gating of cell objects is used to annotatedistinct cell subsets, which are passed into CytoMAP for analysis. CytoMAP segments these spatial datasets into individual neighborhoods and uses clustering algorithms to define similar groups of neighborhoods, or tissue “regions,” which are explored and spatially reconstructed in 2D or 3D space. (C) CytoMAP contains multiple tools to quantify and visualize the tissue architecture, including analysis of spatial correlations between different cell types,investigation of distance relationships of cells with architectural landmarks, analysis of neighborhood heterogeneity within individual tissues or across multiple samples, and quantitative visualization of tissue architecture.

Figure 2.. CytoMAP Identifies Major Features of LN Tissue Structure

(A) Confocal image and zoom-in image of a LN section from a C57BL/6 mouse, adoptively transferred with 10^6 naive OT-II CD4+ T cells and immunized with OVA plus alum. Overview image scale bar, 200 μm; zoom-in scale bar, 30 μm. Only select channels are shown. (B) Histo-cytometry plots of cell object MFI for different channels demonstrating the gating used for identification of the indicated immune cell populations. (C) Positional plot of cell data from (B) (area matches the zoom-in image in A). CytoMAP was used to calculate the number of cells in 30-μm-radius neighborhoods (denoted by the circles in the bottom left), which were raster scanned as denoted by the arrow. (D) Heatmap of the neighborhood composition (percentage of each cell phenotype per neighborhood) after SOM clustering. Individual clusters, or “regions,” aredenoted by the color bar at the top of the graph. Arrowheads at the bottom highlight specific neighborhoods. (E) Region color-coded positional plot of the neighborhoods from (D). (F) Pseudo-space plot with the neighborhoods sorted based on B cell composition (sorted to the left) and T cell composition (sorted to the right). (G) Dimensionality reduction plots of the neighborhoods in which the standardized numbers of cells and total MFI of all channels were used for the dimensionalityreduction. t-SNE, PCA, UMAP, and PHATE were all calculated for the same input neighborhoods, which are color coded based on region type from (D). For this experiment, an imaging volume of 0.03 mm³, 139,399 cells, and 11,328 neighborhoods were analyzed.

Figure 3.. CytoMAP Analysis of a Murine Colorectal Tumor Sample

(A) Confocal image of a representative 20-μm-thick CT26 tumor section isolated 9 days after subcutaneous inoculation. Zoom-in image shows a region in the tumor periphery. Main image scale bar, 500 μm; zoom-in scale bar, 50 μm. (B) Positional plot of the lymphocyte (top) and myeloid cell (bottom) populations. (C) Heatmap of the normalized immune cell composition of regions defined by SOM clustering of 50-μm-radius neighborhoods from all imaged samples. (D) Positional plot of neighborhoods color coded by region defined in (C). (E) Region prevalence plot showing the percentage of the neighborhoods from each sample in each region. (F) t-SNE plots of the neighborhoods from all three samples, color coded based on region type. (G) Pseudo-space plot with TAM-rich neighborhoods sorted to the left and DC-rich neighborhoods sorted to the right. (H) Heatmap of the Pearson correlation coefficients of the number of cells per neighborhood across the imaged tumor samples. For this experiment, 3 tumor samples with a total imaging volume of 1.8 mm³ and a total of 102,354 myeloid cells, 68,801 lymphocytes, and 192,785 neighborhoods were analyzed.

Figure 4.. CytoMAP Analysis of Mtb-Infected Lung Granuloma

(A) Confocal image of a 20-μm-thick section from a Mtb-infected murine lung sample. The left image shows multiple imaged channels overlaid in white. Scale bar, 500 μm. Zoom-in images show separately acquired regions of interest of an uninvolved lung area and the Mtb granuloma. Zoom-in scale bars, 100 μm. (B) Positional plots of the immune cell subsets in the uninvolved (top) versus granuloma (bottom) lung areas shown in (A). (C) Heatmap of SOM-clustered, 50-μm-radius neighborhoods. (D) Positional plots of the neighborhoods defined in (C), color coded by region, within the uninvolved (top) versus granuloma (bottom) lung areas. (E) Pseudo-space plots of the uninvolved (left) versus granuloma (right) lung neighborhoods after sorting for Alv. Macs (sorted to the left) and Mtb⁺ cells (sorted to the right). (F) Rotated half-heatmaps of the Pearson correlation coefficients of the number of cells within the neighborhoods across the uninvolved lung (left) and thegranuloma (right). For the uninvolved lung region, an imaging volume of 0.05 mm³, 36,194 cells, and 4,725 neighborhoods were analyzed. For the granuloma lung region, an imaging volume of 0.07 mm³, 140,453 cells, and 7,350 neighborhoods were analyzed.

Figure 5.. CytoMAP Reveals Diverse Patterns of Myeloid Cell Organization in LNs

(A) Confocal image of a representative steady-state, non-draining LN section from an immunized C57BL/6 mouse. Scale bar, 200 μm; right zoom-in scale bar, 50 μm. (B) Histo-cytometry gating scheme used to annotate cell subsets. (C) Positional plots of the myeloid cell subsets (D) Heatmap of the cellular composition of 30-μm-radius neighborhoods showing regions as defined by SOM clustering (top color bar) and regions defined by manual annotation (bottom color bar). (E) Positional plot of the LN neighborhoods as color coded by the region type (top color bar in D). (F) Positional plot of surfaces generated on the manually annotated regions (bottom color bar in D). (G) Violin plot showing the number of cells as a function of distance to the border of the TZ+IFZ annotated region for the representative LN sample. Distances tothe left of zero represent cells inside the region; distances to the right of zero represent cells outside the region. (H) Plot of the distances of cells to the TZ+IFZ region, in which each dot represents the mean distance of the indicated cell population in an individual LN sample (n = 5). (I) Mean distances of indicated cells to the CMC region. (J) Mean distances of indicated cells to the B follicle region. (K) Plot of the Pearson correlation coefficients between the number of indicated cells per neighborhood using either all tissue neighborhoods (left) or only the TZ+IFZ region neighborhoods (right). Correlations were averaged over the sample cohort (n = 5) from one experiment. For this experiment, 5 LNs with an average imaging volume of 0.023 mm³ and an average of 6,611 myeloid cells, 84,347 spots, and 12,708 neighborhoods were analyzed per sample. Data are representative of at least two independent experiments.

Figure 6.. CytoMAP Analysis of 3D LNs Reveals Preferential Proximity of DCs with Vasculature

(A) Representative confocal image of a 500-μm-thick, Ce3D-cleared, steady-state LN slice. Main image scale bar, 200 μm; zoom-in scale bars, 50 μm. (B) Violin plot comparing the distances of the indicated cells to the nearest blood vessels within a representative LN sample. (C) Mean distances of cell populations to the nearest blood vessel, with each symbol representing an individual LN sample. Circles and squares representsamples from two independent experiments. Cells below the dotted line at 20 μm are considered proximal to vasculature. (D) Percentage of proximal cells for the LN samples shown in (C). (E) Heatmap indicating the Pearson correlation coefficients, averaged across all imaged samples, between the number of different cell or landmark object typesper 50-μm-radius neighborhood. For this experiment, an average imaging volume of 0.93 mm³ and an average of 39,139 myeloid cells, 36,049 blood vessel objects, 98,279 spots, and 36,049 neighborhoods were analyzed per sample. Data represent four samples from two independent experiments.

Figure 7.. Individual Blood Vessel Branches in LNs Demonstrate Selective Associations with Specific DC Subsets

(A) Confocal image demonstrating the CD31-labeled vascular networks in the same representative LN sample as presented in Figure 6A. Scale bar, 200 μm. (B) Zoom-in of the region in (A) (denoted with a rectangle), indicating the centers of CD31⁺ blood vessel objects (white dots). The circles with yellow dots represent spherical-object-centered neighborhoods with a radius of 20 ¼m. Scale bar, 20 μm. (C) Heatmap of the myeloid cell composition for the vessel-centered neighborhoods after SOM clustering. The top color bar indicates color-coded regions ofneighborhoods. The bottom color bar indicates manually annotated regions. (D) t-SNE plot of neighborhoods from all samples color coded by region from (C). (E) Positional plot of the vascular neighborhoods in a 200-μm-thick virtual Z section, color coded based on the annotations in (C). (F) Region prevalence plot showing the percentage of the neighborhoods from each sample in each region. (G) Mean number of blood vessel objects in each neighborhood for each region showing the size of the vessels associated with each cell type. ****p < 0.0001, ascalculated by repeated-measures one-way ANOVA with multiple comparisons. Data represent the same four samples from two independent experiments used in Figure 6.