Characterization of the Zebrafish Cell Landscape at Single-Cell Resolution

Zebrafish have been found to be a premier model organism in biological and regeneration research. However, the comprehensive cell compositions and molecular dynamics during tissue regeneration in zebrafish remain poorly understood. Here, we utilized Microwell-seq to analyze more than 250,000 single cells covering major zebrafish cell types and constructed a systematic zebrafish cell landscape. We revealed single-cell compositions for 18 zebrafish tissue types covering both embryo and adult stages. Single-cell mapping of caudal fin regeneration revealed a unique characteristic of blastema population and key genetic regulation involved in zebrafish tissue repair. Overall, our single-cell datasets demonstrate the utility of zebrafish cell landscape resources in various fields of biological research.


INTRODUCTION
The zebrafish (Danio rerio) is widely accepted as an in vivo model to study vertebrate biological processes for its similarity with human (Torraca and Mostowy, 2018). The combination of optical accessibility and genetic tractability make zebrafish an important system to investigate gene regulation and cell lineage specification. Recently, owing to the advent of high-throughput singlecell RNA sequencing (scRNA-seq), it is now possible to fully define the cell-type composition in complex biological systems (Klein et al., 2015;Macosko et al., 2015). Single-cell data resources are important for understanding cellular regulations and functions (Gupta et al., 2018). The mapping of cell atlases for whole organisms such as human, mouse, Caenorhabditis elegans larvae, planarians, and cnidarians has been achieved (Fincher et al., 2018;Han et al., 2018;Plass et al., 2018;Sebe-Pedros et al., 2018;Packer et al., 2019;Han et al., 2020). Cellular heterogeneity and developmental trajectory in zebrafish embryos have also been analyzed (Tang et al., 2017;Alemany et al., 2018;Farrell et al., 2018;Spanjaard et al., 2018;Wagner et al., 2018;Farnsworth et al., 2019). However, a comprehensive adult zebrafish cell landscape and the related molecular networks have not been fully characterized. As a powerful model, zebrafish are able to fully regenerate many tissues, such as the brain, spinal cord, kidney, heart, liver, and caudal fin (Tanaka, 2016). It remains unclear why mammals usually have limited tissue repair capacity, while lower vertebrates possess stronger ability to regenerate. Powerful molecular tools have increased our understanding of regenerative mechanisms, such as lineage tracing with transgenic lines, live imaging, and scRNA-seq (Lin et al., 2021). Profiling of Xenopus tail regeneration revealed a previously unrecognized cell type at single-cell resolution, which was associated with key regenerative pathways (Aztekin et al., 2019). Single-cell analysis of axolotl limb regeneration showed the characteristic and molecular dynamics of blastema cells (Gerber et al., 2018). Zebrafish caudal fin regeneration is also a good system to explore the process of blastema cell formation. However, the unique signatures of blastema cells and the underlying molecular mechanisms during caudal fin regeneration remain unclear.
Here, we used Microwell-seq to construct an initial compendium of a "Zebrafish Cell Landscape, " comprising more than 250,000 cells isolated from zebrafish embryo and adult tissues. Single-cell analysis of caudal fin regeneration displayed a unique characteristic of blastema population and key signaling pathways involved in zebrafish tissue repair. Our results provide a general and valuable resource for future studies on zebrafish.

Zebrafish Husbandry, Fin Amputation, and BMP Inhibitors Treatment
Zebrafish (D. rerio) wild-type Tübingen strain was raised and maintained in standard zebrafish units at Core Facilities, Zhejiang University School of Medicine. The adult zebrafish aged from 4 to 12 months were studied. The zebrafish research analysis conducted in this study was approved by the Ethics Committee of the Zhejiang University Laboratory Animal Center.
After caudal fins were amputated one or two segments proximal to the origin of bifurcation, adult zebrafish (6 months old) were kept in system water at 28.5 • C and the system water with BMP inhibitors was replaced daily. 10 µM DMH1 (Target Mol), and 10 µM Dorsomorphin (Target Mol) were used as specific BMP inhibitors. All were dissolved in DMSO and final DMSO concentration in fish water was 0.1%. The control zebrafish were kept in fish water with 0.1% DMSO. The lengths of caudal fin were detected after 3 days of BMP inhibitors treatment.

Fabrication of Microwell Device
The diameter and depth of the microwells are 28 and 35 µm, respectively. First, a silicon plate containing 100,000 microwells was manufactured by Suzhou Research Materials Microtech Co. Ltd. (Suzhou, China). The silicon microwell plate was then used as a mold to make a PDMS plate with the same number of micropillars. Prior to experiments, a disposable agarose microwell plate was made by pouring 5% agarose solution onto the surface of the PDMS plate. Both the silicon and the PDMS plates are reusable. One silicon microwell plate allows almost permanent use.

