Research ArticleIMMUNOLOGY

Harnessing lipid signaling pathways to target specialized pro-angiogenic neutrophil subsets for regenerative immunotherapy

See allHide authors and affiliations

Science Advances  30 Oct 2020:
Vol. 6, no. 44, eaba7702
DOI: 10.1126/sciadv.aba7702


To gain insights into neutrophil heterogeneity dynamics in the context of sterile inflammation and wound healing, we performed a pseudotime analysis of single-cell flow cytometry data using the spanning-tree progression analysis of density-normalized events algorithm. This enables us to view neutrophil transitional subsets along a pseudotime trajectory and identify distinct VEGFR1, VEGFR2, and CXCR4 high-expressing pro-angiogenic neutrophils. While the proresolving lipid mediator aspirin-triggered resolvin D1 (AT-RvD1) has a known ability to limit neutrophil infiltration, our analysis uncovers a mode of action in which AT-RvD1 leads to inflammation resolution through the selective reprogramming toward a therapeutic neutrophil subset. This accumulation leads to enhanced vascular remodeling in the skinfold window chamber and a proregenerative shift in macrophage and dendritic cell phenotype, resulting in improved wound closure after skin transplantation. As the targeting of functional immune subsets becomes the key to regenerative immunotherapies, single-cell pseudotime analysis tools will be vital in this field.


Neutrophils have classically been viewed as a largely homogeneous population of myeloid immune cells that do not have the ability to polarize into subsets, unlike other leukocyte populations with multiple phenotypes, such as T lymphocytes or macrophages (1, 2). Traditional neutrophil activities include phagocytosis, degranulation, and release of neutrophil extracellular traps (NETs) aimed at responding to and combating microbial invasions or tissue injury (3). Neutrophils are considered the immune system’s short-lived “first responders” to infection or injury with limited and well-defined responses (4). While traditional neutrophil functions serve a beneficial role in containing infiltrating foreign cells, excessive neutrophil accumulation, degranulation, and release of NETs can contribute to tissue fibrosis, chronic inflammatory responses, and the development of diseases such as diabetes, atherosclerosis, and cancer (5). Moreover, there is now an emerging evidence for the existence of heterogeneous subsets within the larger neutrophil population that are able to carry out functions in both homeostatic and pathogenic immune responses in which they can actively contribute to regenerative processes following tissue injury (6, 7). Specific subpopulations of neutrophils are also implicated in the early establishment of a healthy vasculature. Specifically, pro-angiogenic neutrophils, identified by their high expression of CD49d, VEGFR1 (vascular endothelial growth factor receptor 1), VEGFR2, and CXCR4, are efficiently recruited to nonvascularized tissues and are capable of inducing angiogenesis (8). In addition, aged neutrophils, identified by their small size and increased expression of CD11b and CD49d, are highly reactive phenotypes that are long lived in the circulation compared with classical neutrophils and serve as first responders in inflammatory reactions (9). We believe that a regenerative medicine approach that selectively targets specialized neutrophil subpopulations holds substantial promise.

Advancements in the understanding of the mechanisms that lead to the resolution of inflammation have also led to the discovery of proresolving lipid mediators that regulate key signaling events that are involved in inflammation resolution. One example of such proresolving mediators is resolvins, which are derived from both eicosapentaenoic acid (EPA) and docosahexanoic acid (DHA) (10). Conversion of membrane lipids into specialized bioactive lipid mediators is a critical step in cessation of polymorphonuclear neutrophil infiltration and active resolution of inflammation (11). Aspirin-triggered resolvin D1 (AT-RvD1) is a D-series resolvin, which only differs functionally from resolvin D1 in that it is resistant to rapid conversion and inactivation in vivo. Thus, with a relatively longer in vivo half-life, it is an attractive small molecule for sustained local delivery to reduce chronic inflammation and promote tissue regeneration (11). We and others have shown that treatment with D-series resolvins is able to promote angiogenesis after both ischemic and nonischemic tissue injury while also reducing the damaging aspects of prolonged neutrophil interactions with implanted biomaterials (12, 13). We believe that the application of advanced single-cell analytical techniques on traditional flow cytometry data could elucidate how nontraditional neutrophil functions may be targeted for immunotherapy using locally administered AT-RvD1 or other bioactive lipids.

Assessing the functional heterogeneity of neutrophil subpopulations based on surface marker expression profiles is notoriously challenging. Recent multicenter studies, such as the Immunological Genome Project, show that manual gating is one of the largest confounding factors when analyzing single-cell data-based experiments (14). Innovative methods of dimensionality reduction are emerging to enhance the robustness of surface marker expression profiles of single cells and their association with functional cell states. In this study, we explore neutrophil heterogeneity in the context of sterile inflammation using murine dorsal skinfold window chamber as a model of excisional skin injury (fig. S1A). The window chamber allows for the longitudinal intravital imaging of the injury microenvironment and single-cell tracking of immune cell invasion in the presence or absence of sustained AT-RvD1 presentation (15). The sustained AT-RvD1 delivery is accomplished using a poly(lactic-co-glycolic acid) (PLGA) thin film, which we have characterized in previous studies (13). We used traditional “manually gated” immune populations, where regions of interest (manual gates) were drawn over a series of bivariate plots, for further analysis (fig. S1B). We used several advanced analytical techniques to identify a pro-angiogenic neutrophil subset that is involved in the early angiogenic process, which is important in the successful resolution of inflammation (8). Because of its unique ability to resolve subtly different cell populations while preserving the global structure, we first applied uniform manifold approximation and projection (UMAP) for dimension reduction and unbiased cluster detection based on the canonical markers of innate immune cells (fig. S1C). We then performed a novel pseudotime analysis based on spanning-tree progression analysis of density-normalized events (SPADE) (16, 17) to reconstruct complex cellular hierarchies of immune cell transitions to reveal rare cell states and order single-cell expression profiles in “pseudotime”—a quantitative measure of progress through a biological process (fig. S1D) (18). In addition, we classified each node to an immune cell subset and calculated the percent cell frequency of each population relative to the total number of events in each node. The response of these identified subpopulations to AT-RvD1 could then be analyzed within SPADE. Because every SPADE node carries its own unique surface marker expression combination, these can be compiled into a heatmap including every node on a SPADE dendrogram (fig. S1E). In addition to the single-cell analysis, we performed bright-field imaging and whole mount immunohistochemistry to prescribe functional significance to the neutrophil subpopulations identified using SPADE (fig. S1, F and G). We developed a poly(ethylene glycol)–maleimide (PEG-MAL) hydrogel to deliver AT-RvD1 in the murine syngeneic skin transplant model (fig. S1H). This model was chosen to chronicle wound healing and tissue integration over time (fig. S1K) and validate biomaterial-mediated AT-RvD1 as a therapeutic approach. Macrophage and dendritic cell heterogeneity has been well characterized in the wound healing response, as M2 macrophages and IL10+ dendritic cell subsets specifically are notable for promoting a proregenerative injury milieu (19, 20). Thus, to establish the rigor of our analytical approach, we identify macrophage and dendritic cell subpopulations using SPADE and investigate their response to AT-RvD1 treatment.

