Dynamic cell transition and immune response landscapes of axolotl limb regeneration revealed by single-cell analysis

The axolotl, Ambystoma mexicanum, has extraordinary capability to fully recover multiple tissues after lost, whereas such capability has disappeared in mammals. Thus, deciphering detailed mechanisms underlying axolotl regeneration could provide valuable lessons for regenerative medicine. However, many questions, such as the origin of essential progenitor cells and key responses of individual types of cells for regeneration remain elusive (Haas and Whited, 2017). Newly developed single-cell RNA sequencing (scRNA-seq) method enables researchers to observe cellular and molecular dynamics in axolotl regeneration at the single-cell resolution (Gerber et al., 2018; Leigh et al., 2018), but the reported transcriptome landscapes are only for certain cell types or in certain regenerative stages. A complete overview of the regeneration process for all cell types is still lacking. To this end, we performed axolotl limb amputation followed by scRNA-seq (10× Genomics) on upper forearm tissues at the homeostatic stage (uninjured control, 0 h), and 7 regeneration stages (3 h–33 d) (Fig. 1A, Supplementary Material). In total, we obtained high-quality sequencing data for over 41,000 cells (Table S1). After clustering all cells and reflecting them in uniform manifold approximation and projection (UMAP), we identified 12 putative cell types expressing specific markers (Figs. 1B and S1A; Table S2), which were further confirmed by differentially expressed genes (DEG) analysis (Fig. S1B). By performing integrative comparison of our dataset with ones from Gerber et al. and Leigh et al. (Gerber et al., 2018; Leigh et al., 2018), we found a remarkable similarity of cell-type distribution between our and Leigh et al.’s data (Fig. S2A and S2B), thus further supporting the reproducibility of our experimental procedure. Of note, cells expressing the embryonic limb marker Peroxiredoxin 2 (Prdx2) (Gerber et al., 2018) were identified in all three datasets (Fig. S2C). Located at the center of the UMAP plot, Prdx2 blastema cells were associated with connective tissue (CT) cells from 1 to 7 d (Fig. 1B and 1C) and expressed genes related to stem cell maintenance, cell morphogenesis and multi-tissue development pathways (Table S3), reinforcing the previous study that this cell type may originate from other types of cells by dedifferentiation, such as CT cells, and may serve as progenitor cells for redifferentiation as previously suggested (Gerber et al., 2018). More interestingly, the time-course clustering results revealed a clear inverse correlation of cell populations between epithelial, Prdx2 blastema and CT cells at two phases (first: 1 d to 7 d; second: 7 d to 22 d) (Fig. 1C and 1D). The number of basal epidermis (BE) and intermediate epidermis (IE) cells, dropped rapidly from 1 d to 7 d and then gradually returned from 7 d to 33 d (Fig. 1C and 1D). In contrast, the number of mesenchymal and Prdx2 blastema cells increased in the first phase, and then gradually returned to normal level in the second phase (Fig. 1C and 1D). Such results suggest a potential cell state transformation from epidermal cells into mesenchymal-state cells, and the second phase indicates a reversed process, while it is also possible the rates of cell proliferation/death are distinct between different cell types. To further evaluate such two possibilities, we pooled different cell populations from various regeneration time points for pseudotime trajectory analysis (Figs. 1E and S3A). While the majority of epidermal (IE and BE) and mesenchymal cell populations were located at two ends of the pseudotime trajectory axis, an evident proportion of both the IE and BE was positioned along the conversion route, and their numbers increased toward the mesenchymal-state end from 1 d to 14 d (Fig. 1F). Prdx2 blastema cells were located between the epidermal-state and mesenchymal-state ends with some overlaps (Fig. 1F), suggesting they may possess cellular characteristics of both CT cells and epidermal cells. To further investigate the cellular correlation, we clustered Prdx2 blastema, CT, BE and IE cells into 8 states along the pseudotime axis (Figs. 1G and S3B), and projected them back onto UMAP (Fig. S3C). CT cells were mainly in State 1 (S1) and S2, which are the mesenchymal-like states, while IE and BE mostly located at S7 and S8, which are the epidermal-like states. The intermediate states, S3, S4, S5 and S6, were mixed with CT, epidermal and Prdx2 blastema cells (Figs. 1G and S3B). At individual sampling time point, the cell number in the epidermal-like states continuously


