We selected diseases based on the availability of public scRNA-seq datasets across skin, lung, colon and kidney healthy control tissue, and their respective diseases (encompassing AD, PSO, COPD, IPF, UC, IgAN, LN) and at least two samples per group. We collected and analyzed the following 15 scRNA-seq datasets in our analysis: 1) E-MTAB-8142 (Reynolds et al., 2021) for healthy skin, AD lesional skin, and AD nonlesional skin, PSO lesional skin and PSO nonlesional skin, 2) GSE147424 (He et al., 2020) for healthy skin, AD lesional skin, and AD nonlesional skin, 3) GSE153760 (Rindler et al., 2021) for healthy skin and AD lesional skin, 4) GSE173706 (Merleev et al., 2022) for healthy skin, PSO lesional skin, and PSO nonlesional skin, 5) GSE220116 (Kim et al., 2023) for healthy skin and PSO lesional skin, 6) EGAS00001004344 (Travaglini et al., 2020) for healthy lung, 7) GSE136831 (Adams et al., 2020) for healthy lung, COPD lung and IPF lung, 8) GSE171541 (Huang et al., 2022) for healthy lung and IPF lung, 9) GSE122960 (Reyfman et al., 2018) and 10) GSE135893 (Habermann et al., 2020) for healthy lung and IPF lung, 11) SCP259 (Smillie et al., 2019) and 12) GSE116222 (Parikh et al., 2019), and 13) GSE231993 (Du et al., 2023) for healthy colon, UC inflamed colon, and UC uninflamed colon, 14) GSA: HRA000342 (Zheng et al., 2020) for healthy kidney and IgAN kidney, and 15) LN_Kidney_AMP (from the Accelerating Medicines Partnership (AMP) SLE phase 2 consortium) for healthy kidney and LN kidney (Arazi et al., 2019; Hoover et al., 2023; Horisberger et al., 2024; Izmirly et al., 2024).

We processed the publicly deposited count matrices for these datasets using CellBridge [version 1.0.0 (Nouri et al., 2023)], an automated, Docker-based scRNA-seq data processing pipeline developed by our research group, with empirical filtering parameters for each dataset (specified in Supplementary Table 1) and specific pipeline steps described hereafter. Doublet removal was performed using the Scrublet Python package [version 0.2.3; (Wolock et al., 2019)]. The Seurat R package (version 5.0.1; (Hao et al., 2024)) was used to perform initial processing (“NormalizeData,” “FindVariableFeatures,” “ScaleData,” “RunPCA” with default parameters). To improve our cell type annotation (which relies on dimensionality reduction for smoothing), sample-based batch correction of the PCA was performed using the Harmony R package (version 1.0; (Korsunsky et al., 2019)), followed by the nearest-neighbor graph construction using “FindNeighbors” with 30 dimensions (with default settings except “dims = 1:30”), subsequent cluster identification with “FindClusters” (with default settings except “resolution = 0.7”), and the Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction with “RunUMAP” on 30 dimensions (with default settings except “dims = 1:30”). Cell type annotation was performed by the SignacX R package [version 2.2.5; (Chamberlain et al., 2023)], also produced from our group, that is a classifier designed to predict cellular phenotypes in scRNA-seq data from multiple tissues. Cell classifications from SignacX include: memory B cells (“B.memory” in SignacX nomenclature), naïve B cells (“B.naive”), dendritic cells (“DC”), endothelial, epithelial, fibroblasts, macrophages, classical monocytes (“Mon.Classical”), non-classical monocytes (“Mon.NonClassical”), neutrophils, natural killer cells (“NK”), plasma cells (“Plasma.cells”), CD4⁺ naïve T cells (“T.CD4.naive”), CD4⁺ memory T cells (“T.CD4.memory”), CD8⁺ naïve T cells (“T.CD8.naive”), central memory CD8⁺ T cells (“T.CD8.cm”), CD8⁺ effector memory T cells (“T.CD8.em”), and regulatory T cells (“T.regs”).