In this work, we show that local, sustained delivery of AT-RvD1 is able to selectively enrich injured tissue with pro-angiogenic neutrophils that have high expression of CD49d, VEGFR2, VEGFR1, and CXCR4 at the onset of the inflammatory response. Furthermore, we show that this shift in the early microenvironment increases the downstream accumulation of M2 macrophages, a proregenerative macrophage subset, and IL10+ tolerogenic dendritic cells, a dendritic cell subset responsible for activating regulatory T cells and graft tolerance (20). Our findings demonstrate the importance of applying single-cell dimensionality reduction and pseudotime analysis techniques for the study of rare cellular heterogeneity. Moreover, it suggests that a novel small molecule–based strategy for enhancing early angiogenesis can be achieved through the targeting of a pro-angiogenic neutrophil subset, which is understudied due to the novelty of neutrophil heterogeneity and difficulty in analyzing neutrophil subset response to therapeutic perturbations (21).


Dimensionality reduction analysis reveals neutrophil heterogeneity

To investigate the role of inflammation progression and AT-RvD1 treatment on neutrophil heterogeneity in the blood and in the tissue, we created an array of UMAPs and SPADE dendrograms (Fig. 1). We collected blood and tissue samples at days 1 and 3 after surgery for both vehicle control and AT-RvD1 treatment groups. Using manual gating techniques, we selected CD11b+ cells as representative of myeloid cells and made UMAP and SPADE projections of this group of cells. By overlaying where the neutrophil population (depicted as red dots) lies on these projections, we show that regardless of time point and treatment, the circulating neutrophil population is constrained to a single island and, thus, is composed of a homogeneous group (Fig. 1, A to D). However, when these neutrophils extravasate into the tissue, they show increased diversity, visualized by the spread of neutrophils among distinct clusters in the UMAP and in different trajectories in the SPADE dendrogram (Fig. 1, E to H). In the tissue, neutrophil heterogeneity increases at day 3 and with AT-RvD1 treatment. When compared with day 1, the day 3 UMAP shows a similar number of CD11b+ islands; however, the spread of Ly6G+ neutrophils among these islands is increased at day 3, which suggests an increased neutrophil heterogeneity. Similarly, at 3 days after surgery with sustained AT-RvD1 treatment, the UMAP shows neutrophils scattered among separate islands, and the SPADE pseudotime analysis shows the neutrophils progressing along separate trajectories while the vehicle control shows neutrophils clustered along one main trajectory (Fig. 1H). These results indicate that neutrophil heterogeneity is induced by local inflammatory signals along with lipid signaling pathways modulated by sustained AT-RvD1 treatment in the tissue microenvironment.

Fig. 1 Visualization of neutrophil heterogeneity using UMAP and SPADE.

A traditional gating strategy in which live, single cells were analyzed using CD11b+ as a general myeloid cell marker, and Ly6G+ within the CD11b+ population as a neutrophil marker was used to visualize neutrophil heterogeneity in the myeloid cell compartment. UMAP (left) and SPADE (right) of CD11b+ cells analyzed from the blood (A to D) and tissue (E to H) with an overlay showing the location and frequency of Ly6G+ cells. With a limited set of markers, neutrophil populations are not constrained to a single homogeneous cluster on the UMAP, which would imply one population. Comparison of UMAP and corresponding SPADE-based pseudotime trajectory at day 1 after surgery with a vehicle control (A and E), day 1 after surgery with sustained AT-RvD1 treatment (B and F), day 3 after surgery with a vehicle control (C and G), and day 3 after surgery with sustained AT-RvD1 treatment (D and H). Each UMAP/SPADE was created with the pooled data of four mice.

Pseudotime analysis identifies distinct pro-angiogenic neutrophil trajectories

We then took a closer look at the heterogeneity of neutrophils in the tissue by analyzing the pseudotime dendrograms generated by a SPADE analysis of CD11b+Ly6G+ neutrophils. To interpret and annotate the SPADE trees shown in Fig. 3, we created several versions of each tree, colored according to the median intensity of each measured marker (fig. S2, A to D). These colored trees were then used to identify the subpopulation of neutrophils represented on the different parts of the tree. This method has been shown to identify cell populations that are mischaracterized in manual gating strategies (18). We used selected angiogenic markers—VEGFR1, VEGFR2, and CXCR4—and integrin markers—CD11b and CD49d—to assign phenotype definitions to the SPADE dendrogram trajectories. At day 1, we identified a pro-angiogenic neutrophil subset (red stroke) expressing high levels of CD49d, VEGFR1, VEGFR2, and CXCR4 and an aged neutrophil subset (yellow stroke) with a small size and expressing high levels of CD49d, CD11b, and CXCR4 (Fig. 2, A and B). At day 3, the aged neutrophil subset is not present in the neutrophil population, but the pro-angiogenic subset persists (Fig. 2, C and D). We identify transitioning/intermediate neutrophil subsets (blue, purple, and orange strokes) within the CD49d+ population that lack the expression of one or more of the angiogenic markers (Fig. 2, I to VI). The relative marker expressions of each node on the SPADE tree can also be visualized on the corresponding heatmap showing every node’s unique surface marker expression level above or below the means of the population (fig. S3, A to D). This validates the phenotypes assigned to the SPADE tree as being unique subsets. The ability to separate the transitioning phenotypes from the pro-angiogenic population allows a more specific insight into the effects of AT-RvD1 on the neutrophil population.

