3D deconvolution of human skin immune architecture with Multiplex Annotated Tissue Imaging System

Routine clinical assays, such as conventional immunohistochemistry, often fail to resolve the regional heterogeneity of complex inflammatory skin conditions. We introduce MANTIS (Multiplex Annotated Tissue Imaging System), a flexible analytic pipeline compatible with routine practice, specifically designed for spatially resolved immune phenotyping of the skin in experimental or clinical samples. On the basis of phenotype attribution matrices coupled to α-shape algorithms, MANTIS projects a representative digital immune landscape while enabling automated detection of major inflammatory clusters and concomitant single-cell data quantification of biomarkers. We observed that severe pathological lesions from systemic lupus erythematosus, Kawasaki syndrome, or COVID-19–associated skin manifestations share common quantitative immune features while displaying a nonrandom distribution of cells with the formation of disease-specific dermal immune structures. Given its accuracy and flexibility, MANTIS is designed to solve the spatial organization of complex immune environments to better apprehend the pathophysiology of skin manifestations.


INTRODUCTION
The skin acts as a barrier organ that separates the body from the external environment. Upon inflammation, blood-circulating immune cells are recruited to help orchestrate the cutaneous immunity and are often nested nearby key structural elements (e.g., postcapillary venules, hair follicles, dermal-epidermal junction, etc.) (1,2). In pathological settings, the nature and activation status of the skin immune landscape often represent precious biological information that can help establish an accurate diagnosis, apprehend interpatient heterogeneity, and select the most appropriate treatment. The use of imaging-based approaches to identify cutaneous immune cells is still challenging because of the high level of autofluorescence arising from the tissue itself, the potential spectral spillover when more than four fluorochromes are used simultaneously, and the entanglement of all cells within thick and polarized structural appendages.
The vast majority of microscopic diagnoses of inflammatory skin conditions relies on repeated immunohistochemistry analysis of one or two proteins and/or hematoxylin and eosin (H&E) staining in thin (2 to 5 μm) formalin-fixed, paraffin-embedded (FFPE) specimens (3,4). While such two-dimensional (2D) approaches are reproducible and suitable for routine practice, they do not permit to apprehend the complex topology and heterogeneity of immune cells (5), particularly those nested in-between epidermal appendices. The development of image-based histo-cytometry, which consists of analyzing segmented multicolor images with classical flow cytometry gating strategies, has paved the way toward the development of sophisticated image generation systems coupled to computational imaging (6). Recently, highly multiplexed imaging systems have substantially advanced our understanding of tissue-resident immune subsets and of their spatial distribution with regard to tissue structures, with a strong focus on cancer samples and tumor heterogeneity, such as CODEX (Co-detection by indexing) (7,8), MIBI-TOF (Multiplex Ion Beam Imaging by Time Of Flight) (9), IMC (Imaging Mass Cytometry) (10), MuSIC (Multiplexing using Spectral Imaging and Combinatorics) (11), CyCIF (tissue-based cyclic immunofluorescence) (12), Cell Dive (13), and others (14). While multiplexed imaging has an immense potential, there is a strong need to democratize these methods with the use of inexpensive instrumentation compatible with standard tissue processing and coupled to an analysis interface that is userfriendly enough to be used in routine practice.
Here, we present an integrated framework primarily designed for spatially resolved immune cells phenotyping in FFPE human skin biopsies. We first set up a simple and inexpensive method to acquire 10 fluorescent signals simultaneously and in 3D using a classical confocal microscope. We next designed MANTIS (Multiplex Annotated Tissue Imaging System), an adaptable and interactive analytical system that automatically generates a digitalized version of the skin immune landscape and enables single-cell quantitative data visualization. On the basis of these settings, MANTIS could be implemented in most laboratories coupled to existing confocal equipment to bridge the gap between sophisticated research tools and standard-of-care diagnostic procedures with minimal human intervention.

