Spatial Point Pattern Analysis Identifies Mechanisms Shaping the Skin Parasite Landscape in Leishmania donovani Infection

Increasing evidence suggests that in hosts infected with parasites of the Leishmania donovani complex, transmission of infection to the sand fly vector is linked to parasite repositories in the host skin. However, a detailed understanding of the dispersal (the mechanism of spread) and dispersion (the observed state of spread) of these obligatory-intracellular parasites and their host phagocytes in the skin is lacking. Using endogenously fluorescent parasites as a proxy, we apply image analysis combined with spatial point pattern models borrowed from ecology to characterize dispersion of parasitized myeloid cells (including Man^R+ and CD11c⁺ cells) and predict dispersal mechanisms in a previously described immunodeficient model of L. donovani infection. Our results suggest that after initial seeding of infection in the skin, heavily parasite-infected myeloid cells are found in patches that resemble innate granulomas. Spread of parasites from these initial patches subsequently occurs through infection of recruited myeloid cells, ultimately leading to self-propagating networks of patch clusters. This combination of imaging and ecological pattern analysis to identify mechanisms driving the skin parasite landscape offers new perspectives on myeloid cell behavior following parasitism by L. donovani and may also be applicable to elucidating the behavior of other intracellular tissue-resident pathogens and their host cells.

Introduction

Vector-borne kinetoplastid parasites of the Leishmania donovani complex are the causative agents of zoonotic and anthroponotic visceral leishmaniasis (VL). Leishmania parasites occur in two distinct stages: as extracellular promastigotes in their insect vector, hematophagous phlebotomine sand flies (Diptera: Psychodidae: Phlebotominae) (1), and as intracellular amastigotes in the mammalian host (2). VL is associated with 3.3 million disability-adjusted life years (DALYs) (3) and ~20,000 annual deaths in 56 affected countries (4). There is no available prophylactic vaccine (5) and available drugs are either expensive and/or show significant toxicity (6).

Infection of a naïve host with Leishmania parasites occurs through the bite of infected sand flies into the skin, the only interface between insect vector and host. L. donovani disseminate from the skin to deeper tissues including the spleen, liver, and bone marrow causing systemic pathology, but the skin is commonly asymptomatic (7). Blood parasite load (parasitemia) is considered the main determinant and a surrogate measure of outward transmission potential (i.e., infectiousness) in VL (8). Other studies, however, point to the importance of the skin as an additional reservoir of parasites. In a recent study (9), the infectiousness of L. donovani-infected patients in India, measured by xenodiagnoses, was found to positively correlate with spleen parasite load as well as blood parasitemia determined by qPCR. Of note, out of 7 infectious patients that remained infectious after drug treatment, 6 had no detectable blood parasitemia (9). Furthermore, microbiopsy sampling of skin in a VL-endemic region in Ethiopia indicated that parasites can be detected even in asymptomatic individuals, though multiple biopsies from the same patient were often qualitatively and quantitatively different in terms of parasite load (10). We previously demonstrated in a RAG mouse model of VL that despite high parasite burden in deep tissues, parasitemia was comparatively low (2,000–2,500 parasites/ml blood or equivalent to ~2–3 amastigotes/blood meal) and did not correlate with outward transmission success (11). Conversely, there was a good correlation between outward transmission and mean parasite burden in skin (11). Collectively, these data, together with studies in the hamster model of VL (12), in canine VL (13) and in human post-kala azar dermal leishmaniasis (9, 14, 15) support the need to consider skin parasite load as a potential contributor to infectiousness across the disease spectrum attributable to the L. donovani complex.

A confounding feature associated with sampling in many studies relating host infectiousness to skin parasite residence is the possibility of spatial heterogeneity in parasite distribution. Using an immunodeficient model of L. donovani infection in mice, we recently demonstrated that obligatory intracellular amastigotes accumulate in patches in asymptomatic host skin, and that this patchy landscape is a strong predictor for the infectiousness of an asymptomatic host (11). However, a detailed understanding of the mechanisms involved in the micro- and macro-scale dispersion of Leishmania parasites (and, de facto, of parasitized host cells) in the skin is lacking.

Here, we describe in detail the micro- and macro-scale dispersion of tandem-dimer Tomato (tdTom) L. donovani parasitized cells in the skin of 23 long-term infected RAG mice. The absence of acquired immune pressure in RAG mice allows underlying innate pathways controlling the dispersion of parasitized myeloid cells to be more clearly studied. Due to the intracellular nature of Leishmania amastigotes, their location is equivalent with that of their host cells. Thus, detection of amastigotes serves as a proxy for detection of parasitized myeloid cells in vivo. We used image analysis tools to extract data based on tdTomato signal from whole skin images (ImageJ) and also from microscopy images of punch-biopsy skin sections (StrataQuest; see Methods). We then adopted spatial point pattern methodologies from ecology and other fields to interrogate signal dispersion in silico and to make predictions about modes of dispersal (16). Our results suggest that patches representing infected myeloid cells form “clusters” with a larger patch at the cluster center surrounded by smaller patches. These clusters are consistent with reseeding L. donovani amastigotes from a larger patch within a locally varying radius around it. Our data therefore suggest that in the absence of acquired immunity, these patches are the main driver of their own dispersal. This process would also serve to enhance patch growth, as larger patches consist of a merger of smaller patches. By translating established mathematical models into an immunological context, this study provides new insights into the complex mechanisms underlying the dispersal and dispersion of L. donovani-infected myeloid cells in the skin and exemplifies an approach that may help understand dispersal and dispersion of phylogenetically diverse intracellular pathogens in the skin and other tissues.

Materials and Methods

Ethics Statement

All animal usage was approved by the University of York Animal Welfare and Ethics Review Committee and the Ethics Review Committee at FERA Science Ltd., and performed under UK Home Office license (“Immunity and Immunopathology of Leishmaniasis” Ref # PPL 60/4377).

Mouse, Leishmania, and Sand Fly Lines

C57BL/6 CD45.1.Rag2^-/- (RAG) mice (originally obtained from the National Institute of Medical Research, London, UK) were used. All animals were bred and maintained at the University of York according to UK Home Office guidelines. The Ethiopian Leishmania (Leishmania) donovani strain (MHOM/ET/1967/HU3) had been transfected with a tdTom gene (17) to generate a tdTom-L. donovani line in another study (18). The parasites were maintained in RAG mice by repeat-passage of amastigotes. For both parasite passage and RAG mouse infection, tdTom-L. donovani amastigotes were isolated from the spleens of long-term infected RAG mice and i.v. injected at 8 × 10⁷ amastigotes/uninfected RAG mouse as previously described (11).

Mouse Skin Harvest

All animal data used in this study were repurposed from 23 procedural control animals pooled from 7 different experiments. All animals were infected in the same manner with tdTom-L. donovani (see above), while the applied experimental conditions were comparable and neither adversely affected parasite and myeloid cell behavior nor contributed a significant source of variability between animals (Figures S1A, B). Mouse skin harvest was performed as previously described (11). Briefly, RAG mice were sacrificed in a CO₂ chamber, followed by cervical dislocation. Animals were first shaven with electric clippers and then treated with a depilation cream (Veet by Reckitt, Slough, UK) for up to 3 min. The depilation cream was scraped off and residual cream was washed off under running water. Mice were skinned by an abdominal double-Y incision, removal of ears, and exclusion of paws and tail (Figure S1C). The skin was pulled from the rump over the head and then cut off at eye level. Excess adipose tissue was removed with curved tweezers and skins were kept in complete RPMI on ice until stereomicroscopic imaging and subsequent biopsy punching.

Genomic DNA Extraction

Sixteen or 24 skin biopsies (ø 0.6 cm) were taken from each skin as previously described (11). All genomic DNA (gDNA) extractions were performed with either the DNeasy^® Blood & Tissue spinning column kit (Qiagen, Venlo, Netherlands) or the DNeasy^® 96 Blood & Tissue plate kit (Qiagen) according to supplier’s protocol.

Quantitative Polymerase Chain Reaction

All quantitative real-time PCRs (qPCR) were performed as previously described (11). Briefly, previously characterized Leishmania-specific kinetoplastid DNA primers (Accession number AF103738) (19) were used at a final concentration of 200 nM. Two nanograms of total gDNA extracted from skin biopsies were used/reaction. Fast SYBR^® Green Master Mix (Applied Biosystems) was used according to supplier’s guidelines. Reactions were run in a StepOnePlus™ thermal cycler (Applied Biosystems, Waltham, MA, USA) with a thermal cycle of 95°C for 20 s, a cycling stage of 40 cycles of 95°C for 3 s, 60°C for 5 s, 72°C for 20 s, and 76.5°C for 10 s (data read at final step), followed by the standard melt curve stage. Data were analyzed by StepOne™ Software v.2.3.

Microscopy

Fluorescent stereomicroscopy of whole skins was performed as previously described (11). Briefly, a series of 30–40 images at 12× magnification were taken/skin to image the whole skin area from the hypodermal face. All images were taken under the same conditions (exposure time: 400 ms). Image series were collated in Adobe^® Photoshop^® to render a single image of the whole mouse skin (Abode Inc., San Jose, CA, USA).

For confocal microscopy, fixation and cryo-preservation of skin punch biopsies (ø 0.6 cm) was adapted from Accart et al. (20) to optimize tdTomato signal and skin integrity preservation after skin freezing and cutting (20). Skin biopsies were fixed in 2% formaldehyde in PBS for 15 min at room temperature (RT), rinsed 3× in 15% Sucrose/PBS and incubated overnight at 4°C in a Tris-Zinc fixative (per 100 ml: 1,211.4 mg Tris-base, 501.5 mg ZnCl₂, 500.5 mg zinc acetate dihydrate, 50.6 mg CaCL₂, 15 g sucrose dissolved in dH₂O) on a shaker. The next day, biopsies were rinsed 3x in 15% sucrose/PBS and incubated in 7.5% porcine gelatine (300 bloom) in 15% sucrose/PBS in individual cryo-molds at 37°C for 1 h. Biopsies were sealed upright in the gelatine blocks by congealing at RT. Gelatine/specimen blocks were cut down in size, flash-frozen in isopentane on dry ice, further frozen in OCT in a cryo-mold on dry ice and stored at −80°C. Frozen skin biopsies were longitudinally cut at 20 µm thickness to expose all skin layers and section were placed on super-frosted microscopy glass slides. Gelatine was removed by slide incubation in PBS at 37°C for up to 30 min.

Antibody labeling of skin sections was adapted from Dalton et al. (21). Briefly, slides were rinsed in PBS, permeabilized for 15 min with 0.3% Triton X-100/PBS (perm. buffer), rinsed in PBS, blocked for 30 min with Image-iT™ FX Signal-enhancer (Invitrogen, Waltham, MA, USA), rinsed with PBS, and further blocked for 30 min with 5% serum + FcBlock (1:1000) in perm. buffer. Primary antibodies were diluted in 5% BSA in perm. buffer [CK-10 (1 µg/ml) (Abcam: ab76318)], Mannose Receptor [1 µg/ml] (Abcam: ab64693), CD11c [5 µg/ml] (Abcam: ab33483)] and incubated with specimen for 90–120 min at RT (Abcam PLC, Cambridge, UK). Slides were washed 3 times in 0.5% BSA/PBS. Secondary antibodies were diluted in 5% BSA in perm. buffer [BV421 (1 µg/ml) (BioLegend: 406410); Alexa Fluor 647 (Invitrogen: A-21451)] and incubated for at least 60 min at RT (BioLegend, San Diego, CA, USA). Slides were washed 2 times in 0.5% BSA/PBS and once in PBS at RT. Specimens were counterstained with 0.2 µM YOYO-1 in PBS for 30 min at RT and washed 2 times in PBS. Slides were mounted in ProLong Gold (ThermoFisher Scientific, Waltham, MA, USA) and left overnight at 4°C. Slides were sealed the next day and imaged with a 40× oil-immersion objective (400× magnification) on a fully motorized invert LSM 880 confocal microscope with Airyscan (Zeiss, Jena, Germany) by z-stacked tile-scanning. To overcome high levels of background fluorescence in the skin, we exploited the lambda mode of the Zen black software package (v. 2.3) for linear unmixing of previously recorded fluorescent spectra.

