Single-cell Transcriptional Landscape of Temporal Neutrophil Response to Burn Wound in Larval Zebrafish Free

Neutrophils are the first responders of the innate immune response to tissue damage and infection. As the most abundant circulating leukocyte, neutrophils are important drivers of these early inflammatory responses and can also play a role in chronic inflammation and tissue damage. Complex injuries, such as burns, lead to robust innate immune inflammation and damaging tissue responses. Previous work has identified neutrophil activation in burn patients; however, the role of neutrophils is complicated and likely context-dependent (1, 2). In recent years, single-cell profiling of neutrophils in normal development and disease models has identified diverse neutrophil populations at the molecular level (3–6). Most studies have classified neutrophils based on their maturation status (7). Some studies further identified specific subsets of neutrophils that are positively or negatively associated with disease progression. One recent study focused on neutrophils in the early stage of severe burns and found human neutrophils matching the mouse G_3/4/5a–c subsets in both burned and healthy donors (5, 8). These studies showcase the utility of single-cell transcriptomics to facilitate our understanding of neutrophil involvement in tissue damage contexts.

Although direct profiling of neutrophils from human tissues is feasible, there are limitations to what can be learned from human burns. Neutrophils are often isolated either from local burn wound tissue (eschar) or peripheral blood late after injury, which does not reflect the early state of neutrophils in tissues or provide information about neutrophil progenitors in the bone marrow (9). Zebrafish enable whole organism profiling of neutrophils after burns with unbiased sampling of neutrophils at all tissue locations including the hematopoietic tissue. Previous studies have demonstrated that zebrafish provide a robust model system to study thermal injury that enables the imaging of the temporal and spatial dynamics of the innate immune response to injury (10, 11). Recent single-cell RNA sequencing (scRNA-seq) also confirmed conserved maturation stages and markers of neutrophil subsets in zebrafish and humans (12, 13), justifying the use of zebrafish as a model in studying human neutrophil responses to thermal injury. Using this model, we previously found that both human burn and zebrafish burn tissue showed an increase in IL-6 production (11). In zebrafish, we found that the IL-6 receptor was important for neutrophil but not macrophage recruitment to burn wounds; however, the mechanism for these effects remains unclear.

In this study, we used transcriptional profiling at single-cell resolution of neutrophils and macrophages in larval zebrafish to compare gene expression during homeostasis and following a burn wound. We describe the molecular diversity within neutrophils and macrophages during the early innate immune response to burn. By comparing burned and unwounded conditions, we show subset-level functional changes and subset-to-subset transitions. Furthermore, we present neutrophil–macrophage communications at the subset level and in specific ligand-receptor pairs, including IL-6 to IL-6R signaling. Finally, we identified myeloid subset gene signatures in zebrafish that correspond to burn severity in human patients.

Adult and larval zebrafish are maintained following the animal protocol M005405-R02-A01 approved by the University of Wisconsin–Madison Institutional Animal Care and Use Committee. Embryos were collected from adult Tg(mpx:Dendra2, mpeg1.1:mCherry) fish using natural breeding procedures and kept in E3 medium with 1% methylene blue inside a 28.5°C incubator. At 3 d postfertilization (dpf), larvae were anesthetized with tricaine methanesulfonate (MS222, 200 mg/l) and wounded at the tail fin region using a handheld cauterizer (10). As the cauterizer wire heats up the E3 medium at close vicinity to the tail fin, the affected tissue area changes its transparency and shows as a wave spreading toward the body. The degree of burn wounding is controlled by turning off the cauterizer right before the wave of transparency change hits the tip of the notochord region. After wounding, larvae were washed in E3 medium without MS222 and sent back to the incubator until the time of sample collection.

At 6, 24, and 48 h postburn (hpb), burned and unwounded larvae at matching developmental stages (3, 4, and 5 dpf) were transferred to a 35-mm dish containing calcium-free PBS for 15 min and then anesthetized with tricaine methanesulfonate (MS222, 200 mg/l). A total of 150 fish were used for each time point by condition. Each dish of fish was digested with 2 ml of digestion solution (0.25% trypsin, 1 mM EDTA in PBS) at 28.5°C for 90 min with gentle pipetting every 10 min. Digestion was stopped by adding 200 µl of digestion stop solution (1 mM CaCl₂, 100% FBS). Dissociated cells were filtered with 40-µm cell strainers and centrifuged for 3 min at 3000 rpm at 4°C. Cell pellets were resuspended in PBS with 10% FBS plus DAPI (1 mg/ml).

FACS was performed at University of Wisconsin Carbone Cancer Center Flow Lab using a BD FACSAria. Targeted cells passed gating for cells versus debris, singles versus doublets, and live versus dead and were high in either the 488-nm or 561-nm channel. Both dendra⁺ (488-nm channel) and mCherry⁺ (561-nm channel) cells were sorted into PBS with 10% FBS. After sorting, cells were directly used for library construction in a Chromium controller at the Gene Expression Center of the University of Wisconsin–Madison Biotechnology Center (RRID: SCR_017757). Samples were processed with Chromium single-cell gene expression solution 3′ v2 following the manufacturer’s instructions (10x Genomics).

