Microhyla fissipes breeding adults were collected from Shuangliu, Chengdu, China (30.5825${}^{\circ }$N, 103.8438${}^{\circ }$E) in June 2016. All animal care and treatments were done as approved by Animal Care Committee, Chengdu Institute of Biology (Permit Number: 20150121003). All animals were obtained and maintained as previously described [7]. Developmental stages were determined as described [9].

Tadpoles at stage 18 (S18), S19, S21, and S28 were randomly collected and anesthetized in 0.01% MS222 for tail morphology observation. For gene expression analysis, three biological replicates, each containing 10–30 animal tails, 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 ${}^{\circ }$C.

The detailed method for RNA-Seq and SMRT sequencing, and functional annotation of transcripts were described previously [7]. All raw data from Illumina sequencing have been deposited in the NCBI Short Read Archive (SRA) database with the accession number PRJNA504611.

To compare the transcript levels in the tail among the four developmental stages, FPKM (fragments per kilobase of exon model per million mapped reads) was used to evaluate the transcript levels among the four stages. Then, we used DESeq R package (1.10.1) which is available at https://www.huber.embl.de/users/anders/DESeq/ to identify DETs. The resulting DETs were subjected to GO and KEGG pathway enrichment analyses. The details for these analyses were as reported before [7].

We studied two phases of tail development: tailbud determination and tailbud outgrowth, in M. fissipes (Fig. 1, Ref. [9]). Tailbud determination from S18 (neural tube) to S19 (elongated tailbud) is defined by tailbud elongation while tailbud outgrowth refers to stages from S19 to S28, when a functional tail is formed [5]. At S21, tail fin and medial cloacal tail piece were distinct and the tail starts to outgrow. S28 is the end of embryogenesis, when the tadpole swims freely and start to feed, while at S18–S21, the tadpoles hardly move. Tail development is a fast process in M. fissipes, needing only 15 hours from fertilization to neural tube (tail determination) and 3.5 days to the end of embryogenesis (S28). At each sampling stage, the tailbud was cut from posterior to cloaca as indicated by the red dash line in Fig. 1 and used for RNA isolation.

Fig. 1.

Morphological changes associated with tail development of M. fissipes embryos/tadpoles at different stages used in our study. S18: the neural tube stage when tail bud is small but visible; S19: tail bud I; S21: tail bud III when tail fin and medial cloacal tail piece are distinct; S28: end of embryogenesis and onset of feeding (see more detail at [9]). Red dashed line marks the point of dissection for tissue isolation. The time (hr) from fertilization to each of the four stages are also indicated.

Total RNA used for SMRT sequencing was isolated from heart, liver, spleen, lung, kidney, skin, ovary, and testis from two adults M. fissipes as well as from tadpole tail including at four stages (S38, S40, S41, S43) and tadpole dorsal muscle at three stages (S36, S43, S45). The RNA was sequenced on a PacBio RSII platform. All high-quality short reads from the RNA-Seq from the tadpole tails at eight developmental stages (S18, S19, S21, S28, S39, S40, S41, S43) were used to error-correct the SMRT sequences as previously described [7] and analyzed by using the Proovread correction software [10]. We found that de novo transcripts in the 200–500 bp range account for the most, 60.11% (360,893) of the total reads (Table 1) while only 0.22% (149) of that in SMRT sequencing data [7]. However, 30.78% of full-length transcripts had the length $>$2000 bp while only 3.85% of that in de novo transcripts (Table 1).

Through basic alignments to public databases, in total, 132,414 (22.05%), 82,804 (13.79%), 46,184 (7.69%), 116,599 (19.42%), 118,544 (19.74%), 121,829 (20.29%) and 42,249 (7.04%) de novo transcripts were annotated in the NR, NT, KEGG, SwissProt, PFAM, GO and KOG databases, respectively (Table 1). However, the annotation rate of de novo transcripts in each database was much lower than that of SMRT transcripts, such as 47,112 (69%) and 41,254 (61%) in the NR and GO, respectively [7]. Totally, there are only 32.67% of de novo transcripts were annotated in at least one database while the annotation rate of full-length transcripts were 78% [7]. These data indicated that transcripts obtained from SMRT sequencing were highly improved compared to de novo RNA-seq data.

