Research ArticleFIBROSIS

Interleukin-36γ–producing macrophages drive IL-17–mediated fibrosis

See allHide authors and affiliations

Science Immunology  11 Oct 2019:
Vol. 4, Issue 40, eaax4783
DOI: 10.1126/sciimmunol.aax4783

Sequencing in the matrix

Biological scaffolds that mimic tissue microenvironments can be used to model wound healing. Using two distinct biomaterial environments, one that promotes regeneration and another that promotes fibrosis in conjunction with single-cell RNA sequencing, Sommerfeld et al. have characterized macrophages involved in both fibrosis and regeneration. They have defined populations of macrophages that are involved in both processes, and they have identified a subset of CD9+ interleukin-36γ (IL-36γ)–producing macrophages that participate in IL-17–driven fibrosis. By evaluating wound healing in IL-17–deficient mice, they report IL-17 to be essential for the generation of CD9+ IL-36γ–producing macrophages during fibrosis. Further studies are needed to understand the functional relationship between IL-17– and IL-36γ–producing cells in fibrosis and in other settings.

Abstract

Biomaterials induce an immune response and mobilization of macrophages, yet identification and phenotypic characterization of functional macrophage subsets in vivo remain limited. We performed single-cell RNA sequencing analysis on macrophages sorted from either a biologic matrix [urinary bladder matrix (UBM)] or synthetic biomaterial [polycaprolactone (PCL)]. Implantation of UBM promotes tissue repair through generation of a tissue environment characterized by a T helper 2 (TH2)/interleukin (IL)–4 immune profile, whereas PCL induces a standard foreign body response characterized by TH17/IL-17 and fibrosis. Unbiased clustering and pseudotime analysis revealed distinct macrophage subsets responsible for antigen presentation, chemoattraction, and phagocytosis, as well as a small population with expression profiles of both dendritic cells and skeletal muscle after UBM implantation. In the PCL tissue environment, we identified a CD9hi+IL-36γ+ macrophage subset that expressed TH17-associated molecules. These macrophages were virtually absent in mice lacking the IL-17 receptor, suggesting that they might be involved in IL-17–dependent immune and autoimmune responses. Identification and comparison of the unique phenotypical and functional macrophage subsets in mouse and human tissue samples suggest broad relevance of the new classification. These distinct macrophage subsets demonstrate previously unrecognized myeloid phenotypes involved in different tissue responses and provide targets for potential therapeutic modulation in tissue repair and pathology.

INTRODUCTION

Macrophages play successive roles in the tissue response to trauma and foreign bodies from initial inflammation to innate defense and resolution that result in tissue repair or chronic inflammation and fibrosis depending on environmental signals (1). Biomaterials, medical devices, and prostheses induce a foreign body response (FBR) that is composed, in part, of macrophages and foreign body giant cells. The pathological consequence of the FBR is a capsule of fibrotic tissue surrounding the material. Development of biomaterial scaffolds for regenerative medicine that are designed to promote tissue repair through incorporation of biological cues is shifting and sometimes eliminating the FBR through changes in the macrophage response and phenotype.

Macrophage functional heterogeneity is key to their ability to respond to diverse environments and cues. Yet this complexity in activity and gene expression is not reflected by current phenotyping and nomenclature dogma. Early attempts to classify macrophages resulted in the M1/M2 nomenclature to define the proinflammatory, interferon γ (IFNγ)–activated macrophage (M1) versus the “alternatively IL(Interleukin)-4 activated” macrophage (M2) (2, 3). After this, a spectrum model of macrophage activation was proposed (4). In recent years, a common framework for macrophage activation has been suggested, which considers a set of surface and genetic markers including CD206 as an “M2 marker” and CD86 as an “M1 marker” (5). Recent technological advances, however, including single-cell RNA sequencing (scRNAseq) and mass cytometry suggest that archetypal in vitro phenotypes rarely overlap with those found in physiological conditions (6, 7). These approaches open the door to the refined characterization of macrophage populations that reflects the complexity of their phenotype and function in normal and pathological environments.

Biomaterials generate tissue microenvironments that can reproducibly induce specific macrophage phenotypes. Biological scaffolds, derived from the extracellular matrix (ECM) of tissues, promote a pro-regenerative tissue microenvironment that induces tissue repair through increased expression of IL-4. The repair capacity of these materials correlates with a T helper 2 (TH2) T cell response that directs polarization to a traditionally defined M2 macrophage in combination with a reduction of CD86 expression (8, 9). Synthetic materials induce an FBR that has been associated with inflammatory M1-type macrophages and development of a fibrotic capsule (10, 11). Although the M1-type inflammatory macrophage has been conventionally associated with a TH1 response, we found that the FBR and the associated macrophage phenotype occur in a type 17 immune environment that includes IL-17 production by innate lymphocytes, γδ T cells, and TH17 T cells. IL-17 signaling is required for the fibrosis associated with the FBR, although the level of IL-17 expression may vary depending on the chemical and physical properties of the materials (12). Synthetic and biological materials, therefore, serve as a model for type 2 and type 17 tissue immune microenvironments rich in IL-4 and IL-17, where the associated macrophage phenotype can be studied.

Here, we used scRNAseq to characterize macrophages isolated from a murine model implanted with biomaterials that model tissue repair versus fibrosis. We identified inflammatory macrophages from regenerative and fibrotic environments with varying inflammatory signatures. We also defined phenotypically and functionally a phagocytic macrophage population exclusive to the reparative environments. Last, we identified a CD9hi+IL-36γ+ fibrotic macrophage population in the FBR response with signatures of type 17 immune responses and autoimmunity. We also identified surface markers that defined the distinct macrophage clusters and validated their ability to separate macrophage subsets experimentally by flow cytometry and immunofluorescence. Last, we validated the broader relevance of the macrophage subsets in human pathologies including histiocytosis and fibrosis.

RESULTS

Single-cell profiling of F4/80hi+ macrophages in regenerative and fibrotic tissue microenvironments

To model tissue microenvironments with activated and heterogeneous macrophage populations, we selected biomaterials that are used clinically and induce divergent immune and tissue responses. Biological scaffolds derived from the ECM of tissues promote healing of muscle (13) and induce a type 2 (TH2/IL-4) immune microenvironment (9). We selected a urinary bladder matrix (UBM) that is used clinically in wound healing (14) and hernia repair (15). Polycaprolactone (PCL) was selected as a model synthetic biomaterial that stimulates a fibrotic response (16) and a type 17 (TH17/IL-17) immune microenvironment (12). We sorted differentiated macrophages (CD45+CD64+F4/80hi+) (17) for single-cell analysis 1 week after biomaterial implantation and applied them to the 10× scRNAseq platform (Fig. 1A). We selected this subpopulation of macrophages for analysis because they are used to analyze CD86 (a ligand for CD28 and CTLA4) and CD206 (a scavenger receptor) expression to determine M1 and M2 polarization and infer phenotype.

Fig. 1 Single-cell characterization of macrophages in fibrotic and regenerative microenvironments.