Extraction of single-cell statistics from skin sections by combining conventional confocal laser scanning microscopy and computational imaging
We developed a simple method to generate 3D multiplexed fluorescent images from FFPE (10-μm-thick) or fresh frozen (50-μmthick) skin biopsies that could be implemented in most research or clinical laboratories on existing equipment. Skin sections were first stained with different panels of commercially available fluorochrome-coupled antibodies added simultaneously and then quenched to avoid excessive natural autofluorescence of skin structural elements (Fig. 1A). We acquired 3D fluorescent multiplexed images with a conventional inverted confocal laser scanning microscope equipped with five laser lines, five detectors, and a 40× oil immersion objective, using a strategy of sequential acquisition composed of fast consecutive steps ( Fig. 1B; the detailed description of optical paths and lasers of our 8-year-old Leica SP8 system is provided in Materials and Methods). This setting enabled the acquisition of 8 to 10 fluorescent channels, on a system primarily designed for four colors, over a skin section of the following three dimensions, 0.6 mm (x) by 0.4 mm (y) by 20 μm (z), within 25 min. The obtained 3D images were then deconvoluted and compensated to correct 3D fluorescent spectral spillovers (Fig. 1, C and D) using the Huygens software (Scientific Volume Imaging), a strategy routinely applied in flow cytometry to combine multiple fluorochromes simultaneously (6,15). Compared to classical segmentation strategies based on nucleus expansion (16), which often lead to under-or overestimation of cellular cluster composition, we used the general immune biomarker CD45 as a robust immune staining visualized in most skin-resident immune cells to constitute the core of our cell segmentation strategy for future single immune cell statistics extraction ( Fig. 1, E and F). Using the isosurface algorithm of the Imaris software (Bitplane), we next modeled the 3D fluorescence signal of CD45 for individual immune cells and exported a corresponding single-cell database composed of the mean fluorescence intensity (MFI) of all individual biomarkers and precise x, y, and z tissue coordinates obtained with a resolution of 299 nm by 299 nm by 999 nm per voxel (Fig. 1F). We found that CD45-based segmentation enabled an efficient isolation of single immune cell characteristics, even when those were found aggregated around dermal structural elements. Overall, we demonstrate that it is possible to extract a 10-parameter single-cell database using regular confocal equipment coupled to basic computational imaging steps.

Analysis of the skin immune landscape using MANTIS phenotype attribution matrices
On the basis of the literature, we designed two panels composed of antibodies directed against immune biomarkers suitable to generate a non-exhaustive overview of lymphoid cell (LC) and myeloid cell landscape of the skin, with an average cost of approximately $65 per sample. The combination of CD45, CD3e, CD4, CD8, γδ T cell receptor (TCRγδ), CD20, and CD57 [a terminally sulfated glycan carbohydrate epitope shared by natural killer (NK) and T cells with high cytotoxic potential (17,18)] allows us to identify the following LCs: conventional CD4 and CD8 T cells (being CD57 low or CD57 high ), CD4 + CD8 + double-positive (dp) T cells (19), CD4 − CD8 − double-negative (dn) T cells (20), γδ T cells, B cells, and NK cells (table S1). The combination of CD45, CD207, CD1c, human leukocyte antigen DR [HLA-DR (Human Leukocyte Antigen -DR isotype)], CD123, Siglec8, myeloperoxidase (MPO), and tryptase allows us to identify the following myeloid cells: Langerhans cells, Langerin + (CD207 + ) dermal dendritic cells (dDCs) and Langerin − dDCs, eosinophils, basophils, neutrophils, and mast cells (table S1). The activation status of DC, Langerin + DC, and LC was investigated using levels of HLA-DR expression. A detailed list of excitation/emission/detection strategies is provided in table S2.
We next aimed to develop an adaptable analytical system that could integrate and batch-process extracted single-cell databases and enable an unsupervised phenotyping of immune subsets. To address this latter challenge, we developed MANTIS, an interactive digital tool based on phenotype attribution matrices inspired by the analytical logic of single-cell RNA sequencing that identifies correlations between the single-cell database and the expression profiles of different cell types ( Fig. 2A and fig. S1A). Such an analysis is possible by computing Spearman's ρ correlation, which accommodates nonlinear relationships in the expression values (i.e., in our case, the collected MFI of each biomarker). In practice, MANTIS runs instantaneously a pairwise Spearman's correlation analysis, for each detected single-cell, against selected combinations of biomarkers to identify the immune subsets annotated in the phenotype attribution matrices. The output information is the attribution of specific ρ values per single cell, which then automatically finds the best match of cellular identity and generates associated quantitative statistics (Fig. 2, A to C, and fig. S1, B and C). When analyzed side by side with conventional flow cytometry, we found that the same populations were detected; however, MANTIS enabled the identification of more LCs and CD4 T cells (and less dn T cells), while flow cytometry enabled the identification of more HLA-DR high dDCs and eosinophils ( fig. S2). As a proof of concept, we generated data from two serial sections of an acral lesion from a patient with systemic lupus erythematosus (SLE; i.e., lupus chilblains) stained with a lymphoid and a myeloid panel. The fast 3D acquisition of one region of interest (ROI) composed of six fields of view [i.e., 0.6 mm (x) by 0.4 mm (y) by 20 μm (z)] enabled the annotation of 519 myeloid cells and 708 LCs for a total of 19 different immune subsets identified ( Fig. 2D and fig. S1D). One can then decide to visualize annotated immune populations using either a heatmap, in which the MFI of individual biomarkers is displayed per single cell ( fig. S1E), or a graph-based dimensionality reduction, i.e., t-distributed stochastic neighbor embedding (t-SNE), specifically designed for visualizing clusters of populations and corresponding expression of biomarkers per cluster (Fig. 2E).
A particularly challenging aspect of multiplexed imaging technologies is to circumvent the spatial distribution of immune cells with respect to longitudinal and polarized structural elements (e.g., epidermal appendages of the skin) in thick tissue sections. We developed an interactive software interface that contextualizes the immune topology of the skin by replacing all annotated single immune cells within their 3D spatial context and leverages the natural autofluorescence of keratinocytes to model the epidermal layer to facilitate biopsy orientation ( Fig. 2F and fig. S1, F and G). The algorithm allows us to use two complementary analytical approaches and to switch from one to the other with a simple drawing tool (movie S1). The analysis can start from the visualization of the skin digital immune landscape and then be pursued with the investigation of the immune composition of defined microregions via instantaneous recomputation of drawn ROI. Conversely, it is also possible to start from all annotated immune cells on a t-SNE graph, draw around subsets of interest, and immediately visualize their anatomical distribution in the skin digital immune landscape ( Fig. 2G and movie S1). Together, these data suggest that the MANTIS interactive analytical system can be used to compute the 3D spatial organization of immune and structural elements of inflammatory skin samples from patients.