Libraries were sequenced at the DNA Sequencing Core of the University of Wisconsin–Madison Biotechnology Center (RRID: SCR_017759) on NovaSeq 6000 system in S1 flow cells with read lengths of 29 bp + 90 bp (read1 + read2). Raw sequencing reads were processed by cellranger count in Cell Ranger (v6.1.2, 10x Genomics) using GRCz11 zebrafish genome with an improved gene annotation (14). Filtered cell-by-gene matrices generated were used for downstream analysis in R with Seurat (v4.3.0) (15). We filtered out cells with <200 genes detected or a mitochondrial gene ratio >20%. Doublets were removed using DoubletFinder with an estimated ratio of 2.3% for T1 samples (6 hpb, 3 dpf), 4.6% for 24 hpb, 6.1% for 4 dpf, 6.9% for 5 dpf, and 7.6% for 48 hpb (16).

Individual Seurat objects were log-normalized and merged as one data object for linear dimensional reduction through obtaining the top 2000 highly variable genes (FindVariableFeatures), data scaling and centering (ScaleData), and performing principal component (PC) analysis (PCA; RunPCA). To account for batch effects, the Harmony batch correction method was employed, using the individual sample label as the batch variable (17). Nonlinear dimensional reduction was performed using Harmony-corrected top 100 PCs to generate uniform manifold approximation and projection (UMAP) axes (RunUMAP). Major cell types were identified by Louvain graph clustering at a resolution of 0.1 on the nearest neighbors graph also constructed by corrected PCs (FindNeighbors and FindClusters). Neutrophils, macrophages, and cycling myeloid cells were annotated based on the top differentially expressed genes identified by a Wilcoxon rank-sum test (FindAllMarkers with log₂ fold change threshold of 0.25 and positive only) and subsetted for subclustering.

Isolated myeloid cells were split back into a list of Seurat objects by their original sample identity and normalized through SCTransform. The list was integrated based on the SCT assay (name for the function SCTransform) and then through dimensional reduction by PCA and UMAP with the top 30 PCs. Clustering was performed with the top 30 PCs followed by differential expression analysis as described above.

On the human data side, we received the processed Seurat object provided by the authors from Montaldo et al. (6) and downsampled to keep 2654 cells in each neutrophil maturation stages (precursor, early immature, immature, and mature) (6). On our zebrafish data side, we converted zebrafish gene symbols to their corresponding human orthologs by using the Drosophila RNAi Screening Center Integrative Ortholog Prediction Tool (DIOPT) scoring application programming interface with zebrafish as the source organism and humans as the target organism (18). Only genes with a DIOPT score >6 were kept, ensuring confidence in ortholog calling. For genes with many-to-one orthologs, we used the most highly expressed gene from the zebrafish dataset. All label transfer procedures were conducted with tools in Seurat (v4.3.0). SCTransform was performed for both the received human neutrophil data and our own zebrafish neutrophil data. Anchors between the two datasets were identified by FindTransferAnchors using PCA with 30 PCs as reference reduction space, and SCT as the normalization approach. Label transfer from zebrafish to humans was performed by TransferData with the same dimensions used in anchor identification. The predicted zebrafish cluster name was added to the human neutrophil dataset as a new metadata.