(A) Experimental overview. A virtual aggregate of macrophages in fibrosis and regeneration generated from scRNAseq after sorting of F4/80hi+CD64+. Cells were isolated from murine volumetric muscle injuries at 1 week, treatment with biomaterials UBM (regenerative, IL-4 tissue environment), synthetic (fibrotic, IL-17–rich environment), or saline (wound control). (B) Heatmap of differentially expressed genes. Up to 200 cells per cluster are shown, ordered by cluster, with the top 10 differentially expressed genes. Functionally relevant genes from terminal clusters are annotated. (C) Dimensional reduction projection of cells onto two dimensions using UMAP. Cells are colored by experimental biomaterial condition (top) and computationally determined cluster (bottom). (D) Summary of cluster differentiation trajectories, composition by experimental origin, markers, and biological functions generated by bioinformatics analysis. (E) Flow cytometry strategy informed by computationally determined markers including CD9, CD301b, and MHCII differentiating in vivo macrophage subsets from UBM or synthetic biomaterial. Subsets are colored equivalent to computational clusters, back-gated into tSNE projection.

Sorting and analysis resulted in the capture of ~7200 macrophages with an average read depth of ~50,000 reads per cell across ~13,000 genes, with more than 4000 median unique molecular identifier (UMI) counts per cell (fig. S1 and table S1). By condition, the total number of macrophages captured was 3343 from UBM, 2919 from PCL, and 876 from saline. Unbiased clustering algorithms categorized macrophages into clusters based on global gene expression patterns. We first computationally pooled macrophages from regenerative (UBM), fibrotic (PCL), and control (saline) conditions to create a virtual aggregate. Cells with reduced signal after scaling were removed from the analysis (fig. S2), leaving nine computationally determined clusters. The clusters were largely enriched for regenerative or fibrotic macrophages, with differential expression analysis confirming distinct gene expression profiles (Fig. 1B). UMAP (uniform manifold approximation and projection), a dimensional reduction algorithm, grouped cells by cluster, indicating that UMAP and the clustering algorithm agreed on the similarity of cell phenotypes (Fig. 1C).

The experimental origin of each cell was then superimposed on the UMAP plot so that the enrichment of cells from regenerative or fibrotic microenvironment could be identified across all cell clusters (Fig. 1C). Cell clusters were distinguished by experimental conditions, so we termed macrophages from the UBM tissue microenvironment regenerative-associated clusters and those from the PCL tissue microenvironment fibrosis-associated clusters. Potential batch effects were removed by scaling on each experimental condition. Groups of macrophages isolated from saline-treated wounds were found at the interfaces between fibrotic and regenerative macrophages. This intermediate profile is supported by flow cytometry data, suggesting that macrophages in a muscle wound without a biomaterial have characteristics of fibrotic and regenerative microenvironment (Fig. 2A). Similar cell numbers were sequenced from the fibrotic and regenerative conditions, but two-thirds of the macrophage cell clusters were regenerative-associated clusters, suggesting that there is increased complexity of macrophage phenotypes in the regenerative tissue microenvironment.

Fig. 2 scRNAseq reveals surface markers that discriminate diverse regenerative and fibrotic macrophage subsets.

(A) Histogram of CD86 and CD206 expression of bulk macrophages from regenerative, fibrotic, and saline microenvironments (green, red, and blue) by flow cytometry. (B) Feature plots of scRNAseq expression levels of Cd86 (bottom) and Cd206 (Mrc1, top) expression superimposed on UMAP plots of cells from scRNAseq. Circled borders mark regions enriched for cells from the regenerative or fibrotic experimental condition. (C) Cluster expression of canonical M1 (red) and M2 (green) markers. Expression levels are given as cluster averages normalized to the maximum value per gene. (D) Gene expression (UMI count) of M1 and M2 markers per single cell. Cells are ordered and colored by condition (regenerative as green, saline as blue, and fibrotic as red). Violin plots of cluster gene expression of surface markers identified by differential expression analysis. (E) In vivo flow cytometry strategy using CD9 and CD301b differentiates fibrotic (F1 and F2) and regenerative (R1 and R2) macrophage subsets in vivo. A.U., arbitrary units. (F) Mean fluorescence values indicate subset-specific expression of activation markers CD86 and CD206 (n = 4, **P < 0.005, ****P < 0.0001). (G) Comparison of scRNAseq (SC) and NanoString (NS) expression profiles after FACS of macrophage subsets. Genes are shown as either differentially expressed in both methods with agreement in fold change direction (SC and NS), found differentially expressed in one dataset (SC or NS), not found in either dataset (Not found), or found in both with disagreement in direction of fold change (Disagreement).

To identify relationships between cell clusters and differentiation trajectories, we performed Slingshot pseudotime and RNA velocity analysis on the regenerative- and fibrosis-associated clusters. Regenerative precursor (RP) and fibrotic precursors (FPs)—such as RP1, RP2, and FP1—were selected on the basis of similarities in gene expression in clusters across experimental conditions (fig. S3). RNA velocity, which predicts cell movement on a ~32-hour time scale, confirmed movement of cells from RP2 toward R1 and supported the defined clusters. Pseudotime results indicate a branching lineage in both the regenerative-associated (R1 and R2) and fibrosis-associated (F1 and F2) clusters (Fig. 1D), with two functionally specialized terminal clusters in each condition. R3 was excluded from the pseudotime analysis because of its gene expression profile that included muscle-related genes.

To enable identification of the terminal regenerative and fibrotic macrophages, we determined surface marker combinations from scRNAseq that could differentiate subsets experimentally (Fig. 1D). We performed flow cytometry on cells isolated from the UBM, PCL, and saline treatment conditions using the computationally identified cluster surface markers. The CD45+Ly6cF4/80hi cell populations from all conditions were concatenated together to create a t-distributed stochastic neighbor embedding (tSNE) plot containing a complex mixture of all macrophages. We then identified macrophages expressing the surface markers CD9 (a protein involved in cell adhesion, fusion, and motility), CD301b (a galactose C-type lectin), and major histocompatibility complex class II (MHCII) in the aggregated dataset to represent the computational macrophage clusters. The four terminal clusters F1, F2 (and FP1), R1, and R2 could be separated in the aggregate, suggesting that the subsets can be readily identified experimentally using flow cytometry (Fig. 1E).

Expression of canonical polarization markers CD86 and CD206 is distributed across macrophage clusters

We first explored the correlation of the unbiased single-cell clusters with canonical M1 and M2 polarization markers. Flow cytometry analysis of macrophages confirmed the enrichment of CD206 in the regenerative condition and CD86 in the synthetic condition with saline or untreated wound exhibiting intermediate levels of both markers (Fig. 2A). Histograms were consistent with previous studies that found that UBM treatment down-regulates CD86, whereas CD206 remains constant and PCL slightly decreases CD86 and substantial decreases CD206 compared with saline treatment.

Superimposing Cd206 and Cd86 on the UMAP plot showed enrichment for regenerative and fibrotic macrophages, respectively, but it showed notable heterogeneity and it could not discriminate between phenotypically distinct subsets, as expression levels of neither Cd206 nor Cd86 correlated with the computationally determined clusters (fig. S4). Expression patterns of additional canonical polarization genes across the unbiased clusters found similar disparities.

Comparison of macrophage polarization markers on a per-cell basis in the different experimental conditions also revealed a considerable heterogeneity (Fig. 2C). Retnla was the only type 2 gene not expressed in the fibrotic macrophages. Expression of other type 2 genes, Arg1 and Socs3, did not parallel Cd206 expression on a cell-by-cell basis. At the same time, high levels of Retnla, Arg1, and Socs3 expression were found in cells that did not express Cd206. There was a similar pattern of disparity in expression of Cd86 and type 1 genes. Cd86 expression did not correlate with Nfkbiz and Il1b on a cell-by-cell basis. Many cells expressed high levels of Cd86 expression in parallel with low levels of Nfkbiz and Il1b.

Expression of CD9 and CD301b distinguishes macrophages from fibrotic and regenerative environments

