Blastemal progenitors modulate immune signaling during early limb regeneration

Salamanders, including newts and axolotls, possess the ability to fully regenerate their limbs throughout their lifespan. This regenerative capacity requires the formation of a transient cellular structure distal to the amputation plane, known as the blastema, which comprises progenitors derived from multiple different tissues (Kragl et al., 2009; Sandoval-Guzman et al., 2014; Currie et al., 2016). Although salamanders can perform this process regularly, mammals are incapable of forming a blastema that can give rise to an entire limb. Understanding the mechanisms underlying the initiation of blastema formation may provide important insights into unlocking human regenerative potential.

Blastema formation requires coordinated proliferation and migration of progenitors derived from muscle, bone, dermal fibroblasts, connective tissue and other tissues (Kragl et al., 2009). Researchers have shown that, although the blastema itself is transcriptionally similar to the limb bud during development, regenerating limbs exhibit unique transcriptional profiles during early stages of regeneration before the blastema has formed (Knapp et al., 2013). Elucidating these signals, as well as their specific activities, is crucial for understanding why and how salamanders respond to amputations with blastema formation. Several important signaling molecules and pathways have already been identified during early stages of appendage regeneration. For example, axolotl MARCKS-like protein (axMLP), a molecule that is secreted from the wound epidermis, is required for the induction of tail regeneration and is capable of inducing cell proliferation in intact limbs (Sugiura et al., 2016). Moreover, several core developmental signaling pathways, including FGF and Wnt, are required for appendage regeneration in axolotls and other species (McCusker et al., 2015; Haas and Whited, 2017; Stocum, 2017).

Several studies have investigated bulk transcriptional changes during different stages of axolotl limb regeneration and successfully identified genes that may play important roles in blastema formation and maintenance (Monaghan et al., 2009; Campbell et al., 2011; Knapp et al., 2013; Stewart et al., 2013; Wu et al., 2013; Voss et al., 2015; Bryant et al., 2017; Gerber et al., 2018). More recently, single cell analysis specifically in the connective tissue lineage has elucidated the molecular transitional states during dedifferentiation of mature connective tissue to a progenitor state during limb regeneration (Gerber et al., 2018).

In the present study, we sought to investigate the distinct transcriptional programs active within blastemal progenitors irrespective of lineage, as well as the surrounding tissues, to better understand the genetic programs and signaling interactions that govern the initiation of blastema formation. We reasoned that blastemal progenitors would be among early proliferating cells in stump tissues following amputation. Therefore, to enrich for blastemal progenitors, we transcriptionally profiled dividing cells from regenerating stump tissues during early stages of axolotl limb regeneration. This approach allowed us to examine the transcriptional signature of early dividing cells, inclusive of blastemal progenitors, and to identify potential novel modulators of early blastemal cell induction and maintenance. We additionally profiled non-dividing cells in regenerating stump tissues and the wound epidermis at similar stages for comparison. Pathway analysis of the discrete transcriptional programs of all three subpopulations revealed differential patterns of signaling pathway activation/inhibition and unexpectedly uncovered the strong suppression of several core developmental signaling pathways throughout regenerating stump tissues. Finally, we demonstrate that one dividing cell-enriched candidate, interleukin-8, is a blastemal progenitor-derived regulator of myeloid cell dynamics during the transition from wound healing into blastema formation, and that signaling through its cognate receptor, CXCR-1/2, is necessary for limb regeneration.

Cell cycle re-entry is an integral event for the initiation of limb regeneration prior to blastema formation. We therefore reasoned that we would be able to exploit differences in the cell cycle, such as DNA content, as a means of enriching for blastemal progenitors during early stages of regeneration. We transcriptionally profiled three total cellular fractions at 0, 4 and 5 days post-amputation (dpa): stump-derived dividing (4N) and non-dividing (2N) cells, as well as the whole wound epidermis (Fig. 1A). We developed and optimized a protocol to perform DAPI staining in conjunction with FACS to separate 4N and 2N cells in the regenerating stump during early regeneration (Fig. 1A, Fig. S1A-C). Intact limb tissue is referred to as the 0 dpa timepoint, which was collected to serve as a non-regenerating control.

Principal component analysis (PCA) of the transcriptional profiles of all samples revealed four clusters representing tissue type (stump-derived or epidermis) and regeneration status (Fig. S2A). Differential expression analysis of transcripts between the 4N (proliferating) and 2N (non-dividing) fractions at all three timepoints showed the expected enrichment of cell cycle gene expression in the proliferating fraction (Fig. S2B). Moreover, known blastemal cell markers, including prrx-1 (Satoh et al., 2011; Gerber et al., 2018) and kazd1 (Bryant et al., 2017), were enriched in the stump-dividing cells. Genes known to be upregulated early upon amputation but more likely to be modulating the extracellular matrix degradation rather than marking progenitor cells, such as mmp9 and mmp3 (Vinarsky et al., 2005; Monaghan et al., 2012; Stewart et al., 2013), were enriched in non-dividing stump tissues (Fig. S2C). To investigate the specificity of our strategy to enrich for blastemal progenitors rather than other dividing cell types, such as immune cell infiltrates, we examined the predicted activation and inhibition of immune signaling pathways by applying Ingenuity Pathway Analysis (IPA) software to the differentially expressed transcripts in each cell population relative to the non-regenerating limb. Most innate immune signaling pathways were inhibited within dividing cells (Fig. S3, Table S1), while activated in the other fractions, suggesting that we did not preferentially enrich for dividing immune cells. In all, these observations validated the enrichment of dividing blastemal progenitors using this protocol and allowed us to examine unique gene expression patterns of these cells.