Stereomicroscopic Image Analysis (Macro-Scale)

All whole-skin images were analyzed with custom-build IJ1 macros in Fiji ImageJ (22). Whole RAG mouse skins measured 40.4 cm² (RAG15) to 82.5 cm² (RAG10) in surface area (Table S1). In comparison, C57BL/6 mouse skins are negligibly thin, measuring on average 0.034 cm (female) and 0.037 cm (male) in thickness (23). Thus, we dismiss tissue thickness at the macro-scale and treated whole-skin images as two-dimensional surfaces to simplify the analysis. Patch detection was based on 8-bit grayscale analysis of the L. donovani tdTomato signal in the isolated red channel. The bright fluorescent signal of densely accumulated amastigotes per host cell strongly enhances the fluorescent-signal far beyond what endogenously fluorescent myeloid cells could have achieved. The obligatory intracellular nature of Leishmania amastigotes allowed us to use the parasite signal as a proxy for infected myeloid cells, which allowed their indirect detection.

To determine a reliable detection threshold, we applied a three-level approach of the “Multi-Otsu-Threshold” algorithm plugin for Fiji ImageJ to tdTom-L. donovani infected and naïve RAG skins (24). Naïve skins rendered maximum intensity thresholds <35, while infected skin showed a range of about 60–150 on an 8-bit grayscale. For best detection sensitivity versus specificity balance, we chose a threshold of 70 on an 8-bit grayscale. Furthermore, we had to define the minimum size of a patch. Best image resolution allowed for a pixel size of 21.47 μm × 21.47 μm = 461 μm², which did not allow to resolve for individual L. donovani amastigotes, which are 2–3 μm in diameter (r² * π = 1.5² μm² * 3.1416 ≈ 7.1 μm²). Within a mammalian host, L. donovani amastigotes are intracellular, most commonly residing in macrophages (25). Thus, L. donovani amastigote distribution can be described by their host cell distribution. The cell body shape of macrophages and DCs is plastic and, in tissues like the skin, is highly irregular, making it difficult to calculate the 2D surface area of such cells. Assuming that the overall cell surface is the same between a rounded and irregularly shaped macrophage, accepting an average diameter of 25 μm and assuming a flat cell body shape, we can calculate an approximate 2D cell area (r² * π = 12.5² μm² * 3.1416 ≈ 491 μm²) (26). A single macrophage is therefore defined by at least 4 pixels. We thereby defined the smallest arguable patch size as 2 × 2 macrophages (at least 9 pixels) in 2D. Thus, (3 * 21.47 μm)² = 4,148.65 μm² was set as the smallest possible size to define patches in our stereomicroscopic images. Using these signal detection and patch size thresholds on naïve RAG skins showed false-positive detection of ≤15 small patches (<0.1 mm²), which could largely be attributed to dust particles and other autofluorescent events (Figure S22A). Furthermore, manual patch counting in 1 cm² areas on infected skins was used to assess the false-negative detection rate, showing underestimation of very small size (<0.1 mm²; Figure S22B) and loss of detection sensitivity at larger patch fringes, underestimating patch areas. The macro for density-peak detection was adapted from a publicly available protocol on particle counting in cells (Andrew McNaughton, University of Otago, NZ).

Furthermore, two standard plugins available for Fiji ImageJ were used; a plugin published by Haeri & Haeri (27) for nearest neighbor distance (NND) analysis (27) and a plugin published by Ben Tupper to analyze patch clustering (28). For the latter, a maximum distance of 3.9 mm between patch borders was chosen. This value was based on the idea that an average sand fly measures ~2 mm from proboscis tip to abdominal posterior. If pinned down at its posterior and swung around, the proboscis tip would describe a circle representing the most distal point where a sand fly could bite irrespective of its orientation. Using 1.95 mm as a radius, this circular landing zone has a diameter of 3.9 mm (Figure S1D). Conceptually, any patch reaching into the landing zone may be accessible to a sand fly without it moving away from its landing spot. Even parasite patches bordering the periphery of the landing zone, might allow parasites from these patches to be incorporated in the blood pool, which has an average diameter of ~1 mm (11).

Confocal Image Analysis (Micro-Scale)

Parasite Distance from the Epidermis Measured by Volocity

Volocity (V.6; Quorum Technologies, Laughton, UK) was used for amastigote to epidermis distance measurements. Individual amastigotes were detected by eroding the tdTomato signal until signal segments coincided with separate sets of YOYO-1-labeled small nuclei and kinetoplasts. The epidermis was demarcated by consolidating all Brilliant Violet 421 signals for Cytokeratin-10 (CK-10) positive areas into a single mass under exclusion of skin appendages. CK-10 is a marker of the spinous layer located right above the basal layer. Thus, the epidermis is about twice as thick as indicated by CK-10. Distances were measured by taking the shortest distance from the centroid of the tdTom amastigotes to the edge of the CK-10 positive band. To collate distance measurements from images from the same RAG mouse and to compare between RAG mice, the measurements required normalization due to variations in skin thickness. We estimated the mean thickness of the epidermis, dermis, and hypodermis, respectively, for each image by averaging at least three cross measurements in different places for each layer and then adjusted all measurements by the factor of difference between the mean skin layer thickness and the optical skin thickness as measured by Sabino et al. (23). While a difference in sex did not significantly affect the macroscale analysis, at the micro-scale, skin shows significant difference between sexes. While epidermal thickness is comparable between male and female C57BL/6 mice, the dermis of male C57BL/6 mice is ~1.4 times thicker than that of female mice. Conversely, the hypodermis in female C57BL/6 mice is ~1.84 times thicker than that of male mice. Thus, male and female B6.RAG2^-/- mice are not directly comparable.

Measurements of Parasite Localization in the Skin Measured by StrataQuest

StrataQuest (TissueGnostics, Vienna, Austria) was used to measure L. donovani amastigote co-localization in and distribution of Man^R+ and CD11c⁺ cells in the skin. Amastigotes were detected by segmentation on the tdTomato signal and host nuclei were detected by segmentation on a mask of YOYO-1 signal minus tdTomato signal. Host nuclei were used as a seed and a mask expanded based on Brilliant Violet 421 signal of Man^R and Alexa Fluor 647 signal of CD11c to determine cell type. The epidermis mask was demarcated by consolidating the signal of Brilliant Violet 421 signal of Man^R, Dylight 650 CD11c and YOYO-1 into a grayscale image and segmented using StrataQuest algorithms to distinguish from the dermis. Tissue morphology was used to define masks for the dermis and hypodermis. The tissue and cell masks were used to determine the location of the amastigotes. Host cells were classed as infected when at least one amastigote was found within them. However, for technical reasons, we did not differentiate between infected host cells by number of infecting amastigotes.

