Skip to main content
Download PDF

Identification and differential regulation of microRNAs during thyroid hormone-dependent metamorphosis in Microhyla fissipes

BMC Genomics volume 19, Article number: 507 (2018) Cite this article

Abstract

Background

Anuran metamorphosis, which is obligatorily initiated and sustained by thyroid hormone (TH), is a dramatic example of extensive morphological, biochemical and cellular changes occurring during post-embryonic development. Thus, it provides an ideal model to understand the actions of the hormone and molecular mechanisms underlying these developmental and apoptotic processes. In addition to transcriptional factors, microRNAs (miRNAs) play key roles in diverse biological processes via post-transcriptional repression of mRNAs. However, the possible role of miRNAs in anuran metamorphosis is not well understood. Screening and identification of TH-responding miRNAs are required to reveal the integrated regulatory mechanisms of TH during metamorphosis. Given the specific role of TRs during M. fissipes metamorphosis and the characteristics of M. fissipes as an ideal model, Illumina sequencing technology was employed to get a full scope of miRNA in M. fissipes metamorphosis treated by T3.

Results

Morphological and histological analysis revealed that 24 h T3 treatment M. fissipes tadpoles resembled that at the climax of natural metamorphosis. Thus, small RNA libraries were constructed from control and 24 h T3 treatment groups. A total of 164 conserved miRNAs and 36 predicted novel miRNAs were characterized. Furthermore, 5′ first and ninth nucleotides of miRNAs were significantly enriched in U in our study. In all, 21 miRNAs were differentially expressed between the T3 and control groups (p < 0.01). A total of 10,206 unigenes were identified as target genes of these differentially expressed miRNAs. KEGG pathway analysis indicated that the most overrepresented miRNA target genes were enriched in the “PI3k-Akt signaling pathway”. In addition, a network associated with the TH signaling pathway provides an opportunity to further understand the complex biological processes that occur in metamorphosis.

Conclusions

We identified a large number of miRNAs during M. fissipes metamorphosis, and 21 of them were differentially expressed in the two groups that represented two different metamorphic stages. These miRNAs may play important roles during metamorphosis. The study gives us clues for further studies of the mechanisms of anuran metamorphosis and provides a model to study the mechanism of TH-affected biological processes in humans.

Background

Amphibian metamorphosis is a dynamic post-embryonic development phase that adapts aquatic to terrestrial life. Thyroid hormone (TH: T3 and T4) induces the complete metamorphosis of anuran tadpoles into juvenile frogs, which resembles mammalian postembryonic development around birth in the presence of high levels of TH [1]. Morphological changes induced by TH cover organogenesis such as limb formation, complete organ degeneration such as gill and tail degeneration [2], and comprehensive organ remodeling such as intestine, pancreas, liver, dorsal muscles and brain remodeling [3, 4]. Whereas it is difficult to manipulate uterus-enclosed mammalian embryos, tadpoles are easily manipulated and independent of any maternal influence [5]. Therefore, anuran metamorphosis regulated by TH provides an ideal model for studying the mechanisms underlying hormone-responsive morphological and functional remodeling, tissue-specific responses to hormone, the precise timing of morphological changes, and the development of adult organ-specific stem cells in vertebrates. Furthermore, anuran postembryonic development, which is triggered and orchestrated by TH, serves as a powerful model to identify the effects of environmental chemical contaminants on TH signaling in vertebrates.

Metamorphosis can be divided into three periods based on the sequence of developmental processes: premetamorphosis, where tadpoles grow and develop hindlimb buds; prometamorphosis, where hindlimbs develop and digits differentiate; and metamorphic climax, where forelimbs emerge and the tail is resorbed completely [6]. A significant number of researches have been devoted to study the mechanisms of anuran metamorphosis induced by THs. THs participate in cellular signal transduction by acting as the ligands of thyroid hormone receptors (TRs) [3]. When TH is absent, the unliganded TR/9-cis-retinoic acid receptor (TR/RXR) heterodimers recruit the co-repressors complex that has the histone deacetylase (HDAC) 3 activity [7]. Histone deacetylation leads to a closed chromatin conformation that is inaccessible for transcriptional machinery, thus repressing the expression of downstream genes. When liganded with TH, a conformational change in the TR induces it to release the co-repressor complex and recruit an SRC/p300 co-activator complex that contains histone acetyl transferase activity [6, 8]. These histone modifications result in chromatin opening and the activation of direct TH response genes expression. TRs can also recruit co-activator complexes such as the SWI/SNF complex, which is involved in chromatin remodeling and the Mediator complex directly involved in transcription activation [9]. In turn, the products of these direct target genes affect the expression of downstream genes. A large number of genes and pathways activated by TH in the larval amphibian during metamorphosis have been identified [10,11,12,13,14,15]. Particularly, TH directly upregulates the expression of Shh and induces bone morphogenetic protein-4 (BMP-4) expression in the Xenopus laevis intestine during metamorphosis, which indicates that the Shh/BMP-4 signaling pathway plays key roles in amphibian intestinal remodeling [16]. Notch signaling plays a role in stem cell development by regulating the expression of Hairy genes during intestinal remodeling in metamorphosis [5]. Two key members of the mitochondrial death pathway, XR11 and caspase-9, are found to regulate the proliferative status of neuroblasts in the Xenopus brain during metamorphosis [17]. Therefore, TH alters a number of genes and regulatory pathways leading to morphological changes during metamorphosis.

In addition to transcriptional factors, microRNAs (miRNAs) play key roles in the regulation of post-transcriptional gene expression during diverse biological processes such as cell proliferation, programmed cell death, fat metabolism, organ development, tumorigenesis and behavior [18]. MiRNAs are endogenous, approximately 22 nucleotides (nt), small non-coding RNA molecules that regulate gene expression via complementary binding to the 3′ untranslated regions (UTRs) of target messenger RNAs (mRNAs) and causing mRNA cleavage or translation blockage [19, 20]. Previous studies have shown that miRNAs are involved in Xenopus oocyte, egg, early embryo, liver, and skin development [20,21,22,23]. Given that miRNAs play key roles in regulating gene expression during anuran development, identifying miRNAs and their target genes that are responsible for amphibian metamorphosis is particularly meaningful. Although a series of genes, signaling pathways and even hormones related to metamorphosis have been identified previously, miRNAs involved in anuran metamorphosis, as well as their potential functional networks, are still less reported. Only 75 unique miRNAs belonging to 30 conserved families in Xenopus tropicalis have been identified during anuran metamorphosis so far [24]. Although some effects of miRNAs on metamorphosis have been clarified in these studies, the specific functions of miRNAs and the complicated molecular regulatory networks they may participate in remain unclear. Therefore, our limited knowledge of miRNAs related to anuran metamorphosis prevents us from having a systematic and common understanding of the miRNA:mRNA regulation of anuran metamorphosis.

Microhyla fissipes is a small-sized and widely distributed anuran in southeastern Asia and eastern Asia. Due to its fast development, transparent tadpole phase, biphasic life cycle, diploid, and embryo development in vivo, M. fissipes is a good model to study tissue development and apoptosis, environmental toxicology, and adaptive mechanisms from aquatic to terrestrial life [25].

