Microhyla fissipes breeding adults were collected from Shuangliu, Chengdu, China (30.5825°N, 103.8438°E). All animal care and treatments were done as approved by Animal Care Committee, Chengdu Institute of Biology (Permit Number: 20150121003). Fertilized eggs were obtained from one pair of frogs, incubated in glass petri dishes (diameter, 160 mm) with dechlorinated tap water. After hatching, every 60 tadpoles was transferred to a plastic container (420 × 300 × 230 mm³) with 80 mm depth dechlorinated tap water at 25 ± 0.6 °C, and tadpoles were fed with spirulina powder once daily. All animals were maintained under the 12: 12 h light: dark cycle. Developmental stages were determined as described (21).

Six tadpoles at stage 39 (S39), S40, S41, and S43 each were randomly collected and anesthetized in 0.01% MS222 for tail morphological measurement. Three animals at each stage were used for hematoxylin-eosin (HE) staining and the rest were used for terminal deoxynucleotidyl-transferase mediated dUTP nick end labeling (TUNEL) analysis. For gene expression analysis, three biological replicates, each containing one animal tail, were used for RNA-seq at each developmental stage. To maximize the discovery of M. fissipes gene pools, heart, liver, spleen, lung, kidney, skin, ovary, testis from adult M. fissipes, tails (at S38, S40, S41, S43) and dorsal muscle (at S36, S40, S43, S45) were dissected for SMRT sequencing. All dissected tissue samples were immediately frozen in liquid nitrogen and stored at −80°C.

Apoptotic cells in the tail during metamorphosis were detected by using a TUNEL detection kit (Roche, 11684817910) according to manufacturer's instruction.

Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) and DNase I (Sigma, St. Louis, MO, USA) was used to remove DNA. RNA concentration was measured using Qubit 2.0 Fluorometer (Life Technologies, CA, USA) and RNA purity was checked by using Nanodrop 2000 Spectrophotometer (IMPLEN, CA, USA). Additionally, the integrity of RNA (RIN) was determined using Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) with RIN>8.0. Briefly, mRNA was purified from total RNA by using poly-T oligo-attached magnetic beads and chemically fragmented. First strand cDNA was synthesized by using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H–). Second strand cDNA synthesis was subsequently performed by using DNA Polymerase I and RNase H. Twelve cDNA libraries were generated by using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) following the standard protocols. The libraries were sequenced on the Illumina X-ten platform to obtain 150 nt paired-end reads at NovoGene (Beijing, China). All the raw data from Illumina sequencing were deposited in the NCBI Short Read Archive (SRA) database with the accession number PRJNA504611.

Tissue samples were homogenized and total RNA was isolated. The concentration and quality of RNA were assessed as above for RNA-Seq. The first full length strand cDNA was synthesized from 3 μg total RNA by using the SMARTer PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA). After PCR cycle optimization, the double-strand cDNA was generated with 11 PCR cycles, and the PCR products were separated with agarose gel-based size selection into cDNA fractions of length 1–2, 2–3, and 3–6 kb. Then, large-scale PCR was utilized to amplify the cDNAs. In the end, three SMRTbell template libraries were generated from the amplified cDNA by using Pacific Biosciences template preparation kit (Menlo Park, CA, USA) as per the standard protocol. Then SMRT sequencing was performed on a PacBio RS II platform at NovoGene (Beijing, China). Full-length transcripts were obtained by using Pacific Biosciences' SMRT Analysis Server 2.0 (22). Finally, full-length transcripts were corrected with Proovread Software (23) by using Illumina Hiseq reads. The proofread-corrected sequences after removing the redundant sequences with using CD-HIT-EST (24) were used as the reference sequences for further analyses. All the raw data from Illumina sequencing were deposited in the NCBI Short Read Archive (SRA) database with the accession number PRJNA504611.

The final corrected clean reads above were searched against the databases NR, Swiss-Prot and KOG by using diamond software version 0.8.22 with an E-value threshold of 1.0E-5 and KOG with an E-value threshold of 1.0E-3 to predict the gene identities. NT annotation was done with NCBI blast 2.2.28⁺ software version 2.2.28⁺ with an E-value threshold of 1.0E-5. Hmmer 3.0 package hmmscan was used for PFAM annotation with an E-value threshold of 0.01. GO annotations were determined based on the best BLASTX hit from the NR database using the Blast2GO software version 2.5 (E-value = 1.0E−6). KEGG pathway analyses were performed using the KEGG Automatic Annotation Server (KAAS, E-value = 1.0E−10). Clean reads were mapped to the full length reference transcripts by using Bowtie, and then gene expression level was calculated with RSEM (RNA-Seq by expectation-maximization) (25).