Quantitative validation of MANTIS annotation system using healthy-looking skin and inflammatory pathological lesions
With the constant increase in the number of cases, a large panel of putative skin manifestations of coronavirus disease 2019 (COVID-19) have been observed worldwide (21,22), including an unprecedented high rate of acral lesions, which represent~75% of all cases and commonly named "COVID-toes" (23)(24)(25)(26). Such manifestations ( fig. S3A), compared to non-inflamed healthy-looking skin ( fig.  S3B), are associated with an important immune cell infiltration ( fig. S3C) and tend to develop in young patients with no or very mild respiratory symptoms (26,27). While some pathological features of those lesions have been described (3,28,29), a precise analysis of their spatial immune profile is currently missing, which impairs the development of a clear readout to better diagnose and treat these rare cutaneous lesions. A possible explanation could be a collateral clinical manifestation of an efficient antiviral type 1 interferon response because acral lesions are also commonly observed in patients with interferonopathies, such as the Aicardi-Goutières syndrome (30) and SLE (31). With this in mind, we decided to benchmark the effective performance of MANTIS to resolve the immune topology of skin lesions of similar clinical severity from five patients with COVID-toes, two patients with the multisystem inflammatory syndrome, which is clinically similar to Kawasaki syndrome (i.e., a rare severe systemic inflammatory condition triggered by severe acute respiratory syndrome 2 infection, named hereafter "Kawasaki syndrome"), and three patients with SLE chilblains. Abdominal skin biopsies from five healthy-looking controls were used to set the baseline of a natural steady-state immune environment, albeit from a distant anatomical region.
We validated the quantitative performance of MANTIS to annotate immune cells by calculations of statistical correlations with a supervised approach of histo-cytometry (6, 15) applied on the same datasets for each antibody panel in all skin samples. This last method consists of a manual gating of immune subsets on the same principle used in traditional flow cytometry. A total of 20,464 single CD45 + immune cells were identified with the following distribution per condition: 1670 immune cells in five healthy-looking skin samples (i.e., with 895 myeloid cells and 775 LCs), 1703 in two patients with Kawasaki syndrome (i.e., with 932 myeloid cells and 775 LCs), 5076 in three patients with SLE chilblain (i.e., with 1560 myeloid cells and 3516 LCs), and 12,015 in five patients with COVID-toes (i.e., with 2838 myeloid cells and 9177 LCs). A classical gating strategy based on mutually exclusive biomarkers was used to assess the presence of myeloid cell ( fig. S4A) and LC ( fig. S5A) subsets by histo-cytometry. We identified a total of 19 different immune subsets and found very similar distributions of cell counts by either supervised histo-cytometry or unsupervised MANTIS algorithm (figs. S4, B to E, and S5, B to E). The calculated R coefficients were between 0.75 and 1, regardless of the antibody panel, the patients analyzed, or the disease (figs. S4F and 5F). We observed that all healthy-looking skin samples exhibited a proportion of immune cells aligned with previously described skin-resident immune populations at steady state in human (1,32). However, we noted a slightly higher tendency to detect rare populations of blood-circulating CD45 + CD3 − CD20 + B cells or CD45 + HLA-DR − CD123 + Siglec8 + basophils, only when 3D images were computationally analyzed with MANTIS (figs. S4B and 5B). This is consistent with the fact that the skin is a highly vascularized tissue and that recent studies identified rare B cells in healthy skin (33).
Having validated the quantitative and qualitative performance of MANTIS-based annotation, we next defined a high-level view of the complex immune environment of pathological lesions from all 10 patients. Compared to healthy-looking samples, pathological lesions contained large immune infiltrates, confirming their inflammatory status (Fig. 3, A and B). All three conditions were associated with an infiltration of myeloid cells composed of a large number of neutrophils, eosinophils, mast cells, and conventional CD45 +-CD1c + CD207 − HLA-DR + dDCs (Fig. 3, A and C). While detected in relatively low numbers in all analyzed skin samples, no difference was observed between healthy-looking and pathological samples for CD45 + CD1c − CD207 + HLA-DR + LC or CD45 + CD1c + CD207 + HLA-DR + dDC populations.
Compared to patients with Kawasaki syndrome, patients with SLE and COVID-19 tended to have an increased proportion of LCs (Fig. 3, B and C), with an enrichment in conventional CD4 + or CD8 + T cells and NK cells and, to a lesser extent, in γδ T cells. We also observed dp CD45 + CD4 + CD8 + CD3 + TCRγδ − and dn CD45 + CD4 − CD8 − CD3 + TCRγδ − T cells in all inflamed and some healthy-looking samples, albeit in smaller numbers (Fig. 3, B and C, and fig. S4). Such populations of T cells were often understudied, as CD4 and CD8 biomarkers are thought to be mutually exclusive; however, they have been often reported in autoimmune and chronic inflammatory disorders (19,34), including SLE (35,36).
We next performed an unsupervised clustering of all patients and healthy-looking controls based on the quantitative analysis of their immune signature using both a detailed heatmap based on immune profiles (Fig. 3C) and a principal components analysis (PCA) per patient (Fig. 3D). Healthy-looking skin samples clustered together, with no apparent relationship with the pathological samples (Fig. 3, C and D). Patients with Kawasaki syndrome and COVID-toes had a tendency to form disease-specific clusters, while patients with SLE were distributed between both conditions (Fig. 3, C and D). Although these data were obtained on a restricted number of patients, they suggest that all analyzed pathological lesions displayed common quantitative immune features (Fig. 3C), with nevertheless potential disease-intrinsic characteristics suggested upon analysis with a dimensional reduction PCA (Fig. 3D). To explore this hypothesis further, we refined our analysis by investigating the activation status of conventional CD4 + and CD8 + T cells based on their expression level of CD57, a biomarker classically associated with a high cytotoxic potential (i.e., pro-tissue damage) during viral infections and autoimmune disorders, including COVID-19 (37). We found that, compared to other pathological conditions, three COVID-toes cases were particularly enriched in CD4 + and CD8 + T cells, exhibiting high levels of CD57 (i.e., CD57 high ; calculated as CD57 MFI z score; Fig. 3, E to G).
During inflammatory skin conditions, cytotoxic immune cells can relocate nearby to/in contact with keratinocytes and contribute to severe epidermal damage (38,39). To calculate the anatomical location of all immune cells with respect to the epidermal layers, we acquired the spatial coordinates of the modeled epidermis. We next incorporated into MANTIS a k-dimensional tree algorithm (40,41), which automatically decomposes the structural element coordinates (i.e., as exemplified here with the epidermis) into virtual subspaces and enables us to calculate the nearest neighbor  S6E) were all significantly enriched near the epidermal layer in cases of COVID-toes. CD8 + CD57 low T cells were not found enriched in the epidermis (Fig. 3, H to K), suggesting a biological link between the expression levels of CD57 and epidermal migration in CD8 T cells. Although the number of patients studied is limited, these findings strongly suggest the potential formation of tissue-damaging subepidermal inflammatory clusters composed of cytotoxic T cells and NK cells in COVID-toes.
MANTIS enables topographic exploration of skin lesions by solving the α-shape of in situ immune substructures Inflammatory dermatoses are characterized by the presence of large inflammatory infiltrates composed of specific immune cells and thought to be critical for the development of the pathology (e.g., type 2 immune cells and eosinophils in atopic dermatitis). To better understand the regional heterogeneity of pathological lesions from SLE, Kawasaki syndrome, and COVID-toes, we took advantage of α-shape algorithms that enable us, by tuning the α parameter, to define a precise shape of sets of points by drawing bounding polygons based on the principle of Delaunay triangulation (42). When combined with the digital immune landscapes generated with MANTIS, α-shape algorithms automatically generate polymorphic α-shapes around n-clusters composed of a minimum of 15 cells (i.e., 15 being the minimum number of cells often found in clusters of inflammatory but not in healthy-looking samples; Fig. 4A). This method enables us to automatically detect and quantify the major inflammatory clusters (i.e., named hereafter "αROIs") to provide a high-level view of the in situ immune architecture of the skin lesion for each patient and disease. We generated lymphoid (Fig. 4, B and C) and myeloid (Fig. 4, D and E) αROIs for all the samples. Healthy-looking controls displayed a few lymphoid αROIs, and four of five controls did not show myeloid αROIs. These data indicate that, in human skin at the steady state, LCs have a tendency to form aggregates [i.e., composed of perivascular T lymphocytes (1)], while myeloid cells are more likely to be randomly distributed. In line with the data presented in Fig. 3, we found a higher proportion of both lymphoid and myeloid αROIs in pathological samples as compared to healthy-looking controls (Fig. 4, B to E). Using global unsupervised hierarchical clustering of αROIs per disease, we can generate a high-level view of inflammatory cluster composition and observe trends in disease-specific immune responses (Fig. 4, F and G). Notably, lymphoid αROIs of both Kawasaki syndrome and COVID-toes exhibited a significantly higher proportion of CD4 + CD57 low T cells than that of patients with SLE (Fig. 4, F and H). Both COVID-toes and SLE lesions displayed significant clusters of CD8 + CD57 high cytotoxic T cells, highlighting the cytolytic aspect of pathological lesion microenvironment in these conditions (Fig. 4, F and I). We next analyzed myeloid αROIs for all cases. We found that, COVID-toes had

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E
a particularly high density of clusters enriched in activated HLA-DR high dDCs (Fig. 4, G and J). Conversely, Kawasaki syndrome lesions showed an enrichment in both HLA-DR high LCs (Fig. 4, G and K) and mast cells (Fig. 4, G and L), while SLE lesions displayed large aggregates of eosinophils (Fig. 4, G and M). This finding is consistent with previous reports of strong eosinophilia in SLE (43)(44)(45). No significant differences were observed regarding other immune cell subsets in αROIs. While the precise role played by specific inflammatory clusters of immune cells in each disease remains elusive, these data strongly suggest that combining MANTIS digital maps with α-shape-based algorithms can reveal a significant nonrandom distribution of skin immune cells in skin lesions, with the presence of disease-specific immune structures. MANTIS analytic pipeline can thus enable us to quickly solve the spatial organization of complex immune environments and open interesting perspectives for future investigations in the field of dermatopathology.