In order to identify the gene expression program underlying tail development, we mapped RNA-Seq clean reads at four indicated stages (S18, S19, S21, S28) to proofread-corrected reference sequences. Totally 329,650,058 (50.11%) reads were mapped (Supplementary Table 1). Then, the expression level was determined by using RSEM software [11]. 48,826 expressed transcripts with FPKM $>$0.3 were identified in the tail among the four stages. Among of them, 22,208 transcripts were commonly expressed at all four stages while 2313, 854, 2576 and 5305 transcripts were uniquely expressed at S18, S19, S21 and S28, respectively (Fig. 2). Next, pair-wise comparisons of gene expression among the four stages were carried out, revealing 8807 DETs in the tail during tailbud determination and outgrowth (Fig. 3, Supplementary Table 2). Among the six pair-wise groups, most DETs were found when any of the other three stages were compared with S28 while relatively few DETs were found when comparisons were done between any two stages among S18, S19 and S21 (Fig. 3). Overall, we identified 4284 DETs (3695 upregulated, 589 downregulated), 5311 DETs (4075 upregulated, 1236 downregulated) and 3093 DETs (2624 upregulated, 469 downregulated) for comparisons between S28 vs S18, S28 vs S19 and S28 vs S21, respectively. These results indicated that most of the DETs were upregulated at S28 comparing to other three stages during tail development.

Fig. 2.

Venn diagram analyses reveal stage-specific of gene expression. Venn diagrams showing the number of expressed transcripts (FPKM $>$0.3) at different stages during tail development. Note that totally 48,826 transcripts were expressed, among which 22,208 transcripts were expressed at all four stages, and that each stage had some uniquely expressed transcripts, with S28 having the most (5305). FPKM: fragments per kilobase of exon model per million mapped reads.

Fig. 3.

Pair-wise comparisons reveal that S28 tail has the highest numbers of differentially expression transcripts (DETs, q 0.05) during tail development. Volcano plot of DETs in the tail for six different comparison groups. X axis represents fold-change (log${}_2$FC) of the DETs and the Y axis represents the –log10 padj- (adjusted P- or q-) value of the DET, with 0.05 set as the significance level cutoff. Up-regulated DETs are shown in red while down-regulated ones in green; non-significantly DETs are presented in blue.

Hierarchical clustering of totally 8807 DETs were performed and the results were shown in a heatmap of the expression levels of each DET at four stages in Fig. 4. The DETs fell largely into two groups with highest or lowest expression levels at S28, and there were more DETs in the group with the highest expression level at S28 (Fig. 4). Thus, tail maturation at S28 involves more DETs than the other developmental stages and that most of the DETs involved in tail maturations are upregulated by S28. In addition, the tail gene expression patterns at S18, S19, and S21 were very similar to each other, especially between S18 and S19, consistent with the lack of significant morphology changes in the tail between S18 to S19 but drastic changes at S28 (Fig. 1).

Fig. 4.

A heatmap showing the global expression patterns of the 8807 DETs at the four stages. The intensity of color indicates relative expression levels. Colors from red to blue correspond to high to low levels of expression. Note that the DETs fall into largely two categories with highest or lowest expression at stage 28, respectively. In addition, the expression levels for the majority of the DETs are most different between S18/19 and S28.

To understand the global gene functional categories and biological pathways during tail development, GO and KEGG analyses were performed on the 8807 DETs. These analyses revealed 110 significantly enriched GO categories and 6 highly enriched KEGG pathways (q 0.05) (Supplementary Fig. 1, Supplementary Table 3). Interestingly, 5 of the 6 enriched KEGG pathways were different metabolic processes while the last one was ECM-receptor interaction, highlighting the importance of metabolism during tailbud growth period. In support of this, the most significantly enriched GO terms were related to oxidoreductase activity and metabolism/catabolism.