Because the single-cell analysis confirmed that Cd86 and Cd206 expression did not differentiate phenotypic subsets, we explored alternative surface markers in the scRNAseq dataset. Assessment of surface markers revealed that Cd301b, Cd9, and Cd74 (encoding the MHCII-associated invariant chain) were sufficient to identify each of the macrophage clusters (Fig. 2D). We then tested these surface markers on bulk cell isolates from the tissue environments to confirm that the gene expression correlated with surface protein expression and could separate the macrophage populations using multiparametric flow cytometry. The proposed surface markers were able to discriminate the macrophage subsets corresponding to computationally determined clusters in the regenerative and fibrotic conditions (Fig. 2E) using the gating schematics in figs. S5 and S6. The surface marker CD301b allowed separation of macrophages specific to the regenerative microenvironment, whereas CD9 and CD74 further separated populations into the scRNAseq-predicted subsets with distinct CD206 and CD86 profiles (Fig. 2F). The CD9 and CD301b surface marker paradigm differentiates in vivo macrophage subsets not equivalent to the commonly used cytokine-induced in vitro macrophage phenotypes (IL-4, IFNγ + lipopolysaccharide, and IL-10) that all express high levels of CD9 (fig. S7).

Because the surface markers could separate the macrophage subpopulations, we could further validate the gene expression of the subsets. We sorted the R1, R2, F1, and F2 + FP1 macrophage populations by flow cytometry and compared gene expression using NanoString. Comparison of differential expression from single-cell and NanoString data for 92 genes found in both datasets across the four clusters found 142 cases of agreement in fold change for differentially expressed genes and only 5 cases of disagreement between the two methods (Fig. 2F).

Computational phenotyping reveals inflammatory and phagocytic macrophages in the regenerative environment

We used pseudotime analysis to elucidate relationships between the regenerative-associated clusters and found a branching differentiation trajectory with two terminally differentiated clusters (R1 and R2) stemming from three precursor clusters (RP1 to RP3) (Fig. 3, A and B). Precursor RP3, a direct precursor to only R2, shared a similar but reduced phenotype to R2 based on differential gene expression (fig. S8). Although R2 is composed entirely of macrophages derived from the regenerative condition, RP3 had a relatively high composition of saline-derived macrophages (fig. S9). Flow cytometry analysis of R2 confirmed that it is also present in the saline control wound but appears at later time points (Fig. 3H). The regenerative-associated clusters expressed Il4ra with R2 expressing the highest levels, correlating with the UBM induction of IL-4 and the macrophage response to the tissue environment induced by the biological scaffold.

Fig. 3 Regenerative-associated clusters reveal distinct macrophage phenotypes including inflammatory R1 and phagocytotic R2.

(A) Predicted lineage schematic of regenerative subsets from Slingshot pseudotime analysis and subset surface markers. (B) Slingshot pseudotime trajectory of regenerative-associated clusters shown on a principal component plot (PC1 versus PC2). Cells are colored by cluster. (C) Heatmap of top 20 differentially expressed genes in a comparison of R1 and R2. nES, normalized enrichment score. (D) Violin plots for differentially expressed genes comparing R1 and R2. (E) Gene set enrichment comparing R1 and R2. Plots with higher peaks (red) indicate enrichment of gene sets in R2, whereas plots with negative peaks (blue) indicate enrichment of gene sets in R1. (F) Gene network plots of R1 and R2 generated. Nodes represent genes with connections generated by STRING metadata analysis. Sets of genes associated with specific functions are annotated. (G) Flow cytometry gating scheme validates in vivo protein marker combination for R1 and R2. Macrophages defined as F4/80hi+ from live, CD45+. (H) One- to 6-week time course of R1 and R2 subsets in UCM, PCL, and saline microenvironments (n = 4, biologically independent). Two-way analysis of variance (ANOVA) with subsequent multiple testing P values is presented. ****P < 0.0001. n.s., not significant. (I) Immunofluorescence histology of defect region for CD301b (red) and F4/80 (green) at 1 week post-VML with UBM (scale bars, 25 μm; arrows indicate colocalization).

To characterize the phenotype and potential functional properties of R1 and R2, we compared outcomes of their differential expression, gene set enrichment analysis (GSEA), and gene network analysis. Differential gene expression analysis of R1 showed up-regulation of antigen-presenting capacity, inflammatory activity, and glycolysis, including Cd74, Ccr2, Il1b, and Gapdh (Fig. 3, C and D). Although R1 expression levels of the inflammatory genes Il1b and Tnfa were elevated when compared with other regenerative-associated clusters, they were low when compared with the highly inflammatory fibrosis-associated clusters (fig. S10). Gene set enrichment of R1 also found elevation of inflammatory responses (Fig. 3E). R1 was enriched in leukocyte activation gene sets, suggesting that these macrophages play a role in communication and activation of the adaptive immune system. Last, network analysis supported both the differential expression and gene set enrichment. R1 expressed gene modules associated with glycolysis (Eno1 and Gapdh), antigen presentation (H2 genes and Cd74), and inflammatory cytokines and chemokines (Cxcl1, Ccr2, Ccl5, Tnfa, and Il1b) (Fig. 3F).

In contrast to the inflammatory R1 macrophage subset, R2 expressed multiple genes associated with alternatively activated or anti-inflammatory macrophages (Fig. 3C). Differential expression and network analysis both revealed enrichment for anti-inflammatory genes such as Chil3, Cd163, and Mrc1 (gene encoding CD206) (Fig. 3D). In addition, R2 expressed high levels of Ccl24 (Eotaxin-2), a chemokine for eosinophil attraction that is observed responding to ECM materials (18). Gene set enrichment and gene modules from network analysis also support a metabolic profile associated with oxidative phosphorylation by expression of Cox5a, Uqcrq, Ndufa1, and Ndufc2 (Fig. 3, E and F), in contrast to a glycolytic profile in R1. R2 expression also included gene set enrichment and endocytic gene modules (Cltc, Clta, and Ap2a2) that suggest phagocytic activity in this subset.

The scRNAseq-defined clusters R1 and R2 were validated in vivo by flow cytometry using a combination of CD301b, CD9, and MHCII on cells isolated from the UBM tissue microenvironment. Although CD301b separated the terminal regenerative clusters R1 and R2 from progenitor clusters and R3, we found that the R1 and R2 could be further defined as CD9+MHCII+ and CD9MHCII, respectively. A combination of both markers provided distinct separation of R1 from R2 for analysis and cell sorting (Fig. 3G). As predicted computationally, CD301b+CD9+ R1 and CD301b+CD9 R2 macrophages express different levels of both MHCII and CD11c when quantified by surface marker expression using flow cytometry (fig. S11). R1 and R2 are distinguished by higher CD11c and lower CD206 expression, respectively. These pronounced subset-specific profiles elucidate earlier suggestions of complex phenotypes associated with the regenerative UBM microenvironment (18, 19).

To evaluate the kinetic evolution of the regenerative macrophage subsets in different experimental conditions, we performed flow cytometry on macrophages isolated from regenerative, fibrotic, and control volumetric muscle loss (VML) microenvironments at 1, 3, and 6 weeks (Fig. 3H). Consistent with RNA velocity analysis and prediction (fig. S12), R1 increased significantly in the regenerative condition from 1 to 3 weeks, whereas the R2 population maintained elevated levels with respect to the other conditions. As expected, the fibrotic condition had low expression of both R1 and R2 at all time points. The saline wound control had increased levels of both R1 and R2 at 3 weeks and R2 at 6 weeks, suggesting the UBM tissue microenvironment may be following an accelerated course of regeneration that is observed functionally with respect to an untreated wound. Although UBM and wound controls are significantly enriched in R1 and R2, they are present in low levels in the PCL implants, highlighting the complexity of the wound and biomaterial environment.