Code Availability

All Fiji ImageJ macros and R scripts used in this study were deposited in on GitHub (https://github.com/joedoehl/LeishSkinSPPM) and are freely available for download.

Statistics

All statistics were conducted in GraphPad PRISM v.8, SPSS v.25 (29) or Rstudio v.1.3.959 (30) running R v.4.0.2 (31). Prior to test application, parametric test assumptions (normality by Shapiro–Wilk test, equal variance by Brown–Forsythe test, etc.) were tested for all datasets. If these assumptions were not met even after data transformation, appropriate non-parametric tests were selected. In the case of the Quadrat test, which is based on chi-squared, chi-squared was replaced by a Fisher’s exact test whenever squares of the grid applied to the window of observations had several squares with <5 or at least one square with 0 observations, which was the case for all RAG mice. Prior to statistical analysis, all percent datasets were arcsin transformed. Applied statistical tests are stated throughout the text and in corresponding figure legends. Harvested mouse skins varied considerably in size between RAG mice (40.4–82.5 cm²). For comparison and correlation analyses between RAG mice, some data, like total patch area, patch counts, parasite counts were normalized to area/cm² and counts/cm², respectively. Where data skewness was an issue, we used the median as a more robust location measure. Spatial point pattern analysis was performed in Rstudio based on the work of Baddeley and co-workers (16). The following R packages were used for code execution: abind (32), deldir (33), e1701 (34), ggpubr (35), goftest (36), graphics (31), grDevices (31), polyclip (37), RandomFields (38), rlist (39), rms.gof (40), scales (41), sp (42), sparr (43), spatstat (16), spatstat.utils (44), stats (31), tensor (45), tidyverse (46), and utils (31). Each observation in the analysis was verified by at least two different tests. Running multiple test increases confidence in the result as different tests are not always in agreement.

Spatial Point Patterns Methodology

A “spatial point pattern” is a dataset containing the spatial locations of events or objects as points within a window of observation. The characterization of the spatial arrangement of these points is the main focus of spatial point pattern methodologies. This is generally done by measuring inter-point distances and characterizing the dispersion between points. Such methodologies find application in a myriad of fields, including epidemiology, revealing important spatial features not discernable by eye (47). Of course, not all “point data” permit analysis as a “point pattern”, particularly when sampling location is not relevant to the ultimately studied sample property. Conversely, datasets may represent point patterns that are not immediately recognizable, e.g., when the objects’ physical size exceeds technical classification as points; such limitations can be accommodated by choosing the coordinates of the center of mass of the objects as the point and attaching the objects’ size as an attribute, termed mark, to the respective points (48). This method was used in our study where the size of some observed skin patches clearly exceeds the definition of point objects. Very large skin parasite patches also showed great heterogeneity in tdTomato signal density, containing multiple high-density areas. This suggested that very large skin patches were conglomerates of smaller patches that had coalesced forming continuous, although heterogeneously densely parasitized, patches.

We used functions provided by the spatstat R package by Baddeley and co-workers (16). The following definitions are all taken from their work “Spatial Point Patterns: Methodology and Applications with R” (16). Several functions used in our study were derived from the empirical K-function, originally proposed by Ripley in 1977 (49):

$$\widehat{K}(r)=\frac{\left|W\right|}{n(n-1)}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{{}_{j\ne i}^{j=1}}^{n}1\{d_{ij}\le r\}{e}_{ij}}(r)}$$

where d_ij is the distance between all ordered paired points x_i and x_j (d_ij = ||x_i – x_j||), r is the value of measured distances (r ≥ 0), |W| is the area of the observation window, n(n – 1) is the total number pairs of distinct points, 1{…} is the “indicator” notation, which is 1 when statement “…” is true and 0 when false, and e_ij (r) is an edge correction weight. Riley’s K-function implicitly assumes that the investigated point process (X) has homogeneous intensity and is stationary. If a point process is suspected to be inhomogeneous, then that must be considered in the calculation of the K-function. The standard estimator of K can be extended to produce the inhomogeneous K-function:

$${\widehat{K}}_{inhom}(r)=\frac{1}{D^p\left|W\right|}{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}\frac{1\{d_{ij}\le r\}}{\widehat{\lambda }(x_i)\widehat{\lambda }(x_j)}}}e(x_i,x_j;r)$$

where e(x_i,x_j;r) is an edge correct weight and λ(x_i) and λ(x_j) are estimates of the intensity functions λ(x_i) and λ(x_j), respectively. The constant D^p is the pth power of

$$D=\frac{1}{\left|W\right|}{\displaystyle \sum _i\frac{1}{\widehat{\lambda }(X_i)}}$$

The inhomogeneous K-function assumes that the spatial scale of interaction remains constant, while the intensity is spatially varying. However, the spatial scale may be locally dependent. In this case, a rescaled template process may be applied that is locally stationary and isotropic, but shows a varying rescaling factor between locations within the observation window. Thus, the empirical K-function is adjusted for each pair of data points (X_i,X_j) by calculating the rescaled distance

$$d_{ij}^{\ast }=\frac{\left|x_i-x_j\right|}{s(x_i,x_j)}$$

where the rescaling factor is

$$s(x_i,x_j)=\frac{1}{2}\left(\frac{1}{\sqrt{\widehat{\lambda }(x_i)}}+\frac{1}{\sqrt{\widehat{\lambda }(x_j)}}\right)$$

This means the locally scaled K-function is defined by

$${\widehat{K}}_{scaled}(r)=\frac{\left|W\right|}{n(n-1)}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{{}_{j\ne i}^{j=1}}^{n}1\{d_{ij}^{\ast }\le r\}}}{e}_{ij}(r)$$

The L-function, as proposed by Besag (49), is then derived from the various K-function by

$$L(r)=\sqrt{\frac{K(r)}{\pi }};L_{inhom}(r)=\sqrt{\frac{K_{inhom}(r)}{\pi }};L_{scaled}(r)=\sqrt{\frac{K_{scaled}(r)}{\pi }}$$

respectively. The pair correlation function is derived from the K-function as

$$g_h(r)=\frac{K(r+h)-K(r)}{2\pi rh+\pi h^2}$$

When h becomes very small, then πh² becomes negligible and $\frac{K(r+h)-K(r)}{h}$ becomes the derivative of the K-function with respect to r (K’(r)). Thus, in two dimensions, the pair correlation function can be defined by

$$g(r)=\frac{K\prime (r)}{2\pi r}$$

More details can be found in Baddeley and co-workers book “Spatial Point Patterns: Methodology and Applications with R” (16).

Results

Parasites Accumulate in Phagocytic Cells in the Reticular Dermis

We generated confocal images of 20-µm-thick longitudinal skin biopsy sections from four different RAG mice (RAG19-22). Visual image inspection suggested that L. donovani amastigotes resided primarily in the dermis (Figure 1A). For confirmation, two distinct approaches were used. We calculated the shortest distance measurements of individual amastigotes to the epidermis with Volocity (Figure 1B) and conducted computationally aided amastigote counting in each skin layer with StrataQuest (Figure 1C) as two independent, but confirmatory approaches for assessing amastigote localization in the skin. The former showed that the majority of parasites resided in the reticular dermis, which comprises the deeper 90% of the dermis (50). The latter showed that the dermis harbored on average ~36× and ~9× more amastigotes/mm² than the epidermis and hypodermis, respectively.

FIGURE 1

Figure 1 Evaluation of micro-scale parasite distribution in the skin. (A) Representative tile-scan z-stacked confocal microscope image of a 20-µm-thick skin biopsy section of the RAG19 mouse. In cyan Cytokeratin-10 (CK-10; epithelium), in blue YOYO-1 (nuclei and kinetoplasts) and in red tdTomato (L. donovani amastigotes). The white lines demarcate the hypodermis’s upper and lower edges. Their locations were determined in raw images prior to linear unmixing when the autofluorescence of adipocytes was still visible. (B) Violin plot of normalized L. donovani amastigote distance from the epithelium as measured by the Volocity approach in confocal images. The number of amastigote distances measured are given above each plot. RAG19–21 were female and RAG22 was a male B6.RAG2^-/- mouse. The background color code for skin layer thickness: pink = epidermis, yellow = dermis, green = hypodermis, tan = subcutaneous tissue. The average sand fly proboscis skin penetration depth is marked by a dot-dashed line. (C–H) Evaluation of amastigote localization in the skin by the StrataQuest approach in confocal images. Each circle is representative for a measurement from one of 21 analyzed images across RAG19-22. (C) Parasite counts/skin layer (Friedman test: p < 0.0001; Dunn’s post hoc test: Epi vs. Der: p < 0.0001, Epi vs. Hypo: p = 0.0407, Der vs. Hypo: p = 0.0047). (D) Counts of infected host cell/skin layer (Friedman test: p < 0.0001; Dunn’s post hoc test: Epi vs. Der: p < 0.0001, Der vs. Hypo: p = 0.0021). (E) Correlation plot of total parasite counts/skin section against all infected cells/skin section (Spearman r: p < 0.0001, r = 0.8766, 95% CI: 0.7091–0.9505). (F) Representative image of the Man^R (cyan) and CD11c labelling (yellow) together with the endogenously expressed tdTomato signal from L. donovani amastigotes (red) and YOYO-1-labeled nuclei (gray). (G) Percent of detected infected Mannose Receptor positive cells of all infected cells compared to detected infected CD11c⁺ cells of all detected infected cells (Paired t-test: p < 0.0001, 95% CI: −33.11 to −18.02). (H) Percent of detected infected Mannose Receptor-positive cells of all detected Mannose Receptor cells compared to detected infected CD11c⁺ cells of all detected CD11⁺ cells (paired t-test: p < 0.0001, 95% CI: −11.82 to −4.74).