Given the specific role of TRs during M. fissipes metamorphosis and the characteristics of M. fissipes as an ideal model [6, 25], Illumina sequencing technology was employed to obtain a full scope of miRNAs in M. fissipes metamorphosis induced by T3. From these sequencing data and the subsequent miRNA annotation, we identified differentially expressed miRNAs that may play important roles in regulating metamorphosis initiation and the cellular process during metamorphosis. The identification and characterization of known and novel miRNAs will enable us to better understand the roles of these miRNAs in anuran metamorphosis. Subsequently, we predicted the mRNA targets of M. fissipes miRNAs and correlated their expression with the mRNA expression identified in a previous study [6]. The integrated analysis of miRNA and mRNA expression profiles will explore their possible regulatory patterns at the critical stage of anuran development. This study will advance our understanding of how post-transcriptional regulators may regulate TH dependent metamorphosis and will provide a model to study the mechanisms how TH affects such biological processes in humans.

Methods

Tadpoles

Mature female and male M. fissipes were collected from Shuangliu, Chengdu, China (30.5825o N, 103.8438o E) in June, 2016. The care and treatment of animals in this study were performed according to the Guideline for the Care and Use of Laboratory Animals in China, and the animal experiments were approved by the Experimental Animal Use Ethics Committee of the Chengdu Institute of Biology (Permit Number: 2016036). After acclimatized in our laboratory for 1 week, the male and female were injected LHRHa with 0.3 μg/g body weight resolving dosage. Embryos were obtained from one pair of frogs and were subjected to a 12:12 h light:dark cycle at 25 ± 0.6 °C. Tadpole developmental stage was recorded using the M. fissipes developmental table [26].

Sample collection

Tadpoles at stage 33 (premetamorphosis, oar-shaped limb bud) were exposed to 10 nM 3,3′,5-triiodo-l-thyronine sodium (T3, Sigma). The chemicals were renewed every 12 h of exposure when the medium was refreshed. At each exposure time point (0 h, 12 h, 24 h, 36 h, and 48 h), six tadpoles were randomly collected for morphological and histological observations. Tadpoles treated for 0 h were set as the control group. For the small RNA library construction, samplings of the control group (Control-1, Control-2, and Control-3) and the T3 exposure group (24 h) (T3–1, T3–2, and T3–3) had three biological replicates, respectively, and each replicate was made up of one individual. Tadpoles in this study were anesthetized in 0.01% MS222 (ethyl 3-aminobenzoate methanesulfonate) and were then processed as described below. Each individual from the control and T3 groups for small RNA sequencing was immediately immersed in liquid nitrogen and stored at − 80 °C until the RNA was extracted (within 3 days).

Morphological and histological observations of tadpoles

A stereo microscope (JSZ8T, Jiang Nan Yong Xin, Nanjing, China) with the Mshot Image Analysis system (Mc50-N) was used to observe and take photos of tadpoles and to measure their total length (TOL), tail length (TL), body width (BW), hindlimb length (HLL), hindlimb width (HW), snout length (SL), interocular space (IOS) and mouth width (MW) (Fig. 1). After external morphological observations, intestines were obtained for measuring intestine length (IL). Then, the small intestines were dissected and were fixed in Bouin’s Fluid. After dehydration in a graded series of ethanol and transparency by xylene, intestines were embedded in paraffin and were sectioned in serial transverse sections (7 μm thick) using a YD-1508R rotary microtome (Huaian Kehuai Instrument Ltd., Jiangsu, China). Dewaxed serial sections were stained with Delafield’s hematoxylin and were counterstained with eosin (HE) to show general histological characteristics.

Fig. 1
figure 1

a-e Morphological changes of M. fissipes tadpoles induced by T3. a1-a3: control; b1-b3: T3–12 h; c1-c3: T3–24 h; d1-d3: T3–36 h; e1-e3: T3–48 h. Scale bar = 1 mm. f Morphological parameter of M. fissipes tadpoles. TOL: total length; TL: tail length; BW: body width; HLL: hind limb length; HW: hindlimb width; SL: snout length; IOS: interocular space. g Morphological parameter of M. fissipes tadpoles treated with T3 at different time points. Bars represent the mean ± SD (n = 6). The data were analyzed by One-Way ANOVA and significance is set at p < 0.05. The asterisks indicate a significant difference compared to control, *p < 0.05, **p < 0.01. IL: intestine length

Full size image

The levels of significance in the morphological characteristics (TOL, BW, TL, HLL, HW, IL, SL, IOS and MW) and gene expression were determined by one-way analysis of variance (ANOVA) with SPSS Statistics 13.0 (SPSS Inc., Chicago, IL, USA); a value of p < 0.05 was considered to be statistically significant.

Small RNA sequencing and annotation

Total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The total RNA quantity and purity were analyzed with Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) and RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA) with RIN number > 8.0.

After sRNAs with 15–35 nt in length were isolated from 1 μg total RNA by size fractionation in a 15% TBE urea polyacrylamide gel, the purified sRNAs were then ligated with both 3′ and 5′ adaptors for reverse transcription. After PCR amplification with forward and reverse primers complementary to the 5′ and 3′ adapters and purification with the proper size, the small RNA libraries were constructed. Subsequently, the qualified small RNA libraries were sequenced with the HiSeq 2500 platform, with 50 base-pair end-read sequencing runs (Basebio Tech, Chengdu, China). The raw data were deposited in the NCBI Sequence Read Archive (SRA) database under the study accession SRR5451121, SRR5451123, SRR5451131, SRR5451133, SRR5451135, and SRR5451136.