Overall, a total of 21,077 transcripts (19,417 genes at 4 dpa) and 16,510 transcripts (15,205 genes at 5 dpa) were differentially expressed across all three fractions during early regeneration (Fig. 1B,C). Of these, only a strikingly small percentage were commonly up- or downregulated across all tissues relative to non-regenerating tissue at 0 dpa (1.51% at 4 dpa and 0.58% at 5 dpa), strongly suggesting that all three subpopulations initiate distinct transcriptional programs following amputation. Lists of differentially expressed transcripts for all analyses have been deposited in GEO under accession number GSE111213. The most highly expressed and commonly upregulated transcripts corresponded to enzymes involved in modulating extracellular matrix (ECM) degradation (mmp18, mmp2, timp1, tena and adam8) and to transcription factors involved in limb development and regeneration, including sall4 (Neff et al., 2005; Akiyama et al., 2015; Erickson et al., 2016) and runx-1 (Umansky et al., 2015; Deltcheva and Nimmo, 2017), indicating that they may play a role in orchestrating regeneration across all tissues.

We examined signaling interactions between regenerating subpopulations through IPA analysis of differentially expressed transcripts in each fraction at both timepoints. IPA software uses algorithms that account for the expression levels of signaling components, activators and inhibitors of well-known signaling pathways to predict whether a particular pathway is activated or inhibited. Surprisingly, pathway analysis revealed that several growth factor (e.g. FGF, Notch) (Fig. 1D) and intracellular (e.g. calcium and mTOR) signaling pathways implicated in regeneration were strongly inhibited in regenerating stump tissues during these early timepoints (Table S2). Notable exceptions include a handful of pathways: Hippo, canonical Wnt and TGF-β signaling. Among the repressed pathways were neurovascular signaling pathways (CNTF, NGF, neuregulin, VEGF signaling) that were activated in the wound epidermis and suppressed in stump tissues. These data further suggest that the wound epidermis orchestrates early neurovascular regeneration and that early repression of these pathways in regenerating stump tissues may be required. In all, these results show that exact spatiotemporal modulation of growth factor signaling pathways in specific tissues occurs during early stages of regeneration.

We next used our dataset to identify transcripts enriched in dividing cells at both regenerating timepoints (see Materials and Methods for filtering criteria). A total of 1217 transcripts (1181 genes) were enriched in dividing cells and, of these, only 298 transcripts (265 genes) were annotated. A heatmap representing the top 75 annotated and most highly expressed enriched transcripts, with little to no expression in the wound epidermis, is shown in Fig. 2A. A complete list of all annotated enriched transcripts identified can be found in Table S3. The small fraction of annotated enriched transcripts suggests that many key modulators of early regenerative events might be novel genes or at least may not have identifiable orthologous genes in other species present in existing datasets.

Among the most highly expressed and enriched transcripts in dividing cells were regenerative ECM components, including tenascin (tena), collagens (co5a2, co5a1, co1a2, coba1 and coca1), emilin1 (emil1) and fibrillin 2 (fbn2), suggesting that dividing cells play an early role in building the blastemal niche. Most notably, the transcriptional signatures of early dividing cells suggest that they may be directly regulated by TGF-β signaling. Pathway analysis revealed that TGF-β signaling was specifically activated in early dividing cells (Fig. 1D). Further examination of TGF-β signaling components revealed upregulation of both up- and downstream regulators (tgf-β1, tgf-β1r, smad2, ltbp1 and inhba) as well as direct targets of TGF-β signaling associated with epithelial-to-mesenchymal transition (EMT), such as snail1 and twist1 (Fig. 2B). These data indicate that an autoregulatory TGF-β signaling network is established both intra- and extracellularly in dividing cells, including blastemal progenitors, during early stages of regeneration.

To validate our differential gene expression findings and to learn more about where these transcripts are expressed, we performed time course RNA in situ hybridization on two candidates with high enrichment in dividing cells: transmembrane protein 119 (tm119) and E3 ubiquitin-protein ligase (lin41) (Fig. 2C-H). In situ hybridization confirmed that these transcripts were specific to early dividing cells and likely expressed in blastemal progenitors during early stages of regeneration. Tm119 has been shown to play a role in bone development and osteoblast proliferation (Kanamoto et al., 2009; Mizuhashi et al., 2012, 2015), whereas lin41 has a highly conserved role in stem cell maintenance as well as cellular reprogramming (Slack et al., 2000; Worringer et al., 2014). Both candidates showed little to no expression in uninjured limbs (Fig. 2C,F). Co-expression of tm119 and lin41 with a dividing cell marker, top2a, via double in situ hybridization at 7 dpa validated that these transcripts were expressed within dividing cells (Fig. 2D,G). Moreover, in situ hybridization of tm119 and lin41 at 21 dpa revealed strong expression of tm119 across the blastema, and also lower, but pan-blastemal, expression of lin41 (Fig. 2E,H), suggesting that they were indeed expressed within blastemal progenitors at earlier stages.

Interestingly, we noticed enriched expression of several inflammatory cytokines within early mitotically active cells, suggesting that blastemal progenitors may play a role in immune regulation. Among these was il-8, a secreted cytokine that has well known roles in inflammation, angiogenesis and proliferation (Russo et al., 2014), but no current known role in limb regeneration. We therefore chose to focus on il-8 for further study.

Time course in situ hybridization for il-8 during limb regeneration revealed a strong, transient upregulation during early stages of regeneration (Fig. 3). il-8 was undetectable in non-regenerating limbs (Fig. 3A,A′) and highly induced upon amputation (Fig. 3B-E′). At 1 dpa, il-8 was expressed in cells lining the bone (Fig. 3B,B′) and in the dermis. il-8 expression peaked at 3 dpa in mesenchymal cells within regenerating stump tissues and began to decrease by 7 dpa (Fig. 3C-D′). By early blastemal stages at 14 dpa, il-8 was weakly and sparsely expressed throughout the blastema and basal layers of the apical epithelial cap (AEC) (Fig. 3E,E′).