Leishmania amastigotes are obligatorily intracellular parasites, residing in phagocytic cells. Not surprisingly, therefore, infected host cells were also found predominantly in the dermis (~29× and ~26× more frequent than in epidermis and hypodermis, respectively; Figure 1D). Thus, parasite dispersion was reflected in host cell dispersion. In fact, the number of parasites/mm² skin correlated very strongly with the number of infected cells/mm² (Spearman r: p < 0.0001, r = 0.9364, 95% CI: 0.8432–0.9749; Figure 1E), suggesting that the proportion of infected cells increased with the increase in parasite burden.

In the skin, macrophages and dendritic cells (DCs) are the most abundant phagocytic cell populations in steady state (51). Alternatively, activated tissue-resident macrophages (TRM2s) expressing the mannose receptor (Man^R) are a preferred host for L. major in the mammalian skin (52). CD11c is a generic DC marker (53), although it can be expressed by tissue macrophages under inflammatory conditions (54). Thus, Man^R and CD11c were used to detect these myeloid cell subsets (Figure 1F). We observed that 42% (95% CI: 31%–52.7%) and 7% (95% CI: 4.7%–9.3%) of all infected cells were Man^R+ and CD11c⁺, respectively (Figure 1G). Thus, while Man^R+ TRM2s also appear to be permissive host cells for L. donovani amastigotes, ~51% of all infected cells detected were neither Man^R+ or CD11c⁺. Furthermore, of all Man^R+ and CD11c⁺ cells, 19.1% (95% CI: 12.6%–25.7%) and 9.8% (95% CI: 5.5%–14.1%) were infected, respectively (Figure 1H), leaving a large proportion of these cells uninfected (~80.9% and ~90.2% of all Man^R+ and CD11c⁺ cells, respectively). These data suggest that the observed diverse cellular tropism is not due to saturation of permissive targets such as TRM2s.

Patch Landscape Is Sculpted by Parasite Reseeding and Growth

Whole body skins from 23 untreated, tdTom-L. donovani long-term infected (6–7 months) RAG mice (Rag1–23) were used to analyze the dispersion of skin parasites within their phagocytic host cell. These RAG mice had been used as control mice in seven independent experiments (Exp 1–7) and were comparably treated (Figures S1A–C). We used qPCR-based parasite kDNA detection in multiple skin biopsies from each RAG mouse matched against a standard curve to estimate parasite burden per skin biopsy. Comparing parasite burden from biopsies across the skin of each mouse established whether skin parasite burdens varied geographically across the mouse skin and between mice and the mean parasite burden per biopsy per mouse (Figure 2A and Figure S2A). We found considerable variability in mean parasite burden between RAG mice (Kruskal–Wallis test: p < 0.0001) and in parasite loads between biopsies from the same RAG mouse (Figure S2A), confirming a heterogeneous parasite dispersion in the host skin (11). Furthermore, data on the parasite burden per biopsy showed in most cases a departure from a normal distribution for each RAG mouse (Shapiro–Wilks test: see Table S1) and a lack of equal variance (Brown–Forsythe test: p < 0.0001). Thus, the median parasite burden per biopsy was used to estimate the total parasite burden in skin for each RAG mouse. In turn, this was used to estimate the mean parasite burden per cm² of skin (Table S1). Despite intravenous inoculation of equivalent numbers of parasites, there was considerable variability in total parasite burden in skin between RAG mice, ranging from ~1 × 10⁶ (RAG12) to 1.24 × 10⁹ (RAG3) parasites per mouse (Table S1).

FIGURE 2

Figure 2 Observational description of the skin parasite landscape. All data presented here were generated by the ImageJ-based analysis of stereomicroscopic images. (A) Accumulative qPCR data of estimated parasite loads/biopsy from 432 biopsies (ø 0.6 cm) of 23 RAG mice. (B) Average patch counts/cm² of skin for 23 RAG mice. (I, K) Comparison of ratios of parasitized skin by body region for 23 RAG mice (Friedman Test: p = 0.3113 and p = 0.2931, respectively). (C) Accumulative data of the individual area of 21,904 measured patches from 23 RAG mice. (D) Representative image of detected patches in the RAG4 mouse skin. (E) Arranged data from 21,904 individually detected patches from 23 RAG mice arranged by “A”, the ratio of patch counts/size category, and by “B”, the contribution of size category to the overall patched area in the skin. (F) Correlation plot of total parasite patched skin by area (µm²) against the size of the biggest parasite patch detected in each of the 23 RAG mice (both on a log-scale; Spearman r: p < 0.0001, r = 0.9792, 95% CI = 0.9496–0.9915). (G) Schematic distribution of patches over the skin of RAG4 marked by patch area according to the categories in (E). (H) Correlation plot of 21,904 skin parasite patches by area (µm²) from 23 RAG mice (on a log-scale) against their respective mean pixel density on an 8-bit grayscale (Spearman r: p < 0.0001, r = 0.6554, 95% CI = 0.6476–0.6631). (I) Representative 8-bit grayscale image showing the tdTomato light intensity per pixel (RAG4). In magenta are the outlines of detected skin patches. (J) Accumulative data of mean patch density/patch of 21,904 individually detected patches from 23 RAG mice. (K) X–Y plots of five 5-mm cross sections through patches measuring tdTomato signal intensity/pixel on grayscale [0–255] in the RAG4 skin image. The arbitrary cutoff for patch detection set at 70 is marked in these graphs with black dashed line. (L) Image of the RAG4 image indicating where cross sections through patches were made.

Akin to parasite loads, skin patch counts were also highly variable between RAG mice, ranging from ~1.2 patches per cm² (Rag14) to ~32.9 patches per cm² (Rag8) and 76 (RAG14) to 2,457 (RAG10) patches per mouse (Figure 2B and Table S1). We also measured the area of skin patches. Patches measured ~0.00004 cm² (50× smaller than a 25G needle prick) to 53 cm² (>80% of the mouse skin; Figures 2C, D, Figure S3 and Table S1). Most patches were small (<0.1 mm²) and only few were large (>1 mm²), resulting in a heavily skewed dataset (Figure S2B). More precisely, ~84% of detected patches were 0.004 mm²–0.099 mm² but accounted for merely ~12.7% of the total patch area (Figure 2E). Conversely, only ~3.3% of detected patches were >1 mm² but accounted for ~64% of the total patch area. In an individual mouse, these proportions were even more extreme (Figure S4). In fact, there was a strong correlation between the total patch area and the largest patch in any given mouse (Spearman r: p < 0.0001, r = 0.9792, 95% CI = 0.9496–0.9915; Figure 2F). Interestingly, we identified no preferred site of patch accumulation, suggesting that patch dispersion may be random. However, when we plotted patches by size category, we observed that mid-range (red and blue symbols; 0.01–1 cm²) and large-range (violet and green symbols; ≥1 cm²) patches were usually surrounded by several small patches (Figure 2G, Figure S5). Moreover, when focusing on mid-range patch dispersion, these patches appeared to follow unidentified tracks across the skin in several RAG mice (Figure S5). These observations suggested that, while patch dispersion may look different between RAG mice, dispersal is not a random process, but governed by unknown factors in the skin. Our subsequent analysis aimed to test this assertion.

We observed that patches had higher tdTomato signal at their center compared to their periphery. Since tdTomato amastigotes are approximately of equal brightness, it is reasonable to associate the brighter signal at the patch center with higher amastigote density. A grayscale pixel analysis in the red channel to measure tdTomato signal intensity per pixel (Figure 2I, Figure S6) showed that most mean grayscale values per patch were close to the cutoff of 70 (Figure 2J, Figure S2C). However, when patch area was correlated with their respective average grayscale value, it showed a strong correlation (Figure 2H; Spearman r: p < 0.0001, r = 0.6554, 95% CI = 0.6476–0.6631). This suggested that as patches grew outward, they became more densely populated with amastigotes. Cross-section analysis of tdTomato signal intensity in selected patches showed volcano-like graphs with steep shoulders, suggesting that the amastigote bulk is concentrated at the center of a patch and then thins out quickly at the patch periphery (Figures 2K, L).

We also observed that larger patches contained multiple areas of high tdTomato signal (Figure S6). Correlating patch area and density peaks per patch suggested that all mid-range and larger patches had multiple density peaks (Spearman r: p < 0.0001, r = 0.3831, 95% CI = 0.3714–0.3947; Figure 3A). The statistically weak correlation was explained by the broad base of small patches with only one or two peaks. Accumulation of density peaks per patch only became apparent when patch areas exceeded 1 mm² (>10⁶ µm²). This suggested that smaller expanding patches grew into one another to form bigger patches. Furthermore, large patches indicated areas where initial parasite seeding had been higher and/or earlier.

FIGURE 3