Fig. 2 Analysis of neutrophil pseudotime trajectories using SPADE culminates in a unique marker expression profile of a putative pro-angiogenic subpopulation.

Pseudotime analysis generated SPADE trees consisting of CD11b+Ly6G+ neutrophils from tissue at day 1 after surgery with vehicle control (A), day 1 after surgery with sustained AT-RvD1 treatment (B), day 3 after surgery with a vehicle control (C), and day 3 after surgery with sustained AT-RvD1 treatment (D). Pseudotime trajectories show that distinct neutrophil phenotypes emerge. Day 1 after surgery shows the presence of an aged neutrophil phenotype (A and B). At day 3, the aged neutrophil subpopulation is not present, but the accumulation of pro-angiogenic neutrophils is notable in response to AT-RvD1 (D). The SPADE tree is annotated to distinguish the different subpopulations that emerge through pseudotime analysis. Gray nodes are CD49d classical neutrophils, and black nodes are CD49d+ with color strokes separating these cells into separate neutrophil phenotypes based on functional markers. The expression level, in the form of a heatmap overlay onto the nodes, of these functional markers for every subpopulation is displayed as a call-out box. The following colors and their expression level of a combination of functional markers are the strokes and markers that separate the neutrophil phenotypes: blue, CD49d+VEGFR2Lo(I)VEGFR1Lo(II)CXCR4Lo(III); purple, CD49d+VEGFR2Hi(IV)VEGFR1Lo(V)CXCR4Lo(VI); orange, CD49d+VEGFR2Lo(X)VEGFR1Hi(XI)CXCR4Hi(XII); red, CD49d+VEGFR2Hi(VII)VEGFR1Hi(VIII)CXCR4Hi(IX) pro-angiogenic neutrophil; and yellow, CD49d+CD11bHi(X)CXCR4Hi(XI)FSCLo(XII) aged neutrophil. Each SPADE tree was constructed with pooled data of n = 4.

AT-RvD1 selectively targets pro-angiogenic neutrophils in injured tissue

To assess the effect sustained AT-RvD1 has on the different neutrophil populations in an inflammatory environment, we quantified the cell frequency of each SPADE-identified neutrophil subset (Fig. 3, A to D). At 3 days after injury, there is an increase in nonclassical CD49d+ neutrophils in the local tissue environment in response to sustained AT-RvD1 treatment (Fig. 3A). The driving force behind this increase in CD49d+ neutrophils is a targeted increase in pro-angiogenic neutrophils, as that is the only subset affected by AT-RvD1 delivery (Fig. 3, B to D). We explored this targeted effect further by constructing days 1 and 3 tornado plots (Fig. 3, E and F). To do so, we generated a day 1 and day 3 SPADE tree that includes both cells from AT-RvD1–treated and vehicle control dorsal tissue (fig. S4, A and B). The difference of the frequency of cells from the treated tissue and the frequency of cells from the control tissue is calculated for each node. The greater the number resulting from this calculation, the more dominated the node is by cells from the AT-RvD1–treated tissue. The nodes are assigned a phenotype according to their expression of angiogenic and integrin surface markers (fig. S4, C and D), and the corresponding bar on the tornado plot is colored according to the phenotype. At day 1, there is no specific neutrophil phenotype that is dominated by treatment; however, at day 3, the pro-angiogenic neutrophil nodes are dominated by AT-RvD1 treatment (Fig. 3, E and F). In addition, at day 3, 20 of the 30 nodes found to be dominated by cells from AT-RvD1–treated tissue are a pro-angiogenic neutrophil cluster of cells (Fig. 3F). Thus, the pseudotime analysis of neutrophils using SPADE reveals a selective increase in pro-angiogenic neutrophils due to AT-RvD1–sustained release.

Fig. 3 AT-RvD1 treatment selectively increases the frequency of CD49d+VEGFR2HiVEGFR1HiCXCR4Hi pro-angiogenic neutrophils.

Day 1 (A) and day 3 (B) SPADE dendrograms containing both vehicle control and AT-RvD1–treated CD11b+Ly6G+ neutrophils. Every node on the SPADE tree is assigned a neutrophil phenotype and then whether the cells in that node are dominated by those that come from a treated sample or not are visualized in a tornado plot. The corresponding tornado plots for day 1 and day 3 are shown in (C) and (D), respectively. (C and D) The rank difference between AT-RvD1–treated and vehicle control samples for day 1 (C) and day 3 (D) was calculated and plotted for each node. The nodes in the tornado plot are colored according to their identified neutrophil phenotype from SPADE, showing the targeted effect of AT-RvD1 treatment on pro-angiogenic neutrophils. A difference (AT-RvD1 treated to vehicle control) under 0 reflects a decrease in that particular cell phenotype number with AT-RvD1 treatment, whereas a positive difference suggests that there is an increase in cell number with AT-RvD1 treatment. Frequency of the neutrophil populations identified using SPADE as a percentage of all CD11b+Ly6g+ neutrophils shows the driving force for an increase in frequency of CD49d+ neutrophils due to AT-RvD1 treatment is specifically the pro-angiogenic population (E and F). Statistical analyses were performed using a one-way analysis of variance (ANOVA) with Tukey’s multiple comparisons; *P < 0.05 and ****P < 0.0001; n = 4 animals per group.