To correlate the gene expression changes with the two phase of tail development, we first used K-means algorithm to cluster the DETs based on gene expression profiles across the four developmental stages. The DETs were scattered into eighteen clusters (Supplementary Fig. 2). To explore the gene regulation programs underlying different developmental stages, we focused on cluster 15, cluster 13, and cluster 3, in which the expression of the DETs were upregulated to high levels by S19, S21 and S28, respectively, suggesting their involvement in the respective transitions in the developmental stages. The DETs in each of the three clusters were classified into three GO categories by the Blast2GO, including biological process (BP), cellular component (CC), and molecular function (MF), and directed acyclic graph (DAG) was used to display the GO structures (Fig. 5). Top 10 GO term levels were selected as the master nodes in which only highly enriched GO term levels, such as levels 1, 6, 7, 8, 9, and 10, were displayed, with no significantly enriched GO terms were found in the GO term levels 2–5 (Fig. 5A, note that a node represents a GO term, an arrow represented a parent node, and the branches represent the containment relationships, with the range of function decreasing in size from top to bottom). For cluster 15, only DETs within the BP category were significantly enriched and the DAG showed that all enriched GO terms were eventually linked to three GO nodes: leucine catabolic process, valine catabolic process and tryptophan metabolic process (level 10 GO terms, Fig. 5A), suggesting an important role of amino acid metabolism in the tail development at S19 and later.

Fig. 5.

Directed acyclic graph (DAG) displaying the results of DEG enrichment according to GO analysis in three different DEG clusters. “N” in the left panel represents the number of DEGs in each cluster. For each cluster, the significantly enriched GO terms were arranged according to the GO levels, with the first level being either BP (Biological Process), CC (Cellar component), or MF (Molecular Function). Each box and circle indicate a GO term and the descriptions inside of which, from top to bottom, show GO term id, GO description, q value and the enriched DETs numbers/background gene numbers in the indicated GO term, respectively. Color depth represents the enrichment degree (based on q value, red is the most significantly enriched). Note that some levels had no GO terms shown as there were no GO terms enriched among the DEGs. (A) A cluster where the DETs were upregulated by S19 and remained highly expressed subsequently. These DETs were only significant enriched in biological process categories with leucine catabolic process, valine catabolic process and tryptophan metabolic process as the most significantly enriched. (B) A cluster where the DETs were upregulated by S21 and remained highly expressed subsequently. The GO terms “cytoskeletal anchoring at plasma membrane” of BP and “dystrophin-associated glycoprotein complex” of CC, and “hyaluronic acid binding” of MF are the most significantly enriched. (C) A cluster where the DETs were upregulated by S28. The GO terms “potassium ion transport” of BP and “myosin complex” of CC, and “potassium channel activity” of MF are the most significantly enriched.

For cluster 13, we found significantly enriched GO terms in three GO categories, BP, CC and MF and that there was a single node within each GO category linked to all enriched GO terms (Fig. 5B). These single nodes were GO term “cytoskeletal anchoring at plasma membrane” within BP, “dystrophin-associated glycoprotein complex” within CC, and “hyaluronic acid binding” within MF (Fig. 5B), indicating an important role of these process for the transition from tailbud determination to tailbud outgrowth (S19–S21).

Finally, the DETs upregulated by S28 (cluster 3) also had significantly enriched GO terms in all three GO categories, BP, CC and MF. There was a single node, the GO term “potassium ion (K${}^{+}$) transport” within BP. There were two nodes, the GO terms “myosin complex” and “intermediate filament cytoskeleton” within CC, and two nodes, the GO terms “potassium channel activity” and “fructose 1,6-bisphosphate 1-phosphatase activity” within MF, respectively (Fig. 5C). These processes are likely important for the function of the tail as embryogenesis is complete by S28. Interestingly, a K${}^{+}$-related GO term was a node within both BP and CC, implicating a role of K${}^{+}$-dependent processes in tail maturation and/or function. Notably, the most significant difference between S28 and S18–21 is that tail can move freely at S28 but not at S18–21, and muscle activity was reported to associate with potassium displacements. Importantly, muscle K${}^{+}$ homeostasis contributes to muscle fatigue and recovery [12, 13]. It is intriguing that tadpole tail muscle can keep K${}^{+}$ homeostasis under the almost nonstop contracting. Our finding that the genes in potassium ion processes are involved in tail development may have important implications for understanding muscle fatigue in human.