Figure 3 Correlation plots. All data presented here were generated by the ImageJ-based analysis of stereomicroscopic images. (A) Correlation plot of 21,904 skin parasite patches by area (µm²) from 23 RAG mice against the number of high-density peaks per individual patch (both on a log-scale; Spearman r: p < 0.0001, r = 0.3831, 95% CI = 0.3714–0.3947). (B) Correlation plot of total parasite patched skin by area (µm²) normalized to 1 cm² of skin against parasite burden/cm² for each of the 23 RAG mice (both on a log-scale; Spearman r: p < 0.0001, r = 0.7589, 95% CI = 0.4948–0.8947). (C) Correlation plot of total skin parasite patch counts normalized to 1 cm² of skin against mean parasite burden per punch biopsy for each of the 23 RAG mice (on a log-scale; Spearman r: p = 0.0021, r = 0.5741, 95% CI = 0.1997–0.8022). (D) Correlation plot of total parasite patched skin by area (µm²) normalized to 1 cm² of skin (on a log-scale) against total skin parasite patch counts normalized to 1 cm² of skin for each of the 23 RAG mice (Spearman r: p < 0.0001, r = 0.8656, 95% CI: 0.6984–0.9432). (E) Correlation plot of the median patch area (µm²) against total skin parasite patch counts for each of the 23 RAG mice (Spearman r: p = 0.006, r = −0.5148, 95% CI = −0.7701 to −0.1175). (F) Correlation plot of total parasite patched skin by area (µm²; on a log-scale) against the skewness of parasite patch area data for each of the 23 RAG mice (Spearman r: p < 0.0001, r = 0.8686, 95% CI = 0.7044–0.9445). (G) Correlation plot of total parasite patched skin by area (µm²; on a log-scale) against the test statistic that infers skewness to the population for each of the 23 RAG mice (Spearman r: p < 0.0001, r = 0.9279, 95% CI = 0.8311–0.9701).

Furthermore, we investigated the relationships between mean parasite burden/cm², patch counts/cm² and total patch area/cm² by correlation plot analysis. Mean parasite burden/cm² correlated strongly with total patch area/cm² of skin, suggesting that increases in total skin parasite burden would also increase total patch area and thus may be the driver of patch growth (Figure 3B; Spearman r: p < 0.0001, r = 0.7589, 95% CI = 0.4948–0.8947). Patch counts/cm² correlated only moderately with parasite burden/biopsy (Figure 3C; Spearman r: p = 0.0021, r = 0.5741, 95% CI = 0.1997–0.8022), suggesting that increases in skin parasite burden was more closely related to patch expansion than new patch seeding. Conversely, when total patch area/cm² was correlated with the patch count/cm², it showed a very strong positive correlation (Figure 3D; Spearman r: p < 0.0001, r = 0.8656, 95% CI: 0.6984–0.9432), suggesting that the proportion of parasitized host skin is largely determined by the number of observed patches. Correlating the median of patch areas per skin for each RAG mouse to the total patch counts per skin, we found a moderately negative correlation (Figure 3E; Spearman r: p = 0.006, r = −0.5148, 95% CI = −0.7701 to −0.1175), which suggested that the more patches a RAG mouse skin had the more of these patches were likely to be small (<0.1 mm²).

While data skewness is an underlying characteristic of skin parasite dispersion for any given sample, skewness cannot be directly inferred onto a population. To test whether the observed sample skewness was a general population feature, we divided sample skewness by the standard error of skewness (SES) (55). The resulting test statistic being outside the interval [-2, 2] infers that data skewness was a population property (Table S1). For instance, total patch area per skin was very strongly correlated with the degree of skewness in the patch area data (Figure 3F; Spearman r: p < 0.0001, r = 0.8686, 95% CI = 0.7044–0.9445), suggesting that as patchiness increased in the skin, mid-range patches give way to a few large patches surrounded by many small patches. This appears to be a common population feature, given the correlation of the test statistic to the total patch area/skin was even stronger than compared to the sample skewness (Figure 3G; Spearman r: p < 0.0001, r = 0.9279, 95% CI = 0.8311–0.9701).

The Skin Patch Landscape Is Heterogeneously Clustered

We previously established that the skin parasite dispersion pertained to concepts of landscape epidemiology (11). Thus, we applied methodologies of spatial point pattern analysis commonly used in landscape epidemiology to characterize further the dispersion of parasites and their myeloid host cells in skin and identify potential modes of their dispersal (16). We converted our patch data into spatial point patterns by mapping the center of mass of each patch to a window characterized by the respective RAG skin outline (Figure 4A, Figure S7 and Table S2). After verifying the absence of point duplications, characteristics of patches (e.g., area and density) could optionally be reintroduced to the pattern as a mark (Figure 4B and Figure S8). Focusing on patch area, we calculated the intensity of point distributions (equivalent to dispersion) with or without a mark and graphically represented the measured point pattern intensities (Figures 4C, D and Figures S9, S10). Alternatively, we calculated the intensity per cm², which is akin to the number of patches per cm² and found that both methods corroborated one another (Tables S1, S2). As shown in Figure 4C (without area mark), intensity showed the distribution of patches by center of mass location, whereas in Figure 4D (with area mark), intensity showed the distribution of patches by size in the skin (see also Figures S9, S10). Both intensities showed a clear departure from a homogeneous distribution. To gain further understanding of how patches related to one another, we performed a cluster analysis, choosing 3.9 mm as the maximum distance between patch borders (see Methods). This showed small patches clustering strongly around larger patches (Figure 4E and Figure S11), confirming our previous observation from Figure 2F (see also Figure S5). Furthermore, this analysis showed the formation of networks of closely clustered patches, which may enhance the likelihood of parasite pickup by a biting sand fly without having to densely parasitize the whole skin (Figure 4E and Figure S11).

FIGURE 4

Figure 4 Spatial point pattern analysis of skin parasite patch distribution. All graphs are examples for spatial point pattern analysis based on the RAG4 data. (A) Spatial point pattern of skin parasite patches based on center of mass of each detected patch without respective marks. (B) Spatial point pattern of skin parasite patches based on center of mass of each detected patch with respective patch area mark. (C) Patch intensity representation of skin parasite patches based on center of mass of each detected patch without respective marks. (D) Patch intensity representation of skin parasite patches based on center of mass of each detected patch with respective patch area mark. (E) Representation of a cluster analysis output linking all patches whose patch borders are ≤3.9 mm into a network. (F) Graphical representation of the isotropy test showing the theoretical Poisson process (red dotted line), the RAG4 data (black solid line), and the global envelopes (gray). (G) Graphical representation of correlation-stationary test based on a locally scaled K-function. The two graphs show the separate distribution analysis of the left and right half of the point shown in panel (A). (H) Scaled L-function showing the theoretical Poisson process (red dotted line), the RAG4 data (black solid line), and the global envelopes (gray). (I) Pair correlation function showing the theoretical Poisson process (red dotted line) and the RAG4 data (black solid line). (J) Data points independence test based on K-function of the data residuals, showing the theoretical Poisson process (red dotted line), the RAG4 data (black solid line), and the upper (green) and lower limits (blue) of the confidence intervals of the theoretical Poisson process. (K) Data points independence test based on nearest-neighbor (G-) function of the data residuals, showing the theoretical Poisson process (red dotted line), the RAG4 data (black solid line), and the upper (green) and lower limits (blue) of the confidence intervals of the theoretical Poisson process.

These observations suggested that patch dispersion is likely to be neither homogeneous nor random. To formally test for homogeneity of patch distributions, a quadrat test based on the Fisher’s exact test was applied, which rejected homogeneity for all patch distributions (p = 0.001; Table S2). Thus, an inhomogeneous Poisson point process (IPP) was used as a null hypothesis. Analogous to a homogeneous Poisson point process, an IPP is assumed to be correlation-stationary (the pair correlation between two spots in the window depends only on their relative positions) and isotropic (the distribution is unaffected by the pattern orientation). However, isotropism was rejected for about half of all patch patterns (Figure 4F and Figure S12); RAG14 did not produce a reliable test result due to its sparse number of patches. To assess correlation-stationarity, we split the point patterns into two portions of equal numbers of observations and applied the inhomogeneous K-function, which is derived from Ripley’s K-function, to each half of the pattern. The analysis showed considerable disagreement between the two plots for several mice indicating that the whole point pattern was not correlation-stationary (Figure S13). In addition, the inhomogeneous K-function showed considerable disagreement between different applied border corrections, which suggested that an IPP was a bad fit to describe patch patterns. In addition, we rejected a Poisson distribution for our patch patterns based on a Quadrat Count test (Table S2). These disagreements suggested that the spatial scale of point interactions, which is assumed to be constant by the inhomogeneous K-function (16), was locally variable for our patch patterns. Thus, by assuming a point process that was locally stationary and isotropic but showed a varying rescaling factor between locations in the skin, we could use a locally scaled K-function to rerun the correlation-stationary test. Although not perfect, agreement between the two plots per mouse and between different border corrections was significantly improved (Figure 4G and Figure S14), suggesting that patch distribution was locally scaled.

We proceeded with the scaled L-function, derived from Besag’s L-function, to determine if patches were evenly spaced (regular), randomly spaced (independent), or clustered (Figure 4H and Figure S15). While the fitted data (black solid line) generally departed from the hypothetical random distribution toward a clustered one (black solid line above dotted red line), only RAG11, RAG13, RAG14, and RAG23 showed a clear departure from the global simulation envelops of a hypothetical random distribution (gray shaded area) toward clustering. To confirm these results, we used the pair correlation function (PCF). The scaled PCF analysis showed for most mice, patch patterns departed from the theoretical random distribution (dotted red line) toward a clustered point process (black line above the red dotted line), in particular, at short distances between patches (Figure 4I and Figure S16). With increased distances between patches, some point patterns showed an approximation to more point independence. Applying the K- and nearest-neighbor (G-) function to the point pattern residuals showed a clear departure of independence toward clustering at shorter distances, indicating positive point dependence [black line outside the upper confidence intervals (dotted green line)] for both tests (Figures 4J, K and Figures S17, S18).

