Research Article

Single-cell analysis uncovers convergence of cell identities during axolotl limb regeneration

See allHide authors and affiliations

Science  26 Oct 2018:
Vol. 362, Issue 6413, eaaq0681
DOI: 10.1126/science.aaq0681

How the axolotl makes a new limb

Unlike most vertebrate limbs, the axolotl limb regenerates the skeleton after amputation. Dermal and interstitial fibroblasts have been thought to provide sources for skeletal regeneration, but it has been unclear whether preexisting stem cells or dedifferentiation of fibroblasts formed the blastema. Gerber et al. developed transgenic reporter animals to compare periskeletal cell and fibroblast contributions to regeneration. Callus-forming periskeletal cells extended existing bone, but fibroblasts built new limb segments. Single-cell transcriptomics and Brainbow-based lineage tracing revealed the lack of a preexisting stem cell. Instead, the heterogeneous population of fibroblasts lost their adult features to form a multipotent skeletal progenitor expressing the embryonic limb program.

Science, this issue p. eaaq0681

Structured Abstract

INTRODUCTION

Axolotls (Ambystoma mexicanum) and other salamanders are the only tetrapods that can regenerate whole limbs. During this complex process, changes in gene expression regulate the outgrowth of a new appendage, but how injury induces limb cells to form regenerative progenitors that differentiate into diverse cell fates is poorly understood. Tracking and molecular profiling of individual cells during limb regeneration would resolve distinct differentiation pathways and provide clues for how cells convert from a mature resting state into regenerative cell lineages.

RATIONALE

Axolotl limbs are composed of many different cell types originating from neural, myogenic, epidermal, and connective tissue (CT) lineages. Upon limb amputation, cells from nearby the amputation plane accumulate in a distinctive tissue called the blastema, which serves as a progenitor cell source to build the new limb. Transgenic axolotl strains in which descendants of distinct adult cell types can be labeled, tracked, and isolated during the regenerative process provide an opportunity to understand how particular cell lineages progress during blastema formation and subsequent limb regrowth. Combining transgenic axolotl strains with single-cell RNA sequencing (scRNA-seq) enables the tracking of individual cell types, as well as the reconstruction of the molecular steps underlying the regeneration process for these particular cell lineages. CT cells, descendants of lateral plate mesoderm, are the most abundant lineage contributing to the blastema and encompass bone and cartilage, tendons, periskeleton, and dermal and interstitial fibroblasts. These cells detect the position of the amputation site, leading to the regeneration of appropriate limb parts and making the CT a key cell lineage for deciphering and understanding molecular programs of regeneration.

RESULTS

We used an inducible Cre-loxP fluorescence system to establish genetically marked transgenic axolotl strains for isolating CT cells from adult limb tissue as well as CT descendants in the blastema. We used scRNA-seq to molecularly profile CT cells along a dense time course of blastema formation and the outgrowth of the regenerated arm, as well as stages of embryonic limb development. This profiling indicated that CT cells express adult phenotypes that are lost upon the induction of regeneration. The heterogeneous population of CT-derived cells converges into a homogeneous and transient blastema progenitor state that at later stages recapitulates an embryonic limb bud–like program. Notably, we did not find evidence of CT stem cells or blastema-like precursors in the mature arm. We found that CT subtypes have spatially restricted contributions to proximal and distal compartments in the regenerated limb. Specifically, a particular CT subtype—periskeletal cells—extended the severed skeleton at the amputation site whereas fibroblastic CT cells de novo regenerated distal skeletal segments. By using high-throughput single-cell transcriptomics and Brainbow axolotl-based clonal lineage tracing, we could follow the redifferentiation trajectories of CT lineages during the final stages of regeneration. These findings established the formation of a multipotent skeletal progenitor cell that contributed to tendons, ligaments, skeleton, periskeleton, and fibroblasts.

CONCLUSION

CT cells are a key cell type for understanding regeneration because they form the patterned limb skeleton that guides the regeneration of the other limb tissues, such as muscle. Because of the cells’ heterogeneity and intermingling with other cell types, it had been difficult to study how CT forms regenerative blastema cells. The use of these newly generated transgenic reporter strains combined with single-cell transcriptomic analysis and clonal tracing have allowed us to determine that CT cells with diverse molecular features traverse through a distinctive molecular state as they dedifferentiate into a common, multipotent progenitor resembling an embryonic limb bud cell. In the future, it will be important to test which components of the transition state are required for the dedifferentiation process. Furthermore, this work opens the possibility to examine how regeneration-associated genes and their associated chromatin structure are regulated during this transition. Lastly, the work raises the possibility that the limited regeneration seen among mammals is due to an inability to reprogram CT to such embryonic states.

Regeneration of the axolotl upper arm CT.

Transgenic axolotl strains were used to isolate CT cells at different stages during limb regeneration. Single-cell transcriptomes were generated, and transcriptional signatures were used to track the fate of cells during regeneration. We found that heterogeneous CT cell types in the adult upper arm funnel into a homogeneous multipotent skeletal progenitor state that recapitulates axolotl limb development, before their redifferentiation into the heterogeneous CT cell subtypes.

Abstract

Amputation of the axolotl forelimb results in the formation of a blastema, a transient tissue where progenitor cells accumulate prior to limb regeneration. However, the molecular understanding of blastema formation had previously been hampered by the inability to identify and isolate blastema precursor cells in the adult tissue. We have used a combination of Cre-loxP reporter lineage tracking and single-cell messenger RNA sequencing (scRNA-seq) to molecularly track mature connective tissue (CT) cell heterogeneity and its transition to a limb blastema state. We have uncovered a multiphasic molecular program where CT cell types found in the uninjured adult limb revert to a relatively homogenous progenitor state that recapitulates an embryonic limb bud–like phenotype including multipotency within the CT lineage. Together, our data illuminate molecular and cellular reprogramming during complex organ regeneration in a vertebrate.

Among tetrapods, only salamanders have the capacity to replace a lost limb (1). The adult axolotl (Ambystoma mexicanum) limb is composed of many different cell types originating from neural, myogenic, epidermal, and connective tissue (CT) lineages. Upon limb amputation, cells from areas near the amputation plane accumulate in a tissue called the blastema. Cells within the blastema are capable of fully regenerating the missing limb [reviewed in (2)]. CT cells, descendants of lateral plate mesoderm (LPM), are the most abundant lineage contributing to the blastema (35) and encompass bone and cartilage, tendons, periskeleton, and dermal and interstitial fibroblasts. The CT cells are a key cell type for deciphering the molecular program of regeneration, as they express factors that guide the regeneration of appropriate limb parts (4, 68). Nevertheless, how mature CT produces blastema cells has not been molecularly defined because of the inability to isolate and deconstruct this cell population. It is even unclear whether mature CT cells molecularly “reprogram” into embryonic cell–like limb progenitors or whether the CT harbors preexisting stem cells that selectively seed the blastema.

We generated Cre-loxP reporter lines to track CT compartments and uncover compartment-associated contributions to different limb segments. We used single-cell mRNA sequencing (scRNA-seq) to dissect CT heterogeneity in the blastema, as well as the adult, embryonic, and regenerated forelimb, enabling the characterization of CT cell types and lineage relationships as cells regenerated the limb. Notably, CT cells lose their mature phenotypes to form multipotent limb bud–like progenitors in the blastema. The combination of lineage tracking and single-cell transcriptomics resolves the origin and molecular profiles of redifferentiated CT cell types that emerge from the blastema. Our work provides a molecular view of individual cells that build a blastema and reconstitute a patterned limb skeleton.

A Cre-loxP reporter system tracks CT cells during axolotl limb regeneration