Expression analysis predicted that R2 macrophages were phagocytic, an important function in wound cleanup that is associated with tissue repair. Sorted R1 and R2 macrophages were evaluated for the ability to phagocytose fluorescent microbeads. The R2 (CD301b+CD9) macrophage subset phagocytosed beads, whereas R1 (CD301b+CD9+) macrophages did not internalize any beads (fig. S13). This result confirms the functional properties predicted by the gene expression analysis and provides surface markers that can specifically identify phagocytic macrophages. Flow cytometry sorting studies confirmed that gene expression of surface markers correlated with protein expression. We further validated protein expression in vivo in tissues without cell isolation and found Cd301b expression on macrophages (F4/80) for the regenerative UBM condition, and saline control, but not PCL implants (Fig. 3I and fig. S14).

Last, R3 expressed many genes associated directly with muscle tissue including Mylpf, Myl1, Acta1, Tnnc2, and Tnnt3 as well as genes associated with dendritic cells, Lag3 and Cd11c (fig. S15A) (20). Gene set enrichment and network analysis supported this finding, with gene modules associated with skeletal muscle (fig. S15B) and strong enrichment of myogenic gene sets (fig. S15C). Because the clusters were also enriched for genes associated with endocytosis (fig. S15D), two major possibilities exist for this cluster. It could be a cluster of myeloid cells assisting in muscle production, or the presence of muscle-related transcripts could be the result of endocytosing muscle debris.

Fibrosis-associated macrophages include distinct subsets F1 (MHCIIhi+) and F2 (CD9hi+IL-36γ+)

We also applied pseudotime and RNA velocity analysis to determine fibrosis-associated cluster differentiation and polarization relationships (fig. S12). From pseudotime analysis, one precursor cluster (FP1) led to a branching trajectory with two terminal fibrosis-associated macrophage clusters, F1 and F2 (Fig. 4A). RNA velocity confirmed pseudotime results and continued polarization of F2 from the bulk macrophages. The F1 cluster expressed traditional markers of inflammation, and genes associated with the IFN response including Irf7, Il18, and Tlr2 (Fig. 4B). Gene sets for both IFNα and IFNγ responses were enriched in F1 (Fig. 4C). Gene networks also showed modules associated with IFN response (Stat1, Myd88, Irf7, and Tlr2) and cytokines associated with inflammatory function (Il18, Ccl4, Ccl7, and Cxcl10) (Fig. 4D). Although both macrophage populations have elements of inflammation, F1 and F2 exhibit significantly different expression profiles.

Fig. 4 Fibrosis-associated macrophages include distinct subsets F1 (MHCIIhi+) and F2 (CD9hi+IL-36γ+).

(A) Lineage schematic of fibrotic subsets from Slingshot pseudotime analysis and trajectory including descriptive marker combination. Pseudotime trajectory is shown in a principal component plot (PC1 versus PC2). (B) Fibrotic subsets are distinguished by specific marker. (C) Heatmap of gene set enrichment scores normalized across clusters for gene sets found up-regulated in F1 and running gene set enrichment plots for the IFNγ and IFNα responses. (D) Gene network representation for relationships of differentially expressed genes in F1 (top) and F2 (bottom) by STRING metadata scores. (E) Flow cytometry gating strategy specific to F1 and (F2 + FP1) from F4/80hi+ macrophages using CD9, MHCII, and CD301b. (F) Time course of F1 and (F2 + FP1) subsets in UBM, PCL, and saline microenvironments (n = 4, biologically independent). Two-way ANOVA P values are presented. *P < 0.05; ***P < 0.0001. (G) Immunofluorescence histology for IL-36γ (violet) and F4/80 (green) at 1-week VML with synthetic PCL material. Low power indicates location of PCL material and muscle defect; arrows detail colocalization (scale bars, 500 and 50 μm).

Although F2 was small, it had a distinct gene expression signature that included recently discovered cytokines and genes with limited characterization. Although F2 expressed inflammatory markers Slpi, Hdc, Tlr2, and Il1b (fig. S16), it also expressed genes associated with autoimmunity: Il36γ, Trem1, Asprv1, and Il17ra (Fig. 4B). The unique nature of this subset and possible functional relevance is exemplified in the expression of Il36γ (also known as Il1f9) that is found in lesions of skin psoriasis and participates in a positive feedback loop in type 17 immune responses (21). IL-17, the primary cytokine of type 17 responses, is also associated with fibrotic diseases not yet associated with autoimmunity, including idiopathic pulmonary fibrosis (22, 23), cardiac fibrosis (24), and the FBR, suggesting a type 17 autoimmune connection. The F2 cluster expressed increased Il17ra, further supporting the role of IL-17 in this subset and the macrophage responsiveness in the PCL tissue microenvironment.

To validate experimentally the fibrosis-associated clusters, we performed flow cytometry using the computationally defined surface marker strategy including CD9, CD301b, and MHCII (Fig. 4, E and F). According to the computational analysis, the surface markers differentiate the F1 cluster, but the F2 and FP1 clusters are not separated. Both F1 and F2 + FP1 were elevated in the PCL fibrotic condition at 1 week with a significant reduction in F1 from 1 to 3 weeks. However, levels of F2 + FP1 were consistently elevated, suggesting that F2 may be important in the development of pathological fibrosis. The UBM regenerative condition had low levels of both F1 and F2 at 1 week, which was reduced at 3 and 6 weeks. The control wound condition had high levels of F1 at 1 and 3 weeks with a significant decrease at 3 weeks. There was a low, decreasing number of F2 macrophages in the control wound environment, further suggesting that F2 is associated with fibrotic pathology. Again, the complexity of the wound environment is apparent with all macrophage subsets present at early time points in each condition. However, enrichment of macrophage subsets in specific tissue environments predicted their expression at later time points as seen with F2 that is the primary macrophage in PCL implants at 6 weeks. Last, we performed additional protein-level validation for IL-36γ and CD9 on F4/80+ macrophages with immunofluorescence, and coexpression of these markers representing the F2 subset was visible in tissue environments containing PCL (Fig. 4G and fig. S17).

Generation of F2 macrophages is dependent on IL-17

The F2 macrophage subset was enriched in the PCL condition that induces pathological fibrosis as part of the FBR. Because the F2 macrophages expressed IL-36γ and genes associated with IL-17 immune responses, we further evaluated the F2 population and its dependence on IL-17 signaling. PCL was implanted in mice lacking IL-17A (Il17a−/−) or the IL-17 receptor (Il17ra−/−) to knock down IL-17 signaling. Fibrosis in response to PCL decreased in Il17a−/− mice but was only completely abrogated in Il17ra−/− mice (12). Connecting this functional outcome with macrophage responses, immunostaining for CD9 and F/480 decreased in Il17a−/− and Il17ra−/− mice treated with PCL compared with wild type (WT) after 12 weeks (Fig. 5A). Il36γ, coding for a key cytokine expressed in F2 that links IL-17 and autoimmune responses in the tissue, significantly increased after PCL implantation compared with saline in WT animals but decreased significantly in Il17a−/− and Il17ra−/− mice (Fig. 5B). Over the 6 weeks of PCL biomaterial implantation, Il36γ, along with Trem1 and hallmark genes of fibrosis Col1a1 and Col3a1, significantly increased from 1 to 3 weeks and persisted through 6 weeks (Fig. 5C). Macrophages cocultured with TH17 T cells increased expression of Tnfa, Il36γ, Trem1, and Il17ra (fig. S18A). Cultured fibroblasts exposed to IL-36γ increased expression of Nfkbiz, an early regulator linked to TH17- and IL-36γ–dependent signaling (25) as well as additional Il36γ itself (fig. S18B). Overall, these data support a connection between IL-17 signaling and IL-36γ in F2 and fibrotic macrophage subsets with potential implications toward fibrosis.