To identify potential ligand-receptor interactions taking place between cells, we used CellphoneDB version 5.0, which contains a curated repository of ligand–receptor interactions and a permutation-based statistical framework for inferring enriched interactions between cell types from scRNA-seq data (Efremova et al., 2020). Count matrices and metadata files extracted from the Seurat RDS provided by CellBridge were processed using the default parameters of the statistical method of CellphoneDB v.5.0 (https://github.com/ventolab/CellphoneDB) in Python (version 3.7; Python Software Foundation). We filtered the CellphoneDB statistically significant means for interactions where either the ligand or the receptor is a chemokine gene (i.e., HGNC gene symbol starting with CC or CX) and from a list of genes known to be involved in immune cell extravasation (listed in Supplementary Table 2; (DeGrendele et al., 1997; Johnson-Léger et al., 2002; Ludwig et al., 2005; Farkas et al., 2006; Wegmann et al., 2006; Engelhardt, 2008; Kumar et al., 2012; Azcutia et al., 2013; Rocha et al., 2014; Baeyens et al., 2015; Glatigny et al., 2015; Zhang et al., 2016; Bros et al., 2019; Bui et al., 2020; Arif et al., 2021; Amersfoort et al., 2022; Czubak-Prowizor et al., 2022; Ma et al., 2023; Xu et al., 2023)). For the extravasation-based output, the sender cell (i.e., the cell expressing the ligand) was filtered to Endothelial cells, as this is the cell type involved in mediating immune cell extravasation through blood vessel endothelium (Luster et al., 2005).

We conducted pairwise (disease vs. healthy) comparisons of cell type proportions using the Brunner-Munzel test in the brunnermunzel R package (version 2.0) on each dataset. Analysis was done with the minimum number of subjects per group to hold a comparison (per cell state) set to two, and the minimum number of total cells per subject set to 200. Results with Benjamini-Hochberg-corrected p-values <0.05 were considered statistically significant.

We conducted pairwise (disease vs. healthy) DEG analysis using the zero-inflated negative binomial model (zlm) in the model-based analysis of single-cell transcriptomics (MAST) R package [version 1.26.0; (Finak et al., 2015)] on each cell type in each dataset with sample covaried. For a gene to be considered, it had to be expressed in at least 25% of cells in a given cell type and have an absolute log2-fold change of at least 0.25. Genes with Benjamini-Hochberg-corrected false discovery rate (FDR) of <0.05 were further filtered for our target genes, which included chemokine genes (i.e., HGNC gene symbol starting with CC or CX), immune cell extravasation, proliferation, and activation genes (listed in Supplementary Table 2 as indicated earlier).

After conducting each analysis (cell-cell communication, cell proportion, DEG) per dataset, RStudio (RStudio Team, version 4.2) was used to extract common statistically significant findings across datasets, within each disease group in the case of cell proportion, DEG and cell-cell communication analyses, as well as the corresponding healthy groups (healthy skin, healthy lung, healthy colon, healthy kidney) in the case of cell-cell communication analyses (Figure 1).

For cell-cell communication, interactions that were statistically significant (CellphoneDB’s permutation test, p < 0.05) were selected for further analysis. We analyzed whether an interaction was multi-tissue (i.e., across multiple (healthy or disease) tissues), or disease-enriched (i.e., specific to diseased tissues only), which could then be further subdivided according to specificity, from occurring in multiple diseases across tissues, to occurring only in diseases of specific tissues (e.g., skin), to those only occurring in a particular disease. Visualizations were performed in RStudio (RStudio Team, version 4.2) using the ggplot2, ggsankey, and circlize libraries (Gu et al., 2014).

To validate potential interactions that arose from our cell-cell communication analyses, we chose two findings that had been unexplored in the literature in terms of chemotaxis assays: chemokine (C-X-C motif) ligand 2 (CXCL2) signaling (from various immune cells) to T.CD8.em cells (via dipeptidyl peptidase-4 (DPP4)) and retinoic acid receptor responder 2 (RARRES2) signaling to classical and nonclassical monocytes. First, we measured chemotaxis of pan CD8⁺ T cells (Charles River, cat. PB08NC-2) and pan monocytes (Charles River, cat. PB14-16NC-1) from 2-3 donors towards CXCL2 (R&D Systems, cat. 276-GB-010) and RARRES2 (R&D Systems, cat. 224-CM-025), respectively, using a transmigration chamber assay with 5 μm pore size transwell inserts (Corning, cat. CLS3421) according to the manufacturer’s instructions. Briefly, the transwell inserts were placed into 24-well plates containing 500 ul of Roswell Park Memorial Institute (RPMI) 1,640 media (Gibco, cat. 22400-089) with 1% heat-inactivated fetal bovine serum (FBS; Gibco, cat. A38400-01, lot. 2717806RP) alone (the medium only negative control condition), medium with a positive control, or medium with each target ligand (at multiple dilutions). 2–4 wells were used per condition. For the CD8⁺ T cell experiments, chemokine (C-X-C motif) ligand 12 (CXCL12, 10 ng/mL; R&D Systems, cat. 350-NS-010) was used as the positive control. For the monocyte experiments, two positive controls were used; 5% (v/v) cobra venom activated human complement serum (CAS; Complement Technology Inc, cat. NC1769554), as well as CC motif chemokine ligand 8 (CCL8; R&D Systems, cat. 281-CP-010). Cells were quickly thawed in a 37°C water bath, spun down at 300 g for 10 min and allowed to rest for 1 h in an incubator (37°C, 5% CO2) in RPMI 1640 with 10% FBS, prior to being spun down at 300 g for 10 min, counted with the CellacaMX (Nexcelom Bioscience), and then resuspended in 1% FBS media and loaded at a 100 ul volume, at 200,000 cells/transwell insert for the CD8⁺ T cell experiments and 250,000 cells/transwell insert for the monocyte experiments. Transwell plates were then incubated for 3 h to permit chemotaxis. Cell migration was quantified using CellTiter-Glo (Promega, cat. G9241) relative luminescence (RLU) measured with a plate reader (EnVision, 2104-0010). To find significant changes in cell migration of each concentration of tested ligand compared to the medium control, student t-tests were performed on R studio (version 4.2). For each donor, migration index (MI) was calculated using the following formula and used for plotting. M I w e l l = log R L U w e l l   R L U ¯ m e d i u m   n e g a t i v e   c o n t r o l  

After conducting the initial chemotaxis experiments, we investigated whether knock-out of the receptor DPP4, identified through cell-cell communication analyses in mediating CXCL2 binding in CD8⁺ effector memory T cells, inhibits chemotaxis, and whether the expression of DPP4 changes post-chemotaxis. Firstly, on day 0, a subset of CD8⁺ T memory cells (cluster of differentiation 45 restricted (CD45RO+); STEMCELL Technologies, cat. 200-0168) from two donors underwent the chemotaxis assay as described in the previous section. Transwell-migrated cell solutions (two sets of two wells pooled, per condition) were used for flow cytometry analysis for the receptor genes of interest. These cells along with compensation beads were placed in a 96-well clear round bottom plate (Corning, cat. 3,799), washed twice by spinning down at 350 g for 5 min and resuspended with 1X PBS (Gibco, cat. 10010049), stained in the dark for 15 min at room temperature (RT) with Zombie NIR Fixable viability dye (Biolegend, cat. 423105), washed twice with 1X PBS, stained in the dark for 30 min at 4°C with an extracellular panel of antibodies at 1:100 dilution (Pacific blue CD8 (cat. 344717), BV570 CD45RO (cat. 304225), BV785 CD25 (cat. 302637), FITC DPP4 (cat. 302638), PE CXC receptor 4 (CXCR4, cat. 306505), PE/Dazzle 594 CD62 L-selectin (CD62L; cat. 304841), APC CXCR2 (cat. 320710); all anti-human and Biolegend) with Fc block (Biolegend, cat. 422302), washed twice with cell staining buffer (Biolegend, cat. 420201), permeabilized and fixed (using Fixation/Permeabilization Concentrate, ThermoFisher, cat. 00-5123-43 and Fixation/Permeabilization Diluent, ThermoFisher, cat. 00-5223-56) for 30 min at RT, washed twice with perm/wash buffer (ThermoFisher, cat. 00-8333-56), stained in the dark again with the CXCR4, CXCR2 and DPP4 antibodies and Fc block for intracellular staining, washed twice with cell staining buffer, resuspended in 40 ul cell staining buffer, and quantified with the iQue3 (Sartorius) at a sip time of 37 s. Analysis was conducted on OMIQ (Dotmatics); unfiltered data was gated for cells, followed by live/dead gating, singlets, and then fluorescence minus one (FMO)-based gating of CD8⁺ cells, followed by the target proteins (CXCR4, CXCR2, DPP4, CD25). Percentage of cells positive for each marker (receptors including CXCR4, CXCR2, DPP4 and activation marker CD25) and their fluorescence intensity was compared to the medium negative control with student t-tests in R Studio (version 4.2).

The remaining CD8⁺ memory T cells not used for the chemotaxis assay were incubated in a T25 flask with ImmunoCult-XF T cell expansion medium (ICEM; STEMCELL Technologies, cat. 10981) and 50 IU/mL of IL-2 (STEMCELL Technologies, cat. 78220) at 1 million cells/mL. On day 1, cells were activated with 25 ul/mL CD3/28 activator (STEMCELL technologies, cat. 10971). On day 4, cells were divided amongst the following knock-out groups: wild-type (WT), “electroporation-only,” DPP4 knock-out (KO), and non-target guide control (NTC). The DPP4 KO and NTC groups were washed with 1X PBS and resuspended with P3 primary nucleofection solution (Lonza, cat. V4XP3012) containing ribonucleoprotein (RNP) complexes of 20 pmol cas9 nuclease protein (Horizon Discovery, cat. CAS12205) and 200 pmol guide RNA (gRNA) per 1 million cells, targeting either DPP4 (Horizon Discovery, cat. SQ-004181-01-002) or the NTC locus (Synthego, CRISPRevolution sgRNA EZ kit). The RNPs were assembled by 37°C. incubation for 10 min, followed by RT incubation for 5 min. The cell-RNP mixture was nucleofected using an electroporator (4D-Nucleofector X Unit, Lonza) according to the manufacturer’s instructions, and then incubated with ICEM in a 12-well plate. The “electroporation-only” group underwent the same procedure, but without any introductions to RNPs, and the WT group was simply washed with 1X PBS, resuspended in ICEM and transferred to the 12-well plate. On day 11, 100,000 cells per group and per donor underwent the flow cytometry protocol outlined earlier to verify whether DPP4 was KO’d in the DPP4 KO group and still expressed amongst the control groups (WT, electroporation-only and NTC). On day 12, WT, DPP4 KO and NTC groups were used in the CXCL2 chemotaxis assay (with medium and CXCL12 controls) for chemotaxis and flow cytometry quantification as described earlier.

We curated 15 datasets representing the following: 1) three datasets of AD skin (lesional in all and non-lesional in two of them), 2) three datasets of PSO skin (lesional in all and non-lesional in two of them), 3) two datasets of IPF lung (one of which also contained COPD) and another lung dataset containing COPD, 4) three datasets of UC colon (inflamed and non-inflamed), and 5) two kidney datasets, one of IgAN, and another one of LN. All datasets contained healthy controls as well (Figure 1A).