Pro-angiogenic neutrophils preferentially localize near remodeling vessels

To assess the functional significance of the accumulation of pro-angiogenic neutrophils after AT-RvD1 treatment, we analyzed the vessel remodeling after injury through whole mount immunohistochemistry (Fig. 4A). Representative images of crops around the implant for different time points and treatments show that AT-RvD1 treatment has a vessel remodeling effect (Fig. 4B). We quantified the metrics of vessel remodeling by taking intravital bright-field images of the vascular bed at the time of surgery and 3 days after surgery (Fig. 4, C to E). In accordance with qualitative observations in the confocal images (Fig. 4, A and B), the quantitative vascular metrics show a significant increase in vessel tortuosity and arteriolar diameter because of AT-RvD1 treatment (Fig. 4, D and E). We performed a multiplex cytokine analysis to further investigate possible pathways AT-RvD1 treatment could be activating en route to a pro-angiogenic outcome (Fig. 4, F to I). At both day 1 and day 3, the AT-RvD1–treated and vehicle control samples separate along component 1 (Fig. 4, F and H). Variable correlation plots show a sweep in the cytokines that strongly correlates these cytokines with the separation of the samples on component 1 (Fig. 4, G and I). Furthermore, at both time points, VEGF, SDF-1α (stromal cell–derived factor 1α), and MCP-1 (monocyte chemoattractant protein-1) heavily contribute to principal component 1 and, thus, are key components contributing to the vascular remodeling observed with AT-RvD1 treatment (Fig. 4I and fig. S5).

Fig. 4 Delivery of AT-RvD1 promotes vascular remodeling.

(A) Image (×10) of the whole mount dorsal skin tissue. CD49d was the marker chosen to extrapolate pro-angiogenic neutrophil functionality in response to AT-RvD1 treatment. (B) Representative confocal images of neutrophil infiltration and remodeling vasculature (white circles indicate biomaterial implantation location). (E) Quantification of the vascular metric length density (C), tortuosity (D), and arteriolar diameter. PCA (principal components analysis) on cytokine concentration at day 1 (F and G) and day 3 (H and I) after surgery. PCA plot showing the vehicle control and AT-RvD1–treated samples on two components (F and H). The corresponding PCA loading plot follow the PCA plots for day 1 (G) and day 3 (I). Statistical analyses were performed using two-tailed t tests; *P < 0.05 and **P < 0.01; n = 4 to 6 animals per group. Scale bars, 500 μm for 10× confocal and 100 μm for 20× confocal.

Pro-angiogenic neutrophils have been implicated in playing a role in local vessel remodeling during inflammation (1, 22). Therefore, we sought to validate our finding from the SPADE analysis that AT-RvD1 results in an accumulation of pro-angiogenic neutrophils and gain further insights on the functional effects of AT-RvD1. We rendered the whole mount confocal images (Fig. 4, A and B) into surfaces using Imaris (Fig. 5A). This allows us to count the number of Ly6g+ neutrophils along with the number of Ly6g+ neutrophils also expressing CD49d while calculating the distance of these cells to the surrounding vessels. In line with the SPADE analysis showing an increase in CD49d+ neutrophils due to AT-RvD1 treatment, using immunohistochemistry, we see a significant increase in these cells due to AT-RvD1 treatment (Fig. 5, B and D). In addition to a higher number of the CD49d+ neutrophils, this analysis reveals that this subset of neutrophils specifically localizes closer to vessels when treated with AT-RvD1 compared with a vehicle control (Fig. 5, C and E). The finding that AT-RvD1 has targeted, direct effects on this pro-angiogenic subset of neutrophils is supported by the higher expression of its receptor—FPR-2 (formyl peptide receptor 2)—on CD49d+ neutrophils compared with CD49d neutrophils (Fig. 5F). In addition, we used SPADE to identify potential pro-angiogenic subsets of the general Ly6G myeloid cell pool at days 1 and 3 after injury (fig. S6, A to D). We identify different subpopulations of myeloid cells using VEGFR1, VEGFR2, and CXCR4 and show that AT-RvD1 has no significant effect on any of these subgroups (fig. S6, E and F).

Fig. 5 Recruitment and infiltration of CD49d+ neutrophils into the dorsal tissue.

CD49d was the marker chosen to extrapolate pro-angiogenic neutrophil functionality in response to AT-RvD1 treatment. (A) Imaris renderings from the confocal images of neutrophil infiltration into tissue at day 1 and day 3. (B to E) Quantification of Imaris renderings for total CD49d+ neutrophils and distance to CD31+ vessel at days 1 and 3. (F) Quantification of FPR-2 receptor expression on CD49d and CD49d+ neutrophils. Statistical analyses were conducted using two-tailed Mann-Whitney test and two-tailed t test; *P < 0.05 and ***P < 0.001; n > 100 cells, across three to four animals per group. Scale bars, 30 μm for Imaris.

AT-RvD1 treatment has a lasting immunoregenerative effect after skin transplantation

