Zhong J et al. (2023), Multi-species atlas resolves an axolotl limb de...

XB-ART-60861

Nat Commun 2023 Oct 10;141:6346. doi: 10.1038/s41467-023-41944-w.

Show Gene links Show Anatomy links

Multi-species atlas resolves an axolotl limb development and regeneration paradox.

Zhong J , Aires R , Tsissios G , Skoufa E , Brandt K , Sandoval-Guzmán T , Aztekin C .

???displayArticle.abstract???
Humans and other tetrapods are considered to require apical-ectodermal-ridge (AER) cells for limb development, and AER-like cells are suggested to be re-formed to initiate limb regeneration. Paradoxically, the presence of AER in the axolotl, a primary model organism for regeneration, remains controversial. Here, by leveraging a single-cell transcriptomics-based multi-species atlas, composed of axolotl, human, mouse, chicken, and frog cells, we first establish that axolotls contain cells with AER characteristics. Further analyses and spatial transcriptomics reveal that axolotl limbs do not fully re-form AER cells during regeneration. Moreover, the axolotl mesoderm displays part of the AER machinery, revealing a program for limb (re)growth. These results clarify the debate about the axolotl AER and the extent to which the limb developmental program is recapitulated during regeneration.

???displayArticle.pubmedLink??? 37816738
???displayArticle.pmcLink??? PMC10564727
???displayArticle.link??? Nat Commun

Species referenced: Xenopus laevis
Genes referenced: arid3b axin2 bambi bmp2 bmp4 cdh1 dlx2 dlx3 dlx5 dpt dsp epcam fgf7 fgf8 fgf9 fn1 hpse jag2 kazald1 krt17 lef1 lgr5 lum msx2 notch1 prdm1 sp6 sp8 tcf7 tfap2c tp63 vwa2 wnt3a wnt5a