Previous studies have shown that il-8 is highly expressed in invading monocytes in wound healing and tumorigenic contexts (Fu et al., 2015; Williams et al., 2016). In contrast, double in situ hybridization of il-8 with two blastemal cell markers prrx-1 and kazd1, as well as csf1-r, a monocyte marker in mammals and amphibians (Grayfer et al., 2014; Stanley and Chitu, 2014), confirmed that il-8 was indeed expressed within blastemal progenitors and not monocytes. At both 3 and 7 dpa, il-8 was co-expressed in a subpopulation of prrx-1⁺ cells (Fig. 3F), and exhibited little to no expression in csf1-r⁺ monocytes (Fig. 3G). Quantification of prrx-1⁺ and il-8⁺ populations revealed high concordance between prrx-1 and il-8 expression. At 3 dpa, 99.0% of prrx-1⁺ cells expressed il-8 and 96.6% of il-8⁺ cells expressed prrx-1. At 7 dpa, only 67.9% of prrx-1⁺ cells expressed il-8, but 97.5% of il-8⁺ cells expressed prrx-1, consistent with a decrease in il-8 expression over time (Fig. 3H). In contrast, at 3 dpa, only 3.8% of il-8⁺ cells were csf1r⁺ and 1.7% of csf1r⁺ cells were il-8⁺, whereas at 7 dpa, 4.4% of il-8⁺ cells were csf1r⁺ and 1.5% of csf1r⁺ cells were il-8⁺ (Fig. 3I). il-8 was also expressed in a subset of kazd1⁺ blastemal progenitors at 7 dpa (Fig. 3J). Altogether, these data provide evidence that il-8 is strongly expressed within a subset of blastemal progenitors primarily during early stages of limb regeneration.

Because il-8 was strongly expressed in blastemal progenitors, we examined whether il-8 expression was sufficient to induce cell behaviors characteristic of the initiation of blastema formation, such as immune cell recruitment or cellular proliferation. We designed a myc-tagged il-8 overexpression vector (pCMV-IL8myc-T2A-tdTomato) and control vector (pCMV-T2A-tdTomato), and validated secretion of myc-tagged IL-8 protein in 293T cells (Fig. S4A-E). We injected and electroporated either the control tdTomato vector or the il-8 overexpression vector into intact limbs of non-regenerating animals, briefly pulsed the animals with EdU at 3 days post-injection, and collected samples 24 h later (Fig. 4A).

As IL-8 is a well-known inflammatory cytokine (Zeilhofer and Schorr, 2000), we first examined whether ectopic expression of il-8 induced monocyte and granulocyte recruitment by performing α-Naphthol Acetate (NSE) staining and Naphthol AS-D Chloroacetate (NCAE) staining, respectively. A 2.3-fold increase in NSE⁺ monocytes (117.53 vs 51.94 cells/mm², P=0.0041) and a 2.0-fold increase in NCAE⁺ granulocytes (121.54 versus 61.97 cells/mm², P=0.0353) was observed in limbs expressing il-8 relative to controls (Fig. 4B-D′), indicating that IL-8 is sufficient to recruit myeloid cells.

We also assessed whether IL-8 could induce cellular proliferation by comparing the percentage of total EdU⁺ cells between tdTomato and il-8-expressing limbs. We observed a modest, yet significant, increase in the total number of EdU⁺ cells in limbs overexpressing il-8 compared with control tdTomato limbs (13.92% versus 11.27%, P=0.0015) (Fig. 5A,C). Quantification of dividing bone/perichondrial, epidermal, endothelial and satellite cells revealed that this total increase in dividing cells was primarily due to an increase in the percentage of dividing bone/perichondrial cells (17.78% versus 8.71%, P=0.004), although there is a modest increase in the percentage of dividing epidermal cells (24.60% versus 19.70%, P=0.0146) (Fig. 5B-E). No increase in proliferation was observed in CD34⁺ endothelial cells or pax7⁺ muscle satellite cells. These data suggest that ectopic il-8 expression is sufficient to promote proliferation of bone/perichondrial and epidermal cells in a non-regenerative context, i.e. without amputation.

We next examined whether loss of function of il-8 may impair limb regeneration. To this end, we performed a dual pulse morpholino experiment (Fig. 6A) using translation-blocking il-8 targeted or control morpholinos, injected into either the right or left forelimb, at 2 days pre-amputation and 3 dpa. This treatment produced a strong and transient knockdown of endogenous IL-8 protein levels (Fig. S4F). il-8 morpholino-treated limbs displayed delayed blastemal outgrowth (Fig. 6B). Quantification at 21 dpa revealed a 37% decrease in blastemal length (572.89 vs 915.40 µm, P<0.0001) and a 32% reduction in blastemal area (491.29 vs 723.90 µm, P=0.0005) (Fig. 6C-F). In addition, control limbs reached digit stages of differentiation before il-8 morphant limbs (Fig. 6B), indicating that transient knockdown of IL-8 delays limb regeneration. At this time point, there was no defect in cell proliferation in either blastemal cells or the wound epidermis (Fig. S5) and no visible phenotypes with the blastema aside from the difference in size. These observations in conjunction with the strong expression of il-8 during early stages of limb regeneration suggested that il-8 morphant limbs were likely exhibiting defects at earlier timepoints in regeneration.