Having established the targeted effects of AT-RvD1 on pro-angiogenic neutrophils, we used a genetically identical donor-to-host tail skin tissue graft model to establish the rigor of our approach and assess whether the early proregenerative, AT-RvD1–mediated shift in the neutrophil population has effects on cell populations downstream in the inflammatory cascade. We developed PEG-MAL hydrogels capable of locally delivering therapeutic doses of AT-RvD1 in the first 24 hours (fig. S7) and injected the hydrogel into the wound immediately after skin transplantation. We collected tissue samples at 3, 7, and 14 days after injury from both vehicle control and AT-RvD1–treated groups for single-cell analysis. Using manual gating techniques, we selected CD11b+CD64+Mertk+ macrophages to analyze using SPADE. We pooled all the tissue samples (across all days and treatments) and constructed a macrophage-specific SPADE dendrogram. Using SPADE, we are able to identify a CD206loLy6Chi M1 macrophage group, a CD206hiLy6clo M2 macrophage group, and a CD206hiLy6chi M1/M2 hybrid macrophage group (Fig. 6A). With this SPADE dendrogram, we are able to isolate specific tissue samples and identify how AT-RvD1 treatment shifts the frequency of macrophages over time (Fig. 6, B to G). At day 3, both vehicle control and AT-RvD1–treated macrophage populations show a similar composition of macrophage subpopulations (Fig. 6, B and C). However, from day 3 to day 7, we see that while in the vehicle control case, the macrophages shift further toward an M1 phenotype, and in the AT-RvD1 treatment case, the macrophages shift toward an M2 phenotype (Fig. 6, D and E). This shift continues from day 7 to day 14 as the vehicle control macrophages settle into the profibrotic hybrid phenotype, while the AT-RvD1–treated macrophages have a predominant M2 phenotype (Fig. 6, F and G). By quantifying these shifts in the macrophage population visualized in SPADE, we confirm that there is a decrease in M1 macrophages at day 7 and an increase in M2 macrophages at day 14 with AT-RvD1 treatment (Fig. 6, H to J).

Fig. 6 AT-RvD1 promotes a proregenerative shift of the macrophage population.

(A) CD64+MERTK+ macrophages pooled from all samples from both treatments and three time points (days 3, 7, and 14) were used to construct a SPADE dendrogram. This SPADE dendrogram is annotated according to the Ly6C and CD206 surface markers expression. The red annotation corresponds to a hybrid phenotype with high expression of both Ly6C and CD206 (I and II), the blue annotation corresponds to an M2 phenotype with low expression of Ly6C and high expression of CD206 (III and IV), and the green annotation corresponds to an M1 phenotype with high expression of Ly6C and low expression of CD206 (V and VI). (B to G) Cells corresponding to a specific time point and treatment combination are isolated and mapped out on the SPADE dendrogram to analyze where these cells are highly populated (gray nodes) and where they are less densely populated (white nodes). The movement of these cells across time is depicted with dashed arrows that show where the cells are visualized to move (B to E). (H to J) The visualization method shows an AT-RvD1–mediated shift in cell frequency toward an M2 phenotype that can also be seen in the corresponding quantifications. Statistical analyses were performed using a two-way ANOVA with Tukey’s multiple comparisons; *P < 0.05 and ***P < 0.001; n = 4 animals per group.

Using an identical approach, we manually selected CD11b+CD11c+ dendritic cells to analyze using SPADE. We pooled all the tissue samples (across all days and treatments) and constructed a dendritic cell–specific SPADE dendrogram. Using SPADE, we are able to identify a IL10hi tolerogenic dendritic cell group and a IL10lo classical dendritic cell group (fig. S8A). With this SPADE dendrogram, we are able to isolate specific tissue samples and identify how AT-RvD1 treatment shifts the frequency of dendritic cells over time (fig. S8, B to G). From day 3 to day 7, we see a shift in dendritic cells toward an IL10hi phenotype with the AT-RvD1 treatment but see an opposite shift in the vehicle control case (fig. S8, B to E). By quantifying this shift in the dendritic cell population visualized in SPADE, we confirm that there is an increase in IL10hi tolerogenic dendritic cells at day 7 with AT-RvD1 treatment. Combined, these results establish the utility of SPADE and indicate that there is a lasting proregenerative effect on macrophages and dendritic cells after AT-RvD1 treatment.

To further demonstrate the effectiveness of biomaterial-mediated AT-RvD1 treatment as a regenerative immunotherapy, we performed a multiplex cytokine analysis and investigated the wound healing and integration of skin tissue grafts using macroscopic gross images of the grafted syngraft over time. As expected, both control and AT-RvD1 groups show no signs of transplant rejection. However, at day 10, the AT-RvD1 hydrogel treatment group appears to be closer to full wound closure than the vehicle control. After the principal components analysis of the multiplex cytokine panel, at days 3, 7, and 10, AT-RvD1–treated and vehicle control samples separate along component 1 (fig. S9, A to C). Similar to the cytokine analysis performed in the dorsal skinfold window chamber model (Fig. 4), VEGF and SDF-1α heavily contribute to the separation (fig. S9, D to F). These cytokines also have a significantly higher concentration with AT-RvD1 treatment compared with control (fig. S10, E and P).


The immune system is able to achieve a wide range of responses in health and disease due to the variety and heterogeneity of specialized immune cells. However, classical analytical techniques greatly underestimate the functional diversity of the immune system by providing an imperfect systems-level view of bulk immune cell populations. This is a major hurdle in the regenerative medicine field, where there is an increasing focus on the recruitment or polarization of specific immune cell subsets toward a therapeutic outcome (2325). Recently, new analytical and visualization techniques have been developed that allow scientists to gain new insights into immune system diversity (26). As these techniques are applied in the field of regenerative medicine, significant mechanistic discoveries can be made that aid in developing the most effective therapeutic strategies. However, there are complex environmental cues and timing factors that contribute to high cell-to-cell variation in the expression of surface markers, which presents a unique challenge. Such high variation complicates the analysis of these experiments since a population of cells captured at the same time includes many distinct intermediate differentiation states, and considering only their average properties masks trends occurring across individual cells (27, 28). We overcome this challenge using SPADE, a dimensionality reduction and pseudotime analysis technique that allows us to explore neutrophil heterogeneity in response to a chronic inflammatory environment in an in vivo dorsal skin wound model and explore macrophage and dendritic cell response in a skin transplant model to the sustained delivery of AT-RvD1, a proresolving lipid mediator.