To label and track axolotl CT cells, we generated a germline transgenic line that expresses the tamoxifen-inducible Cre-ERT2 (Cre-ER) gene under the control of the Prrx1 limb enhancer element (Prrx1:TFPnls-T2A-Cre-ERT2;CAGGs:lp-GFP-3pA-lp-Cherry, hereinafter referred to as Prrx1:Cre-ER;Caggs:lp-Cherry animals) (figs. S1 and S2 and tables S1 and S2). Prrx1, a paired-related homeobox gene, is expressed in CT precursors in the developing limb bud and in the limb blastema (911). Immunofluorescence staining of PRRX1 protein confirmed specific expression in axolotl limb bud cells (Fig. 1A), blastema cells (Fig. 1B), and adult CT (fig. S2). Administration of tamoxifen to Prrx1:Cre-ER;Caggs:lp-Cherry animals at the limb bud stage resulted in an efficient (>80%) genetic labeling of adult limb CT (Fig. 1, C and D, and fig. S1E). Notably, after limb amputation, we found that Prrx1-expressing blastema cells express mCherry, showing that the transgenic reporter efficiently marks the adult precursors to the blastema cells (Fig. 1B). Examination of regenerates 25 days postamputation (dpa) revealed mCherry-expressing cells in upper and lower arm CT (Fig. 1D and fig. S1, C to F), showing that CT gives rise to new CT during regeneration. Therefore, this transgenic line provides a system to track CT cells during limb regeneration.

Fig. 1 Tracking and molecular profiling of axolotl limb CT.

(A) Longitudinal section of a stage 47 limb bud stained with anti-PRRX1 antibody (Ab) (red) identifying PRRX1 as a pan-CT marker during limb development. Arrowheads indicate the absence of PRRX1 staining in the epidermis. The boxed area in the top panel corresponds with the higher-magnification views in the bottom three panels. (B) Longitudinal section of a blastema 11 dpa, stained with anti-PRRX1 Ab (green). Red, converted cells; blue, Hoechst-stained nuclei. Scale bars, 500 μm. (C) Embryos after induction of Prrx1:Cre-ER;CAGGs:lp-Cherry by 4-OHT showing expression of mCherry only in limb mesenchyme. (D) Fluorescence image of converted cells in an uninjured limb and a regenerated limb (conversion at limb bud stage) indicating stable labeling of CT before and after regeneration. The arrowhead indicates the amputation plane. (E) (Left) tSNE plot visualizing scRNA-seq data for 2379 single cells (circles) from the adult axolotl upper arm. Dashed lines outline related cell types. (Right) mCherry expression is detected exclusively in CT cell types. (F) Bar plots showing the mean expression of marker genes in each cluster. The x axis represents cell clusters identified in (E). Error bars indicate SD. UMI, unique molecular identifier.

We used a high-throughput droplet-based scRNA-seq method (10× Genomics) to sample the cellular diversity in the uninjured adult limb and further validate this transgenic line. We converted cells at the limb bud stage and performed scRNA-seq on the dissociated uninjured adult limb tissue containing labeled and unlabeled cells (2379 cells) (table S3). Using unbiased clustering and the expression of marker genes, we identified endothelial, epidermal, immune, muscle, red blood, and CT cells (Fig. 1E). mCherry mRNA from Prrx1:Cre-ER;Caggs:lp-Cherry converted cells was detected only in the CT cluster, which included periskeletal, tendon, dermal, and fibroblastic cell subpopulations as identified on the basis of canonical marker expression (Fig. 1F). To specifically examine CT heterogeneity, we used t-distributed stochastic neighbor embedding (tSNE) clustering to analyze 2375 single-cell transcriptomes after the isolation of labeled Prrx1:Cre-ER;Caggs:lp-Cherry line–derived CT cells by fluorescence-activated cell sorting (FACS) from the adult upper forelimb (Fig. 2, A and B, and table S5). We identified eight distinct clusters that we assigned on the basis of marker gene expression as tenocytes (identified by the expression of marker gene Tnmd), periskeletal cells (with marker gene Col8a2), actively cycling cells (with marker gene Ccnb1), and five fibroblastic CT subpopulations (fCT I to V) (Fig. 2, B and C, and table S6). The two fibroblastic CT populations fCT IV and fCT V could be identified as components of the dermis on the basis of Twist2 and Ptgds expression. The cycling population contained multiple different cell types, including tenocytes, fCT I, fCT III, and fCT IV (fig. S4A). We detected only very few skeletal cells (fig. S4B), likely because of their inability to dissociate from the skeletal matrix. Because of the low cell number, skeletal cells fall within the periskeletal cluster in the tSNE plot. Previous live imaging data (12) and our own cell tracking data (fig. S9) show that skeletal cells do not contribute to the regeneration. Since regeneration of bone progresses via endochondral ossification, “skeletal cells” here refers to osteogenic and chondrogenic cells, the proportions of which differ between mature and regenerating bone. In summary, these profiling data validate the specificity of our Prrx1:Cre-ER;Caggs:lp-Cherry reporter animals, provide a cell atlas and a marker set for cell types of the uninjured adult axolotl limb (table S4), and characterize the heterogeneity of the upper arm CT (table S6).

Fig. 2 Blastema formation from axolotl upper arm CT cells involves molecular funneling during regeneration.

(A) Schematic of CT scRNA-seq experiments performed on mCherry+ CT cells of the uninjured axolotl upper arm (0 dpa) and the regenerating upper arm blastema (3, 5, 8, 11, and 18 dpa) from Prrx1:Cre-ER;Caggs:Lp-Cherry animals (conversion at 1-cm size). Cells were sorted by FACS. Amp., amputation. (B) A tSNE projection of 2375 single-cell transcriptomes visualizes the cellular heterogeneity of the uninjured upper arm CT and reveals eight clusters that refer to seven CT subtypes and one cluster of cycling cells (with expression of Ccnb1) (inset). (C) Violin plots showing distribution of expression for selected tSNE cluster marker genes (B). Colors refer to tSNE clusters (B). (D) Diffusion map projection (17) describing lineage relationships between uninjured CT cells, blastema-derived cells, and cells from a regenerated upper arm. DC, diffusion component. (E) DC 3 captures the cell type heterogeneity in the uninjured CT that is lost in the blastema. (F) Expression of cell type marker genes (groups i to vi) identified for the uninjured CT, shown as a heat map for uninjured CT and blastema time points, with genes in rows and cells hierarchically clustered in columns. Transcript levels are scaled across rows. Perisk., periskeleton. (G) Mean pairwise correlation (Pearson) between genes of each gene group (F) across all cells, calculated for each time point. Mean correlation coefficients decrease over the course of blastema formation. Error bars, SD. (H) Heat-map visualization of time point–specific marker genes (columns) with cells (rows) ordered by diffusion pseudotime (see also fig. S4I). GO enrichments (bottom) and exemplary genes (top) are shown (see also fig. S5A). Colored sidebar, time points. resp., response; act., activity; pos. reg., positive regulation; prol., proliferation; trans., transcription. (I) Pseudotemporal expression of different gene signatures across all cells from uninjured upper arm CT to 18-dpa blastema. Conditional means smoothed by using LOESS (local regression) are presented.

Blastema formation from axolotl upper arm CT cells involves molecular funneling during regeneration

To understand the molecular pathways involved in CT regeneration, we used a high-transcriptome-coverage scRNA-seq method (Fluidigm C1) to analyze mCherry-positive (mCherry+) cells isolated by FACS (Prrx1:Cre-ER;Caggs:lp-Cherry) along a dense time course that captures the major transitions during regeneration. Time points were determined on the basis of the average blastema cell cycle length (53 to 103 hours) (13, 14) and previous bulk transcriptional dynamics (15) and therefore included the uninjured upper forelimbs (108 cells), as well as blastema stages (3 dpa, 108 cells; 5 dpa, 167 cells; 8 dpa, 121 cells; 11 dpa, 163 cells; and 18 dpa, 135 cells) and a fully regenerated upper forelimb (3 to 12 months postamputation, 128 cells) (Fig. 2A, fig. S3, and table S7). We first analyzed the uninjured upper forelimb data and could confirm the presence of the CT subpopulations described above, which allowed identification of additional marker genes for each cell type (fig. S4, C and D). We next combined fibroblastic and periskeletal cells from the uninjured CT, which are the major populations contributing to the blastema (12), with all blastema states and the regenerated CT and used a diffusion pseudotime estimate (16, 17) to reconstruct the lineage path from uninjured CT through the blastema state to the regenerated CT (Fig. 2D and fig. S4E). This analysis revealed a general time-dependent progression, with some intermixing between time points. Diffusion component 3 visualizes how heterogeneous differentiated CT cells funnel into the relatively homogeneous early blastema state (Fig. 2E) and regenerated CT cells funnel out of the blastema to reestablish the initial cellular heterogeneity (fig. S4F). We confirmed this observation with two alternative analyses. When focusing on the cell type–specific expression patterns found in the uninjured tissue (Fig. 2F), we did not observe a comparable heterogeneity within the blastema cell populations and instead found that heterogeneity was diminished at all blastema time points (Fig. 2G) up until the redifferentiation of CT cell types (fig. S4G). Further, we calculated for each time point the number of unique genes detected within the top 100 most expressed genes per cell as a measure of intercellular heterogeneity per time point. We detected a decrease in the number of unique genes from the adult uninjured upper limb through the early blastema stages up to the 11-dpa blastema, followed by an increase at 18 dpa and in the regenerate (fig. S4H). This is consistent with a progressive decrease in intercellular heterogeneity during blastema formation until 11 dpa, with the highest heterogeneity found in the uninjured adult upper arm, the regenerate, and the 18-dpa blastema when cells begin to differentiate. Notably, we did not find blastemalike cells within the uninjured tissue because all cells from the adult upper forelimb cluster separately from blastema-derived cells in a tSNE analysis (fig. S4I). Also, when the cycling cells in the uninjured tissue are analyzed together with cycling blastema cells, the cells derived from the uninjured tissue are distinct because they express signatures of differentiated CT cell types (e.g., Mfap5 and Fbn1) and lack expression of blastema markers (e.g., Nrep and Prdx2) (fig. S4, J and K).