The raw data were processed by summarizing data production, evaluating sequencing quality, filtering low quality reads, 5′ adapter contaminants, reads without 3′ adapter, reads without insert fragment, reads containing poly (A) stretches, and reads less than 15 nt or longer than 35 bp. Length distribution analysis was performed on small RNA reads. For cataloging clean reads in terms of genome annotation to separate out rRNA, tRNA, snRNA and snoRNA, the clean RNAs were mapped to the X. tropicalis genome (http://www.ensembl.org), one of M. fissipes closely related species, by Novoalign with the maximum penalty score of alignment = 60. Subsequently, the remnant RNA sequences were re-queried against the miRNA database, miRBase 21.0 (ftp://mirbase.org/pub/mirbase/21) by using Novoalign to the further identification of conserved miRNAs. The remaining sequences that did not match known miRNAs were mapped to the X. tropicalis genome to identify potentially novel miRNAs by using miRDeep2. Novel miRNAs were predicted if the extended sequences at the mapped positions have the propensity of forming hairpin structures. Subsequently, we analyzed the base bias on each position of all identified conserved miRNAs in our study.

Analysis of differentially expressed miRNAs

To determine the abundance of miRNAs and the differentially expressed miRNAs, the read counts of miRNAs were transformed into TPM (transcript per million) through the following criteria: Normalization formula: Normalized expression = (Actual miRNA count/Total count of clean reads)*1,000,000 [27]. The fold-change and p-value (p adjust) were calculated from the normalized expression using the DESeq R packages. The p value < 0.01 and |log2 (foldchange)| > 1 were set as the threshold for significantly differential expression.

Target prediction and function analysis

The available RNA transcriptome data (GECV01000000) [6] was downloaded and used for the prediction of potential target mRNA candidates for significantly differently expressed miRNAs between the control and T3 groups using miRanda 3.3a [28]. To determine the function of target genes, Gene Ontology (GO) (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (http://www.genome.jp/kegg/) enrichment were conducted. GO annotation based on categories of biological processes, molecular functions and cellular components was performed by Blast2GO 3.1 analysis [29]. KOBAS 2.0 (http://kobas.cbi.pku.edu.cn) [30] was used to test the statistical enrichment of the target gene candidates in the KEGG pathways. In addition, because of the importance of TH signaling in anuran metamorphosis, the interaction analysis of miRNA-mRNA pairs were performed on this pathway based on the targets prediction and the expression levels of mRNAs and miRNAs. As we got the mRNA data from the previous study [6], we performed the pairwise comparison (MC vs. PM) using the DEGseq R package for screening differentially expressed genes. Unigenes were considered as differentially expressed genes at a normalized fold change > ± 2 at p value < 0.005 after adjustment for the false discovery rate. Afterwards, the TH pathway was downloaded from the KEGG and integration of miRNA-seq with mRNA-seq was achieved by integrating the expression profiles of differently expressed miRNAs and differently expressed mRNAs in the TH pathway with the addition of differently expressed miRNA-targeting information.

Quantitative reverse transcriptase PCR for miRNAs

Quantitative stem-loop RT-PCR (qRT-PCR) was performed to profile the relative expression of 8 miRNAs, including 7 known and 1 novel miRNAs in the control and T3 group. Total RNA was isolated using the same reagents as described above. The expression of miRNA was detected using U6 snRNA as an endogenous control. All the primers used in the qRT-PCR are listed in Additional file 1: Table S1. The relative amount of miRNA to U6 was calculated using the 2−ΔΔct method, and all quantitative data presented were the mean ± SEM.

Results

Morphological and histological observations in T3-induced metamorphosis

In the presence of exogenous T3, the morphology of M. fissipes tadpoles changed dramatically, similar to the natural metamorphosis. For example, after 48 h of T3 treatment, the tail shortened; the hindlimb buds differentiated and formed the toes; angiogenesis was observed in the hindlimb; the skin remodeled from transparent to black with pigment; the mouth became bigger; and the eyes have moved ahead (Fig. 1a). TOL, BW, TL, SL, and IOS increased after exposure to T3 at 12 h and then decreased from 24 h (Fig. 1c). Meanwhile, HLL and HW increased during T3-induced metamorphosis. MW decreased from 0 h to 36 h and then increased from 36 h to 48 h. Furthermore, TOL significantly decreased by 9.98% from 12 h to 24 h (p < 0.05). TOL and TL significantly decreased by 21.27 and 21.84% (p < 0.01) from 24 h to 36 h. Significant decreases were detected in SL (p < 0.01), BW (p < 0.05) and IOS (p < 0.01) at 48 h. Thus, the altered trends of these morphological characteristics after T3 treatment were similar to the natural metamorphosis (Additional file 2: Figure S1), and tadpoles exposed to 24 h of T3 treatment resembled the tadpoles at the climax of metamorphosis.

In addition to external morphological changes in response to T3, we examined internal metamorphic changes. The intestine is one of the best-studied organs that remodels during metamorphosis including apoptosis of larval epithelial cells and proliferation of adult precursor epithelial cells [31]. Additionally, T3 treatment of premetamorphic X. laevis and X. tropicalis tadpoles are known to induce similar changes in the intestine to those during natural metamorphosis [32]. In our study, IL significantly decreased by 33.37, 59.86 and 72.84% at 12 h (p < 0.05), 36 h (p < 0.01) and 48 h (p < 0.01), respectively. To further examine the changes in the intestine during metamorphosis, the anterior small intestines from M. fissipes tadpoles after T3 treatment were cross-sectioned and HE stained. The premetamorphic (S33) intestine was a simple tube and contained mostly a monolayer of larval epithelial cells, with thin layers of muscles (Fig. 2a). After 12 h of T3 treatment, metamorphosis had begun, and the thin larval muscle in the intestine had begun to increase in thickness (Fig. 2b). At 24 h, larval epithelial cell death and adult epithelial cell proliferation had begun, and there was also an evident increase in the muscle and connective tissue of the intestine (Fig. 2c). At 36 h, there were numerous proliferating adult epithelial cells, apoptotic larval epithelial cells, and connective tissues (Fig. 2d). At 48 h, the epithelium developed into the multiply folded adult structure, and the connective tissue and outer muscle were abundant (Fig. 2e). T3-induced intestinal remodeling mimicked its remodeling observed during natural metamorphosis, as reported in X. trpocalis [32]. In particular, the M. fissipes tadpoles that were treated for 24 h with T3 resembled the natural tadpoles at the climax of metamorphosis. Thus, morphological and histological analysis revealed that M. fissipes tadpoles treated for 24 h with T3 would be selected for further miRNA-seq.

Fig. 2
figure 2

Histological changes of M. fissipes intestinal remodeling induced by T3. a control; b T3–12 h; c T3–24 h; d T3–36 h; e T3–48 h. Ct: connective tissue; Ep: epithelium; Mu: muscle. The red arrow indicates the apoptosis larval epithelial cells, while the blue arrow indicates the proliferating adult epithelial cells. Scale bar = 10 μm

Full size image

Overview of small RNAome in M. fissipes metamorphosis

To survey miRNAs involved in M. fissipes metamorphosis, six small RNAs libraries from the control (Control-1, Control-2, Control-3) and 24 h T3 treatment (T3–1, T3–2, and T3–3) groups were constructed and sequenced. After removal of the adapters, low-quality reads, reads with unknown nucleotides, reads smaller than 15 nucleotides, reads larger than 35 nucleotides, reads without a 3′ adaptor, and reads with polyA/T, 17.2, 19.1, 14.2, 18.6, 15.3, and 18.8 million clean reads were generated in the Control-1, Control-2, Control-3, T3–1, T3–2, and T3–3 libraries, respectively (Additional file 3: Table S2). An overview of the reads for sRNA-seq from the raw data to high-quality data and with quality filtering is shown in Additional file 3: Table S2.

The length distributions of sRNAs were similar between two groups in which 22 nt RNAs were the most abundant (Additional file 4: Figure S2a). After genome annotation to separate out rRNA, tRNA, snRNA, snoRNA and mRNA, more than 6 million sequences corresponding to 70.11 and 65.40% of all clean reads were annotated to known miRNAs in the control and T3 treatment group, respectively (Additional file 5: Figure S3). These results showed that our sRNA libraries were highly enriched with mature miRNAs.

Identification of conserved and novel miRNAs in M. fissipes metamorphosis

Many miRNAs have been conserved during evolution, and the sequence identity is extremely high even among evolutionarily distant species [33, 34]. To investigate the conserved miRNAs in M. fissipes metamorphosis, the sequences of the miRNAs from the six libraries were compared to X. tropicalis miRNA sequences in miRBase 21.0. In total, 164 conserved miRNAs belonging to 82 miRNA families in the libraries (Additional file 6: Table S3) were identified. The details families of the member numbers of conserved miRNAs are summarized in Additional file 6: Table S3. A total of 34 conserved miRNA families contained more than one member.

The conserved miRNAs have a broad range of expression levels, varying from 1 count to 978,287 counts, confirming the thoroughness of our miRNA screen strategy. For T3 and control libraries, the majority of unique sequences for conserved miRNAs exhibited high abundance (> 1000 counts), although a few unique sequences were present at less than 10 counts (Additional file 4: Figure S2c). Among the miRNAs with high abundance (more than 100,000 counts), 12 miRNAs (mfi-miR-192, mfi-miR-26, mfi-miR-143, mfi-miR-148a, mfi-miR-205a, mfi-miR-22-3p, mfi-miR-181a-5p, mfi-miR-182-5p, mfi-miR-194, mfi-miR-200a, mfi-miR-92a, and mfi-let-7f) were most highly expressed in M. fissipes metamorphosis. In addition, 7 members of the let-7 family (including mfi-let-7a, − 7b, − 7c, −7e, −7f, −7 g, and -7i) were expressed in metamorphosis and their responses to T3 were different.

The mature miR-#-5p and miR-#-3p are pairwise miRNAs that align to the 5′ and 3′ end regions of the same precursors, respectively. In M. fissipes, 15 pairs of known miR-#-5p and miR-#-3p were detected. Among these, 10 pairwise miRNAs (mfi-miR-142, mfi-miR-17, mfi-miR-18a, mfi-miR-20a, mfi-miR-181a, mfi-miR-182, mfi-miR-199a, mfi-miR-30a, mfi-miR-9a, and mfi-miR-9b) were found to have relatively lower expression levels of miR-#-3p than their miR-#-5p counterparts, while the other 5 pairwise miRNAs (miR-29c, miR-22, miR-363, miR-24a, and miR-126) showed relatively higher expression levels of miR-#-3p.

Furthermore, one of the greatest advantages of high-throughput sequencing is that this technology can be used to identify novel miRNAs. In this study, 36 novel miRNAs were predicted. In addition, all novel miRNAs were found in all libraries.

The composition patterns of miRNAs are closely correlated with the biological characteristics of different development stages. The nucleotide bias on each position of all identified miRNAs is shown in Additional file 4: Figure S2b. The majority of miRNAs tended to start with 5’-U but not 5′-G in M. fissipes. Interestingly, the 5′ ninth nucleotide was also significantly enriched in U in our study.

Identification and validation of miRNAs differentially expressed in M. fisssipes metamorphosis

A total of 21 differentially expressed miRNAs were identified between the control group and T3 treated group (|log2foldchange| > 1, p < 0.01) (Additional file 4: Figure S2d), including 20 known miRNAs (mfi-miR-10a, mfi-miR-10b, mfi-miR-133c, mfi-miR-148b, mfi-miR-15a, mfi-miR-15b, mfi-miR-17-5p, mfi-miR-181b, mfi-miR-199a-5p, mfi-miR-216, mfi-miR-222, mfi-miR-25, mfi-miR-301, mfi-miR-30d, mfi-miR-30e, mfi-miR-31a, mfi-miR-363–3p, mfi-miR-429, mfi-miR-9a-5p, and mfi-miR-9b-5p) and one novel miRNA (novel-miR-10).

Among these, the expression levels of 8 and 13 miRNAs were suppressed and promoted by T3 treatment, respectively (Fig. 3a). To validate the expression patterns of miRNAs by deep sequencing, 8 miRNAs with significant expression changes were randomly selected for qRT-PCR analysis. Except novel-miR-10, the expression patterns obtained by qRT-PCR for the control and T3 groups were highly consistent with the sequencing results (Fig. 3b). These results indicated that small-RNA sequencing is a sensitive and reliable method to identify differentially expressed miRNAs in M. fissipes metamorphosis.

Fig. 3
figure 3

Differential expression levels of miRNAs in the control and T3 groups. a Significantly different expression miRNAs in the two groups based on deep sequencing data. b Validation of miRNAs by qRT-PCR, * means a statistically significant difference (p < 0.05) between the T3 and control groups; y-axis means the Fold change with control group as 1

Full size image

Target prediction of differentially expressed miRNAs and function analysis in M. fissipes metamorphosis

The functions of miRNAs on gene regulation mainly rely on the interaction between miRNAs and their target genes via complementation [35]. Therefore, identification of miRNA target genes is essential to understand the regulatory function of miRNA. Target prediction was performed against the unigenes from the M. fissipes transcriptome using miRanda (3.3a). A total of 10,206 unigenes (out of 26,323) were predicted as potential target genes of miRNAs in this study. It indicated that each miRNA always corresponded to hundreds of target genes. For example, mfi-miR-181b could bind to collagen genes, including col5a3, col6a3, and col12a1; ECM proteolysis related genes, including matrix metalloproteinase 2 and matrix metalloproteinase 14; and the hedgehog target gene twist2. All of these genes were involved in natural metamorphosis during the climax of metamorphosis stage for patterning [6]. Additionally, a protein coding gene may also be targeted by multiple miRNAs at different targeting sites. For example, TRα, which is essential for metamorphosis initiation in M. fissipes, was predicted as the target of 7 miRNAs including mfi-miR-10a, mfi-miR-10b, mfi-miR-148b, mfi-miR-15b, mfi-miR-31a, mfi-miR-9a-5p, and mfi-miR-9b-5p. Therefore, many predicted targets were likely to be regulated by multiple miRNAs at different targeting sites, while miRNAs targeted to different genes at different sites.

To obtain an overview of the pathways in which the miRNAs were involved, the predicted target genes for the differentially expressed miRNAs were subjected to GO analysis and KEGG pathway analysis. GO enrichment analysis showed that 216 cellular component ontology terms, 134 molecular function ontology terms and 704 biological process ontology terms were highlighted (Fig. 4). Similarly, all predicted target genes were also blasted against the KEGG database, and 251 pathways were highlighted (p < 0.01), suggesting that each pathway included several miRNAs. The top twenty enriched pathways are shown in Fig. 5. Interestingly, the most overrepresented miRNA target genes were enriched in the “PI3k-Akt signaling pathway”. In addition, RNA transport, ECM-receptor interaction, Wnt signaling pathway, mTOR signaling pathway, ubiquitin mediated proteolysis and lysosome were also significantly enriched.

Fig. 4
figure 4

Gene ontology (GO) classification of differentially expressed target genes

Full size image
Fig. 5
figure 5

KEGG pathway analysis of differently expressed target genes. The top 20 significant enriched pathways based on KEGG analysis. The y-axis shows 20 different enriched pathways in which the target gene productions are involved in. The x-axis shows the number of genes enriched in each pathway

Full size image

Network of putative interactions in the TH pathway

Since a single miRNA can target multiple mRNAs and one mRNA can also be targeted by multiple miRNAs, a context-specific network of miRNA:mRNA interactions can be illuminated. Based on the data from RNA-seq and the importance of the TH signaling pathway in the metamorphosis, differentially expressed miRNAs and their target genes that showed inverse expression patterns were co-mapped in the TH signaling pathway (Fig. 6). MiRNAs and their target genes with opposite counts during metamorphosis were considered potential interactions.

Fig. 6
figure 6

A network of putative interactions between differently expressed miRNAs and target genes in the TH signaling pathway. The red box indicates the genes that were upregulated by T3 treatment, while the green box indicates the genes that were downregulated by T3. The blue line presented the relative expression of miRNAs targeting the gene

Full size image

Among the 10 significantly differential expressed target genes involved in the TH signal pathway, 4 genes (HSP90, PTEN, RXRβ and TRβ) were upregulated from premetamorphosis to the climax of metamorphosis, while TRα, RXRα, type I iodothyronine deiodinase (D1) and transthyretin were downregulated. It was reported that TRα-deficient tadpoles develop faster, however, with a smaller body size at the same stage than their wild-type siblings, suggesting that TRα plays an important role in controlling the timing of metamorphosis [36, 37]. In addition, mfi-miR-10a, mfi-miR-10b, mfi-miR-148b, mfi-miR-9a-5p and mfi-miR-9b-5p, which are predicted to target TRα, were up-regulated. RXRβ and TRβ, which initiated the process of metamorphosis, were up-regulated, while their predicted interaction miRNAs (mfi-miR-133c, mfi-miR-25, mfi-miR-31a and mfi-miR-363–3p; mfi-miR-15a, mfi-miR-15b and mfi-miR-216) were down-regulated. We also obtained the interactions of miRNAs and genes from the PI3k-Akt signaling pathway which has been enriched in the KEGG analysis. The expression of Akt decreased during metamorphosis, while the expression of mfi-miR-9a-5p and mfi-miR-9b-5p, which target Akt increased. Casp9-mfi-miR-181b interaction was also detected, indicating that miRNA regulated anuran metamorphosis by posttranscriptional regulation of genes in PI3k-Akt signaling pathway and in the TH signaling pathway. Therefore, synergistic regulatory pathways among multiple miRNAs and mRNAs play important roles in anuran metamorphosis.

Discussion

In anuran, the transition from larval to juvenile is controlled directly by steadily rising levels of TH. The rising level of TH initiated the transcription of immediate response genes at varying thresholds in each responding tissue and led to extensive morphological and physiological remodeling of tissues. Treated with T3, M. fissipes prematurely accelerated development toward the metamorphic climax similar to that of X. tropicalis [31]. Treatment of M. fissipes tadpoles for 24 h with T3 resulted in the similar of morphological and histological characteristics of natural tadpoles at the climax of metamorphosis. Therefore, the control group and the 24 h T3 treatment group were selected for miRNA-seq to screen and identify TH-responding miRNAs and to reveal the integrated regulation of TH during M. fissipes metamorphosis. Furthermore, due to the near total dependence of anuran metamorphosis on TH, M. fissipes provides an excellent model to develop an assay for TH system disruption across all vertebrates.

Few miRNAs have been reported to be involved in anuran metamorphosis [24]. Using the high-throughput sequencing approach, we identified 164 conserved miRNAs and 36 novel miRNAs and analyzed the expression profiles of these miRNAs during M. fissipes metamorphosis. The most-abundant size class of sRNAs was 22 nt, followed by 23 and 21 nt. This was consistent with the known size of Dicer-processed small RNAs in X. tropicalis somatic libraries and other organisms [20, 38, 39]. The expression levels of all miRNA sequences ranged from 1 to 978,287 counts, indicating that not only miRNAs with high expression but also miRNAs with low expression were detected by Illumina sequencing. Therefore, Illumina sequencing is a more efficient and accurate approach than the traditional cloning method to identify conserved and novel miRNAs.

The expression levels of mfi-miR-192, mfi-miR-26, mfi-miR-143, mfi-miR-148a, mfi-miR-205a, mfi-miR-22-3p, mfi-miR-181a-5p, mfi-miR-182-5p, mfi-miR-194, mfi-miR-200a, mfi-miR-92a, and mfi-let-7f were highest in this study, implying their potential significant functions in M. fissipes metamorphosis. Intriguingly, the number of miR-143a clones was almost one third of the total number of cloned miRNAs during X. tropicalis metamorphosis [24], suggesting its certain conserved role in supporting metamorphosis. MiR-143, miR-181a-5p and miR-22-3p were highly enriched in skeletal muscle tissue [40, 41], which was the major tissue of M. fissipes tadpoles for small RNA-seq. MiR-92 was reported to play a regulated role in the proliferation of myeloid cells [42], while miR-192 suppressed cell proliferation and induced apoptosis. High expression of these two miRNAs in M. fissipes tadpoles further confirmed that their roles in proliferation, differentiation and apoptosis during anuran metamorphosis [40, 41, 43]. MiR-182 was significantly upregulated in anaplastic thyroid cancer (ATC) tissues and cells [44]; its high abundance in this study may suggest its role in thyroid cell activation during metamorphosis. Different expression levels of the let-7 family members implied that the let-7 family likely conducted precise regulation of different target genes, and let-7 miRNAs have been reported to widely participate in tissue development and metabolism during development and the regulation of temporal transitions associated with cell proliferation and differentiation during Japanese flounder metamorphosis [45]. Furthermore, let-7 family miRNAs have indeed been reported to be involved in numerous cellular processes, including as master regulators of the cell cycle [46, 47], indicating that the let-7 family was likely critical to morphological reorganization during metamorphosis by regulating the cell cycle.

Different expression levels between miR-#-5p and their miR-#-3p counterparts were detected in 15 pairs of miRNAs during M. fissipes metamorphosis. Although both strands of the miRNA duplex were equally produced after the transcription of the miRNA genes, their stability in the the single-stranded form was different [48]. As a result, the expression levels of miR-#-3p/5p mainly relied on their respective degradation degrees and degradation rates [27]. Some miR-#-3 ps were reported as mature functional miRNAs with abundant expression, and miR-#-3p/miR-#-5p ratios may vary dramatically in different stages of development [49]. On the other hand, it was found that most of the miR-#-3 ps were detected at the same level or at relatively lower expression levels than miR-#-5p [27, 50, 51]. In this study, more miR-#-5 ps with higher expression levels indicated their importance in regulating gene expression.

The 3′ and 5′ edges of the ‘seed region’, which were known to have a critical role in targeting miRNA to mRNA for translational inhibition or mRNA cleavage were both flanked by U, which has been detected in other vertebrates [46]. This phenomenon suggested that the U base may have effects on miRNA action, such as miRNA biogenesis, mRNA target recognition and miRNA binding to targets for gene regulation. High conservation was also observed at the first position upstream and downstream of the seed match in the target gene. This nucleotide was often a conserved A, which could pair to the nucleotide of the miRNA with 5′ and 3′ edges of its ‘seed region’ were U [52].

In total, 20 known miRNAs (mfi-miR-10a, mfi-miR-10b, mfi-miR-133c, mfi-miR-148b, mfi-miR-15a, mfi-miR-15b, mfi-miR-17-5p, mfi-miR-181b, mfi-miR-199a-5p, mfi-miR-216, mfi-miR-222, mfi-miR-25, mfi-miR-301, mfi-miR-30d, mfi-miR-30e, mfi-miR-31a, mfi-miR-363–3p, mfi-miR-429, mfi-miR-9a-5p and mfi-miR-9b-5p) and one novel miRNA (novel-miR-10) were found to be differentially expressed between the control and T3 groups, indicating that these miRNAs might play important roles in diverse biological processes that resulted in the transition from larval to juvenile M. fissipes. MiR-222 and miR-181b are potential molecular markers used to differentiate papillary thyroid cancer from normal thyroid tissue [53], suggesting their potential roles in regulating thyroid cell division and differentiation during metamorphosis. Interestingly, T3 treatment activated the expression of these two miRNAs in M. fissipes tadpoles. It seems that miR-222 and miR-181b are involved in a positive feedback cycle in amplifying the TH signal to achieve the dramatic metamorphosis. MiR-133 is a muscle-specific miRNA reported to positively regulate differentiation and proliferation of myoblasts [54]. Our study showed that the expression of mfi-miR-133c was suppressed by T3 treatment. Given that replacement and conversion between different types of myofilaments was one of the most dramatic changes during anuran metamorphosis [43], down-regulation of mfi-miR-133c may be necessary to suppress the expression of tadpole-type myofilaments.

The impact of these 21 microRNAs is amplified by a larger number of predicted target genes. Therefore, identification of miRNA targets is an important step to thoroughly describing the function of miRNAs. Further investigation of these differentially expressed miRNAs would be helpful to explain the complex genetic network that controls processes of metamorphosis and to provide insight into metamorphic changes. Based on the principle of high miRNA-target complementarity, 10,206 target genes were predicted. It was currently estimated that the human miRNAs targeted ∼ 60% of all protein-coding genes [55], indicating that each miRNA always corresponded to hundreds of target genes. GO and KEGG pathway analysis showed that the predicted target genes were involved in diverse biological processes, especially cellular operations. These results showed that miRNAs participated in metamorphosis by targeting genes in these pathways. Intriguingly, the most significantly enriched pathway targeted by the differentially expressed miRNAs was the PI3k-Akt signaling pathway, which is the central signaling pathway responding to hormones, extracellular signal molecules and nutrients to regulate numerous cellular events such as metabolic homeostasis, proliferation, growth, and survival [56]. Pathways associated with RNA transport, the Wnt signaling pathway, the mTOR signaling pathway, ubiquitin-mediated proteolysis and lysosome-related pathways were also significantly enriched, indicating the role of differentially expressed miRNAs in the regulation of cell growth, cell proliferation, and apoptosis. Moreover, the enrichment of cell-cell and cell-ECM interactions suggested that the differentially expressed miRNAs participated in communication between cells and the extracellular matrix, and then affected cellular changes, including apoptosis, cell proliferation, migration and differentiation [57]. Ubiquitin -mediated proteolysis and lysosomal degradation occur during metamorphosis, where massive breakdown of larval tissue occurs [58], thus explaining the enrichment of miRNA-targeted genes in the ubiquitin-mediated proteolysis and lysosome degradation pathway. Furthermore, miRNAs enriched in the Wnt signaling pathway are necessary for proper formation of important tissues [59], such as the hindlimb, adult dorsal muscle, adult intestine and adult liver. Therefore, the enrichment results indicated that these differentially expressed miRNAs participated in regulating cell growth, proliferation, communication and apoptosis during M. fissipes metamorphosis.

Based on the integrated miRNA and mRNA analysis together with target prediction and the importance of the TH signaling pathway during metamorphosis, a network of differentially expressed miRNAs and differentially expressed genes in the TH pathway was constructed. The expression of miR-10a was switched on by RARα, a hormone receptor modulated by RXRα [60]. Because RXRα and TRα were upregulated by T3 during metamorphosis and TRα was the target gene of miR-10a, T3-RXRα-RARα-miR-10a-TRα played the negative feedback role in M. fissipes metamorphosis. Elevated expression of TRα and decreased mfi-miR-10a, mfi-miR-10b, mfi-miR-148b, mfi-miR-9a-5p and mfi-miR-9b-5p in the pathway indicated that these miRNAs were essential for the initiation of metamorphosis in M. fissipes. It also implied that TRα is not only the receptor mediated by TH signals but is also a potential target gene regulated by miRNA. It was important for future studies to functionally validate the predicted function of these differentially expressed miRNAs in this pathway.

Conclusion

Our study is the first study on the miRNAs in M. fissipes metamorphosis and the integrated analysis of mRNA-seq and miRNA-seq in amphibians. This study provides an overview of the miRNAs expressed in TH-dependent metamorphosis in M. fissipes through the use of high-throughput Illumina sequencing technology. In total, 164 conserved and 36 predicted novel miRNAs were identified in the M. fissipes metamorphosis. Of those, 21 differentially expressed miRNAs were detected in two metamorphic groups, and functional annotation showed that the differentially expressed miRNAs along with their target genes might be essential to cell growth, cell proliferation, communication and apoptosis. This study provided the basis for future analysis of miRNA function in M. fissipes metamorphic development and gave us clues for further studies of the mechanisms of M. fissipes metamorphosis. Furthermore, using M. fissipes as the model organism, the miRNA detected in this study may provide a marker for the identification of environmental chemical contaminants that perturb TH signaling in vertebrates.

Abbreviations

BW:

body width

GO:

gene ontology

HLL:

hindlimb length

HW:

hindlimb width

IL:

intestine length

IOS:

interocular space

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MW:

mouth width

qRT-PCR:

quantitative reverse transcriptase PCR

RXR:

9-cis-retinoic acid receptor

SL:

snout length

TH:

thyroid hormone

TL:

tail length

TOL:

total length

TR:

thyroid hormone receptors

UTR:

untranslated region

References

  1. Shi YB, Hasebe T, Fu L, Fujimoto K, Ishizuya-Oka A. The development of the adult intestinal stem cells: insights from studies on thyroid hormone-dependent amphibian metamorphosis. Cell & Bioscience. 2011;1(1):30.

    Article  CAS  Google Scholar 

  2. Brown DD, Cai L. Amphibian metamorphosis. Dev Biol. 2007;306(1):20–33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  3. Das B, Cai L, Carter MG, Piao Y-L, Sharov AA, Ko MSH, Brown DD. Gene expression changes at metamorphosis induced by thyroid hormone in Xenopus laevis tadpoles. Dev Biol. 2006;291(2):342–55.

    Article  PubMed  CAS  Google Scholar 

  4. Shi YB. Amphibian metamorphosis: from morphology to molecular biology. New York: Wiley; 1998.

    Google Scholar 

  5. Hasebe T, Fujimoto K, Kajita M, Fu L, Shi Y-B, Ishizuya-Oka A. Thyroid hormone-induced activation of notch signaling is required for adult intestinal stem cell development during Xenopus Laevis metamorphosis. Stem Cells. 2017;35(4):1028–39.

    Article  PubMed  CAS  Google Scholar 

  6. Zhao L, Liu L, Wang S, Wang H, Jiang J. Transcriptome profiles of metamorphosis in the ornamented pygmy frog Microhyla fissipes clarify the functions of thyroid hormone receptors in metamorphosis. Sci Rep. 2016;6:27310.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Shi Y-B, Matsuura K, Fujimoto K, Wen L, Fu L. Thyroid hormone receptor actions on transcription in amphibia: the roles of histone modification and chromatin disruption. Cell & Bioscience. 2012;2(1):42.

    Article  CAS  Google Scholar 

  8. Morvan-Dubois G, Demeneix BA, Sachs LM. Xenopus laevis as a model for studying thyroid hormone signalling: from development to metamorphosis. Mol Cell Endocrinol. 2008;293(1):71–9.

    Article  PubMed  CAS  Google Scholar 

  9. Huang S, Cao X, Tian X, Wang W. High-throughput sequencing identifies micrornas from posterior intestine of loach (Misgurnus anguillicaudatus) and their response to intestinal air-breathing inhibition. PLoS One. 2016;11(2):e0149123.

  10. Tata JR. Amphibian metamorphosis as a model for studying the developmental actions of thyroid hormone. Biochimie. 1999;81(4):359–66.

    Article  PubMed  CAS  Google Scholar 

  11. Hasebe T, Fujimoto K, Kajita M, Ishizuya-Oka A. Thyroid hormone activates Wnt/β-catenin signaling involved in adult epithelial development during intestinal remodeling in Xenopus laevis. Cell Tissue Res. 2016;365(2):309–18.

    Article  PubMed  CAS  Google Scholar 

  12. Ishizuya-Oka A, Kajita M, Hasebe T. Thyroid hormone-regulated wnt5a/ror2 signaling is essential for dedifferentiation of larval epithelial cells into adult stem cells in the Xenopus laevis intestine. PLoS One. 2014;9(9):e107611.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Sun G, Heimeier RA, Fu L, Hasebe T, Das B, Ishizuya-Oka A, Shi Y-B. Expression profiling of intestinal tissues implicates tissue-specific genes and pathways essential for thyroid hormone-induced adult stem cell development. Endocrinology. 2013;154(11):4396–407.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. Nakai Y, Nakajima K, Robert J, Yaoita Y. Ouro proteins are not essential to tail regression during Xenopus tropicalis metamorphosis. Genes Cells. 2016;21(3):275–86.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Grimaldi A, Buisine N, Miller T, Shi Y-B, Sachs LM. Mechanisms of thyroid hormone receptor action during development: lessons from amphibian studies. Biochim Biophys Acta Gen Subj. 2013;1830(7):3882–92.

    Article  CAS  Google Scholar 

  16. Ishizuya-Oka A, Hasebe T, Shimizu K, Suzuki K, Ueda S. Shh/BMP-4 signaling pathway is essential for intestinal epithelial development during Xenopus larval-to-adult remodeling. Dev Dyn. 2006;235(12):3240–9.

    Article  PubMed  CAS  Google Scholar 

  17. Coen L, Le Blay K, Rowe I, Demeneix BA. Caspase-9 regulates apoptosis/proliferation balance during metamorphic brain remodeling in Xenopus. Proc Natl Acad Sci. 2007;104(20):8502–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Krol J, Loedige I, Filipowicz W. The widespread regulation of microRNA biogenesis, function and decay. Nat Rev Genet. 2010;11(9):597–610.

    Article  PubMed  CAS  Google Scholar 

  19. Bartel DP. MicroRNAs: Genomics, biogenesis, mechanism, and function. Cell. 2004;116:281–97. Cell 2007, 131(4):11–29

    Article  PubMed  CAS  Google Scholar 

  20. Armisen J, Gilchrist MJ, Wilczynska A, Standart N, Miska EA. Abundant and dynamically expressed miRNAs, piRNAs, and other small RNAs in the vertebrate Xenopus tropicalis. Genome Res. 2009;19(10):1766–75.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Faunes F, Almonacid LI, Melo F, Larrain J. Characterization of small RNAs in Xenopus tropicalis gastrulae. Genesis. 2012;50(3):260–70.

    Article  PubMed  CAS  Google Scholar 

  22. Harding JL, Horswell S, Heliot C, Armisen J, Zimmerman LB, Luscombe NM, Miska EA, Hill CS. Small RNA profiling of Xenopus embryos reveals novel miRNAs and a new class of small RNAs derived from intronic transposable elements. Genome Res. 2014;24(1):96–106.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Ahmed A, Ward NJ, Moxon S, Lopez-Gomollon S, Viaut C, Tomlinson ML, Patrushev I, Gilchrist MJ, Dalmay T, Dotlic D, et al. A database of microRNA expression patterns in Xenopus laevis. PLoS One. 2015;10(10):e0138313.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Hikosaka A, Takaya K, Jinno M, Kawahara A. Identification and expression-profiling of Xenopus tropicalis miRNAs including plant miRNA-like RNAs at metamorphosis. FEBS Lett. 2007;581(16):3013–8.

    Article  PubMed  CAS  Google Scholar 

  25. Liu L-S, Zhao L-Y, Wang S-H, Jiang J-P. Research proceedings on amphibian model organisms. Dong wu xue yan jiu = Zoological research. 2016;37(4):237–45.

    PubMed  PubMed Central  Google Scholar 

  26. Wang S, Zhao L, Liu L, Yang D, Khatiwada JR, Wang B, Jiang J. A complete embryonic developmental table of Microhyla fissipes (Amphibia, Anura, Microhylidae). Asian Herpetological Research. 2017;8(2):108–17.

    Google Scholar 

  27. Yi S, Gao Z-X, Zhao H, Zeng C, Luo W, Chen B, Wang W-M. Identification and characterization of microRNAs involved in growth of blunt snout bream (Megalobrama amblycephala) by Solexa sequencing. BMC Genomics. 2013;14:754.

  28. Enright A, John B, Gaul U, Tuschl T, Sander C, Marks D. MicroRNA targets in Drosophila. Genome Biol. 2003;4(11):P8.

    Article  Google Scholar 

  29. Gotz S, Garcia-Gomez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talon M, Dopazo J, Conesa A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36(10):3420–35.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li C-Y, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39:W316–22.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Sterling J, Fu L, Matsuura K, Shi Y-B. Cytological and morphological analyses reveal distinct features of intestinal development during Xenopus tropicalis metamorphosis. PLoS One. 2012;7(10):e47407.

  32. Heimeier RA, Das B, Buchholz DR, Fiorentino M, Shi Y-B. Studies on Xenopus laevis intestine reveal biological pathways underlying vertebrate gut adaptation from embryo to adult. Genome Biol. 2010;11(5):R55.

  33. Loh Y-HE, Yi SV, Streelman JT. Evolution of micrornas and the diversification of species. Genome Biology and Evolution. 2011;3:55–65.

    Article  PubMed  CAS  Google Scholar 

  34. Juanchich A, Bardou P, Rue O, Gabillard J-C, Gaspin C, Bobe J, Guiguen Y. Characterization of an extensive rainbow trout miRNA transcriptome by next generation sequencing. BMC Genomics. 2016;17:164.

  35. Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, Weigel D. Specific effects of micrornas on the plant transcriptome. Dev Cell. 2005;8(4):517–27.

    Article  PubMed  CAS  Google Scholar 

  36. Wen L, Shi Y-B. Unliganded thyroid hormone receptor alpha controls developmental timing in Xenopus tropicalis. Endocrinology. 2015;156(2):721–34.

    Article  PubMed  CAS  Google Scholar 

  37. Choi J, K-iT S, Sakuma T, Shewade L, Yamamoto T, Buchholz DR. Unliganded thyroid hormone receptor alpha regulates developmental timing via gene repression in Xenopus tropicalis. Endocrinology. 2015;156(2):735–44.

    Article  PubMed  CAS  Google Scholar 

  38. Yao Y, Ma L, Jia Q, Deng W, Liu Z, Zhang Y, Ren J, Xue Y, Jia H, Yang Q. Systematic characterization of small RNAome during zebrafish early developmental stages. BMC Genomics. 2014;15:117.

  39. Wang B, Gan Z, Cai S, Wang Z, Yu D, Lin Z, Lu Y, Wu Z, Jian J. Comprehensive identification and profiling of Nile tilapia (Oreochromis niloticus) microRNAs response to Streptococcus agalactiae infection through high-throughput sequencing. Fish & Shellfish Immunology. 2016;54:93–106.

    Article  CAS  Google Scholar 

  40. Kim JY, Park Y-K, Lee K-P, Lee S-M, Kang T-W, Kim H-J, Dho SH, Kim S-Y, Kwon K-S. Genome-wide profiling of the microRNA-mRNA regulatory network in skeletal muscle with aging. Aging-Us. 2014;6(7):524–44.

    Article  CAS  Google Scholar 

  41. Chu WY, Liu LS, Li YL, Chen L, Wang KZ, Li HH, Du SJ, Zhang JS. Systematic identification and differential expression profiling of micrornas from white and red muscles of Siniperca chuatsi. Curr Mol Med. 2013;13(8):1397–407.

    Article  PubMed  CAS  Google Scholar 

  42. Manni I, Artuso S, Careccia S, Rizzo MG, Baserga R, Piaggio G, Sacchi A. The microRNA miR-92 increases proliferation of myeloid cells and by targeting p63 modulates the abundance of its isoforms. FASEB J. 2009;23(11):3957–66.

    Article  PubMed  CAS  Google Scholar 

  43. Nishikawa A: Cell interaction during larval-to-adult muscle remodeling in the frog, Xenopus laevis. In: Gowder S, editor. Cell interaction. Rijeka: InTech; 2012: Ch. 01,4-8.

  44. Liu Y, Zhang B, Shi T, Qin H. miR-182 promotes tumor growth and increases chemoresistance of human anaplastic thyroid cancer by targeting tripartite motif 8. OncoTargets and therapy. 2017;10:1115–22.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Fu Y, Shi Z, Wang G, Zhang J, Li W, Jia L. Expression of let-7 microRNAs that are involved in Japanese flounder (Paralichthys olivaceus) metamorphosis. Comparative Biochemistry and Physiology B-Biochemistry & Molecular Biology. 2013;165(2):106–13.

    Article  CAS  Google Scholar 

  46. Mi X, Wei Z, Zhou Z, Liu X. Identification and profiling of sex-biased microRNAs from sea urchin Strongylocentrotus nudus gonad by Solexa deep sequencing. Comparative Biochemistry and Physiology D-Genomics & Proteomics. 2014;10:1–8.

    Article  CAS  Google Scholar 

  47. Roush S, Slack FJ. The let-7 family of microRNAs. Trends Cell Biol. 2008;18(10):505–16.

    Article  PubMed  CAS  Google Scholar 

  48. Okamura K, Ishizuka A, Siomi H, Siomi MC. Distinct roles for argonaute proteins in small RNA-directed RNA cleavage pathways. Genes Dev. 2004;18(14):1655–66.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Guo L, Lu Z. The fate of miRNA* strand through evolutionary analysis: implication for degradation as merely carrier strand or potential regulatory molecule? PLoS One. 2010;5(6):e11387.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Chi W, Tong C, Gan X, He S. Characterization and comparative profiling of miRNA transcriptomes in bighead carp and silver carp. PLoS One. 2011;6(8):e23549.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Xu Z, Chen J, Li X, Ge J, Pan J, Xu X. Identification and characterization of micrornas in channel catfish (Ictalurus punctatus) by using solexa sequencing technology. PLoS One. 2013;8(1):e54174.

  52. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120(1):15–20.

    Article  PubMed  CAS  Google Scholar 

  53. Aherne ST, Smyth P, Freeley M, Smith L, Spillane C, O'Leary J, Sheils O. Altered expression of mir-222 and mir-25 influences diverse gene expression changes in transformed normal and anaplastic thyroid cells, and impacts on MEK and TRAIL protein expression. Int J Mol Med. 2016;38(2):433–45.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. Townley-Tilson WHD, Callis TE, Wang D. MicroRNAs 1, 133, and 206: critical factors of skeletal and cardiac muscle development, function, and disease. International Journal of Biochemistry & Cell Biology. 2010;42(8):1252–5.

    Article  CAS  Google Scholar 

  55. Sapre N, Selth LA. Circulating MicroRNAs as biomarkers of prostate Cancer: the state of play. Prostate cancer. 2013;2013:539680.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. Zhang L, Zhou F, ten Dijke P. Signaling interplay between transforming growth factor-β receptor and PI3K/AKT pathways in cancer. Trends Biochem Sci. 2013;38(12):612–20.

    Article  PubMed  CAS  Google Scholar 

  57. Damjanovski S, Amano T, Li Q, Ueda S, Shi YB, Ishizuya-Oka A. Role of ECM remodeling in thyroid hormone-dependent apoptosis during anuran metamorphosis. Ann N Y Acad Sci. 2000;926:180–91.

    Article  PubMed  CAS  Google Scholar 

  58. Ciechanover A, Orian A, Schwartz AL. Ubiquitin-mediated proteolysis: biological regulation via destruction. Bioessays. 2000;22(5):442–51.

    Article  PubMed  CAS  Google Scholar 

  59. Nayak L, Bhattacharyya NP, De RK. Wnt signal transduction pathways: modules, development and evolution. BMC Syst Biol. 2016;10(Suppl 2):44.

  60. Lee DY, Lin TE, Lee CI, Zhou J, Huang YH, Lee PL, Shih YT, Chien S, Chiu JJ. MicroRNA-10a is crucial for endothelial response to different flow patterns via interaction of retinoid acid receptors and histone deacetylases. Proceedings of the National Academy of Sciences. 2017;114(8):2072–7.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank Yan Wang for M. fissipes rearing.

Funding

This study was supported by Important Research Project of Chinese Academy of Sciences (KJZG-EW-L13) and Western Light Talent Culture Project of the Chinese Academy of Sciences (2015XBZG_XBQNXZ_B_011). The funding bodies have no role in study design, data collection and analysis, or prepration of the manuscript.

Availability of data and materials

The raw reads were deposited in the NCBI Sequence Read Archive (SRA) database under the study accession SRR5451121, SRR5451123, SRR5451131, SRR5451133, SRR5451135, SRR5451136.

Author information

Authors and Affiliations

  1. Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu, 610041, China

    Lusha Liu, Wei Zhu, Jiongyu Liu, Shouhong Wang & Jianping Jiang

  2. University of Chinese Academy of Sciences, Beijing, 100049, China

    Shouhong Wang

Authors
  1. Lusha Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Wei Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jiongyu Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Shouhong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Jianping Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LSL and JPJ conceived and designed the study. LSL performed the molecular experiments, analyzed the data and wrote the manuscript. WZ and JYL assisted with the bioinformatics analysis. SHW collected of samples. All authors helped to revise the manuscripts and approved the final manuscripts.

Corresponding authors

Correspondence to Lusha Liu or Jianping Jiang.

Ethics declarations

Ethics approval and consent to participate

Mature female and male M. fissipes were collected from Shuangliu, Chengdu, China (30.5825o N, 103.8438o E). As prescribed by “Law of People’s Republic of China on the Protection of Wildlife” and “Regulations for the Implementation of the People’s Republic of China on the Protection of terrestrial Wildlife” (State Council Decree [1992] No. 13), the specimen collection did not ask for any ethical or institutional approvals, because M. fissipes is not endangered or protected by any current law. The care and treatment of animals in this study were performed according to the Guideline for the Care and Use of Laboratory Animals in China. The animal experiments (including LHRHa and T3 treatment, morphological, histological and molecular biology experiments) were approved by the Experimental Animal Use Ethics Committee of the Chengdu Institute of Biology (Permit Number: 2016036).

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Primers sequences for qRT-PCR. (DOC 38 kb)

Additional file 2:

Figure S1. Morphological characteristics of M. fissipes tadpoles during natural metamorphosis. TOL: total length; TL: tail length; BW: body width; SL: snout length; IOS: interocular space; Pre-: premetamorphosis; Pro-: prometamorphosis; climax: the climax of metamorphosis; End: end of metamorphosis. (TIF 761 kb)

Additional file 3:

Table S2. Overview of readcounts for sRNA-seq from the raw data to high quality reads, and quality filtering. (DOC 42 kb)

Additional file 4:

Figure S2. a. Length distribution and abundance of small RNA sequences in M. fissipes, as determined by Illumina small-RNA deep sequencing. b. Nucleotides bias on the specific position of miRNAs in M. fissipes c. Count distribution and abundance of unique small RNA sequences in M. fissipes. d. Scatter plot map for miRNAs expression in the control and T3 groups. Each plot represented an individual miRNA, while the red plot indicated the significantly differentially expressed miRNA (p < 0.01 and |log2 (foldchange)| > 1). (JPG 1976 kb)

Additional file 5:

Figure S3. Pie charts of different abundance of small RNA in control group and T3 treated group. (TIF 579 kb)

Additional file 6:

Table S3. Conserved miRNAs in M. fissipes. (XLS 61 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, L., Zhu, W., Liu, J. et al. Identification and differential regulation of microRNAs during thyroid hormone-dependent metamorphosis in Microhyla fissipes. BMC Genomics 19, 507 (2018). https://doi.org/10.1186/s12864-018-4848-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4848-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us