Neutrophils, until recently, were considered to be a relatively homogeneous cell population with a set of conserved functions—phagocytosis, degranulation, and release of nuclear material. However, current research in the field of neutrophil biology has revealed that neutrophils have a much more diverse repertoire of functional responses (1, 29). This emergence of functional neutrophil phenotypes opens up the possibility to harness the therapeutic potential of neutrophil subsets. Using UMAP and SPADE to investigate neutrophil heterogeneity, we found that in the presence of an inflammatory injury, there is a much larger degree of neutrophil heterogeneity in the injured tissue than in circulation. This result speaks to the long-standing question of whether defined phenotypes are programmed in the bone marrow or whether mature neutrophils in circulation can be reprogrammed by local external stimuli (1). Thus, our results show that the local inflammatory microenvironment provides cues that induce the neutrophils to acquire specialized functions, which is depicted as different subpopulations spread across distinct islands on the UMAP and among different trajectories and nodes on the SPADE dendrograms. We locally deliver AT-RvD1 into this microenvironment to skew this neutrophil reprogramming to a more proregenerative state. More specifically, because we and others have shown the pro-angiogenic effects of AT-RvD1, we sought to use AT-RvD1 to locally induce the pro-angiogenic polarization of neutrophils (12, 13).

Pro-angiogenic neutrophils are a distinct neutrophil subset recently found to display tropism for angiogenic foci and contribute to rapid and early vascular growth (8). Impaired tissue vascularization is associated with chronic inflammation, and early angiogenesis helps restore appropriate tissue structure and function. Resolvins, a class of specialized proresolving lipid mediators, have been extensively characterized in their ability to inhibit neutrophil infiltration and promote their timely apoptosis (11, 30). However, the use of such lipid signaling pathways to target specific neutrophil subsets for regenerative immunotherapy has yet to be explored. We use SPADE to track the tissue neutrophil lineage in pseudotime, which allows us to look at the progression of pro-angiogenic marker expressions that results in transitioning neutrophil subsets that culminate in a pro-angiogenic neutrophil subpopulation. With this, we are able to analyze the effects that sustained AT-RvD1 delivery on each subpopulation. We found that by day 3 after AT-RvD1 delivery, there is an overwhelming increase in the pro-angiogenic subset of neutrophils, while there is no effect on the transitioning nonclassical neutrophils identified. SPADE identifies a separate trajectory for aged neutrophils at day 1 that is not present at day 3 after injury. This could be explained by the idea that some neutrophil subtypes may derive from aging neutrophils migrating into tissues given their activated and highly reactive state (31). The specific accumulation of pro-angiogenic neutrophils is a novel mechanism of action for this proresolving lipid mediator and is the first small-molecule approach that is able to selectively target a therapeutically functional neutrophil subpopulation. Our results suggest that the local delivery of AT-RvD1 results in the local reprogramming of neutrophils that have migrated into the tissue; however, the mechanism behind this reprogramming by an alteration of lipid signaling pathways requires further studies.

Furthermore, using whole mount immunohistochemistry of the injured dorsal skin tissue, we see an increase in Ly6G+CD49d+ cells. We are able to use this simplified marker staining to imply the function and behavior of the pro-angiogenic subset of neutrophils since we show in our previous results that the pro-angiogenic neutrophil subset is the overwhelming driving force behind any changes in this CD49d+ parent neutrophil population. Not only does this accumulation of Ly6G+CD49d+ cells validate our cytometry results, but it also allows us to gain unique insights into the functional changes that AT-RvD1 delivery is having on pro-angiogenic neutrophils. Our results show that these cells are positioned in closer proximity to the remodeling vasculature of the tissue in response to AT-RvD1. This accumulation and activation of pro-angiogenic neutrophils with sustained AT-RvD1 delivery coincide with an increase in arteriolar diameter and vessel tortuosity obtained from intravital bright-field imaging. Since the AT-RvD1 receptor FPR-2 is shown to be enriched on CD49d+ neutrophils, the delivery of this small molecule could be acting directly on this neutrophil subset. Linking the vessel remodeling that occurs with AT-RvD1 delivery directly to this neutrophil subset is still in need of further studies and is an exciting avenue for the field.

Evidence suggests that the early stages, right after the onset of inflammation, initiate an active, controlled process of inflammation resolution (32). We use an AT-RvD1–loaded PEG-MAL hydrogel in a syngeneic skin graft model not only to further establish the rigor of our approach but also to show that modulating the early neutrophil response, as shown in the dorsal skinfold model, leads to a cascade of proregenerative immune cell accumulation. Specifically, our results show that sustained AT-RvD1 delivery in the first 24 hours after skin transplant is able to increase the accumulation of M2 macrophages and IL10+ tolerogenic dendritic cells at day 7 after transplantation. Recent studies have established a connection between neutrophil activity and macrophage or dendritic cell reprogramming (33, 34). Thus, a strategy in which modulating the neutrophil response to injury leads to a proregenerative shift in downstream immune cell populations is an exciting avenue to explore. Whether AT-RvD1 delivery causes a change in the production of proresolution lipid mediators to have this downstream effect is also in need of further studies. Overall, the ability for AT-RvD1 delivered with a hydrogel to increase the wound closure after skin transplant further exhibits the therapeutic potential of this approach.

This study uses a novel method in which a nonlinear dimensionality reduction technique, UMAP, is used to visualize neutrophil heterogeneity on a two-dimensional map, and SPADE is used to make cellular hierarchy inferences among subpopulations of similar neutrophils. This method allowed us to uncover a previously unknown mechanism in which the sustained delivery of the proresolving lipid mediator AT-RvD1 leads to a targeted accumulation of pro-angiogenic neutrophils, increased early vessel remodeling, and enhanced wound closure after a skin transplant. Furthermore, as the targeting of key immune subsets becomes vital to regenerative immunotherapies, our studies can be used as a model to explore the subtle responses of different cellular populations to various therapeutic perturbations.


Fabrication of polymeric thin films

PLGA thin films were made, as previously described. Briefly, PLGA (50:50; DLG 5E, Evonik Industries) was solubilized in dichloromethane using a sonicator at 37°C until dissolved. Ten micrograms (100 μl) of AT-RvD1 (Cayman Chemical; solution in ethanol) was added to make AT-RvD1–loaded films. Solutions were cast in Teflon molds and stored at −20°C until full organic solvent evaporation was observed. Films were then lyophilized for 24 hours. A 1-mm biopsy punch was used to produce films used in studies.

