A spatially coordinated keratinocyte–fibroblast circuit recruits MMP9⁺ myeloid cells to drive type I interferon-driven inflammation in photosensitive autoimmunity

A spatially coordinated keratinocyte–fibroblast circuit recruits MMP9⁺ myeloid cells to drive type I interferon-driven inflammation in photosensitive autoimmunity

Study design

Study participants were recruited under IRB-approved protocols (H00021295, H-14848) at the University of Massachusetts Medical School. Written informed consent was obtained from all participants (or their legally authorized representatives, when applicable) before enrollment. Individuals with vitiligo, psoriasis, CLE or DM were diagnosed by board-certified dermatologists, and age-matched and sex-matched healthy volunteers without autoimmune or inflammatory skin disease served as controls. Skin samples included healthy skin, lesional skin and non-lesional skin (clinically normal tissue ≥10 cm from lesions). Samples were used for scRNA-seq, flow cytometry, proteomics and spatial transcriptomics as appropriate (Supplementary Table 1).

Suction blister biopsies and processing

Suction blisters were generated as previously described43 using 10–15 mmHg negative pressure at 40 °C for 30–60 min. Blister fluid was aspirated and centrifuged at 350g for 10 min at 4 °C; the supernatant was then flash-frozen for protein analysis. The cell pellet was processed for inDrop scRNA-seq and flow cytometry.

InDrop scRNA-seq