Modeling the Point Process Behind Patch Patterns Suggests Influence of Covariates

To gain further mechanistic insights into patch distribution, we used a model that accounted for positive point (patch) dependence (clustering) and was tolerant to a locally varying distribution scale. We used different cluster processes (Cauchy, Matérn, Thomas and Variance-Gamma) as described by Baddeley et al. [“spatstat” (16)] that assume cluster formation based on a parental point distribution that seeded the observed points around themselves within a radius, which is variable between clusters due to the locally scaled nature of cluster distribution. This assumption is reflected in our data showing that small patches cluster around mid-range and large patches (Figures 2F, 4E and Figures S5, S11).

Thirteen of the 23 patch distributions were best fitted by a Matérn cluster process, while the other 10 were best fitted by a Thomas cluster process. The Matérn cluster process assumes a uniform distribution of offspring around the parental point, while the Thomas cluster process assumes offspring to be randomly displaced from the parental point (16). Checking the parameters selected by the model showed that they were within a reasonable range, confirming model validity (Table S3). This was further confirmed by the scaled L-function of the selected fitted models for each patch pattern (Figure 5A and Figure S19). The L-functions of the fitted models closely resembled the L-functions from the point pattern analyses (Figure 4H and Figure S15). Furthermore, we used the selected fitted models to run 4 simulations each to see if the outcome would resemble the distributions in the original data (Figure 5B and Figure S20). With a few exceptions, like RAG3, the simulation outputs provided a good representation of the original patch dispersion data (Figure 4A and Figure S7).

FIGURE 5

Figure 5 Basic model evaluation of skin parasite patch distribution. (A) Scaled L-function showing the theoretical Poisson process (red dotted line), the RAG4 predicted model data (black solid line), and the global envelopes (gray). (B) Four model simulations based on the Matérn cluster process for RAG4.

Based on these findings, we can deduce how the dispersal of parasites and their myeloid host cells may operate. Both the Matérn and Thomas cluster processes assume points of origin that seed new points within a set radius to create an observed clustered point dispersion. In a locally scaled setting, as it is the case for our patch distribution, this radius is locally varying in length and thus distinct for each patch cluster. Analogously, after initial skin infection, the dispersal of patches may be explained by existing patches locally seeding new patches around themselves, while growing outward, creating a self-propagating network of patch clusters in the skin (Figure 6 and Figure S21).

FIGURE 6

Figure 6 Proposed mechanism for the dispersal of patches of L. donovani-infected myeloid cells. (Step 1) After an initial seeding event of Leishmania amastigote-infected host cells (in red) in the skin from the circulation, uninfected phagocytic cells (in purple) in the skin are attracted to the infected cells. (Step 2) These recruited myeloid cells form a type of innate granuloma in which suitable uninfected myeloid cells are gradually infected. This results in the formation of a primary patch of infected myeloid cells. (Step 3) As these patches keep growing, they seed new patches within a radius “r” around themselves, potentially by escape of infected myeloid cells or via transient release of free amastigotes. This process then keeps repeating itself, forming self-propagating networks of patch “clusters”. The image was created with BioRender.com.

Discussion

It is known that heterogeneity in skin parasite dispersion provides a mechanism to enhance host infectiousness to sand fly vectors (11). Here, we have provided evidence for the underlying mechanisms driving this heterogeneity. To our knowledge, this is the first study to make use of spatial point pattern methodologies (16) in this way. The Matérn and Thomas cluster processes, which were identified to best describe how dispersion of patches of parasites and their infected myeloid host cells occurred, both assume points of origin that seed new points within locally varying radii clustering around them. Based on this analysis, we therefore propose that existing patches locally seed new patches around themselves, giving rise to patch clusters with a central large patch. Since patches are densely parasitized, these patch networks would increase the chance of high parasite loads being acquired by a sand fly and eliminate the need for parasitism of the entire skin. Indeed, our conservatively chosen tdTomato signal and patch size thresholds may underestimate very small patches leading to patch networks that are denser than reported here, further enhancing outward transmission potential.

In vivo, Leishmania amastigotes reside in phagocytic cells in the dermis, primarily macrophages (56). Therefore, parasite dispersion is linked to host cell dispersion, a subject on which scant information is available. In the steady-state dermis, macrophages and DCs are the primary phagocytic cell populations (51). Interestingly, macrophages and DCs have very particular distributions in steady-state skin. Branched DCs are predominantly interstitially localized directly underneath the dermal–epidermal junction (57, 58). Macrophages are commonly found below these DCs, mostly interstitially. Deeper in the reticular dermis, both cell types associate perivascular to blood vessels rather than lymphatic vessels and DCs here are more rounded. Among the different macrophage populations in the skin, it was shown that Man^R+ TRM2s are important as long-term L. major host cells (52). In addition, the rate of parasite infection in Man^R+ TRM2s (~45%) correlated strongly with the severity of skin pathology and parasite survival in the host. The rate of infection in DCs was much lower (22.5%), indicating a preference of L. major for Man^R+ TRM2s (52). In our study, the observed infection rates for both cell types with L. donovani were less than half of this (~19.1% and 9.8%, respectively), but L. donovani similarly showed a preference for Man^R+ TRM2s over DCs (5:1). These cells constituted almost 50% of all infected host cells. Conversely, 50% of host cells belonged to neither group; presumably, other myeloid-derived phagocyte populations (56). Interestingly, the majority of Man^R+ TRM2s and CD11c⁺ DCs (~80.9% and ~90.2%, respectively) were uninfected in the analyzed skin sections, despite the high number of parasites present in some of the sections. This raised the question of why the small proportion of infected host cells should be clustered close together and also indicates that observed cellular tropism is not due to saturation of preferred targets.

Furthermore, our analysis suggested that the process of forming patches may be influenced by yet-to-be-identified covariates. One covariate may be proximity to blood vessels as Man^R+ TRM2s often colocalize to them in the dermis, where they can pick up macromolecules from the blood (52, 59). Given that parasitized phagocytic cells in circulation may have properties associated with apoptotic cells (60), it is possible that TRMs capture parasitized host cells by efferocytosis, thus seeding new skin sites for patch formation. Alternatively, infected phagocytes in circulation could exit the vasculature by extravasation at sites of sub-clinical inflammation. Either mechanism could explain why mid-range patch distribution appeared to follow unidentified lines in the skin and why there was no preferred location for parasite accumulation in the skin. Furthermore, proximity to blood vessels may also enhance access to nutrients from the circulation for parasitized phagocytic cells and thus indirectly promote parasite proliferation. Also, parasite proximity to blood vessels would improve their chances to be found by a probing sand fly. Future analysis could potentially confirm blood vessel association of parasite patches by proximity measurements. Three-dimensional (e.g., light sheet) microscopic imaging data would be needed to achieve this. This and other covariates that may affect the establishment of parasite niches can be integrated in the future into the Matérn and Thomas cluster process models to determine their influence on the dispersal and dispersion of parasites and their myeloid host cells.

Independent of the means of initial skin seeding of parasites, our data indicate that the accumulation of infected myeloid cells into patches could be analogous mechanistically to phagocyte migration during granulomatous inflammation. For example, early in the process of L. donovani-dependent granuloma formation in the liver, uninfected Kupffer cells (KCs) migrate towards a granuloma-initiating infected KC resulting in a heterogeneous distribution of KCs and their condensation and subsequent infection at the granuloma core (18, 61). Conversely, studies in zebrafish embryos infected with Mycobacterium marinum highlight the innate capacity of phagocytes to exit granulomas (62), an innate feature of the phagocytic cell response to infection that appears to be obscured in fully immunocompetent mouse models of granulomatous inflammation (61). Our RAG mouse model of skin infection similarly lacks the confounding effects of an adaptive immune response, providing a similarly unique window on processes driving innate myeloid cell aggregation and a platform on which to build increasing levels of immune complexity. It remains to be determined to what extent patches form because of the inherent migratory properties of phagocytic cells and/or modification of such properties by intracellular amastigotes. Manipulation of phagocyte function by intracellular Leishmania has been well-described (63) though few studies have addressed phagocyte migratory potential. Similarly, while phagocyte heterogeneity in skin is well described (see above), the mechanisms regulating mobility of these cells in their natural 3D-environment remain poorly understood.

The notion that new patches are seeded primarily from existing patches is supported by the observation that larger patches were surrounded by smaller patches that emerged within limited distances from the larger patches, resulting in patch clustering. If, akin to parasite dispersion, parasite dispersal is also linked to its host cells, then, mechanistically, patch clustering could occur by escape of infected host cells from these “innate granuloma”-like parasite patches as was demonstrated in M. marinum-infected zebrafish embryos (64). Conversely, parasites could be redistributed to uninfected host cells within a limited radius from patches by amastigote release from ruptured host cells or by intercellular passage. Macrophages are known to form a variety of open-ended tunneling nanotubes (TNT) between one another at a distance and such structures have been implicated in the spread of several respiratory viruses and HIV-1 (65–67). We know of no formal evidence, however, to show that Leishmania parasites are transported between macrophages via TNTs.

The cores of patches are very dense with parasites and heavily infected host cells, but the patch fringes are much less heavily parasitized, suggesting that parasites proliferating strongly within a patch are then pushed outward. Mandell and Beverley (68) observed that not all amastigotes in the host skin proliferate at the same rate and described two distinct amastigote populations; one that was dormant, the other proliferative (68). Thus, patches may represent areas where amastigotes are proliferative, while inter-patch space could mark areas of dormant amastigotes. Whether host cell metabolism stimulates/supports amastigotes to proliferate in the mammalian host skin (69) remains to be determined.