To compare the transcript expression levels among the four stages, the reads of each transcript were normalized to yield the number of sequenced fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) for the four developmental stages. Pairwise comparisons of transcript expression levels among the four stages were carried out by using the DESeq R package (1.10.1) to identify DETs. Transcripts were considered as DETs at q value < 0.05 after adjustment for the false discovery rate (FDR). For heatmap analysis, log10(FPKM+1) values were used for each tested transcript. Gene expression profile and log2 (ratios) values were used for K-means clustering.

DETs were subjected to GO enrichment analysis by using the GOseqR package based on the Wallenius non-central hyper-geometric distribution (26). KEGG pathway enrichment analysis of the DETs was done with the KOBAS software (27).

Studies on X. laevis indicate that tail resorption takes place mainly after stage 62 (S62) when rapid reduction in tail length occurs (28). This process involves apoptosis in essentially all tail tissues, especially, the epidermis, muscle, and notochord (29–31). To determine if M. fissipes undergoes a similar tail resorption process, we first compared the morphology of M. fissipes tadpoles at different stages of metamorphosis to that of X. laevis tadpoles (28). As shown in Figure 1, M. fissipes tadpoles at S39 resemble X. laevis tadpoles at S56, considered to be the early phase of metamorphosis (with stage 54/55 typically considered the onset of metamorphosis in X. laevis) (28). Under comparable rearing temperatures, M. fissipes tadpoles at S39 complete metamorphosis (reaching S45) in about 10 days, about ½ of the time that it takes for X. laevis tadpoles at S56 to reach the end of metamorphosis when tail resorption is complete (Figure 1). In particular, tail resorption, as defined from the onset of tail length reduction to complete resorption, requires only 3 days for M. fissipes tadpoles (from S42–S45) but 10 days for X. laevis tadpoles (from S61–S66) (Figure 1), indicating that M. fissipes tadpoles undergo a faster T3-dependent metamorphosis.

We next isolated tails from animals at 4 different metamorphic stages, S39, S40, S41, and S43 for further analysis. M. fissipes S39 is equivalent to X. laevis S56 (Figure 1). The only noticeable metamorphic change that takes place at this stage is the hindlimb development (28). For this reason, X. laevis at S56 is often referred to as a premetamorphic stage and for simplicity, we will refer M. fissipes tadpoles at S39 as premetamorphic animals as well. S40–S43 are metamorphic climax stages, although tail length reduction occurs only after S41. When tail sections at these different metamorphic stages were analyzed, it was clear that significant muscle atrophy was observed only at S43 (Figure 2A). In addition, TUNEL assay showed that while some apoptotic cells were observed by S40/41, drastic increase in apoptotic cells were found at S43 (Figure 2B), consistent with rapid tail length reduction. These findings indicate that M. fissipes undergoes a similar tail resorption program as that during X. laevis tail metamorphosis (29–31).

To facilitate the genome-wide analysis of the gene expression program underlying tail resorption, we first carried out SMRT sequencing of different M. fissipes organs at different stages. From 108,614,218 nucleotide sequences thus obtained, a total of 67,939 transcripts were assembled. The lengths of these assembled transcripts ranged from 103 to 8,847 bp with an average length of 1,599 bp and N50 of 1,956 bp (Table S1). Annotations through NR, NT, KEGG, Swiss-Prot, PFAM, GO and KOG databases showed that 52,881 transcripts (78% of the 67,939 transcripts) were annotated in at least one of the databases (Table S1).

To determine the gene expression program underlying tail resorption, we carried out Illumina RNA-Seq on the tail at the four metamorphic stages (S39, S40, S41, S43). The gene expression level was determined from the sequence data by using RSEM (25). 50,577 expressed transcripts with FPKM >0.3 were thus identified in the tail among these four metamorphic stages and the vast majority of transcripts (33,264 transcripts) were commonly expressed at all four stages (Figure 3).