After conducting each form of analysis (cell-cell communication via CellphoneDB, cell proportion, DEGs) per dataset, we extracted findings that were present in all datasets of a particular disease group (AD nonlesional skin, AD lesional skin, PSO nonlesional skin, PSO lesional skin, COPD lung, IPF lung, UC uninflamed colon, UC inflamed colon, IgAN kidney, LN kidney) in the case of cell-proportion, DEG and cell-cell communication analyses, as well as healthy group (healthy skin, healthy lung, healthy colon, healthy kidney) in the case of cell-cell communication analyses (Figure 1B). Statistical results for the cell proportion data are reported in Supplementary Table 3 and for the DEG data reported in Supplementary Table 4.

For cell-cell communication analysis, after extracting findings that were present in all datasets of a particular disease group and healthy group, the next process involved four steps: 1) comparing CellphoneDB findings from each disease tissue with those from all healthy tissues to extract ‘disease-enriched’ interactions, 2) further comparing these ‘disease-enriched’ findings against each other to extract disease-specific interactions (e.g., interactions only present in COPD lung but in no other healthy or disease group), 3) comparing findings from each healthy tissue against those from every other healthy tissue to extract healthy-tissue-specific interactions (e.g., interactions only present in healthy lung, but not in healthy skin, colon or kidney), and 4) extracting disease interactions that overlapped with healthy tissue interactions. These CellphoneDB cell-cell communications and their interaction values [the mean of the gene expression level for cells in a pair (Efremova et al., 2020)] found in healthy tissues are reported in Supplementary Table 5, interactions that are disease-enriched are reported in Supplementary Table 6, and disease tissue interactions that also occurred in healthy tissue are reported in Supplementary Table 7.

