RNA-Seq Reveals Dynamic Changes of Gene Expression in Key Stages of Intestine Regeneration in the Sea Cucumber Apostichopus japonicas

Background Sea cucumbers (Holothuroidea; Echinodermata) have the capacity to regenerate lost tissues and organs. Although the histological and cytological aspects of intestine regeneration have been extensively studied, little is known of the genetic mechanisms involved. There has, however, been a renewed effort to develop a database of Expressed Sequence Tags (ESTs) in Apostichopus japonicus, an economically-important species that occurs in China. This is important for studies on genetic breeding, molecular markers and special physiological phenomena. We have also constructed a library of ESTs obtained from the regenerative body wall and intestine of A. japonicus. The database has increased to ∼30000 ESTs. Results We used RNA-Seq to determine gene expression profiles associated with intestinal regeneration in A. japonicus at 3, 7, 14 and 21 days post evisceration (dpe). This was compared to profiles obtained from a normally-functioning intestine. Approximately 5 million (M) reads were sequenced in every library. Over 2400 up-regulated genes (>10%) and over 1000 down-regulated genes (∼5%) were observed at 3 and 7dpe (log2Ratio≥1, FDR≤0.001). Specific “Go terms” revealed that the DEGs (Differentially Expressed Genes) performed an important function at every regeneration stage. Besides some expected pathways (for example, Ribosome and Spliceosome pathway term), the “Notch signaling pathway,” the “ECM-receptor interaction” and the “Cytokine-cytokine receptor interaction” were significantly enriched. We also investigated the expression profiles of developmental genes, ECM-associated genes and Cytoskeletal genes. Twenty of the most important differentially expressed genes (DEGs) were verified by Real-time PCR, which resulted in a trend concordance of almost 100% between the two techniques. Conclusion Our studies demonstrated dynamic changes in global gene expression during intestine regeneration and presented a series of candidate genes and enriched pathways that contribute to intestine regeneration in sea cucumbers. This provides a foundation for future studies on the genetics/molecular mechanisms associated with intestine regeneration.


Introduction
Regeneration, which can be described as the regrowth, or repair, of cells, tissues and organs, is a widespread phenomenon found among certain metazoans [1,2]. Besides its obvious advantages in terms of survival, regeneration also enhances the adaptive capacities of a species within its natural environment [3]. A considerable amount of research has focused on the fundamental biology of regeneration, including cell dedifferentiation, transdifferentiation, proliferation and migration. There is, however, a need for further studies into key aspects of regeneration, such as the target genes associated with the regenerative processes, and the molecular mechanisms involved [2,4].
Echinodermata possess spectacular regenerative capacity [1,3]. The sea cucumbers (Holothuroidea) are capable of regenerating most tissues and organs, including the intestine, respiratory tree, gonads and the body wall. They are thus considered as excellent models for organ regeneration studies [1,5,6,7]. Regeneration has been the topic of considerable research, with an emphasis on visceral (intestine) regeneration, including the histological changes (for example, tissue layer changes in the intestine wall and changes in the radial nerve cord) and cellular events (cell origin, migration and proliferation) associated with such regeneration [8,9,10,11,12,13,14]. Nevertheless, due to the lack of information on the genome of sea cucumbers, only relatively few genes -such as Ependymin, Wnt9, Bmp1 and Serum amyloid A etc. -have been screened, recognized and analyzed, using traditional methods [11,15,16]. A greater effort on the study of the molecular/genetic mechanisms involved in intestine regeneration in the Holothuroidea is therefore needed.
Pablo et al. constructed an EST library (,5000 ESTs) and used a microarray technique to analyze gene expression profiles during intestine regeneration in the cucumber Holothuria glaberrima [7,10]. However, approximately 5000 ESTs generated by Sanger sequencing were not nearly sufficient to cover the transcriptome of sea cucumber. Moreover, because of limitations due to lower throughput, high background noise and lower sensitivity of microarray, it's hard to obtain high-quality global gene expression profiles. The development of an ultrahigh-throughput sequencing RNA-Seq technique, however, provided a powerful alternative technology for gene expression profiling [17,18]. The RNA-Seq technique has the potential to overcome microarray limitations and provide an expression profile with a greater and reproducible dynamic range. RNA-Seq has been successfully applied to many research projects, including the transcriptional landscape of the yeast genome, research on sulfur-deprived Chlamydomonas cell survival and mapping, as well as the quantification of mammalian transcriptomes [19,20,21].
In a previous study, our team constructed a transcriptome of the regenerative body wall (after 4 days of regeneration) and intestine (7dpe) in A. japonicus by 454 life sequencing which, to a large extent, increased the EST information required for a molecular mechanism study on regeneration [22]. This research was, however, of a preliminary nature, in contrast with the present study on intestine regeneration in A. japonicus, which includes five key stages: wound healing (0-3dpe), blastema formation (3-7dpe), lumen formation (7-14dpe), intestine differentiation (14-21dpe) and growth (21dpe-) [23]. The following processes took place during different stages: Stage I: activation of damage repair capacity in response to injury, so as to accumulate energy in preparation for regeneration [24]; Stage II: formation of the blastema, which encompasses cell migration, dedifferentiation and transdifferentiation; Stage III: luminal epithelium growth by cell division and migration, which results in lumen formation of new intestine; Stage IV: gradual development of the intestine to form a complete structure in which the digestive and absorptive functions are restored. Stage V: enlargement of the intestine to reach its original normal size. This process can thus be described as a continuous dynamic change in the gene expression that underlies the molecular mechanisms involved in regeneration.
In the present study, we used RNA-Seq to determine the global dynamic changes in the gene expression profile, occurring during the intestine regeneration in A. japonicus. For robustness, the reference library, which integrated nine transcriptomes, (in all ,30,000 isotigs), included all tissues of sea cucumbers at every developmental stage (embryo, larva, white juvenile and black juvenile) and under different physiological conditions, i.e. active, regeneration and aestivation conditions [22,25]. Based on this premise, we compared the expression profiles of the intestine in a Normal situation (with intact intestine) with those recorded at 3, 7, 14, 21 days post evisceration (dpe). The global dynamic changes in gene expression during each of these stages were analyzed. The sampling time covered all stages of intestine regeneration. As expected, a total of up to 5119 DEGs associated with intestine regeneration were identified. Furthermore, gene ontology (GO) and pathways enrichment analysis were conducted for investigating the main function of DEGs and providing an overview of the gene regulation process during intestine regeneration. For example, ''organic substance transport'' was only observed at 3dpe, which indicated that during the early stage of the process sufficient nutritional and energy resources were assembled to enable the commencement of regeneration. Besides the expected pathways, the enriched ''Notch signaling pathway'', ''ECMreceptor interaction'' and ''Cytokine-cytokine receptor interaction'' were substantially activated during intestine regeneration. It should be noted that little overlap was observed between the screened top DEGs when comparing our results with those of Pablo et al. To sum up, our work may provide a more representative basis for the future study of molecular mechanisms associated with intestine regeneration.