DISCUSSION
Here, we propose a general framework for 3D quantitative and spatial analysis of skin immune cells at the cellular level. We first describe a simple method to perform a fast 3D acquisition of up to 10 biomarkers simultaneously and extract a single-cell database containing the biological identity (including spatial coordinates) of skin lymphoid and myeloid cells. We then analyze the extracted databases using an automated and interactive analytic pipeline composed of phenotype attribution matrices coupled with cell-tostructure distance calculations and α-shape algorithm-based detection of major inflammatory clusters. Our analysis was focused on FFPE samples, as it is still the most easily available source of pathological tissues and can enable analysis of patients' skin-sampled in routine clinical practice. However, the use of cryopreserved samples is compatible with the approach that we describe here and enables the analysis of thicker tissue sections ( fig. S2).
We identified that the first step of the process, which consists of the generation of good-quality 3D multiplexed images with a great ratio signal over noise, is critical for the rest of the study. This is why we emphasized the capability of a conventional 8-year-old (noncustom-built) confocal laser scanning system to acquire 10 different fluorochromes simultaneously. This method of acquisition can be democratized to most academic/clinical facilities because it involves a conventional equipment coupled to basic spectral spillover compensation and single-cell data extraction strategies, via the use of commercially available software (Huygens and Imaris; described in detail in Materials and Methods).
Single-cell segmentation is also very critical, as it will constitute the very core of the future analysis of immune subpopulations and expression of biomarkers. Possible mistakes made at this step, e.g. the inability to separate immune cells in large infiltrates, would then result in misinterpretation of MANTIS-generated results. We tested different approaches to automatically segment healthy-looking and inflammatory skin samples, including random forest-based classifiers (e.g., ilastik machine learning). While such a method was suitable to segment healthy-looking images with a low concentration of immune cells, it failed to segment complex inflammatory lesions, where large and packed immune clusters were present. We thus opted for a semisupervised segmentation of single immune cells using the software Imaris, in which the segmentation of each inflammatory cluster was quality-controlled manually and was in 3D (Fig. 1F). While this approach is probably more time-consuming, we could ensure a precise 3D segmentation and further extraction of an accurate single-cell database to be processed with MANTIS. A recent study has reported the use of an analysis pipeline, including a new segmentation strategy, adapted from the field of astronomy named "AstroPath" (46). Using this approach and only six biomarkers, they could identify important pathological features in biopsies from patients with melanoma. These results, in line with our findings, highlight the importance of carefully selecting a list of biomarkers to be studied and of having the right analytic pipeline to draw reliable insights.
Because immunologists are more commonly used to identifying immune cell populations with manual gating of populations based on flow cytometry, we validated the quantitative performance of MANTIS by analyzing the extracted single-cell databases for the 15 patients analyzed with the conventional flow cytometry software, FlowJo. We found that the MANTIS-based analysis on 3D images generated with two different panels and just a minimal number of 10 antibodies per panel were sufficient to distinguish 19 immune subsets and identify disease-specific trends in skin lesions.
Because MANTIS attribution matrices can be quickly adjusted to any set of markers, they could be compatible with single-cell databases generated using technology with high multiplexing capabilities such as CODEX (7,8), MIBI-TOF (9), IMC (10), MuSIC (11), CyCIF (12), Cell Dive (13), and others (14). We included in MANTIS the α-shape algorithm that enables us to define the precise shape of the inflammatory immune clusters based on the principle of Delaunay triangulation (42). When applied to digital immune landscapes, the α-shape algorithm automatically identifies and quantifies dermal and epidermal inflammatory clusters (i.e., αROIs) composed of a minimum of 15 immune cells (Fig. 4A). This method enables one to directly analyze the major αROIs to provide a fast high-level view of the skin immune architecture in a given lesion. Using this method, we could quickly identify that the immune subset (e.g., mast cells, HLA-DR high dDCs, CD8 + CD57 high cytotoxic T cells, eosinophils, etc.) was specifically enriched in dermal areas only in some patient subgroups. These preliminary observations suggest that, depending on their etiology, pathological lesions could be due to distinct pathological mechanisms. This type of analysis opens interesting perspectives for the 3D cartography of complex inflammatory skin lesions and should be pursued by additional studies on a larger number of patients. While 2D immune landscapes are represented here to facilitate the visual assessment in the figures (3D graphs are hardly perceivable on static pictures), the single-cell segmentation and extraction of cellular spatial coordinates were all performed in 3D.
On the basis of CODEX high multiplexing capacity, previous studies (8,47) have shown that it is possible to generate a highlevel view of the cell-to-cell interaction landscape based on the principle of the Delaunay neighborhood graph (48). The MANTIS αshape algorithm is complementary, as it automatically identifies major immune structures while deciphering their cellular composition. Combining α-shape and neighborhood approaches could help quickly solve the biology of major inflammatory clusters in the skin, by drawing the ligand-receptor interactome of immune and structural cells within the identified cluster. Such a high-dimensional analysis of the skin immune architecture could provide a promising avenue for understanding the complexity of inflammatory skin manifestations with potential benefits for patient stratification and/ or diagnosis.
There is a strong need to design new tools to assist clinical decision-making and/or better apprehend the complexity of inflammatory dermatoses. While very promising processes have been made in the field of spatial biology (49)(50)(51), there is an unmet need for a nonexpensive and standardized multiplexed imaging analytical framework capable of automatically resolving the immune architecture of an inflamed skin. Here, we show that the MANTIS analytical system is uniquely positioned to examine numerous questions in the fields of skin immunobiology and should lay the foundation for a fast and automated analysis pipeline of relevant in situ inflammatory environments in both research and clinical facilities.