Fig. 5 Profibrotic CD9hi+IL-36γ+ macrophages are dependent on IL-17 signaling, and terminal clusters are relevant in human pathologies.

(A) Immunofluorescence staining for mouse macrophage marker F4/80 and CD9 in WT, Il17ra−/−, and Il17a−/− mice 12 weeks after implantation with PCL (scale bars, 50 μm). (B) Il36γ gene expression in WT, Il17ra−/−, and Il17a−/− mice with PCL normalized to saline controls (n = 4, biologically independent; ANOVA with multiple comparison, ***P < 0.001). (C) Time course of Il36γ signaling and fibrosis-related gene expression over 6 weeks, normalized to saline controls (n = 4, biologically independent; ANOVA with multiple comparison, ***P < 0.001). (D) Immunofluorescent staining for CD64-, CD9-, and IL-36γ–positive macrophages in human breast implant tissue capsules (scale bars, 50 μm), juvenile xanthogranuloma, and Langerhans cell histiocytosis (scale bars, 200 μm). (E) Gene expression correlations of human IL17RA with IL36γ, IL17RA, and CD9 with MSR1 in human breast implant fibrotic capsules. (F) Network diagrams and similarity heatmaps for terminal fibrotic and regenerative macrophage clusters to clusters from repository single-cell RNA datasets for murine models of cancer (sarcoma ± immunotherapies aCTLA-4, aPD-1), lung fibrosis (± bleomycin induction), and human liver. Circles represent percent compositions of clusters by condition.

The relevance of the fibrosis-associated F2 profile was explored in human fibrotic pathologies. Positive immunofluorescence staining for F2-specific markers (CD64+CD9+IL-36γ +) in human breast implant tissue capsules as well as histiocytosis (juvenile xanthogranuloma and Langerhans cell histiocytosis) suggest that these macrophages are also relevant to human disease (Fig. 5D). Furthermore, IL17RA expression is correlated to IL36γ and the expression of the macrophage marker MSR1, supporting a role of the IL-17/IL-36γ feedback loop in human macrophages (Fig. 5E).

To determine the broader relevance of the model biomaterial-induced macrophages, we used the SingleCellNet program (26) and trained the scRNAseq cell classification algorithm using the terminal macrophage datasets (fig. S19). We then compared the results with macrophages computationally extracted from publicly available datasets of healthy human liver (27), a mouse model of idiopathic pulmonary fibrosis (28), and a mouse model of sarcoma with cancer immunotherapy (6). After clustering macrophages within each dataset, we applied SingleCellNet to quantify similarity with the terminal biomaterial macrophages (Fig. 5F). The R1 and R3 macrophage subsets were found in all of the conditions that we evaluated. Although the expression of muscle markers in R3 initially suggests that they may be unique to the muscle tissue environment, their presence in all the conditions studied suggests that even this subset has broader relevance in multiple tissue types and pathologies. We found that only liver, a strongly regenerative tissue, contained macrophages similar to R2. Most of the macrophages in the sarcoma were similar to F1, but there was also a small cluster similar to F2. The F2 cluster was also present in the lung tissue.

DISCUSSION

Here, we applied scRNAseq to identify and characterize macrophage phenotypes associated with tissue microenvironments modeled using biomaterials that promote divergent tissue responses, immune profiles, and functional outcomes. The unbiased classification and characterization from the single-cell analysis provide distinct phenotypic profiles that can be identified using a combination of surface markers for the differentiated macrophages that are typically used to characterize M1 and M2 polarization. Terminal clusters discovered from single-cell analysis provide refined phenotypic and functional macrophage characterization in different tissue environments. We identified a previously unknown macrophage population that links a fibrotic tissue environment with type 17 immune responses and signatures of autoimmunity, which was abrogated with loss of IL-17 signaling.

The macrophages associated with UBM and IL-4 in the tissue are heterogeneous and distributed in phenotypically distinct clusters. IL-4 is a cytokine recognized for promoting repair of muscle (9), liver (29), and cartilage (30) and is critical for macrophage polarization in a healing wound (31). The UBM environment induced greater macrophage heterogeneity with two primary terminal subsets with phenotypes relevant to tissue repair. Expression analysis of the R1 cluster suggests that it is important for mobilizing and educating immune cells through the expression of chemokines and increased antigen presentation that is required during the early wound healing process. The R2 macrophages, with the highest level of Il4ra, expressed genes relevant to stimulation of other cell types important for type 2 responses and regeneration including Ccl24 (Eotaxin-2), coding for a protein that attracts and activates eosinophils. Both macrophage populations may be acting, in part, by mobilizing and instructing the immune response. This finding is supported by Heredia et al. (32), who demonstrated that the IL-4–secreting eosinophils are critical to muscle repair. The metabolic profiles of R1 (glycolysis) and R2 (oxidative phosphorylation) correlate with distinct functions of antigen presentation and adaptive-related chemokine expression versus phagocytosis, which was validated experimentally in sorted R2 macrophages. Glycolysis has been associated with inflammatory macrophages (33) and oxidative phosphorylation with alternatively activated macrophages (34). In vitro studies of conventional M2 macrophages required inhibition of both metabolic pathways to inhibit IL-4–induced STAT6 (signal transducer and activator of transcription 6) phosphorylation (35).

The distribution of macrophages isolated from the PCL-treated wounds was less heterogeneous than the ECM-treated tissues and included the functional subsets F1 and F2. The F1 cluster expresses many genes associated with inflammation including IFN-related cytokines and activation of the innate and adaptive immune system. The R1 cluster also expressed markers of inflammation and mobilization, but the magnitude of expression and types of inflammatory markers were significantly different. This difference in the F1 and R1 inflammatory profile suggests the importance of the early inflammatory response in directing the subsequent tissue repair or development of a foreign body capsule or fibrosis. The time course of flow cytometry revealed that R1 increased with ECM treatment. Because ECM treatment improves tissue repair, R1 may represent an inflammatory phenotype that can be targeted to enhance tissue development.

The F2 cluster associated with PCL treatment expressed genes that connected type 17 immunity and markers of autoimmune disease. Type 17 immune responses are associated with autoimmunity in diseases such as psoriasis, irritable bowel syndrome, and inflammatory arthritis (21, 3638). The F2 macrophages express IL-36γ, a cytokine that is found clinically in the skin of patients with psoriasis and in inflammatory arthritis (39). This cytokine was also recently identified as a target in tumors that, when blocked, enhances responsiveness to immunotherapy (40). It is also implicated in a positive feedback loop with IL-17 (21). In other work, we demonstrated that IL-17 is produced by innate lymphocytes, γδ, and CD4+ T cells in response to PCL implantation in mice and in the fibrous capsule surrounding human breast implants (12). IL-17 is implicated in fibrotic disease in the lung (23), heart (24), and liver (41) in addition to the FBR (42). The F2 macrophage population expressing IL-17 receptor A links IL-17 signaling, fibrosis, and autoimmune disease.