When we carried out KEGG pathway analysis on these three clusters (${}^{\mathrm{\#}}$15, 13, 3), we found 2, 7 and 11 KEGG pathways highly enriched (q 0.05) for the three clusters, respectively (Supplementary Fig. 3, Supplementary Table 4; Supplementary Fig. 4, Supplementary Table 5;Supplementary Fig. 5, Supplementary Table 6). The two KEGG pathways for cluster 15 were “Pyruvate metabolism” and “Citrate cycle (TCA cycle)” (Supplementary Fig. 3, Supplementary Table 4), both of which are related to the energy metabolism, likely important during tailbud outgrowth because energy metabolism is essential for the biosynthesis of amino acids, nucleotide and lipid during development. Both clusters 13 and 3 highly enriched (q 0.05) KEGG pathways “ECM-receptor interaction” and “Focal adhesion”, implicating a critical involvement of cell-cell and cell-EMC interactions during late stages of tail development as genes in these clusters were upregulated by stage 21 and 28, respectively.

In contrast to the above three clusters, in which the DETs were upregulated to highest levels at S28, two clusters 4 and 11, with a combined 1485 DETs, were DETs whose expression was downregulated to lowest levels by S28 (Supplementary Fig. 2). Analyses of these 1485 DETs revealed the enrichment of 31 GO terms (Supplementary Fig. 6A, Supplementary Table 7) and 9 KEGG pathways (Supplementary Fig. 6B, Supplementary Table 7) with q 0.05. Within the GO category Biological Process, there were 20 organ development GO terms highly enriched (Supplementary Table 7). To further investigate the potential function of the DETs involved in tail development, we identified the common DETs present in each of these 20 GO terms and found 29 DETs, including Xld3, HDAC1, dkk1, 14-3-3, Ruvb-like2, Foxc1, Xiro3, Arih2, Hes1, etc. (Fig. 6, Supplementary Table 8). All these genes were highly expressed at S18 and/or S19, but were gradually repressed from S21 to S28, suggesting that these genes are important for early tail development, i.e., tail determination period (Fig. 6), but not for tail outgrowth or that their repression is important for tail outgrowth.

Fig. 6.

A heatmap showing the developmental expression patterns of the 29 EDTs present in total 20 organ development GO terms enriched among genes downregulated between stage 18 and stage 28. Most of the genes are highly expressed at S18 and S19 while downregulated by S28. The colors red to blue corresponds to high to low levels of expressions.

The KEGG pathways enriched among these 29 downregulated DETs included TGF-$\beta$ signaling pathway, signaling pathways regulating pluripotency of stem cells and Notch signaling pathway. Common among these pathways was the regulation of smad1/5/8, Id, Hes1/5 and HDAC genes.

It is of interest that of the 9 highly enriched (q 0.05) KEGG pathways among downregulated DETs in clusters 4 and 11 were Cell cycle and Hippo signaling pathway (Supplementary Fig. 6B, Supplementary Table 7) (Fig. 7). Both Cell cycle and Hippo pathways are very important for the organ development. The downregulation of genes in these pathways, including Skp2, CDK2-cyclinA, Mps1, HDAC, Cdh1 and Yap/TAZ, Slug, etc. (Fig. 7), suggest that they are involved in tailbud determination and perhaps early tailbud outgrowth but are no longer needed as tail development approaches completion by S28.

Fig. 7.

Selected pathways enriched among DETs that are highly expressed between S18 and S19 but downregulated by S28. (A) Cell cycle pathway. Note that Skp2, CDK2-cyclinA, HDAC, Mps1 and Cdh1, etc., were among the regulated genes. (B) Hippo signaling pathway. Note that BMPs, Wnt and YAP/TAZ were among the regulated genes. The red color indicates the genes that significantly enriched of the pathways.