Synthesis of Barcoded Beads
Magnetic beads (diameter 20-25 µm) coated with carboxyl groups were provided by Suzhou Knowledge & Benefit Sphere Tech. Co., Ltd. (Suzhou, China 1 ). The barcoded oligonucleotides on the surfaces of the beads were synthesized by three rounds of split-pool. All the sequences used are the same as those reported previously (Han et al., 2018).
For each batch of bead synthesis, 300-350 µl of carboxyl magnetic beads (50 mg/ml) was washed twice with 0.1 M MES (2-[N-morpholino]ethanesulfonic acid). The beads were then suspended in a final volume of 635 µl 0.1 M MES. EDC [1ethyl-3-(3-dimethylaminopropyl) carbodiimide hydrochloride] (3.08 mg) was added to the beads, and 6.2 µl of beads was then placed in each well of a 96-well plate. Amino-modified oligonucleotides (2.5 µl, 50 µM in 0.1 M MES) were then added to each well. After vortexing the mixture and incubating it for 20 min at ambient temperature, 0.5 µl of mix (6 mg of EDC in 100 µl of 0.1 M MES) was distributed into each well. After an additional round of vortexing and incubation for 20 min at ambient temperature, an additional 0.5 µl of mix (6 mg of EDC in 100 µl of 0.1 M MES) was distributed into each well. After vortexing and incubation for 80 min at ambient temperature, the beads were collected in 1 ml of 0.1 M PBS containing 0.02% Tween-20. After centrifugation, the supernatant was carefully removed. The beads were then washed twice in 1 ml of TE (pH 8.0).
In the second split-pool, the beads were washed with water and divided among the wells of another 96-well plate containing PCR mix (1 × Phanta Master Mix, Vazyme) and 5 µM oligonucleotides. The oligonucleotides in each tube encoded a sequence with reverse complementarity to linker 1, a unique barcode and a linker 2 sequence. The PCR program was as follows: 94 • C for 5 min; 5 cycles of 94 • C for 15 s, 48.8 • C for 4 min, and 72 • C for 4 min; and a 4 • C hold. The third split-pool procedure was the same as the second one. The PCR program was as follows: 94 • C for 5 min, 48.8 • C for 20 min, 72 • C for 4 min and a 4 • C hold. Beads were mixed sufficiently between denaturation (95 • C) and primer annealing (48.8 • C) in every cycle. The oligonucleotides used in each tube encoded a linker 2 reverse-complementary sequence, a unique barcode, a unique molecular identifier (UMI) sequence and a poly-T tail. All the oligonucleotides were synthesized by Sangon Biotech Co. Ltd. with HPLC purification. To remove the chains without the third barcoded sequence, beads were collected and suspended in 200 µl of exonuclease I mix [containing 1× exonuclease I buffer and 1 U/µl exonuclease I (NEB)], and incubated at 37 • C for 15 min (beads were mixed by rotary mixer). After being washed with 200 µl of TE-TW and 200 µl of 10 mM Tris-HCl pH 8.0, beads were resuspended in 1 ml of ddH 2 O. To remove complementary chains, the beads were placed in a 95 • C water bath for 6 min and separated using a magnet, removing the supernatant quickly, 2 times. The beads could be stored in TE-TW (10 mM Tris pH 8.0, 1 mM EDTA, 0.01% Tween-20) for 4 weeks at 4 • C.

Cell Preparation
Dissociation of embryonic tissues was performed similarly as previously described (Manoli and Driever, 2012). In brief, 50-100 zebrafish embryos were grown to the indicated times and chorions were removed by incubating in 1mg/ml Pronase (Sigma) for 3-4 min followed by washing in 0.5× Danieau Buffer.
[10× Danieau Buffer = 174 mM NaCl, 2.1 mM KCl, 1.2 mM MgSO4, 1.8 mM Ca(NO3)2, 15 mM HEPES, pH 7.6]. Yolk were removed by blowing in deyolking buffer [55 mM NaCl, 1.8 mM KCl, 1.25 mM NaHCO3] for 10 times followed by washing in 0.5× Danieau Buffer. Embryo tissues were triturated to homogeneity in 1-5 ml FACSmax cell dissociation solution (AMS Biotechnology) and incubated for 4-5 min at room temperature. Then embryonic single cells were collected after passage through a 40-µm strainer (Biologix). Zebrafish tissues were carefully dissected, put to cold DPBS and minced into ∼1-mm pieces on ice using scissors. The tissue pieces were transferred to a 15-ml centrifuge tube, rinsed twice with cold DPBS and suspended in 5 ml of a solution containing dissociation enzymes. The samples were treated with various enzymes at 37 • C for different amounts of time (Supplementary Table 1). During the dissociation, the tissue pieces were pipetted up and down gently several times until no tissue fragments were visible. The dissociated cells were centrifuged at 300 × g for 5 min at 4 • C and then re-suspended in 3 ml of cold DPBS. After passage through a 40-µm strainer (Biologix), the cells were washed twice, centrifuged at 300 × g for 5 min at 4 • C, and re-suspended at a density of 1 × 10 5 cells/ml in cold DPBS containing 2 mM EDTA.

