JmjC Domain-Encoding Genes Are Conserved in Highly Regenerative Metazoans and Are Associated with Planarian Whole-Body Regeneration.

The capacity for regeneration varies greatly among metazoans, yet little is known about the evolutionary processes leading to such different regeneration abilities. In particular, highly regenerative species such as planarians and cnidarians can regenerate the whole body from an amputated fragment; however, a common molecular basis, if any, among these species remains unclear. Here, we show that genes encoding Jumonji C (JmjC) domain-containing proteins are associated with high regeneration ability. We classified 132 fully sequenced metazoans into two groups with high or low regeneration abilities and identified 118 genes conserved in the high regenerative group that were lost in species in the low regeneration group during evolution. Ninety-six percent of them were JmjC domain-encoding genes. We denoted the candidate genes as high regenerative species-specific JmjC domain-encoding genes (HRJDs). We observed losses of HRJDs in Helobdella robusta, which lost its high regeneration ability during evolution based on phylogenetic analysis. By RNA sequencing analyses, we observed that HRJD orthologs were differentially expressed during regeneration in two Cnidarians, as well as Platyhelminthes and Urochordata, which are highly regenerative species. Furthermore, >50% of the head and tail parts of amputated planarians (Dugesia japonica) died during regeneration after RNA interference of HRJD orthologs. These results indicate that HRJD are strongly associated with a high regeneration ability in metazoans. HRJD paralogs regulate gene expression by histone demethylation; thus, HRJD may be related to epigenetic regulation controlling stem cell renewal and stem cell differentiation during regeneration. We propose that HRJD play a central role in epigenetic regulation during regeneration.

Introduction

Regeneration is a process wherein lost body parts are replaced, and unveiling the complicated biological processes associated with regeneration remains a major challenge in biology. The regeneration abilities of metazoans have been classified into five levels, including regeneration of the whole body, structural elements (e.g. limbs, tail, and fins), internal organs (e.g. heart and liver), tissues (e.g. gut lining), and cells (e.g. axon and muscle fiber regeneration) (Bely and Nyberg 2010). Regeneration is an obvious beneficial trait; however, large differences in regeneration capability are found among metazoans (Bely and Nyberg 2010; Agata and Inoue 2012). For example, planarians and sea anemones can regenerate a new individual from a piece of the body, whereas some animals such as birds and nematodes cannot regenerate any structures (Sánchez Alvarado 2000; Bely and Nyberg 2010; Agata and Inoue 2012).

It is critical to understand the molecular basis leading to the different regeneration abilities among metazoans, and some common genetic factors related to regeneration have already been reported. Wnt/β-catenin, fibroblast growth factor (FGF), and Bone morphogenetic proteins (BMPs) signaling are commonly associated with regeneration in metazoans and are important regulators of cell differentiation or self-renewal during regeneration (Holstein 2012). For example, Wnt genes are expressed in the regenerating heart of mammals (Ozhan and Weidinger 2015), fins of zebrafish (Poss et al. 2000), tail of amphibians (Caubit et al. 1997), and head of hydra (Bode 2003). Furthermore, BMP cooperates with FGF signaling to induce nerve cells to regulate limb regeneration in amphibians (Makanae et al. 2016). BMP genes regulate the symmetric regeneration in the early stage of regeneration in planarians (Reddien et al. 2007). These genes are often shared in several metazoan regeneration processes including structure and whole-body regeneration; however, a shared genetic basis (if one exists) truly conferring a high regeneration ability, such as whole-body regeneration, is still unknown.

New techniques such as next-generation sequencing and high-performance supercomputing enable high-throughput, genome-wide data analysis for large numbers of species. A recent report revealed a common set of differentially expressed genes (DEGs) responding to the regeneration process among distantly related species (Fumagalli et al. 2018). It may be possible to identify conserved genes that are associated with regeneration capacity in highly regenerative metazoans by performing comparative-genome and whole-transcriptome analyses. Basal lineages including Cnidaria, Placozoa, Ctenophora, and Porifera can regenerate their whole body, although the evolution of genetic background conferring the ability has been poorly understood yet (Sánchez Alvarado 2000; Bely and Nyberg 2010; Agata and Inoue 2012; Bely et al. 2014). Based on the phylogenetic distribution of high regenerative species in metazoans, the highly regenerative ability seemed to originate from early animals after they acquired multicellularity as an epiphenomenon of development (Bely and Nyberg 2010). We hypothesized that a common ancestor of all metazoans had a high regeneration ability and that particular genes were lost independently in multiple phyla, resulting in some species in individual phyla showing low regeneration abilities. To test this hypothesis, we defined species that regenerate whole, anterior, or posterior body parts as highly regenerative species, and other species that regenerate only appendages such as limbs, tails, fins, or lower structures as low regenerative species. According to this definition, we classified metazoans with fully sequenced genomes as highly regenerative species (Urochordata, Cephalochordata, Echinodermata, Hemichordata, Annelida, Platyhelminthes, Cnidaria, Placozoa, Ctenophora, and Porifera) or low regenerative species (Vertebrata, Arthropoda, Mollusca, and Nematoda) (supplementary table S1, Supplementary Material online). The purpose of this study was to identify genes that are differentially expressed during regeneration in four highly regenerative metazoans (Schmidtea mediterranea, Nematostella vectensis, Hydra vulgaris, and Ciona intestinalis) whose orthologs were conserved only in the genomes of species in the highly regenerative group during evolution. We provide evidence that the identified candidate genes are associated with regeneration in planarians by RNA interference (RNAi) knockdown experiments (Rouhana et al. 2013).