The F2 macrophages also expressed multiple forms of Trem (triggering receptor expressed on myeloid cells) and Trem ligand that are associated with autoimmune diseases such as inflammatory bowel disease and psoriasis (43, 44). TREM integrates and broadly modifies inflammatory signals across the innate and adaptive immune system. The presence of F2 surface markers and related cytokines in human tissues provides evidence that the tissue immune environment created by PCL and the mechanisms of response may be broadly relevant to various pathological conditions. Further supporting the broader relevance of the macrophage subsets, we found macrophage clusters similar to F2 in publicly available datasets of idiopathic lung fibrosis and sarcoma. Multiple genes in the F2 subset have functions that remain unknown. On the basis of the potential importance of this subset in disease pathology, further studies into these unknown genes may be warranted. Although small in number, the F2 macrophage subset is highly differentiated and may play a critical role in pathologies associated with the proinflammatory macrophages in tissue fibrosis.

Further studies to investigate the macrophage subsets that we identified by single-cell analysis will determine whether the surface markers and their associated subpopulations and respective expression profiles are stable and broadly relevant. We characterized only a subset of macrophages in this study and further analysis into expression profile and function of other macrophage and myeloid populations in the tissue environment that affect healing and fibrosis outcomes. Additional spatial profiling of the macrophage subsets will enable understanding of their communication with other immune cell types and proximity to cells involved directly in healing or fibrosis.

MATERIALS AND METHODS

Experimental model and subject details

Surgical procedures and implantation

All animal procedures were performed in adherence to approved Johns Hopkins University Institutional Animal Care and Use Committee protocols. Animals included female mice WT C57BL/6J (Jackson Laboratories) and transgenic Il17a−/− (Y. Iwakura, University of Tokyo, Tokyo, Japan) and IL17ra−/− (T. Mustelin, Amgen, Seattle) at ages from 8 to 10 weeks. The bilateral traumatic muscle defect was created as previously described (9). The defects were filled with 30 mg of a synthetic material or biological scaffold material. PCL was used as a synthetic material (particulate, Mn = 50,000 g/mol, mean particle size <600 μm, Polysciences). In turn, as a biological scaffold material, decellularized UBM (MatriStem, ACell) was implanted from 0.05 ml of a suspension (400 mg/ml) in phosphate-buffered saline (PBS). Control surgeries were injected with 0.05 ml of PBS as a no implant control. All materials were ultraviolet-sterilized before use. Mice were given subcutaneous carprofen (Rimadyl, Zoetis) at 5 mg/kg for pain relief. For the sample harvest, mice were euthanized at 1, 3, 6, and 12 weeks after surgery.

Specimen harvest

We obtained murine samples by dissecting the quadriceps femoris muscle followed by fine dicing. Tissues were digested for 45 min at 37°C with Liberase TL (1.67 Wünsch units/ml) (Roche Diagnostics) and deoxyribonuclease I (0.2 mg/ml) (Roche Diagnostics) in RPMI 1640 medium (Gibco). The digested tissues were ground through 70-μm cell strainers (Thermo Fisher Scientific), rinsed with PBS + 0.05% bovine serum albumin, and then washed twice with 1× PBS. The enriched single-cell suspension was washed and stained with the following antibody panels, respective to the application. There were no noticeable differences in single-cell isolation from different material environments.

Flow cytometry and fluorescence-activated cell sorting

For cell isolation using fluorescence-activated cell sorting (FACS), suspensions of single cells from digested muscles were stained for 20 min at 4°C using Viability Dye eFluor 780 (eBioscience). For single-cell analysis, staining in an antibody panel was conducted for 30 min at 4°C including F4/80 phycoerythrin (PE)–Cy7 (BioLegend), CD11b Alexa Fluor 700 (BioLegend), CD64 PerCP-Cy5.5 (BioLegend), MHCII (I-A/I-E) Alexa Fluor 488 (BioLegend), CD3 allophycocyanin (APC) (BioLegend), Ly6c Brilliant Violet 510 (BioLegend), CD45 Brilliant Violet 605 (BioLegend), and Fc Block TruStain fcX (anti-mouse CD16/32) (BioLegend) (table S2). Sorting of murine macrophages (CD45+F4/80hi+CD64+) was performed from live, CD45+CD3 using a BD FACSAria II (BD Biosciences). The post-sort purity of macrophages used for scRNAseq was >98%. For subset analysis, an Attune NxT Flow Cytometer (Thermo Fisher Scientific) after viability stain, a panel comprised of F4/80, CD9, CD11c, MHCII, CD301b, CD206, CD86, Ly6c, CD45, and Fc BLock TruStain FcX (table S3). Data analysis was performed in FlowJo Flow Cytometry Analysis Software (TreeStar). For tSNE projection of multidimensional flow cytometry datasets, downsampling of macrophages (CD45+F4/80hi+) was performed to 50,400 cells total of samples from regenerative, fibrotic, saline control conditions (n = 4, biologically independent). Then, samples were concatenated and tSNE projection was computed (iterations = 1000, perplexity = 30) (fig. S20). For visualization of separation macrophage, subtypes were back-gated onto two-dimensional projection.

Histopathology and immunofluorescence microscopy

After harvest, implanted tissues at 1, 3, 6, and 12 weeks after implantation were fixed in 10% neutral buffered formalin for 48 hours before ethanol and xylene dehydration. Archival formalin-fixed paraffin-embedded tissues from patients with conditions known to involve increased tissue macrophages, such as Langerhans cell histiocytosis, juvenile xanthogranuloma, and dermal scar, were obtained from the Johns Hopkins Hospital surgical pathology archives. Human tissue samples including silicone breast implants were acquired from patients undergoing implant exchange or replacement surgery exemption IRB00088842. The average age was 56 years old (range of 41 to 70 years old), and the average implant residence time was 41 months (range of 1 to 360 months). Paraffin-embedded sections of 5 μm thickness were produced using a microtome (Leica RM2255 microtome). Immunofluorescence staining for CD9, IL-36γ, F4/80, and CD301b was performed with a tyramide signal amplification method using consecutive staining with Opal-520, Opal-570, and Opal-650, respectively. Histological slides were rehydrated in an ethanol to water gradient. After post-fixation in 10% neutral buffered formalin for 30 min, antigen retrieval was conducted in 1 × AR6 buffer (Perkin-Elmer) at 95°C for 15 min. After cooling, endogenous peroxidases were quenched in 3% (v/v) aqueous peroxide (Sigma-Aldrich). All staining was conducted at room temperature after blocking with antibody block/diluent (Perkin-Elmer). For each staining round, the primary antibody in block/diluent (Perkin-Elmer) was incubated at room temperature for 60 min, followed by 10 min of incubation with horseradish peroxidase (HRP) polymer–conjugated secondary antibody, and 10 min of Opal reagent (1:150) in 1× plus amplification diluent (Perkin-Elmer). The antibody HRP complex was stripped at 95°C 1× AR6 buffer (Perkin-Elmer) for 15 min. After quenching of autofluorescence in 0.04% (w/v) Sudan Black (Sigma-Aldrich) dissolved in 70% (v/v) ethanol, slides were then counterstained with 4′,6-diamidino-2-phenylindole (DAPI) for 5 min before being mounted using Dako mounting medium. Imaging of the histological samples was performed on Zeiss Axio Observer with Apotome.2 and Zeiss Zen Blue software version 2.5.

Phagocytosis assay