Hydrogel fabrication

Four-arm PEG (10-kDa molecular weight) end functionalized with maleimide (>95% purity; Laysan Bio) at 4.5% (w/v) was used for all hydrogel formulations. PEG macromers were functionalized with RGD peptide (GRGDSPC), cross-linked with the cysteine-flanked peptide VPM (GCRDVPMSMRGGDRCG) (AAPPTec) in 0.5 M MES buffer (pH 5.5). The final concentration of RGD was 1.0 mM. Gels were also loaded with AT-RvD1 (4 μg/ml) (Cayman Chemical). The cross-linker concentration was based on the concentration of nonreacted maleimide groups remaining on PEG macromers. For hydrogels used in animal studies, all components were filtered through a spin column after pH measurements and kept under sterile conditions until injection into the animals.

Dorsal skin fold window chamber surgery

Animal experiments were performed using sterile techniques in accordance with an approved protocol from the Georgia Institute of Technology Institutional Animal Care and Use Committee. Male C57BL/6 mice (the Jackson laboratory) aged 6 to 12 weeks were anesthetized by inhaled isofluorane and surgically fitted with sterile dorsal skinfold window chambers (APJ Trading Co.), as previously described. Briefly, the dorsal skin was shaved, depiliated, and sterilized via three washes with 70% ethanol and chlorhexidine. The dorsal skin was drawn away from the back of the mouse, and one side of the titanium frame was attached to the underside of the skin. Sterile surgical microscissors were then used to expose the microvasculature through the removal of the epidermis and dermis in a 12-mm-diameter circle. Mice were implanted with two of the same films (either empty PLGA vehicle film or AT-RvD1–loaded PLGA film) placed on opposite sides of the window chamber. Before implantation, the films were washed in 70% ethanol for 30 s, followed by washing with sterile Ringer’s solution for 30 s. Exposed tissue was then sealed with a sterile glass coverslip. Mice were administered with sustained-released buprenorphine intraperitoneally (0.1 to 0.2 mg/kg) and allowed to recover in heated cages. All mice received standard laboratory diet and water ad libitum throughout the course of the experiment.

Skin transplant graft experiment

Female C57BL/6 mice aged 5 to 7 weeks (Charles River) were transplanted with full-thickness tail skin grafts from age- and sex-matched Balb/c or C57BL/6 mice. Briefly, donor mice were euthanized under terminal anesthesia via cervical dislocation. Donors were placed tail first into a sterile decapicone with the tail exposed. The tail skin was sterilized via three washes with 70% ethanol and chlorhexidine. The tail skin was removed with surgical scissors and kept in ice-cold sterile saline until use. Recipient mice were anesthetized with 50 μl of ketamine-xylazine-saline cocktail (ratio of 5:1:4) and shaved and sterilized via three washes with 70% ethanol and chlorhexidine. A full-thickness defect was cut into the dorsal skin of the recipient mouse. Hydrogels were mixed and loaded into a syringe and injected in depots in the defect area. Skin grafts were cut from donor tails to fit into the defect area. Wounds were bandaged, and mice were administered with sustained-released buprenorphine intraperitoneally (0.1 to 0.2 mg/kg) and allowed to recover in heated cages. Bandages covering grafts were removed 7 days after graft placement. All mice received standard laboratory diet and water ad libitum throughout the course of the experiment.

Tissue harvest and flow cytometry

To collect samples for flow cytometry analysis, mice were euthanized via CO2 asphyxiation. The dorsal tissue was excised and digested with collagenase type 1A (1 mg/ml; Sigma-Aldrich) at 37°C for 30 min and further separated with a cell strainer to create a single-cell suspension. Single-cell suspensions were stained for flow cytometry analysis using standard methods and analyzed on a FACSAria III flow cytometer (BD Biosciences). Dead cells were identified using the Zombie Green fixable viability stain (BioLegend). The antibodies used for identifying neutrophils and neutrophil subsets were as follows: BV421-conjugated anti-CD11b (BioLegend), BV510-conjugated anti-Ly6C (BioLegend), APC-Cy7–conjugated anti-Ly6G (BioLegend), PE-conjugated anti-CD49d (BioLegend), PerCP-Cy5.5–conjugated anti-CXCR4 (BioLegend), APC-conjugated anti-VEGFR1 (BioLegend), and PE-Cy7–conjugated anti-VEGFR2 (BioLegend). Staining using BV (Brilliant Violet) dyes was performed in the presence of Brilliant Stain Buffer (BD Biosciences). Positivity was determined by gating on fluorescence minus one controls.

High-dimensional analysis of flow cytometry data