As ectopic il-8 expression was sufficient to recruit immune cells, we hypothesized that il-8 morphant limbs may exhibit deficiencies in myeloid cell recruitment. Surprisingly, we observed no significant deficiency in recruitment of either NCAE⁺ granulocytes or NSE⁺ monocytes to the amputation site at either 1 dpa or 5 dpa (Fig. 7A,D). However, the mean density of granulocytes was lower overall than that of controls at both of these timepoints, suggesting that il-8 may play a minor role in recruiting granulocytes during wound healing. More interestingly, there was a significant decrease in both granulocytes (374.16 versus 483.24 cells/mm², P<0.01) and monocytes (198.20 versus 290.92 cells/mm², P<0.05) during later stages of wound healing at 7 dpa (Fig. 7). We did not observe a difference in overall cell proliferation or cell death at this timepoint. Therefore, the data suggest that il-8 likely plays a role in retaining both granulocytes and monocytes during the transition from wound healing to blastema formation.

To determine whether IL-8 was directly regulating myeloid cell behaviors, we examined the expression pattern of its cognate receptors, cxcr-1 and cxcr-2, during early stages of regeneration. Although other organisms, including mammals and zebrafish, have both receptors, analysis of publicly available transcriptomes for the axolotl (Bryant et al., 2017; Nowoshilow et al., 2018) revealed only one receptor with homology to both cxcr-1 and cxcr-2 in different species. We therefore refer to the detectable cognate receptor as cxcr-1/2. Double in situ hybridization of cxcr-1/2 with csf1r or mpo, a neutrophil marker (Walters et al., 2010), at 3 dpa revealed co-expression of cxcr-1/2 with subpopulations of csf1r⁺ and mpo⁺ cells (Fig. 8), suggesting that IL-8 likely directly acts on subpopulations of monocytes and granulocytes, including neutrophils during wound healing. Altogether, these data provide one of the first examples of blastemal progenitors as immunomodulators, specifically as the source of canonical chemokines, during wound healing, and further suggests that retention of myeloid cells during the transition to blastema formation is crucial.

Finally, we asked whether inhibition of CXCR-1/2 signaling would also impair limb regeneration. We perturbed the pathway by treating regenerating axolotls with a small molecule inhibitor of CXCR-1/2: SB-225002 (White et al., 1998) (Fig. 9). Inhibitor treatment beginning from the time of amputation (0 dpa) completely inhibited limb regeneration, while DMSO-treated limbs regenerated normally (Fig. 9B-E). Owing to the well-documented role of CXCR-1/2 signaling in immediate wound-associated inflammatory responses (Ha et al., 2017), we hypothesized that failure to regenerate in this context was mainly due to prevention of immediate inflammatory responses stimulated by limb amputation. Therefore, we treated regenerating animals with SB-225002, beginning at later stages of wound healing (3 dpa) and found that this also prevented regeneration (Fig. 9F-I), suggesting that prolonged CXCR-1/2 signaling is essential. Yet, treatment beginning after blastema formation (15 dpa) led to normal limb regeneration (Fig. 9J,K). In all, these data suggest that CXCR-1/2 signaling is necessary during early stages of regeneration, but dispensable for later stages.

Since il-8 morphants displayed defective retention of myeloid cells during late wound healing, we examined whether inhibition of CXCR-1/2 signaling affected myeloid cell behavior similarly during early stages of regeneration. Unexpectedly, monocytes and granulocytes in SB-225002-treated limbs (from 0 dpa) displayed unhealthy morphologies at 7 dpa. The majority of NSE⁺ and NCAE⁺ cells at the distal amputation site had formed apoptotic bodies, suggestive of myeloid cell death (Fig. 10A). Quantification of healthy NCAE⁺ and NSE⁺ myeloid cells revealed statistically significant decreases in both populations in SB-225002-treated limbs at 7 dpa compared with DMSO controls (Fig. 10B,C, NCAE⁺ granulocytes: 128.9 versus 362.3 cells/mm², P<0.0001, NSE⁺ monocytes: 94.71 versus 240.6 cells/mm², P=0.0022). In concurrence with these observations, TUNEL staining of inhibitor-treated limbs revealed higher levels of cell death near the amputation plane at 7 dpa (6.03% versus 1.47%, P=0.0154) (Fig. 10D), suggesting that myeloid cells were likely undergoing apoptosis. Interestingly, increased cell death was not observed in il-8 morphant limbs, suggesting that other cytokines in addition to il-8 may synergistically modulate myeloid cell behaviors. Altogether, these data newly identify CXCR-1/2 signaling as an immunomodulatory pathway that is, in part, regulated by IL-8, and is crucial for successful blastema formation.

Elucidating the molecular mechanisms underlying the initiation of regeneration is key to understanding the difference in responses to tissue loss between regenerative and non-regenerative organisms. Here, we transcriptionally profiled distinct regenerating subpopulations during early pre-blastemal stages of limb regeneration to differentiate gene expression changes that occur in early dividing cells, inclusive of blastemal progenitors, from those of surrounding tissues. Using this dataset, we have gained insights into the patterns of suppression and activation of signaling pathways present within subsets of early regenerating limb tissues. Notably, our examination of the expression profiles of early dividing cells revealed an immunomodulatory role for blastemal progenitors during early stages of regeneration.