Cells from blister fluid were processed using the inDrop protocol44, with an extended 8 bp unique molecular identifier (UMI). Libraries were sequenced on Illumina NovaSeq or NextSeq, and raw FASTQ files were processed using the inDrop pipeline (https://github.com/garber-lab/indrop_processing). Sequences were aligned to GRCh38 with STAR45, and digital gene expression matrices were generated using the End Sequencing Analysis Toolkit (ESAT)45,46. A custom GTF file was applied to merge genes with overlapping 3′ ends. Additional details are provided in Supplementary Methods.

scRNA-seq analysis

The UMI count tables were processed using Seurat (v.5)47. Cells with 200–4,000 UMIs and <20 hemoglobin UMIs (HBA1, HBA2, HBB, HBD) were retained. Mitochondrial and ribosomal genes (RPL*, RPS*) were excluded from the 2,000 variable genes before performing principal component analysis (PCA). The top 20–30 PCs were used for uniform manifold approximation and projection (UMAP) and clustering (FindNeighbors/FindClusters). Cells were first grouped into major lineages, followed by subtype reclustering. Doublet clusters, identified by inappropriate marker co-expression, were removed before final clustering of lymphoid and myeloid compartments.

For lymphocytes, 130 cell cycle genes were excluded from PCA. Interferon-response scores were derived from 50 high-confidence interferon-induced genes (log2(fold change)-ranked, FDR < 0.01) and regressed out during scaling. Visualizations were constructed using ggplot2 (ref. 48).

Differential expression was performed on pseudo-bulk profiles normalized to 1 M UMIs per sample and analyzed with DESeq2 (ref. 49) (design: ~skin; comparisons: healthy control versus lesional and healthy control versus non-lesional within disease and lineage). Differentially expressed genes required P < 0.01 and |log2(fold change)| > 1. Subtype-level differential expression used edgeR50 adapted for Seurat3,50,51, with significance assessed at FDR < 0.05. Heatmaps were generated with ComplexHeatmap52.

Compositional differences were assessed with scCODA53, which models dependencies among cell type proportions. KC_suprabasal was used as the reference because of its relative abundance stability. Proportion heatmaps show disease-level percentages scaled by cell type.

Visualization functions are packed in the R package AddOns, available on GitHub (https://github.com/garber-lab/AddOns).

10× scRNA-seq of three DM samples analysis

Raw count matrices were normalized using the NormalizeData function and the LogNormalize method in Seurat (v.5). Data were integrated with the CCAIntegration method on the top 30 PCs, using the top 2,000 variable genes selected with the FindVariableFeatures function, with default parameters. UMAP and clustering were performed on the integrated top 30 PCs.

GSE179633 scRNA-seq analysis

Data from a previous publication24 (GSE179633) were obtained as raw UMI count matrices and analyzed using Seurat (v.5)47and AddOns. In total, we retained 13 epidermal samples (three healthy controls, five DLE, five SLE) and 16 dermal samples (four healthy controls, five DLE, seven SLE).

Each sample underwent standard Seurat quality control. Doublets were identified with both scDblFinder54 and DoubletFinder55, and any cell flagged by either tool was removed. Cell cycle effects were regressed out, and Harmony was applied to integrate dermal and epidermal datasets. Cells were first grouped into major lineages, followed by subtype reclustering.

MC_CD14 cells were extracted and embedded separately. All dermal and epidermal samples were combined, and the top 500 variable genes (vst method) were used for PCA. The first five PCs were integrated with Harmony (k.weight = 30) and used to compute a diffusion map (destiny56). The first diffusion component defined the LYVE1-to-MMP9 continuum; cells were ranked along this axis. Heatmaps were generated using ComplexHeatmap52. For trajectory heatmaps, MC_CD14 cells were divided into 50 equal bins, and mean expression per bin was calculated. Expression values were scaled per gene across all bins and reference cell types.

Proteomic analysis

The NULISA (Nucleic Acid Linked Immuno-Sandwich Assay) Inflammatory Panel (Supplementary Table 2) quantified 250 proteins in blister fluid samples. Only samples with NULISA protein quantification values above the level of detection were included in visualization and differential expression analysis.

Cytokine concentrations in cell culture media were measured using Olink’s Proximity Extension Assay (Target 48 Cytokine/Immune Surveillance Bundle; Supplementary Table 6). For moDCs, secretion levels were calculated by subtracting baseline protein levels in the input media from post-stimulation concentrations.

Statistical analyses were performed using two-sided t-tests. Differential protein abundance was assessed using the olink_ttest function from the OlinkAnalyze package (https://doi.org/10.32614/CRAN.package.OlinkAnalyze) with unpaired comparisons. The BoxPlot.nulisa function in the AddOns package was used for visualization (https://github.com/garber-lab/AddOns).

In vivo UVB photoprovocation

UVB photoprovocation was performed under IRB approval (H00021295), using a standardized protocol57. On day −1, six sites on the lower back were irradiated with UVB (25–150 mJ cm−2) using a SolRx 100-Series handheld device (285–350 nm, Philips PL-S 9W/12/2P). Then, 24 h later (day 0), the minimal erythema dose was determined for each participant.

A dose of 1.5× the minimal erythema dose was used for photoprovocation. If this dose fell within the initial testing range (≤150 mJ cm−2), samples were obtained 24 h (suction blister) or 72 h (punch biopsy) after threshold testing without additional irradiation. For participants whose 1.5× minimal erythema dose exceeded 150 mJ cm−2, a second visit was scheduled to administer the full dose, followed by sampling at the same time points.

UVB-treated and matched untreated sites were sampled by suction blistering or punch biopsy. In a subset of participants, photoprovocation was performed before and after three monthly 300 mg intravenous doses of anifrolumab to assess treatment-associated changes in UVB responsiveness. Cutaneous disease activity was scored using CLASI58.

T cell isolation, activation and cytokine treatment

CD4+ T cells were isolated from healthy donor peripheral blood mononuclear cells (PBMCs) using the CD4+ T Cell Isolation Kit (Miltenyi Biotec). The cells were resuspended in T cell culture media (RPMI supplemented with HEPES, penicillin–streptomycin and FBS) and plated in 24-well Nunc cell culture plates.

For T cell activation, wells were pre-coated overnight at 4 °C with anti-CD3 (10 ng ml−1; BioLegend, Ultra-LEAF Purified anti-human, Clone OKT3) and anti-CD28 (2 ng ml−1; BioLegend, Ultra-LEAF Purified anti-human, Clone CD28.2) antibodies dissolved in PBS.

After plating, the cells were treated with cytokines—either IFNβ (50 ng ml−1; R&D Systems) or IFNγ (100 ng ml−1; R&D Systems)—and cultured for 3 days. Following the incubation, cells were stained for viability and surface markers CD2 and CD4. Subsequently, the cells were fixed, permeabilized and stained intracellularly for granzyme B (GZMB, BioLegend, Clone QA16A02) before being analyzed by flow cytometry.

Cytokilling assay

The cytokilling assay was performed by measuring impedance of the cells using the xCELLigence Real Time Cell Analysis technology (RTCA, Agilent). A total of 20,000 BT-20 cells (ATCC, HTB-19) were seeded in an RTCA E-96 well plate and left to adhere for 6 h. Then, 2 × 105 freshly isolated CD4 T cells from PBMCs were added. Solitomab (Sigma-Aldrich), a bispecific antibody anti-CD3 and anti-epithelial cell adhesion molecule, was used to stimulate the cells at a concentration of 100 ng ml−1. IFNβ (R&D) at 10, 100 or 1,000 ng ml−1 or IFNγ (R&D) at 100 ng ml−1 were added to the culture. Cells were left in the RTCA analyzer until maximum cell death was reached (~40 h). The percentage of cytolysis was calculated with the measurement of the impedance, using the software provided with the RTCA technology.

Formalin-fixed, paraffin-embedded scRNA-seq and data analysis

Three lesional biopsies from three patients with DM were processed for scRNA-seq. For each sample, two formalin-fixed, paraffin-embedded curls (25 μm each) were dissociated using the Miltenyi Biotech FFPE Tissue Dissociation Kit (CG000632 RevA, 10× Genomics). The resulting cell suspension was divided equally into four centrifuge tubes, each containing approximately 300,000 cells. Cells in each tube were hybridized with a unique Probe Barcode, as per the instructions in the ‘Chromium Fixed RNA Profiling Reagent Kits for Multiplexed Samples’ user guide (CG000527, 10× Genomics). Post hybridization, cells from the four tubes were washed, counted and pooled in equal proportions. Approximately 40,000 cells from the pooled suspension were loaded onto a Chromium Q chip (PN-1000422, 10× Genomics). Sequencing libraries were prepared and sequenced on an Illumina NovaSeq platform using paired-end dual-indexing (28 cycles for read 1, ten cycles for i7, ten cycles for i5 and 90 cycles for read 2). The sequencing data were demultiplexed using bcl2fastq (Illumina). The resulting FASTQ files were processed with Cell Ranger using the multi-pipeline and the GRCh38-2020-A reference genome. Data analysis was performed using the Seurat package. Quality control and preprocessing steps included doublet removal using scDblFinder54 and ambient RNA correction using SoupX59. The data were then normalized, and variable features were identified. To correct for technical variation, we applied Harmony60 integration using the first 20 principal components, followed by UMAP. Clustering was conducted using the Leiden algorithm with a resolution parameter of 1.2. Cell type annotation was performed using established marker genes. For subsequent analyses, we focused specifically on the myeloid cell population.

Keratinocyte and fibroblast isolation and culture

Pannus samples were obtained under protocol H11731 (UMass Memorial Hospital) and stored in 70% ethanol, followed by sterile saline before processing. Four to five 4 mm punch biopsies were incubated in dispase I (25 mg ml−1; 04942086001, Roche) for 30 min and held overnight at 4 °C.

The following day, epidermis was separated from dermis and incubated in 0.25% trypsin-EDTA (25200056, Gibco) at 37 °C for 15 min with gentle agitation. Trypsin was neutralized with DMEM + 10% FBS to generate keratinocyte suspensions.

For fibroblast isolation, dermis was digested in collagenase (1 mg ml−1; 9001-12-1, Sigma-Aldrich) for 1 h at 37 °C, neutralized with DMEM + 10% FBS and dissociated by gentle pipetting. Cells were pelleted and cultured in DMEM with 10% FBS and 1% penicillin–streptomycin, with media changes every other day and passaging at sub-confluence (~5–7 days).

Immortalized N/TERT2G keratinocytes61 were generously provided by H. Rheinwald under a Material Transfer Agreement. Unless otherwise specified, N/TERT2G cells were used in an unprimed and undifferentiated state. N/TERT2G cells were cultured in Keratinocyte Serum-Free Media supplemented with epidermal growth factor and bovine pituitary extract (17005042, Gibco) and 1% penicillin–streptomycin (15140122, Gibco). Cultures were fed every other day and passaged upon reaching sub-confluence, approximately every 5–7 days.

Fibroblast stimulation

Fibroblasts were plated at 105 cells per well in a 24-well plate and incubated for 2–3 days at 37 °C. Cells were then stimulated for 24 h with either keratinocyte-conditioned medium (generated from NTERT keratinocytes treated and subsequently exposed to UVB), recombinant human IFNβ (50 ng ml−1; R&D Systems, 8499-IF-010/CF), recombinant human TNF (50 ng ml−1; BioLegend, 570102) or recombinant human IL-1α (100 ng ml−1; PeproTech, 200-01A-10UG). After 24 h, cells were lysed in RLT buffer for RNA library preparation.

Dendritic cell differentiation and stimulation

PBMCs from de-identified healthy donor leukopaks (UMass Leukocyte Core) were isolated by Ficoll-Paque PLUS centrifugation (17-1440-02, GE Healthcare), washed in PBS–EDTA and resuspended in RPMI 1640. Monocytes were enriched by CD14⁺ magnetic selection (EasySep Human CD14 Positive Selection Kit II; 17858, StemCell Technologies).

Monocytes were plated at 1 × 106 cells per ml in RPMI + 10% FBS, 1% penicillin–streptomycin, 2 mM L-glutamine, 50 ng ml−1 GM-CSF (572904, BioLegend) and 50 ng ml−1 IL-4 (574006, BioLegend). Medium was half-replenished with cytokine-supplemented RPMI on day 3. On day 7, moDCs were stimulated for 24 h with keratinocyte-conditioned media (from UVB-treated N/TERT2G keratinocytes), ultrapure LPS (100 ng ml−1; E.coli O111:B4, Invitrogen), recombinant IFNβ (50 ng ml−1; 8499-IF-010/CF, R&D Systems) or TNF (50 ng ml−1; 570102, BioLegend).

For selected studies, conditioned media were pre-treated with proteinase K (200 µg ml−1; RP107B-1, QIAGEN) at 55 °C for 1 h, then heat-inactivated DNase I (100 U ml−1, 37 °C, 30 min; 79254, QIAGEN) or passed through 0.22 µm filters (SLGVR33RS, MilliporeSigma) to assess the biochemical nature and size of active mediators. For blocking assays, moDCs were pre-incubated with adalimumab (20 µg ml−1, 1 h; HY-P9908, MedChemExpress) or anifrolumab (50 µg ml−1, 24 h; HY-P99168, MedChemExpress) before stimulation.

After incubation, cells were collected for RNA extraction or flow cytometry, and supernatants were collected for proteomic analysis. For T cell chemotaxis assays, moDCs were washed 12 h after exposure to keratinocyte supernatants and placed in fresh media; supernatants were collected 6 h later.

T cell chemotaxis with keratinocyte-conditioned DC supernatants

Fresh PBMCs from healthy donors were incubated overnight in chemotaxis medium (FBS-free RPMI + 1% ultrapure BSA, penicillin–streptomycin and L-glutamine) to enhance chemokine responsiveness. Non-adherent cells were then plated at 100,000 cells per well in the upper chamber of a 5 µm Boyden system and allowed to migrate for 2 h toward moDC-conditioned media. These moDC supernatants were generated after exposure to keratinocyte supernatants from UVB-irradiated ± IFNβ-primed N/TERT2G cells. Supernatants from untreated or LPS-stimulated moDCs served as negative and positive controls. Migrated cells were analyzed by flow cytometry (CD4, CD8, CD3, CD19, CD14) to identify responding T cell subsets.

In vitro UVB treatment

N/TERT2G keratinocytes (0.8–1 × 105 per well) were seeded in 12-well plates 3 days before irradiation. Cells were washed once with PBS and exposed to 50 or 100 mJ cm−2 UVB using a SolRx 100-Series handheld device (285–350 nm, Philips PL-S 9W/12/2P). PBS was then replaced with keratinocyte medium, and cells were incubated for 24 h at 37 °C with 5% CO2.

Membrane permeability assays

Cells were incubated with Hoechst 33342 (1:30,000; P3566, Invitrogen) and propidium iodide (1:1500; H3570, Invitrogen). Images were acquired every 2 h for 24 h using a Lionheart FX automated microscope (BioTek) and analyzed with Gen5 software. Membrane integrity was quantified as the ratio of propidium iodide-positive to Hoechst-positive cells.

ATP production assay

Keratinocyte ATP levels were quantified using the CellTiter-Glo 2.0 Assay (G9242, Promega). Keratinocytes were plated in 12-well dishes, treated as indicated and assayed 24 h later by adding an equal volume of CellTiter-Glo reagent. After a 10 min incubation at room temperature (22 °C (range 20–25 °C)) to ensure lysis and signal stabilization, luminescence was measured on a microplate reader. Values were normalized to untreated controls and used as a readout of metabolic activity.

Flow cytometry

The cells obtained through suction blister skin sampling were treated with Fc receptor blocking solution Human TruStain FcX (422302, Biolegend) and stained with the following dye and fluorescent antibodies according to the manufacturer’s instructions: Zombie NIR Fixable Viability Kit (423106, Biolegend), PerCP anti-human CD45 (368506, Biolegend; Clone 2D1, Monoclonal, Mouse IgG1, κ), Brilliant Violet 510 anti-human CD3 (344828, Biolegend; Clone SK7, Monoclonal, Mouse IgG1, κ), Spark NIR 685 anti-human CD19 (302270, Biolegend; Clone HIB19, Monoclonal, Mouse IgG1, κ), BUV737 anti-human CD56 (612766, BD Biosciences; Clone NCAM16.2, Monoclonal, Mouse BALB/c IgG2b, κ), Brilliant Violet 570 anti-human CD16 (302036, Biolegend; Clone 3G8, Monoclonal, Mouse IgG1, κ), Spark Blue 550 anti-human CD14 (367148, Biolegend; Clone 63D3, Monoclonal, Mouse IgG1, κ), PE/Fire 810 anti-human HLA-DR (307683, Biolegend; Clone L243, Monoclonal, Mouse IgG2a, κ) and eFluor 450 anti-human CD11c (48011642, Invitrogen eBioscience; Clone 3.9, Monoclonal, Mouse IgG1, κ). For flow cytometry analyses for T cell chemotaxis, Zombie Violet Fixable Viability Kit (423113, Biolegend), CD4 (344622, Biolegend; clone SK3, Monoclonal, Mouse IgG1, and R7-20041, Cytek; clone SK3, monoclonal, mouse IgG1, k), CD8 (344742, Biolegend; clone SK1, Monoclonal, Mouse IgG1, κ), CD3 (317332, Biolegend; clone OKT3, Monoclonal, Mouse IgG1, κ), CD19 (302206, Biolegend; Clone HIB19, monoclonal, mouse IgG1, k) and CD14 antibody (367148, Biolegend; Clone 63D3, Monoclonal, Mouse IgG1, κ) were used antibodies according to the manufacturer’s instructions to stain the migrated cells. For analysis of activated CD4 T cells, Zombie NIR viability dye, CD2 (300206, Biolegend; Clone RPA-2.10, Monoclonal, Mouse IgG1, κ), CD3 and CD4 antibodies, as mentioned above, were used for surface staining, and intracellular staining was performed using BD Cytofix/Cytoperm Fixation/Permeabilization Kit (554714, BD Biosciences) and Pacific Blue or PE/Dazzle 594 anti-human/mouse Granzyme B Recombinant Antibody (372216 and 372218, Biolegend; Clone QA16A02, Monoclonal, Mouse IgG1, κ) according to the manufacturer’s instructions. Cultured DCs were detached using ice-cold PBS, and prepared and stained with Fc receptor blocking solution, live dead zombie green (423111, Biolegend), CD11b (101226, Biolegend; Clone M1/70, Monoclonal, Rat IgG2b, κ), CD11c (301608, Biolegend; Clone 3.9, Monoclonal, Mouse IgG1, κ), CD14 (367148, Biolegend; Clone 63D3, Monoclonal, Mouse IgG1, κ), HLA-DR (562804, BD; Clone L243, Monoclonal, Mouse IgG2a, κ), CD83 (BD, 612823; Clone HB15e, Monoclonal, Mouse IgG1, κ), CD86 (305442, Biolegend; Clone IT2.2, Monoclonal, Mouse IgG2b, κ) and CD80 (563315, BD Biosciences; Clone L307.4, Monoclonal, Mouse IgG1, κ). An Aurora cytometer from Cytex Biosciences was used for sample analysis, and FlowJo Software (v.10.9.0) was used for data processing. All antibodies and reagents were tested and validated per the respective producers; information is available at biolegend.com, bdbiosciences.com, cytekbio.com and thermofisher.com.

Bulk RNA-seq

Total RNA was extracted from frozen samples in RLT buffer using the RNeasy Mini Plus Kit (QIAGEN). To eliminate residual genomic DNA, the RNA was treated with RNase-free DNase I for 15 min at room temperature. RNA-seq libraries were constructed using an optimized CellSeq2 protocol for bulk RNA libraries62. Each library was prepared by pooling 24 samples, with 50 ng of input RNA per sample. Library construction primers are provided in Supplementary Data 1.

The RNA libraries were sequenced on the NovaSeq 6000 platform. Raw FASTQ files were processed using the following pipeline, https://github.com/shuo-shan/viafoundry-celseq2-bulk-pipeline, with a modified GRCh38 genome as the reference (details in Supplementary Information), to generate the final count matrices that were then analyzed in R (v.4.3.3). When required, batch effects owing to sequencing run and inter-donor variability were adjusted using the ComBat-seq63 function from the sva package. Differential expression analysis was performed using DESeq2 (ref. 49), with a design formula in which stimulation was encoded as the condition factor. Genes were considered differentially expressed if they met all the following criteria: mean normalized expression (baseMean) ≥ 10, P value < 0.05 and absolute log2(fold change) ≥ 1.

Spatial genomics gene positioning system (GenePS)

A custom seqFISH panel containing 521 genes, including 21 genes to be imaged one-at-a-time, was designed with Spatial Genomics (Primary Probe Kit, v500; no. 30021131). Fresh-frozen human skin biopsy samples obtained at UMass Memorial Hospital under protocol H00021295 were cryo-sectioned 10 µm thick and mounted on specialized coverslips (10200003, Spatial Genomics). Immediately after collection, the sections were fixed with freshly prepared 4% paraformaldehyde (28908, Thermo Scientific) in 1× PBS for 15 min at room temperature. Following fixation, the sections were washed three times with 1× PBS for 5 min each and dehydrated in 70% ethanol for 30 s at room temperature. After air drying, the sections were stored at −80 °C until further preparation for seqFISH experiments. Tissue processing was performed in accordance with Spatial Genomics’ Sample Processing for RNA Profiling Guide (ASY-053).

Tissue sections were imaged using the Gene Positioning System (GenePS; Spatial Genomics), an automated platform for image acquisition, reagent delivery and data processing. The system was loaded with the sample, decode plate and buffers following the GenePS Instrument User Guide (INS-205) and guided by on-screen instructions. Regions of interest within human skin sections were selected using DAPI-stained nuclei imaged with a ×10 air objective. Once regions of interest were selected, the automated imaging workflow was initiated, using a ×60 oil objective with standard image acquisition settings and capturing either three or six z-planes at a z-step size of 1.5 µm. The experimental protocol included iterative rounds of decode probe hybridization, imaging and signal removal, continuing until all hybridization rounds were completed.

Raw image files were aligned across hybridization rounds and processed directly on the GenePS instrument to localize RNA signals. Subsequent analysis using Spatial Genomics’ Navigator software enabled transcript decoding and cell segmentation. The segmentation algorithm used machine learning to delineate nuclei based on DAPI staining, with segmentation boundaries expanded by 15 pixels (~1.5 µm) past the nuclei borders. RNA transcripts were decoded with a maximum allowed FDR of 5%, and the transcripts were assigned to individual cells, generating cell-by-gene count matrices for each region of interest. The processed outputs included a cell-by-gene matrix, cell coordinates, transcript list, DAPI image and segmentation mask, providing high-resolution spatial and transcriptional data for downstream analysis.

seqFISH+ analysis

Cell-by-gene matrices, cell coordinates and transcript locations were used to create the Seurat47 object. Cells with over 50 transcripts or over ten transcripts at the skin surface were retained. IntegrateLayers – RPCA was used for integration on the top 30 PCs, followed by UMAP and clustering (FindNeighbors/FindClusters). Cells were first grouped into major lineages, followed by subtype reclustering. Doublet clusters, identified by inappropriate marker co-expression, were removed before final clustering of lymphoid and myeloid compartments.

To order MC_CD14 cells along the MMP9 transition, the top 100 variable features were selected by FindVariableFeatures and then used for PCA. The top five PCs were used for CCAIntegration and the following DiffusionMap56. The first diffusion component defined the transition. For heatmaps, MC_CD14 cells were divided into 30 equal bins, and mean expression per bin was calculated. Expression values were scaled per gene across all bins and reference cell types.

Spatial images were generated using software Spatial Genomics Navigator Lite.

Differential expression analysis was performed using the Model-based Analysis of Single-cell Transcriptomics (MAST) framework64, with default parameters wrapped by Seurat. Differential analysis along the CD14+ myeloid cell transition axis was performed using tradeSeq65 with cellWeights = 1.

Colocalization and visualization functions are packed in the scSpatial R package available on GitHub (https://github.com/garber-lab/scSpatial). The colocalization strength of each query cell within a reference cell type was calculated using the createField function (with s.d. = 400 pixels) and the getValueInField with default parameters. The summation of colocalization strengths of all query cells was used to represent the colocalization score between two cell types. For the heatmap, colocalization scores were scaled so that the total score of all cells in the sample was zero, and the score of the query cell type itself was set to 100. Heatmap visualization and clustering were performed using the ComplexHeatmap R package52,66. See the Supplementary Information and Methods for more details on colocalization analysis. Gene expression trends along colocalization scores or distances were assessed using the FeaturePlot.cont.rollmean function, with window.prop = 1.

The chemokine expression levels around a given CD14+ cell were summed within a radius of 400 pixels (41.2 μm). Local expression of chemokines as a function of CD14+ cells transition score was visualized with geom_smooth (method = “loess”, span = 0.2).

The smoothed expression level of TNF and IFNG over spatial embedding was visualized with ImageFeaturePlot.contour (sigma = 400, threshold = 0.9, n_levels = 6).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.