Human skin samples
The pathological biopsies were performed as part of routine care for diagnosis purposes in Lyon, Reims, and Toulouse university hospitals. All patients have given written informed consent to the publication of their case details. Skin biopsies from patients with lupus erythematosus were obtained from Toulouse University Hospital. Biopsies of COVID-toes were obtained from Toulouse, Reims, and Lyon university hospitals. Skin biopsies of multisystem inflammatory syndrome were obtained from Reims University Hospital. Anonymized healthy human skin samples were obtained from donors that underwent abdominoplasty procedures and had given their written informed consent. Donors did not have any record of allergies or dermatological disorders and did not use corticosteroids. Control healthy skin biopsies were age-and gendermatched with pathological samples and obtained from Genoskin SAS (www.genoskin.com/). Genoskin has obtained all legal authorizations necessary from the French Ministry of Higher Education, Research and Innovation (AC-2017-2897) and the Personal Protection Committee (2017-A01041-52). All studies were conducted according to Declaration of Helsinki protocols.

Skin section preparation, histology, and staining
Human skin samples were either frozen in optimal cutting temperature compounds (OCT, Tissue-Tek) or formalin-fixed and paraffin-embedded. FFPE-tissue sections (10 μm) were heated at 95°C for 20 min. Sections were subsequently immersed into xylene for 30 min, washed in a graded series of ethanol (100, 95, 70, 50, and 30% for 5 min each), and abundantly washed with distilled water. They were then treated using a heat-induced epitope retrieval method as previously described (52).
FFPE-tissue sections were blocked and permeabilized with phosphate-buffered saline (PBS), 0.5% (w/v)% bovine serum albumin (BSA; Sigma-Aldrich), and 0.3% Triton X-100 (Merck) for 30 to 60 min at room temperature and then incubated with fluorophore-coupled antibodies or unconjugated antibodies overnight at 4°C in the dark. The sections were then washed three times in PBS 0.5% (w/v)% BSA, 0.3% Triton X-100 and incubated, if needed, with secondary antibodies in PBS, 0.5% (w/v)% BSA, and 0.3% Triton X-100 for 2 hours at room temperature in the dark. Last, samples were treated with an autofluorescence quenching solution named TrueVIEW (Vector Laboratories) for 5 min. The slides were mounted in Mowiol medium (Sigma-Aldrich) and sealed with a coverslip. All conjugated and unconjugated antibodies used in this study were validated in single immunostainings of human skin and tonsils and are listed in table S1.
Acquisition Z-stack images (512 × 512 pixel; 1 μm) were acquired using an 8year-old confocal microscope SP8 (Leica Microsystems) equipped with an HC PL APO CS2 with 40× numerical aperture 1.3 oil objective, an ultraviolet diode (405 nm), and four lasers in visible range wavelengths (405, 488, 532, 552, and 635 nm). The setup was made up of five detectors [three hybrid detectors with high quantum yield compared to classical photomultiplier (PMT) detectors and two PMTs]. Mosaic sequential images were acquired using the between-stack configuration to simultaneously collect individual seven or eight channels and tiles before merging them to obtain one single image. Use of the between-stack configuration and the modulation of the detectors' detection windows help reduce the leaking of fluorophores. Last, a digital zoom of 1.9 was applied during the acquisition and a mosaic multicolor image was obtained and exported into a .lif format. Detection windows and microscope configuration used in our study are listed in table S2.
Image deconvolution and correction of spectral spillover 3D mosaic images were then imported into Huygens SVI software to correct the signal by applying deconvolution and cross-talk correction. Two deconvolution methods were used: the express deconvolution (theoretical and fast) or the deconvolution wizard (possibility to use experimental or theoretical parameters and to adjust the background value). Automatic cross-talk correction estimation was obtained, and the coefficients were slightly adjusted manually, if needed, for optimal spillover correction.