Our molecular data taken together with previous cell lineage observations (12) and digital tracking data (18) show little numerical bias in cell contribution during blastema formation, allowing us to conclude that blastema formation does not involve the selection of a preexisting stem cell population. With the cell cycle length (13, 14) and cell counts pointing to an expansion of cell numbers by an average of six times from the uninjured to the 11-dpa time point (supplementary materials and methods), a “preexisting” blastema cell would need to constitute at least 20% of the mature CT population to account for the cell expansion. In summary, these data indicate that mature CT cells dedifferentiate into a relatively homogeneous pool of progenitor cells when forming the blastema.

Next, we explored the transcriptional changes that ensue during the regenerative process within the blastema (Fig. 2H and fig. S5, A to C). The uninjured CT state is characterized by the expression of many extracellular matrix (ECM) genes that are down-regulated during blastema formation (Fig. 2, H and I). We observed an inflammation response (e.g., expression of Il11 and Cxcl2) in the early blastema (3 to 5 dpa) coinciding with the expression of various matrix metalloproteinase genes (Mmp8 and Mmp13), indicating extensive ECM remodeling induced by injury. By 8 dpa, cells express genes involved in proliferation (e.g., Ccnb3), and by 11 dpa, there is positive regulation of transcription (Hmga2 and Lbh expression) and evidence of repatterning (Hoxa13 and Hoxd12 expression). Reestablishment of ECM (via Matn4) is initiated by 11 dpa but peaks at 18 dpa. Notably, genes necessary for skeletal development (Hoxd13, Col2a1, and Sox9) are detected in the blastema 18 dpa. In addition, the expression of genes for various signaling molecules (Dkk1, Bambi, and Lep) associated with regeneration are induced throughout blastema development (15, 19, 20). These data constitute a molecular map of CT cells as they transition through the limb-regenerative process.

CT reprogramming progresses through a blastema-specific state before recapitulating embryonic limb bud development

We next explored the extent to which CT regeneration recapitulates limb development. To cover different stages of limb development, we examined the single-cell transcriptomes of stage 28 limb field cells (82 cells), as well as stage 40 (76 cells) and stage 44 (121 cells) limb bud cells (21) (Fig. 3A and table S7) and identified CT cells on the basis of Prrx1 expression (279 cells). By stage 44, the limb bud contains approximately 1500 cells. We observed that CT precursors in the developing limb are distinct from and less differentiated than adult CT (Fig. 3B and fig. S6A). Although cells from the two limb bud stages showed major similarities, cells from the stage 28 limb field were distinct from stage 40 and 44 limb bud cells, as shown by the unique expression of genes such as Wnt8a or Hoxd1 (fig. S6B). Variation in the expression of spatial patterning genes constituted the major source of heterogeneity in the limb buds (stages 40 and 44) (fig. S6C), whereas patterning genes were generally not yet expressed in stage 28 limb field cells (Fig. 3C). This variation enabled reconstruction of the proximal-distal and anterior-posterior axes through the correlated expression of patterning genes, confirming that our data captured the representation of the limb bud cells (Fig. 3D) (22). We investigated when patterning axes start to be established during limb regeneration and found that anterior-posterior markers (Fgf8 and Shh, respectively) and proximal-distal markers [Meis2 (proximal) and Hoxa11 and Hoxa13 (distal)] are established between 8 and 11 dpa (Fig. 3E). In a manner similar to that for the limb bud, patterning genes enabled the spatial reconstruction of cell positions in a blastema 11 dpa (Fig. 3F and fig. S6D).

Fig. 3 CT reprogramming progresses through a blastema-specific state before recapitulating embryonic limb bud development.

(A) Overview of scRNA-seq experiments on three axolotl limb bud stages. Limb bud CT cells (279) were in silico identified based on Prrx1 expression and compared with blastema cells on the transcriptome level. (B) Heat map showing expression of genes (columns) that distinguish mature limb CT cells from limb bud CT cells (rows, hierarchically clustered). The expression of CT cell type marker genes is shown on the right. (C) Bar graphs showing the fraction of cells per embryonic stage that express proximal-distal (Meis2, Hoxa11, and Hoxa13) or anterior-posterior (Fgf8 and Shh) patterning genes. (D) Intercellular correlation network of stage 44 limb bud cells (circles) placing cells in a hypothetical position within an imaginary limb bud on the basis of the expression of five known patterning genes (see also fig. S6C). (E) Limb bud patterning genes are reactivated during blastema formation. Bar graphs show the fraction of cells per blastema time point that express proximal-distal (Meis2, Hoxa11, and Hoxa13) or anterior-posterior (Fgf8 and Shh) patterning genes. (F) Intercellular correlation network of 11-dpa blastema cells (circles) placing cells in a hypothetical position within an imaginary limb blastema on the basis of the expression of five patterning genes (see also fig. S6D). (G and H) Box plots showing distributions of scaled correlation between single-cell transcriptomes at any given time point and the mock bulk transcriptome of stage 44 limb bud (G) or stage 28 limb field (H) CT cells. Limb bud and limb field progenitors are most similar to 11-dpa blastema cells. (I) Scatter plot showing differential correlation of single-cell transcriptomes (dots; color represents time point) with limb bud versus mature CT transcriptomes (y axis) and with 5-dpa blastema versus 11-dpa blastema transcriptomes (x axis). (J) Dot plot visualizing the expression of genes shared between the 11-dpa blastema and limb bud progenitor cells. The circle size represents the fraction of cells at each time point expressing the gene, and color represents the average expression level.

We next quantified the similarity between different blastema stages and limb development. We created mock bulk transcriptomes for each uninjured, regenerate, embryonic, and blastema time point and calculated the correlation (Spearman’s r) of each single cell with each bulk transcriptome. We found that the correlation of blastema cells with stage 40 and stage 44 limb buds peaks at 11 dpa, with a progressively weaker correlation at earlier blastema time points (Fig. 3G and fig. S6, E and F). Notably, we observed the same trend when we compared blastema cells with the stage 28 limb field (Fig. 3H). Additionally, we found that cells in the 3- and 5-dpa blastemas are similar neither to uninjured adult cells nor to limb bud or limb field cells, suggesting that a cell state emerges in the blastema that is distinctive for limb regeneration (Fig. 3I and fig. S6E). Within our dataset, we identified many genes that were expressed only in the early blastema stages, whereas 11-dpa cells largely share expression patterns with the cells derived from the developing limb (Fig. 3J and fig. S6, G and H). Our data suggest that the former CT cells within the blastema initiate reprogramming by using genetic mechanisms that are distinct from embryonic limb bud development but arrive at a limb bud–like state by day 11.

Tracking different CT subpopulations in the blastema reveals distinct cell sources for proximal versus distal limb regenerate tissue