Macrophage subsets R1 and R2 were sorted according to the gating strategy described in fig. S5. Around 100,000 cells were collected in RPMI 1640 + 5% (v/v) fetal bovine serum (Gibco). Macrophages were seeded at 60,000 cells/cm2 and incubated at 37°C for 12 hours. Then, phagocytosis was conducted for 4 hours at a concentration of 1.0 × 106 cells/ml with Fluoresbrite flash red (Polysciences) with an average diameter of 2.00-μm particles at 2.0 × 108 particles/ml. Last, DAPI and Alexa Fluor 594 anti-phalloidin were used to counter stain macrophages for fluorescent microscopy of bead phagocytosis of fluorescent beads.

Cell culture

For macrophage and T cell coculture experiments, bone marrow–derived monocytes (BMDMs) were isolated and differentiated into macrophage as previously described. At day 0, BMDMs were seeded at 60,000 cells/cm2 and incubated at 37°C for 7 days in medium supplemented with macrophage colony-stimulating factor (100 ng/ml) (PeproTech) into a 12-well insert of a transwell plate (Corning). Naïve CD4+ T cells were isolated from murine lymph nodes at day 4 using a magnetic bead negative selection kit (Miltenyi). Then, T cells were seeded at 1 million cells/ml in 2 ml into six-well plates (Corning) in naïve condition supplemented with IL-2 (100 ng/ml) (PeproTech) or skewed toward the TH17 phenotype using a commercially available kit (R&D Systems). At day 7, differentiated macrophages were exposed to naïve and TH17 T cells at 250,000 cells/cm2 or IL-17A, IL-17F (PeproTech, 100 ng/ml), and control medium for 72 hours.

Fibroblasts were isolated using previously reported methods (45), seeded to T-75 tissue culture plates (Corning), and passaged once before reseeding at 50,000 cells/cm2 for 24 hours and exposed to IL-36γ (R&D Systems) or IL-17A (PeproTech) at 100 ng/ml for 24 hours. Cells were lysed using RLT plus buffer (Qiagen) reagent for quantitative reverse transcription polymerase chain reaction (qRT-PCR).

qRT-PCR gene expression assay

For total and enriched mRNA, lysis was conducted on whole tissues using TRIzol reagent or in vitro cell culture using RLT plus buffer containing 1% (v/v) β-mercaptoethanol. RNA purification was performed with RNeasy Plus Micro and Mini kits (Qiagen). All qRT-PCR was performed using TaqMan Gene Expression Master Mix (Applied Biosystems) according to the manufacturer’s instructions. Briefly, 2 μg of enriched mRNA was used to synthesize complementary DNA (cDNA) using Superscript IV VILO Master Mix (Thermo Fisher Scientific). The cDNA concentration was set to 50 to 100 ng/well (in a total volume of 20 μl of PCR) to match the manufacturer recommendations. The qRT-PCRs were performed on the StepOne Plus Real-Time PCR System (Thermo Fisher Scientific), as TaqMan single-plex assays, using the manufacturer recommended settings for quantitative and relative expression. For tissue samples, β2m was used as the reference gene and samples were normalized to PBS-treated controls. All qRT-PCR data were analyzed using the Livak method, where ΔΔCt values are calculated and reported as relative quantification (RQ) values calculated by 2−ΔΔCt (46). RQ values are represented by the geometric means, with error bars representing the geometric standard deviation. Low-expressing mRNA transcripts were preamplified using the TaqMan Pre-Amp System (Thermo Fisher Scientific) according to the manufacturer recommendations with 14 cycles of amplification with the primer probes of interest (tables S4 and S5).

Single-cell encapsulation and library generation

After sorting of macrophages (CD45+F4/80hi+), single cells were encapsulated in water-in-oil emulsion along with gel beads coated with unique molecular barcodes using the 10× Genomics Chromium Single-Cell Platform. For single-cell RNA library generation, the manufacturers’ protocol was performed (10× Single Cell 3’ v2). Sequencing was performed using an Illumina HiSeq2500 Rapid Mode with 310 million reads per sample and a sequencing configuration of 26 × 8 × 98 (UMI × Index × Transcript read). The Cell Ranger pipeline software was used to align reads and generate expression matrices for downstream analysis.

NanoString gene expression analysis and single-cell data comparison

Macrophage subsets R1, R2 and F1, and F2 + FP1 were sorted according to the gating strategy described in figs. S4 and S5 (table S2). After lysis in RLT plus buffer containing 2-mercaptoethanol, RNA was purified using the RNeasy Plus Micro Kit (Qiagen). The gene expression analysis for n = 3 biologically independent samples was conducted with the nCounter Mouse Myeloid Innate Immunity V2 Panel (XT-CSO-MMII2-12, NanoString Technologies). After quantification of RNA using a Qubit RNA fluorometric assay high-sensitivity kit (Thermo Fisher Scientific), 50 ng of RNA per sample was added to a barcoded probe set reagent and hybridized for 18 hours at 65°C. NanoString data were processed using the nSolver 4.0 software kit. The differential expression data were used for a direct comparison with single-cell clusters. All differentially expressed genes (P < 0.05 and positive fold change) for each cluster were taken for further analysis. Undetected genes in the single-cell dataset were removed. Differential expression of the sorted clusters was performed on the single-cell dataset using the gene set generated from NanoString. If the NanoString and single-cell differential expression both showed significance (P < 0.05) in the same direction for a given gene and cluster, the two techniques were considered in agreement. If there was significance in the opposite direction, the techniques were considered in disagreement. If only one technique found significance, it was noted as “found in one,” and if neither found significance, “not found” was noted for that gene and cluster.

Computational analysis

Sequence alignment, filtering, normalization, and scaling

Alignment was performed with STAR (47) through the Cell Ranger pipeline. Filtering, normalization, and scaling were performed using Seurat (7, 48). Cells with UMI counts for fewer than 200 genes and genes with expression in less than 0.1% of cells were both dropped from analysis. Data were normalized by Enorm = log(UMI*10,000/UMItotal), where UMItotal is total UMI expression for a given gene. Scaling was performed to remove unwanted effects correlated to batch, total UMI count, percent of mitochondrial genes, and cell cycle. To determine cell cycle scores for use with scaling, Seurat’s CellCycleScoring function was used with a previously determined set of genes (49) correlated with G2M or S phase. Data scaling was performed by first fitting a linear model with the parameters to scale out as independent variables (batch, total UMI count, percent mitochondrial genes, G2M score, and S score) (fig. S21). For each gene, the residuals from the fit were Z-scored and used for downstream analysis. Last, principal components were calculated and the top 50 were selected on the basis of leveling of variance per principal component as determined by an elbow plot (fig. S22).

Clustering analysis

SC3 (50) and Seurat’s unsupervised clustering algorithms using the top 50 principal components were compared using silhouette coefficients. For each clustering algorithm, a reasonable range of resolution or k parameters was determined a priori. Silhouette coefficients were determined for each result and compared. Seurat was found to have consistently higher silhouette coefficients and was used. A resolution value was selected on the basis of a combination of high silhouette coefficient and reasonable number of clusters for biological consideration. After clustering, a small population of contaminant fibroblasts was found and removed (fig. S23).

Removal of low signal clusters

Downstream analysis also found evidence of clusters that had little signal remaining after scaling. To attempt to quantify this, differential expression test was calculated using Seurat’s negative binomial statistical test with a minimum log fold change of 1. On the basis of these results together with quality control information such as total UMI count (fig. S2), we decided to remove clusters O1 to O4 from analysis to avoid diluting the signal from clusters with stronger signal. The resulting data matrix was reprocessed identically to find nine clusters of macrophages with strong enough signal to identify surface markers and suggest biological function.

Differential expression and gene set enrichment