The transcriptional signatures of early dividing cells and blastemal progenitors suggest that the formation of the early blastemal niche may be regulated by canonical Wnt, Hippo and TGF-β signaling. Dividing cells showed enriched expression of regenerative ECM components, many of which are important for regeneration (Calve et al., 2010; Godwin et al., 2014), suggesting that they are drivers of the transition to a regenerative ECM. In addition, Hippo, Wnt and TGF-β signaling pathways were activated within dividing cells. As downstream effectors of all three pathways interact in development and tumorigenesis (McNeill and Woodgett, 2010; Attisano and Wrana, 2013), it is likely that synergy between these pathways is essential for early blastemal cell establishment and maintenance. Interestingly, TGF-β signaling, which is necessary for axolotl limb regeneration (Levesque et al., 2007; Denis et al., 2016), was specifically active in dividing cells. Our data further suggest that TGF-β signaling is sustained through autocrine feedback in early dividing cells, which exhibit exclusive upregulation of tgf-β1, tgf-β1r and smad-2, as well as regulators and direct targets, including ltbp1, twist1 and snail-1. snail-1 directs epithelial-to-mesenchymal transition (EMT) behaviors (Fuxe et al., 2010), and activates expression of both regenerative ECM components and twist family members, which are expressed in limb blastemal cells (Kragl et al., 2013; Bryant et al., 2017). Moreover, dividing cells highly expressed ECM components such as emilin-1 and fibrillin-2, which modulate TGF-β signaling through interactions with ltbp1 (Neptune et al., 2003; Randell and Daneshtalab, 2017), indicating that both an intra- and extracellular TGF-β signaling network is established.

Surprisingly, we observed the strong repression of many growth factor signaling pathways in early regenerating stump tissues. FGF, Notch, IGF-1, PDGF and non-canonical Wnt signaling pathways (PCP and Wnt/Ca⁺) all appear to be inhibited, yet many of these signaling pathways are necessary and/or sufficient for appendage regeneration in amphibians as well as zebrafish (Poss et al., 2000; Yokoyama et al., 2000, 2001; Stoick-Cooper et al., 2007; Chablais and Jazwinska, 2010; Satoh et al., 2011; Grotek et al., 2013; Makanae et al., 2014; Rodrigo Albors et al., 2015; Currie et al., 2016; Nacu et al., 2016; Shibata et al., 2016). Furthermore, signaling pathways involved in neuronal (e.g. neuregulin) and vascular (e.g. VEGF) regeneration (Yu et al., 2014; Farkas et al., 2016; Farkas and Monaghan, 2017; Ritenour and Dickie, 2017) were active in the wound epidermis, yet suppressed in regenerating stump tissues, signifying that early neurovascular regeneration is coordinated by the wound epidermis and that repression of these pathways within regenerating stump tissues may be important. A potential explanation for this apparent dichotomy lies in the timing of activation. Premature activation of signaling pathways such as Notch (Grotek et al., 2013) or non-canonical Wnt signaling (Stoick-Cooper et al., 2007) inhibits blastemal growth during zebrafish fin regeneration. Therefore, our findings suggest that repression of these pathways in early stages of regeneration may be necessary and that precise timely release of inhibition ensures successful regeneration. It is interesting to note that the analysis did not reveal strong predictions for activation or repression of other pathways essential for limb regeneration, such as sonic hedgehog (shh) signaling (Singh et al., 2015). These other pathways may be downstream and act at later stages of regeneration (or are regulated post-transcriptionally). Nevertheless, these findings could provide targetable insights for improving regenerative outcomes.

Most notably, we provide one of the first examples that blastemal progenitors play an early paracrine immunomodulatory role in appendage regeneration. Others have shown that IL-8 recruits neutrophils and macrophages in zebrafish organ and appendage regeneration (de Oliveira et al., 2013; Xu et al., 2018); however, the source of IL-8 was not examined. Contrary to other injury and tumorigenic contexts, where IL-8 is secreted from macrophages (Arango Duque and Descoteaux, 2014), we found that the primary source of  IL-8 in early stages of limb regeneration is a subpopulation of blastemal progenitors. Furthermore, the high concordance of expression between il-8 and prrx-1, recently demonstrated as a connective tissue blastemal cell marker (Gerber et al., 2018), suggests that il-8 is primarily expressed in early regenerating connective tissue. Meanwhile, its receptor cxcr1/2 is expressed in myeloid cells. Our functional data further shows the importance of signaling between blastemal progenitors and the immune system. il-8 knockdown delayed regeneration and CXCR-1/2 inhibition prevented limb regeneration. We also demonstrate that il-8 can recruit both myeloid cell-types in non-regenerating limbs. Thus, whereas leukocytes have been shown to regulate blastemal cells in other models (Nguyen-Chi et al., 2017), our findings newly suggest that blastemal progenitors direct immune support and that this bi-directional signaling is important for early stages of regeneration.

Failure to retain myeloid cells during the transition from wound healing to blastema formation may account for the regenerative delay in il-8 morphants. Monocyte-derived macrophages are necessary for successful appendage regeneration in the axolotl and other model systems, playing roles that include directing blastemal outgrowth, clearing senescent cells and modulating blastemal cell proliferation (Li et al., 2012; Godwin et al., 2013; Petrie et al., 2014; Yun et al., 2015; Nguyen-Chi et al., 2017; Simkin et al., 2017). Others have demonstrated the importance of pro-regenerative (M2) rather than pro-inflammatory (M1) macrophage subtypes in regeneration (Pei et al., 2016; Simkin et al., 2017). Therefore, it is possible that IL-8 may retain M2 macrophages during the initiation of blastema formation. The extent of the role of granulocytes in regeneration, however, seems to be more context dependent (Nakayama et al., 2011; Li et al., 2012; Kurimoto et al., 2013; Paris et al., 2016; Lindborg et al., 2017). Granulocytes, including neutrophils, clear debris during wound healing (Wang, 2018). Therefore, low levels of granulocytes in il-8 morphants may have slowed wound healing, delaying regeneration. Furthermore, in the zebrafish, IL-8 can act as a potent neutrophil chemoattractant (de Oliveira et al., 2013, 2016) or a chemokinetic molecule, aiding in neutrophil reverse migration from the sites of injury through Cxcr2 signaling (Powell et al., 2017). In contrast, our findings suggest that IL-8 may retain granulocytes, including neutrophils, during late wound healing rather than facilitating their exit strategy as in the zebrafish, highlighting potential species-specific differences in immune responses during regeneration. As IL-8 morphant limbs eventually regenerate, it is clear that other pathways compensate to ensure wound healing resolution.