We next sought to analyze how different CT cell subpopulations contribute to redifferentiated cell types in the regenerated limb. We first generated a new transgenic line by using the Col1a2 promoter (Col1A2:TFPnls-T2a-ERT2-Cre-ERT2;CAGGs:lp-GFP-3pA-lp-Cherry, hereinafter referred to as Col1A2:ER-Cre-ER;Caggs:lp-Cherry) in which only a subset of CT cells are labeled and tracked (fig. S7). Specifically, conversion of 3-cm larvae leads to limbs with genetic labeling primarily of the skeleton, periskeleton, and tendons but very few dermal fibroblasts and no interstitial fibroblasts, thus constituting a subset of cell types labeled in the Prrx1:Cre-ER;Caggs:lp-Cherry line (Fig. 4A and fig. S7A). scRNA-seq on 36 sorted, labeled cells from the uninjured adult upper arm revealed periskeletal cells and tenocytes (bone cells were not recovered in the dissociation, but they do not contribute to the blastema), confirming the histological analysis of the labeled cells (Fig. 4B and table S8). We next performed scRNA-seq on labeled Col1a2:ER-Cre-ER;Caggs:lp-Cherry (referred to as Col1a2) descendants in the 11-dpa blastema (349 cells) and could identify undifferentiated progenitors as well as immature skeletal and nonskeletal cell lineages (Fig. 4C and table S8). A pseudotemporal trajectory displayed two branching paths: one transitioning from progenitors (expressing Mycl, Matn4, and Nrep) to nonskeletal cells (expressing Tnmd, Fcn2, and Aspn; likely precursors to tenocytes and periskeletal cells) and the other to skeletal cells (expressing Cnmd/Lect1 and Otos) (Fig. 4D, fig. S6I, and table S10).

Fig. 4 Tracking CT subpopulations in the blastema reveals distinct cell sources for proximal versus distal limb regenerate tissue.

(A to G) Col1a2:ER-Cre-ER;Caggs:lp-Cherry–labeled animals were used. (A) (Top) Bright-field image with fluorescent overlay showing limbs at 0 dpa. (Bottom) Upper arm (UA) limb cross section revealing periskeletal and tendon cell labeling (Hoechst, blue; mCherry, red). Scale bars, 2 mm. (B) Heat map showing marker expression for 36 mature limb periskeletal and tendon cells (hierarchically clustered; Pearson). (C) Heat map showing scaled expression (exp.) of genes identified by PCA (table S10) in 349 Col1a2 line–derived 11-dpa blastema cells. (D) Diffusion map projection of Col1a2 line–derived 11-dpa blastema cells with signature scores shown. (E) (Top left) Image overlay of 25-dpa limbs (scale bar, 2 mm). The arrow shows the site of amputation. (Top right, bottom left, and bottom right) Limb cross sections along the proximal-distal axis. Samples are stained with Hoechst (blue), and images are overlaid with DIC and the mCherry fluorescence signal (red) of converted cells. The arrowhead indicates callus formation. Scale bars, 200 μm. (F) Fraction of converted cells in five CT subtypes at different proximal-distal positions after regeneration (n > 8000 cells, three limbs). reg., regenerate; fibro., fibroblasts; iF, interstitial fibroblasts (fCT I to III). (G) mCherry-labeled humerus transplantation into an unlabeled host. (Top) Images of live animals. Scale bars, 1 mm. (Center) Longitudinal sections of 15-dpa blastema showing converted cells (red), Hoxa11 expression (green), and Hoechst staining (blue). Scale bar, 500 μm. (Bottom) Colocalization (arrowheads) of mCherry+ and Hoxa11+ cells. (H to L) Prrx1:Cre-ER;Caggs:lp-Cherry–labeled animals were used. (H) mCherry-labeled humerus transplantation into an unlabeled host. (Left to right) Images of live animals (scale bars, 1 mm) and limb cross sections from UA-callus, UA-regenerate, and LA-regenerate positions (15 dpa; converted cells, red; Hoechst, blue; scale bars, 200 μm). (I) Quantification of converted cells among periskeletal and skeletal cells and an aggregate of both subtypes (All) along the proximal-distal axis. cal., callus. (J) Unlabeled humerus transplantation into an mCherry converted host. (Left to right) Images of live animals (scale bars, 1 mm) and limb cross sections from UA-callus, UA-regenerate, and LA-regenerate positions (15 dpa; converted cells, red; unconverted cells, green; Hoechst, blue; scale bars, 200 μm). The magnification shows mCherry+ periskeletal (Ps) and skeletal (Sk) cells and interstitial fibroblasts (iF) (fCT I to III). (K) Quantification of converted cells among periskeletal and skeletal cells and an aggregate of both subtypes along the proximal-distal axis. Error bars in (F), (I), and (K) indicate SD. (L) Diffusion map projection of Prrx1 line–derived 18-dpa blastema cells with signature scores shown. The map was created by using genes identified by PCA (table S10).

Microscopic examination of regenerates revealed a spatial bias in cell contribution. Labeled cells were found primarily in the extension to the existing bone at the amputation site (humerus), with few cells in more distal segments (Fig. 4, E and F). We also found a similar spatial restriction of Col1a2 descendants to extending bone in lower arm amputations (fig. S8B). This is in contrast to the Prrx1:Cre-ER;Caggs:lp-Cherry line in which all proximal and distal positions of the regenerate are labeled. These data provide direct evidence supporting a previous hypothesis that distinct cell sources are used for bone extension versus de novo segment formation (23). The Col1a2 line–derived blastema cells could be intrinsically or extrinsically limited to extending existing bone. We transplanted a humerus from a Col1a2:ER-Cre-ER;Caggs:lp-Cherry donor into an unlabeled host and found that, after upper arm amputation, some rare mCherry+ cells were found in the distal, Hoxa11+ region of the blastema 15 dpa and showed Hoxa11 staining (Fig. 4G). These observations lead us to conclude that the Col1a2 line–derived blastema cells are not intrinsically limited in their segmental identity but rather are strongly associated with the callus and therefore spatially biased toward extending their bone of origin.

We performed multiple transplantation experiments to confirm the above-mentioned observations and resolve the source of the distal CT. To exclude the possibility that the Col1a2 driver marked only a subset of periskeleton and tendon with spatially restricted potential, we grafted a humerus from Prrx1:Cre-ER;Caggs:lp-Cherry (referred to as Prrx1) converted limbs into an unlabeled host limb, replacing the host humerus before upper arm amputation (Fig. 4, H and I). As in the Col1a2 tracing, mCherry+ cells were found in the callus and the regenerated periskeleton and skeleton just beyond the amputation site, with progressive depletion toward the lower arm regenerate. We next performed a complementary graft of unlabeled humerus into Prrx1 converted hosts to label nonbone CT. We found nearly complete labeling of tendon, skeleton, and periskeleton as well as dermal and interstitial fibroblasts in the lower arm regenerate, indicating that nonskeletal CT cells regenerate the distal CT (Fig. 4, J and K). We separately grafted a muscle fiber bundle containing labeled interstitial fibroblasts from embryonic Prrx1 converted animals into unlabeled hosts and found extensive labeling of the distal CT (fig. S10). These data indicated that the fibroblastic CT cells were the major contributors to the distal limb regenerate. To analyze the regeneration of distal limb CT on the molecular level, we performed a lineage reconstruction analysis on the scRNA-seq data of labeled Prrx1 descendants in the 18-dpa blastema (Fig. 4L, fig. S6J, and table S10). This analysis revealed a bifurcated path where uncommitted blastema progenitors branch off into a nonskeletal and a skeletal lineage. At the 18-dpa time point, differentiated dermal and interstitial fibroblasts were not yet identifiable, consistent with live imaging data showing the late differentiation of these lineages (12). Notably, both Col1a2 descendant blastema cells (Fig. 4C, mostly periskeleton derived) and Prrx1 descendant blastema cells (Fig. 4L, mostly fibroblastic CT derived) had an expression profile resembling that of limb bud progenitors. Taking these results together, we conclude that cells from multiple CT compartments funnel into an undifferentiated progenitor, with the tissue of origin biasing the spatial contribution of the cells.

CT lineages reemerge through multipotent progenitors