Overall, we identified 39 unique chemokine-based L-R pairs by cell-cell communication analysis (Figure 2A). The majority of these interactions (28 or 72%) were “disease-enriched,” i.e., they did not occur in any healthy tissue. Of these disease-enriched L-R pairs, 11 (28%) were specific to a particular disease (referred to as “disease-specific”), 12 (31%) were specific to a particular organ but co-occurring in multiple diseases of that organ (referred to as “disease organ-specific”), five (13%) occurred in various diseases (referred to as “disease multi-tissue”). Interestingly, among the disease multi-tissue L-R pairs, the sender and receiver cell types implicated in their expression were entirely distinct depending on the diseased organ type. For example, the CCL3L1 to CCR1 interaction was found in IgAN kidney, LN kidney and UC uninflamed colon; however, in UC uninflamed colon, the CCL3L1 to CCR1 interaction occurred between Neutrophils to Macrophages, whereas in IgAN or LN kidney, other cell types were involved in this communication.

Unlike the disease-enriched interactions, there were no L-R pairs that were “healthy-enriched,” i.e., only found in healthy tissues; of the remaining 11 (28%) interactions, all occurred in both disease and healthy tissue. Of those, six (16%) were specific to a particular healthy tissue, while also occurring in various diseases (referred to as “healthy tissue-specific and disease multi-tissue”), three (8%) were occurring in varieties of both healthy and disease tissues (referred to as “healthy and disease multi-tissue”), and two (5%) were occurring in various healthy tissues and ubiquitously occurring in all the diseases examined in the study (referred to as “healthy multi-tissue and disease ubiquitous”).