Cell Collection and Lysis
Cell concentration should be carefully controlled during Microwell-seq. Both cell and bead concentrations were estimated using a haemocytometer. The proper cell concentration is ∼100,000/ml (with 10% of the wells occupied by single cells). The proper bead concentration is ∼1,000,000/ml (with every well occupied by single beads). An evenly distributed cell suspension was pipetted onto the microwell array, and extra cells were washed away. To eliminate cell doublets, the plate was inspected under a microscope. Cell doublets were reduced by pipetting over the region of high cell density. The bead suspension was then loaded into the microwell plate, and the plate was placed on a magnet. Excess beads were washed away slowly. Cold lysis buffer (0.1 M Tris-HCl pH 7.5, 0.5 M LiCl, 1% SDS, 10 mM EDTA, and 5 mM dithiothreitol) was pipetted over the surface of the plate and removed after 12 min of incubation. The beads were then collected, transferred to an RNase-free tube, and washed once with 1 ml of 6 × SSC, once with 500 µl of 6 × SSC and once with 200 µl of 50 mM Tris-HCl pH 8.0. Finally, ∼50,000 beads were collected in a 1.5-ml tube.

Exonuclease I Treatment
The beads were washed with 200 µl of TE-TW and 200 µl of 10 mM Tris-HCl pH 8.0, resuspended in 100 µl of exonuclease I mix containing 1× exonuclease I buffer and 1 U/µl exonuclease I (NEB), and incubated at 37 • C for 60 min with mixing on a rotary mixer to remove oligonucleotides that did not capture mRNA. The beads were then pooled and washed once with TE-SDS, once with 1 ml of TE-TW and once with 200 µl of 10 mM Tris-HCl pH 8.0.

cDNA Amplification
The beads were distributed into 4 PCR tubes. To each tube, 12.5 µl of PCR mix [1 × HiFi HotStart Readymix (Kapa Biosystems) and 0.1 µM TSO_PCR primer] was added. The PCR program was as follows: 98 • C for 3 min; 4 cycles of 98 • C for 20 s, 65 • C for 45 s, and 72 • C for 6 min; 10-14 cycles of 98 • C for 20 s, 67 • C for 20 s, and 72 • C for 6 min; 72 • C for 10 min; and a 4 • C hold. After pooling all PCR products, AMPure XP beads (Beckman Coulter) were used to purify the cDNA library.

Transposase Fragmentation and Selective PCR
The purified cDNA library was fragmented using a customized transposase that carries two identical insertion sequences. The customized transposase was included in the TruePrep DNA Library Prep Kit V2 for Illumina (Vazyme). The fragmentation reaction was performed according to the instructions provided by the manufacturer. We replaced the index 2 primers (N5××) in the kit with our P5 primer to specifically amplify fragments that contain the 3 ends of transcripts. Other fragments will form selfloops, impeding their binding to PCR primers. The PCR program was as follows: 72 • C for 3 min; 98 • C for 30 s; 5 cycles of 98 • C for 15 s, 60 • C for 30 s, and 72 • C for 3 min; 72 • C for 5 min; and a 4 • C hold. The PCR product was purified using AMPure XP beads. Then, 25 µl of PCR mix (1× HiFi HotStart Readymix and 0.1 µM 2100 primer) was added to each sample. The PCR program was as follows: 95 • C for 3 min; 8 cycles of 98 • C for 20 s, 60 • C for 15 s, and 72 • C for 15 s; 72 • C for 5 min; and a 4 • C hold. To eliminate primer dimers and large fragments, AMPure XP beads were then used to purify the cDNA library. The size distribution of the products was analyzed on an Agilent 2100 bio-analyzer, and a peak in the 400-700 bp range was observed. Finally, the samples were subjected to sequencing on an Illumina Hiseq system.

Processing of the Microwell-Seq Data
Standard procedures for processing the Microwell-seq datasets were performed using the protocols described in the previously published paper (Han et al., 2018). Reads from zebrafish cell landscape data were aligned to the D. rerio GRCz10 genome using STAR (Dobin et al., 2012) and the digital gene expression (DGE) data matrixes were obtained using the Dropseq Core Computational Protocol 2 with default parameters. For quality control, we filtered out cells with detection of less than 500 transcripts. Cells with high proportion (>20%) of transcript counts derived from mitochondria-encoded genes were also excluded.

Clustering of Single Cell Data Matrix
Seurat (V3.0)  was used to perform clustering analysis of single cell data from different tissues. The DGE data was used as inputs. Cells from the pre-processed data and genes expressed in more than 3 cells were selected for further analysis. The filtered data was log 2 (TPM/100 + 1) transformed, then the number of UMI and the percentage of mitochondrial gene content were regressed out according to published method (Buettner et al., 2015). About 2,000 various genes with average expression more than 0.5 and a dispersion greater than 0.125 were used as inputs for initial principal component analysis (PCA) and number of principal components (PCs) used for Non-linear Dimensional Reduction (t-SNE) analysis is chosen according to the PCElbowPlot function and JackStrawPlot function. For clustering, we set different resolution parameters between 0.6 and 4 in FindAllCluster function and narrowed down to certain cluster numbers by distinguishing differential genes among clusters. The heatmap, produced by DoHeatmap function is one of basis for judging the quality of clustering. These parameters, including resolution and number of principal components, were adjusted on per-tissue basis. The default wilcox rank sum test was used by running FindAllMarkers function in Seurat to find DEGs in each cluster. Gene ontology (GO) analysis of DEGs was performed with ToppGene 3 . Finally, we annotate each cell type by extensive literature reading and searching for the specific gene expression pattern.
For processing of the complete zebrafish tissue dataset (258,902 cells), we used Scanpy  in python environment to perform the analysis. Background-removed DGE data with cells analyzed in each tissue and genes expressed in at least 20 cells was the used as inputs for Scanpy. Then the DGE data was log 2 (TPM/10 + 1) transformed. We selected about 3,000 highly-variable genes according to their average expression and dispersion, regressed out UMI and gene numbers, scaled each gene to unit variance with clip values exceeding a standard deviation of 10. We chose about 50 PCs for the PCA, and we computed the neighborhood graph of cells. We then used the Louvain clustering to cluster cells with resolution = 2.5 and k = 10. Finally, 63 clusters for the zebrafish landscape were produced and marker genes were calculated by Wilcoxon rank sum test.

