m6A promotes planarian regeneration

Abstract Regeneration is the regrowth of damaged tissues or organs, a vital process in response to damages from primitive organisms to higher mammals. Planarian possesses active whole‐body regenerative capability owing to its vast reservoir of adult stem cells, neoblasts, providing an ideal model to delineate the underlying mechanisms for regeneration. RNA N 6‐methyladenosine (m6A) modification participates in many biological processes, including stem cell self‐renewal and differentiation, in particular the regeneration of haematopoietic stem cells and axons. However, how m6A controls regeneration at the whole‐organism level remains largely unknown. Here, we demonstrate that the depletion of m6A methyltransferase regulatory subunit wtap abolishes planarian regeneration, potentially through regulating genes related to cell–cell communication and cell cycle. Single‐cell RNA‐seq (scRNA‐seq) analysis unveils that the wtap knockdown induces a unique type of neural progenitor‐like cells (NP‐like cells), characterized by specific expression of the cell–cell communication ligand grn. Intriguingly, the depletion of m6A‐modified transcripts grn, cdk9 or cdk7 partially rescues the defective regeneration of planarian caused by wtap knockdown. Overall, our study reveals an indispensable role of m6A modification in regulating whole‐organism regeneration.


| INTRODUCTION
Regeneration is a process of replacing or restoring damaged or missing cells and tissues to full function, which widely exists in the animal kingdoms. 1 For example, hydra can regenerate from tiny body fragments to entire organisms, 1 and zebrafish can regenerate large portions of the heart. 2 For humans, some tissues also have regenerative potential throughout life, including the liver and muscle. 3 Planarians are being considered a popular model system to study regeneration at the whole-organism level. 4,5 The discovery of the planarian stem cell (neoblast) marker piwi (smedwi-1), 6 the development of neoblast isolation 7 and RNA interference (RNAi) knockdown 8 methods have opened up opportunities to unravel the mechanisms of planarian regeneration. A number of pathways and regulatory factors have been identified to be essential for planarian regeneration, such as Wnt 9-14 and EGFR signalling. 15 Especially along with the advances of third-generation DNA 16 and large-scale single-cell RNA sequencing technologies, 17,18 accumulating evidence suggests that the pluripotent neoblasts are cellular sources for regeneration. 19,20 In addition, some studies have revealed that epigenetic modifications also exert functions in planarian regeneration. For instance, the COMPASS family MLL3/4 histone methyltransferases are essential for the differentiation and regeneration of the planarian. [21][22][23] In addition, the CREBbinding protein (CBP) and p300 family of histone acetyltransferases homologues Smed-CBP2 and Smed-CBP3 displayed distinct roles in stem cell maintenance and functions. 24,25 N 6 -methyladenosine (m 6 A), the most abundant dynamic internal mRNA modification, has been shown to be an epi-transcriptomic marker playing diverse regulatory roles under physiological and/or pathological conditions. 26 Studies in various model systems have revealed that m 6 A regulates stem cell self-renewal and differentiation. [27][28][29][30] Depletion of either m 6 A methyltransferases or its demethylases dramatically affects gene expression profiles. [31][32][33][34] High-throughput sequencing technologies enable the detection of the m 6 A location within the transcriptome. 31,35 Subsequently, m 6 A has been identified in various RNA species and has been implicated in stem cell biology, developmental and cancer biology. 26 For example, two separate studies showed that m 6 A regulates the transition of embryonic stem cells from a pluripotent state to a differentiated state. In this case, m 6 A selectively marks transcripts that code for key transcription factors involved in differentiation. 27,36 This demonstrates that m 6 A is a molecular switch that regulates stem cell differentiation, a fundamental mechanism in development and stem cell biology. Moreover, recent work showed that m 6 A regulates haematopoietic stem cell regeneration 37 and axon regeneration in the mouse nervous system. 38 While those previous studies have provided important insights into the underlying molecular mechanism of m 6 A-mediated regulation of regeneration at the cellular level, whether and how m 6 A is implicated in regeneration in tissues and entire organisms is unknown.
Here, we employed planarian Schmidtea mediterranea as the model to investigate the role of m 6 A in stem cell function and wholebody regeneration. The results illustrated that depletion of the m 6 A methyltransferase regulatory subunit wtap abolished planarian regeneration, which is mainly mediated by wtap-dependent m 6 20 These timepoints display a unique composition of piwi-1 (smedwi-1) positive population and smedwi-1 negative population, 20 resembling the combination of both processes of tissue regeneration from proliferating neoblasts (epimorphosis) and remodelling of the existing tissues (morphallaxis). 39 In order to gain an overview of the gene expression dynamics from our own perspective during planarian regeneration, we performed RNA-seq at five different time points after amputation (0 hpa, 6 hpa, 3 dpa, 7 dpa, 11 dpa; 0 hpa means sample collected immediately after amputation; Figure 1A). To investigate the gene expression changes of the planarian during regeneration, we grouped all the expressed genes into five clusters according to their distinct expression pattern, with each cluster containing groups of genes upregulated specifically at a certain timepoint after amputation ( Figure 1B). Since m 6 A modification is catalysed by three key components of the methyltransferases, including mettl3, mettl14 and indispensable regulatory component wtap, 40,41 we found that all these three genes displayed upregulated expression at 3, 7 and 11 dpa ( Figure 1C), and all belong to the clusters 3 and 4. Based on the gene ontology (GO) functional analysis, several genes in cluster 4 were observed to be related to the regulation of system process, signalling and cell communication ( Figure 1D, Table S1), including growth factors and transcription factors such as wnt2 and egfr ( Figure 1D). Furthermore, GO enrichment analysis revealed that genes from cluster 4 also involve in nervous system development and differentiation ( Figure 1E, Table S1). Since the components of m 6 A enzymes have a similar expression trend as genes in cluster 4, these findings suggest the potentially important roles of m 6 A in regulating planarian regeneration through the aforementioned pathways.