There were 30 unique extravasation-based L-R pairs identified by cell-cell communication analysis (Figure 2B). Unlike among the chemokine-based L-R pairs, a minority (33%) of extravasation-based L-R pairs were “disease-enriched’; of these interactions, nine (30%) were disease-specific, one (3%) was “disease organ-specific,” and one (3%) was in the “disease multi-tissue” set. The only “disease multi-tissue” L-R pair was SELE to GLG1, which occurred in inflamed UC colon and nonlesional AD and PSO skin.

As with the chemokine-based L-R pairs, none of the extravasation-based L-R pairs were found exclusively in healthy tissue. However, contrary to the chemokine-based L-R pairs, the majority of the extravasation-based interactions (20; 67%) occurred in both healthy and disease tissue. Of the 20, seven (23%) were specific to a particular healthy tissue and specific to a diseased organ (referred to as “healthy tissue-specific and disease organ-specific”; interestingly, all of them were specific to the kidney), six (20%) were specific to a particular healthy tissue, while also occurring in various diseases (referred to as “healthy-tissue specific and disease multi-tissue”), and six (20%) were occurring in varieties of both healthy and disease tissues (referred to as “healthy and disease multi-tissue”). Of note, none of the extravasation-based L-R pairs were healthy multi-tissue and disease ubiquitous, as found in the chemokine-based L-R pairs.

Taking into account all four components of cell-cell communication—the sender cell, its ligand, the receiver cell, and its receptor—we proceeded to examine interactions specific to each of the healthy tissues studied (skin, lung, colon, and kidney) and their corresponding diseases. Intriguingly, in each of the healthy tissues (lung, skin, colon and kidney), sender-receiver pairs were unique to each tissue type regardless of the chemokine-based ligands and receptors. For example, in healthy lung, there was a Macrophage to Endothelial cell interaction, and this cell pair was not found in any other tissue group. To illustrate these findings, we used bolded chords for all the healthy tissue chemokine interactions in healthy skin (Figure 3A), healthy lung (Figure 3C), healthy colon (Figure 3E) and healthy kidney (Figure 3G). In the disease datasets, we observed both tissue-exclusive and non-exclusive sender-receiver pairs. Tissue specificities for sender-receiver pairs involved in cell-cell communication are also indicated in Supplementary Table 5 (for the healthy tissues) and Supplementary Table 6 (for the disease-tissues). The subsequent sections will elaborate on the tissue-specific interactions in more detail.