Cross-Species Transcriptome Comparison
To make the gene expression profiles of cross-species cell types comparable, we downloaded the homology correspondences between zebrafish and human, mouse, axolotl provided by dmodENCODE (Celniker et al., 2009). The gene expression profiles for zebrafish (this study) and other species (Gerber et al., 2018;Han et al., 2018Han et al., , 2020 were normalized to the total number of transcripts and multiplied by 100,000. To attenuate the effects of noise and outliers, we used pseudo-cells (Tosches et al., 2018) for further analysis; each pseudo-cell was an average of 20 cells randomly selected from the same cell type. To compare cross-species transcriptomes, we performed MetaNeighbor (Crow et al., 2018) analysis between zebrafish and mouse. MetaNeighbor was formalized through neighbor voting based on cell-cell similarities. The mean area under the receiver operator characteristic curve (AUROC) scores was used to measure the similarity of cell types, and we chose 0.8 as the threshold. The Circlize package was used to view in the similarity of cell types among different species (Gu et al., 2014). The network was visualized using Cytoscape with the "edge-weighted spring embedded" layout (Shannon et al., 2003).

Weighted Correlation Network Analysis
WGCNA (V1.69) was performed with default parameters (Langfelder and Horvath, 2008). In brief, high variable genes among the caudal fin regeneration single cells of interest were calculated using the distance to the median metric. An adjacency matrix, representing a "unsigned" gene network, was built setting the soft power parameter to 4 (calculated from the pick Soft Threshold function). Modules were identified using the function block wise Modules, with a minimum module size of 30 and the save TOMs parameter set to TRUE. The TOM dissimilarity measure (1-TOM) to the fourth power was then used to cluster genes with the function TOM plot. And the correlation of modules and gene expressed in each cell type was calculated by Pearson correlation coefficient.
Virtual Inference of Protein-Activity by Enriched Regulon Analysis VIPER (V1.20.0) was performed to computational inference of TF activity from gene expression profile data with default parameters (Alvarez et al., 2016). Briefly, an appropriate cell context-specific regulatory network was generated with Frontiers in Cell and Developmental Biology | www.frontiersin.org ARACNE algorithm 4 from caudal fin regeneration and nonregeneration single-cell gene expression profiles (Margolin et al., 2006). We then performed master regulator analysis and estimated its significance, including P-value and normalized enrichment score (NES), by comparing each regulon enrichment score to a null model using msVIPER algorithm 5 . Master TFs and cofactors regulated each other and their downstream target genes were visualized by Cytoscape (V3.8.2).

Gene Set Enrichment Analysis
Signature enrichment of regeneration relevant clusters against others was performed by gene set enrichment analysis (GSEA) (V1.0). The GSEA was implemented using JAVA downloaded from the Broad Institute 6 . D. rerio pathway database for enrichment was downloaded from https://www.wikipathways. org/index.php/Download_Pathways.

Receptor-Ligand Pairing Analysis
Analysis of potential receptor-ligand pairings was performed using the method CellPhoneDB (V1.0.0) (Vento-Tormo et al., 2018). First, we aggregated the gene expression levels of 20 cells from each cluster in the caudal fin regeneration. To eliminate the effect of variable cell numbers in each cluster, we randomly sampled three pseudo-cells for analysis. Only receptors and ligands expressed in more than 10% of the cells in the specific cluster were considered. By permuting cluster labels randomly 1,000 times to calculate the mean expression values of ligands and receptors, interaction was constructed as a receptor-ligand pairing matrix. Then, we used pairwise comparisons between all cell types and obtained a likelihood of P value to filter the falsepositive interaction. The cut off was set with the mean expression greater than 0.1 and P values smaller than 0.1. We used the sum of the number of receptor-ligand pairs in each cell-cell pairing to indicate the strength of the cell-cell interactions.

Statistical Analysis
All the results were presented as mean ± SEM. Two-tailed Student's t-test was performed to compare the differences between two groups. All quantitative experiments were repeated at least 3 times independently.

Constructing a Zebrafish Cell Landscape Using Microwell-Seq
Single-cell RNA sequencing technologies enable identification of cell identity. To comprehensively characterize different cell types in zebrafish, we performed mRNA-seq expression profiling in individual cells isolated from major tissues using Microwell-seq (Han et al., 2018) (Figure 1A). We collected the pharyngula stage (24 h post-fertilization, hpf) and hatching stage (72 hpf) zebrafish embryos, as well as blood, brain, caudal fin, eye, gill, heart, intestine, kidney, liver, muscle, ovary, pancreas, skin, spleen, swim bladder, and testis samples from adult zebrafish. Tissues were obtained, carefully dissected and prepared into single-cell suspension (Supplementary Table 1). Single cells were then processed with single-cell mRNA-seq.
The single-cell transcriptomics data were processed by published pipelines (Macosko et al., 2015;Satija et al., 2015). In general, we analyzed >250,000 single cells with an average 537 genes and ∼1,180 UMIs per cell from zebrafish embryos and 16 adult tissues (Supplementary Figure 1A). To detect relationships between cells from different organs, we visualized all cells with t-distributed stochastic neighbor embedding (t-SNE) and grouped them with unbiased, graph-based clustering ( Figure 1B). In a global view, we identified tissue cell types in 63 major clusters (Supplementary Table 1). Experimental batches from the same tissues were well controlled ( Figure 1C). We elucidated the cell type identity of clusters by examining marker genes and comparing them to those in previous studies. For example, cluster 7 (C7) and C9 were defined as hepatocyte with high expression of leg1.1 and bhmt. C50 partly derived from zebrafish embryos and was identified as primitive hepatocyte with additional high level of dao.1 and acot17. Compared with C7, C9, and C50 showed a closer spatial relationship in the t-SNE map. Similar to human and mouse liver cell atlas, the heterogeneity of zebrafish hepatocyte may imply zonation features at the single-cell level (Aizarani et al., 2019). Epithelial cells (C6, C8, C14, C28, and C35) were composed of multiple tissues such as caudal fin, eye, gill, muscle and skin, suggesting the consistency of transcriptome signatures of epithelial cell in different tissues. Conversely, innate immune cells of skin (C58) and intestine (C63) were identified as two independent subgroups, indicating specific tissue residency and function execution (Figures 1B,C and Supplementary Figure 1B). We then performed subclustering analysis for each of the 63 major clusters and generated a hierarchy of more than 600 cell-type sub-clusters ( Figure 1D). Together, we clustered the expression profiles of the individual cells and constructed a comprehensive cell type landscape for zebrafish.

Cellular Heterogeneity in Embryo and Adult Tissues
t-SNE analysis and differential gene expression analysis for each specific organ type were performed (Figure 2, Supplementary  Figure 2, and Supplementary Table 2). The brain is the most complex organ that serves as the center of the nervous system in all vertebrates, controlling over the other organs of the body (Raj et al., 2018). Similar to the previous singlecell analysis of adult zebrafish brain (Alemany et al., 2018;Spanjaard et al., 2018), our data identified four major cell groups. Two neuron sub-clusters (C7 and C15) expressed high level of snap25a, syt1a and kiss1. C6 and C12 were defined as radial glia expressing cx43 and fabp7a. Gene expression of C1 and C19 encoded the classic myelin associated protein, a major constituent of the myelin sheath of oligodendrocyte in the nervous system. C3 and C14 were identified as microglia with high expression of apoeb (Figures 2A,B). Notably, our data additionally identified neuro precursor cells, such as C9 with high level of neurod1, hes6 (Weng et al., 2019) and C10 as quiescent radial glia expressing ascl1a, her15.2 and her4.1 (Chapouton et al., 2011). Mural lymphatic endothelial cells (C20) were defined with high level of lyve1b (Supplementary Figure 3A). Functional experiments verified its biological significance of regulating meningeal angiogenesis in zebrafish brain (Bower et al., 2017). By analyzing adult zebrafish eye, we defined 12 distinct clusters with specific molecular markers ( Figure 2C). Similar to previous single-cell analysis of adult eye (Alemany et al., 2018), retinal rod cells (C2) and retinal cone cells (C11) were identified. Additionally, we defined C10 as müller glia, a type of retinal glial cell with high expression of ascl1a, glula and rlbp1a (Supplementary Figure 3C), which would be induced to dedifferentiate and produce multipotent neuronal progenitor cells in damage retina models (Raymond et al., 2006). Epithelial cells (C1 and C7), stromal cells (C6 and C8) and immune cells (C5 and C12) were also detected in zebrafish eye single-cell data ( Figure 2D).

Cross-Species Analysis of Cell-Type Similarity
Single-cell transcriptomics offers an opportunity for comprehensive cross-species and cross-tissues analysis of cell types. The lung is the primary organ for gas exchange in mammals. The swim bladder and gill are specialized organs in teleosts that regulate respiration. Whether the lung evolves from the gill or the swim bladder is still a controversy based on morphological evidences (Perry and Sander, 2004;Zheng et al., 2011). A recent study reported that the cell types could be proposed as "evolutionary units" in comparative cell biology (Wang et al., 2021). To infer the evolutionary relationship of respiratory system at the single-cell level, we extracted orthologous genes and performed cross-species clustering analysis of zebrafish gill (Z_G), swim bladder (Z_SB) and mouse adult lung (M_AL), fetal lung (M_FL) (Han et al., 2018) (Figures 3A,B and Supplementary Table 3). As revealed by the heatmap and circos plot, the gene expression patterns of major cell types showed strong correlations between mouse and zebrafish, such as immune cells, stromal cells and proliferating cells (Figure 3C and Supplementary Figure 4A). Notably, based on transcriptome distances that were evaluated by correlation coefficient, we found that adult mouse lung alveolar type 1 (AT1) cells specialized for gas exchange displayed strong correlations with zebrafish swim bladder epithelial cells, while other epithelial cells shared transcriptional similarity with gill ionocytes (Figure 3D and Supplementary Figure 4B). Recent studies identified Foxi1+ pulmonary ionocytes in mouse lung, a rare cell type in the conducting airway at single-cell resolution (Plasschaert et al., 2018). We propose that a gene's function is strongly predictive of conservation in gene expression program; the evolution of cell-type regulatory network may be independent of tissue evolution. Taken together, our single-cell transcriptomes demonstrated the conservation and divergence of cell types between zebrafish swim bladder, gill and mouse lung, providing potential resources to infer the cell-type evolutionary relationship.
Unlike mammals, zebrafish adult tissues exert powerful regenerative ability (Gemberling et al., 2013). To understand the potential cellular mechanism, we calculated the correlation of cell-type between the zebrafish cell landscape and human cell landscape, mouse cell atlas, axolotl limbs regenerative landscape (Han et al., 2018(Han et al., , 2020Leigh et al., 2018). scRNA-seq datasets revealed a well-conserved cellular architecture that enables matching of homologous cell types between zebrafish and other species, such as embryonic cell, oligodendrocyte, enterocyte, epithelial cell and granulocyte etc. (Supplementary Table 3). Notably, we found that zebrafish adult stromal cells showed a strong correlation with both human and mouse fetal stromal cells, resembling zebrafish embryonic stromal cells ( Figure 3E and Supplementary Figure 4C). In cross-species clustering analysis, adult zebrafish stromal cells were clearly close to fetal and neonatal mammalian stromal cells ( Figure 3F). Moreover, zebrafish adult stromal cells were strongly associated with fibroblast-like blastema in axolotl limbs regenerative landscape (Supplementary Figure 4D). Taken together, adult zebrafish stromal cells are intrinsically different from the adult human and mouse stromal cells. They possess a progenitor-like phenotype that is only seen in fetal mammalian tissues. This may help to explain stronger regenerative potentials in the adult zebrafish tissues when compared to higher organisms.

Characteristic of Blastema Cells During Caudal Fin Regeneration
To investigate the regenerative capability of the zebrafish system, we then focused on zebrafish caudal fin ontogeny and  its regeneration process. scRNA-seq was performed to profile uninjured fin and outgrowth of regenerated blastema cells at 3 days post-amputation (dpa) from two independent biological replicates termed replicate 1 and replicate 2 (Poss et al., 2003) ( Figure 4A). An independent analysis of the two biological replicates revealed that the results of both experiments strongly overlapped. Transcriptomic profiling of uninjured (0 dpa) and regenerated fin (3 dpa) contributed to nearly all cell clusters, suggesting that the cell identity was unaffected by batch effect (Figure 4B). Three major cell types (osteoblast, epidermal and blastema-like cell) were identified based on gene expression signatures (Supplementary Table 4). Among them, blastema-like cells were mainly composed of four subpopulations. In replicate 1, C4 was identified as stromal-related cells with high expression of col1a1a and col1a1b (C7 in replicate 2), C8 expressed high level of fibroblast growth factor binding protein (fgfbp2a) (C10 in replicate 2), C7 and C13 differed mainly in higher expression of proliferation genes (C4 and C8 in replicate 2), such as hmgn2, pcna ( Figure 4C).
To further reveal the unique characteristic of blastemalike cells, we next analyzed the differentially expressed genes (DEGs) between uninjured (0 dpa) and regenerated fin (3 dpa). Upregulated markers in blastema-like cells during caudal fin regeneration were identified ( Figure 4D). Extracellular matrix genes, such as col1a1a and col1a2, showed significantly enhanced expression, indicating that stromal cells mainly participated in blastema formation. The result is in accordance with previous findings that connective tissues transit to a blastema state during adult axolotl and frog limb regeneration (Gerber et al., 2018;Lin et al., 2021). Besides, blastema were enriched in cells expressing c1qtnf5, ecrg4a, and1, clu, and fgfbp2a, with confirmation of their localization by in situ hybridization (Figures 4E,F and Supplementary Figures 5A,B). We found that c1qtnf5 and ecrg4a were only expressed in outgrowth blastema involved in regulation of cell proliferation, which is consistent with zebrafish mantle cells (hair-cell progenitors) that express high levels of c1qtnf5 and ecrg4a during hair cell regeneration (Steiner et al., 2014). clu, a marker of mouse intestine revival stem cell (Ayyaz et al., 2019), exhibited potential regulation function during zebrafish tissue repair. Together, our data illuminated molecular characteristics of blastema cells during zebrafish caudal fin regeneration at singlecell resolution.

Genetic Regulation During Tissue Regeneration
Transcription factors (TFs) directly interpret the genome, exerting control over processes that specify cell types and controlling specific pathways (Lambert et al., 2018). To understand the genetic regulation during zebrafish caudal fin regeneration, we next focused on the activity of TFs by enriched regulon analysis. We performed weighted correlation network analysis (WGCNA) to look for clusters (modules) of highly correlated genes (Langfelder and Horvath, 2008). We constructed a gene co-expression network in zebrafish caudal fin regeneration single-cell dataset. Blastema-like cells (C4, C7, C8, and C13 in replicate 1) were integrated into one regeneration gene module, and other clusters were defined as non-regenerative module (Supplementary Figure 6A). Virtual inference of protein-activity by enriched regulon analysis (VIPER) was then performed to infer the relative activity of a regulatory TF based on the enrichment of its most closely-regulated targets on a given gene expression signature (Alvarez et al., 2016). Comparing the regeneration and non-regenerative gene modules, we identified high activity TFs regulating regeneration progress, such as twist1a, prrx1a, msx1b, fosab, etc., as well as low activity TFs such as tp63, elf3 and klf2b (Figures 5A,B, Supplementary  Figure 6B, and Supplementary Table 5). Previous studies found that twist1a and prrx1a play a vital role in axolotl limb regeneration, msx1b and fosab are essential for cell proliferation and differentiation (Gerber et al., 2018;Hirsch et al., 2018).
Next, we sought to look for key signaling pathways that regulate tissue regeneration. In the gene enrichment analysis, bone morphogenetic protein (BMP) signaling pathway was found to be strongly associated with the zebrafish caudal fin regenerative process ( Figure 5C and Supplementary Figure 6C), exhibiting significant difference between regenerative and non-regenerative module ( Figure 5D and Supplementary Figure 6D). The BMP signaling pathway was mainly regulated by blastema cells, in which bmp4 was only expressed in 4 clusters of blastema, bmp2a was enriched in C4 and C8, and bmp1a was highly expressed in C7 and C13 ( Figure 5E). Besides, ligand-receptor map was constructed to reveal cell-cell interactions during caudal fin regeneration (Supplementary Figure 6E). We found that the 4 clusters of blastema cell were at the center of the network, regulating the regeneration process. Previous studies have reported that BMP signaling pathway contributed to zebrafish cardiomyocyte regeneration (Wu et al., 2016) and whole body regeneration in acoels (Srivastava et al., 2014). To further examine the function of BMP signaling during caudal fin repair, adult zebrafish were treated with two specific BMP inhibitors, DMH1 and dorsomorphin (Dor), 12 h before amputation (−12 h) to 3 dpa ( Figure 5F). The results showed that DMH1 and Dor significantly inhibited caudal fin regeneration compared to DMSO, suggesting that the BMP signaling pathway plays crucial roles during zebrafish caudal fin repair (Smith et al., 2006) (Figures 5G,H). Summarily, insights into genetic regulation and signaling pathways involved in caudal fin regeneration have clear implications for future prospects in tissue regenerative engineering.

DISCUSSION
Here, we used Microwell-seq to generate a zebrafish cell landscape. At the current stage, although the sequencing relatively shallow, the zebrafish cell landscape data certainly allows for the separation of major cell types in the zebrafish system. After the cell clustering, reads from the same cell type can then be aggregated for deeper investigation of genetic regulation. Single-cell transcriptome analysis has already been  applied in zebrafish embryo, larvae, and in brain, kidney, eye, caudal fin, heart, liver and pancreas of adult fish (Tang et al., 2017;Spanjaard et al., 2018;Farnsworth et al., 2019). When compared to available single-cell data of zebrafish tissues, we have covered a more comprehensive tissue types in terms of the breath of the analysis. In pharyngula stage (24 hpf) singlecell analysis, neural cell, hatch gland, epithelial cell, muscle, ionocyte, macrophage, lens and otic were identified in both our study and a previously published work . Moreover, our single-cell data additionally identified xanthophore and hepatocyte which is missing in Wagner's research. In zebrafish kidney, Tang et al. (2017) identified 14 clusters of kidney and hematopoietic cells, while we identified 18 clusters in our microwell-seq dataset of zebrafish kidney. Among these clusters, nephron epithelial cell, proximal tubular cell, distal tubular cell, vascular endothelial cell, mucin cell, HSC, macrophage, erythrocyte, and blood progenitors are highly comparable. The proliferating cell and innate immune cell identified in our study correspond to the kidney progenitor cell and lymphoid cell in Tang's research, respectively. Interestingly, we found a fraction of immune-active epithelial cells expressing high level of mal in zebrafish caudal fin, skin, and swim bladder, indicating an important role of epithelial cells involved in immune responses (Supplementary Figure 3D) (Schleimer et al., 2007;Han et al., 2020). Our work is by no mean a complete representation of all zebrafish cell types, but we constructed an initial draft to create an organism-wide cellular hierarchy for the adult zebrafish.
Single-cell transcriptomics offers an opportunity for comprehensive cross-species and cross-tissues analysis of cell types. A recent study reported that the cell-type could be proposed as "evolutionary units" in comparative cell biology (Wang et al., 2021). In the current study, we propose that the different epithelial cells in mammal lung may evolve from different organs in zebrafish including swim bladder and gill. In other words, the evolution of cell-type regulatory network may be independent and superior to tissue evolution.
Unlike mammals, zebrafish is a model system with an amazing capacity for regeneration (Marques et al., 2019). To understand the potential cellular mechanism, we performed cross-species analysis of cell-type similarity between the zebrafish cell landscape and human cell landscape, mouse cell atlas. The results showed that adult zebrafish stromal cells shared strong similarity to human and mouse fetal stromal cells. This may help to explain stronger regenerative potentials in the adult zebrafish tissues when compared to higher organisms.
During caudal fin regeneration, cells near the amputation plane accumulate into a distinctive tissue called the blastema. Blastema formation is a comprehensive process, comprising various different cell dedifferentiation, proliferation and redifferentiation to rebuild the missing fin structures (Pfefferli and Jaźwińska, 2015). In this study, our single-cell analysis of blastema population showed unique transcriptional signatures and key signaling pathways involved in caudal fin repair. The two biological replicates show high similarities in both blastema gene expression and TF regulation. Previous study revealed that zebrafish mantle cells expressed high level of c1qtnf5 and ecrg4a during hair cell regeneration (Steiner et al., 2014). Similarly, we found c1qtnf5 and ecrg4a were involved in the formation of blastema during caudal fin regeneration, and verified it by in situ hybridization ( Figure 4F). Focused on the genetic regulation, we found previously identified key blastema marker genes, including the muscle segment homeobox family member msx1b and bone development and regeneration TF twist1a (Hou et al., 2020). During axolotl limb regeneration, there is a population of relatively homogenous progenitor cells similar to the embryonic state (Gerber et al., 2018). Hmgb is a well-conserved nuclear protein and plays vital role in the development of zebrafish pectoral fin buds and mouse forelimb buds, which also regulates caudal fin regeneration found in our study (Itou et al., 2011). However, cx43, a gap junction protein required to build the right fin length, was not identified in this study ( Figure 5B). We believe that combination of other molecular tools, such as lineage tracing and transplantation-based functional assays, will further broaden our knowledge on the origin of blastema cells, the initial signals of cell regeneration, and the impact of microenvironment in the near future (Lin et al., 2021). A deep knowledge on the single-cell data of zebrafish tissue repair would inspire new strategies for controlling tissue regeneration in mammals.

CONCLUSION
We present a zebrafish cell landscape with single-cell composition for many tissues that have not been well characterized. We reveal a unique molecular and cellular phenotype during caudal fin regeneration. Our single-cell datasets improve our understanding of the zebrafish and provide a valuable reference for future studies.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by the Ethics Committee of the Zhejiang University Laboratory Animal Center.

AUTHOR CONTRIBUTIONS
GG, XH, and JP were responsible for the conception, design, and study supervision. MJ, HC, CG, and YL were responsible for the development of methodology, analysis, and experiments. YX, WE, LM, JW, and QG performed the bioinformatics analysis and constructed the website of zebrafish cell landscape. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We would like to thank Zebrafish Resource Center (Hangzhou) for supporting this project. We thank Vazyme (Nanjing) for supplying customized enzymes in the study. We thank for technical support by the Core Facilities of Zhejiang University School of Medicine and computational support by Center of Cryo-Electron Microscopy at Zhejiang University. We thank JP and LM for their guidance. GG is a participant of the Human Cell Atlas Project (International), the Alliance for Atlas of Blood Cells (China), and the Cell Atlas Project (Zhejiang University Stem Cell Institute).

SUPPLEMENTARY MATERIAL
The Sheet 3: The number of sub-cluster in 63 zebrafish cell types.
Sheet 4: The oligonucleotide sequences used in single cell mRNA-seq library construction.
Supplementary Table 2 | Cellular heterogeneity in embryo and adult tissues. Differentially expressed genes detected in each cell type for all tissues in zebrafish cell landscape datasets. Yellow labels indicate specific marker genes of cell clusters. Genes are selected by log foldchange > 0.25, Bonferroni-adjusted p-value < 0.1, expressed in at least 15% of cells in either population (Seurat FindAllMarkers). Log fold change is calculated as arithmetic mean of log10 cpm values of one population minus the arithmetic mean of log10 cpm values of the second, and fold change is 10log_foldchange. P-values were calculated by the Wilcoxon rank sum test. The last sheet showing cell barcodes of each cell in zebrafish tissues.
Supplementary Table 3 | Cross-species analysis of cell-type similarity.
Sheet 1: Cell type annotation of zebrafish gill, swim bladder and mouse lung.
Sheet 2: The correlation coefficient of cell types among zebrafish gill, swim bladder and mouse lung.
Sheet 3: The correlation of cell-type among zebrafish cell landscape and human, mouse cell landscape.
Sheet 4: The correlation of cell-type among zebrafish cell landscape and axolotl limbs regenerative landscape.
Supplementary Table 4 | Characteristic of blastema cells during caudal fin regeneration. Differentially expressed genes detected in each cell type for caudal fin regeneration. Genes are selected by log foldchange > 0.25, Bonferroni-adjusted p-value < 0.1, expressed in at least 15% of cells in either population (Seurat FindAllMarkers). Log fold change is calculated as arithmetic mean of log10 cpm values of one population minus the arithmetic mean of log10 cpm values of the second, and fold change is 10log_foldchange. P-values were calculated by the Wilcoxon rank sum test.
Supplementary Table 5 | Genetic regulation during tissue regeneration.
Sheets 1 and 2: The transcription factor-activity of caudal fin regeneration module.
Sheet 3: Activated and repressed transcription factors in caudal fin regeneration module, as well as its target genes.
Sheets 4 and 5: Gene set enrichment analysis of caudal fin regeneration module.
Frontiers in Cell and Developmental Biology | www.frontiersin.org