| m 6 A changes dynamically during planarian regeneration
To determine the regulatory role of m 6 A in planarian regeneration, we first measured the relative level of m 6 A in planarian mRNA, which shows a much higher level of over 0.5% relative to another modification on adenosine (N 1 -methyladenosine, m 1 A) by using UHPLC-MRM-MS/MS (ultra-high-performance liquid chromatography-triple quadrupole mass spectrometry coupled with multiple-reaction monitoring; Figure 2A).
To obtain m 6 A landscape during planarian regeneration, we performed transcriptome-wide m 6 A mapping on samples of five time points (0 hpa, 6 hpa, 3 dpa, 7 dpa, 11 dpa) of regenerative process ( Figure 1A). We carried out methylated RNA immunoprecipitation sequencing (MeRIP-seq) to investigate the features and distribution dynamics of m 6 A in mRNA. In total, we identified 9057-13,362 m 6 A peaks over all five different regeneration time points. We found that most genes have one MeRIP peak ( Figure S1A). The total number of m 6 A peaks at each stage (especially at 6 hpa and 7 dpa) of (A) Percentage of modified A to total A (%)  Figure S1C. (D) Barplot showing the significant gene ontology (GO) terms for 12,016 mRNAs that with m 6 A modification of at least one timepoint during regeneration showing in Figure S1C. (E) Barplot showing the significant GO terms for 1532 mRNAs that with m 6 A modification of four different regeneration timepoints (6 hpa, 3 dpa, 7 dpa and 11 dpa) showing in Figure S1C. (F) Line chart showing one of the trends (fourth category) of mRNAs m 6 A level during regeneration, which with gradual decreased m 6 A level from 0 hpa to 11 dpa. mRNAs with different expression pattern were defined by MEV with parameter-distance-metric-selection = Pearson correlation, number of cluster = 4, maximum-iterations = 50. (G) Barplot showing the significant GO terms for mRNAs shown in Figure 2F. See also Figure S1.
regeneration is slightly higher than that at 0 hpa ( Figure S1B). Consistent with previous observations in mammals, 31 we found that m 6 A peaks were markedly enriched near the stop codon. Also, the distribution pattern of m 6 A is similar among different regeneration time points ( Figure 2B). While transcripts of 6385 genes contain m 6 A peaks at all regeneration time points including the one before amputation (0 hpa), 1532 genes show transcript being modified only after amputation (i.e., only in timepoints other than 0 hpa). We can also find m 6 A-modified transcripts that are unique to each particular time point during regeneration ( Figure S1C, Table S2).
GO enrichment analysis revealed that the 6385 shared m 6 Amodified gene transcripts are associated with various essential biological processes, including regulation of signalling, regulation of cell communication and cell differentiation ( Figure 2C, Table S1). The  F I G U R E 3 Legend on next page. m 6 A-modified genes during regeneration were related to signal transduction, cell differentiation, and nervous system development ( Figure 2D, Table S1), indicating an important role of m 6 A in these regeneration-related processes. The overlap of all m 6 A-modified transcripts with five different clusters indicates that the cluster 4 genes with upregulation at 7 and 11 dpa show the highest overlap mostly with m 6 A-modified gene transcripts ( Figure S1D). GO analysis shows that the genes with their transcripts exclusively modified during the regeneration period are related to response to stimulus, cell communication and signal transduction pathways ( Figure 2E, Table S1).  Table S1). These results may imply a potential role of m 6 A in cell-cell communication and neuron development.