In the IgAN kidney dataset, we found cell-cell communications which can be considered atypical chemokine ligand to receptor bindings that have not yet been investigated in chemotaxis assays based on our review of the literature, to the best of our knowledge. This included: 1) CXCL2 signaling, known as a neutrophil chemoattractant (Li et al., 2016), to T.CD8.em cells via DPP4, also known as CD26 and 2) RARRES2, an adipokine (Tang et al., 2023), signaling to classical and nonclassical monocytes via chemokine receptor CCRL2. RARRES2 has been gaining attention in recent years due to its roles in metabolic and inflammatory processes (Kwiecien et al., 2020; Ma et al., 2024). First, we measured chemotaxis of pan CD8⁺ T cells and pan monocytes towards their respective target chemokines in a transwell assay and found that CXCL2 did indeed chemoattract CD8⁺ T cells; specifically, 1 ng/mL of CXCL2 had a significantly higher migration index compared to the medium-only control (p < 0.05 for both donors, Supplementary Figure 2A). CXCL12, a chemokine well-known to chemoattract T cells (Hara et al., 2006), also led to significantly higher migration (p < 0.001 for both donors). However, RARRES2 did not chemoattract monocytes (Supplementary Figure 2B), unlike CCL8 [also known as monocyte chemoattractant protein-2 (MCP-2)] which led to significantly higher migration for both donors, particularly at the 0.1 ng/mL concentration (p < 0.01, Supplementary Figure 2C), compared to the negative controls. The data and t-test results for the transwell assays are reported in Supplementary Tables 8, 9, respectively.

After establishing which of the tested ligands chemoattracted the predicted cell type(s), we asked two follow-up questions. First, we aimed to determine whether receptor expressions would change in the presence of the chemoattracting ligand, and specifically, a) if the expression of DPP4 and CXCR2, the canonical receptor of CXCL2 (Zhang et al., 2017), increases in the presence of CXCL2 post-chemotaxis and b) if CXCR4, the canonical receptor of CXCL12 (Shi et al., 2020), increases when CXCL12 is expressed, and c) if the expression of CD25, an activation marker (Bajnok et al., 2017), increases post-chemotaxis. Second, we tested the involvement of DPP4 in the CXCL2-induced chemoattraction of CD8⁺ effector memory T cells by knocking out DPP4.

To answer the first question, we conducted the CXCL2 transwell assay with the CD8⁺ effector memory T cells and then collected the migrated cells for flow cytometry quantification (Figure 5A). Ultimately, there were no significant changes in the levels of expression of DPP4, CXCR2, CXCR4, or CD25. However, we did find that memory CD45RO+ CD8⁺ T cells, most of which were effector memory cells (CD62L-CD45⁺; Figure 5B), migrated towards CXCL2 (Figure 5C). Like the pan CD8⁺ T cell findings, 1 ng/mL CXCL2 led to significantly higher migration compared to the medium-only condition (p < 0.05 for both donors). Next, we assessed if knocking out DPP4 would affect the ability of the CD8⁺ memory T cells to migrate to CXCL2 (Figures 5D–F). We determined that DPP4 KO cell’s CXCL2-induced migration was not significantly different from the medium-only condition, contrasting with wild-type (WT) and NTC cells in donor 1 and 2 respectively, that showed significantly higher migration towards CXCL2 compared to the medium-only control (p < 0.05). In contrast, the CXCL12 positive control showed similar levels of significantly higher migration compared to the medium-only control regardless of treatment group (p < 0.001 for WT, KO, NTC for both donors) Figure 5G. These findings suggest that the CXCL2-DPP4 axis might indeed play a role in the chemoattraction of CD8⁺ effector memory T cells.