Last, we show that CXCR-1/2 signaling is necessary during early, but not late, stages of limb regeneration. CXCR-1/2 inhibitor-treated limbs beginning during early and late stages of wound healing failed to form a blastema. The high level of myeloid cell death, which was not observed in il-8 morphants, during early regeneration was likely the main cause of failure to regenerate, suggesting that CXCR-1/2 signaling may serve as an important survival pathway for monocytes and granulocytes in early limb regeneration. Additionally, the difference in severity of phenotypes between il-8 morphant and CXCR-1/2 inhibitor-treated limbs is likely attributed to the fact that CXCR-1/2 binds to other interleukins, including IL-1 and IL-6 (Baggiolini et al., 1997). Therefore, it is possible that other cytokines are acting in concert with IL-8 to control CXCR-1/2-specific myeloid cell behaviors during early limb regeneration. Altogether, these data highly suggest that CXCR-1/2 signaling may be key in bridging communication between early blastemal cells and the immune system during the initiation of regeneration.

In conclusion, our approach reveals differential patterns of activation/suppression of core developmental signaling pathways within separate cellular subsets during early regeneration and a role for early dividing cells, inclusive of blastemal progenitors, in shaping the regenerative niche. We demonstrate that blastemal progenitors play an early immunomodulatory role through IL-8 and that signaling through its cognate receptor, CXCR-1/2, is necessary for limb regeneration. These findings highlight the importance of further characterizing the complexity of bi-directional crosstalk between the immune system and blastemal cells on a broader scale.

Axolotl (Ambystoma mexicanum) husbandry and surgeries were performed in accordance with the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC) and Institutional Animal Care and Use Committee (IACUC) guidelines at Harvard University. Sub-adult white and albino axolotls (15-18 cm) provided by the Ambystoma Genetic Stock Center (AGSC, University of Kentucky) were used for the initial RNA-sequencing experiment. Juvenile white axolotls (5-8 cm) were used for all morpholino and overexpression experiments. For CXCR-1/2 inhibitor experiments, animals (3-5 cm) were immersed in either 500 nM SB-225002 (Tocris Bioscience) or DMSO (Sigma) beginning at either 0, 3 or 15 dpa and solutions were changed daily.

Briefly, animals were anesthetized in 0.1% Tricaine (Sigma-Aldrich) and all four limbs were amputated at the mid-radius/ulna level. Bone was trimmed back to facilitate wound closure and regeneration. Approximately 2-3 mm of tissue directly proximal to the amputation plane was collected at either 0, 4 or 5 dpa. We chose to transcriptionally profile regenerating limbs at 4 and 5 dpa because preliminary experiments showed a distinct increase in cellular proliferation in the regenerating stump beginning at 3 dpa (data not shown here) and we wanted to capture the transcriptional signatures of dividing cells during the duration of cell-cycle re-entry. In order to obtain enough material for the sorting and sequencing protocol, we chose to perform the experiments at 4 and 5 dpa. Twelve limbs were pooled for each biological replicate per timepoint and the experiment was performed in biological triplicate.

To prep the tissue for FACS, the intact epidermis or wound epithelium was micro-dissected off and placed into 0.25% Trypsin-EDTA for 15 min with agitation at room temperature to dissociate epithelial cells from the dermal layer. The remaining stump tissues were micro-dissected further into small pieces with dissecting scissors and chemically dissociated in a solution composed of 5 mg/ml collagenase (Sigma-Aldrich), 7.3 mg/ml dispase II (Roche Diagnostics) and 1.36 mg/ml D-glucose in 80% PBS (Kumar and Brockes, 2007) for 15 min with agitation at room temperature. In order to obtain a representation of dermal cells in the stump tissue fraction, we also chemically dissociated micro-dissected intact epidermis or wound epithelia along with the stump tissues. For regenerating samples, the wound epidermis region, which is visibly transparent, was carefully removed to isolate adjacent full-thickness skin with both epidermal and dermal layers, and then dissociated with the stump tissue fraction. Dissociated cells were then transferred to a new tube (non-dissociated chunks of tissue were left behind) and the dissociation was serum inactivated. The following protocol for staining, FACS and RNA isolation was optimized and adapted from Hrvatin et al. (2014). The cell suspension was pelleted, re-suspended, passed through a 70 µM cell strainer, washed twice with 80% PBS and fixed in 4% paraformaldehyde/0.1% saponin with a 1:50 dilution of RNasin plus RNase Inhibitor (Promega) for 30 min at 4°C. Fixed cells were then washed twice with a 1% BSA/ 0.1% saponin (1:40 dilution RNasin plus RNase Inhibitor, Promega) in PBS wash buffer and DAPI staining (10 µg/ml) of cells was performed for 30 min at 4°C in a 0.1% saponin solution (1:20 RNasin plus RNase Inhibitor, Promega). The stump fraction of DAPI-stained cells were immediately sorted into 2N and 4N fractions using FACS, whereas the corresponding epithelial fraction was subjected to the same staining treatment, but not sorted. Cells were sorted into RNAlater solution (Invitrogen) and RNA was isolated with the RecoverAll Total Nucleic Acid Isolation Kit for FFPE (Invitrogen). The quality of the RNA samples was assessed using the Agilent RNA 6000 Pico kit on the Bioanalyzer 2100 (Agilent Technologies).