Single-cell RNA-seq (scRNA-seq)
The scRNA-seq library was constructed using the Chromium single-cell 3 prime v2 reagent kit (10× Genomics) according to the manufacturer's instructions (https://support.10xgenomics.com/single-cell-gene-expression/index/doc/user-guide-c hromium-single-cell-3-reagent-kits-user-guide-v2-chemistry). All libraries were further prepared based on the requirements of the BGISEQ-500 sequencing platform manufactured by MGI ® . The DNA concentration was determined by a Qubit (Invitrogen ® ). Then, samples with 2 pmol of nucleotides were pooled to generate single-strand DNA circles (ssDNA circles). DNA nanoballs (DNBs) were generated with the ssDNA circles by rolling the circles during replication to significantly increase the fluorescent signals during the sequencing process. The DNBs were loaded into the patterned nanoarrays and sequenced on the BGISEQ-500 sequencing platform with a paired end read length of 26-100 bp.

Gene reannotation
To obtain the comprehensive genetic annotation information, we reannotated the transcriptome sequence using the protein sequence from the NCBI and UniProt. The human (GRCh38), mouse (GRCm38), gallus gallus domesticus (GRCg6a), zebrafish (GRCz11) and xenopus laevis (V2) protein sequences were downloaded from the NCBI, while the SwissProt and human marker databases were downloaded from UniProt. We then aligned all the public data against the axolotl transcriptome by using blastall (Altschul et al., 1990) (version 2.2.25) with the main parameter of "-e 1e-5". The best hit result from alignments with a Z-score≧ 200 was used to perform the homologous annotation.

scRNA-seq data processing
The average number of reads per sample of scRNA-seq data was 667,250,032, and median number of genes per cell was 1,433 (Table S1). The data were aligned to the axolotl genome(Nowoshilow et al., 2018) (https://www.axolotl-omics.org/assemblies) with STAR(Dobin et al., 2013), which was implemented in the 10× Genomics Cell Ranger-2.0.1 software with default parameters to generate the absolute Unique Molecular Identifier (UMI) counts. For the alignment, the ribosomal protein genes were excluded from all datasets. The resulting UMI counts generated the gene-cell matrices used for the downstream analysis; the matrices initially contained 23,106 unique genes from all samples and a mean of ~ 103,260 reads per cell in each sample.

Unbiased clustering with Seurat
The cell-gene matrices were filtered to retain genes expressed in at least five cells from each sample. Cells with low gene expression (< 500) were also removed. Filtered data were imported into R and integrated by using Seurat(Stuart et al., 2019) (version 3.0.3). Downstream analyses, such as normalization, shared nearest neighbor graph-based clustering, and differential expression analysis and visualization, were performed by using the R package Seurat (version 3.0.3). To enhance the identification of common cell types and enable comparative analyses of different stages, we integrated datasets from all stages by using the 'FindIntegrationAnchors' function in Seurat, which was implemented in the Seurat workflow (https://satijalab.org/seurat/v3.0/immune_alignment.html).
The clusters were identified using a community identification algorithm that was implemented in the Seurat 'FindClusters' function. The resolution parameter used to find the resulting number of clusters was set to 0.3 to produce an appropriate number of clusters that were large enough to capture most of the biological variability. The clustering results were presented using a uniform manifold approximation and projection (UMAP) based on plots generated by the 'RunUMAP' function. The 'FindAllMarkers' function was used to identify differentially expressed genes with default parameters.

DEG analysis of the regeneration stages
Each regenerative stage was compared with the control to evaluate the gene expression changes over time, and the differentially expressed genes (DEGs) were identified by the 'FindMarkers' function in Seurat (version 3.0.3). The DEGs were filtered by using |log2(foldchange)| > 1 and adjusted with a P value < 0.01. The combined mean expression of the DEGs in each cell type at different stages was calculated, and the results were shown in Figure S1B.

Single-cell trajectory construction
For EMT (epithelial to mesenchymal transition) and MET (mesenchymalto epithelial transition) related cells (intermediate epidermis, basal epidermis, connective tissue, blastema, satellite cells, and sclerotome cells), we used Monocle(Trapnell et al., 2014) (version 2.10.0) to construct the single-cell trajectory for each defined cell type with known markers. Highly variable genes with q value < 0.01 were selected, and a discriminative dimensionality reduction (DDRTree) was performed with regression according to the UMI counts to eliminate unwanted variation introduced by variation in the sequencing depth between samples, and the cell states were ordered along the trajectory according to their pseudotime values.

UMAP_2
BGI_Macrophage   Adult Axolotl Limb Regeneration. The first column is DEGs cluster number. The second column provides gene names or codes. The third to eighth columns provide average fold change (avg_logFC) and adjust P value (p_val_adj) data of listed genes between cells at different stages (1 d, 3 d, 7 d, 14 d). Stages with less than 50 Prdx2 + blastema cells were removed. The first stage with more than 50 Prdx2 + blastema cells was used as baseline stage for comparison. Data is saved as csv file.