Ethics Statement
Not applicable. Our research did not involve human participants or samples.

Animals
Adult sea cucumbers, A. japonicus (70-100 g), were collected from the coast of Qingdao, Shandong Province and cultured in a laboratory for 1 week in sea water at 15-17uC, prior to setting up the experiment. Evisceration was induced by injecting about 2 ml 0.35 M KCl into the coelom [6,10,16]. The non-eviscerated sea

RNA preparation and cDNA synthesis
Total RNA from regenerative (from animals at 3 dpe, 7 dpe, 14 dpe, 21 dpe) and normal intestines (control) was extracted and DNase-treated using RNeasy Mini Kit and RNase-Free DNase Set (Qiagen, Germany), following the manufacturer' s instructions. The quality and concentration of RNA were measured by NanoDrop 1000 (Thermo). Total RNA from 15 individuals per stage were pooled. The starting amount of RNA was 1ug per pool (including 15 individuals). After that, the mRNA in total RNA was enriched by using the oligo(dT) magnetic beads. After the addition of fragmentation buffer, the mRNA was interrupted and formed short fragments (,200 bp). The first strand cDNA was then synthesized by the random hexamer-primer, using the mRNA fragments as templates. Buffer, dNTPs, RNase H and DNA polymerase I were added to synthesize the second strand. The double strand cDNA was purified using the QiaQuick PCR extraction kit, and washed with EB buffer so as to facilitate end repair and the addition of a single nucleotide A (adenine). Finally, sequencing adaptors were ligated to the fragments. The required fragments were purified by means of agarose gel electrophoresis and enriched by PCR amplification. The library products were then ready for sequencing analysis via Illumina HiSeq TM 2000 (BGI, Shenzhen).