???attribute.lit??? ???displayArticles.show???

	Fig. 1 Multi-species limb atlas reveals developing axolotl limbs have cells with apical-ectodermal-ridge (AER) characteristics. a Schematic to generate a multi-species limb atlas using publicly available single-cell RNA-Seq (scRNA-Seq) datasets. b UMAP plot of Seurat-integrated multi-species limb atlas. Individual datasets from each species with different developmental stages are integrated. Dots are colored by cell identities. c UMAP plot of the ectodermal lineage. Dots are colored by cell identities. d UMAP plot of species contribution to the AER cluster. Dots are colored by species. e Dot plot showing AER marker expressions in AER cells from different species. The dot color indicates the mean expression that was normalized to the max of each dataset and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. f Heatmap showing the MetaNeighbor score for pair-wise similarities of basal ectoderm and AER clusters. X- and Y-axes indicate species and developmental stages. Asterisks (*) denote the pairs with scores above 0.9. Source data provided as a Source data file. g Heatmap showing signaling ligands gene set enrichment analysis scores for AER, and non-AER basal ectoderm clusters. The basal ectoderm represents the transcriptome-wide most similar population to the AER, and is used for comparison. Colored dots in the Y-axis indicate different species. Source data provided as a Source data file. h Single optical section of z-stacks of confocal images of Stage 46 axolotl forelimb buds stained for Dr999-Pmt21178 (referred to as Dr999) mRNA via hybridization-chain-reaction (HCR). Different z-stacks were shown from left to right, representing different levels of the dorsal-ventral axis. (Top) Gray, Hoechst; Bottom Gray, Dr999 mRNA. Scale bar: 100 μm. i Max-projection confocal image of Stage 53 axolotl forelimb digit tips stained for Dr999 mRNA via HCR. Green, Dr999 mRNA; Gray, Hoechst. Scale bar: 250 μm. j Zoomed single optical section image of the axolotl limb bud from (h) stained for Dr999 mRNA. Red arrows show Dr999+ squamous cells, and yellow arrows show outer layer peridermal cells. The basement membrane is labeled with a dashed line. Green, Dr999 mRNA; Gray, Hoechst. Scale bar: 10 μm.
	Fig. 2 Spatial organization of axolotl AER cells during limb development. a Max-projection confocal images of ST 46 (top), 50 (mid), and 53 (bottom) axolotl forelimbs correspond to ST50, 52, 54 hindlimb scRNA-Seq data, respectively, stained for Dr999 and Vwa2 mRNA via HCR. Green, Dr999 mRNA; Red, Vwa2 mRNA; Blue, Hoechst. Yellow * is added to the bottom images to indicate autofluorescence. Scale bar: 100 μm (top and mid) and 200 μm (bottom). b Max-projection confocal image of axolotl ST 53 hindlimb buds stained for Dr999 and Epcam mRNA via HCR. Gray, Hoechst; Magenta, Dr999 mRNA; Cyan, Epcam mRNA. Enlarged views of corresponding regions indicated by the numbered yellow boxes are shown side by side. Scale bar: 20 μm for zoomed images, 100 μm for the merged image and 200 μm for the others. c Msx2 expression in axolotl limb buds based on the scRNA-Seq datasets. d (Left) Max-projection confocal image of axolotl Stage 53 hindlimb buds stained for Dr999 and Msx2 mRNA via HCR. Scale bar: 100 μm. Red arrows show skin cells, and the green arrows show non-skin mesodermal cells to highlight Msx2 staining profile. (Right) Close-up images of the dashed labeled region of Dr999 and Msx2 double-positive skin cells. Scale bar: 10 μm.
	Fig. 3 Axolotl AER cells are not entirely re-formed during limb regeneration. a Schematics describing examples of the used scRNA-Seq datasets of axolotl limb development and regeneration are illustrated. b UMAP plot of the basal ectoderm cells of the integrated axolotl limb development and regeneration datasets. Subclustered cell identities are labeled by different colors and text. c Sample contribution to the integrated AER/AEC cluster from b. Red dots indicate cells from the selected sample; gray dots indicate the other cells in the AER/AEC cluster. hpa: hours-post amputation; dpa: days-post amputation. d Dot plot showing AER marker expressions in (left) the scRNA-seq datasets of axolotl limb development and regeneration AER or AEC clusters, respectively, and (right) spatial transcriptomics (Visium) AEC cluster (g). The dot color indicates the mean expression that was normalized to the max of each gene; the dot size represents the percentage of cells with non-zero expressions. Please note that the Visium and scRNA-seq datasets were normalized separately within the dataset. Significant differentially expressed genes between AER and AEC were labeled in bold (two-sided Wilcoxon rank-sum test; P-values < 0.05). e Dot plot showing signaling ligand expressions in (left) scRNA-seq datasets of axolotl limb development and regeneration AER or AEC cells, respectively, and (right) Visium AEC cluster (g). The dot color indicates the mean expression that was normalized to the max of each gene; the dot size represents the percentage of cells with non-zero expressions. Please note that the Visium and scRNA-seq datasets were normalized separately within the dataset. Significant differentially expressed genes between AER and AEC were labeled in bold (two-sided Wilcoxon rank-sum test; P-values < 0.05). f Hematoxylin and eosin stained 5 dpa axolotl limb regeneration tissue section. The tissue is oriented with the anterior to the top, the posterior to the bottom, proximal to the left, and distal to the right. g Clustering of the Visium spots identified known tissue types, including the AEC (pink) and the blastema (green). For the full clustering results, see Supplementary Fig. 15. h Expression profiles of selected markers in the AEC and the blastema clusters.
	Fig. 4 Axolotl limbs contain mesodermal cells showing part of the epithelial AER transcriptional program. a Schematics of the strategy to evaluate the mesodermal cells showing the AER transcriptional program. AEC and connective tissue (CT) cells are clustered based on differentially expressed genes (DEGs) in the AER cluster of the multi-species limb atlas. b UMAP plots of the clustering based on the DEGs in the AER cluster of the multi-species limb atlas (Fig. 1b) for development datasets. Mouse E10.5, chicken E4.5, human CS13, frog NF51, and axolotl ST52 UMAPs are shown. Light blue, limb bud mesoderm cells; pink, AER cells. c UMAP plots of the clustering based on the DEGs in the AER cluster of the multi-species limb atlas for the regeneration datasets. Dark blue, CT cells; dark pink, AEC cells. CT cells gathered with the AEC and the AEC are circled with a dashed line. d Heatmap showing shared gene expressions in the axolotl CT and AEC populations. Dendrogram based on gene expression profile indicates a subset of CT cells shows similarity to the AEC population and is highlighted with black dendrograms. Note that these shared genes were identified by K-means clustering (Supplementary Fig. 22a) of the whole AER gene list used in Fig. 4b, c, except for the manually added Fgf8 due to its high relevance to AER. e The expression profile of the example epithelial and fibroblast genes in the UMAP plot in Fig. 4c axolotl regeneration dataset. f Expression profiles of the example epithelial and fibroblast genes are visualized on the Visium dataset. Please note that the Krt17 data is the same as in Fig. 3h. g Violin plot showing expression levels of the example epithelial and fibroblast genes in the CT, the blastema, and the AEC clusters of the Visium dataset. h Schematics illustrating the AER transcriptional program in developing (axolotl, frog, human, mouse, chicken) (left) and regenerating (axolotl, frog tadpole) limbs (right). Red to orange colors are used to indicate changes in the AER transcriptional program, relevant explicitly in the axolotl regeneration schematics. Dash rectangle indicates the zoom-in view of the axolotl regenerating limb with selected AER genes.
	Supplementary Fig. 1 \| Quality assessment of the re-analyzed publicly available limb development datasets of five species. a) Barplot showing the used cell numbers after filtering in each dataset is visualized. Two human datasets of the same developmental stage from independent studies were used and denoted as Human_CS13 (1) and Human CS13 (2). Please see Supplementary Data 1 for full details. b) Boxplot showing the transcript number per cell in each dataset is visualized. Please see Supplementary Data 1 for full details. c) Boxplot showing the gene numbers per cell in each dataset is visualized. Please see Supplementary Data 1 for full details. d) Boxplot showing the percentage of mitochondrial genes per cell in each dataset is visualized. Please see Supplementary Data 1 for full details.
	Supplementary Fig. 2 \| UMAP visualization of the re-analyzed individual datasets of limb development. UMAP visualization of individual limb development datasets for (a) human, (b) chicken, (c) mouse, (d) frog, and (e) axolotl are visualized. Two human datasets of the same developmental stage from independent studies were used and denoted as Human_CS13 (1) and Human CS13(2). Lineage and cell type annotations are based on marker genes in Supplementary Fig. 3, and labeled by different colors and text, except the axolotl AER cells, which are identified based on Fig. 1d and re-plotted in this figure. Note that to highlight the identified AER-like cells in the axolotl stage 50 sample, a zoom- in view of the ectoderm population is shown in the upper left corner.
	Supplementary Fig. 3 \| Seurat-integrated multi-species limb atlas. a) UMAP plot of Seurat-integrated multi-species limb atlas is shown. Individual datasets from each species and developmental stage (full list in Supplementary Fig. 2 and Supplementary Data 1) are integrated. Cell lineage and cell type identities are labeled by different colors and text. Please note that this plot is the same as in Fig. 1b. b) UMAP plots of species and developmental stage contribution to the Seurat-integrated multi- species limb atlas are shown for each species. Species are color-coded and developmental stages are indicated with a shading of the same color. c) Violin plots of the number of detected genes (top) and the percentage of expressed mitochondrial genes (bottom) in each cell type in the Seurat-integrated multi-species limb atlas. Cell types were defined in Fig. 1a. d) Dotplot showing used marker genes to annotate clusters. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expressions.
	Supplementary Fig. 3 \| Seurat-integrated multi-species limb atlas. e) A confusion matrix is plotted to determine the annotation accuracy for cells between the Seurat- integrated multi-species atlas and individual maps from Supplementary Fig. 2. f) Pie chart showing the percentage of cells from each species and developmental stages in the Goblet cells cluster in the Seurat-integrated UMAP in Fig. 1a. g) Barplot showing the ratio of identified AER cluster to the whole basal ectoderm cluster in indicated species and developmental stages. h) UMAP plot of Seurat-integrated multi-species limb atlas without the chicken dataset. Cells are colored by their lineages and cell type identities. Inserted box: the species contribution to the AER cluster is visualized. Cells from different species are color-coded.
	Supplementary Fig. 4 \| SAMap-integrated multi-species limb atlas a) UMAP plot of SAMap -integrated multi-species limb atlas is shown. Schwann Endo Immune Percent 100 75 50 25 0 Ectoderm Limb bud mesoderm Limb bud mesoderm Coarse annotation Individual clustering 17.5 0.2 2.4 0.2 2.9 0.3 8.5 1.1 9.5 0.3 0.5 0 0.2 3.6 27.1 44.9 0 5.9 0 0.2 0 88.1 7.3 0 0 1.9 0 0.7 89.5 0.1 0 0 5.5 1 0.4 89.8 0 0.1 0 0 0 0 90 0 0.1 0 0 0.1 0 0 0.4 99.6 0 0 0 0 0 0 0 0 99 0.9 0 0 0.1 0 0 0 0 0 SAMap integration Limb bud mesoderm Ectoderm Muscle Blood Immune Endo Pericyte Schwann Immune_Blood Unknown Individual datasets from each species and developmental stage (full list in Supplementary Fig. 2 and Supplementary Data 1) are integrated. Clustering and annotation based on marker gene expressions (please see Supplementary Fig. 3) are indicated lineage and cell type identities, and labeled by different colors and text. b) UMAP plots of species and developmental stage contribution to the SAMap-integrated multi- species limb atlas are shown for each species. Species are color-coded and developmental stages are indicated with a shading of the same color.
	Supplementary Fig. 4 \| SAMap-integrated multi-species limb atlas. c) UMAP plot of species contribution to the SAMap-integrated multi-species limb atlas AER cluster is visualized. Cells from different species are color-coded and the developmental stage is indicated. d) A confusion matrix is plotted to determine the annotation accuracy between the SAMap- integrated multi-species atlas and individual maps from Supplementary Fig. 2.
	Supplementary Fig. 5 \| AER marker expressions in AER cells of the analyzed species. Dotplot showing an extended list of AER markers. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. Please note that some of the genes in this figure and their values are the same as in Fig. 1e.
	Supplementary Fig. 6 \| Basal ectoderm and AER clusters similarities of the analyzed species. Heatmaps showing the pair-wise similarity scores of basal ectoderm and AER clusters across species. X and Y axes indicate different species with their developmental stages. Similarity scores are calculated by using (a) the top 3000 highly variable genes, and (b) the top 50 principal components. Asterisks (*) denote the pairs with similarity scores above 0.9. AER clusters are compared to basal ectoderm clusters to highlight specific transcriptome profiles.
	Supplementary Fig. 7 \| Signaling ligand gene set enrichment scores for lineages and cell types in the analyzed species. Heatmap showing the (a) calculated and (b) max-normalized signaling ligand gene set enrichment scores for each annotated population in the multi-species limb atlas (Fig. 1b). Grey indicates missing cell cluster. Please note that the same data for the basal ectoderm and the AER clusters are shown in Fig. 1g.
	Supplementary Fig. 8 \| Signaling pathways-associated gene expressions in AER cells in the analyzed species. Dotplots showing expression profiles of various genes from FGF (a), WNT (b), BMP (c), NOTCH (d), and TGFb (e) pathways. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression.
	Supplementary Fig. 9 \| Spatial organization of axolotl AER cells during limb development. a) Heatmap showing the expression of 14 shared AER-specific differentially expressed genes in different cell types in axolotl developing limbs. Stages are indicated on the top. b) Violin plots showing expression levels of the putative axolotl AER markers are shown. Please note that the analyzed scRNA-Seq datasets are for hindlimbs, and morphologies of corresponding developing limbs are indicated in the figure. c) Co-expression of AER-like genes in the ectoderm in axolotl limb buds. Scatter plot (lower left) shows the co-expression pattern of Dr999-Pmt21178 and Vwa2. Each dot represents a cell whose location is defined by its expression of Dr999-Pmt21178 (x-axis) and Vwa2 (y-axis). AER-like cells are colored in shades of green while non-AER ectoderm in grey. Density plots representation of the distribution of the Dr999-Pmt21178 (top) and Vwa2 (lower right) expression in AER (in green) and non-AER ectoderm (in grey). d) (Left) Merged image of the top row in Fig. 2a. Scale bar 100 μm. (Right) Enlarged views of corresponding regions indicated by the numbered yellow boxes. Scale bar 20 μm.
	Supplementary Fig. 10 \| Vwa2 expression in humans, mice, chickens, and frogs. a) Violin plots showing Vwa2 expression across all cell types in indicated species and developmental stages based on the analyzed scRNA-Seq datasets. b) Max-projection confocal image of mouse E10.5 hindlimb limb buds stained for Fgf8, and Vwa2 mRNA via HCR. Gray, Hoechst; Magenta, Fgf8 mRNA; Cyan, Vwa2 mRNA. Scale bar: 100 μm for both. c) Max-projection confocal image of chicken HH22 hindlimb buds stained for Fgf8, and Vwa2 mRNA via HCR. Gray, Hoechst; Magenta, Fgf8 mRNA; Cyan, Vwa2 mRNA. Scale bars: (Top) 250 μm and (Bottom) 25 μm.
	Supplementary Fig. 10 \| Vwa2 expression in humans, mice, chickens, and frogs. d) Merged image of (b). Zoom-in views of yellow boxes regions are shown by side. Scale bar: 100 μm and 10 μm for enlarged images. e) Merged image of (c). Zoom-in views of yellow boxes regions are shown by side. Scale bar: 100 μm for full image and 10 μm for enlarged images.
	Supplementary Fig. 11 \| AER-specific Fgf expression in developing axolotl limbs. a) Max-projection confocal image of axolotl Stage 53 hindlimb buds stained for Fgf7 and Fgf18 mRNA via HCR. Gray, Hoechst; Cyan, Fgf7; Magenta, Fgf18 mRNA; Scale bars: (Top) 100 μm and (Bottom) 20 μm. The red arrow indicates the target signals. b) Zoomed region of (a) indicated by yellow box number 1 shows the target signals. Scale bars: 20 μm. c) Zoomed region of (a) indicated by yellow box number 2 shows high background signals. Scale bars: 20 μm.
	Supplementary Fig. 12 \| Quality assessment of the re-analyzed publicly available Xenopus limb regeneration dataset. a) Barplots showing the used cell numbers, the transcript number per cell, the gene numbers per cell, and the percentage of mitochondrial genes per cell after filtering in Xenopus laevis 5 days-post amputation (dpa) limb regeneration dataset are visualized.
	Supplementary Fig. 12 \| Quality assessment of the re-analyzed publicly available Xenopus limb regeneration dataset. b) UMAP plot of the identified and annotated clusters. Clusters were annotated based on marker genes listed in Supplementary Fig. 12c. c) Heatmap showing the expression profile of marker genes that are used to annotate clusters. d) Dotplot showing signaling ligand expressions in the AER or AEC clusters during Xenopus laevis limb development and regeneration, respectively. Please note that a similar figure was generated in Aztekin et al, 2021, but in this figure, NF Stage 50, 51, and 52 from Lin et al. 2021 were used.
	Supplementary Fig. 13 \| Quality assessment of re-analyzed publicly available axolotl limb regeneration dataset. a) Barplots showing the used cell numbers, the transcript number per cell, the gene numbers per cell, and the percentage of mitochondrial genes per cell after filtering in axolotl limb regeneration datasets are visualized. b) UMAP plot of the identified and annotated clusters. Clusters were annotated based on marker genes listed in Supplementary Fig. 13d. The AEC cells were identified by the analysis described in the text and in Fig. 3a-c, then re-plotted and colored dark pink in this figure. c) Sample contribution to the axolotl regeneration dataset is visualized. Red dots indicate cells from the selected sample; gray dots indicate the other cells. hpa: hours-post amputation; dpa: days-post amputation. d) Heatmap showing the expression profile of marker genes that are used to annotate clusters, except AEC cells and their marker expressions indicated following the analysis indicated in text and Fig. 3a-c.
	Supplementary Fig. 14 \| Differentially expressed gene analysis between axolotl AER and AEC cells. a) Volcano plot showing the differentially expressed genes (DEGs) between axolotl AER and AEC populations. Red and blue dots indicate genes significantly enriched in the AER or the AEC, respectively. Gray dots indicate statistically not significant genes. The top 10 DEGs are labeled. b) Barplots showing enriched GO terms based on the top 200 DEGs (ordered by fold change) in AER (red) and AEC (blue) cells.
	Supplementary Fig. 15 \| Quality assessment and clustering results of axolotl 5 days-post- amputation 10X Visium spatial transcriptomics (Visium) dataset. a) Barplots showing used Visium spot numbers, transcript number per spot, gene numbers per spot, and the percentage of mitochondrial genes per spot after filtering are visualized. b) 5 dpa axolotl limb regeneration tissue section that is used for the Visium is shown and stained for hematoxylin and eosin. Please note that this image is the same as Fig. 3f but added in this figure for easier visual comparison to Supplementary Fig. 15d, which contains full annotation. c) UMAP plot of the Visium clusters is shown. Clusters were annotated by marker gene expressions (Supplementary Fig. 15e) and confirmed by morphologies in hematoxylin and eosin staining on the tissue section in Fig. 3f and Supplementary Fig. 15b. The identified clusters are color-coded and tissue annotation is labeled. d) Visium spots were colored by tissue types identified in Supplementary Fig. 15c.
	Supplementary Fig. 15 \| Quality assessment and clustering results of axolotl 5 days-post- amputation 10X Visium spatial transcriptomics (Visium) dataset. e) Heatmap showing the top 20 DEGs (ordered by fold change) for the identified clusters is shown.
	Supplementary Fig. 16 \| Consensus nonnegative matrix factorization (cNMF) to identify AER- related modules. a) The strategy to perform consensus negative matrix factorization (cNMF) is described. cNMF was used on AER and CT clusters to identify cell identity and cell activity modules. b) Stability (blue line) is measured by the Euclidean distance silhouette score of the clustering and Frobenius error of the consensus solution (red line) for each tested k value (module number) shown in the X axis. Please see Methods for more detail. c) Violin plots showing activity scores for identified cNMF modules in CT and AER clusters. Modules specific to AER, CT, or shared are color-coded: pink, AER-specific modules; blue, CT- specific modules; green, shared modules.
	Supplementary Fig. 16 \| Consensus nonnegative matrix factorization (cNMF) to identify AER- related modules. d) Enriched GO terms are determined for the top 100 genes (ordered by gene expression program scores) in the identified AER-specific modules. Barplot showing the top GO is visualized, and GO terms are labeled.
	Supplementary Fig. 17 \| Single cell gene set enrichment analysis (scGSEA) using differentially expressed AER gene sets for individual datasets. UMAP visualization of the scGSEA of indicated datasets of limb development (a) and regeneration (b) based on the top 500 AER DEGs in multi-species limb atlas (Fig. 1b). The enrichment scores are indicated in shades of red. Cells failed to pass the enrichment threshold are colored in grey.
	Supplementary Fig. 18 \| Clustering based on differentially expressed AER gene sets for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 AER DEGs in multi-species limb atlas (Fig. 1b) for individual datasets. In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line. Please note that Mouse E10.5, chicken E4.5, human CS13, frog ST51, and axolotl ST52 development datasets, and the axolotl and frog tadpole regeneration datasets are also shown in Fig 4b and c.
	Supplementary Fig. 19 \| Clustering based on cNMF module 11 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 11 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 20 \| Clustering based on cNMF module 12 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 12 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 21 \| Clustering based on cNMF module 13 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 13 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 22 \| AER gene expression profile in AER, AEC, limb bud mesoderm, and connective tissue cells during limb development and regeneration. Heatmaps showing the expression profile of 500 AER DEGs in regenerating limbs (a) and developing limbs (b). For axolotl regenerating limbs (Axolotl Reg), genes were grouped into three by the K- means algorithm, and the first group is shown in Fig. 4d. Please note that part of the axolotl regeneration dataset is also presented in Fig. 4f. Colored squares above the heatmap represent cells from indicated populations: light blue, limb bud mesoderm cells in the development dataset; light pink, AER cells in the development dataset; dark blue, connective tissue cells in the regeneration dataset; dark pink, AEC cells in the regeneration dataset.
	Supplementary Fig. 23 \| Mesodermal cells showing AER transcriptional program are present in intact limbs and express signaling ligands. a) (Top) UMAP plot from Fig. 4c and Supplementary Fig. 18b axolotl regeneration dataset is used to distinguish mesodermal cells showing part of the AER program and labeled by green. The rest of the connective tissue is labeled dark blue, and AEC is labeled dark pink. (Bottom) Dotplot showing signaling ligands expressions in the populations distinguished in Supplementary Fig. 23a. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. b) The expression profile of previously reported regeneration-associated genes (18–20) in the UMAP plot of clustering in Fig 4c axolotl regeneration dataset. c) Sample contribution to UMAP plot of re-clustered regeneration dataset is visualized. Red dots indicate cells from the selected sample; gray dots indicate all the other cells.

	Fig. 1 Multi-species limb atlas reveals developing axolotl limbs have cells with apical-ectodermal-ridge (AER) characteristics. a Schematic to generate a multi-species limb atlas using publicly available single-cell RNA-Seq (scRNA-Seq) datasets. b UMAP plot of Seurat-integrated multi-species limb atlas. Individual datasets from each species with different developmental stages are integrated. Dots are colored by cell identities. c UMAP plot of the ectodermal lineage. Dots are colored by cell identities. d UMAP plot of species contribution to the AER cluster. Dots are colored by species. e Dot plot showing AER marker expressions in AER cells from different species. The dot color indicates the mean expression that was normalized to the max of each dataset and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. f Heatmap showing the MetaNeighbor score for pair-wise similarities of basal ectoderm and AER clusters. X- and Y-axes indicate species and developmental stages. Asterisks (*) denote the pairs with scores above 0.9. Source data provided as a Source data file. g Heatmap showing signaling ligands gene set enrichment analysis scores for AER, and non-AER basal ectoderm clusters. The basal ectoderm represents the transcriptome-wide most similar population to the AER, and is used for comparison. Colored dots in the Y-axis indicate different species. Source data provided as a Source data file. h Single optical section of z-stacks of confocal images of Stage 46 axolotl forelimb buds stained for Dr999-Pmt21178 (referred to as Dr999) mRNA via hybridization-chain-reaction (HCR). Different z-stacks were shown from left to right, representing different levels of the dorsal-ventral axis. (Top) Gray, Hoechst; Bottom Gray, Dr999 mRNA. Scale bar: 100 μm. i Max-projection confocal image of Stage 53 axolotl forelimb digit tips stained for Dr999 mRNA via HCR. Green, Dr999 mRNA; Gray, Hoechst. Scale bar: 250 μm. j Zoomed single optical section image of the axolotl limb bud from (h) stained for Dr999 mRNA. Red arrows show Dr999+ squamous cells, and yellow arrows show outer layer peridermal cells. The basement membrane is labeled with a dashed line. Green, Dr999 mRNA; Gray, Hoechst. Scale bar: 10 μm.
	Fig. 2 Spatial organization of axolotl AER cells during limb development. a Max-projection confocal images of ST 46 (top), 50 (mid), and 53 (bottom) axolotl forelimbs correspond to ST50, 52, 54 hindlimb scRNA-Seq data, respectively, stained for Dr999 and Vwa2 mRNA via HCR. Green, Dr999 mRNA; Red, Vwa2 mRNA; Blue, Hoechst. Yellow * is added to the bottom images to indicate autofluorescence. Scale bar: 100 μm (top and mid) and 200 μm (bottom). b Max-projection confocal image of axolotl ST 53 hindlimb buds stained for Dr999 and Epcam mRNA via HCR. Gray, Hoechst; Magenta, Dr999 mRNA; Cyan, Epcam mRNA. Enlarged views of corresponding regions indicated by the numbered yellow boxes are shown side by side. Scale bar: 20 μm for zoomed images, 100 μm for the merged image and 200 μm for the others. c Msx2 expression in axolotl limb buds based on the scRNA-Seq datasets. d (Left) Max-projection confocal image of axolotl Stage 53 hindlimb buds stained for Dr999 and Msx2 mRNA via HCR. Scale bar: 100 μm. Red arrows show skin cells, and the green arrows show non-skin mesodermal cells to highlight Msx2 staining profile. (Right) Close-up images of the dashed labeled region of Dr999 and Msx2 double-positive skin cells. Scale bar: 10 μm.
	Fig. 3 Axolotl AER cells are not entirely re-formed during limb regeneration. a Schematics describing examples of the used scRNA-Seq datasets of axolotl limb development and regeneration are illustrated. b UMAP plot of the basal ectoderm cells of the integrated axolotl limb development and regeneration datasets. Subclustered cell identities are labeled by different colors and text. c Sample contribution to the integrated AER/AEC cluster from b. Red dots indicate cells from the selected sample; gray dots indicate the other cells in the AER/AEC cluster. hpa: hours-post amputation; dpa: days-post amputation. d Dot plot showing AER marker expressions in (left) the scRNA-seq datasets of axolotl limb development and regeneration AER or AEC clusters, respectively, and (right) spatial transcriptomics (Visium) AEC cluster (g). The dot color indicates the mean expression that was normalized to the max of each gene; the dot size represents the percentage of cells with non-zero expressions. Please note that the Visium and scRNA-seq datasets were normalized separately within the dataset. Significant differentially expressed genes between AER and AEC were labeled in bold (two-sided Wilcoxon rank-sum test; P-values < 0.05). e Dot plot showing signaling ligand expressions in (left) scRNA-seq datasets of axolotl limb development and regeneration AER or AEC cells, respectively, and (right) Visium AEC cluster (g). The dot color indicates the mean expression that was normalized to the max of each gene; the dot size represents the percentage of cells with non-zero expressions. Please note that the Visium and scRNA-seq datasets were normalized separately within the dataset. Significant differentially expressed genes between AER and AEC were labeled in bold (two-sided Wilcoxon rank-sum test; P-values < 0.05). f Hematoxylin and eosin stained 5 dpa axolotl limb regeneration tissue section. The tissue is oriented with the anterior to the top, the posterior to the bottom, proximal to the left, and distal to the right. g Clustering of the Visium spots identified known tissue types, including the AEC (pink) and the blastema (green). For the full clustering results, see Supplementary Fig. 15. h Expression profiles of selected markers in the AEC and the blastema clusters.
	Fig. 4 Axolotl limbs contain mesodermal cells showing part of the epithelial AER transcriptional program. a Schematics of the strategy to evaluate the mesodermal cells showing the AER transcriptional program. AEC and connective tissue (CT) cells are clustered based on differentially expressed genes (DEGs) in the AER cluster of the multi-species limb atlas. b UMAP plots of the clustering based on the DEGs in the AER cluster of the multi-species limb atlas (Fig. 1b) for development datasets. Mouse E10.5, chicken E4.5, human CS13, frog NF51, and axolotl ST52 UMAPs are shown. Light blue, limb bud mesoderm cells; pink, AER cells. c UMAP plots of the clustering based on the DEGs in the AER cluster of the multi-species limb atlas for the regeneration datasets. Dark blue, CT cells; dark pink, AEC cells. CT cells gathered with the AEC and the AEC are circled with a dashed line. d Heatmap showing shared gene expressions in the axolotl CT and AEC populations. Dendrogram based on gene expression profile indicates a subset of CT cells shows similarity to the AEC population and is highlighted with black dendrograms. Note that these shared genes were identified by K-means clustering (Supplementary Fig. 22a) of the whole AER gene list used in Fig. 4b, c, except for the manually added Fgf8 due to its high relevance to AER. e The expression profile of the example epithelial and fibroblast genes in the UMAP plot in Fig. 4c axolotl regeneration dataset. f Expression profiles of the example epithelial and fibroblast genes are visualized on the Visium dataset. Please note that the Krt17 data is the same as in Fig. 3h. g Violin plot showing expression levels of the example epithelial and fibroblast genes in the CT, the blastema, and the AEC clusters of the Visium dataset. h Schematics illustrating the AER transcriptional program in developing (axolotl, frog, human, mouse, chicken) (left) and regenerating (axolotl, frog tadpole) limbs (right). Red to orange colors are used to indicate changes in the AER transcriptional program, relevant explicitly in the axolotl regeneration schematics. Dash rectangle indicates the zoom-in view of the axolotl regenerating limb with selected AER genes.
	Supplementary Fig. 1 \| Quality assessment of the re-analyzed publicly available limb development datasets of five species. a) Barplot showing the used cell numbers after filtering in each dataset is visualized. Two human datasets of the same developmental stage from independent studies were used and denoted as Human_CS13 (1) and Human CS13 (2). Please see Supplementary Data 1 for full details. b) Boxplot showing the transcript number per cell in each dataset is visualized. Please see Supplementary Data 1 for full details. c) Boxplot showing the gene numbers per cell in each dataset is visualized. Please see Supplementary Data 1 for full details. d) Boxplot showing the percentage of mitochondrial genes per cell in each dataset is visualized. Please see Supplementary Data 1 for full details.
	Supplementary Fig. 2 \| UMAP visualization of the re-analyzed individual datasets of limb development. UMAP visualization of individual limb development datasets for (a) human, (b) chicken, (c) mouse, (d) frog, and (e) axolotl are visualized. Two human datasets of the same developmental stage from independent studies were used and denoted as Human_CS13 (1) and Human CS13(2). Lineage and cell type annotations are based on marker genes in Supplementary Fig. 3, and labeled by different colors and text, except the axolotl AER cells, which are identified based on Fig. 1d and re-plotted in this figure. Note that to highlight the identified AER-like cells in the axolotl stage 50 sample, a zoom- in view of the ectoderm population is shown in the upper left corner.
	Supplementary Fig. 3 \| Seurat-integrated multi-species limb atlas. a) UMAP plot of Seurat-integrated multi-species limb atlas is shown. Individual datasets from each species and developmental stage (full list in Supplementary Fig. 2 and Supplementary Data 1) are integrated. Cell lineage and cell type identities are labeled by different colors and text. Please note that this plot is the same as in Fig. 1b. b) UMAP plots of species and developmental stage contribution to the Seurat-integrated multi- species limb atlas are shown for each species. Species are color-coded and developmental stages are indicated with a shading of the same color. c) Violin plots of the number of detected genes (top) and the percentage of expressed mitochondrial genes (bottom) in each cell type in the Seurat-integrated multi-species limb atlas. Cell types were defined in Fig. 1a. d) Dotplot showing used marker genes to annotate clusters. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expressions.
	Supplementary Fig. 3 \| Seurat-integrated multi-species limb atlas. e) A confusion matrix is plotted to determine the annotation accuracy for cells between the Seurat- integrated multi-species atlas and individual maps from Supplementary Fig. 2. f) Pie chart showing the percentage of cells from each species and developmental stages in the Goblet cells cluster in the Seurat-integrated UMAP in Fig. 1a. g) Barplot showing the ratio of identified AER cluster to the whole basal ectoderm cluster in indicated species and developmental stages. h) UMAP plot of Seurat-integrated multi-species limb atlas without the chicken dataset. Cells are colored by their lineages and cell type identities. Inserted box: the species contribution to the AER cluster is visualized. Cells from different species are color-coded.
	Supplementary Fig. 4 \| SAMap-integrated multi-species limb atlas a) UMAP plot of SAMap -integrated multi-species limb atlas is shown. Schwann Endo Immune Percent 100 75 50 25 0 Ectoderm Limb bud mesoderm Limb bud mesoderm Coarse annotation Individual clustering 17.5 0.2 2.4 0.2 2.9 0.3 8.5 1.1 9.5 0.3 0.5 0 0.2 3.6 27.1 44.9 0 5.9 0 0.2 0 88.1 7.3 0 0 1.9 0 0.7 89.5 0.1 0 0 5.5 1 0.4 89.8 0 0.1 0 0 0 0 90 0 0.1 0 0 0.1 0 0 0.4 99.6 0 0 0 0 0 0 0 0 99 0.9 0 0 0.1 0 0 0 0 0 SAMap integration Limb bud mesoderm Ectoderm Muscle Blood Immune Endo Pericyte Schwann Immune_Blood Unknown Individual datasets from each species and developmental stage (full list in Supplementary Fig. 2 and Supplementary Data 1) are integrated. Clustering and annotation based on marker gene expressions (please see Supplementary Fig. 3) are indicated lineage and cell type identities, and labeled by different colors and text. b) UMAP plots of species and developmental stage contribution to the SAMap-integrated multi- species limb atlas are shown for each species. Species are color-coded and developmental stages are indicated with a shading of the same color.
	Supplementary Fig. 4 \| SAMap-integrated multi-species limb atlas. c) UMAP plot of species contribution to the SAMap-integrated multi-species limb atlas AER cluster is visualized. Cells from different species are color-coded and the developmental stage is indicated. d) A confusion matrix is plotted to determine the annotation accuracy between the SAMap- integrated multi-species atlas and individual maps from Supplementary Fig. 2.
	Supplementary Fig. 5 \| AER marker expressions in AER cells of the analyzed species. Dotplot showing an extended list of AER markers. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. Please note that some of the genes in this figure and their values are the same as in Fig. 1e.
	Supplementary Fig. 6 \| Basal ectoderm and AER clusters similarities of the analyzed species. Heatmaps showing the pair-wise similarity scores of basal ectoderm and AER clusters across species. X and Y axes indicate different species with their developmental stages. Similarity scores are calculated by using (a) the top 3000 highly variable genes, and (b) the top 50 principal components. Asterisks (*) denote the pairs with similarity scores above 0.9. AER clusters are compared to basal ectoderm clusters to highlight specific transcriptome profiles.
	Supplementary Fig. 7 \| Signaling ligand gene set enrichment scores for lineages and cell types in the analyzed species. Heatmap showing the (a) calculated and (b) max-normalized signaling ligand gene set enrichment scores for each annotated population in the multi-species limb atlas (Fig. 1b). Grey indicates missing cell cluster. Please note that the same data for the basal ectoderm and the AER clusters are shown in Fig. 1g.
	Supplementary Fig. 8 \| Signaling pathways-associated gene expressions in AER cells in the analyzed species. Dotplots showing expression profiles of various genes from FGF (a), WNT (b), BMP (c), NOTCH (d), and TGFb (e) pathways. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression.
	Supplementary Fig. 9 \| Spatial organization of axolotl AER cells during limb development. a) Heatmap showing the expression of 14 shared AER-specific differentially expressed genes in different cell types in axolotl developing limbs. Stages are indicated on the top. b) Violin plots showing expression levels of the putative axolotl AER markers are shown. Please note that the analyzed scRNA-Seq datasets are for hindlimbs, and morphologies of corresponding developing limbs are indicated in the figure. c) Co-expression of AER-like genes in the ectoderm in axolotl limb buds. Scatter plot (lower left) shows the co-expression pattern of Dr999-Pmt21178 and Vwa2. Each dot represents a cell whose location is defined by its expression of Dr999-Pmt21178 (x-axis) and Vwa2 (y-axis). AER-like cells are colored in shades of green while non-AER ectoderm in grey. Density plots representation of the distribution of the Dr999-Pmt21178 (top) and Vwa2 (lower right) expression in AER (in green) and non-AER ectoderm (in grey). d) (Left) Merged image of the top row in Fig. 2a. Scale bar 100 μm. (Right) Enlarged views of corresponding regions indicated by the numbered yellow boxes. Scale bar 20 μm.
	Supplementary Fig. 10 \| Vwa2 expression in humans, mice, chickens, and frogs. a) Violin plots showing Vwa2 expression across all cell types in indicated species and developmental stages based on the analyzed scRNA-Seq datasets. b) Max-projection confocal image of mouse E10.5 hindlimb limb buds stained for Fgf8, and Vwa2 mRNA via HCR. Gray, Hoechst; Magenta, Fgf8 mRNA; Cyan, Vwa2 mRNA. Scale bar: 100 μm for both. c) Max-projection confocal image of chicken HH22 hindlimb buds stained for Fgf8, and Vwa2 mRNA via HCR. Gray, Hoechst; Magenta, Fgf8 mRNA; Cyan, Vwa2 mRNA. Scale bars: (Top) 250 μm and (Bottom) 25 μm.
	Supplementary Fig. 10 \| Vwa2 expression in humans, mice, chickens, and frogs. d) Merged image of (b). Zoom-in views of yellow boxes regions are shown by side. Scale bar: 100 μm and 10 μm for enlarged images. e) Merged image of (c). Zoom-in views of yellow boxes regions are shown by side. Scale bar: 100 μm for full image and 10 μm for enlarged images.
	Supplementary Fig. 11 \| AER-specific Fgf expression in developing axolotl limbs. a) Max-projection confocal image of axolotl Stage 53 hindlimb buds stained for Fgf7 and Fgf18 mRNA via HCR. Gray, Hoechst; Cyan, Fgf7; Magenta, Fgf18 mRNA; Scale bars: (Top) 100 μm and (Bottom) 20 μm. The red arrow indicates the target signals. b) Zoomed region of (a) indicated by yellow box number 1 shows the target signals. Scale bars: 20 μm. c) Zoomed region of (a) indicated by yellow box number 2 shows high background signals. Scale bars: 20 μm.
	Supplementary Fig. 12 \| Quality assessment of the re-analyzed publicly available Xenopus limb regeneration dataset. a) Barplots showing the used cell numbers, the transcript number per cell, the gene numbers per cell, and the percentage of mitochondrial genes per cell after filtering in Xenopus laevis 5 days-post amputation (dpa) limb regeneration dataset are visualized.
	Supplementary Fig. 12 \| Quality assessment of the re-analyzed publicly available Xenopus limb regeneration dataset. b) UMAP plot of the identified and annotated clusters. Clusters were annotated based on marker genes listed in Supplementary Fig. 12c. c) Heatmap showing the expression profile of marker genes that are used to annotate clusters. d) Dotplot showing signaling ligand expressions in the AER or AEC clusters during Xenopus laevis limb development and regeneration, respectively. Please note that a similar figure was generated in Aztekin et al, 2021, but in this figure, NF Stage 50, 51, and 52 from Lin et al. 2021 were used.
	Supplementary Fig. 13 \| Quality assessment of re-analyzed publicly available axolotl limb regeneration dataset. a) Barplots showing the used cell numbers, the transcript number per cell, the gene numbers per cell, and the percentage of mitochondrial genes per cell after filtering in axolotl limb regeneration datasets are visualized. b) UMAP plot of the identified and annotated clusters. Clusters were annotated based on marker genes listed in Supplementary Fig. 13d. The AEC cells were identified by the analysis described in the text and in Fig. 3a-c, then re-plotted and colored dark pink in this figure. c) Sample contribution to the axolotl regeneration dataset is visualized. Red dots indicate cells from the selected sample; gray dots indicate the other cells. hpa: hours-post amputation; dpa: days-post amputation. d) Heatmap showing the expression profile of marker genes that are used to annotate clusters, except AEC cells and their marker expressions indicated following the analysis indicated in text and Fig. 3a-c.
	Supplementary Fig. 14 \| Differentially expressed gene analysis between axolotl AER and AEC cells. a) Volcano plot showing the differentially expressed genes (DEGs) between axolotl AER and AEC populations. Red and blue dots indicate genes significantly enriched in the AER or the AEC, respectively. Gray dots indicate statistically not significant genes. The top 10 DEGs are labeled. b) Barplots showing enriched GO terms based on the top 200 DEGs (ordered by fold change) in AER (red) and AEC (blue) cells.
	Supplementary Fig. 15 \| Quality assessment and clustering results of axolotl 5 days-post- amputation 10X Visium spatial transcriptomics (Visium) dataset. a) Barplots showing used Visium spot numbers, transcript number per spot, gene numbers per spot, and the percentage of mitochondrial genes per spot after filtering are visualized. b) 5 dpa axolotl limb regeneration tissue section that is used for the Visium is shown and stained for hematoxylin and eosin. Please note that this image is the same as Fig. 3f but added in this figure for easier visual comparison to Supplementary Fig. 15d, which contains full annotation. c) UMAP plot of the Visium clusters is shown. Clusters were annotated by marker gene expressions (Supplementary Fig. 15e) and confirmed by morphologies in hematoxylin and eosin staining on the tissue section in Fig. 3f and Supplementary Fig. 15b. The identified clusters are color-coded and tissue annotation is labeled. d) Visium spots were colored by tissue types identified in Supplementary Fig. 15c.
	Supplementary Fig. 15 \| Quality assessment and clustering results of axolotl 5 days-post- amputation 10X Visium spatial transcriptomics (Visium) dataset. e) Heatmap showing the top 20 DEGs (ordered by fold change) for the identified clusters is shown.
	Supplementary Fig. 16 \| Consensus nonnegative matrix factorization (cNMF) to identify AER- related modules. a) The strategy to perform consensus negative matrix factorization (cNMF) is described. cNMF was used on AER and CT clusters to identify cell identity and cell activity modules. b) Stability (blue line) is measured by the Euclidean distance silhouette score of the clustering and Frobenius error of the consensus solution (red line) for each tested k value (module number) shown in the X axis. Please see Methods for more detail. c) Violin plots showing activity scores for identified cNMF modules in CT and AER clusters. Modules specific to AER, CT, or shared are color-coded: pink, AER-specific modules; blue, CT- specific modules; green, shared modules.
	Supplementary Fig. 16 \| Consensus nonnegative matrix factorization (cNMF) to identify AER- related modules. d) Enriched GO terms are determined for the top 100 genes (ordered by gene expression program scores) in the identified AER-specific modules. Barplot showing the top GO is visualized, and GO terms are labeled.
	Supplementary Fig. 17 \| Single cell gene set enrichment analysis (scGSEA) using differentially expressed AER gene sets for individual datasets. UMAP visualization of the scGSEA of indicated datasets of limb development (a) and regeneration (b) based on the top 500 AER DEGs in multi-species limb atlas (Fig. 1b). The enrichment scores are indicated in shades of red. Cells failed to pass the enrichment threshold are colored in grey.
	Supplementary Fig. 18 \| Clustering based on differentially expressed AER gene sets for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 AER DEGs in multi-species limb atlas (Fig. 1b) for individual datasets. In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line. Please note that Mouse E10.5, chicken E4.5, human CS13, frog ST51, and axolotl ST52 development datasets, and the axolotl and frog tadpole regeneration datasets are also shown in Fig 4b and c.
	Supplementary Fig. 19 \| Clustering based on cNMF module 11 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 11 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 20 \| Clustering based on cNMF module 12 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 12 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 21 \| Clustering based on cNMF module 13 for individual datasets. UMAP visualization of the clustering of indicated datasets of limb development (a) and regeneration (b) based on the top 500 genes from AER-specific cNMF module 13 (Supplementary Fig. 16 and Supplementary Data 4). In development datasets: light blue, limb bud mesoderm cells; pink, AER cells. In regeneration datasets: dark blue, CT cells; dark pink, AEC cells. Connective tissue cells gathered with the AEC population are highlighted with a dashed line.
	Supplementary Fig. 22 \| AER gene expression profile in AER, AEC, limb bud mesoderm, and connective tissue cells during limb development and regeneration. Heatmaps showing the expression profile of 500 AER DEGs in regenerating limbs (a) and developing limbs (b). For axolotl regenerating limbs (Axolotl Reg), genes were grouped into three by the K- means algorithm, and the first group is shown in Fig. 4d. Please note that part of the axolotl regeneration dataset is also presented in Fig. 4f. Colored squares above the heatmap represent cells from indicated populations: light blue, limb bud mesoderm cells in the development dataset; light pink, AER cells in the development dataset; dark blue, connective tissue cells in the regeneration dataset; dark pink, AEC cells in the regeneration dataset.
	Supplementary Fig. 23 \| Mesodermal cells showing AER transcriptional program are present in intact limbs and express signaling ligands. a) (Top) UMAP plot from Fig. 4c and Supplementary Fig. 18b axolotl regeneration dataset is used to distinguish mesodermal cells showing part of the AER program and labeled by green. The rest of the connective tissue is labeled dark blue, and AEC is labeled dark pink. (Bottom) Dotplot showing signaling ligands expressions in the populations distinguished in Supplementary Fig. 23a. The dot color indicates the mean expression that was normalized to the max of each cell type and to the max of each gene; the dot size represents the percentage of cells with non-zero expression. b) The expression profile of previously reported regeneration-associated genes (18–20) in the UMAP plot of clustering in Fig 4c axolotl regeneration dataset. c) Sample contribution to UMAP plot of re-clustered regeneration dataset is visualized. Red dots indicate cells from the selected sample; gray dots indicate all the other cells.