To reconstruct the reestablishment of each CT cell type in the upper arm regenerate, we performed high-throughput droplet-based scRNA-seq of labeled Prrx1:Cre-ER;Caggs:lp-Cherry descendants in the late blastema at three time points, including 18 dpa (9939 cells), 25 dpa (9019 cells), and 38 dpa (2861 cells) (Fig. 5A and table S9). We inferred differentiation trajectories and pseudotemporal cell relationships (24) from diffusion map embedding (17) and identified a trifurcated path, where multipotent blastema progenitors expressing embryonic limb and cell cycle markers (e.g., Prdx2, Nrep, and Ccnb1) branch off into a nonskeletal lineage or a skeletal lineage that then bifurcates into either cartilage or bone (Fig. 5, B and C). We observed temporal differences in the lineage progression, as progenitor cells are present only at two earlier time points (18 and 25 dpa), whereas cells on the nonskeletal branch are found at all three time points. A skeletal precursor state is found mainly at 18 dpa, before cartilage and bone start to differentiate at 25 dpa and 38 dpa, respectively (see also fig. S11B). We next analyzed the heterogeneity in the nonskeletal branch and found that the main nonskeletal lineages (periskeletal cells, tenocytes, and fibroblastic CT cells) present in the adult uninjured tissue had started to reemerge by 38 dpa (Fig. 5D and fig. S11, C and D). These data suggest that, at the transcriptome level, the progenitor pool between 18 and 25 dpa is still relatively homogeneous and that these progenitor cells diversify sometime after 25 dpa into diverse CT lineages.

Fig. 5 CT lineages reemerge through multipotent progenitors.

(A) High-throughput scRNA-seq was performed on mCherry+ CT cells of the uninjured axolotl upper arm (0 dpa) and the upper arm blastema (at 18, 25, and 35 dpa) from Prrx1:Cre-ER;Caggs:Lp-Cherry converted animals. Cells were sorted by FACS. (B) Diffusion analysis (17) and tree construction (24) of blastema time points (18, 25, and 38 dpa) identify four main branches. Pie charts show time point proportions per branch. (C) Pseudotemporal marker gene expression for each branch. (D) SPRING analysis (48) of the nonskeletal blastema branch cells (blue inset) together with mature CT cells (a total of 3151 cells) reveals the reemergence of CT subpopulations identified in the mature tissue. (E) SPRING plots colored by early blastema (Nrep), cell cycle (Ccnb1), and CT subtype (Tnmd, Col4a2, and Twist2) marker expression show the cell differentiation during late phases of regeneration.

We further sought to validate our molecular analysis by clonal lineage tracing using a Brainbow transgenic animal, which allowed us to determine the spectrum of cell types formed from single CT cells. To establish clonal tracing conditions, recombination was induced in mature limb cells and examined after 7 days. To identify color combinations showing appropriately low occurrence for clonal analysis, clone pairs were identified as infrequently occurring adjacent sisters derived from a cell division, and their distribution in color and hue-saturation (HS) space was determined (fig. S11E). From these data, a rare recombinant type that mapped in the blue color range was identified (fig. S11F). Examination of nine different source zone equivalents in the mature tissue identified 2, 1, 0, 1, 1, 0, 0, 2, and 0 blue cells per zone (fig. S11G). Amputated limbs were allowed to regenerate for 25 days (Fig. 6A) and then examined for the presence of blue clones. Three out of 10 limbs contained a blue clone, whereas the other seven limbs showed no blue clone, confirming its rare occurrence (Fig. 6, B and D). The frequency distribution in color space of the blue cells in each regenerate was consistent with clonal origin (Fig. 6, C and E; for more examples, see fig. S11, H and I). Assignment of the blue cells from within each limb to CT subtypes revealed contributions of clonal descendants to skeletal, periskeletal, fibroblastic, and tendon cells (Fig. 6B). These lineage tracing data confirm the molecular profiling conclusions that limb blastema cells acquire a limb bud progenitor identity and form a multipotent CT progenitor (25, 26).

Fig. 6 Brainbow clonal analysis confirms multilineage potential of CT cells upon limb regeneration.

(A) Image of a regenerated Brainbow axolotl limb. A presumptive clone of blue cells is observed throughout the limb, from the digit tip (A′), the elbow (A′′), and the amputation plane (solid line) at the injection site (arrowhead) (A′′′). Scale bars, (A) 300 μm; (A′ to A′′′), 100 μm. (B) HS color distribution of cells from a representative image, including the presumptive clone of blue cells (white circle). Multiple cells of each CT subtype are represented. Similar distributions were observed in 3 of 10 analyzed samples (for examples, see fig. S11H). (C) Frequency distribution in HS color space, calculated by using the formula in fig. S11E, for known clonally related cells (fig. S11, F and G) (triangles), presumptive clonally related cells [the blue cells in the white circle in (B)] in a regenerated limb (circles), and non–clonally related cells in a regenerate (squares). The frequency distribution of the presumptive clone is indistinguishable from that of known clones (Kruskal-Wallis analysis and Dunn’s multiple comparison). (D) Example of HS color distribution of cells from a representative image lacking a discreet cluster of blue cells (white circle). Similar distributions were observed for 7 of 10 analyzed samples (for examples, see fig. S11I). (E) Frequency distribution as in (C) for the sample shown in (D). Because of the lack of identification of a clonally related subset, no presumptive clone could be mapped.

Summary

The molecular understanding of blastema formation was previously limited by the inability to identify and isolate blastema precursor cells in the adult tissue. We have demonstrated the importance of genetically marked transgenic axolotl strains for isolating blastema precursor cells from adult limb tissue, and we have molecularly profiled these cells by using single-cell transcriptomic methods. This profiling has indicated that CT cells express adult phenotypes that are lost upon the induction of regeneration and funnel into an embryonic limb bud–like phenotype that includes multipotency within the CT lineage (25, 26). The molecular reprogrammability of adult cells to cells of embryonic limb potential capable of orchestrating complex limb morphogenesis has clear implications for future prospects in regenerative engineering.

Materials and methods

Axolotl husbandry, transgenesis, and 4-OH tamoxifen treatment

Axolotls (Ambystoma mexicanum) were bred in MPI-CBG, CRTD, and IMP facilities. All animal handling and surgical procedures were carried out in accordance with local ethics committee guidelines. Animal experiments were performed as approved by the State Authorities Saxony and the Magistrate of Vienna. “White” refers to a nontransgenic d/d strain of axolotl that has white skin due to the absence of melanocytes. Animals were anesthetized in 0.03% benzocaine (Sigma) before amputation and surgery. Reference (27) describes axolotl husbandry, transgenesis, and 4-hydroxy tamoxifen (4-OHT) treatment in detail.

For generating transgenic animals, Prrx1 (10), Col1A2 (28), and Col2A1 (29) driver elements were cloned at the 5′ end of the TFPnls-T2A-ERT2-Cre-ERT2 (ER-Cre-ER) or TFPnls-T2A-Cre-ERT2 (Cre-ER) cassette with flanking SceI sites. TFPnls refers to nuclear teal fluorescent protein. The Col1A2 1.5-kb promoter and the Prrx1 2.4-kb enhancer/promoter from the mouse are kind gifts from George Bou-Gharious and Malcolm Logan, respectively. Constructs used for generating transgenic axolotl driver lines are available at Addgene (111150, 111151, 111152). Caggs:LoxP-eGFP-3polyA-LoxP-Cherry (Caggs:lp-Cherry) transgenic animals and the Col2A1 promoter were described in previous publications (29, 30). Transgenic driver animals were generated as described previously by using Sce1 meganuclease (27). General information on the founder animals can be found in table S1. In order to perform tracing, driver lines were crossed with Caggs:lp-Cherry reporter animals to obtain double transgenic progeny. We believe that our Caggs:lp-Cherry reporter line carries multiple insertion copies of the cassette and thus produces offspring with either one or more than one copy of Caggs:lp-Cherry. Occasionally, we obtained animals that showed better conversion efficiency than the rest of the animals. Thus, for each experiment, animals were carefully chosen to ensure that they had a similar level of green fluorescent protein (GFP) intensity before conversion. All tracing and transplantation experiments involving transgenic animals were performed with F1 or later generations.

4-OHT treatment was done as described previously (27) and as presented in table S2. Briefly, small animals were treated with 1 to 2 μM 4-OHT by bathing, and large animals were intraperitoneally injected with 5 μl of 10-mg/ml 4-OHT per gram of body weight of the animals.

Immunohistochemistry and microscopy