| m 6 A methyltransferase wtap depletion impairs regeneration in planarians
To further study the functional roles of m 6 A modification in planarians, we first performed whole-mount in situ hybridization of m 6 A methyltransferase regulatory subunit, wtap, and found that wtap specifically colocalizes with the neoblast marker smedwi-1, especially at the posterior region, suggesting a regulatory role of m 6 A in neoblast ( Figure 3A,B). To study whether m 6 A modification is indispensable for planarian regeneration, the key regulatory unit, wtap, was knocked down by RNAi, and then the worm was amputated at both anterior and posterior regions ( Figure 3C). Protein expression and mRNA levels of wtap were significantly reduced compared to the control, as shown by whole-mount immunostaining ( Figure 3D), western blot ( Figure S2A) and quantitative reverse-transcription PCR (qRT-PCR) ( Figure S2B). Accordingly, the mRNA m 6 A level was reduced by nearly 80% at 7 dpa ( Figure S2C).
Most importantly, wtap-deficient planarians failed to fully regenerate the amputated tissue parts compared to the controls, which regrow head and tail tissues at 9 dpa ( Figure 3E). Moreover, newly regenerated tissue in the wtap-deficient planarian area was significantly smaller at 3 dpa than that of the control planarians ( Figure S2D). The blastema formed at the wound region of wtapdeficient planarians was also smaller and failed to grow the missing anterior and posterior tissues in comparison with the control planarians ( Figure 3F). Especially at 7 dpa, when control planarians had regrown photoreceptors, the wtap knockdown planarians appeared to be still unable to give rise to any new tissues from their healed wounds ( Figure S2D). We noticed that smedwi-1 RNA expression did not change obviously in wtap knockdown planarians ( Figure 3G). Since These results indicate that m 6 A is essential for planarian regeneration, and its depletion leads to abnormal cell proliferation accompanied by increased apoptosis.

| WTAP-mediated m 6 A controls cell cycle-and cell-cell communication-related factors essential for regeneration
To investigate the underlying molecular mechanism for the defective planarian regeneration mediated by a decreased m 6    25 knockdown planarians ( Figure 1A). Even though the distribution of m 6 A along the CDS only slightly differs in wtap knockdown compared to controls, m 6 A peaks near stop codons become visibly reduced in wtap knockdown compared to the control worms at each time point during regeneration ( Figure 4A). We found that among the genes having upregulation trends in cluster 4 ( Figure 1B), many genes showed a disordered expression pattern after wtap knockdown ( Figure 4B). Further, GO functional analysis indicated that these genes are related to the regulation of signalling, cell communication and signal transduction-associated pathways ( Figure 4B, Table S1). We then analysed the differentially  Table S3).
In order to explore the regulation of m 6 (Table S4). GO analysis revealed that during the regeneration, after wtap knockdown, the upregulated genes in epidermal were enriched in the nucleic acid metabolic process and cell cycle-related functions ( Figure S4B, Table S7), and related to cell cycle regulation functions in neoblast ( Figure S4C). In addition, the   Figure S4D).
The functions of downregulated genes in the neoblast are associated with RNA splicing, DNA replication initiation, protein folding, and so forth ( Figure S4E). Neuronal cells are associated with protein folding, neurotransmitter secretion and signal release from synapse functions ( Figure S4F). These data suggest that wtap knockdown-induced regeneration retardation might through affecting the cell cycle and signal transduction pathways of neoblast and neuronal cells.  Table S5).
To further investigate the functions of NP-like cells in planarian regeneration, we analysed the effect of wtap knockdown on the expression of NP-like cell markers. First, the IGV and enrichment analyses showed that the levels of m 6 A modification on grn and ubiqp transcripts were significantly decreased upon wtap knockdown at both 7 and 11 dpa (Figures 6A and S5C). Consistently, WTAPmediated m 6 A depletion leads to a drastic increase in grn and ubiqp expression in wtap-deficient planarians ( Figure S5B). Furthermore, through both whole-mount fluorescent in situ hybridization and immunofluorescence assays, we observed an increase in grn mRNA and protein levels ( Figure 6B-D), as well as increased colocalization of two markers of NP-like cells, grn and dph1 ( Figure 6B,C). In addition, the expression of dph1 has been restored by the double knockdown of wtap and grn ( Figure S6A). It has been reported that GRN is a secretory protein important for neuronal functions [48][49][50] and can exert function over a distance. 48,51 Furthermore, GRN can stimulate cell division by promote G1 phase to M phase directly. 52 We found that the dou-