Processed ZebraHub transcriptome data object covering all developmental time points (zf_atlas_full_v4_release.h5ad) was downloaded (https://zebrahub.ds.czbiohub.org/data) and converted into the H5Seurat format to be used in the R environment with Seurat (v4.3.0). Only 5 dpf cells were used for visualizing general clustering and localization of myeloid populations on UMAP axes. Myeloid cells from all developmental time points were extracted by CellSelector based on the myeloid marker expression pattern and used for coexpression visualization using FeaturePlot with blend = true.

Single-cell transcriptome data from the adult zebrafish tessellated lymphoid network were processed as previously described (19). Neutrophil was annotated based on the top differentially expressed genes identified by Wilcoxon rank-sum test (FindAllMarkers with log₂ fold change threshold of 0.25 and positive only) and subsetted for subclustering. After neutrophil-only dimensional reduction, they were used for coexpression visualization using FeaturePlot with blend = true.

We followed an analysis procedure adapted from Arora et al. (20) with details listed below.

For each cluster in neutrophils, we generated a pseudobulk profile by calculating the average expression level for each gene using the RNA assay. We performed FindVariableFeatures on neutrophils and filtered the average expression matrix to include only the top 2000 highly variable genes. Expression levels of these genes were used for Pearson correlation calculation and hierarchical clustering.

Pseudobulk profiles were generated similarly as above. We added one step before performing PCA to bring the average expression matrix back to log scale, which was missing in the original pipeline (20). These profiles were used as input for PCA with scale = true. Cell clusters were plotted onto two selected PC axes wrapped either by time point of sample collection to show burn-unwounded difference or by condition to show time course change. Top 500 genes ranked by positive or negative loading were extracted for Gene Ontology (GO) enrichment analysis using Metascape v3.5.20230501 (21).

Gene ranking for each cluster was generated by conducting Wilcoxon rank sum test using wilcoxauc from presto v1.0.0 comparing between neutrophil subset (I. Korsunsky, A. Nathan, N. Millard, and S. Raychaudhuri, manuscript posted on bioRxiv, DOI: 10.1101/653253). By-cluster ranking was then taken as input for gene set enrichment analysis of GO terms from clusterProfiler v4.6.2 with gene set sizes between 15 and 500 (22). GO terms were ranked by p value and with the top five of them selected for visualization.

For each GO term of interest, all genes associated with it in each myeloid subset were pulled out. Average expression level of each gene in each subset by condition and time was normalized by the mean value from all cells and summed up to represent the expression level of the corresponding GO term. These GO expression levels were converted to z scores for visualization.

Cell trajectories were constructed using Slingshot v2.6.0 (23) or Monocle3 v1.3.1 (24). From neutrophil-only visualization, we found outliers not clustered with the rest of the neutrophils and removed them from trajectory analysis. UMAP embedding and neutrophil subset identities were used as input for Slingshot calculation.

For Monocle, we performed independent analysis from filtered matrices generated by cellranger count. A cell dataset from each sample was combined into one data object for preprocessing, dimensional reduction, clustering, and graph identification using default parameters. Cells clustered with neutrophils labeled in the Seurat object were chosen for cell ordering with the root of trajectory manually selected at the Precursor-like_Neuts side.

RNA velocity was calculated by velocyto v0.17.17 (25) using the default parameters from the alignment generated by cellranger count for each sample collected by condition and time. The loom files generated by velocyto were used as input in complement to the Seurat object containing only neutrophils for velocity estimation with scVelo v0.2.5 in Python (26). In brief, Seurat object and loom files were merged and filtered to contain only the top 3000 genes with at least 20 counts shared between unspliced and spliced measurements. The resulting datasets were processed through pca and neighbors and then for velocity estimation using the stochastic mode.

The Seurat object containing all myeloid cells was converted to AnnData format then loaded in Python for network analysis with CellOracle v0.14.0 (27). The AnnData object was normalized by the total number of unique molecular identifier counts per cell and filtered to include only the top 3000 genes with at least one count. After filtering, the AnnData object was renormalized and used for perform_PCA followed by knn_imputation. For Gene Regulatory Network (GRN) construction, we used the zebrafish promoter base GRN (danRer11_CisBPv2_fpr2) provided by the package. Cluster-specific GRN was constructed with the top 2000 source-target links kept based on edge strength. All transcription factors (TFs) in the top 2000 links from all myeloid subsets were merged and clustered based on their degree of centrality. Top 30 TFs were included for visualization. The expression level of genes targeted by these TFs were extracted using the same approach as in Functional annotation for zebrafish myeloid subsets.

Cell–cell interaction analysis was performed using the CellChat v1.6.1 R package with customized database (28). All myeloid populations were included for analysis. Functions from CellChat were used to identify changes in interactions between populations across conditions at comparable timepoints. Visualizations were constructed from built-in visualization functions in CellChat or through export of CellChat data and visualization with ggplot2.

Burn microarray data from the Systems Biology Coagulopathy of Trauma (SYSCOT) study group was first stratified based on patient total burn surface area (TBSA) (29). Patients with 0–10% TBSA were classified as low TBSA, 11–20% were moderate TBSA, 21–30% were high TBSA, and 31%+ were very high TBSA. Using the human gene converted Seurat object constructed (see Human–zebrafish scRNA-seq label transfer), we identified zebrafish myeloid subset markers in the form of their human orthologues with FindAllMarkers implementing the model-based analysis of single-cell transcriptomics (MAST) hurdle model. Differentially expressed genes were further filtered based on a log₂ fold change of 0.5 and an adjusted p value of 0.01.

Microarray datasets were then evaluated at a per-sample level by averaging the top 50 gene expression values from each subset-specific signature. In the case where a subset gene signature contained fewer than 50 genes, all of the genes were used. The distribution of gene signature scores was visualized as raincloud plots with their means compared using a Wilcoxon rank-sum test as implemented in the ggpubr library function stat_compare_means. Null datasets were generated through 10000 permutation tests where TBSA labels were shuffled, and gene signature averages were calculated. The Spearman correlation between individual genes within the found signatures and the subset gene expression signature averages were calculated with the base R cor function with the method parameter equal to the Spearman correlation.

Peripheral blood was obtained without identifiers from the TSB BioBank, which maintains an IRB-approved protocol (UW HS IRB# 2016-0934) and honest broker system for collection and distribution of samples from consented subjects, meeting the definition of not human subjects research. Peripheral blood was collected from deidentified healthy donors and burn patients (Supplemental Table I). Burn samples (blood) were collected within 1 wk of hospital admission. All incubation and washing steps were performed at room temperature. Samples were incubated with RBC lysis buffer (BioLegend) diluted 1:10 in deionized H₂O for 10 min within 1 h of sample collection and then washed twice with cell staining buffer (BD Biosciences). For cell surface staining, human TruStain FcX (BioLegend) was added at 1:20 to block Fc receptors for 10 min, followed by incubation in Zombie NIR Live/Dead stain at 1:1000 (BioLegend) for 10 min. Protein-targeted staining was done with the following Abs in cell staining buffer: Alexa Fluor 700–conjugated anti-human CD16 at 1:400, PE/Fire 640–conjugated anti-human CD66b at 1:400, PerCP-conjugated anti-human CD3 at 1:200, PerCP-conjugated anti-human CD19 at 1:200, PerCP-conjugated anti-human CD56 at 1:400, PerCP-conjugated anti-human CD203c at 1:100 (BioLegend), BUV395-conjugated mouse anti-human CD45 at 1:100 (BD Biosciences), and PE-conjugated GPR84 at 1:200 (Alomone Labs). After staining for 30 min, cells were washed with cell staining buffer and resuspended in fixation buffer (BioLegend) for 15 min. After fixation, cells were washed and resuspended in 300 μl of cell staining buffer and filtered through a 70-μm filter. Flow cytometry data were collected on the Cytek Aurora flow cytometer. FCS files were exported and analyzed using FlowJo software v.10.9.0. Neutrophils were identified as CD3/CD19/CD56/CD203c⁻, CD45/CD66b/CD16⁺ cells with granulocyte scatter. These cells were then evaluated for the expression of GPR84. Fluorescence minus one of GPR84-PE was used as a control to quantify cells expressing GPR84.

All raw and processed sequencing data generated in this study have been submitted to the National Center for Biotechnology Information’s Gene Expression Omnibus (http://ncbi.nlm.nih.gov/geo/) under accession number GSE250610. Source codes and processed data are provided at https://github.com/huydinhlab/HouKhatriZebrafishBurn.

To characterize the molecular heterogeneity of myeloid cells in larval zebrafish, we used transgenic zebrafish Tg(mpx:Dendra2, mpeg1.1:mCherry) with fluorescently labeled neutrophils and macrophages. Cells were collected from both burned and unwounded larval zebrafish along a time course of 3 d (Fig. 1A; see Materials and Methods). We chose three time points to capture burn wound healing (6, 24, and 48 hpb) that cover the myeloid cell recruitment and inflammation resolution stages (10). In addition, we included the corresponding time points in normal development as unwounded controls (3, 4, and 5 dpf). We obtained transcriptomes of all cells at single-cell resolution using a droplet-based 3′ mRNA capturing approach from 10x Genomics. Initial processing revealed six known major cell populations, including neutrophils (mpx, lyz), macrophages (mpeg1.1, marco), lateral line cells (prox1a, f11r.1), epithelial cells (cldnb, krt4), erythrocytes (hbae1.1, hemgn), and photoreceptors (opn1mw1, gnb3b), and one cycling myeloid population with both neutrophil and macrophage signatures and markers labeling cell cycle (mki67, top2a). One identified population has not been previously characterized, with higher expression of cebpd and her6 (Supplemental Fig. 1A, 1B). We focused on neutrophils and macrophages (a total of 19,118 cells) for further clustering analysis that identified six macrophage subtypes and four neutrophil subtypes (Fig. 1B, Supplemental Figs. 1C, 2). Two of the six macrophage subsets showed enriched expression of cell cycle–related genes, including Cycling_Mac1, which had high level scoring for both S phase and G₂M phase markers (Supplemental Fig. 3).

Current molecular annotation for myeloid subsets in zebrafish is still limited, although there has been a recent report from adult zebrafish bone marrow (13). Therefore, we named clusters both by known annotations and their representative top expressed genes. For macrophage subsets (Fig. 1C), we examined the expression distribution of known M₁/M₂ subtype markers and could not match zebrafish subsets with this polarization nomenclature (Supplemental Fig. 1D). Based on gene signatures, we named the macrophage clusters as ccl34a.4_Macs, lgals3bpb_Mac1, lgals3bpb_Mac2, and cx43_Macs (Fig. 1B, Supplemental Figs. 1C, 2). Among these, ccl34a.4⁺ and lgals3bpb⁺ macrophages were also defined in a recent macrophage scRNA-seq dataset from adult zebrafish under homeostasis (30). For neutrophils, we performed cross-species label transferring between our zebrafish data and a comprehensive scRNA-seq mapping of human neutrophils, including four stages: precursors and early immature, immature, and mature neutrophils (6). We mapped zebrafish neutrophil subsets to human neutrophil stages using homologs of marker genes (see Materials and Methods; Supplemental Fig. 4A, 4B) and identified three out of four zebrafish neutrophil subsets that align well with human neutrophil stages, named as Precursor-like_Neuts, Immature_Neuts, and Mature_Neuts, respectively (Fig. 1D, Supplemental Fig. 4C). The il6r- and padi2-expressing neutrophil subset could not be mapped to any of the human neutrophil states, thus keeping its marker-based name as il6r_padi2_Neuts. We confirmed the existence of this population in scRNA-seq data from the ZebraHub database (Supplemental Fig. 5A–C), which covers up to 5 dpf, and in the adult zebrafish tessellated lymphoid network (Supplemental Fig. 5D, 5E) (19).

All subsets exist in both burned and unwounded conditions at each time point (Fig. 1E). This consistent distribution suggests that as early as 3 dpf, larval zebrafish have established heterogeneous populations of myeloid cells that are sustained through 5 dpf. In addition, burn wounding does not eliminate or induce any subset. In the burn healing process, Precursor-like_Neuts and lgals3bpb_Mac1 showed more apparent increases at an early time point compared with the unwounded control (Supplemental Fig. 1E). On the contrary, Mature_Neuts and cx43_Macs showed a more apparent decrease in burn. Taken together, we demonstrate the existence of heterogeneity within myeloid populations in both normal development and during burn wound healing in larval zebrafish. Hereafter, we focus more on dissecting the molecular features of neutrophil subsets and their interactions with macrophages.

We next aimed to understand the function of neutrophil subsets in normal development and burn wounds. First, to identify the global differences across subsets, conditions, and time points, we summarized expression in each subset and performed hierarchical clustering among subsets (Fig. 2A). Precursor-like_Neuts are grouped together, regardless of their condition or time point of collection. However, Immature_Neuts and il6r_padi2_Neuts at the earliest time point, either in the burn or unwounded condition, are grouped with the Mature_Neuts subset. To further dissect the variance between subsets, we performed PCA on the pseudobulk profiles. The top two PCs explained 24% of variance (Fig. 2B). By examining the grouping of subsets based on the first six PCs, we found that PC2 correlates with time progression and separated burned from unwounded at T1, while PC5 better separated different subsets (Supplemental Fig. 6).

Next, we asked what drives the separation on PC axes. We performed GO enrichment analysis on the top 500 genes with positive or negative correlation with each PC loading for PC annotation (21) (Supplemental Fig. 7). Genes with negative loading of PC2 showed enrichment for both active protein synthesis (dre04141 protein processing in endoplasmic reticulum, R-DRE-532668 N-glycan trimming) and degradation functions (GO:0030163, protein catabolic process; GO:0006517, protein deglycosylation). This suggests a higher turnover rate for neutrophils in T1 (3 dpf/6 hpb) than later time points and in the burn versus unwounded condition at T1. Genes associated with PC5 enrich for terms that reflect differences between subsets, with the positive direction associated with the Precursor-like_Neuts, showing cell cycle features (GO:0006260, DNA replication, dre03410 base excision repair, R-DRE-69206 G₁/S transition), whereas the negative direction showed more mature cell functions (GO:0045321, leukocyte activation; GO:0050900, leukocyte migration) and the signaling pathways that mediate them (GO:0007229, integrin-mediated signaling pathways; GO:0034097, response to cytokine; GO:0007179, TGF-β receptor signaling pathways). The data suggest that from Precursor-like_Neuts, to Immature_Neuts, and then to il6r_padi2_Neuts, the cell cycle–related signature drops while immune activity increases.

While PCA provides a global separation between subsets, it is not enough for comprehensive functional annotation for each. To increase the granularity of functional annotation, we performed gene set enrichment analysis for GO terms for all subsets (Fig. 2C) (22). Compared to overrepresentation in GO enrichment analysis, this approach allowed for both the existence and ranking of genes to be taken into consideration. Zebrafish GO terms were used as gene sets for rank testing. We further complemented the enrichment analysis with the average expression level of the genes associated with each GO term in each subset at individual time points (see Materials and Methods). We found subset-specific enrichment in certain biological processes, including translation-related terms in Precursor-like_Neuts, Immature_Neuts, and il6r_padi2_Neuts (GO:1990904, ribonucleoprotein complex; GO:0043043, peptide biosynthetic process). Precursor-like neutrophils are also highly enriched in steroid metabolism-related biological processes (GO:0008202, GO:0006694, GO:0016125), but without significant expression level changes between burned and control conditions. Both Immature_Neuts and il6r_padi2_Neuts enrich for the collagen catabolic process (GO:0030574). Genes associated with this term are also highly expressed in them, with highest levels at the 6 hpb condition (T1 in burn condition). This suggests potential early involvements of Immature_Neuts and il6r_padi2_Neuts in the wound healing processes. Immature_Neuts also uniquely enriched for neuropeptide signaling (GO:0007218), suggesting that zebrafish might have a subset of neutrophils that are regulated by neuropeptides similar to humans (31). Mature_Neuts enriched for GPCR activity (GO:0008227), protein autophosphorylation (GO:0046777), oxygen carrier activity (GO:0005344), Ig complex (GO:0019814), and cyclic nucleotide degradation (GO:0004114). Among these terms, GPCR activity and cyclin nucleotide degradation activity have their associated genes more highly expressed in the 48 hpb condition (T3 in burn). cAMP, as one type of cyclic nucleotides, is known to be at high levels in the urine of severe burn patients that gradually drops to a level closer to that in mild burn patients (32). It acts downstream of PGE₂ to suppress the activation and growth of T lymphocytes and granulocyte-macrophage progenitors (33, 34). Mature_Neuts in zebrafish could be involved in the attenuation of the cAMP-mediated immunosuppression effect at a later stage in burn wound healing. These results indicate that zebrafish neutrophil subsets play distinctive roles during normal development and could selectively enhance specific functions during burn wound healing.

From zebrafish–human label transfer (Fig. 1), we learned that zebrafish neutrophils follow the same path as human neutrophils through maturation. However, we also saw that il6r_padi2_Neuts do not match with human neutrophil states. What is the relationship between this unique subset and the conserved differentiation path? We first examined the state transition paths through two independent approaches: trajectory reconstruction with pseudo-time analysis (slingshot, monocle) and RNA velocity estimation (see Materials and Methods). Both slingshot- and monocle-inferred trajectories consistently predict a branched trajectory, with one branch matching the states of human neutrophils developing from Precursor-like_Neuts to Immature_Neuts and then to Mature_Neuts (Fig. 3A, Supplemental Fig. 8A). Il6r_padi2_Neuts are on the other branch of the trajectory that connects with Immature_Neuts, but with unclear directionalities. We examined the expression of mmp9, a marker associated with zebrafish neutrophil maturity, and found that il6r_padi2_Neuts have a higher level than do Mature_Neuts (Supplemental Fig. 8B) (13). However, human homologs of the top-ranking signatures of il6r_padi2_Neuts are mainly expressed at the early immature or immature stage in human neutrophils (Supplemental Fig. 8C).

To further understand the direction of the transitions among subsets, we performed RNA velocity analysis in neutrophils (Fig. 3B, Supplemental Fig. 8D; see Materials and Methods). RNA velocity incorporates the dynamics of the spliced and unspliced forms of mRNA across genes to infer the future state of a cell (25, 26). While the differentiation path matches the predicted directionality and remains stable, the path connecting il6r_padi2_Neuts with Immature_Neuts could be lost. This suggests that the il6r_padi2_Neuts might not be the terminal state of the path similar to the Mature_Neuts. When comparing between conditions, we found a strong velocity trend from Precursor_Neuts to Immature_Neuts and Mature_Neuts in burn, especially at T2 (Fig. 3B). This suggests a potential trauma-induced emergency hematopoiesis in zebrafish in response to burn wounding. Connecting this to the proportional changes shown in Fig. 1E, we see that the transition resulted in expansion of Precursor-like_Neuts and il6r_padi2_Neuts and a shrinkage of Mature_Neuts.

We next aimed to identify what factors drive the state transition of neutrophil subsets. We performed network analysis using CellOracle, which takes advantage of both the open chromatin signature and gene expression at single-cell resolution to build cell state–specific GRNs (27). We used the built-in promoter base GRN models in CellOracle to compensate for our lack of matching chromatin accessibility profiles. We constructed a TF to target gene link and compared the degree centrality for TFs in different subsets across conditions and time points (Fig. 3C, Supplemental Fig. 9A). The degree of centrality level reflects the number of downstream targets a TF might control and is used in this study to infer TF activity. Macrophages were included for elucidating cell type–specific TF activity or excluded to focus on the between-subset comparisons in neutrophils. We found that raraa and ved are more active in neutrophils, whereas foxg1c, sox9b, and rxraa are more active in macrophages (Supplemental Fig. 9A). As expected, between-subset distinctions are less pronounced compared with those from between-cell types. To further understand neutrophil subset-specific TFs across conditions and time, we calculated the average expression level of predicted downstream genes in neutrophil subsets for each candidate TF (Fig. 3C). In general, TFs with higher inferred activity in a subset also presented a relatively higher expression level of downstream target genes in that subset. This suggests that active TFs are not only controlling more genes, but also enhancing their expression level. Mature_Neuts and il6r_padi2_Neuts are the two branches in the trajectory analysis, so we compared the two and identified TFs with their activities associated with these populations. We found TFs associated with one branch independent of wounding conditions (msx2b versus prrx1a, hmx1, hoxb2a, and rx2), and TFs that showed condition-dependent patterns (Supplemental Fig. 9B). Most of the candidates do not have known neutrophil-associated functions but are associated with retinoic acid (RA) signaling, including suppression (alx4a, dlx4a, and sox9a/b) (35–37), activation (hmx1, sox4b, msx2b, hoxb2a, and nr6a1a) (38–42), and serving as coreceptors (nr4a2b) (43). Altogether, we identified the state transition paths in zebrafish neutrophils compared with human neutrophils and identified TFs that potentially drive the deviations.

Cells in a living organism communicate with one another to adjust their function to their cellular microenvironment. In search of cell–cell communication pathways of early innate immune responses that are activated or suppressed in burn wounding conditions, we performed CellChat analysis using all myeloid subsets (28). We converted a curated human ligand–receptor pair database for use in zebrafish to infer outgoing and incoming signals in each population in a context-dependent manner as described in Materials and Methods. At the global level, the total number of interactions between myeloid subsets in both burned and unwounded conditions first drops, then increases as time goes from T1 to T3 (Fig. 4A, 4B). However, the average strength of these interactions, measured by coexpression of ligand and receptor genes, presents a condition-specific pattern, where it mostly remains high in burned samples, but rapidly decreases in the unwounded condition. Zooming in on the subset level, we found a differential contribution from myeloid subsets to the overall level of outgoing and incoming signals in the burned versus unwounded comparison (Fig. 4C). At T1, Precursor_Neuts, Immature_Neuts, and Mature_Neuts followed the same enrichment pattern, with relatively higher signals in both outgoing and incoming directions in burn. At T2, all neutrophil subsets picked up this pattern. In addition, all subsets except for lgals3bpb_Mac2 showed relatively higher incoming signals in burn conditions. At T3, all myeloid subsets except for Mature_Neuts and lgals3bpb_Mac2 showed incoming enrichment in the burn condition, whereas all except for Cycling_Mac1 and lgals3bpb_Mac1 showed outgoing enrichment in the burn condition. Overall, myeloid cells in burn had extended signaling interactions until 48 hpb, with neutrophils showing more striking differences between conditions.

We next determined the key ligand–receptor pairs that contribute to the burn-induced interactions. We calculated information flow and ranked each ligand–receptor pair by the differences between burn and unwounded conditions and the statistical significance of the comparison (Fig. 4D). We looked for common and unique ligand–receptor pairs that are relatively enriched in burn conditions across time points (Fig. 4E). Then, we narrowed down our search for candidate pathways by focusing only on the ones with high levels of overall information flow. Among the candidate pathways, THY1 and IL16 signaling both showed a high burn-to-unwounded ratio only at early time points, with the former in the incoming direction and the latter in the outgoing direction for neutrophils (Fig. 4D, Supplemental Fig. 10A). Ligands of THY1 bind integrin and play a role in mediating endothelial adhesion and migration of myeloid cells (44–46). Binding of recombinant human THY1 induces MMP9 expression in neutrophils, which in turn enhances their migration ability (47). This aligns with the collagen metabolic function enrichment pattern in Immature_Neuts and il6r_padi2_Neuts, suggesting that their migration regulation is conserved with humans and mice. The IL6 pathway presents a trend-switching pattern, where it shows higher information flow in burn conditions in T1 and T3, but lower conditions in T2. By examining the networks formed by IL6–IL6R interactions among myeloid subsets, we found that the IL6 signal is mostly sent out by macrophages and received by neutrophils (Fig. 4F, Supplemental Fig. 10B). Previous work using a zebrafish il6r mutant has shown that the il6/il6r axis is essential for neutrophil but not macrophage recruitment to the burn wound area (11). This matches with the high-in-burn pattern for IL6–IL6R information flow at T1 and T3. T2 (24 hpb) is the time point when most of the neutrophils would have migrated out of the wound site in wild-type larval zebrafish after burn wounding. Therefore, we think that the reverse pattern in T2 might suggest the need for il6–il6r signal attenuation for inflammation resolution. This finding points to the hidden complexity of the signaling pathways at the myeloid subset level.

To test whether our findings in zebrafish burn models could translate to human burns, we examined the association between human orthologs of zebrafish myeloid subset signatures and burn severity in human patients using published microarray data (Fig. 5A) (29). We extracted top expressed genes from myeloid subsets in zebrafish and converted them to human orthologs using the DIOPT score (18). To gain a refined range of human burn severity measured by TBSA, we divided the TBSA levels into four categories: 0–10% (low), 11–20% (medium), 21–30% (high), and 31–40% (very high) with >100 samples per each category. First, we asked whether any myeloid subset signatures are associated with the TBSA levels and found that Immature_Neuts signatures positively associate with TBSA, whereas lgals3bpb_Mac1 negatively associates with TBSA (p < 2.22e−16 for both subset comparing between low and very high TBSA; Fig. 5B). Intriguingly, the trend is preserved with only samples collected at day 1 (with sample size reduced by half; both p < 1e−3 comparing between low and very high TBSA), indicating this might be a signature of early innate response in human burns. To ensure that the correlation was not observed due to the large sample size, we permuted TBSA levels for 10,000 times and found no correlation (p > 0.1), as expected. Furthermore, we aimed to identify whether there are any individual genes that correlate with the TBSA levels and found that GYG1, GPR84, CAPG, and MMP8 positively correlated with TBSA levels (Fig. 5C), whereas CD74, RUNX3, and PSAP negatively correlated with TBSA levels among the top 50 genes in each subset.

Among these burn severity–associated marker genes, GPR84, CD74, and PSAP all locate on the plasma membrane, making them directly detectable through flow cytometry. We used GPR84, a G_αi-coupled receptor, with its activation by medium-chain fatty acids inducing proinflammatory response in neutrophils (48), as the top candidate to test whether the RNA level in whole blood association with burn severity is reflected by protein level on the cell surface. We performed flow cytometry using peripheral blood collected from both burn patients and healthy donors with generally matching demographics (see Materials and Methods; Supplemental Table I). We found that there is a higher frequency of GPR84⁺ neutrophils in the blood neutrophils of burn patients compared with that of healthy controls (Fig. 5D, paired sample t test, p = 0.022), suggesting its involvement in the human burn wound healing process. Taken together, these findings demonstrate the potential of translating zebrafish burn wounding research to better understand innate immune responses in human burn patients.

Neutrophil heterogeneity at both the phenotypic and molecular level has been suggested by recent studies. However, to what degree such heterogeneity exists in tissue damage, in particular, in burn injury, has not been well characterized. In this work, we focused on the myeloid diversity in the burn wound healing response using larval zebrafish as a model. By performing transcriptional profiling at single-cell resolution, we identified four neutrophil subsets and six macrophage subsets in both normal development from 3 to 5 dpf and burn wound healing from 6 to 48 hpb. These subsets are present at the earliest time point we captured in development, suggesting that fish have developed this diverse group of myeloid cells as early as 3 dpf. In addition, the general population structure remains the same in burn wound healing compared with the unwounded control. This consistent distribution of myeloid subsets suggests that burn wounds do not induce the generation of new myeloid populations.

Owing to the limited molecular definition of neutrophils in zebrafish, we took advantage of published human and mouse transcriptomes for subset identity confirmation. Through label transfer between human and fish neutrophils, we found that three out of the four neutrophil subsets in fish could map to human neutrophil stages. Among these, Mature_Neuts showed more identity overlap with human neutrophils not in the mature stage, suggesting that the zebrafish neutrophils could be less mature than the human counterparts. Trajectory analysis further confirmed that these three subsets follow a maturation path similar to that in humans. This subset-level conservation makes fish-based wounding models for studying innate immunology more relevant to humans.

Il6r_padi2_Neuts, the neutrophil subset that does not score well in human neutrophil stages, lies on a branch separated from the trajectory of maturation. RNA velocity analysis also suggests that this population might not be an end state, as we found with the Mature_Neuts. To further understand the factors driving the separation between il6r_padi2_Neuts and Mature_Neuts, we constructed gene regulatory networks for each and identified key TFs associated respectively. Most of the candidates or their homologs are known to be involved in RA signaling. RA is essential for neutrophil maturation, as it has been studied in the context of acute promyelocytic leukemia patients, cell lines, and knockout mouse models (49). The subset-dependent RA signatures suggest that RA signaling might be involved in finer control of neutrophil maturation in zebrafish. Further work is needed to explore the il6r_padi2_Neuts subset and its development in zebrafish.

To understand the roles of these neutrophil subsets in burn wound healing, we performed functional annotation and compared between burned and unwounded conditions. Although general expression of GO term–associated genes present a similar subset enrichment pattern in both conditions, a time point–dependent burn response could be observed in certain terms. For instance, genes associated with collagen catabolic processes and cyclic nucleotide degradation are expressed at much higher levels in the burned condition at 6 and 48 hpb, respectively. This suggests that burn could induce time-sensitive functional enhancement in neutrophil subsets.

Burn wounding not only affects functions of individual neutrophil subsets, but it also changes the way they communicate. Cell–cell interaction analysis revealed more outgoing and incoming signaling from neutrophils in burned than in unwounded conditions. In particular, THY1, IL16, and IL6 stood out as candidate ligand–receptor pairs that possess higher information flow in the burned condition. Their enrichment patterns aligned well with functional study from humans and mice, suggesting their roles in myeloid recruitment and migration. In our work, we relied on the coexpression patterns in myeloid subsets to infer the communication between them. However, other cell types could express the same ligands or receptors and be involved in the communication as well. Although other cell types were captured in the cell sorting, we chose to leave them out due to the potential bias in our sorting strategy to more likely be obtaining cells with autofluorescence from populations other than the labeled myeloid cells. To better characterize cell–cell interaction between neutrophil subsets and other cell types in burn wound healing, an all-inclusive profiling that samples other cell types in appropriate ratios is needed in future work.

As an end goal of this study, we evaluated the translational value of the larval zebrafish model for human burn wounding. We confirmed associations between zebrafish neutrophil subset signatures and human burn severity using a human microarray dataset. Although other factors such as presence of inhalation injury, comorbidities of the patient, and depth of burn also contribute, TBSA is still the major factor in determining severity of burn. These findings suggest the potential of using zebrafish as a model to study human burn responses at the molecular level. In fact, neutrophils carrying these signatures could be found in a much higher level in the peripheral blood of burn patients than healthy donors. With more in-depth characterization aligning burn wound healing time courses between humans and zebrafish, the zebrafish model will become more helpful in facilitating our understanding of cell type–specific responses in burn in humans.

The authors have no financial conflicts of interest.

We thank all members of the Huttenlocher lab, the Dinh lab, and the Gibson lab for discussion and support in completing this work. We thank the Gene Expression Center at the University of Wisconsin–Madison for single-cell library construction and sequencing assistance, and the University of Wisconsin Carbone Cancer Center Flow Cytometry Laboratory for use of its facilities and services. We thank the University of Wisconsin Carbone Cancer Center BioBank for use of its facilities and services. We thank Dr. Montaldo Elisa and Dr. Ostuni Renato for providing us with processed data for comparing human and zebrafish neutrophils.