Next, pair-wise comparisons of gene expression among the four stages were carried out, leading to the identification of 4,555 DETs in the tail during natural metamorphosis (Figure 4, Table S2). Among the 4 developmental stages, most DETs were found when comparing gene expression at S43 to any one of the other three stages, while relatively few DETs were found when the comparison was done among S39, S40, S41 (Figure 4). Consistently, Venn diagrams of the DETs between S39 and any one of the other three stages (Figure 5A) or the DETs between consecutive two stages (Figure 5B) showed very few commonly regulated genes. In contrast, Venn diagram of the DETs between S43 and any one of the other three stages showed that most of the DETs were common (Figure 5C). Finally, a heatmap of all 4,555 DETs showed that the vast of majority of the genes fell into two categories: highest and lowest expression at S43, respectively, but with similar expression levels among S39, S40, and S41 (Figure S1). These findings indicate that the tail at S43 has a very different gene expression program compared to those at the other three stages, while the tail gene expression programs at S39, S40, and S41 were very similar to each other, consistent with the lack of significant changes in the tail from S39 to S41 but drastic tail resorption by S43 (Figure 1).

To identify the critical gene regulation programs governing tail resorption, we further analyzed the 4,555 DETs in the tail during natural metamorphosis. First, K-means clustering was performed to divide the DETs into the 26 possible expression patterns (clusters) (Figure 6). In agreement with the findings above, most of the DETs fell into two groups of clusters: one group of clusters (the up-clusters #5, 8, 13, 21, 23, 26) in which the DETs were gradually upregulated by S43 and another group of clusters (the down-clusters #3, 9, 10, 12, 20, 22, 25) in which the expression of the DETs changed little from S39 to S41 but were downregulated significantly between S41 and S43. Since metamorphosis begins around S39 and drastic tail resorption starts between S41 and S43, we focused on these two groups of clusters. Analyses of gene functional categories and biological pathways with GO and KEGG, respectively, on the 1,030 genes in the up-clusters revealed no pathways or GO terms was significantly enriched in the DETs based on the more stringent criterion of q < 0.05. On the other hand, a number of significantly enriched KEGG pathways and GO terms were found based on a value of p < 0.01 (Tables S3, S4, Figure S2). Most significant among them was the proteolysis GO term with 60 DETs (Figure S2), essentially all of which were upregulated by S43 but changed little from S39 to S41 (Figure S3, Table S5). The upregulation of this GO term is consistent with the rapid protein degradation associated with tail resorption after S41.

Similar analyses on the 1,596 genes in the down-clusters revealed the enrichment of 21 KEGG pathways (Figure 7A, Table S6) and 499 GO terms (Figure 7B, Table S7) even based on the more stringent statistical criterion of q < 0.05. Many of the most significantly enriched KEGG pathways were related to metabolism, consistent with the tissue resorption process as the cells undergo apoptosis and reduce their own metabolism. The top 30 most significantly enriched GO terms were related to three higher level GO terms “Biological Process” (BP), “Cellular Component” (CC), and “Molecular Function” (MF) (Figure 7B), and were strongly associated with the protein import into nucleus, intermediate filament, myosin complex, motor activity and nucleotide binding (Figure 7B, boxed terms). When we carried out the analysis with directed acyclic graph (DAG), we found that the enriched GO terms within BP were all linked to a single node: the GO term “protein import into nucleus” (Figure 7Ca). On the other hand, the enriched GO terms within CC and MF could be linked to two GO nodes: the GO terms “myosin complex” and “intermediate filament” for CC (Figure 7Cb), and the GO terms “cytoskeletal components” and “muscle contraction” for MF (Figure 7Cc), respectively. The enrichment of these GOs among the down-regulated genes is again consistent with the tissue resorption program as the cells undergo apoptosis.

The discoveries of many KEGG pathways and GO terms associated with tail resorption above were based on genes that were coordinately regulated in up- or down-clusters. To determine if all DETs in a given KEGG pathway or GO term were coordinately regulated, we compared the gene regulation profiles of all DETs in selected pathways or GO terms that are likely important for tail resorption. It is well-known that matrix metalloproteinases (MMPs) are important for apoptosis and associated ECM degradation. A heatmap of all DETs encoding MMPs showed that all 9 DETs were coordinately upregulated dramatically at S43 (Figure 8A). Similarly, over 80% of the DETs belonging to the GO term “muscle function” (Figure 8B) or “mitochondrial function” (Figure 8C) coordinately down-regulated by S43, consistent with the apoptotic process taking place at S43.

A similar result was obtained on DETs in the KEGG pathways “glycolysis/gluconeogenesis” (Figure 9A) or “Alzheimer's disease” (Figure 9B). In the former, all enriched DETs were down regulated between S41 and S43 while in the latter, 9 out of 10 DETs were downregulated. Thus, DETs in KEGG pathways or GO terms that are associated with apoptosis and tissue resorption are coordinately regulated during tail resorption.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.