Segmentation
3D mosaic images were imported into Imaris software to separate objects (cells) using a 3D surface segmentation. Before creating the surface objects in Imaris, classical image processing was required. For instance, defining a threshold, adding a median filter, and/or normalizing the layers were sometimes applied to clean the background. Images were either cleaned using the CD45 surface objects or other channels by applying appropriate masks for each channel. Then, segmentation was applied on the CD45 channel surface. Statistics were exported into .csv format.

Epidermis modeling
The epidermis was identified using the natural autofluorescence of the tissue. On the basis of the autofluorescence found in multiple appropriate channels, an epidermis surface was created using the Object Creation semi-automated tool of the Imaris software. The coordinates of the epidermis were then exported into .csv format. The modeling of the epidermis shape showed on the digital maps was obtained using the α-shape algorithm (described below).

Segmentation troubleshooting
In some cases, the surface creation parameters were not efficient in automatically obtaining good object creation, or the module was not sensitive enough to detect low-intensity objects. In this case, the creation of small objects was done manually, and the threshold selection was also reduced. If the detected object was below 1 μm, then a manual object unification with surrounding objects of the same intensity was performed.

Statistical data exportation
Statistical properties of each segmented object (cell) in the processed 3D Imaris Multiplex image were automatically calculated. Object volume, sphericity, area, xyz position, and MFI in all channels were exported as a .csv table.