Materials and Methods

Identification of Genes Specific to Highly Regenerative Species

We downloaded the complete protein sequences of 126 metazoan species with fully sequenced genomes deposited in Ensembl (release 87) and Ensembl Metazoa (release 33). We also downloaded the protein sequences for highly regenerative species, Acanthaster planci v1.0 (Hall et al. 2017) and Ptychodera flava v3.0 (Simakov et al. 2015), from OIST (http://marinegenomics.oist.jp), H. vulgaris from NCBI (Chapman et al. 2010), Branchiostoma floridae v1.0 (Putnam et al. 2008) from JGI Genome Portal (https://genome.jgi.doe.gov), S. mediterranea (Robb et al. 2015) from SmedGD (http://smedgd.neuro.utah.edu), and Exaiptasia pallida v1.0 (Baumgarten et al. 2015) from REEFGENOMICS (http://aiptasia.reefgenomics.org). To identify genes conserved only among highly regenerative species, we classified the species into two groups (low and highly regenerative species), based on their regenerative abilities according to the work of Bely and Nyberg (2010) and related papers (supplementary table S1, Supplementary Material online). We did not use species for which a high regeneration ability has not been documented (Bely and Nyberg 2010), even though their genomic sequences were available. In most phyla, the representative regenerative ability of a phylum reported by Bely and Nyberg (2010) is consistent with the regenerative abilities of species in the phylum, but the degree of the regeneration ability does not match between C.teleta (high regenerative species) and H.robusta (low regenerative species) in Annelida. We found that 116 species belonging to 4 phyla and 14 species belonging to 10 phyla were in the low and highly regenerative groups, respectively (supplementary table S1, Supplementary Material online). Cephalochordata was not classified in the highly regenerative group by Bely and Nyberg (2010), although a recent study reported the extensive regeneration ability of Branchiostoma lanceolatum (Somorjai et al. 2012); therefore, we classified Cephalochordata into the highly regenerative group. We conducted BlastP searches using protein sequences of all potential highly regenerative species as queries against protein sequences from each species used in this study to obtain homologs of potential highly regenerative species (threshold: BlastP score >100). Genes conserved among all species in the highly regenerative group that were lost in all species of the low regenerative group were regarded as candidate genes possibly associated with highly regeneration ability. We also identified genes for which a phylum of the low regenerative group had an ortholog as candidate genes (table 1).

Table 1

Number of Genes Specifically Conserved in High Regenerative Metazoans

Phyla	Species	Regeneration Group	Number of Genes Conserved in High Regenerative Species
			With JmjC Domain	Without JmjC Domain
Porifera	A. queenslandica	High	1	0
Placozoa	Trichoplax adhaerens	High	16	1
Ctenophora	Mnemiopsis leidyi	High	3	0
Cnidaria	N. vectensis	High	9	1
Cnidaria	H. vulgaris	High	2	0
Cnidaria	E. pallida	High	8	0
Annelida	H. robusta	Lowa	1	0
Annelida	C. teleta	High	3	0
Mollusca	Lottia gigantea	Lowb	2	0
Mollusca	Octopus bimaculoides	Lowb	1	0
Mollusca	Crassostrea gigas	Lowb	3	1
Platyhelminthes	S. mediterranea	High	2	0
Echinodermata	A. planci	High	5	0
Echinodermata	Strongylocentrotus purpuratus	High	3	0
Hemichordata	P. flava	High	7	1
Cephalochordata	B. floridae	High	21	0
Urochordata	Ciona savignyi	High	18	0
Urochordata	C. intestinalis	High	13	1
Total			118	5

Helobdella robusta is a low regenerative species having a JmjC domain-encoding gene.

Mollusca having JmjC domain-encoding genes is in the low regenerative group.

RNA Sequencing Analyses of the Regeneration Process

We used published RNA sequencing (RNA-Seq) data sets for N. vectensis (Schaffer et al. 2016), H. vulgaris, C. intestinalis (Spina et al. 2017), and S. mediterranea (Zeng et al. 2018) for our analyses. During the regeneration process in N. vectensis, individuals were divided into two parts, and total RNA from oral and physa blastemas was extracted at 0, 8, 24, and 72 h post-amputation (hpa) (Schaffer et al. 2016). Total H. vulgaris blastema RNA was extracted from the head blastemas at 0, 2, 4, 6, 12, and 48 hpa. Cionaintestinalis RNA samples were also extracted from tissues of non-regenerated normal oral siphons immediately after amputation as controls and from tissues at different times during oral regeneration (1, 3, and 8 days post-amputation [dpa]). The RNA samples of S. mediterranea were extracted from regenerating individuals at various times (0, 1, 3, 6, 9, 12, 24, 36, 48, 60, 72 hpa and 4, 5, 6, 7, 10, 14 dpa) (Zeng et al. 2018). The single-end libraries from H. vulgaris, N. vectensis, C. intestinalis, and S. mediterranea samples were sequenced using the Illumina HiSeq 2500, HiSeq 2000, and Ion Torrent Proton platforms, respectively. We downloaded RNA-Seq data sets derived from both normal and regenerative conditions for three highly regenerative species from the Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra), which were previously deposited under accession numbers PRJNA270225 (H. vulgaris), PRJNA330595 (N. vectensis), and PRJNA421768 (S. mediterranea). We did not use raw C. intestinalis reads from the SRA because the number of mapped short reads from C. intestinalis could be downloaded from the Gene Expression Omnibus website (https://www.ncbi.nlm.nih.gov/geo). The numbers of RNA-Seq sample replicates for each stage were two, two, four, and three for H. vulgaris, N. vectensis, S. mediterranea, and C. intestinalis, respectively.

We checked the quality of the downloaded RNA-Seq reads for H. vulgaris and N. vectensis by FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Contaminating Illumina adapters in the reads and low-quality reads (quality value <30) were excluded using the FASTX_Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Then, we aligned short RNA-Seq reads for H. vulgaris, N. vectensis, and S. mediterranea to their reference genomes from NCBI, Ensembl Metazoa (http://metazoa.ensembl.org), and SmedGD (http://smedgd.neuro.utah.edu), respectively, using Tophat2, version 2.1.1 (Kim et al. 2013). Next, we estimated the read counts using HTSeq, version 0.6.1 (Anders et al. 2015). Read counts for C. intestinalis were estimated according to our previously published procedure (Spina et al. 2017). Finally, we identified DEGs under normal and regenerative conditions in these three species using the “TCC” package, version 1.20.0 (Sun et al. 2013) in R.

Identification of Orthologous Gene Groups for Four Highly Regenerative Species

We conducted an all-against-all BlastP search (BlastP score >100) with the protein sequences for N. vectensis, H. vulgaris, S. mediterranea, and C. intestinalis and found the best hit (highest BlastP score) for each species. Genes sharing the same best hit were clustered for the four highly regenerative species. We considered clustered genes as genes in the same orthologous gene group. We determined the number of orthologous gene groups depending on the presence of genes for all four species in an orthologous gene group.

Phylogenetic Relationships of Jumonji C Domain-Encoding Genes in Metazoans

To understand the phylogenetic relationships of high regenerative species-specific Jumonji C (JmjC) domain-encoding genes (HRJDs) in metazoans, we analyzed all HRJDs in the highly regenerative species shown in table 1 against outgroup genes. The inferred HRJD orthologs identified by BlastP were included in the highly regenerative species and the Mollusca species. When searching for outgroup genes, we found a human gene (ENSG00000155666, KDM8) that was the best hit with respect to the protein sequences of our candidate genes in N. vectensis by BlastP searching (BlastP score >100). Then, we collected orthologs found in highly regenerative species and in several model organisms (human, mouse, chicken, zebrafish, fruit fly, and worm) that were the best hits with respect to the human KDM8 protein sequence by BlastP searching. We aligned the protein sequences using MAFFT v7.310 (Katoh et al. 2002) and generated a maximum-likelihood phylogenetic tree using RAxML, v8.2.10 (Stamatakis 2014).

Whole-Mount In Situ Hybridization

The planarians (Dugesiajaponica, strain SSP) were cultured in an incubator at 22–24 °C. The planarians used in this study (5–7 mm in length) were starved for approximately 1 week before conducting the experiments. We conducted whole-mount in situ hybridization for two homologs (DjHRJDa and DjHRJDb) of our candidate genes in the planarian D. japonica, which is one of most well-characterized highly regenerative species (table 1). Amplified DNA fragments of DjpiwiA (positive control) were inserted into plasmids according to methods described in a previous study (Shibata et al. 2016). DjHRJDa and DjHRJDb were cloned into the pCR II-TOPO and pBluescript SK (+/–) vectors, respectively. The DjpiwiA plasmid was linearized by digestion with SamI, and the DjHRJDa and DjHRJDb plasmids were linearized by digestion with BamHI. T7 RNA polymerase was used to prepare digoxigenin-labelled antisense RNA probes using linearized plasmids as the templates. The expression signals were detected using antisense RNA probes. The method followed was described in detail previously (Shibata et al. 2016). The sequences of primers used are shown in supplementary table S4, Supplementary Material online.

Synthesis of DjHRJDa, DjHRJDb, EGFP, and DjpiwiB dsRNA

We prepared dsRNA as described previously (Alvarado and Newmark 1999; Rouhana et al. 2013). EGFP (negative control) and DjpiwiB (positive control) (Hayashi et al. 2010; Shibata et al. 2016) were inserted into the plasmid pBluescript SK, DjHRJDa was cloned into vector pCR II-TOPO, and DjHRJDb was cloned into pBluescript SK (+/–). Each gene was amplified using Ex Taq polymerase to add a T7 promoter site on both ends of the target sequences. The PCR products were gel-purified using the Gel/PCR Extraction Kit (Genetics) and used as a template for synthesizing dsRNA. The primer and dsRNA sequences are shown in supplementary table S4, Supplementary Material online.

RNAi Knockdown by the Feeding Method

Twenty-five microliters of chicken liver solution, 6 μl of 2% agarose, and 6 μl of dsRNA (32 μg/μl) were mixed and fed to a group of 15 planarians. This food mixture was prepared in small aliquots (∼6 μl each) and frozen at −30 °C for 30 min. Then, we conducted three successive feedings at 2-d intervals. siRNA against EGFP was used as a negative control and DjpiwiB siRNA was used as a positive control.

RNAi Knockdown by Microinjection

Planarians (D. japonica, strain SSP) were cultured in an incubator at 22–24 °C. The individuals used in this study (5–7 mm in length) were starved for approximately 1 week before microinjection. The dsRNA was injected 3–5 times into planarians (32 nl/injection) with two successive treatments given 2 days apart using a Drummond Scientific Nanoject injector (Broomall, PA, USA). Control planarians were injected with EGFP dsRNA, which did not target any mRNA encoded in the planarian genome. At 6 h following the last set of injections, the planarians were divided into three parts (the head, trunk, and tail) along the anteroposterior axis using sterile surgical blades. We incubated the fragments at 24 °C in the dark for regeneration. After 2 weeks of regeneration. We counted the surviving fragments in each group to calculate the survival rate.

Quantitative Real-Time PCR Analysis to Determine the RNAi Efficiency

We examined the relative expression levels of target genes in four groups (EGFP RNAi, DjHRJDa RNAi, DjHRJDb RNAi, and double knockdown [DKD]) by real-time PCR. We extracted total RNA from the whole bodies of 15 individual D. japonica samples in each group, and cDNA was synthesized using a QuantiTect Transcription Kit (Qiagen). We carried out quantitative analysis as previously described in Nishimura et al. (2012). The sequences of primers used for all experiments are shown in supplementary table S2, Supplementary Material online.

Results

Regeneration-Related Genes Conserved in Highly Regenerative Species

To identify genes conserved only in the genomes of species with a high regeneration ability, we conducted a BlastP search against a database of protein sequences from 132 metazoans (supplementary table S1, Supplementary Material online) using the protein sequences of 14 species from 10 phyla with the potential for high regeneration as queries. No genes were specifically found in the genomes of species in the highly regenerative group (BlastP score >100; table 1). Then, we used a less stringent condition in which the candidate genes included not only those in highly regenerative phyla, but also those in a phylum with low regeneration, to screen for genes that are largely specific to highly regenerative species. In this manner, we found 123 genes conserved in all species of the highly regenerative phyla and one phylum of species with low regeneration (Mollusca) (table 1). As discussed below, we found the ortholog from the low regenerative species Helobdella robusta in the highly regenerative phylum Annelida (Bely and Nyberg 2010).

JmjC Domain-Encoding Genes Conserved in Highly Regenerative Species

To investigate the protein domains found in genes conserved in highly regenerative species, we analyzed their protein sequences using SMART (http://smart.embl-heidelberg.de/). Remarkably, 96% of the conserved genes encoded JmjC domain-containing proteins (table 1). Only two (EDO46114 and ENSCINT00000030781) of the genes conserved in highly regenerative species did not have any known protein domains. We designated the JmjC domain-encoding genes as highly regenerative species-specific JmjC domain-encoding genes (HRJDs). In general, JmjC domain-containing proteins are related to epigenetic factors and belong to a large gene family (Klose et al. 2006). Although many JmjC domain-containing proteins have other domains such as an FBOX domain or PHD domain (Klose et al. 2006), HRJDs are in a group that contains only a JmjC domain. It has been reported that JmjC domain-containing proteins play an essential role in regulating stem cell renewal (Loh et al. 2007; Shen et al. 2009; Xiao et al. 2017). JmjC domain-containing proteins have conserved residues within the predicted cofactor-binding sites (Klose et al. 2006; Xiao et al. 2017). Three amino acid residues bind to the Fe(II) cofactor, and two additional residues bind to αKG within the JmjC domain. There are amino acid variations in the conserved residues representing active or abrogated enzymatic activities. Each of the identified HRJDs has at least one ortholog with active amino acid residues in highly regenerative species, except for HrHRJD1 in H. robusta, which only has a different amino acid variation in the conserved residue, which may abolish protein function (fig. 1). We examined whether the substituted amino acid in HrHRJD1 was deleterious by Provean (Choi et al. 2015) (http://provean.jcvi.org/). The Provean score was −6.350, which was deleterious. The result indicates that HrHRJD1 in H. robusta may have a divergent function compared with other HRJD orthologs.

. 1.

—Amino acid alignments in highly conserved sites of the JmjC domain for high regeneration species and H. robusta. Gray and black arrows/boxes represent amino acid residues associated with Fe(II)-binding and αKG-binding sites, respectively. The amino acid of HrHRJD1 marked by a solid white character in a green box is different compared with other highly regenerative species.

Phylogenetic Relationship of Genes in the HRJD Family

To examine the phylogenetic relationships of HRJDs in the fully sequenced genomes of metazoans, we analyzed all HRJDs shown in table 1. We also used the sequence for the human KDM8 protein, which is the closest human paralog for HRJDs, as well as KDM8 orthologs in highly regenerative species and several model species (chicken, mouse, zebrafish, fruit fly, Caenorhabditis elegans, and Xenopus tropicalis) as outgroup genes. The phylogenetic tree showed that HRJD orthologs in the highly regenerative species clustered together and were completely separated from the KDM8 outgroup genes, as supported by high bootstrap values (100; fig. 2). Note that we did not find any HRJD orthologs in low regenerative species except for the leech H. robusta and the Mollusca. Helobdellarobusta had only one HRJDsHelroP165856 (HrHRJD1), and another Annelida (Capitella teleta) had three HRJDs (CtHRJD1, CtHRJD2, and CtHRJD3), as shown in figure 2. These results suggest that H. robusta lost HRJDs during evolution. We propose that low regenerative species lost HRJDs independently during metazoan evolution. The tree topology of HRJDs does not follow that of the species tree (fig. 2) because multiple duplication (and/or loss) events of HRJDs are ancient and their protein sequences are quite divergent. Thus, it is hard to detect duplication or loss events in HRJDs by comparing the gene tree to the species tree.

. 2.

—The phylogenetic relationships of metazoans in this study. (A) The phylogenetic relationships among HRJDs in highly regenerative species. Phylogenetic tree of HRJDs in the highly regenerative group aligned using MAFFT software and built using the maximum-likelihood algorithm of RAxML. The gene IDs corresponding to the names on the tree are shown in supplementary table S5, Supplementary Material online. Branch lines shown in red and black indicate HRJD and KDM8 orthologs, respectively. The candidate genes of highly regenerative species were phylogenetically separated from the outgroup genes (KDM8 orthologs). The solid and hollow arrows represent HRJD orthologs for C. teleta and H. robusta, respectively. This figure was made by FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). (B) The phylogenetic tree of highly regenerative groups and Mollusca conserved HRJDs. Abbreviations of species names (Ap, Acanthaster planci; Aq, Amphimedon queenslandica; Bf, Branchiostoma floridae; Ce, Caenorhabditis elegans; Ct, Capitella teleta; Ci, Ciona intestinalis; Cs, Ciona savignyi; Cg, Crassostrea gigas; Dr, Danio rerio; Dm, Drosophila melanogaster; Ep, Exaiptasia pallida; Gg, Gallus gallus; Hr, Helobdella robusta; Hs, Homo sapiens; Hv, Hydra vulgaris; Lg, Lottia gigantea; Ml, Mnemiopsis leidyi; Mm, Mus musculus; Nv, Nematostella vectensis; Ob, Octopus bimaculoides; Pf, Ptychodera flava; Sm, Schmidtea mediterranea; Sp, Strongylocentrotus purpuratus; Ta, Trichoplax adhaerens; Xt, Xenopus tropicalis).

Changes in Expression Levels for HRJDs during Regeneration

To identify genes expressed in response to the regeneration process, we used RNA-Seq data sets from four highly regenerative metazoans (S. mediterranea, N. vectensis, H. vulgaris, and C. intestinalis). Before analyzing the RNA-Seq data sets, we identified 3,070 orthologous gene groups shared by the four high regenerative species (see Materials and Methods section). RNA was extracted from oral and physa blastemas at 0, 8, 24, or 72 h post-amputation (hpa) in N. vectensis (Schaffer et al. 2016); from the head blastemas of H. vulgaris (0, 2, 4, 6, 12, and 48 hpa); and from the blastemas of C. intestinalis oral siphons (Spina et al. 2017) (0, 1, 3, and 8 dpa). RNA samples from S. mediterranea were extracted from regenerating planarians at 0, 1, 3, 6, 9, 12, 24, 36, 48, 60, and 72 hpa, and 4, 5, 6, 7, 10, and 14 dpa (Zeng et al. 2018). We identified DEGs in N. vectensis, H. vulgaris, S. mediterranea, and C. intestinalis, respectively, by comparing their gene expression levels with those observed under normal conditions (supplementary table S2, Supplementary Material online).

Although not all HRJDs responded to the regeneration process, we found that four of nine, one of two, one of two, and four of 13 HRJDs overlapped with DEGs specific to N. vectensis, H. vulgaris, S. mediterranea, and C. intestinalis, respectively. The expression levels of 10 DEGs are shown in supplementary figure S1, Supplementary Material online. To test the enrichment of DEGs in HRJDs, we identified 238 orthologous gene groups (containing 6,304 genes) in which all four species had at least one DEG, out of 3,070 orthologous gene groups (containing 39,510 genes) shared by the four highly regenerative species (supplementary fig. S2, Supplementary Material online). DEGs were significantly enriched in genes of the HRJD orthologous group (10/26 vs. 6,304/39,510; P = 0.0042, χ2 test) during the regeneration of the four highly regenerative species.

DjHRJDa and DjHRJDb Expression in Neoblasts, and Differentiated Cells of D. japonica

The planarian is a typical highly regenerative species and a model species for examining the function of regeneration-related genes (Nishimura et al. 2012; An et al. 2018). Dugesia japonica and S. mediterranea each have two HRJD orthologs, which are DjHRJDa and DjHRJDb in D. japonica and SmHRJDa and SmHRJDb in S. mediterranea (fig. 2 and supplementary table S3, Supplementary Material online). Based on the RNA-Seq analysis in S. mediterranea, we found that SmHRJDb was a DEG during regeneration, whereas SmHRJDa was not (fig. 3). The expression level of SmHRJDb was negatively regulated at the beginning of regeneration, but that of SmHRJDb was upregulated after 6 days of regeneration (fig. 3). We observed that the expression level of smedwi-1, smedwi-2, and smedwi-3 continuously increased after amputation, but they started to decrease after 72 h (fig. 3). Interestingly, the expression level of SmHRJDb negatively correlated with that of smedwi-1 (Pearson product-moment correlation, r = −0.87, P = 0.0047) and smedwi-3 (Pearson product-moment correlation, r = −0.79, P = 0.012) after 72 h (fig. 3). This result indicates that there is direct or indirect interaction between piwi and HRJDs during the regeneration process. Amino acid sequences similarity between DjHRJDb and SmHRJDb, particularly those encoding the JmjC domain, was higher than that between DjHRJDa and SmHRJDa (fig. 3).

. 3.

—Gene expression of HRJD orthologs in planarians. (A) and (B) Expression changes of two HRJDs (SmHRJDa and SmHRJDb) and three piwi genes (smedwi-1, smedwi-2, and smedwi-3) during regeneration in S. mediterranea, respectively. The x axis and y axis indicate the time after amputation and FPKM (fragments per kilobase of exon per million mapped reads) values in S. mediterranea, respectively. The error bars indicate the standard error. The pink area represents the different expression pattern between SmHRJDb and semdwi genes, SmHRJDa was not a DEG during regeneration under any condition. (C) Schematic structures and the identity of amino acid sequences between HRJD orthologs in two planarians. The percentages represent the identity of amino acid sequences for each segmental region of a HRJD protein in two planarians. (D–F) The expression patterns of DjpiwiA (control), DjHRJDa, and DjHRJDb in nonirradiated planarians by in situ hybridization. (G–I) The expression patterns of DjpiwiA (control), DjHRJDa, and DjHRJDb in X-ray irradiated planarians by in situ hybridization.

We speculated that the regeneration-related genes might be expressed specifically in neoblasts which are planarian adult pluripotent stem cells, because neoblasts play a central role in regeneration (Agata and Watanabe 1999; Newmark and Alvarado 2001; Agata et al. 2006). Therefore, we examined the expression patterns of DjHRJDa and DjHRJDb in D. japonica by in situ hybridization (fig. 3). The positive control DjpiwiA, which is a reliable marker of neoblasts (Hayashi et al. 2010; Shibata et al. 2016), was expressed only in neoblasts, and thus we observed that DjpiwiA expression disappeared after specifically killing neoblasts with X-ray irradiation (fig. 3). In contrast, we observed high DjHRJDa and DjHRJDb expression in the whole body even after X-ray irradiation (fig. 3), although the signals for DjHRJDa and DjHRJDb expression without irradiation were higher than those in irradiated planarians. These results indicate that DjHRJDa and DjHRJDb were highly expressed in both neoblasts and differentiated cells of D. japonica.

Regenerative Failure after RNAi-Mediated Knockdown of Genes Encoding the JmjC-Domain Proteins DjHRJDa and DjHRJDb in D. japonica

To investigate whether HRJDs play an important role in high regeneration ability, we conducted RNAi-mediated knockdown of DjHRJDa and DjHRJDb in D. japonica. We used knockdown of DjpiwiB (a reliable marker of neoblasts) (Hayashi et al. 2010; Shibata et al. 2016) and EGFP as RNAi positive and negative controls, respectively. It took approximately 5 days for planarians to digest the double-stranded RNA (dsRNA) after last feeding. We divided eight of 15 individuals into three pieces (the head, trunk, and tail) for each knockdown experiment, and the remaining seven individuals were in the intact group. Note that we did not distinguish among the amputated body parts (the head, trunk, and tail) in the experiment using the dsRNA feeding method. We did not observe abnormal phenotypes in intact individuals treated with dsRNA of DjHRJDa, DjHRJDb, or EGFP (negative control), even after 7 days. This result indicates that the HRJDs were not essential for survival. In contrast, following DjpiwiB knockdown (positive control), 85.7% of intact individuals and 87.5% of individuals with amputation died (fig. 4). The effect of DjpiwiB knockdown was lethal to the planarians, whereas RNAi against DjHRJDa, DjHRJDb, and EGFP did not show any abnormal phenotypes after the first amputation. Then, we conducted a second amputation in which the target was the blastema of amputated individuals (except for individuals with DjpiwiB knockdown). As a result, 28.6% and 46.4% of DjHRJDa and DjHRJDb RNAi knockdown planarians died (including one planarian showing an abnormal phenotype) during the 7-day regeneration process (fig. 4).

. 4.

—Survival rate in knockdown planarians after amputation. (A) The survival rate of amputated pieces during the planarian regeneration process by feeding dsRNA are shown. Individual planarians started to die after the second amputation (n = 28). (B–D) Pictures of dead head fragments after knocking down the expression of DjHRJDa, DjHRJDb, or DjpiwiB (positive control). (E) The survival rates of amputated pieces (head, trunk, and tail) during the planarian regeneration process by injecting dsRNA are shown. Individual planarians in the DjHRJDa RNAi (n = 14), DjHRJDb RNAi (n = 16), and DKD RNAi (n = 17) groups started to die at 1 week post-amputation, and we estimated their survival rates after 2 weeks. Asterisks represent a significant difference in survival rates between knockdown and control individuals (Fisher's exact test, P value <0.05).

Blastemas are important cells for regeneration in planarians (Tasaki et al. 2011). The distribution of blastemas is not uniform in the planarian body. In particular, neoblasts are enriched in the trunk compared with the head and tail (Orii et al. 2005). Thus, we examined the survival rate for each amputated part after RNAi knockdown. To observe more severe phenotypes, we conducted RNAi knockdown by injecting instead of feeding dsRNA. It took 6 h for the planarians to digest the dsRNA after the last microinjection, after which time we amputated all individuals into three pieces (e.g. the head, trunk, and tail) for each knockdown group. To estimate the efficiency of RNAi knockdown, we measured the relative expression levels of DjHRJDa and DjHRJDb after knockdown by quantitative PCR. We observed that the relative expression levels of DjHRJDa and DjHRJDb were reduced by 46.6% and 82.1%, respectively, when each gene was silenced individually (supplementary fig. S3, Supplementary Material online). In addition, DKD experiments (where both genes were targeted by RNAi) reduced the relative expression levels of DjHRJDa and DjHRJDb by 33.5% and 80.5%, respectively (supplementary fig. S3, Supplementary Material online). The knockdown planarians started to die 7 days after the amputation. Death following knockdown likely paralleled the observed expression increase in SmHRJDb (a DjHRJDb ortholog) in S. mediterranea (figs. 3).

The abnormal phenotypes observed in knockdown individuals occurred faster using the microinjection method than with the feeding method. We estimated the survival rates of each group at 14 days after treatment. Almost all amputated pieces in the negative control group (EGFP knockdown) survived (only one tail fragment died). Although the individual trunk fragments in knockdown individuals did not die, the amputated head and tail fragments died in the DjHRJDa (RNAi), DjHRJDb (RNAi), and the DKD groups (fig. 4). Specifically, 50% of the head fragments and 58.8% of the tail fragments died after DjHRJDa and DKD RNAi knockdown, respectively (fig. 4). This finding indicates that regeneration of the head and tail fragments, which include fewer neoblasts compared with the trunk, is likely to be influenced by RNAi knockdown of DjHRJDa and DjHRJDb. The death rate of DKD (DjHRJDa and DjHRJDb) individuals was higher than that for the groups where only DjHRJDa or DjHRJDb was knocked down. The expression level of the DjHRJDb ortholog in S. mediterranea significantly increased during the regeneration process, whereas that of the DjHRJDa ortholog did not change (fig. 3). This result indicates that DjHRJDb expression is needed for regeneration.

Discussion

HRJDs Associated with the Regeneration Process

Common ancestors of metazoans likely had a high regeneration ability, and many species appear to have lost that ability during evolution (Bely and Nyberg 2010). We found that HRJDs were conserved specifically in highly regenerative species. The genes were lost in Vertebrata, Nematoda, and Arthropoda during evolution, although they were conserved in Mollusca. Bely and her colleagues pointed out that regeneration data for Mollusca are limited when compared with other phyla (Bely 2014). Although Mollusca species cannot regenerate many body parts, some of them can regenerate the head (Bely 2014). In addition, it has been reported that Echinoderm larvae can regenerate the whole body, whereas Echinoderm adults cannot (Vickery et al. 2001; Reinardy et al. 2015). Mollusca possessing HRJDs may have also a high potential for regeneration at some spatiotemporal stages (e.g. the larva stage). Further studies are needed to understand the regeneration capability of Mollusca. HRJDs may have multiple functions, with one function related to regeneration ability potentially being lost in Mollusca during evolution (fig. 5).

. 5.

—Overview of losses of high regeneration ability and HRJDs in metazoans during evolution. Phylogenetic relationships of metazoan phyla used in this study. Red and blue characters indicate highly and low regenerative phyla (or species), respectively. Black crosses and gray lines indicate losses of HRJDs and regeneration ability during evolution, respectively.

Members of the Annelida phylum are classified as highly regenerative species (Bely and Nyberg 2010; Bely et al. 2014), although some species (such as the leech) lost the ability to regenerate during evolution (Bely et al. 2014). We surveyed HRJDs of C. teleta and H. robusta in the Annelida phylum and identified three genes conserved in C. teleta, but only one gene conserved in H. robusta (table 1). Leeches lost two HRJDs and the regeneration ability during evolution (fig. 2). Highly regenerative species have multiple HRJDs, whereas H. robusta and the primitive multicellular aquatic metazoan Amphimedon queenslandica have only one HRJD (table 1). These findings suggest that the maintenance of multiple HRJDs may be associated with an increased capacity for regeneration.

Duplication of Ancestral Genes for HRJD and KDM8

There are seven groups (JHMD1, PHF2, JARID, JHDM3, UTX, JHDM2, and JmjC domain only) in the large family of JmjC domain-encoding genes (Klose et al. 2006). HRJDs and their closest paralogs KDM8 are in the JmjC domain only group (fig. 2). It has been reported that KDM8 regulates pluripotency in human embryonic stem cells (Zhu et al. 2014). Amino acid sequences of HRJDs apparently diverged from those of KDM8 and thus they are distinguishable (fig. 2). Although plants and fungi do not have HRJD orthologs, they have KDM8 orthologs (Jones et al. 2010; Gacek-Matthews et al. 2015), metazoans have both HRJD and KDM8 orthologs (fig. 2). Thus, duplication of KDM8 and HRJD would have occurred before the origin of metazoans and after metazoan-fungi split.

Functions of JmjC Domain-Containing Proteins

The JMJD1C protein contains a JmjC domain that is a paralog of HRDJ genes and regulates self-renewal for embryonic stem cells in mice (Peng et al. 2009). Thus, HRJDs may be involved in the regulation of stem cells during regeneration. We propose the following model regarding the molecular functions of HRJDs, based on the known functions of their paralogs that are conserved in metazoans. Members of the family of JmjC domain-containing proteins are epigenetic factors that remove methylation markers on histones to regulate target gene expression and that are part of the epigenetic memory system regulating cell fate and identity (Klose et al. 2006; Tsukada et al. 2006). It has been reported that lysine (K) and arginine (R) residues of histones can be methylated by paralogs of JmjC domain-containing proteins and that histone methylation correlates with transcriptional activation or repression based on the position of K, such as in H3K4 (Eissenberg and Shilatifard 2010), H3K9 (Park et al. 2011), H3K27 (Agger et al. 2007; Hamada et al. 2015; Wiles and Selker 2017), H3K36 (Hsia et al. 2010; Suzuki et al. 2016), or based on the state of demethylation, as with mono- (me1), di- (me2), or tri-demethylation (me3) (Lachner et al. 2003; Martin and Zhang 2005; Chen et al. 2011).

JmjC domain-containing proteins also serve essential roles in regulating stem cell renewal. JMJD1C protein containing a JmjC domain in low regenerative vertebrate species, which is one of HRJD paralogs, demethylates H3K9 and thereby control ERK/MAPK signaling and the epithelial-to-mesenchymal transition, as well as regulate mouse embryonic stem cell self-renewal (Xiao et al. 2017). Other mouse paralogs (e.g. the Jmjd1a and Jmjd2c proteins) functioned as H3K9Me2 and H3K9Me3 demethylases to positively control stem cell renewal, and their depletion lead to stem cell differentiation (Loh et al. 2007). The mouse paralog JmjC domain-containing protein UTX, which is an H3K27 trimethylase, regulates pluripotency, as somatic cells lacking UTX fail to reprogram back to pluripotency induced by transcription factors (Oct4, Sox2, Klf4, and Myc) (Mansour et al. 2012). The human paralogs UTX and JMJD3 are H3K27 trimethylases that decrease the expression level of H3K27me3 associated with HOX genes and delocalized polycomb proteins in vivo during differentiation (Agger et al. 2007). Moreover, the cricket paralog UTX protein promotes joint formation through histone H3K27 methylation to regulate leg regeneration (Agger et al. 2007). We propose that HRJDs influence cellular processes through epigenetic regulation.

Although an ortholog of HRJD (Mina53-like1) in C. intestinalis (table 1) was proposed to promote mesenchymal cell differentiation during embryonic development (Tokuoka et al. 2004), a possible association of this gene with regeneration in C. intestinalis has not been examined. The expression level of Mina53-like1 in C. intestinalis was upregulated by embryonic overexpression of the FGF genes FGF9, FGF16, and FGF20 during development (Tokuoka et al. 2004). FGF genes are important for regeneration processes in many species, such as planarians, zebrafish, and Xenopus laevis (Lee et al. 2005; Fukui and Henry 2011; Shibata et al. 2016). Furthermore, previous data indicated that many genes expressed during embryogenesis are also expressed during regeneration (Christen et al. 2010; Hayashi et al. 2010). Therefore, we speculate that Mina53-like1 in C. intestinalis also plays a role in regeneration.

HRJDs Associated with the High Regeneration Ability of Planarians

We identified two HRJD orthologs (DjHRJDa and DjHRJDb) in planarians and examined their biological functions by RNAi. We observed that many individuals with knocked down DjHRJDa and DjHRJDb expression did not regenerate their lost parts after 2 weeks. Reddien (2013) proposed two models for cell specification in planarian regeneration, namely the naïve-neoblast and specialized-neoblast models. Note that an intermediate process of differentiation achieved by undifferentiated blastema cells was included in both models. The regenerative process in planarians involves two main processes: neoblast proliferation and neoblast differentiation (Agata et al. 2006). The knockdown planarians contained some progenitor cells even though they were not supplied by neoblasts. These progenitor cells normally became differentiated cells during regeneration after the first amputation in the DjHRJD knockdown group; thus, they could survive for 1 week. After 1 week, the progenitor cells were exhausted, which caused the planarians to die. This may explain why none of the planarians demonstrated abnormal phenotypes following RNAi knockdown of DjHRJDs during the first week, whereas many of them died later.

Conclusions

We have presented evidence of strong associations between HRJDs and a high regeneration ability in metazoans. Evolutionary conservation of the HRJDs appears to have been important for maintaining a high regenerative ability in metazoans, and species lacking this ability in Vertebrata, Arthropoda, and Nematoda have lost HRJDs during evolution. We propose that loss of the HRJDs during evolution is associated with the loss of a high ability for regeneration. These findings shed new insights into understanding the common molecular basis of high regeneration capability in metazoans.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.