cDNA was synthesized from RNA samples using the NuGEN Ovation RNA-seq System V2 protocol (Integrated Sciences) according to the manufacturer's instructions using 10 ng of RNA as starting material. The quality and concentration of cDNA preps was then assessed using the Agilent DNA 1000 kit on the Bioanalyzer 2100 (Agilent Technologies). RT-PCR for cell cycle markers ccnb3, ccna2 and cdk1b was then performed on the cDNA generated from 4N and 2N cells to check for proper enrichment of dividing cells. Primer sequences are as follows: ccnb3-For, 5′-CACAAGAATCCAGTGCCACA-3′; ccnb3-Rev, 5′-CCTCCTTTGCAACAGTGTCC-3′; ccna2-For, 5′-GAACGTACAGCCTGGCAAG-3′; ccna2-Rev, 5′-CTGACGGCTGCTCCTTTG-3′; cdk1b-For, 5′-GCCAAACAACGAAATCTGGC-3′; cdk1b-Rev, 5′-AGGGTGGTTCAATGCCTCTT-3′.

To prepare cDNA sequencing libraries, 200 ng of cDNA was first sheared to a peak size of 200 bp using the Covaris S220 according to the manufacturer's protocol. Good quality and correct size distribution of sheared DNA fragments was assessed with the Agilent DNA High Sensitivity kit on the Bioanalyzer 2100 (Agilent Technologies). Sequencing libraries were then generated using the Wafergen PrepX Complete ILMN DNA Library kit (Takara Bio) protocol on the Apollo 324 NGS Library Prep System (Takara Bio). The quality of the DNA sequencing libraries was performed with the Agilent DNA High Sensitivity kit on the Bioanalyzer 2100 (Agilent Technologies). Concentration of the DNA libraries was doubly confirmed using the Qubit dsDNA HS Assay kit (ThermoFisher Scientific) and Kapa Illumina Library Quantification kit (Kapa Biosystems). Libraries were multiplexed and sequenced on either the Illumina Hiseq 2500 system (125 bp reads) or Nextseq 500 (150 bp reads) at the Harvard Bauer Core Sequencing Facility.

Reads were trimmed using Trimmomatic (Bolger et al., 2014) to a minimum length of 100 bp and poor quality reads were removed from the sequencing analysis. Alignment was performed on the trimmed reads using Kallisto (Bray et al., 2016) and the previously published well-annotated axolotl transcriptome (Bryant et al., 2017). Raw read data and the processed data matrix containing TMM-normalized TPM values for each sample have been deposited in GEO under accession number GSE111213. Differential expression analysis of genes and transcripts relative to the non-regenerating timepoint (0 dpa) for each fraction was performed using DESeq2 (Love et al., 2014) with an adjusted P-value cutoff of 0.05. Core and comparison analysis of differentially expressed transcript lists for each cellular population was performed using Ingenuity Pathway Analysis software (Qiagen). Uniprot IDs of transcript blast hits (extracted using the associated Trinotate file by Bryant et al., 2017) were converted to human IDs for IPA analysis. We focused on strongly activated or inhibited growth factor signaling pathways, i.e. absolute value of the Z-score>1 and signal detected across at least three differentially expressed analyses out of a total of six.

To identify transcripts enriched within dividing cells at both timepoints, we first compared dividing and non-dividing cells averaged at both regenerating timepoints (4 and 5 dpa) and identified 6834 differentially expressed transcripts (3510 of which were upregulated within dividing cells). Transcripts that were normally differentially expressed in non-regenerating dividing and non-dividing cells (at 0 dpa), such as cell cycle-associated transcripts, were filtered out. Of these, 628 total transcripts (583 genes) were regeneration specific, annotated and had at least a twofold change between dividing and non-dividing cells; only 298 transcripts in this list were upregulated in dividing cells.

Tissue was collected at 0, 7 and 21 dpa and fixed in 4% paraformaldehyde overnight at 4°C, washed in PBS, brought up a sucrose gradient to 30% sucrose, and embedded in OCT. The blocks were serially sectioned and 16 µm sections were collected. Custom RNAscope probes for the axolotl orthologs of top2a, cxcr-1, tm119, lin41, il-8, kazd1, mpo and prrx-1 were generated (Advanced Cell Diagnostics) in either the C1 or C2 channels. Double chromogenic section in situ hybridization was performed on frozen cryosections using the RNAscope 2.5 HD Duplex Detection Kit (Advanced Cell Diagnostics) according to the manufacturer's instructions.

The total numbers of prrx-1⁺only, il-8⁺ only, csf1r⁺ only or double-positive cells were quantified by counting cells that exhibited at least 10 puncta for either probe. These are categorized as strongly-expressing cells, according to RNAscope standards (score of 4+). These numbers were used to quantify the percentages of co-positive and single-positive cells.

To generate an il-8 overexpression construct, the open reading frame for il-8 was amplified out of cDNA from regenerating limbs at 7 dpa with primers that added on a NheI site, and a kozak sequence directly upstream of the ATG start codon and a HindIII site directly downstream of the coding sequence. This PCR product was cloned into a pCMV-TVA-T2A-tdTomato backbone (a gift from J. Whited, Harvard University, MA, USA) in place of TVA. An empty pCMV-T2A-tdTomato vector was used as a control.