Skin enzymatic digestion and gradient separation
Immune cells were isolated from fresh skin biopsies from Genoskin. Briefly, skin was harvested in predigestion medium, fragmented into pieces, and incubated at 37°C on a rotating plate to remove epithelial cells and other impurities. The supernatant was discarded, and the samples were digested for 45 min on a rotating plate with 1.25 mg of Liberase (Sigma Aldrich, #5401127001) and 2.5 mg of deoxyribonuclease I (Sigma Aldrich, #10104159001) to disaggregate the tissue. Samples were further dissociated with the Miltenyi gentleMACS Dissociator. Cells were then enriched with a Percoll gradient, washed, and used for staining with the same panels used for MANTIS. Data were acquired on a FACSymphony (BD) cytometer and were analyzed using FlowJo (Tree Star) software.

Implementation of MANTIS reference panels
To enable cell identification, we built a binary table containing a literature-based theoretical signature of biomarkers expressed in each cell population identified by the used panel (naturally depending on the used set of antibodies), known as the reference attribution panel. If a cell population is positive for a marker, then the value is set to 1; otherwise, it is set to 0. If a cell population can be positive for a marker, then there are two columns, one with the value set to 1 and the other with the value set to 0 (i.e., γδ T cells can express CD4 or not). Two reference tables were implemented and designated by lymphoid and myeloid reference attribution matrices.

Dynamic adaptation of reference matrices
Sample heterogeneity led to different acquisition parameters (laser power, gain, etc.). To standardize data processing, we scaled the reference tables and dynamically adapted, for each sample, the table values according to the MFI values. In practice, the value "1" in the binary table was replaced by the maximum MFI value acquired in the corresponding channel from the tested sample.

Automatic cell type identification
To annotate the segmented objects, a correlation matrix between the MFI table and the adapted reference panel was generated by performing a pairwise Spearman's rank correlation using the R software (2021). Each object was then phenotypically assigned to the cell type having the highest correlation coefficient. Objects with multiple highest correlation coefficients were assigned as "Other" cell types.

Accuracy validation
The accuracy of MANTIS automatic cell identification was verified by comparing quantification results to classical histo-cytometry (6). Briefly, linear regression of cell type density was computed between both attribution methods, and regression coefficients were calculated. Regression coefficients ranging between 0.75 and 1 reflect MANTIS technique robustness.
Activation status detection MANTIS panels were designed to include not only discriminant markers for cell attribution but also nondiscriminant and informative markers, for instance, activation markers. The cell populations of interest (CD4 + and CD8 + T cells in the lymphoid panel and DCs, LCs, and CD207 + DCs in the myeloid panel) and the activation markers that reflect the activation status of these populations (CD57 in the lymphoid panel and HLA-DR in the myeloid panel) were defined in the MANTIS algorithm. The latter automatically computes the MFI density curve associated with the activation markers within the selected populations. Subsequently, the MFI corresponding to the first peak of the density curve is defined as the MFI value above which the cell is considered positive for the activation marker.
α-Shape calculation α-Shape was calculated using the alphashape Python package. Briefly, Delaunay triangulation of a given set of points formed a bounding polygon that contains all the points of the set. The α parameter was defined by the value α, and a circle with 1/α radius was drawn in such a way that two points of the dataset are located on the boundaries of the circle and the circle is empty. For each empty circle found, the line between the two points formed a side of the bounding polygon, i.e., the α-shape. As α decreased, the α-shape changed from a convex hull (e.g., epidermis α-shape, α = 0.4) to a more tightly fitting bounding box, resulting in more refined αshapes [e.g., ROI α-shape (αROI), α = 0.1].

Cell-to-structure distance calculation and nearestneighbor search
x-y coordinates of epidermis α-shape contours were stored using the k-dimensional tree method, which allows data ranking and structuration. Briefly, data points were classified on the basis of nodes and branches space-partitioning, allowing a fast nearestneighbor calculation. For a given point (cell) of the dataset, the nearest neighbor in the epidermis α-shape was found, and the distance defined by r was calculated using the scipy.spatial Python package (40). The distance of cells contained in the epidermis αshape was set to 0. Data clustering and αROI analysis ROIs (αROI, i.e., inflammatory cell clusters) were identified using the α-shape algorithm with a tuned α parameter (α = 0.1), allowing correct detection of high-cell density areas. αROI with less than 15 cells were removed from the analysis. For each selected αROI, specific characteristics were calculated and extracted, such as area, total number of cells, cell number, proportion by cell type, and αROI center coordinates.

Data visualization
Visualization charts were obtained using the ggplot2, Pigengene, and ComplexHeatmap R packages and matplotlib and seaborn Python packages. t-SNE was computed with Rtsne.

Statistics
Statistical tests were performed using Prism 8 (GraphPad Software) and the Rstats and rstatix R packages. One-way analysis of variance (ANOVA) with Tukey's test for multiple comparisons or Mann-Whitney test was performed on samples as noted in the respective figure legends. A P value of less than 0.05 was considered statistically significant.

Supplementary Materials
This PDF file includes: Figs. S1 to S6 Tables S1 to S4 Legend for movie S1 Other Supplementary Material for this manuscript includes the following: Movie S1 View/request a protocol for this paper from Bio-protocol.