Sequence annotation, assessment and gene expression levels
The original image data was transferred into sequence data by base calling (using CASAVA version 1.5+, defined as 'raw reads') and saved as fastq files that include the detailed reads sequences and the reads quality information. Specifically, if the sequencing error rate is denoted as E, Illunima HiSeq TM 2000, and the base quality value is denoted as sQ, the relationship is as follows: sQ = 210lgE. To obtain clean reads, it is necessary to remove the dirty raw reads prior to data analysis. This process encompasses the following procedures: (1) remove reads with adaptors; (2) remove reads in which unknown bases are greater than 10%; (3) remove low-quality reads (the percentage of the low quality bases of quality value #5 is more than 50% in a read).
A comparison of ''clean reads'' by SOAPaligner/soap2 was carried out using the reference databases from large-scale transcriptome profiling [22,25]. Assessment of sequencing was conducted in five steps: sequence reads quality assessment, Figure 1. The change level of global differential expression genes during intestine regeneration in sea cucumber A. japonicus. The principle "FDR#0.001 and the absolute value of log2Ratio$1'' was used as a threshold to screen DEGs. The figure showed that not only a large number of genes were differentially expressed, but the change fold of differentially expression was at a high level, especially at the early stage of intestine regeneration. doi:10.1371/journal.pone.0069441.g001 statistics of alignment analysis, sequencing saturation analysis, distribution of reads on reference genes, and gene coverage analysis. Mismatches of no more than two bases were allowed in the alignment. The number of clean tags was calculated and normalized, using the RPKM method (Reads Per kb Million reads) [21].

Screening of differentially expressed genes (DEGs)
According to the Audic et al.'s algorithm, a strict algorithm to identify differentially expressed genes between two samples was developed [26]. The P value was used to detect the difference in gene expression in two different samples [27]. If the P value is less than1E-238 (P,1E-238 ), the P value is given the value ''0''. The FDR (False Discovery Rate) was estimated, to determine the threshold of P-value. The ''FDR#0.001 and the absolute value of log 2 Ratio$1'' was used as the threshold to judge the significance of gene expression difference. All DEGs were screened by program ''stastistic'' written using Language C. Previous studies have demonstrated that 48 genes differentially expressed responding to the injection [28]. However, the intestine regeneration may involve the regulation of immune genes. Hence, These genes were not subtracted from the regeneration profile.

Gene ontology and Pathway enrichment analysis
Gene Ontology (GO) is an international standardized gene functional classification system which offers a dynamic-updated controlled vocabulary and a strictly defined concept to comprehensively describe properties of genes and their products in any organism. GO has three ontologies: molecular function, cellular component and biological process. The basic unit of GO is GOterm. Pathway-based analysis helps to further understand genes biological functions. KEGG is the major public pathway-related database. The first stage is to map all DEGs to GO terms in the database (http://www.geneontology.org/) using program Blas-t2GO and KEGG (http://www.genome.jp/kegg/) using KEGG Automatic Annotation Server (KAAS) and calculated gene numbers, after which the hypergeometric test is used to find significantly-enriched GO terms and pathways in DEGs. GO terms conforming to p-value through Bonferroni Correction#0.05 were defined as significantly enriched GO terms and pathways [29].

Real-time PCR validation
To validate RNA-seq results, some significant DEGs were chosen to carry out Real-time PCR. According to the sequence information in the transcriptome database, primers were designed     for optimal performance using the primer3 (Table S1). The RNA that was used to synthesize cDNA was the same as that used to RNA-seq. The first strand cDNA was synthesized in 25 ml reaction system as follows. Firstly, 4 ml RNA and 1ml oligodT18 was degenerated at 70uC for 5 min. Then, 1 ml M-MLV reverse transcriptase (Promega), 5 ml M-MLV buffer (25 mmol/L KCl, 10 mmol/L Tris-HCl, 0.6 mmol/L MgCl2, and 2 nmol/L DTT, pH 8.3), 5 ml dNTP, 1 ml ribonuclease inhibitor and 8 ml RNasefree water were added at 42uC for 1 h. The synthesized cDNAs were diluted with RNase-free water and stored at 280uC for subsequent quantitative real-time PCR. The mRNA expression levels were determined using the SYBR GreenH real-time PCR assay with an Eppendorf MastercyclerHep realplex (Eppendorf, Hamburg, Germany). The amplification volume was 25 mL, containing 12.5 mL of SYBR GreenMasterMix (Takara), 0.5 mL (each) of forward and reverse primer (10 mM), 1 mL of diluted cDNA, and 10.5 mL of RNase-free water. Thermal cycling was as follows: (1) 95uC for 5 s; (2) 40 cycles at 95uC for 10 s, 60uC for 20 s and 72uC for 30 s. The melting curve analysis of the amplification products was done to demonstrate the specificity of the PCR products. NADH dehydrogenase (NADHF, NADHR, Table S1) was used as a housekeeping gene for internal standardization [7,16,22,30]. The 22ggCT method was used to analyze the expression level. Five biological pools (N = 5), with each pool (3 dpe, 7 dpe, 14 dpe, 21 dpe and Normal) being mixed with three different sea cucumber individuals (n = 3 animals) were analyzed, because of the small quantity of regenerative intestine. All data were given as mean 6 S.E. (N = 5) and the level of statistical significance was set at P,0.05. Analysis was carried out using SPSS16 software.

Reads sequencing, quantification and assembly
A total of 4 868 208, 4 727 453, 5 030 570, 4 715 682 and 4 877 984 reads were sequenced using RNA-Seq technique in Normal, 3dpe, 7dpe, 14dpe and 21dpe libraries, respectively, which have been submitted to NCBI (accession NO. GSE44995; Table 1). After trimming the dirty raw reads (reads with adaptors, more than 10% unknown bases and over 50% bases of quality value #5), 4 848 595, 4 708 854, 5 011 237, 4 696 545 and 4 855 872 clean reads were obtained in Normal, 3dpe, 7dpe, 14dpe and 21dpe libraries, respectively. (Table 1). The map reference transcriptome from 454 life sequencing contains nine libraryembryos (4 h, 23 h); larvae (30 h, 6d,8d, 10d); white juveniles (16d, 22d); black juveniles (32d, 37d); female gonads; male gonads; intestine, respiratory trees and coelomic fluids from active adults; intestine, respiratory trees and coelomic fluid from aestivating adults; and regenerating intestine and body wall. Therefore there are abundant expressed sequence tags (29 667 ESTs) in the reference transcriptome to cover the whole transcriptome [22,25]    Sequencing Saturation Analysis Saturation of the library is related to the depth of sequencing. When barely any new genes are detected, sequencing reaches saturation. With the number of reads increasing, the number of detected genes increased. The more reads were sequenced, the better results were achieved. Figure S1 indicates that, the percentage of detected genes was less than 2% when the number of reads increased from 4 to 5M. ,5 M reads sequenced in each library provide sufficient coverage depth to achieve good results. Sequencing saturation analysis also presented the transcriptome size in the five libraries: R3 = R7.R14.R21.Normal.

Differentially expressed genes (DEGs)
The RNA-seq results showed sequential large-scale genes expression profiles at 3dpe, 7dpe, 14dpe and 21dpe during intestine regeneration in A. japonicus. The transcripts detected with at least two-fold differences (log 2 Ratio$1) are classified as DEGs (FDR#0.001) (Figure 1). Scatterplot was applied, to graphically describe the expression profile of global differential genes. A large number of genes were differentially expressed at a high level between normal and regenerative libraries, especially at the early stage of intestine regeneration (Figure 1). For example, when compared to Normal, 2415 up-regulated genes (10.64%, 2415/ 22691) and 1099 (4.84%, 1099/22691), the down-regulated genes were observed at 3dpe; and 2520 (11.11%, 2520/22691) upregulated genes and 1084 (4.78%, 1084/22691) down-regulated genes were screened at 7dpe. The number and expression level of DEGs gradually decreased in the 21dpe library, and only 1012 (4.46%, 1012/22691) up-regulated genes and 554 (2.44%, 554/ 22691) down-regulated genes were observed. In addition, we also developed a strategy to focus on key regeneration genes by comparing gene expression between two libraries (3dpe and 7dpe) during the early stages of regeneration. Only 280 up-regulated genes and 256 down-regulated were observed at 3dpe, in comparison with 7dpe (Table S2).
There were occasional No Reads observations in some libraries, resulting in the high change fold of DEGs between two libraries. In such cases, though, it was difficult to assess the real difference of DEGs. In addition, some mapped genes with no reference information in the NCBI database could not be annotated. These unknown genes were hard to analyze and in such situations we here showed the representative top 10 DEGs, based on the criteria (Tables 2,3,4,5 and 6). When the number of reads mapping to one gene is 0 in one library, the counterpart must be more than 100 in the other library. As shown on the top DEGs list, most DEGssuch as low density lipoprotein-related protein 2-like, orthodenticle, cyclin B3 Hu/elav isoform 7, Rp2 Lipase, proprotein convertase subtilisin/kexin type 9 -were significantly differentially expressed throughout the intestine regeneration ( Figure 2, Tables 2,3,4 and 5). Some special DEGs-kelch-like ECHassociated protein 1, HORMA domain containing 1-like and Human Fc gamma BP et al -were screened at 3dpe Vs 7dpe (Table 6). Unknown DGEs with high change fold were also notable, since they might correspond to novel sea cucumbersspecific sequences associated with their striking regenerative capacities. Unknown DEGs with high change fold can be found in Table S3, S4, S5, S6 for further comprehensive and specific research.
The genes unique to a certain stage are very important. 386 genes that were unique to 3dpe were found on Table S8, 117 of which could be annotated, such as f-box only protein 47-like, Notch homolog Scalloped wings-like and heat shock protein 70  Table S8 were unique to 7, 14 and 21 dpe respectively (Table S8). Though so many genes were unique to a time point, these genes were low expressed. only one or a few unique reads could be mapped to them.

Validation of Results by real-time PCR
To validate the expression profiles, top 10 up-regulated genes (abbreviation: LDP, GL2, speA, RAP, OTD, SCF, CB, TFP, HE7 and His1) and top 10 down-regulated genes (abbreviation: Rp2, CRBP, FCBP, PCSK, AMY, GAP, LOC, TAG, LYS, CBG), selected from Table 2, were applied to real-time PCR at four stages of regeneration (normal, 3dpe, 7dpe and 21dpe) ( Table 7; Figure 3). NADH dehydrogenase was taken as a reference gene to normalize gene expression data, based on a previous investigation [7,16,22,30]. All primers were tested to be effective on PCR amplication. One peak in the melting curve was detected in all real-time PCR, which indicated that the all PCR product was specifically amplified (data not shown). Real-time PCR results indicated that the expression level of these selected genes reached its highest at 3dpe. The overall trend of Real-time PCR-based expression patterns among these selected genes was similar to those obtained by RNA-Seq based detection. However, the changes fold of most data measured by Real-time PCR was smaller than those by RNA-seq (Table 7).

Gene ontology analysis of DEGs
DEGs, that were involved in certain biological functions during intestine regeneration, were evaluated by GO (gene ontology) analysis. GO terms were assigned to 1111(31.6%, 1111/3514), 1230 (34.1%, 1230/3604), 905 (35.3%, 905/2561) and 498 (31.8%, 498/1566) DEGs for 3, 7, 14 and 21dpe, respectively ( Figure 4). The DEGs associated with cellular and metabolic processes in the category 'biological process'; the DEGs associated with cell and cell part in the category 'cellular component' and the DEGs associated with binding and catalytic activity in the category 'molecular function' were notably abundant throughout the intestine regeneration. The GO distribution of the DEGs was highly similar at different regeneration stages. Moreover, we screened the significantly enriched GO terms classified in biological function category. As indicated in Table 8 , some GO terms -'pancreas development', 'gene expression', 'translation' and 'reproduction' -were significantly enriched at 3, 7 and 14dpe. At 3dpe, only six GO terms (pancreas development, gene expression, translation, reproduction, organic subtance transport and reproductive process) were significantly enriched, which were mainly associated with development (4) and transport (1). GO term 'organic subtance transport' was only enriched at 3dpe. 15 GO terms (ribonucleoprotein complex biogenesis, cellular component biogenesis at cellular level, ribosome biogenesis and primary metabolic process et al)were significantly enriched at 7dpe, which were associated with development (5) and metabolism (10) GO term 'organ development' was only enriched at 7dpe. In particular, as much as 28 GO terms were significantly enriched at 14dpe, which were associated with development (4), metabolism (14), cell cycle (6), signal transduction (2) and transport(2). 14 GO term (cell cycle, proteasomal protein catabolic process and signal transduction in response to DNA damage et al) were only enriched at 14dpe (Table 8, corrected p-value#0.05). No significantly enriched GO terms in DEGs were detected at 21 dpe. Besides, two unique GO terms-'' serine family amino acid metabolic process'' and ''polyol metabolic process'' were significantly enriched at 3dpe Vs 7dpe (Table 8).

Pathway enrichment analysis of DEGs
Intestine regeneration associated biological pathways were evaluated by enrichment analysis of DEGs. Significantly enriched metabolic pathways and signal transduction pathways were identified, and listed in Table 9 (Q value ,0.05). Seven significantly enriched pathways for the up-regulated DEGs were screened. ''Cytokine-cytokine receptor interaction'' and ''Notch signaling pathway'', ''Ribosome'', ''Spliceosome'', ''RNA transport'', ''DNA replication'' and ''Ribosome biogenesis in eukaryotes'', were related to the biogenesis of protein, DNA and RNA. There were eight significantly enriched pathways for the downregulated DEGs, which were involved in digestion and absorption of vitamin fats and carbohydrates, the Renin-angiotensin system, Pancreatic secretion, and metabolism of Glycerophospholipid, alpha-Linolenic and Starch and sucrose. In addition, ''Retinol metabolism'', ''Protein digestion and absorption'', ''PPAR signaling pathway'', ''Glutathione metabolism'' and ''ECM-receptor interaction'', including some up-regulated and down-regulated DEGs, were also enriched. For example, DEGs for adipocyte differentiation were up-regulated (12) and DEGs for lipid metabolism were down-regulated (15) in the PPAR signaling pathway. Only one unique pathway'' Glycine, serine and threonine metabolism'' were found at 3dpe Vs 7dpe. The DEGs in enriched pathways were listed in Table S7.

Individual key DEGs
Key genes associated with the regenerative process were classified into three groups: developmental genes, extracellular matrix (ECM) associated genes, and cytoskeletal genes.
a. Developmental genes. Regeneration and development share similar mechanisms, so developmental genes are excellent candidates for future study on molecular mechanisms of regeneration [4]. The expression profiles of well-known developmental genes are therefore summarized in Table 10. Most genes (Wnt, Hox, BMP and syndecan) were up-regulated; while krueppel like6 were down-regulated during regeneration. In addition, genes from the same family showed different expression patterns. For example, Hox1 and Hox3 were up-regulated at early stage 3 ,7 14 dpe and reached a peak at 3dpe. No changes were observed for Hox9/10 at 3dpe, but its expression remained at a high level from 7dpe to 21dpe. In the case of Hox11/13, the up-regulating trends occurred from 14dpe. b. ECM-associated genes. Quinones et al. have found that, in the sea cucumber H. glaberrima, the ECM content undergoes significant changes during intestine regeneration and that those changes are closely related to MMPs activity [13]. In this study six family genes (collagen, tenascin, laminin, MMPs, spondin and fibropellin) associated with ECM were screened. All of these genes, except MNP14, up-regulated at one or several stages of intestine regeneration (Table 11).  c. Cytoskeletal genes. Changes in cytoskeletal protein synthesis during regeneration have been demonstrated in many regeneration tissues [31,32,33]. All cytoskeletal genes showed significant changes at 3 or 7dpe (Table 12). Some (alpha-tubulin, beta-tubulin and actin) were up-regulated; while others (gamma tubulin, myosin and gelsolin) were down-regulated.

Discussion
In this study, an RNA-Seq technique was used to construct the large-scale dynamic expression profile during intestine regeneration in A. japonicus. Real-time PCR was used to verify gene expression profiles. We found that the results of the verifications were comparable to the gene expression profiles. They had similar trends of up-and down-regulation. However, the change scales of these gene expressions, especially for down-regulated genes, were larger for RNA-Seq, compared with those for Real-time PCR. The RNA-Seq was probably more sensitive in terms of determining gene expression levels, particularly for low-abundance transcripts [17]. The individual variation may be another reason for the difference in change fold. For RNA-Seq assay, 15 individuals were pooled to run; while for Real-time PCR, 5 samples were used separately.
This was to be expected, as intestine regeneration is a very complicated process that includes epithelial cell migration and proliferation, contraction of injured area, changes in cellular function, and communication and interaction of multiple cell types. [13,34,35]. This result is not completely consistent with a report by Pablo et al in which ,2000 ESTs -a third of the total -showed a different expression in H. glaberrima. [7]. Besides which, there was almost no overlap in top DEGs between these two studies. The difference scales of top DEGs were much larger in our results than Pablo et al. The number of ESTs matched to known genes in the reference libraries was also different: only ,1500 ESTs (,20%) in the H. glaberrima library and ,13000 ESTs (47%) in the A. japonicus library. The reasons for these differences might be due to three factors. First, most H. glaberrima ESTs were only obtained from regenerative intestine, so the original data was biased in favour of the genes associated with regeneration, while in our studies, the ESTs were obtained from transcriptomes of all tissues at multiple developmental stages and varied physiological conditions found in Table 11. Expression profile of ECM associated genes during intestine regeneration.   [17]. In a previous study we found 324 up-regulated genes and 80 down-regulated genes during regeneration (body wall at 4dpe and intestine at 7dpe) [22], which is much lower than the number of up-regulated genes (2520) and down-regulated genes (1084) in intestine regeneration at 7dpe. This may be due to the sequencing depth of RNA-Seq (,5M reads) being larger than that of 454 sequencing (182,473 reads). Moreover, single transcriptome data was insufficient to cover the whole transcriptome, so fewer genes were assembled and recognized. It thus becomes clear that our study is a useful supplement to research carried out on a molecular basis as well as being a more representative collection for studying intestine regeneration. Gene expression profiles changed with tissue and cellular events during intestine regeneration. During the early stage of intestine regeneration a large number of genes, such as seawi, presenilin, matrix metalloproteinase and Tbx2/3, were abundantly expressed, which contributed to many cellular events, such as cellular migration and differentiation (Table S3, 4). As regeneration proceeded, genes such as helicase, growth arrest-specific 8-like, DnaJ and cyclin-dependent kinase were up-regulated to promote cell proliferation at 7-and14-dpe (Table S4, S5). Lost tissues were then restored and regeneration entered into the intestine growth stage (which included intestinal wall thickening and other processes) so that the number and change fold of DEGs gradually decreased. We have listed the top significant DEGs as a reference for further studies. In our previous transcriptome comparative study [22], genes in the top significant DEGs list, such as proprotein convertase subtilisin/kexin type 9 and cyclin B3, also showed a high degree of differential expression. Orthodenticle, Hu/elav, regeneration-associated protein and speedy A, that were deemed to be associated with regeneration or development in some studies, are good candidate genes for studying intestine regeneration [7,36,37]. Nevertheless, many genes in the top DEGs list, such as low-density lipoproteinrelated protein, cellular retinol-binding protein, and proprotein convertase subtilisin/kexin type 9, were mainly related to intestinal digestion, absorption, and metabolism. It was thus obvious that their expression levels changed when the intestine lost its function [38,39,40]. When this is taken into consideration, it is clear that a search for important regeneration-associated genes should not be strictly confined to the top DEGs as this was not sufficient to search important regeneration-associated genes. To focus on key genes, we compared gene expression between 3dpe and 7dpe stages. Previous studies have demonstrated that the blastema for developing a new intestine forms at 7dpe [8]. The early stage is therefore crucial for regeneration. The method that we used was successful even though only 536 DEGs were observed. Besides H3.3 histone, ribosomal protein L23a and other genes associated with DNA replication, other previously-unnoticed genes attracted our interest (Table 6). For example, it has been demonstrated that annexin could regulate regeneration by affecting the synthesis of EGF [41,42]. Attractin, a member of the CUB family, plays a crucial role in myelination in the central nervous system, so its over-expression could relate to the reconstruction of the nervous system of regenerative intestine [43,44]. In addition, over-expression of the solute carrier family and the tetraspanin gene were observed in regenerating rat liver [45,46]. The genes unique to a certain stage are very important. For example, f-box only protein and heat shock protein 70 were unique to 3dpe (Table S8). F-box protein only mediates the ubiquitination and subsequent proteasomal degradation of target proteins, which play an important role in degrading protein to promote cell migration. Heat shock protein 70 may response stimuli and help sea cucumbers to resist abnormal state. Besides, about two thirds unique genes were ''unknown genes'' which need further specific research.

Function classification of DEGs
It is to be expected that the significantly enriched GO terms are mainly associated with development, gene expression, metabolism and cell cycle. In many studies, a number of key developmental regulators have been found to be up-regulated during the regenerative process [2]. Genes mapped to ''gene expression'' terms, such as Zinc finger protein, ribosomal protein, and Dap5 protein, were activated to produce mature gene products during intestine regeneration. The results also demonstrated some marked trends in gene functioning during the different stages of regeneration. During the early stages, nutrients were rapidly accumulated by genes executing ''substance transport'', to provide energy for regeneration and a large number of genes associated with development and cellular events were over-expressed to initiate regeneration. At 7dpe, genes associated with metabolism were activated during the rapid blastema formation stage, when anabolism and catabolism activities (macromolecule biosynthesis and metabolism, protein synthesis and degradation) were dominant [8]. At this stage, intestine wall undergo differentiation, when undifferentiated cells gradually differentiate into different cell types. So anabolism and catabolism activities were dramatic. At the 14dpe stage, there was an obvious change in the expression of genes mapped to ''cell cycle'', which was consistent with a previous study in which it was noted that, during this stage, a large number of dividing cells were found in serosa/muscle and mucosa [8]. Besides, at this stage, DGEs were mapped to signal transduction term (signal transduction in response to DNA damage and DNA damage response, signal transduction by p53 class mediator). At the 14dpe stage, the rate of cell proliferation was high, so DNA replication activities was dramatic. Hence, signal transduction was started to response to DNA damage during DNA replication.

Pathways that play a key role in regeneration
Pathway enrichment analysis identified the most significantlyaffected pathways during intestine regeneration. It is no surprise that ''Ribosome'', ''Spliceosome'', ''RNA transport'' and''DNA replication'' pathways were observed, because protein synthesis is needed for initiating and maintaining regeneration, and the digestion and absorption of vitamin, fat and carbohydrate were inevitably weakened because of the loss of intestine function. The metabolism of Glycerophospholipid, alpha-Linolenic acid and Starch and sucrose were also affected due to hypometabolism (associated with high oxygen consumption) during intestine regeneration [24]. The ''Notch signaling pathway'' has been reported to play a key role during the regeneration of human intestine epithelia [47], avian retina [48], rat liver [49], zebrafish heart [50] and newt retina [51] by controling multiple cell differentiation, proliferation and apoptotic programs [52]. In the preset study, many genes involved in this pathway -such as delta, notch, histone deacetylase and presenilin -were identified and upregulated during intestine regeneration. Another noticeable pathway is the ''ECM-receptor interaction'' pathway. In this pathway, genes -for example, collagen, laminin, tenascin, axonemal dynein and CD36 -that are required for cellular activities such as cell adhesion, migration and differentiation, were also affected [53,54]. In addition, genes belonging to the FDGF (fibroblast-derived growth factor), TNF (Tumor necrosis factor) and TGF (Transforming growth factor) family were observed to be up-regulated in the ''Cytokine-cytokine receptor interaction'' pathway. This phenomenon has also been reported during liver regeneration [55,56].
Special types of genes associated with regeneration a. Developmental genes. Many developmental genes were differentially expressed during intestine regeneration, Which is not surprising, since regeneration and development are closely related. It appears that developmental genes are reactivated during the regenerative process [4]. Therefore, the developmental genes that we screened were excellent candidates for regeneration study.
The Wnt gene family, which encodes a large group of secreted cysteine-rich proteins, is one of the most intensively studied gene families and is a critical mediator of key cell-cell signaling events during development and regeneration [57,58]. We found that the expression patterns of wnt4, wnt6 and wnt8 were not exactly the same during intestine regeneration in A. japonicus. As previously reported, distinct Wnt genes have dual roles in regeneration, with some genes promoting regeneration while others inhibited regeneration [59]. The wnt4 and wnt6 genes were up-regulated during intestine regeneration of A. japonicus and have been shown to be involved in epithelial remodeling, myogenesis, epithelialmesencheymal transformation and the endomesoderm genes regulatory network [60,61,62,63]. Only a very low level of transcript of wnt8 was, however, detected during intestine regeneration in A. japonicus, indicating that this gene may not contribute to this process. Results from another study have also indicated that wnt8 transcripts were undetectable by in situ hybridization during zebrafish fin regeneration [64].
The Hox genes that encode transcription factors are an important gene family that determines positional identity along the anteroposterior axis in a wide range of metazoans and guides tissue differentiation [65,66,67,68]. These genes determine the progenitor cell fate and positional identity during wound healing and regeneration [66,67]. As is the case for Wnt genes, the four hox genes also show different expression patterns. The hox1 and hox3 genes are activated at the early stage of intestine regeneration, while hox9/10 and hox 11/13 are stable until a later stage, which is consistent with results from research on Xenopus limb regeneration. This research indicated that Xhoxc10 was up-regulated when cells dedifferentiated in the blastema,whereas XHoxa13 was slightly reexpressed during later regeneration, while XHoxd13 was not expressed until a later regeneration phase [69]. This could be because different Hox genes play distinct roles in regeneration, including wound healing, setting up regeneration, and specifying positional identity [70,71].
Bone morphogenetic proteins (BMPs) are multi-function growth factors that regulate cell proliferation, survival, differentiation, and apoptosis [72]. An increased level of attention has recently been focused on the role of BMPs in regeneration, such as tail and limb regeneration in Xenopus tadpoles [73], newt lens regeneration [74], digit regeneration in fetal mice [75], and bone regeneration in rats [76]. Our results have shown that BMP1 was up-regulated during intestine regeneration, which is consistent with Mashanov et al.' s study on intestine regeneration in H. glaberrima [11]. They suggested that Bmp1 may be involved in regulating the folding of the luminal epithelium and gut looping by modulating morphogenetic movements [11].
In our study we found that two other two genes, syndecan and krueppel-like6, were affected by regeneration. These are reported to be associated with muscle and axon regeneration [77,78]. Syndecans, a member of the transmembrane heparin sulfate proteoglycans (HSPGs) family, are implicated in muscle differentiation and myoblast development to direct muscle regeneration [78,79,80]. Krueppel-like factors, a family of zinc finger transcription factors, were deemed to regulate the intricate gene programs associated with axon regeneration [77,81].
b. Extracellular matrix (ECM) associated genes. ECM remodeling related to organ morphogenesis and regeneration, has been demonstrated in many animals: liver regeneration in rats [82], skeletal muscle regeneration in mouse [83], limb regeneration in newts [84], skin regeneration in mammals [85] and intestine regeneration in sea cucumber H. glaberrima [13]. Collagen, tenascin, laminin, spondin and fibropellin were similarly observed to be up regulated during intestine regeneration. ECM remodeling is regulated by the proteolytic activities of matrix metalloproteases (MMPs), which degrade ECM components [13,86]. Many studies have suggested that MMPs are required for tissue remodeling and regeneration, such as newt limb regeneration [87], zebrafish fin regeneration [88], mice epithelial regeneration [89] and skeletal muscle regeneration [90]. Our results also identified up-regulated expression of all MMP genes except MMP14. We have concluded that ECM genes contributed to intestine regeneration in A. japonicus.
c. Cytoskeletal genes. Myogenesis is one of the most remarkable characteristics of intestine regeneration [35]. Many myocyte events occurred during intestine regeneration of sea cucumbers, including the formation of spindle-like structures, muscle component dedifferentiation, and muscle precursor cell division [35,91,92]. Our results showed a dramatic expression variation in cytoskeletal genes. The Cytoskeletal genes-''tubulin'' [93,94,95], ''actin'' [96], ''myosin'' [97] and ''gelsolin'' [98] have been proved to be related to regeneration. An interesting finding is the downregulated expression of myosin and gelsolin during intestine regeneration. A similar result was also reported in H. glaberrima [7]. The expression of gelsolin was also observed to be down-regulated in regenerating skeletal muscle [99]. This may relate to the observation that there is no muscle layer in the newborn intestine at the early stage of regeneration, so myosin was not needed to regulate muscle motility [8]. Gelsolin has been demonstrated to be a regulator and effector of apoptosis, which would explain why it was down regulated during Cytoskeletal remodeling [100].

Conclusion
The intestine regeneration of sea cucumbers is a complex process where thousands of genes (more than 15%, ,3500 ESTs) showed significantly differential expression during intestine regeneration. RNA-Seq can be used for analyzing variantion in gene expression between different samples. Analysis gained by this study will be useful for future regeneration studies on sea cucumbers. Moreover, we demonstrate that intestine regeneration in the sea cucumber is a good model for studies to identify and characterize the molecular basis of organ regeneration. Figure S1 Sequencing saturation analysis in Normal, 3dpe, 7dpe, 14dpe and 21dpe libraries.