To validate the expression of the vector, 293T cells were transfected with 50 µg of either the il-8 overexpression or tdTomato control construct using Lipofectamine 2000. The media and cells were separately collected at 48 h post-transfection (hpt). Media was concentrated with a 15 ml Amicon filter (3 K) (MilliporeSigma). Protein was extracted from cells using Trizol (ThermoFisher Scientific). Western blotting was performed on both il-8 overexpression and control transfected cell lysates and media using a mouse monoclonal anti-GAPDH (MAB374, MilliporeSigma, 1:2000) and polyclonal rabbit anti-myc tag antibody (ab9106, Abcam, 1:2000). Cells were additionally immunostained at 72 hpt for the myc tag (same antibody as above) and tdTomato using a polyclonal goat anti-tdTomato antibody (LS-C348313, LSBio, 1:1000).

Intact limbs of axolotls were injected and electroporated with myc tagged il-8 overexpression construct (1 µg/µl in PBS). At 3 days post-injection, animals were injected intraperitoneally with EdU (Life Technologies) at a concentration of 8 mg/kg 4 hours prior to tissue collection. Tissue was prepped and embedded in OCT in a similar manner to previously described above. The blocks were serially sectioned at 16 µm. EdU staining was performed using the Click-it EdU Alexa Fluor 488 Imaging Kit (ThermoFisher Scientific). For immunostaining, mouse anti-chick Pax7 (DSHB, 1:200) and rabbit anti-human CD34 (ab81289, Abcam, 1:200) antibodies were used.

The imaging analyses were all conducted blinded. For the quantification of il-8 overexpression limbs, three or four sections were imaged per limb and quantified. Total numbers of EdU⁺, CD34⁺, pax7⁺, CD34⁺ EdU⁺ and pax7⁺ EdU⁺ cells were counted, and percentages were calculated using DAPI⁺ nuclei totals for EdU percentages, pax7⁺ nuclei totals for dividing satellite cells and CD34⁺ cell totals for dividing endothelium. Dividing bone/perichondrial or epidermal cells were counted based on EdU⁺ cells of the cell type out of DAPI⁺ total of the cell type (counted by morphology). A two-tailed paired t-test was used for statistical analysis.

Staining of monocytes and granulocytes was performed using the α-Naphthyl Acetate (Non-specific Esterase) (NSE) kit or Naphthol AS-D Chloroacetate (Specific Esterase) (NCAE) Kit (Sigma-Aldrich), respectively, according to the manufacturer's instructions. The only modification to the protocol was a 10 min fixation step in citrate-acetone-methanol fixative. For the overexpression experiment, the area of each cross-section was measured using ImageJ analysis software. For characterization of the il-8 morphant limbs, tissue was collected and prepared at 1, 5 and 7 dpa. The area of the section 500 μm from the amputation plane was measured and the quantification was performed blinded. Every effort was made to ensure quantification around the same area in all limbs using the humerus and ulna as landmarks. A two-tailed paired t-test was used for all statistical analyses and all of the quantification was performed blinded (averaging two or three sections/limb). For characterization of DMSO- or SB-225002-treated limbs, limbs from animals treated from time of amputation were collected at 7 dpa, stained and analyzed as described above.

Axolotls were anesthetized in 0.1% tricaine and intact forelimbs were injected and electroporated with either control (left limb) or il-8-targeted morpholino (right limb). Approximately 3-4 µl of morpholinos were injected at a final concentration of 5 µM in PBS. At 2 days post-injection, both forelimbs of axolotls were amputated at the mid-radius/ulna level and control or il-8 morpholinos were injected again at 3 days post-amputation into regenerating stump tissue. Both translation-blocking and five-point mutation control morpholino antisense oligonucleotides were designed and generated by GeneTools against the il-8 ORF with the following sequences: control, 5′-CCGATCTTGATGCTCACCTCCTG-3′; il-8, 5′-CCGATGTTCATGGTGACCTGCTG-3′. To validate morpholino knockdown, il-8-targeted and control morpholino injected limbs from the same animal were collected and protein was extracted using Trizol (ThermoFisher Scientific). Western blotting was performed on protein extracts using the anti-GAPDH antibody described above, as well as a cross-reacting polyclonal mouse anti-chick IL-8 antibody (MBS2018201, MyBioSource, 1:400).

For analysis at the blastemal stages, animals were injected intraperitoneally with EdU (Life Technologies) at a concentration of 8 mg/kg 4 hours prior to tissue collection. Tissue was collected at 21 dpa and prepared as described above. Picro-mallory staining was performed on sections for histological analysis to analyze blastema length and area. Cell proliferation was assessed using the Click-it EdU Alexa Fluor 488 Imaging Kit (ThermoFisher Scientific). All imaging analysis was carried out in ImageJ and conducted blinded. A two-tailed paired t-test was used for statistical analysis.

TUNEL staining was performed using the In Situ Cell Death Detection Kit (Roche) as described previously (Zhu et al., 2012). The percentage of TUNEL⁺ nuclei (out of DAPI total) was quantified from two sections per limb and averaged.

Alcian Blue staining of DMSO or SB-225002-treated 40 dpa limbs was performed as described previously (Whited et al., 2013).

We thank J. Whited for providing the overexpression vector (pCMV-TVA-T2A-tdTomato), and for providing scientific guidance and helpful comments on the manuscript. We acknowledge I. Adatto, L. Krug and the Harvard Office of Animal Resources (OAR) for their dedicated animal care in the axolotl facilities, as well as S. Ionescu and J. Lavecchio for performing all of the FACS in the SCRB department flow cytometry core. In addition, we thank the Ambystoma Genetic Stock Center (AGSC), which provided many of the animals for the experiments conducted here. Last, we greatly thank J. Davis, J. Rivera-Feliciano, E. Rosado-Olivieri, N. Sharon, C. Kayatekin, all of the members of the Melton lab and also R. Amamoto for helpful scientific discussion and review of the manuscript.