| DISCUSSION
Regeneration is a dynamic and tightly controlled process involving a complex yet well-orchestrated gene regulation network. However, the potential significance of RNA m 6 A modification in regeneration is not well understood. Recently, one study reported that kiaa1429 knockdown resulted in failed regeneration and decreased m 6 A level as well as accordingly changed gene expression pattern. 53 In our study,     RNA m 6 A modification has been reported to regulate stem cell differentiation and tissue regeneration. 27,29,[36][37][38]54 It is required for axon regeneration by promoting protein translation of regenerationassociated genes in the rodent model. 38 Meanwhile, deletion of m 6 A reader YTH domain-containing family protein 2 (YTHDF2) in mouse haematopoietic stem cells (HSC) was shown to facilitate HSC regeneration by enhancing the stability of mRNAs related to both Wnt signal and survival. 37 In this study, using planarian Schmidtea mediterranea model system, we found the m 6 A peaks increase at CDS and 3 0 -UTR junction during the regeneration progression. Moreover, knockdown of m 6 A regulatory unit, wtap, resulted in a significant reduction in mRNA m 6 A level and abolished the regenerative capability of planarian to regrow missing head and tail. These findings suggest that wtapmediated m 6 A is imperative for the whole-organism regeneration of planarian.

c o n tr o l ( R N A i) w ta p ( R N A i) w ta p ( R N A i) + g r n ( R N A i) w ta p ( R N A i) + c d k 7 ( R N A i) w ta p ( R N A i) + c d k 9 ( R N A i)
Regeneration involves cell proliferation and the consequent generation of new tissues, which requires tight control of the cell cycle. 55 In our study, we found that cell cycle-related genes such as cdk7, cdk9 and ccnt1 are mainly expressed in neoblasts and upregulated upon wtap knockdown. Both CDK7 and CDK9 are members of the cyclindependent protein kinase family with effector CDK activity to phosphorylate Pol II and other targets within the transcriptional machinery, and also serve as a CDK-activating kinase for at least one other essential CDK involved in transcription. [56][57][58] Both CDK7 and CDK9 have been shown to be required for the cell cycle regulation, such as the blockage of all cell divisions in C. elegans by Cdk7 loss 43 and G 1 cell cycle arrest of D.melanogaster cells by RNAi-mediated Cdk9 silencing. 47 Moreover granulin (grn)/progranulin has been shown to stimulate cell division, 52 and grn inhibition decreases gata1 expression and further the differentiation of erythroid in zebrafish. 59 In support, we found an induced expression of grn, cdk7 and cdk9 post wtap knockdown accompanied by impaired regeneration, suggesting that m 6 A modification are critical for planarian regeneration through destabilizing the transcripts of these three genes.
Importantly, accumulative studies revealed that regeneration is a complicated process involving cell-cell communication in various tissue types. 60 The flow of long-range patterning information in regenerative morphogenesis is crucial to control stem cell behaviour in vivo. 61 For example, muscle regeneration in response to injury, is a nonspecific inflammatory response to trauma, involving interaction between muscle and the immune system. 60  Among all major cell types identified in planarian by scRNAseq, 17,18,20 the most important cell type among them is the neoblasts, which can provide the cellular source for robust planarian regeneration. This totipotent stem cell can differentiate into all other cell types after injury or stimulation. Soon after amputation, neoblast starts to proliferate and migrate to wound site to form regenerative blastemal tissue, which becomes the foundation for re-growing the missing tissue. 5 However, both blastema formation and remodelling of preexisting tissue (morphallaxis) are required for the planarian regeneration. 55 These intricate processes require correct cell-cell communications and precise coordination between stem cells and pre-existing tissue to ensure the regeneration is finished on the right track. 5 In our study, we identified a neural progenitor-like cell cluster accumulated during the regeneration upon wtap knockdown, and further experiments demonstrated that the aberrant high expression of m 6 A modified gene, grn, secreted from this cell cluster disrupts the regeneration process. It is likely that such an abnormal cell cluster is supposed to be a group of neural progenitor cells at a special state during normal neural regeneration, which at the same time is also critical for subsequent events of regeneration of other tissue types to proceed. Therefore, we speculate that this special state of neural progenitor cell cluster, as we amplified its cell number by wtap knockdown, triggers a critical checkpoint regulated during the entire regeneration process that involves cell-cell communication mediated by GRN to ensure the correct path for regeneration, and wtap knockdown induced aberrant grn expression will stop the regeneration at this checkpoint.
As depicted in our model ( Figure S7F), during planarian homeostasis, the expression level of grn is held at a moderate level to inhibit overgrowth, as equivalent to grn expression level before amputation at 0 hpa ( Figure S5B). Upon injury, m 6 A modification selectively targets several important gene transcripts, including grn for degradation, manifested as its reduced expression level after 6 hpa ( Figure S5B it is uncertain how the expression levels of grn, cdk7 or cdk9 regulate the cell cycle and cell-to-cell communication, which eventually lead to regeneration failure. Also, it is worth noting that the wtap knockdown phenotype can produce not only regeneration-specific phenotype, but also lysis phenotype under homeostatic conditions (data not shown).
Therefore, it would be interesting to dissect further which wtapmediated m 6 A genes are regeneration specific. Other model systems which can harness overexpression systems may be able to answer these questions in depth.
Collectively, our study uncovered an essential role of WTAPmediated m 6 A modification in the regeneration of an entire organism.
We further discovered that WTAP-mediated m 6 A modification is particularly important for the neural progenitor cell population. So far, this process has not been well understood. Our study may have uncovered a novel mechanism in regulating planarian regeneration that utilizes GRN to monitor and permit the regeneration to happen correctly. Therefore, our study highlights the potential importance of cell-cell communication during the planarian regeneration and provides important insights for future studies in regenerative medicine.

| Planarian culture
Animals Clonal asexual (CIW4) and sexual strains of Schmidtea mediterranea were maintained in Montjuïch salts as previously described

| Replication, size estimation and randomization
For every assay, at least three independent replicates with a minimum of three animals per experiment were performed. For RNAi phenotype characterization, numbers of animals used are indicated in each panel. No sample size estimation was performed. Animals for all experiments were randomly selected. All animals were included in statistical analyses, and no exclusions were done. Images were randomized before quantifications.

| RNA interference
Double-stranded RNA (dsRNA) was synthesized by in vitro transcription as previously described. 63 In summary, whole worm cDNA was generated and subcloned into a vector. In vitro transcription was performed using T7 polymerase (Promega, P2075). Planarians were fed by dsRNA mixed with pig liver paste three times every 4 days and were amputated into three fragments pre-and post-pharyngeal 24 h after the last feeding. For double knockdown of two genes (double RNAi), dsRNAs of two genes were mixed in liver paste each at the same final concentration as single gene knockdown assay, and delivered three times at a 4-day interval. The C. elegans unc-22 gene was used as a negative RNAi control for single knockdown experiments, and was mixed with wtap dsRNA as control for double knockdown. The knockdown efficiency of double RNAi or single RNAi was verified by RT-qPCR. Primers used for the dsRNA synthesis are listed in Table S6. ImageJ (V1.5) were used to measure the pixel area of regenerated blastema divided by that of the whole tissue fraction, and the same size of the amputated tissue fractions were chosen for quantification.

| BrdU labelling and whole-mount immunofluorescent staining
Animals were fed with 20 mg/mL BrdU (Sigma, 19160) as described before. 6 Specimens were harvested after labelling for 6 h and fixed in 4% formaldehyde (FA). The rat anti-BrdU (Abcam, ab6326) antibody was used for detection. Immunofluorescence staining was performed as described before. 63  The primers used are listed in Table S6. Quantification of fluorescent signals was analysed by ImageJ software (V1.5) as previously described. 64

| Total RNA extraction
Total RNA was extracted from planarians with 1 mL TRIzol ® reagent (Invitrogen, 15596018) through homogenizing on a tissue disruptor. mRNAs were purified from total RNAs using Dynabeads ® mRNA Purification Kit (Ambion, 61006) and subjected to TURBO™ DNase (Invitrogen) treatment at 37 C for 30 min and ethanol precipitation.
After centrifugation and extensive washing with 75% ethanol, the mRNA was dissolved and quantified using Qubit 3.0 (Thermo Fisher).

| Western blot
Planarians were dissolved and homogenized in RIPA buffer (Cell Signalling Technology, 9806s) on a tissue disruptor. After incubation on ice for 10 min, the samples were centrifuged at 13,000 rpm for

| RNA-seq data analysis
Adaptor sequences were trimmed off for all raw reads using the Cutadapt software (version 1.2.1). 69 Reads that were less than 18 nt in length or contained an ambiguous nucleotide were discarded by Trimmomatic (version 0.30). 70 The remaining reads were aligned to the Schmidtea mediterranea transcriptome smed_20140614 using Bowtie2 (v2.2.9) 71 default parameters. Counts were calculated for the sum of reads of each transcript. R package 'DEseq2' 72 was used to identify differentially expressed transcripts. Genes with different expression patterns ( Figure 1B and Figure 4B) were defined by the K-means algorithm in MEV software. In detail, genes were clustered by MEV with parameter-distance-metric-selection = Pearson correlation, number of cluster = N, maximum iterations = 50 and the gene expression level of each time point was compared to that of the previous one. 73

| MeRIP-seq data analysis
For MeRIP-seq, m 6 A-enriched peaks in each m 6 A immunoprecipitation sample were identified by MACS2 peak-calling software (version 2.0.10) 74 with the corresponding input sample serving as control.
MACS2 was run with default options except for '-nomodel' to turn off fragment size estimation. A stringent cut-off threshold for a pvalue of 0.001 was used to obtain high-confidence peaks. In each stage, peaks with over 50% overlap in three biological replicates were used in the downstream analysis. The overlapping peaks were reannotated by the highest peak within the region of all three peaks.
The m 6 A level with different patterns was defined by the K-means algorithm in MEV software (version 4.4). In detail, genes were clustered by MEV with parameter-distance-metric-selection = Pearson correlation, maximum-iterations = 50 and m 6 A level of each timepoint was compared to that of the previous one. 73 4.15 | Single-cell RNA-seq data processing Reads were processed using the Cell Ranger 3.0.0 pipeline 75 with default and recommended parameters. FASTQ files generated from Illumina sequencing were aligned to the Schmidtea mediterranea transcriptome smed_20140614 using the STAR algorithm. 76 Next, Gene-Barcode matrices were generated for each individual sample by counting unique molecular identifiers and filtering non-cell associated barcodes. Only genes that can be translated into proteins were retained. We generated a gene-barcode matrix containing barcoded cells and gene expression counts. This output was then imported into the Seurat (v3.1.2) 77 R toolkit for quality control and downstream analysis of single-cell RNA-seq data. All functions were run with default parameters, unless specified otherwise. Low-quality cells (<500 genes/cell, >6000 genes/cell, <15 cells/ gene and >20% mitochondrial genes) were excluded. In order to exclude multiple captures, which is a major concern in microdroplet-based experiments, Double-tFinder (version 2.0.2) 78 was employed to remove top N cells with the highest pANN score for each library separately, where N represents the doublet rates. Then all the datasets were merged using the 'merge' function in Seurat.

| Identification of cell types and subtypes by nonlinear dimensional reduction
The Seurat package implemented in R was applied to identify major cell types. 77 Highly variable genes were generated and used to perform PCA. Significant principle components were determined using JackStraw 77 analysis and finally focusing on PCs 1-20. We grouped cell types based on their expression profile and matched them to known markers. To identify the subtype of neuronal cells, we used the reported markers of neuronal cells and also using 'FindTransferAnchors' and 'TransferData' function 77 which provided by Seurat to predict the similarity between the subcluster and the reported cell types that identified by Plass. 18 The gene name equivalencies beyond SMED to the transcriptome that used by Plass were transferred using the correspondence table which was downloaded from https://planosphere. stowers.org/chado/analysis?name=&program=&sourcename.

| GO analysis
To enrich the GO annotation of planarian, we aligned the planarian transcriptome with the uniport database through blastx programme, and selected the best comparison result. According to the comparison results and the related GO annotation, the GO annotations of planarians were enriched. Then, GO analysis was performed using topGO, 79 with all genes as background.

| Cell-cell communication analysis
In order to explore cell-cell communication networks via ligandreceptor interactions, we initially identified the homologous genes of planarian and human genes by blastx (version 2.7.1). 80

| Quantification and statistical analysis
All statistical analyses of qPCR and imaging were performed at least three biological replicates. Student's two-tailed unpaired t-test was used for statistical comparisons and data were shown as mean ± SD.
Values of p were used for significance analysis; ns, no significance.
All statistical analyses were performed using GraphPad Prism (version 7.0).

This work was supported by grants from the Strategic Priority Research
Program of the Chinese Academy of Sciences, China (XDA16010501,