Differential expression was calculated using edgeR’s exactTest (51). A list of differentially expressed genes by cluster is provided in table S3. Unless otherwise specified, differential expression for a specific cluster was determined by comparison against all other clusters. For heatmaps comparing expression of multiple different genes by cluster, gene expression values were normalized against the highest expression value for each given gene. Feature plots and violin plots for gene expression were generated with Seurat’s FeaturePlot and VlnPlot, respectively, using log-normalized expression values. Ranked gene lists were used to calculate gene set enrichment scores, and GSEA was performed as previously described using fgsea (52). Briefly, all genes in the analysis were ranked by P value with the exception of the direct comparison between R1 and R2, where genes were ranked by log fold change. For a given gene set, a running enrichment score is calculated by stepping through the ranked gene list and adding to the score when encountering a gene in the gene set or subtracting from the score otherwise. The enrichment score for a gene set is maximum (or minimum if negative) of the running enrichment score with a higher enrichment score, indicating higher concentration of genes in the gene set toward the beginning of the ranked gene list. Last, a null distribution for the gene set is calculated by repeatedly permuting the ranked genes and used to calculate a P value. Gene ontology cellular components, biological process, and hallmark gene sets were obtained from the Broad Institute for GSEA.

Protein network analysis

Functional network analysis of gene products in macrophage clusters was conducted using a databank-based query of the STRING consortium (53, 54). Three-letter codes of the top 300 differentially expressed genes per cluster were mapped to the Mus musculus protein products. Connectivity of proteins in the network was score-ranked from 0 to 1, with 1 being the highest possible confidence weighed metadata comparison based on experimental protein-protein interaction evidence, PubMed text mining, and curated databases. Gephi network visualization tool was used to plot the cluster protein networks (55). Network nodes represent proteins with the size proportional to the number of interactions. Edges represent protein-protein associations, with the thickness proportional to the strength of the meta-evidence supporting the interaction.

Pseudotime analysis

Pseudotime analysis was performed using Slingshot (56). Precursor clusters were determined using differential expression and gene set enrichment for clusters. In particular, we found that precursor clusters shared gene set enrichment patterns (fig. S3). R3 was removed as an outlier due to highly unique gene expression. To determine genes associated with pseudotime progression, genes were regressed on pseudotime as determined by Slingshot using a general additive model with the gam R package. The top 10 genes by P value were then used to generate a heatmap, ordering cells by pseudotime. Z score–scaled residual expression values were used for the heatmaps.

RNA velocity analysis

RNA velocity analysis was performed as previously described (57) using the velocyto.py python package for annotating transcripts as spliced or unspliced, followed by the velocyto.R R package to perform velocity estimation. Briefly, transcripts are marked as either spliced or unspliced based on the presence or absence of intronic regions in the transcript. For each gene, a simple model of RNA dynamics is then fit to the data. Last, the RNA velocity is estimated for each cell by looking for over- or underrepresentation of spliced to unspliced ratios. RNA velocity is visualized on a UMAP plot, with vector fields representing the averaged velocity of nearby cells.

Analysis of publicly available datasets

Publicly available scRNAseq datasets were obtained from Gene expression omnibus [GSE115469 (27), GSE111664 (28), and GSE119352 (6)]. For each dataset, we scaled the data and then used SingleCellNet (26) trained on the tabula muris dataset (58) to quantify similarity to known cell types. We then selected only macrophages and alveolar macrophages and completed analysis up through clustering. After determination of optimal clustering resolution, we applied the SingleCellNet program to quantify similarity to the biomaterial wound macrophages, training the algorithm on the terminal clusters R1, R2, R3, F1, and F2. This generated five similarity scores for each public dataset macrophage, one for each terminal cluster. Last, we calculated cluster level similarity scores by taking the mean of all cells in the cluster.

Supplementary Materials

immunology.sciencemag.org/cgi/content/full/4/40/eaax4783/DC1

Fig. S1. Quality control of scRNAseq data.

Fig. S2. Differential expression and quality control metrics show clusters grouped on nondesirable traits.

Fig. S3. Precursor clusters show similarities across condition.

Fig. S4. scRNAseq reveals that canonical M1 and M2 markers do not correlate with single-cell clusters.

Fig. S5. Macrophage gating scheme for regenerative subsets.

Fig. S6. Macrophage gating scheme for fibrotic subsets.

Fig. S7. Murine macrophage CD9 and CD301b profiles using flow cytometry.

Fig. S8. Differential gene expression heatmap for RP3 and R2.

Fig. S9. Cell composition of clusters by experimental condition.

Fig. S10. Expression of inflammatory genes in fibrotic cluster F2.

Fig. S11. Quantification of expression levels of CD11c and CD206.

Fig. S12. RNA velocity analysis.

Fig. S13. Regenerative subsets R1 and R2 phagocytosis assay.

Fig. S14. Immunofluorescence microscopy of F4/80+CD301b+ macrophages.

Fig. S15. Cluster R3 expresses tissue-specific genes.

Fig. S16. Expression of inflammatory genes in fibrotic cluster F2.

Fig. S17. Immunofluorescence microscopy on F4/80+CD9+ macrophages.

Fig. S18. In vitro TH17 dependence of Il36γ expression in macrophages and IL-36γ induction of fibrosis phenotypes in cultured fibroblasts.

Fig. S19. Analysis pipeline for publicly available datasets.

Fig. S20. Scheme to generate for tSNE projection for virtual aggregate flow cytometry data.

Fig. S21. Scaling metric for the scRNAseq raw data.

Fig. S22. Variance by principle component.

Fig. S23. Expression of fibroblast signature genes in small contaminant cluster.

Table S1. Quality control metrics for 10× scRNAseq.

Table S2. Macrophage FACS panel.

Table S3. Macrophage subtype panel.

Table S4. Murine TaqMan gene expression assay probes.

Table S5. Human TaqMan gene expression assay probes.

Table S6. Sorted cell populations for NanoString Gene Expression Profiling.

Table S7. Differentially expressed genes for scRNAseq clusters.

REFERENCES AND NOTES

Acknowledgments: We would like to thank the Cindy Sears laboratory for providing the original breeding mice strain for transgenic Il17ra−/− and Il17a−/− and the Sidney Kimmel Comprehensive Cancer Center (SKCCC) Immune Monitoring Core for NanoString studies. Funding: We acknowledge funding from the Morton Goldberg Chair, the Bloomberg~Kimmel Institute for Cancer Immunotherapy, and ACell Inc. Author contributions: S.D.S., C.C., and J.H.E. conceptualized and drafted figures and the manuscript. S.D.S., C.C., J.H.E., D.M.P., P.C., and F.H. contributed to experimental design and interpretation of results. S.D.S., C.C., L.C., D.R.M., R.M.S., A.T., J.E.S., and J.M.T. conducted experimental procedures and analyzed data. All authors participated in editing and revising the manuscript text and figures. Competing interests: J.H.E. is an inventor on intellectual property related to biological scaffolds and inhibiting fibrosis. J.H.E. is a consultant to ACell and a founder of Aegeria. D.M.P. is an inventor on intellectual property related to inhibiting fibrosis and immunotherapies for cancer. D.M.P. is a consultant and founder to multiple immunotherapy companies. Data and materials availability: Single cell RNA sequencing data is available from the Gene Expression Omnibus under accession number GSE138027. All other needed to support the conclusions in the paper are available in the paper and/or the Supplementary Materials. Requests regarding protocols and resources should be directed to and will be fulfilled by the lead contact J.H.E. (jhe{at}jhu.edu).
View Abstract

Navigate This Article