UMAP is a nonlinear dimensionality reduction algorithm. UMAP is able to embed high-dimensional data into a space of two or three dimensions. Cells are visualized in a scatter plot, where points that are closer together can be considered more similar. Before UMAP dimensional reduction, each flow cytometry sample was manually pregated to select live, single, CD11b+ myeloid cells. The gated FCS (fluorescence correlation spectroscopy) files were imported into Python 3.7 using fcsparser ( and Pandas 2.5. Each channel except for FSC (forward scatter) and SSC (side scatter) was normalized by applying arcsinh transformation with a cofactor of 150 to transform fluorescence data into a fold-channel scale. A composite UMAP plot that used data points from all samples (eight total samples: four vehicle and four AT-RvD1) was performed on MATLAB ( Composite UMAPs were generated for day 1 and day 3. Each cell was then identified as a neutrophil by overlaying the manually pregated live, single-cell, CD11b+, Ly6G+ neutrophils onto the UMAP projection using MATLAB.

SPADE is a visualization tool that was designed to map heterogeneous single-cell populations into two dimensions on the basis of similarities across defined markers. SPADE creates a tree where nodes represent clusters of cells with similar marker expression. The size and color of each node are relative to the number of cells present and the median marker expression. SPADE was performed through MATLAB, and the source code is available at SPADE automatically constructs the tree by performing density-dependent downsampling, agglomerative clustering, linking clusters with a minimum spanning-tree algorithm, and upsampling based on user input. The SPADE tree was generated by exporting uncompensated pregated live, single-cell, CD11b+ myeloid cells or pregated live, single-cell, CD11b+, Ly6G+ neutrophils. The markers used to build the SPADE tree were SSC, FSC, CD11b, Ly6G, CD49d, CXCR4, VEGFR1, and VEGFR2. The following SPADE parameters were used: apply compensation matrix in FCS header, Arcsinh transformation with cofactor of 150, neighborhood size of 5, local density approximation factor of 1.5, maximum allowable cells in pooled downsampled data of 50,000, target density of 20,000 cells remaining, and number of desired clusters of 100.

Tissue whole mount immunohistochemistry and confocal imaging

Following euthanasia, mouse vasculature was perfused with warm saline and then with 4% paraformaldehyde until tissues were fixed. The dorsal tissue was excised and permeabilized overnight at 4°C with 0.2% saponin. The tissues were blocked overnight in 10% mouse serum at 4°C. Tissues were incubated at 4°C overnight in a staining solution containing 0.1% saponin, 5% mouse serum, 0.5% fatty acid–free bovine serum albumin, and the following fluorescently conjugated antibodies: Alexa Fluor 594 anti-CD31 antibody (1:100 dilution; BioLegend) for blood vessel visualization, Alexa Fluor 488 anti-Ly6G (1:200 dilution; BioLegend) for visualization of neutrophils, and Alexa Fluor 647 anti-CD49d (1:200 dilution; BioLegend) for visualization of pro-angiogenic neutrophils. Tissues were washed four times for 30 min with 0.2% saponin and once with PBS (phosphate-buffered saline) and then mounted in 50/50 glycerol/PBS. Mounted samples were imaged on a Zeiss LSM 710 NLO confocal microscope. Crops of 332 × 332 μm at ×20 magnification in the peri-implant area were taken for image analysis in Imaris (Bitplane). Images were then blinded and rendered in Imaris by a third party. Cells expressing Ly6G or double positive for Ly6G and CD49d were identified in Imaris using the surface tool. Surfaces were identified by smoothing with a 1-μm grain size and automatic thresholding on absolute intensity. Touching objects were split using a seed points diameter of 10 μm. CD31+ vasculature was identified in Imaris using the same surface method as described above, also applying a 1-μm grain size, but instead manually selecting the threshold value optimized for each image and manually applying the volume filter to remove small debris. Distance to vasculature calculations between cells and CD31+ vasculature was made by applying a distance transformation to the CD31+ surface and recording the median position of each cell surface relative to CD31+ vessels.

Vascular metrics

Vascular metrics were analyzed, as described previously (13). Briefly, the mouse was mounted on a microscope stage and imaged noninvasively at ×5 magnification on a Zeiss imager D2 microscope with AxioCam MRC 5 color digital camera (Zeiss). Images were acquired on day 0 immediately following film implantation and again on day 3. Regions of interest (ROI) measuring 2000 × 2000 pixel2 were traced around the implant for further analysis. Vessels within these ROI measurements were traced, and the total vessel length per unit area was quantified via ImageJ. Arteriolar diameter measurements were measured via ImageJ, and day 3 diameters were normalized to day 0. Vascular tortuosity measurements were made within the ROIs by measuring the distance metric—the path length of a meandering curve divided by the linear distance between end points in ImageJ.

Cytokine measurements

For cytokine measurements, 4-mm biopsy punches of tissue centered on each biomaterial implant were harvested after euthanasia. Tissue biopsy punches were combined for each animal, digested for 30 min at 37°C in collagenase type 1A (1 mg/ml), and disaggregated through a cell strainer. Protein was isolated from the single-cell suspension in RIPA (radioimmunoprecipitation assay) buffer containing Halt Protease and Phosphatase Inhibitor Cocktail (diluted to 1×; Thermo Fisher Scientific) for 45 min on ice. Following cell lysis, total protein was obtained by centrifugation for 15 min at 14,000g and 4°C. To determine the total protein concentration in each sample, a bicinchoninic acid assay (BCA assay) was carried out using a Pierce BCA protein assay kit (Thermo Fisher Scientific) according to kit instructions. Cytokine measurements were made using the Mouse Magnetic Luminex Screening Assay (catalog number LXSAMSM, R&D Systems) according to kit instructions. Kit analytes included CXCL12/SDF-1α, FGF-β, MMP-9, and VEGF. Cytokine results were normalized to average vehicle expression at day 1 and visualized via heatmap.

Statistical analysis

All statistical analyses were performed using GraphPad Prism version 7.0 (La Jolla, CA). Results are presented as means ± SEM. For pairwise comparisons, unpaired two-tailed t tests with Welch’s correction, if variance was significantly different, were used. For multiple comparisons, one-way analysis of variance (ANOVA) was used with Tukey’s multiple comparisons when relevant. For analysis of cellular distance to vasculature, data reflect cells counted from three ROIs acquired across three to four animals per group, and statistical comparisons were made using a two-tailed Mann-Whitney test.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the core facilities at the Parker H. Petit Institute for Bioengineering and Bioscience at the Georgia Institute of Technology for the use of their shared equipment, services, and expertise. Funding: Research reported in this publication was supported by the NIH under award numbers R01AR056445 and R56AR071708, the NIH Cell and Tissue Engineering training grant under award number GM008433, the National Science Foundation under award number CCF1552784, the National Science Foundation Graduate Research Fellowship under award number DGE-1650044, and the Immuno-Engineering Seed Grant. Author contributions: T.C.T., M.C.P.S., L.A.H., E.A.B., P.Q., and G.A.K. designed the research, analyzed the data, and wrote the paper. T.C.T., M.C.P.S., L.A.H., W.Y.Y., S.V., H.S.L., F.S.P., and Q.D.M. performed the research, analyzed the data, and reviewed the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article