In summary, through spatial point pattern analysis combined with systematic image analysis, we provide a novel mechanistic framework to explain how the dispersal of L. donovani-infected myeloid cells occurs in the absence of the constraints imposed by acquired immunity.

Data Availability Statement

All Fiji ImageJ macros and R scripts used in this study were deposited on GitHub (https://github.com/joedoehl/LeishSkinSPPM) and are freely available for download.

Ethics Statement

All animal usage was approved by the University of York Animal Welfare and Ethics Review Committee and the Ethics Review Committee at FERA Science Ltd., and performed under UK Home Office license (“Immunity and Immunopathology of Leishmaniasis” Ref # PPL 60/4377 and P49487014).

Author Contributions

JD and PK designed the study. JD, HA, NB, and AR conducted the experimental studies. JD conducted the spatial point pattern analysis, modeling, and all the statistical analysis, which were revised by SC and JP. The manuscript was written by JD and PK, and all authors contributed to the final version. PK supervised the study and has senior authorship. PK conceived the study. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by a Wellcome Senior Investigator Award (# WT104726; http://www.wellcome.co.uk) and a MRC Programme Grant (#G1000230; http://www.mrc.ac.uk) to PMK.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank Andrew McNaughton from the Confocal/µCT Facility at University of Otago, NZ, for the permission of adapting his protocol and invaluable discussion, Jesus Valenzuela from the VMBU-LMVR at the NIAID/NIH, USA, for encouragement and time made available to complete the manuscript, the Imaging and Cytometry lab and staff, in particular Graeme Park and Jo Marrison, in the Bioscience Technology Facility at the University of York, UK, for use of their equipment and technically knowhow that supported the image acquisition, and the staff of the Biological Services Facility at the University of York, UK, for animal husbandry.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.795554/full#supplementary-material

References

1. Ready PD. Biology of Phlebotomine Sand Flies as Vectors of Disease Agents. Annu Rev Entomology (2013) 58(1):227–50. doi: 10.1146/annurev-ento-120811-153557

CrossRef Full Text | Google Scholar

2. Bates PA. Revising Leishmania's Life Cycle. Nat Microbiol (2018) 3(5):529–30. doi: 10.1038/s41564-018-0154-2

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Murray CJL, Vos T, Lozano R, Naghavi M, Flaxman AD, Michaud C, et al. Disability-Adjusted Life Years (DALYs) for 291 Diseases and Injuries in 21 Regions, 1990–2010: A Systematic Analysis for the Global Burden of Disease Study 2010. Lancet (2012) 380(9859):2197–223. doi: 10.1016/S0140-6736(12)61689-4

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Alvar J, Vélez ID, Bern C, Herrero M, Desjeux P, Cano J, et al. Leishmaniasis Worldwide and Global Estimates of Its Incidence. PloS One (2012) 7(5):e35671. doi: 10.1371/journal.pone.0035671

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Srivastava S, Shankar P, Mishra J, Singh S. Possibilities and Challenges for Developing a Successful Vaccine for Leishmaniasis. Parasites Vectors (2016) 9(1):277. doi: 10.1186/s13071-016-1553-y

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Singh N, Kumar M, Singh RK. Leishmaniasis: Current Status of Available Drugs and New Potential Drug Targets. Asian Pacific J Trop Med (2012) 5(6):485–97. doi: 10.1016/S1995-7645(12)60084-4

CrossRef Full Text | Google Scholar

7. Ready P. Epidemiology of Visceral Leishmaniasis. Clin Epidemiol (2014) 6:147–54. doi: 10.2147/CLEP.S44267

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Miller E, Warburg A, Novikov I, Hailu A, Volf P, Seblova V, et al. Quantifying the Contribution of Hosts With Different Parasite Concentrations to the Transmission of Visceral Leishmaniasis in Ethiopia. PloS Neglected Trop Dis (2014) 8(10):e3288. doi: 10.1371/journal.pntd.0003288

CrossRef Full Text | Google Scholar

9. Singh OP, Tiwary P, Kushwaha AK, Singh SK, Singh DK, Lawyer P, et al. Xenodiagnosis to Evaluate the Infectiousness of Humans to Sandflies in an Area Endemic for Visceral Leishmaniasis in Bihar, India: A Transmission-Dynamics Study. Lancet Microbe (2021) 2(1):e23–31. doi: 10.1016/S2666-5247(20)30166-X

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kirstein OD, Abbasi I, Horwitz BZ, Skrip L, Hailu A, Jaffe C, et al. Minimally Invasive Microbiopsies: A Novel Sampling Method for Identifying Asymptomatic, Potentially Infectious Carriers of Leishmania Donovani. Int J Parasitol (2017) 47(10-11):609–16. doi: 10.1016/j.ijpara.2017.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Doehl JSP, Bright Z, Dey S, Davies H, Magson J, Brown N, et al. Skin Parasite Landscape Determines Host Infectiousness in Visceral Leishmaniasis. Nat Commun (2017) 8(1):57. doi: 10.1038/s41467-017-00103-8

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Valverde JG, Paun A, Inbar E, Romano A, Lewis M, Ghosh K, et al. Increased Transmissibility of Leishmania Donovani From the Mammalian Host to Vector Sand Flies After Multiple Exposures to Sand Fly Bites. J Infect Dis (2017) 215(8):1285–93. doi: 10.1093/infdis/jix115

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Laurenti MD, Rossi CN, da Matta VL, Tomokane TY, Corbett CE, Secundino NF, et al. Asymptomatic Dogs Are Highly Competent to Transmit Leishmania (Leishmania) Infantum Chagasi to the Natural Vector. Vet Parasitol (2013) 196(3-4):296–300. doi: 10.1016/j.vetpar.2013.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Molina R, Ghosh D, Carrillo E, Monnerat S, Bern C, Mondal D, et al. Infectivity of Post-Kala-Azar Dermal Leishmaniasis Patients to Sand Flies: Revisiting a Proof of Concept in the Context of the Kala-Azar Elimination Program in the Indian Subcontinent. Clin Infect Dis (2017) 65(1):150–3. doi: 10.1093/cid/cix245

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Mondal D, Bern C, Ghosh D, Rashid M, Molina R, Chowdhury R, et al. Quantifying the Infectiousness of Post-Kala-Azar Dermal Leishmaniasis Toward Sand Flies. Clin Infect Dis (2019) 69(2):251–8. doi: 10.1093/cid/ciy891

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Baddeley A, Rubak E, Turner R. Spatial Point Patterns: Methodology and Applications With R. London: Chapman and Hall/CRC Press (2015).

Google Scholar

17. Shaner NC, Campbell RE, Steinbach PA, Giepmans BNG, Palmer AE, Tsien RY. Improved Monomeric Red, Orange and Yellow Fluorescent Proteins Derived From Discosoma Sp. Red Fluorescent Protein. Nat Biotechnol (2004) 22(12):1567–72. doi: 10.1038/nbt1037

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Beattie L, Peltan A, Maroof A, Kirby A, Brown N, Coles M, et al. Dynamic Imaging of Experimental Leishmania Donovani-Induced Hepatic Granulomas Detects Kupffer Cell-Restricted Antigen Presentation to Antigen-Specific CD8+ T Cells. PloS Pathog (2010) 6(3):e1000805. doi: 10.1371/journal.ppat.1000805

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bezerra-Vasconcelos DR, Melo LM, Albuquerque ÉS, Luciano MCS, Bevilaqua CML. Real-Time PCR to Assess the Leishmania Load in Lutzomyia Longipalpis Sand Flies: Screening of Target Genes and Assessment of Quantitative Methods. Exp Parasitology (2011) 129(3):234–9. doi: 10.1016/j.exppara.2011.08.010

CrossRef Full Text | Google Scholar

20. Accart N, Sergi F, Rooke R. Revisiting Fixation and Embedding Techniques for Optimal Detection of Dendritic Cell Subsets in Tissues. J Histochem Cytochem (2014) 62(9):661–71. doi: 10.1369/0022155414539963

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Dalton JE, Maroof A, Owens BM, Narang P, Johnson K, Brown N, et al. Inhibition of Receptor Tyrosine Kinases Restores Immunocompetence and Improves Immune-Dependent Chemotherapy Against Experimental Leishmaniasis in Mice. J Clin Invest (2010) 120(4):1204–16. doi: 10.1172/JCI41281

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: An Open-Source Platform for Biological-Image Analysis. Nat Methods (2012) 9(7):676–82. doi: 10.1038/nmeth.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Sabino CP, Deana AM, Yoshimura TM, da Silva DF, Franca CM, Hamblin MR, et al. The Optical Properties of Mouse Skin in the Visible and Near Infrared Spectral Regions. J Photochem Photobiol B (2016) 160:72–8. doi: 10.1016/j.jphotobiol.2016.03.047

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Liao P-S, Chen T-S, Chung P-C. A Fast Algorithm for Multilevel Thresholding. J Inf Sci Eng (2001) 17(5):713–27. doi: 10.6688/JISE.2001.17.5.1

CrossRef Full Text | Google Scholar

25. Podinovskaia M, Descoteaux A. Leishmania and the Macrophage: A Multifaceted Interaction. Future Microbiol (2015) 10(1):111–29. doi: 10.2217/fmb.14.103

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Guertin DA, Sabatini DM. Cell Size Control. Encyclopedia Life Sci (2006). doi: 10.1038/npg.els.0003359

CrossRef Full Text | Google Scholar

27. Haeri M, Haeri M. ImageJ Plugin for Analysis of Porous Scaffolds Used in Tissue Engineering. J Open Res Software (2015) 3:e1. doi: 10.5334/jors.bn

CrossRef Full Text | Google Scholar

28. Tupper B. Graph Plugin and Demo Nih2010. Available at: https://imagej.nih.gov/ij/plugins/graph/index.html.

Google Scholar

29. Corp I. IBM SPSS Statistics for Windows. 25.0. Armonk NY, editor. IBM Corp (2017).

Google Scholar

30. RStudio Team. RStudio: Integrated Development for R. 1.4. Boston MA, editor. RStudio, PBC (2020).

Google Scholar

31. Core Team R. R: A Language and Environment for Statistical Computing. 4.0.2. Vienna, Austria: R Foundation for Statistical Computing (2020).

Google Scholar

32. Plate T, Heiberger R. Abind: Combine Multidimensional Arrays. 1.4-5. CRAN (2016).

Google Scholar

33. Turner R. Deldir: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation. 0.1-29. CRAN (2020).

Google Scholar

34. Meyer D, Dimitriadou E, Hornik K, Weingessel A, Leisch F. E1071: Misc Functions of the Department of Statistics. Wien TU, editor. CRAN (2020). In: E1071) PTGF.

Google Scholar

35. Kassambara A. Ggpubr: 'Ggplot2' Based Publication Ready Plots. 0.4.0. CRAN (2020).

Google Scholar

36. Faraway J, Marsaglia G, Marsaglia J, Baddeley A. Goftest: Classical Goodness-Of-Fit Tests for Univariate Distributions. 1.2-2. CRAN (2019).

Google Scholar

37. Johnson A, Baddeley A. Polyclip: Polygon Clipping. 1.10-0. CRAN (2019).

Google Scholar

38. Schlather M, Malinowski A, Oesting M, Boecker D, Strokorb K, Engelke S, et al. RandomFields: Simulation and Analysis of Random Fields. CRAN (2020).

Google Scholar

39. Ren K. Rlist: A Toolbox for Non-Tabular Data Manipulation. CRAN (2016).

Google Scholar

40. Mukherji S. Rms.Gof: Root-Mean-Square Goodness-of-Fit Test for Simple Null Hypothesis. CRAN (2013).

Google Scholar

41. Wickham H, Seidel D. Scales: Scale Functions for Visualization. CRAN (2020).

Google Scholar

42. Bivand RS, Pebesma E, Gomez-Rubio V. Applied Spatial Data Analysis With R. New York: Springer (2013). 405 p.

Google Scholar

43. Davies TM, Marshall JC, Hazelton ML. Tutorial on Kernel Estimation of Continuous Spatial and Spatiotemporal Relative Risk. Stat Med (2018) 37(7):1191–221. doi: 10.1002/sim.7577

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Baddeley A, Turner R, Rubak E. Spatstat.Utils: Utility Functions for 'Spatstat'. 1.17-0. CRAN (2020).

Google Scholar

45. Rougier J. Tensor: Tensor Product of Arrays. 1.5. CRAN (2012).

Google Scholar

46. Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. Welcome to the Tidyverse. J Open Source Software (2019) 4(43):1686. doi: 10.21105/joss.01686

CrossRef Full Text | Google Scholar

47. Lessler J, Salje H, Grabowski MK, Cummings DA. Measuring Spatial Dependence for Infectious Disease Epidemiology. PloS One (2016) 11(5):e0155249. doi: 10.1371/journal.pone.0155249

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Ohser J, Muecklich F. Arrangement of Objects. Statistical Analysis of Microstructures in Materials Science. Chichester, England: John Wiley (2000) p. 263–312.

Google Scholar

49. Ripley BD. Modeling Spatial Patterns. J Roy Stat Soc B Met (1977) 39(2):172–212. doi: 10.1111/j.2517-6161.1977.tb01615.x

CrossRef Full Text | Google Scholar

50. Iglesias-Guitian JA, Aliaga C, Jarabo A, Gutierrez D. A Biophysically-Based Model of the Optical Properties of Skin Aging. Comput Graphics Forum (2015) 34(2):45–55. doi: 10.1111/cgf.12540

CrossRef Full Text | Google Scholar

51. Dupasquier M, Stoitzner P, van Oudenaren A, Romani N, Leenen PJ. Macrophages and Dendritic Cells Constitute a Major Subpopulation of Cells in the Mouse Dermis. J Invest Dermatol (2004) 123(5):876–9. doi: 10.1111/j.0022-202X.2004.23427.x

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lee SH, Charmoy M, Romano A, Paun A, Chaves MM, Cope FO, et al. Mannose Receptor High, M2 Dermal Macrophages Mediate Nonhealing Leishmania Major Infection in a Th1 Immune Environment. J Exp Med (2018) 215(1):357–75. doi: 10.1084/jem.20171389

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Merad M, Sathe P, Helft J, Miller J, Mortha A. The Dendritic Cell Lineage: Ontogeny and Function of Dendritic Cells and Their Subsets in the Steady State and the Inflamed Setting. Annu Rev Immunol (2013) 31:563–604. doi: 10.1146/annurev-immunol-020711-074950

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Mai LT, Smans M, Silva-Barrios S, Fabie A, Stager S. IRF-5 Expression in Myeloid Cells Is Required for Splenomegaly in L. Donovani Infected Mice. Front Immunol (2019) 10:3071. doi: 10.3389/fimmu.2019.03071

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Duncan C. Fundamental Statistics for Social Research. Routledge (1997).

Google Scholar

56. Liu D, Uzonna JE. The Early Interaction of Leishmania With Macrophages and Dendritic Cells and Its Influence on the Host Immune Response. Front Cell Infect Microbiol (2012) 2:83. doi: 10.3389/fcimb.2012.00083

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Wang XN, McGovern N, Gunawan M, Richardson C, Windebank M, Siah TW, et al. A Three-Dimensional Atlas of Human Dermal Leukocytes, Lymphatics, and Blood Vessels. J Invest Dermatol (2014) 134(4):965–74. doi: 10.1038/jid.2013.481

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Tong PL, Roediger B, Kolesnikoff N, Biro M, Tay SS, Jain R, et al. The Skin Immune Atlas: Three-Dimensional Analysis of Cutaneous Leukocyte Subsets by Multiphoton Microscopy. J Invest Dermatol (2015) 135(1):84–93. doi: 10.1038/jid.2014.289

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Barreiro O, Cibrian D, Clemente C, Alvarez D, Moreno V, Valiente I, et al. Pivotal Role for Skin Transendothelial Radio-Resistant Anti-Inflammatory Macrophages in Tissue Repair. Elife (2016) 5:e15251. doi: 10.7554/eLife.15251

CrossRef Full Text | Google Scholar

60. Martinez-Lopez M, Soto M, Iborra S, Sancho D. Leishmania Hijacks Myeloid Cells for Immune Escape. Front Microbiol (2018) 9:883. doi: 10.3389/fmicb.2018.00883

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Egen JG, Rothfuchs AG, Feng CG, Winter N, Sher A, Germain RN. Macrophage and T Cell Dynamics During the Development and Disintegration of Mycobacterial Granulomas. Immunity (2008) 28(2):271–84. doi: 10.1016/j.immuni.2007.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Davis JM, Clay H, Lewis JL, Ghori N, Herbomel P, Ramakrishnan L. Real-Time Visualization of Mycobacterium-Macrophage Interactions Leading to Initiation of Granuloma Formation in Zebrafish Embryos. Immunity (2002) 17(6):693–702. doi: 10.1016/S1074-7613(02)00475-2

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Kaye P, Scott P. Leishmaniasis: Complexity at the Host-Pathogen Interface. Nat Rev Microbiol (2011) 9(8):604–15. doi: 10.1038/nrmicro2608

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Davis JM, Ramakrishnan L. The Role of the Granuloma in Expansion and Dissemination of Early Tuberculous Infection. Cell (2009) 136(1):37–49. doi: 10.1016/j.cell.2008.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Onfelt B, Nedvetzki S, Benninger RK, Purbhoo MA, Sowinski S, Hume AN, et al. Structurally Distinct Membrane Nanotubes Between Human Macrophages Support Long-Distance Vesicular Traffic or Surfing of Bacteria. J Immunol (2006) 177(12):8476–83. doi: 10.4049/jimmunol.177.12.8476

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Cifuentes-Munoz N, Dutch RE, Cattaneo R. Direct Cell-to-Cell Transmission of Respiratory Viruses: The Fast Lanes. PloS Pathog (2018) 14(6):e1007015. doi: 10.1371/journal.ppat.1007015

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Bracq L, Xie M, Benichou S, Bouchet J. Mechanisms for Cell-To-Cell Transmission of HIV-1. Front Immunol (2018) 9:260. doi: 10.3389/fimmu.2018.00260

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Mandell MA, Beverley SM. Continual Renewal and Replication of Persistent Leishmania Major Parasites in Concomitantly Immune Hosts. Proc Natl Acad Sci USA (2017) 114(5):E801–10. doi: 10.1073/pnas.1619265114

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Saunders EC, McConville MJ. Immunometabolism of Leishmania Granulomas. Immunol Cell Biol (2020) 98:832–44. doi: 10.1111/imcb.12394

PubMed Abstract | CrossRef Full Text | Google Scholar