Limb cryosections of 10 μM thickness were rehydrated and permeabilized by using phosphate-buffered saline (PBS) with 0.3% Tween-20 and then blocked for 1 hour with PBS with 0.3% Triton-X100 and 2% normal horse serum (Vector Laboratories # S-2000). Slides were then incubated overnight with primary antibodies in blocking buffer in a humidified chamber. The next day, slides were washed with PBS–0.3% Tween-20 and then incubated with secondary antibodies for 2 hours, washed with PBS–0.3% Tween-20, and mounted with mounting medium (Vector Laboratories # H-1000). The cell nuclei were stained with Hoechst 33342 (Sigma; final concentration, 0.5 μg/ml). Mosaic images of sections were acquired on a Zeiss Axio-observer.z1 by using a 20× apochromatic lens with 0.8 numerical aperture (NA) and Axiovision software or Zen2 (Zeiss).

Antibodies

The primary antibodies used in these studies were anti-COL1A2 (DSHB, Sp1.D8), PAX7 (MPI-CBG antibody facility), MHC (monoclonal antibody 4A1025, a kind gift from S. Hughes), β-3TUB (R&D#MAB1195), MBP (Genetex#GTX761141), anti-HOXA11 (MPI-CBG antibody facility), and anti-PRRX1 antibodies (MPI-CBG antibody facility). Secondary antibodies used in these studies were procured from Molecular Probes.

Anti-PRRX1 antibodies

For raising polyclonal anti-PRRX1 antibodies, a DNA fragment of the axolotl Prrx1 coding for the N-terminal amino acids 1 to 101 was cloned to generate a glutathione S-transferase (GST) fusion protein (GST-AxPRRX1-N) and a maltose binding protein (MBP) fusion protein (MBP-AxPRRX1-N). GST-AxPRRX1-N and MBP-AxPRRX1-N fusion proteins were expressed in Escherichia coli (BL21-DE3-pLysS) and purified by using glutathione Sepharose 4B (GE) and amylose resin (NEB), respectively. GST-AxPRRX1-N was injected into the rabbit to generate polyclonal antibodies. The serum was affinity purified against MBP-AxPRRX1-N that was immobilized on NHS-activated Sepharose 4 Fast Flow (GE). Antibodies were dialyzed and concentrated by using Amicon Ultra-15 10K MWCO (Millipore) before usage.

To validate anti-PRRX1 antibodies, we first immunostained sections of a developing axolotl limb bud. As expected for a transcription factor, we noticed strong staining in the nuclei of limb mesenchyme and complete absence of staining in the epidermal cells (Fig. 1A). To test whether anti-PRRX1 can stain CT cells of the mature limb, we performed an LPM transplantation from GFP donor embryos (discussed below). Staining of such LPM-transplanted limbs showed complete colocalization among GFP+ cells and PRRX1+ stained nuclei (fig. S2A). To further validate the specificity of PRRX1 antibodies to CT, we co-stained PRRX1 with markers of muscle fibers (MHC), muscle satellite cells (PAX7), Schwann cells (MBP), and neurons (β-3TUB) in mature limb sections (fig. S2, B to E). All of these markers showed no colocalization with anti-PRRX1 antibodies, suggesting that anti-PRRX1 antibodies specifically label the CT population in the limb bud and the mature limb. Further, as expected, we found an increase in anti-PRRX1 staining in both upper arm and lower arm blastema upon amputation (Fig. 1A). This suggests that PRRX1 is a pan-CT marker for the limb.

Identification of five subtypes of CT

Collagen 1 alpha 2 (COL1A2) is an ECM protein, and in axolotl, it is expressed in dermal fibroblasts and skeleton cells (28). Staining limb sections with COL1A2 antibodies showed specific labeling of basal lamina, tendons, and skeletal and periskeletal structures. Thus, on the basis of co-staining of PRRX1 and COL1A2, along with differential interference contrast (DIC) imaging (which delineated the tissue morphology), we identified the following five distinct CT subtypes (fig. S8): (i) dermis/dermal fibroblasts (dF), cells underneath the basal lamina that are COL1A2+ and PRRX1med; (ii) interstitial fibroblasts (iF) (fCT I to III), cells scattered across the limb in or around muscle fibers, nerve fibers, and blood vessels, etc., that are COL1A2 and PRRX1high or PRRX1med; (iii) tenocytes (Tn), cells found as thick bundles of CT cells with elongated nuclei in the interstitial space that are COL1A2+ and PRRX1med; (iv) periskeletal cells (Ps), cells at the periphery of the skeleton that are COL1A2+ and PRRX1high or PRRX1med; and (v) skeletal cells (Sk), cells that reside inside a cartilaginous structure and are COL1A2+ and PRRX1low.

Characterization of transgenic animals

Prrx1:Cre-ER animals showed expression of TFPnls in limb buds, with fluorescence progressively declining as the limb developed (fig. S1B). In fully developed limbs, TFPnls is no longer visible under a stereoscope. 4-OHT treatment of Prrx1:Cre-ER;CAGGs:lp-Cherry animals at stage 44 (21) (just after hatching) resulted in conversion in limb bud mesenchyme (Fig. 1C). In the developed limbs, these converted cells contributed to the entire CT lineage (Fig. 1D and fig. S1).

In Col2A1:ER-Cre-ER animals, TFPnls fluorescence was visible in cartilage structures of the primary body axis, mainly the jaw skeleton, but not in the secondary body axis (limb skeleton). However, 4-OHT treatments of Col2A1:ER-Cre-ER;CAGGs:lp-Cherry animals at an early hand morphogenesis stage showed conversion not only in the skeleton of the primary axis but also in the limb (fig. S9).

In Col1A2:ER-Cre-ER animals, TFPnls fluorescence was visible in skin and cartilage structures (fig. S7A). In a longitudinal section of a Col1A2:ER-Cre-ER limb, TFPnls showed complete colocalization with anti-COL1A2 antibodies in dermal fibroblasts, tendons, and periskeletal and skeletal cells but not in the interstitial fibroblasts (fig. S7B). Similarly, conversion of Col1A2:ER-Cre-ER;CAGGs:lp-Cherry animals with 4-OHT showed conversion in dermal fibroblasts, tenocytes, periskeletal cells, and skeletal cells but not in the interstitial cells (fig. S7, B and C).

LPM transplantation and skeleton-LPM transplantation

For LPM transplantation, LPM cells from stage 18 Caggs:eGFP embryos were transplanted onto white embryos (4). Clean LPM transplantation results in CT labeling. Animals with skeleton only were obtained in this process, probably by an incomplete LPM transplant.

Humerus transplantation and interstitial CT transplantation

For humerus transplantation, 5- to 7-cm-long Prrx1:Cre-ER;CAGGs:lp-Cherry embryonic converted transgenic animals and white animals were used. Limbs of both transgenic and white animals were amputated at the elbow. All soft tissues including muscle, nerve, blood vessels, and skin were gently pushed back by using forceps to free the humerus from any attachment, and then the humerus was pulled out. Humeri of transgenic and white animals were swapped and were placed into the corresponding cavity. Later, tissues were trimmed to the distal metaphysis of the upper arm to mimic an upper arm amputation and produce a flat amputation surface. Two kinds of humerus-transplanted animals were obtained: (i) Prrx1:Cre-ER;CAGGs:lp-Cherry hosts with unlabeled white humeri and (ii) unlabeled white hosts with Prrx1:Cre-ER;CAGGs:lp-Cherry–labeled humeri.

Similarly, Col1A2:ER-Cre-ER;CAGGs:lp-Cherry–labeled humeri were grafted into the nonlabeled white host to obtain preferential labeling in the skeleton and periskeletal cells. Such limbs were harvested 15 days later for the analysis of HOXA11 staining.

For interstitial CT transplantation, 5-cm-long Prrx1:Cre-ER;CAGGs:lp-Cherry limb bud converted transgenic animals were used as donors. From the amputated limbs from donors, the full thickness of skin was removed and discarded. Next, humeri were separated from soft tissues and discarded. Remaining soft tissues, including labeled interstitial CT cells, were placed in 0.8× PBS. In parallel, white hosts underwent amputation near the distal metaphysis, and humeri were trimmed to produce a flat surface. Some muscle tissues were carefully removed and replaced with the soft tissue obtained from the donor to obtain interstitial CT–transplanted animals.

Five-centimeter-long animals were used for all three transplantation procedures. After transplantation, the hosts were kept anesthetized for another 2 hours in a humidified chamber with 0.005% benzocaine to allow healing. Animals were then carefully transferred to fresh water.

Tracing experiments and analysis

Tracing experiments were carried out with 4- to 5-cm-long animals (table S2). For tracing experiments, animals were amputated at the distal humerus near the metaphysis, and humeri were trimmed to produce a flat amputation surface. Animals were imaged on an Olympus SZX16 fluorescent stereoscope every 5 days to monitor regeneration. At 25 dpa, regenerated limbs were harvested by amputation near the shoulder and fixed in 4% MEMFA overnight at 4°C. The next day, limbs were washed three times with PBS, dehydrated in 30% sucrose, and embedded in Tissue-tek O.C.T. compound (Sakura). Limbs were sectioned and stained as described earlier. Cross-section analysis and identification of each subtype were done on the basis of DIC morphology and immunostaining of PRRX1 and COL1A2 as described earlier. In all cases, limbs were analyzed at four proximal-distal positions: namely, UA-mature (upper arm, mature), UA-callus (upper arm, transition zone), UA-regenerate (upper arm, regenerate), and LA-regenerate (lower arm, regenerate).

For Col2A1:ER-Cre-ER;CAGGs:lp-Cherry tracing, converted cells were counted only in the skeletal compartment because no converted cells were observed in any other CT subtype. For humerus-transplanted animals, counting was done in periskeletal and skeletal compartments. For Col1A2:ER-Cre-ER;CAGGs:lp-Cherry and Prrx1:Cre-ER;CAGGs:lp-Cherry (conversion in the mature limb) tracing, converted cells in all five subcompartments were counted. For Col1A2:ER-Cre-ER;CAGGs:lp-Cherry and Prrx1:Cre-ER;CAGGs:lp-Cherry (conversion in the mature limb) analysis, cross sections were made on sister slides and stained with COL1A2 and PRRX1 antibodies to identify the presence of total CT cells and converted CT cells in each compartment. For Prrx1:Cre-ER;CAGGs:lp-Cherry (conversion in the mature limb) analysis, fold enrichment was calculated as a ratio of percentages of labeled cells in the regenerate over the percentage of soft CT in the mature upper arm limb (preamputation zone). Statistical analyses were performed by using GraphPad Prism 6.0 (GraphPad Software).

Axolotl cell cycle kinetics

The somatic axolotl cell cycle has been measured in the range of 53 to 103 hours (13, 14), which is related to the massive genome (31). This means that by 11 days of regeneration, a maximum of 2.5 to 5 cycles have occurred, which translates to a 6- to 32-fold expansion in cell number (see below for actual cell counts). After limb amputation, cells within several millimeters respond to the injury by initiating S phase. Previous bulk molecular profiling comparing wounding with limb amputation across very fine time points determined that the common injury response lasts for 3 days (15).

Blastema founding population

Two important pieces of information came from previous live imaging of bulk CT cells during digit tip regeneration with the use of Brainbow transgenics (12):

(i) The majority of CT cells (except cartilage) migrate into the blastema, excluding the selection of a rare cell type and validating the FACS approach.

(ii) Cells within 500 μm migrate into the blastema, allowing us to calculate the founding population.

Cell numbers in the axolotl blastema and the founding population

The 11-day blastema has on average 45,000 PRRX1+ cells, whereas a 500-μm slice of the mature limb harbors on average 7500 PRRX1+ cells, excluding cartilage. This implies, on average, a sixfold expansion of cell number during regeneration, which is consistent with an average 103-hour cell cycle. The cell cycle numbers combined with these cell counts imply that an embryonic cell–like precursor would have to be present conservatively at 20 to 100% of the mature cell population.

Brainbow clonal analysis

For clonal tracing, 5-cm-long Brainbow axolotls (12) were injected with Tat-Cre protein to indelibly label cell lineages with different fluorescent proteins. Limb amputations were performed 7 days postinjection, just in front of converted cells. Limbs were allowed to regenerate for 25 days and then fixed by using MEMFA. Limbs were cleared by using a 1-propanolpH9/Eci clearing approach (32). Z-stacks were acquired on an inverted Zeiss LSM780 equipped with a 10×/0.3 EC Plan-Neofluar objective (Carl Zeiss Microscopy GmbH, Germany). Zen 2.3 SP1 (Zen black/64 bit) was used for image acquisition and automatic stitching of images. Representative RGB slices were acquired by using FIJI (33) and processed as described previously (34) with minor changes. Hue and saturation values were extracted by using the 5 by 5 Average option for the Magic Wand tool in Adobe Photoshop CS6. Clonal analysis was performed essentially as described previously (35) with modifications. To determine the distance in color space between known clones, control limbs were left unamputated and surveyed for distinct pairs of cells with similar color, location, and morphology. These cells are assumed to be clonally derived from each other and provide a measure for the expected variance in color space for clonally related cells. To determine the expected prevalence of clonal distribution in color space, unamputated control samples were divided into 350-μm regions and analyzed for clonal diversity. These regions represent the 500-μm area that contributes all the cells required for limb regeneration after taking into account the clearing-induced shrinkage of ±30 to 40%. Regenerated limbs were analyzed by selecting both presumptive clones and non–clonally related cells, while recording cell type, hue, and saturation. Kruskal-Wallis analysis and Dunn’s multiple comparisons test were used to determine differences in clonal analysis. Frequency distributions were used to graphically display samples. Outliers were removed by using ROUT (Q=1%). Three unamputated control samples and 10 regenerated samples were used for analysis. Statistical analysis was performed using Prism7 (version 7.0c).

Dissociating and preparing axolotl upper arm tissue for scRNA-seq experiments

Six to 10 blastemas, limb buds, or uninjured upper forelimbs were dissociated as a pool for each scRNA-seq experiment. Blastema samples and uninjured tissues were dissociated in 500 μl of 1× Liberase (Roche) diluted in 0.8× PBS (PBS without Mg2+ or Ca2+) supplemented with 0.5-U/μl deoxyribonuclease I (DNase I) (Thermo Scientific) at room temperature for ~30 min. Mechanical disruption by forceps as well as smooth pipetting was applied during incubation. Then, the cell suspension was filtered through a 30-μm-diameter strainer to generate a single-cell suspension. The filter was washed three times with 500 μl of 0.8× PBS, and either the obtained cell suspension was immediately subjected to FACS or cells were collected by centrifugation at 300 relative centrifugal force (rcf) for 5 min and again suspended/washed in 1 ml of 0.8× PBS for the 10× Genomics experiments. Centrifugation and washing were repeated once more, and cells were afterward again filtered through a 30-μm-diameter strainer. Cells were collected in 50 μl of supernatant and counted manually. For C1 experiments, cells were sorted by FACS in a single tube and centrifuged at 300 rcf for 5 min, after which they were resuspended in ~15 μl of 0.8× PBS and then kept on ice. Cells were counted manually. Limb bud samples were dissociated in 100 μl of 1× Liberase (Roche) in 0.8× PBS with disrupting and smooth pipetting at room temperature for 30 min. Afterward, the cell suspension was diluted with 200 μl of 0.8× PBS and filtered two times through a 30-μm-diameter strainer to generate a single-cell suspension that was manually counted and directly used for single-cell experiments. Limb field tissues were dissociated in 0.8× PBS while EDTA was added until an easy dissociation with pipette tips could be achieved.

Design of scRNA-seq experiments

As an uninjured reference state and to illuminate general axolotl limb heterogeneity, scRNA-seq experiments using the 10× Genomics platform were performed without any cell type enrichment. In contrast, to obtain deeper insights into CT heterogeneity, scRNA-seq experiments with the uninjured mature upper forelimb were performed by using FACS enrichment. To investigate the process of blastema formation, we collected independent scRNA-seq blastema data along a time course at 3, 5, 8, 11, and 18 dpa of axolotl upper arm forelimbs. In all CT-specific experiments, FACS was used to enrich for mCherry+ CT cells from converted Prrx1:Cre-ER;CAGGs:lp-Cherry animals. Around 40,000 cells were sorted for each single-cell experiment into 300 μl of 0.8× PBS, and at least two microfluidic C1 chips (Fluidigm, 96 capture sites) were used per experiment. For stage 40 and stage 44 limb bud experiments, dissociated cells were directly used for single-cell isolation on microfluidic chips without FACS enrichment. Wild-type as well as Prrx1:Cre-ER;CAGGs:lp-Cherry animals were used for limb bud experiments, and at least two microfluidic C1 chips (Fluidigm, 96 capture sites) were used for each limb bud stage. For the second set of scRNA-seq experiments where we used converted Col1A2:ER-Cre-ER;CAGGs:lp-Cherry animals, individual mCherry+ cells were directly sorted by FACS in single wells of a 96-well plate, and SmartSeq2 (36) was performed.

Capturing of single cells and preparation of cDNA

Prrx1 converted animals were used for all 10× Genomics experiments. FACS enrichment of the CT was performed for two of the uninjured CT experiments as well as for the blastema samples (18, 25, and 38 dpa). Two experiments on whole upper arm sections were performed without CT enrichment to analyze the complete axolotl limb heterogeneity. Cells were loaded onto the 10× cartridge at various concentrations, with a target of 2000 to 6000 cells per sample depending on the input cell concentration and the total cell number obtained after tissue dissociation. Cell encapsulation, cDNA generation, and preamplification as well as library preparation were performed by using the Chromium Single Cell 3′ v2 reagent kit according to the instruction manual. The isolation of single CT cells after FACS enrichment and the generation of cDNA for Prrx1 converted samples (uninjured samples and blastemas at 3 to 18 dpa) and limb bud experiments were performed by using the Fluidigm C1 system. Single cells were captured on a medium-sized (10- to 17-μm cell diameter) or large-sized (17- to 25-μm cell diameter) integrated fluidic circuit RNA-seq chip (Fluidigm). Cells were loaded on the chip at a concentration of 200 to 300 cells/μl and imaged by phase-contrast microscopy to assess the number of cells per capture site. Cell capture, cell lysis, reverse transcription, and cDNA amplification were performed on the chip. Reverse transcription of mRNA was performed by using the SMARTer Ultra Low RNA kit v2 for Fluidigm C1 (TaKaRa Clontech), where ERCC (External RNA Controls Consortium) RNA spike-in mix (Ambion) was added to the lysis reaction and processed in parallel to cellular mRNA. The Adventage2 polymerase chain reaction (PCR) kit (TaKaRa Clontech) was used for cDNA preamplification. For Col1a2 converted samples (uninjured samples and blastema at 11 dpa), single cells were isolated and cDNA was generated by performing SmartSeq2 (37). Briefly, individual cells were sorted in 4 μl of lysis mix [0.1% Triton X-100, 2.5 mM (each) deoxynucleoside triphosphate, and 2.5 μM dT30 reverse transcription primer] in single wells of a 96-well plate and snap frozen on dry ice. Generation of cDNA and its amplification exactly followed SmartSeq2 instructions, with application of 19× PCR cycles for blastema cells and 22× for uninjured mature cells. Similarly, stage 28 limb field samples were processed; however, 1 μl containing a single cell was pipetted into 3 μl of lysis mix, and snap freezing was not applied. For preamplifying cDNA, 18× PCR cycles were applied.

RNA-seq library construction and cDNA sequencing

A high-throughput electrophoresis-based fragment analyzer (Fragment Analyzer, Advanced Analytical Technologies) was used to assess the single-cell cDNA fragment size distribution and its concentration of every single cell for C1 and SmartSeq2 experiments. Illumina libraries were constructed by using the Illumina Nextera XT DNA sample preparation kit according to the protocol supplied by Fluidigm. Up to 192 libraries were pooled (3 μl each) and purified with SPRI beads. Library concentration and size distribution were assessed on an Agilent Bioanalyzer and with Qubit double-stranded DNA high-sensitivity assay kits and a Qubit 2.0 fluorometer. Each cell was paired-end sequenced (100 base reads) on an Illumina HiSeq 2500 to a depth of 2 to 5 million reads, and base calling, adaptor trimming, and de-multiplexing were performed as described previously (38, 39). 10× Genomics sample libraries were sequenced on individual lanes on an Illumina HiSeq 2500, and base calling, adaptor trimming, and de-multiplexing of single cells were performed by using 10× Genomics Cell Ranger 2.1 software.

Processing and analysis of scRNA-seq data

Raw reads for C1 and SmartSeq2 experiments were processed by using a custom script and pseudoaligned to the axolotl transcriptome (40) by using Kallisto with default settings (41). Transcript levels were quantified based on the pseudoalignment as transcripts per million reads (TPM) generated by Kallisto. Single-cell reads obtained with 10× Genomics were aligned to the axolotl transcriptome by using STAR (42) implemented in the 10× Genomics Cell Ranger 2.0 software, which generated absolute unique molecular identifier counts. Expression values (TPM/counts) of different isoforms of the same gene were summed afterward by using custom Perl scripts. We excluded cells for which fewer than 1000 genes were detected, and ribosomal protein genes were excluded from all datasets. For limb bud samples, cells were excluded if they did not express Prrx1 (TPM > 0). Expression levels were converted to the log space by taking the log2(TPM). R studio (43) was used to run custom R scripts for performing principal components analysis (PCA) (FactoMineR package) and hierarchical clustering (stats package) and for constructing heat maps, correlation plots, scatter plots, violin plots, dendrograms, bar graphs, or histograms. Generally, ggplot2 and gplots packages were used to generate data graphs. The Seurat package (44) implemented in R was used to identify cell clusters and perform differential gene expression analysis on the basis of tSNE. Diffusion maps (17) (destiny/dpt packages) were used to analyze cell lineage relationships. Covariance network analysis and visualizations were done by using igraph implemented in R (45). Gene ontology (GO) enrichment analyses were performed by using DAVID Bioinformatics Resources 6.7 (46). A list of transcription factors was obtained from the online animal transcription factor database AnimalTFDB (47).

Supplementary Materials

www.sciencemag.org/content/362/6413/eaaq0681/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S11

Tables S1 to S10

Reference (49)

References and Notes

Acknowledgments: We thank B. Borgonovo and D. Drechsel for technical advice; P. Keller for help with anti-PRRX1 antibodies; M. Blaschek, T. Torrecilla Lobos, H. Andreas, B. Gruhl, S. Mögel, and A. Wagner for excellent animal care and breeding; M. Corriero for the axolotl sketch; T. Engelmaier for help with hematoxylin and eosin staining; Z. He, G. Parra, and J. Kelso for computational support; A. Weihmann and B. Schellbach for sequencing; and members of the Tanaka lab and the Treutlein lab and S. Pääbo for helpful discussions. We thank the Core Unit Durchflusszytometrie (CUDZ) at the College of Veterinary Medicine, University of Leipzig, for providing access to FACS. Animal experiments were performed as approved by the State Authorities Saxony and the Magistrate of Vienna. Funding: We acknowledge financial support from the Max Planck Society (B.T.), DFG grant TA 274/3-3 (P.M., E.M.T.), DFG TA 274/11-1 (D.K.), DFG/ANR TA274/13-1 (P.M., E.M.T.), an ERC Advanced Investigator award (P.M., S.H., and E.M.T.), an ERC starting grant (B.T.), and FZ111 and institutional support from the IMP (E.M.T.). Author contributions: T.G., P.M., D.K., J.G.C., E.M.T., and B.T. designed the study and wrote the manuscript. P.M. generated transgenic axolotl lines with assistance from M.S. and S.K. P.M. performed lineage tracing and immunohistochemistry with help from S.H. P.M. and D.K. performed grafting experiments. P.M. and W.M. performed Brainbow clonal tracing. T.G. performed single-cell experiments with assistance from J.G.C., D.K., and P.M. M.G.-S. performed FACS. T.G. analyzed scRNA-seq data with help from B.T., J.G.C., J.K., P.M., D.K., and E.M.T. S.N. provided assistance with the axolotl transcriptome. J.D.C. generated the HOXA11 antibody. Competing interests: The authors declare no competing interests. Data and materials availability: Single-cell count matrices and fastq files are deposited in NCBI’s Gene Expression Omnibus under the accession number GSE106269.
View Abstract

Navigate This Article