Comparative expression analysis identifies the respiratory transition-related miRNAs and their target genes in tissues of metamorphosing Chinese giant salamander (Andrias davidianus)

Comparative expression analysis identifies the respiratory transition-related miRNAs and their... Background: Chinese giant salamander (Andrias davidianus) undergoes a metamorphosis from aquatic larvae to terrestrial adults, with concomitant transfer of respiration from gills to lungs prior to metamorphosis. These two tissues, as well as skin, were sampled to identify the differentially expressed miRNAs. Results: High-coverage reference transcriptome was generated from combined gill, lung and skin tissues of metamorphosing juveniles, and lung tissue of adults: 86,282 unigenes with total length of approximately 77,275,634 bp and N50 of 1732 bp were obtained. Among these, 13,246 unigenes were assigned to 288 pathways. To determine the possible involvement of miRNAs in the respiratory transition, small RNA libraries were sequenced; 282 miRNAs were identified, 65 among which were known and 217 novel. Based on the hierarchical clustering analysis, the twelve studied samples were classified into three major clusters using differentially expressed miRNAs. We have validated ten differentially expressed miRNAs and some of their related target genes using qPCR. These results largely corroborated the results of transcriptomic and miRNA analyses. Finally, an miRNA-gene-network was constructed. Among them, two miRNAs with target genes related to oxygen sensing were differentially expressed between gill and lung tissues. Three miRNAs were differentially expressed between the lungs of larvae and lungs of adults. Conclusions: This study provides the first large-scale miRNA expression profile overview during the respiration transition from gills to lungs in Chinese giant salamander. Five differentially expressed miRNAs and their target genes were identified among skin, gill and lung tissues. These results suggest that miRNA profiles in respiratory tissues play an important role in the regulation of respiratory transition. Keywords: Deep sequencing, Respiratory system, miRNA, Andrias davidianus * Correspondence: 362205379@qq.com; yuanxh@ffrc.cn Shengyan Su, Yuheng Wang and Huiwei Wang contributed equally to this work. Department of Animal Husbandry & Veterinary Medicine, Jiangsu Polytechnic College of Agriculture and Forestry, Zhenjiang 212400, People’s Republic of China Key Laboratory of Genetic Breeding and Aquaculture Biology of Freshwater Fishes, Ministry of Agriculture; Freshwater Fisheries Research Center, Chinese Academy of Fishery Sciences, Wuxi 214081, People’s Republic of China Full list of author information is available at the end of the article © The Author(s). 2018 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. Su et al. BMC Genomics (2018) 19:406 Page 2 of 15 Background pattern of the expression of these target mRNAs, micro- The endangered Chinese giant salamander (Andrias RNAs (miRNAs) play an important role by specifically davidianus, Cryptobranchidae family), endemic to binding to the 3’UTR of mRNA and (usually) promoting China, is the world’s largest amphibian. It has been mRNA degradation [15]. miRNAs are endogenous ~ 22 called a “living fossil”, and it is referred to locally as nucleotide-long non-coding RNAs, which regulate a wawayu (baby fish) because its vocalization is somewhat range of essential cellular and biological processes by resemblant of a baby’s cry [1]. The Chinese government targeting mRNA for cleavage or translational repression has declared this species a class II protected species [2]. [16–23]. There are indications that miRNAs modulate Its ability to switch from aquatic to aerial respiration pulmonary development and maintain lung homeostasis and the associated morpho-functional adjustments [3, 4] [24, 25]. Misexpression of some miRNAs can contribute further contribute to its status of a species of great inter- towards neurodevelopmental disorders, which can also est to biologists. Studies of this species might provide affect the lung ventilation [26]. key information about body adaptions (anatomical, With two key developmental stages, Chinese giant physiological and molecular) that enable organisms to salamander is a typical amphibian, undergoing a transi- transition from water to land. tion from aquatic larvae to terrestrial adults, with con- Piiper (1982) divided the transition from water to land comitant transfer of respiration from gills to lungs prior into three sub-stages dependent on different organs: fish to the metamorphosis [27]. In order to unveil the role of gills, amphibian gills, and lungs and skin. The features of miRNAs in the adjustment of respiratory system to the respiration in transitional states showed that the adjust- transition from water to land, we have studied the ex- ment for gas exchange style in different breathing pression of miRNAs and their corresponding target models is related to the metabolic rate [5]. Burggren and genes in metamorphosing larvae of the Chinese giant Doyle (1986) reported that the transition from gill res- salamander. Real time quantitative polymerase chain re- piration to lung respiration during the development of action (qPCR) was used to study expression in gill, lung bullfrog larvae is accompanied by a decline in the regu- and skin tissues. We also assembled a high-coverage ref- lation of branchial ventilation frequency reflex and at- erence transcriptome, and determined the putative regu- tendant increasing dependence on breathing air [6]. latory network of miRNAs associated with the studied Based on the neurorespiratory pattern of gill and lung morpho-functional adjustments. This study contributes ventilation, physiologists found that spinal nerve II is to the understanding of the role of miRNA-mediated crucial for lung ventilation [7]. However, these studies regulatory networks in the transition from water to land merely identified the differences among different respir- of Chinese giant salamander. ation styles. In transitional forms, depending on both gill and lung ventilation, gill and lung breathing can be re- Methods placed or complemented by skin and or buccopharyn- Sample collection and RNA isolation geal breathing [8]. This adaptation is necessary for the Healthy Chinese giant salamander’s larvae (5 months-old) adaptive changes of respiratory properties of blood and and adults (20 months-old) were obtained from the ventilatory and circulatory flow rates [8]. Nitric oxide Jurong Chinese Giant Salamander Breeding Center in synthase (NOS)/nitric oxide (NO) levels in gills and Zhenjiang, Jiangsu, China. Specimens were anaesthetized lungs of Protopterus annectens during aestivation and using 2-phenoxyethanol, euthanized using tricaine metha- arousal changed in a tissue-specific manner in parallel nesulfonate (Sigma, USA), and promptly dissected. Three with organ readjustment [9]. However, despite numerous different tissues were collected from three larvae: skin efforts, signal transduction mechanisms that underlie (skin5m), gills (gill5m), and lungs (lungs5m); and only the cyclic branchial and pulmonary remodeling remain lungs from three adults (lungs20m) (Fig.1;Additional file 1: incompletely explained. Assunta et al. (2016) studied the Table S1). This ensured triplicate samples for four differ- water-land transition from the perspective of gene ex- ent tissue groups (skin5m, gill5m, lungs5m, lungs20m). pression patterns; they used lungfish transcriptome se- RNA was extracted from all samples using Trizol Reagent quence (which is a time-effective and economical (Invitrogen, USA) according to the manufacturer’sinstruc- method [10–13]) and phylogenetic relationship analysis tions. Purity of RNA was checked by NanoPhotometer® [14]. They identified 226 concatenated protein sequences spectrophotometer (IMPLEN, CA, USA). The final con- in all vertebrates and 59, 951 informative amino acid po- centration was determined using Qubit® RNA Assay Kit in sitions in both lungfish/tetrapod and the coelacanth/ Qubit®2.0 Fluorometer (Life Technologies, CA, and USA). lungfish + tetrapod nodes. They also identified several The integrity of RNA was assessed using the RNA Nano proteins in the tetrapod lung which were not found in 6000 Assay Kit and Agilent Bioanalyzer 2100 system the actinopterygian fish (pulmonary surfactants A, C, D (Agilent Technologies, CA, and USA). The handling of SFTPA, SFTPC, and SFTPD). Regarding the regulation animals was conducted in accordance with the guidelines Su et al. BMC Genomics (2018) 19:406 Page 3 of 15 Fig. 1 Chinese giant salamander respiratory system transition from gills, lungs and skin to respiration mainly using lungs (a) A transitioning Chinese giant salamander with gills and lungs; (b) Chinese giant salamander with lungs and without gills for the care and use of animals for scientific purposes set cBot Cluster Generation System using TruSeq PE Cluster by the Institutional Animal Care and Use Committee of the Kit v3-cBot-HS (Illumina) according to the manufacturer’s Freshwater Fisheries Research Center, Chinese Academy of instructions. Finally, the libraries were sequenced using Fishery Sciences, Wuxi, China. paired-end reads on an Illumina Hiseq 4000 platform. Library preparation and transcriptome sequencing De novo transcriptome assembly and prediction of Sequencing libraries were constructed using 3 μgof coding sequence RNA per sample using NEBNext®Ultra™ RNA Library In order to obtain high-quality clean reads, the raw data Prep Kit for Illumina® (NEB, USA), following the manu- were trimmed to remove reads containing adapter, reads facturer’s instructions. Index codes were added to attri- containing ploy-N and low-quality reads. Q20, Q30, GC- bute sequences to each sample. Briefly: mRNA was content and sequence duplication level of the clean data purified from total RNA using poly-T oligo-attached were calculated. Clean data were then assembled de magnetic beads and fragmented using divalent cations novo as described before [28]: the sequenced ‘left’ and under elevated temperature in NEBNext First Strand ‘right’ files (read1 and read2 files respectively) from all Synthesis Reaction Buffer (5X). First strand cDNA was samples were pooled into left.fq and right.fq files, re- transcribed using random hexamer primer and M-MuLV spectively. Transcriptome assembly was finalized using Reverse Transcriptase (RNase H). Second strand cDNA these two large files, and the longest fragments per gene synthesis was then performed using DNA polymerase I were designated as unigenes. All these operations were and RNase H. Remaining overhangs were converted into conducted using Trinity [29], with all parameters set to blunt ends using exonuclease/polymerase. After adenyla- default. tion of 3′ ends of DNA fragments, NEBNext adaptors with hairpin loop structure were ligated to prepare for Gene function annotation, SSR detection and primer hybridization. Library fragments were purified with design AMPure XP system (Beckman Coulter, Beverly, USA) to These unigenes were then queried against public databases select cDNA fragments with the length of 150~ 250 bp. for functional annotation: NCBI non-redundant protein In order to generate the products, we used 3 μl USER database using BLASTX searches [30], Protein Families Enzyme (NEW ENGLAND BioLabs Inc., USA) with (Pfam) [31], orthologous groups of genes (EggNOG) [32] size-selected, adaptor-ligated cDNA at 37 °C for 15 min, and Swiss-Prot database [33]. The cut-off of value for anno- followed by 5 min at 95 °C. PCR was performed with tation in NR, Swissport, KEGG databases was set E to 1-e5 Phusion High-Fidelity DNA polymerase, universal PCR using blastx, whereas for the Pfam database the value was primers and index (X) primer. Ultimately, we purified set E to 10 using hmmscan. In addition, the predicted CDS PCR products with AMPure XP system and assessed the were classified into EggNOG categories using HMMER library quality on the Agilent Bioanalyzer 2100 system. version 3.1 with an E-value cutoff of 1E-05 [34]. Kyoto We performed clustering of the index-coded samples on a Encyclopedia of Genes and Genomes database (KEGG) Su et al. BMC Genomics (2018) 19:406 Page 4 of 15 [35] was used to associate unigenes with metabolic path- the Agilent Bioanalyzer 2100 system. We performed the ways. Gene Ontology database (GO) [36] was used to clas- clustering of index-coded samples on a cBot Cluster sify the unigenes functionally into biological process, Generation System using TruSeq PE Cluster Kit v4- molecular functions and cellular components categories. cBot-HS (Illumina) according to the manufacturer′s in- GO term analysis was conducted by blast2go using the Nr structions. After cluster generation, library preparations annotation information. TransDecoder was used to predict were sequenced on an Illumina HiSeq X Ten platform the sequence of the Unigene coding region and its corre- and single-end reads were generated. sponding amino acid sequence, which was recommended by Trinity and Cuffinks. MIcroSAtellite identification tool Small RNA annotation, new miRNAs, and prediction of (MISA) [37] was used to detect SSRs in the transcriptome their target genes with the following parameters: definition (unit size - min. For small RNA sequencing reads, all sequences with repeats) = 1–10 2–63–54–55–56–5, and interruption adapter or poly-N contaminants and low-quality reads (max. Difference between two SSRs) = 100. Six types of SSRs were flittered. After that, sequences < 15 and > 35 nt were investigated:mono-,di-,tri-, tetra-, penta-, and hexa- were removed. Q20, Q30, GC-content and sequence du- nucleotide repeats. Primers for SSRs were designed using plication level of the clean data were calculated as well. Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/). Using Bowtie software, the clean reads were aligned with Rfam, Silva, GtRNAdb and Repbase databases to filter Quantification of gene expression levels the transfer RNAs (tRNA), small nuclear RNAs (snRNA) The read count for each gene was obtained by mapping , small nucleolar RNAs (snoRNA), ribosomal RNAs clean reads to assembled transcriptome using RSEM (rRNA) and other ncRNAs, as well as repeats. By com- [38]. Expression levels of unigenes were calculated as paring the remaining reads with known miRNAs from fragments per kb per million reads (FPKM). Differen- the miRBase, known miRNAs were detected and new tially expressed genes (DEGs) were identified using miRNAs predicted. Non-conserved miRNAs were pre- DESeq R package (1.10.1) [39], employing a model based dicted using miRdeep2 [41] with the following run pa- on the negative binomial distribution. Benjamini and rameters: “-g -1 -b 0”. Secondary structure prediction of Hochberg’s approach was used to adjust the resulting P new miRNAs was performed using Randfold tools [42] value, which was determined by false discovery rate with the parameters: s = simple mononucleotide shuf- (FDR) and fold change (FC) [40]. In the pairwise analysis fling, and number of randomizations = 99. As we did not of the four groups (12 samples with 3 replicates) de- have the whole genome sequence at disposal, RNA-seq scribed before, any gene with an adjusted P < 0.05 and data were used to predict new microRNAs. FC ≥ 2 was considered to be a differentially expressed gene (DEG). miRNA expression levels and functional annotation of their target genes Library preparation and sRNA sequencing To estimate the miRNA expression levels for each sam- A total amount of 1.5 μg RNA per sample was used as ple, sRNAs were mapped back onto the precursor se- input material for the RNA sample preparations. Se- quence [43]. From the mapping results we know the quencing libraries were generated using NEBNext®Ultra™ readcount, which represents the number of reads that small RNA Sample Library Prep Kit for Illumina® (NEB, are mapped to a certain miRNA. Then the expression of USA) following manufacturer’s recommendations, and miRNA in each sample was normalized as transcripts index codes were added to attribute sequences to each per million (TPM): TPM = readcount× 1,000,000/ sample. First of all, 3’ SR Adaptor for Illumina was Mapped Reads. Potential target genes of differentially mixed with RNA and nuclease-free water. This was expressed miRNAs were identified using miRnada v3.3a followed by incubation for 2 min at 70 °C in a preheated (run parameters: -sc 150.0 -en − 30 -scale 4.0 -go − 2.0 thermal cycler. Tube was then transferred to ice, and 3’ -ge − 8.0) and RNA hybrid v2.1.1 (run parameters: -m Ligation Reaction Buffer (2X) and 3’ Ligation Enzyme 50,000 -d 1.9,0.28 -b 1 -e − 30) programs. Intersection of Mix. The mixture was incubated for 1 h at 25 °C in a the target set was regarded as the target gene of a micro- thermal cycler. To prevent adaptor-dimer formation, the RNA. After these step, intersection of predicted targets SR RT Primer hybridizes to the excess of 3′ SR Adaptor were verified by target accessibility (PITA) [44].Target that has remained free after the 3′ ligation reaction, and genes of miRNAs were subjected to BLASTX searches transforms the single-stranded DNA adaptor into a (E-value threshold of 1.00e-5) and annotated against the double-stranded DNA molecule. After ligating the 5′ SR NCBI’s Nr and Pfam databases, as well as the EggNOG Adaptor, reverse transcription and PCR amplification and Swiss-Prot databases again. Target genes were also were performed. At last, PCR products were purified subjected to GO analysis, and classified into biological (AMPure XP system) and library quality was assessed on process, cellular component and molecular function Su et al. BMC Genomics (2018) 19:406 Page 5 of 15 categories. Metabolic pathways were analyzed using the assemble all the trimmed reads. About 77.83% of clean KEGG database. reads were mapped to the newly assembled transcrip- tome. Finally, a total of 86,282 unigenes with total length Validation of differentially expressed miRNAs and their of approximately 77,275,634 bp and N50 of 1732 bp corresponding target genes were obtained (Table 1). These values are higher than Expression levels of the following randomly selected the number of unigenes and N50 in previously reported miRNAs were determined using predesigned TaqMan Chinese giant salamander transcriptomes [47, 48]. MicroRNA assay (Life Technologies, Foster City, CA, Among these, 30,060 unigenes (34.84%) were < 300 bp in USA): aca-miR-142-5p, bta-miR-142-5p, gga-miR-34c-3p, length, 33,789 (39.16%) were > 301 and < 1000 bp, and aca-let-7a-5p, aca-miR-203-3p, aca-miR-203-5p, aja-miR- 22,433 (26.00%) were > 1000 bp (Additional file 2). The 142, dre-miR-203a-5p, efu-miR-223 and ipu-miR-142. overall GC ratio of the unigenes was 49.19%. Most of cDNA was synthesized from 100 ng of total RNA using the CDS were 200~ 300 bp in length (19,941, 31.01%), TaqMan® MicroRNA Reverse Transcription kit (Life Tech- followed by 100~ 200 bp (15,984, 24%), and 300~ 400 bp nologies). qPCR was performed using 7900HT Fast Real- (7652, 11.9%) (as predicted by Trans-Decoder). No CDS Time PCR System (Life Technologies), with primers listed with length < 100 were found, but 1130 (1.76%) of CDS in Additional file 1: Tables S2 and S3. Relative expression were > 3000 bp (Additional file 1: Table S5). A total of of miRNAs was determined using the 2-ΔΔCT method 13,199 SSRs were identified from 22,433 unigenes (> 1 kb) [45], with mamm U6 small nuclear 1 (U6), which belongs using MISA Perl script, with 816 (3.64%) unigenes con- to the class of small nuclear RNAs, used as an en- taining more than one SSR. Among the different types of dogenous reference [46]. No significant differences in SSRs, mono-nucleotide repeats were the most abundant the U6 expression were found among the four groups (10,332), followed by di-nucleotide (1233), tri-nucleotide (Additional file 1: Tables S4). All reactions were per- (687), tetra-nucleotide (88) and penta-nucleotide repeats formed in triplicate [37], data expressed as mean ± SE, (4) (Additional file 1: Table S6). Frequency distribution of and subjected to one-way analysis of variance these putative cSSRs was also calculated (Additional file 3). (ANOVA) using SAS8.0 program. On average, one cSSR locus was found for about every 5.9 kb of a unigene sequence. Results BLASTX and BLASTn searches against GO, KEGG, Assembly and annotation of Andrias davidianus eggnog and Nr databases were conducted to validate transcriptome and annotate the assembled unigenes; this produced To obtain a reference transcriptome for the study of 6390 (21.27%), 9236 (30.74%), 13,246 (44.09%), 16,180 microRNA regulation patterns during the transition (53.85%), 27,181 (90.48%), 28,247 (94.02%) unigene from gill respiration to lung respiration, samples from matches respectively (E-value ≤1.0E-5) (Additional file 1: individuals which only had lung tissue and from those Table S7). A set of 30,044 unigenes were annotated in possessing both gill and lung tissues at the same time these public databases via Trans-Decoder annotation. were used for RNA-Seq analysis. Their mix was used to Functional annotation of all unigenes against protein da- obtain the transcriptome sequence. A total of tabases showed that a total of 18,435 (61.36%) unigenes 41,878,756 bp of clean data were obtained after quality exhibit significant similarity to known proteins in the filtering (Table 1). As there are no available assembled Pfam database, and a total of 13,341 (44.40%) in the and annotated genomic sequences of the Chinese giant Swissprot database. Four species with most BLASTX salamander that could be used for transcriptome assem- hits against the Nr protein database accounted for bly, Trinity de novo assembler program was used to 42.86% of the identified unigenes: Xenopus tropicalis (14.16%), Chrysemys picta (13.89%), Chelonia mydas (7.76%) and Latimeria chalumnae (7.05%). These re- Table 1 Assembly quality statistics of the Chinese giant salamander sults are comparable to those reported in the larvae transcriptome of tiger salamander and the amphibian species usu- Parameter Value ally used for the respiratory transition studies - Xen- Raw reads 41,878,756 opus tropicalis [49]. The remaining 57.14% were High quality bases 12,482,332,994 distributed as follows: Pelodiscus annectens (5.38%), High quality reads % (≥Q30) 89.39% Xenopus laevis (3.37%), Anolis carolinensis (3.22%), Oncorhynchus mykiss (3.16%), Alligator mississippiensis Unigene number 86,282 (2.90%), Alligator sinensis (2.82%), and other species Unigene total length 77,275,634 (36.28%) (Additional file 4). In total, 9236 (10.70%) uni- Average length (bp) 895.62 genes were allocated into three main GO categories Unigene N50 length (bp) 1732 (Additional file 5): within the cellular component category, Su et al. BMC Genomics (2018) 19:406 Page 6 of 15 cell (5172) and cell part (5171) terms were most might have effect on the miRNA target recognition [50]. prominent; within the molecular function category, To identify the conserved miRNAs in A. davidianus,small GO terms were predominantly associated with bind- RNA sequences were compared with the currently avail- ing activity (5022) and catalytic activity (3622); the able mature miRNAs in the miRBase. The miRNAs that biological process category was dominated by cellular perfectly matched with conserved miRNA sequences or process (5471) and single-organism process (5039). had stem-loop precursors were identified as conserved These results are almost identical to Huang et al.’sre- miRNAs, which indicates their conserved function across port [48] and similar to be found by Eo et al. [49]. In various species. Finally, a total of 65 small RNAs in A. the the latter study [49] these terms were associated davidianus were identified as orthologs of known miR- with transcripts differentially expressed between gills NAs in other species (Additional file 1: Table S10). For and lungs. conserved miRNA, most of them ranged from 21 to 23 nt The unigenes were also queried (by BLASTX) in in length; sorted by abundance: 22 nt > 21 nt > 23 nt. KEGG database to further explore their biological func- These results are fully consistent with Yang et al.’s results tions and interactions. A total of 13,246 unigenes were in Chinese giant salamander [51]. Abundance of known assigned to 288 pathways, among which ribosome (345), miRNA family members ranged from 1 to 13. The most MAPK signaling pathway (311) and focal adhesion (303) abundant were let-1 and mir-142 (13 and 12 members re- were the largest three groups. In addition, 189 unigenes spectively). Among the remaining, two families had 9 were associated with oxidative phosphorylation, 61 with members, and four families had > 3 members and < 9 glutathione metabolism, 46 with metabolism of xenobi- members. To uncover novel miRNAs sequences, all non- otics by cytochrome P450 and 45 with drug metabolism- annotated small RNAs were searched against the newly- cytochrome P450. assembled transcriptome. We identified 46 non-conserved miRNAs belonging to 32 miRNA families, five of which Small RNA sequencing data and annotation (families) contained more than one member. Length dis- To assess the involvement of miRNAs in the respiratory tribution of non-conserved miRNAs was identical to that transition, small RNA libraries were constructed for se- of all miRNAs. Thus, length distribution of miRNAs quencing. A total of 183,422,095 clean reads were gener- might provide a window into the role of miRNAs compo- ated from 12 samples after removing low quality reads, nents in morpho-functional adjustments, because length adapter sequences, and sequences with lengths < 15 nt variation is related to enzymatic modifications, such as and > 35 nt (Additional file 1: Table S8). Q30 percentage RNA editing and exonuclease actives [52]. per individual ranged from 98.83% to 99.13%. Reads from all libraries were annotated using GenBank and Differential expression of miRNAs in gills, lungs and skin Rfam databases (Additional file 1: Table S9). By annota- DESeq program was used to determine whether miRNA tion, Rrnas, tRNAs, snRNAs and snoRNAs were sepa- profiles were different between salamanders depending rated from each other. tRNAs and rRNAs were the most predominantly on gills and skin for respiration and those abundant small RNA types. Length of mature miRNAs depending on lungs for respiration (Table 2). Between was mostly between 20 and 22 nt (27.65%), which were skin and gills in larvae, we found 77 differentially typical dicer-processed miRNA products (Additional file 6). expressed miRNAs, 46 of which were up-regulated and The most abundant small RNAs in all of the libraries were 31 were down-regulated. Among them, four miRNAs 21 nt in length, and a significant bias towards U was ob- were unique to skin and two to gills in larvae. All of served in the first nucleotides of all putative miRNA se- these miRNAs belong to novel miRNAs. Between skin quences across all libraries (Additional file 7). U also had and lungs in larvae, we identified 69 differentially high frequency (51%) in positions 9 and 12. Moreover, A expressed miRNAs, 37 of which were up-regulated and + U were dominant at the start position. These features 32 down–regulated. Among them, only one miRNA, Table 2 miRNA count among different tissues Type DEG Number Up-regulated Down-regulated skin_5 m vs gills_5m 77 46 31 skin_5 m vs lungs_5m 69 37 32 skin_5 m vs lungs_20m 82 55 27 gills_5 m vs lungs_5m 32 16 16 gills_5 m vs lungs_20m 45 30 15 lungs_5 m vs lungs_20m 3 1 2 5 m and 20 m indicate the age of sampled specimens Su et al. BMC Genomics (2018) 19:406 Page 7 of 15 named unconservative_c58354.graph_c0_367043, was identified by intersection of miRnada and RNA hybrid unique to skin in larvae, and seven miRNAs were unique prediction programs. The target intersection set was to gills in larvae. Between gills in larvae and lungs in verified by PITA program [44]. Among the 15,739 genes adults, we identified 82 differentially expressed miRNAs, predicted as possible targets, 8358 (53.10%) were suc- 55 up-regulated and 27 down–regulated. Among these, cessfully annotated against the NR, Swiss-Prot, GO, only one miRNA, named unconservative_c63346.graph_ KEGG and Pfam databases (Additional file 1: Table S11). c0_552043, was unique to skin in larvae and three miR- A total of 352, 291, 1016, 117 and 862 differentially NAs were unique to gills in larvae. Between gills and expressed miRNA target genes were found in pairwise lungs in larvae, we identified 32 differentially expressed comparisons of skin5 m vs gill5m, skin5 m vs lungs5m, miRNAs, half of which were up-regulated and half skin5 m vs lungs20m, gill5 m vs lungs5 m and gill5 m vs down–regulated. Here, six miRNAs were unique to lungs20m, respectively, but no differentially expressed lungs in larvae. Between gills in larvae and lungs in genes were found in the lungs5 m vs lungs20m compari- adults, we identified 45 differentially expressed miRNAs, son. Expression levels of ten target genes related to the 30 of which were up-regulated and 15 down–regulated. validated 12 miRNAs were also determined (Fig. 4). Among them, only one miRNA, named unconservative_ Among them, four genes had higher expression in skin, c60271.graph_c0_422386, was unique to gills in larvae while two genes had higher expression in gills (both and four miRNAs were unique to lungs in adults. Be- were compared to the lung tissue). One gene, semialde- tween the adult and larva lungs, three miRNAs were dif- hyde dehydrogenase, had higher expression in both gills ferentially expressed including one up-regulated and two and skin, while only centrosomal protein of 57 kDa iso- down-regulated. Results of all pairwise comparisons in- form X1 exhibited higher expression level in lungs. dicate that aca-miR-203-5p, dre-miR-203a-5p, ipu-miR- The top-scoring BLASTX hits against the Nr protein 203a and rno-miR-203a-5p were expressed only in the database indicated that the top four species were Chrys- lungs of larvae. Distribution of significant changes de- emys picta (17.12%), Xenopus (Silurana) (13.82%), Chelo- tected is shown as a volcano plot in the Additional file 8. nia mydas (8.22%) and Latimeria chalumnae (7.69%), According to the hierarchical clustering analysis, the 12 which accounted for nearly half of the miRNA targets studied samples were classified into three major clusters, identified against the Nr database. The remaining which corresponds to the three studied tissues (Fig. 2). miRNA targets were distributed among Pelodiscus Therefore, miRNA expression patterns indicate that sinensis (4.81%), Xenopus laevis (4.55%), Anolis caroli- miRNAs are involved in the regulation of adaptation of nensis (3.61%), Alligator mississippiensis (3.51%), A. respiratory tissues during the transition from aquatic lar- sinensis (3.20%), Python bivittatus (1.68%), and other vae to terrestrial adults. species (31.83%). Among the large number of differentially expressed miRNAs, we selected ten for further validation. These Construction of the miRNA-gene-network were chosen based on two criteria, their strong differen- The ultimate goal of our study was to test the hypothesis tial expression among the above four tissues (> 2-fold in- that miRNAs and their target genes are associated with crease) and their strong correlation with ventilation the transition from aquatic to aerial respiration. Thus, features, either strong positive or strong negative: a) no we identified five candidate miRNAs expressed in the expression in some tissues, while expression was ob- specific expression patterns in gill, skin and lung. In served in other tissue(s) in differential expression dataset order to identify putative functional modules, we (using DESeq); b) FDRs were lower than 0.05 for all constructed a miRNA-gene network (Fig. 5). Among pairwise comparisons. Based on these criteria, aca-miR- thefivemiRNAsinthe network, threeweredown- 142-5p, bta-miR-142-5p, gga-miR-34c-3p, aca-let-7a-5p, regulated during the respiratory transition: unconserva- aca-miR-203-3p, aca-miR-203-5p, aja-miR-142, dre-miR- tive_c62709.graph_c3_520410, unconservative_c41959. 203a-5p, efu-miR-223 and ipu-miR-142 were selected for graph_c0_140263 and unconservative_c60463.graph_ further validation experiments (Fig. 3). These results c1_429066; whereas one exhibited increased expression largely corroborated the results of transcriptomic during the respiratory transition: unconservative_ analysis. c60430.graph_c0_428023. The miRNAs exhibiting the widest and narrowest degrees in the miRNA-Gene Identification of target genes of differentially expressed network were unconservative_c60430.graph_c0_428023 miRNAs and their functional annotation and unconservative_c57205.graph_c0_339408, which The ultimate goal of our study was to test the hypothesis had 1361 and 3 target genes respectively. Regarding the that miRNAs and their target genes are associated with total differentially expressed genes, most of them were the transition from aquatic to aerial respiration. Poten- targeted by two and more miRNAs in the network, tial targets of differentially expressed miRNAs were whereas 18 genes were regulated by two miRNAs. Su et al. BMC Genomics (2018) 19:406 Page 8 of 15 Fig. 2 (See legend on next page.) Su et al. BMC Genomics (2018) 19:406 Page 9 of 15 (See figure on previous page.) Fig. 2 Hierarchical cluster analysis of differentially expressed miRNAs in the respiration tissue transfer process. Each row represents a different miRNA, and each column represents a different sample. Clustering based on parameter similarity is shown on the left (miRNAs) and above (samples) the figure. Increasing green intensity denotes decreasing gene expression, and increasing red intensity denotes increasing gene expression, with color scale shown on the right. The heatmap was constructed using expression profiling of each sample, with the data calculated from log10 (TPM + 1) GO categorization of these network targets indicated a 140263 and unconservative_c60463.graph_c1_429066 significant enrichment in the cellular process, single- were highly expressed in lungs of larvae. Unconserva- organism process, metabolic process, biological regulation, tive_c60430.graph_c0_428023 subsequently up-regulated cell, cell part, organelle, binding and catalytic activity the expressions of secretory phospholipase A2 receptor (Additional file 9). These targets were also aligned to the isoform X2 gene, which is functionally associated with KEGG database: a total of 203 differentially expressional oxidative stress-induced premature senescence. Two miRNA target genes were involved in 185 pathways, with down-regulated genes were pulmonary surfactant- ‘focal adhesion’ and ‘MAPK signaling pathway’ being the associated protein B precursor and pulmonary most abundant terms (Additional file 10). surfactant-associated protein B-like isoform X1, which are necessary for vertebrate lung development [56]. Discussion These results are consistent with Eo et al.’s report [49], Our understanding of the mechanisms involved in the and illustrate that miRNAs contribute to the observed transition from gill respiration to lung respiration re- morpho-functional adjustments. In amphibians, at least mains incomplete. In previous reports, candidate pro- 30% of oxygen exchange and most of the carbon dioxide teins [14], nitric oxide synthase (NOS)/nitric oxide [9], elimination occurs through the skin [57]. The skin of and conserved molecular networks [53] were reported to lungless salamanders even serves as their only respira- be involved in this complex process. Herein, we focused tory organ [58]. Among the five genes associated with on the expression profiles of miRNAs. In order to asso- respiration identified in Chinese giant salamander [27], ciate the miRNAs with their target genes, a reference NADH-ubiquinone oxidoreductase chain 5, ATPase sub- transcriptome of the Chinese giant salamander (A. davi- unit 6 and cytochrome b were differentially expressed be- dianus) was generated. As a result, along with already tween skin and other tissues in the present study. Huang assembled transcriptomes of the larval-stage tiger sala- et al. [59] found that one of the most highly expressed mander [49] and Chinese salamander [54], scientists genes related to the respiration in the skin of seven dif- aiming to study the expression of genes in the transition ferent anurans was NADH-ubiquinone oxidoreductase from aquatic to aerial respiration now have transcrip- (GO term: “catalytic activity”), which is consistent with tomes of three salamander species at disposal. our result. Therefore, all of these results illustrate that In the present study, miRNA expression profiles were gills, lungs and skin play an important role in the detected in three different tissues. Characteristics of the morpho-functional adjustments associated with respir- miRNAs detected, in particular 21 nt as the most com- ation transition via oxidative stress, lung development mon length and a significant bias towards U in the first and catalytic activity mechanisms. nucleotide, are considered to be typical for miRNAs Secondly, biological functions of miRNAs and their [55]. Based on these identified miRNAs, two steps were target genes involved in respiration regulation were clas- taken to better understand the association between miR- sified. Five known miRNAs were related to hypoxia- NAs and their target genes during the transition from induced respiration regulation. The let-7 miRNA family aquatic to aerial respiration. members are induced by hypoxia [60], and increased Firstly, we focused on the specific expression of miR- level of let-7a-5p targets major pulmonary arterial NAs in different tissues and different functions of these hypertension pathways [61]. Higher expression level of tissues during the respiration transition process. Expres- let-7a-5p, and its target gene sal-like protein 1 (SALL1), sion profiles of gill and skin tissues could be distin- in gills in comparison to lungs illustrated that gill respir- guished from the expression profile of lung tissue using ation might be a reflection hypoxia induced by the 115 miRNAs. Among these, we validated and confirmed change of environment: from oxygen-abundant terres- the differential expression of 12 miRNAs. Among them, trial to oxygen-scarce aquatic environment. Some stud- three new miRNAs differentially expressed between ies have reported that sal-like protein 1 plays essential lungs of larvae and adults show specific expression pat- roles in maintaining self-renewal and pluripotency of terns: unconservative_c60430.graph_c0_428023 was embryonic stem cells (ESCs) [62] and promotes the highly expressed in lungs of adults compared to lungs of intrahepatic cholangiocarcinoma cell migration [63]. larvae, while both unconservative_c41959.graph_c0_ These results indicate that cell differentiation and Su et al. BMC Genomics (2018) 19:406 Page 10 of 15 Fig. 3 Validation of differentially expressed miRNAs in the respiration tissue transfer process. Different letters above bars indicate that two means are significantly different (P < 0.05) migration may be one possible way of promoting adjust- induced apoptosis and fibrosis of cardiomyocytes. For ments in the respiration process by sal-like protein 1. miR-142, Lu et al. [65] found that hypoxia-induced cell Other three miRNAs have been associated with the in- proliferation and invasion could be inhibited by the hibition of hypoxia-induced adaptive changes: Wang overexpression of miR-142 or hypoxia inducible factor1 et al. [64] reported that miR-142-3-P inhibits hypoxia- (HIF-1α). Down-regulation of pulmonary miR-223 and Su et al. BMC Genomics (2018) 19:406 Page 11 of 15 Fig. 4 Validation of differentially expressed mRNAs in the respiration tissue transfer process. Different letters above bars indicate that two means are significantly different (P < 0.05). Note: c82457.graph_c0: sal-like protein 1; c54346.graph_c0: Semialdehyde dehydrogenase, NAD binding domain; c63687.graph_c2: sorting nexin-33; c53027.graph_c0: Ribosomal protein L11; c41045.graph_c0: Carbon-nitrogen hydrolase; c40974.graph_c0: Iron/zinc purple acid phosphatase-like protein C; c59945.graph_c0: glypican-5-like isoform X3; c60274.graph_c0: reverse transcriptase; c71119.graph_c0: centrosomal protein of 57 kDa isoform X1; c76183.graph_c0: trimethyllysine dioxygenase, mitochondrial isoform X3 upregulation of insulin-like growth factor 1 receptor ex- metabolism, as well as ribosomal structure and biogen- pression can be initiated by pulmonary hypoxia in the esis, this illustrates that metabolic ways were coordi- right heart failure [66]. In our study, the three miRNAs nated with adaptive changes in the skin respiratory mentioned above exhibited higher expression in lungs patterns. For two miRNAs, unconservative_c57205. than in skin and gills, which indicates that pulmonary graph_c0_339408 and unconservative_c62709.graph_c3_ breathing is susceptible to adaptive air respiration by 520410, we did not observe a significant difference in inhibiting the hypoxia-induced gene expression changes. expression among the studied tissues, but one of their However, miR-142 and miR-223 targets, Ribosomal pro- targeted genes exhibited a significant difference in tein L11 and Iron/zinc purple acid phosphatase-like pro- expression: trimethyllysine dioxygenase gene (validated tein C respectively, had higher expression in skin. As by qPCR) had higher expression in skin compared to they are involved in carbohydrate transport and other tissues. As its mitochondrial isoform X3 is known Su et al. BMC Genomics (2018) 19:406 Page 12 of 15 Fig. 5 miRNA-Gene-Network for the respiratory transition in Chinese giant salamander. The miRNA-Gene-Network was built using gene expression data and predicted interactions in the TargetScan miRNA database. The squares represent miRNAs and the circles represent their target genes, where red indicates up-regulation, and blue - down-regulation, purple – no significant regulation. Group1 is skin5 m and Group2 is gill5m, Group3 is lungs5 m and Group4 is lungs20m. There are four pairwise group comparisons (skin5 m vs gill5m, gill5 m vs lungs5m, gill5 m vs lungs20m and lungs5 m vs lungs20m) for differential miRNA expression regulation pattern analysis. So different proportions of different colors in squares and cycles correspond to the differential expression levels among the four groups (pairwise comparisons); for example, if half of the square is filled with red, it means that this miRNA was up-regulated in two of the four groups to be down-regulated in response to hypoxia (10% O2) components, qualitative and quantitative (immuno- [67], this result illustrates that these two miRNAs facili- fluorescence microscopy and Western blotting) pat- tate the respiratory system transition from water to land terns of hypoxia inducible factor1 (HIF1) expression by regulating the oxygen-related genes in lungs and skin. correlated inversely to that of nitric oxide synthase in Eo et al. [49] used larval tiger salamander and pyrose- both gills and lungs [9]. Synaptotagmin-like 4, associ- quencing to identify genes differentially expressed be- ated with the GO terms “phospholipid binding” and tween gill and lung tissues. They also found that “neurexin family protein binding”, is involved in intra- ‘pulmonary-associated protein isoform cra_a’ was exclu- cellular trafficking, secretion, and vesicular transport. sively expressed in the lungs. Among these differentially This gene indicates that lipid metabolism and nervous expressed miRNA-targeted genes, Carbon-nitrogen system may be related with the respiration transition. hydrolase (annotation in pFam) is associated with the Zabelinskii et al. [68] studied the lipid and fatty acid “nitrogen compound metabolic process” GO term. This composition of gills and lungs in 18 fish species and gene was involved in signal molecule changes in the gills seven mammalian species, and found different molar and lungs of African lungfish during the maintenance ratios between phospholipids and cholesterol: 2:1 in and arousal phases of aestivation [9]. They also reported fish and 3:1 in mammals. In vitro, cranial nerve (CN) localization and expression of nitric oxide synthase roots and spinal nerve (SN) roots had different fre- (NOS), Protein Kinase B, phospho-Akt (Akt), heat shock quency, and their amplitude bursts patterns corre- protein 90 (Hsp-90), hypoxia inducible factor α (HIF-1α) sponded to gill and lung ventilatory burst patterns, in a tissue-specific manner in parallel with organ re- respectively [69]. So, lipid metabolism and neural ac- adjustment in the gills and lungs of P. annectens during tivity also play a key role in the process of morpho- aestivation and arousal [9]. Among these molecular functional adjustments. Su et al. BMC Genomics (2018) 19:406 Page 13 of 15 Conclusions Additional file 9: Figure S8f. Go classification of differentially expressed This study provides the first large-scale report of miRNA targets for the miRNA-gene-network. Differentially expressed miRNA targets were allocated into three main GO categories: biological miRNA expression profiles and their predicted target process, cellular component and molecular function. Left vertical axis genes during the respiration transition from gills to represents the percentage of genes, while right vertical axis indicates lungs in Chinese giant salamander. A high-coverage the number of genes. (PNG 1185 kb) reference transcriptome was also generated and as- Additional file 10: Figure S9f. KEGG pathway analysis of differentially expressed miRNA targets for the miRNA-gene-network. The horizontal sembled. Two miRNAs with target genes related to axis represents the number of genes, and the percentage of this number oxygen sensing were differentially expressed between in the total number of genes in a given pathway. On the vertical axis are gill and lung tissues. Three miRNAs were differen- the names of KEGG pathways. (PNG 84 kb) tially expressed between lungs of larvae and lungs of adults. Differential expressions of several miRNAs were validated and confirmed; e.g., let-7a-5p was tar- Abbreviations geting major pulmonary arterial hypertension path- Akt: Phospho-Akt; CN: Cranial nerve; DEGs: Differentially expressed genes; EggNOG: Orthologous groups of genes; ESCs: Embryonic stem cells; FC: Fold ways. Their target gene annotation indicates that the change; FDR: False discovery rate; FPKM: Fragments per kb per million reads; observed miRNA regulation pattern is a response to gill5m: Larvae (5 months-old) gill; GO: Gene Ontology database; HIF-1α:Hypoxia the change in the levels of oxygen in the environ- inducible factor1; Hsp-90: Heat shock protein 90; KEGG: Kyoto Encyclopedia of Genes and Genomes database; lungs20m: Adults (20 months-old) lungs; ment. These results suggest that expression profiles of lungs5m: Larvae (5 months-old) lungs; miRNAs: MicroRNAs; MISA: MIcroSAtellite miRNAs in three respiratory tissues are correlated identification tool; NO: Nitric oxide; NOS: Nitric oxide synthase; PITA: Target with oxygen responses during the transition from accessibility; qPCR: Real time quantitative polymerase chain reaction; RNase H: M-MuLV Reverse Transcriptase; rRNA: Ribosomal RNAs; SALL1: Sal-like protein water to land. 1; SFTPA, SFTPC, and SFTPD: Pulmonary surfactants A, C, D; skin5m: Larvae (5 months-old) skin; SN: Spinal nerve; snoRNA: Small nucleolar RNAs; snRNA: Small nuclear RNAs; SSR: Microsatellite; TPM: Transcripts per million; tRNA: Transfer RNAs; U6: mamm U6 small nuclear 1 Additional files Acknowledgments We thank the students and staff of the Aquatic Genetic Laboratory, Freshwater Additional file 1: Supplementary Tables S1 to S11. (DOC 122 kb) Fisheries Research Center of Chinese Academy of Fishery Sciences for their kind Additional file 2: Figure S1. Size distribution of the assembled assistance in the study. The authors are grateful to Dr. Du for valuable technical unigenes. Horizontal axis gives different size intervals of the assembled support. unigenes, vertical axis gives the number of unigenes located in the specific size interval. (PNG 44 kb) Funding Additional file 3: Figure S2. Frequency distribution of the putative cSSRs This study was supported by grants from Central Public-interest Scientific observed. C is the number of SSRs present in compound formation; C* is Institution Basal Research Fund, Chinese Academy of Fishery Sciences the number of sequences containing more than one SSR; p1 – p5 represent (CAFS) (2016RC-LX03), New Aquaculture Project projects of Jiangsu Province the numbers of mono-, di-, tri-, tetra- and penta-nucleotides respectively. (Y2013–51) and the 2016 annual Jiangsu University “Qing Lan Project” young (PNG 21 kb) academic leaders training project (Jiangsu teachers [2016]15), the key Additional file 4: Figure S3. Distribution of the top BLASTX hits for technology programs of Zhenjiang (NY2016004) and the Key Projects in unigenes in the Nr databases. (PNG 64 kb) the National Science & Technology Pillar Program during the twelfth Additional file 5: Figure S4. GO classification of Chinese giant Five-Year Plan Period (2012BAD26B02), the National Nonprofit Institute salamander unigenes. 9236 (10.70%) unigenes were allocated into three Research Grant of CATAS-TCGRI (2013JBFM14). The funders had no role main GO categories: biological process, cellular component and in the design of the study, collection, analysis and interpretation of data, molecular function. Some of these unigenes were annotated with and in writing the manuscript. multiple GO terms. (PNG 109 kb) Additional file 6: Figure S5. Length distribution of novel miRNAs. Availability of data and materials Horizontal axis gives the length of novel miRNAs, vertical axis gives the The datasets generated during and/or analysed during the current study number of novel miRNAs. (PNG 110 kb) (assembled transcriptome and miRNA sequences) are available in NCBI’s Additional file 7: Figure S6. miRNA nucleotide bias. A: Nucleotide Transcriptome Shotgun Assembly (TSA) database (ftp://ftp.ncbi.nlm.nih.gov/ distribution at the first position of miRNA. Each row represents a genbank/tsa) under the accession numbers: PRJNA401072 and PRJNA398145 miRNA of a different length, and each column represents the (SRR5944935, SRR5944934, SRR5944939, SRR5944938, SRR5944937, distribution of bases in the first nucleotide position. B: Overall SRR5944936, SRR5944931, SRR5944930, SRR5944933, SRR5944932, nucleotide distribution. Each row represents a miRNA of a different SRR5944929, SRR5944928). length, and each column represents the overall distribution of the bases. (PNG 31 kb) Authors’ contributions Additional file 8: Figure S7. Volcano plots of differentially expressed SS, CJ, XJ and XP designed the experiments, supervised the study and miRNA for different respiration-related tissues. Red dots represent significantly modified the manuscript. SS, WY, WH, ZY, HW and HC performed the up-regulated genes and green dots significantly down-regulated genes in research, analyzed the data and wrote the manuscript. YX participated in the pairwise tissue comparisons (blue dots are genes with unchanged expression). design of the study, lab work, data analysis, data visualisation, molecular lab The x-axis represents the log2-transformed gene expression. The y-axis is the p work; SS carried out sequence alignments and contributed custom–made value (−log2) adjusted by Benjamini-Hochberg correction. All pairwise tissue software and codes; CJ and HW participated in sample collection; XJ comparisons revealed differentially expressed miRNAs. A: skin5 m vs gill5m; B: conceived and coordinated the study; XP co-wrote the manuscript, and skin5 m vs lungs5m; C: skin5 m vs lungs20m; D: gill5 m vs lungs5m; E: gill5 m reviewed the manuscript. All authors contributed comments to the final vs lungs20m; F: lungs5 m vs lungs20m. (PNG 742 kb) version of the manuscript, read and approved it. Su et al. BMC Genomics (2018) 19:406 Page 14 of 15 Ethics approval and consent to participate 15. Vila-Casadesus M, Gironella M, Lozano JJ. MiRComb: an R package to Study design and all experimental procedures involving animals were analyse miRNA-mRNA interactions. Examples across five digestive cancers. reviewed, approved and supervised by the Institutional Animal Care and Use PLoS One. 2016;11:e0151127. Committee of the Freshwater Fisheries Research Center of Chinese Academy 16. Xu Y, Chu L, Jin Q, Wang Y, Chen X, Zhao H, Xue Z. Transcriptome-wide of Fishery Sciences. Various measures were taken to minimize the suffering identification of miRNAs and their targets from Typha angustifolia by RNA- of the animals. Seq and their response to cadmium stress. PLoS One. 2015;10:e0125462. 17. Liu X, Zhu L, Liao S, Xu Z, Zhou Y. The porcine microRNA transcriptome Competing interests response to transmissible gastroenteritis virus infection. PLoS One. 2015;10: The authors declare that they have no competing interests. e0120377. 18. Coolen M, Bally-Cuif L. MicroRNAs in brain development and physiology. Curr Opin Neurobiol. 2009;19:461–70. Publisher’sNote 19. Ivey KN, Srivastava D. MicroRNAs as regulators of differentiation and cell fate Springer Nature remains neutral with regard to jurisdictional claims in decisions. Cell Stem Cell. 2010;7:36–41. published maps and institutional affiliations. 20. Bailey SG, Sanchez-Elsner T, Stephanou A, Cragg MS, Townsend PA. Regulating the genome surveillance system: miRNAs and the p53 super Author details family. Apoptosis. 2010;15:541–52. Key Laboratory of Genetic Breeding and Aquaculture Biology of Freshwater 21. Dykxhoorn DM. MicroRNAs and metastasis: little RNAs go a long way. Fishes, Ministry of Agriculture; Freshwater Fisheries Research Center, Chinese Cancer Res. 2010;70:6401–6. Academy of Fishery Sciences, Wuxi 214081, People’s Republic of China. 22. Skalsky RL, Cullen BR. Viruses, microRNAs, and host interactions. Annu Rev Department of Animal Husbandry & Veterinary Medicine, Jiangsu Microbiol. 2010;64:123–41. Polytechnic College of Agriculture and Forestry, Zhenjiang 212400, People’s 23. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. Republic of China. Wuxi Fisheries College, Nanjing Agricultural University, 2004;116:281–97. Wuxi 214081, People’s Republic of China. 24. Nana-Sinkam SP, Karsies T, Riscili B, Ezzie M, Piper M. Lung microRNA: from development to disease. Expert Rev Respir Med. 2009;3:373–85. Received: 4 January 2018 Accepted: 12 April 2018 25. Rupani H, Sanchez-Elsner T, Howarth P. MicroRNAs and respiratory diseases. Eur Respir J. 2013;41:695–705. 26. Nowak JS, Michlewski G. miRNAs in development and pathogenesis of the References nervous system. Biochem Soc Trans. 2013;41:815–20. 1. Wang X-M, Zhang K-J, Wang Z-H, Ding Y-Z, Wu W, Huang S. The decline of 27. Geng X, Wei H, Shang H, Zhou M, Chen B, Zhang F, Zang X, Li P, Sun J, Che the Chinese giant salamander Andrias davidianus and implications for its J, et al. Proteomic analysis of the skin of Chinese giant salamander (Andrias conservation. Oryx. 2004;38:197–202. davidianus). J Proteome. 2015;119:196–208. 2. Zhang L, Jiang W, Wang Q-J, Zhao H, Zhang H-X, Marcec RM, Willard ST, 28. Chen L, Fan J, Hu L, Hu Z, Xie Y, Zhang Y, Lou Y, Nevo E, Fu J. A Kouba AJ. Reintroduction and post-release survival of a living fossil: the transcriptomic analysis of bermudagrass (Cynodon dactylon) provides novel Chinese Giant salamander. PLoS One. 2016;11:e0156715. insights into the basis of low temperature tolerance. BMC Plant Biol. 2015; 3. Simons RS, Bennett WO, Brainerd EL. Mechanics of lung ventilation in a 15:216. post-metamorphic salamander, Ambystoma Tigrinum. J Exp Biol. 2000; 29. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis 203:1081–92. X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly 4. Ijspeert AJ, Crespi A, Ryczko D, Cabelguen JM. From swimming to from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29: walking with a salamander robot driven by a spinal cord model. Sci. 644–52. 2007;315:1416–20. 30. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ, 5. Piiper J. Respiratory gas exchange at lungs, gills and tissues: mechanisms Gapped BLAST. PSI-BLAST: a new generation of protein database search and adjustments. J Exp Biol. 1982;100:5–22. programs. Nucleic Acids Res. 1997;25:3389–402. 6. Burggren W, Doyle M. Ontogeny of regulation of gill and lung ventilation in 31. Finn RD, Coggill P, Eberhardt RY, Eddy SR. The Pfam protein families the bullfrog, Rana catesbeiana. Respir Physiol. 1986;66:279–91. database: towards a more sustainable future. Nucleic Acids Res. 2016;44: 7. Gdovin MJ, Torgerson CS, Remmers JE. Neurorespiratory pattern of gill and D279–85. lung ventilation in the decerebrate spontaneously breathing tadpole. Respir 32. Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walter MC, Rattei Physiol. 1998;113:135–46. T, Mende DR, Sunagawa S, Kuhn M, et al. eggNOG 4.5: a hierarchical 8. Lenfant C, Johansen K. Gas exchange in gill, skin, and lung breathing. Respir orthology framework with improved functional annotations for eukaryotic, Physiol. 1972;14:211–8. prokaryotic and viral sequences. Nucleic Acids Res. 2016;44:D286–93. 9. Garofalo F, Amelio D, Icardo JM, Chew SF, Tota B, Cerra MC, Ip YK. Signal 33. Bairoch A, Apweiler R. The SWISS-PROT protein sequence database and its molecule changes in the gills and lungs of the African lungfish Protopterus supplement TrEMBL in 2000. Nucleic Acids Res. 2000;28:45–8. annectens, during the maintenance and arousal phases of aestivation. Nitric 34. Finn RD, Clements J, Arndt W, Miller BL, Wheeler TJ, Schreiber F, Oxide. 2015;44:71–80. Bateman A, Eddy SR. HMMER web server: 2015 update. Nucleic Acids 10. Hu XG, Liu H, Jin Y, Sun YQ, Li Y, Zhao W, El-Kassaby YA, Wang XR, Mao JF. Res. 2015;43:W30–8. De novo transcriptome assembly and characterization for the widespread 35. Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto and stress-tolerant conifer Platycladus orientalis. PLoS One. 2016;11: encyclopedia of genes and genomes. Nucleic Acids Res. 1999;27:29–34. e0148985. 36. Chen L, Chu C, Lu J, Kong X, Huang T, Cai YD. Gene ontology and KEGG 11. Xiong H, Li Q, Liu S, Wang F, Xiong Z, Chen J, Chen H, Yang Y, Tan X, Luo pathway enrichment analysis of a drug target-based classification system. Q, et al. Integrated microRNA and mRNA transcriptome sequencing reveals PLoS One. 2015;10:e0126492. the potential roles of miRNAs in stage I endometrioid endometrial 37. Thiel T, Michalek W, Varshney RK, Graner A. Exploiting EST databases for the carcinoma. PLoS One. 2014;9:e110163. development and characterization of gene-derived SSR-markers in barley 12. Vidal EA, Moyano TC, Krouk G, Katari MS, Tanurdzic M, McCombie WR, (Hordeum vulgare L.). Theor Appl Genet. 2003;106:411–22. Coruzzi GM, Gutierrez RA. Integrated RNA-seq and sRNA-seq analysis 38. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data identifies novel nitrate-responsive genes in Arabidopsis thaliana roots. BMC with or without a reference genome. BMC Bioinf. 2011;12:323. Genomics. 2013;14:701. 39. Anders S, Huber W. Differential expression analysis for sequence count data. 13. Wang Y, Brahmakshatriya V, Lupiani B, Reddy SM, Soibam B, Benham AL, Gunaratne P, Liu HC, Trakooljul N, Ing N, et al. Integrated analysis of Genome Biol. 2010;11:R106. microRNA expression and mRNA transcriptome in lungs of avian influenza 40. Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery virus infected broilers. BMC Genomics. 2012;13:278. rate in behavior genetics research. Behav Brain Res. 2001;125:279–84. 14. Biscotti MA, Gerdol M, Canapa A, Forconi M, Olmo E, Pallavicini A, Barucca 41. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 M, Schartl M. The lungfish transcriptome: a glimpse into molecular accurately identifies known and hundreds of novel microRNA genes in evolution events at the transition from water to land. Sci Rep. 2016;6:21571. seven animal clades. Nucleic Acids Res. 2012;40:37–52. Su et al. BMC Genomics (2018) 19:406 Page 15 of 15 42. Bonnet E, Wuyts J, Rouze P, Van de Peer Y. Evidence that microRNA 64. Wang Y, Ouyang M, Wang Q, Jian Z. MicroRNA-142-3p inhibits hypoxia/ precursors, unlike other non-coding RNAs, have lower folding free energies reoxygenationinduced apoptosis and fibrosis of cardiomyocytes by than random sequences. Bioinformatics. 2004;20:2911–7. targeting high mobility group box 1. Int J Mol Med. 2016;38:1377–86. 43. Lin YL, Ma LT, Lee YR, Lin SS, Wang SY, Chang TT, Shaw JF, Li WH, Chu FH. 65. Lu Y, Ji N, Wei W. MiR-142 modulates human pancreatic cancer proliferation MicroRNA-like small RNAs prediction in the development of Antrodia and invasion by targeting hypoxia-inducible factor 1 (HIF-1alpha) in the cinnamomea. PLoS One. 2015;10:e0123245. tumor microenvironments. Biol Open. 2017;6:252–9. 66. Shi L, Kojonazarov B, Elgheznawy A, Popp R, Dahal BK, Bohm M, Pullamsetti 44. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. The role of site SS, Ghofrani HA, Godecke A, Jungmann A, et al. miR-223-IGF-IR signalling in accessibility in microRNA target recognition. Nat Genet. 2007;39: hypoxia- and load-induced right-ventricular failure: a novel therapeutic 1278–84. approach. Cardiovasc Res. 2016;111:184–93. 45. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using 67. Ganfornina MD, Perez-Garcia MT, Gutierrez G, Miguel-Velado E, Lopez-Lopez real-time quantitative PCR and the 2(−Delta Delta C (T)) method. Methods. JR, Marin A, Sanchez D, Gonzalez C. Comparative gene expression profile of 2001;25:402–8. mouse carotid body and adrenal medulla under physiological hypoxia. 46. Wang X, Li M, Wang Z, Han S, Tang X, Ge Y, Zhou L, Zhou C, Yuan Q, Yang J Physiol. 2005;566:491–503. M. Silencing of long noncoding RNA MALAT1 by miR-101 and miR-217 68. Zabelinskii SA, Brovtsyna NB, Chebotareva MA, Gorbunova OB, Krivchenko inhibits proliferation, migration, and invasion of esophageal squamous cell AI. Comparative investigation of lipid and fatty acid composition of fish gills carcinoma cells. J Biol Chem. 2015;290:3925–35. and mammalian lungs. A model of the membrane lipid component areas. 47. Li F, Wang L, Lan Q, Yang H, Li Y, Liu X, Yang Z. RNA-Seq analysis and gene Comp Biochem Physiol B Biochem Mol Biol. 1995;111:127–40. discovery of Andrias davidianus using Illumina short read sequencing. PLoS 69. Torgerson CS, Gdovin MJ, Remmers JE. Fictive gill and lung ventilation in One. 2015;10:e0123730. the pre- and postmetamorphic tadpole brain stem. J Neurophysiol. 1998;80: 48. Huang Y, Xiong JL, Gao XC, Sun XH. Transcriptome analysis of the Chinese 2015–22. giant salamander (Andrias davidianus) using RNA-sequencing. Genom Data. 2017;14:126–31. 49. Eo SH, Doyle JM, Hale MC, Marra NJ, Ruhl JD, DeWoody JA. Comparative transcriptomics and gene expression in larval tiger salamander (Ambystoma tigrinum) gill and lung tissues as revealed by pyrosequencing. Gene. 2012; 492:329–38. 50. Ji Z, Wang G, Xie Z, Zhang C, Wang J. Identification and characterization of microRNA in the dairy goat (Capra hircus) mammary gland by Solexa deep- sequencing technology. Mol Biol Rep. 2012;39:9361–71. 51. Huang Y, Yang YB, Gao XC, Ren HT, Sun XH. Identification and characterization of the Chinese giant salamander (Andrias davidianus) miRNAs by deep sequencing and predication of their targets. 3 Biotech. 2017;7:235. 52. Li M, Xia Y, Gu Y, Zhang K, Lang Q, Chen L, Guan J, Luo Z, Chen H, Li Y, et al. MicroRNAome of porcine pre- and postnatal development. PLoS One. 2010;5:e11541. 53. Cushing L, Jiang Z, Kuang P, Lu J. The roles of microRNAs and protein components of the microRNA pathway in lung development and diseases. Am J Respir Cell Mol Biol. 2015;52:397–408. 54. Che R, Sun Y, Wang R, Xu T. Transcriptomic analysis of endangered Chinese salamander: identification of immune, sex and reproduction-related genes and genetic markers. PLoS One. 2014;9:e87940. 55. Pontes O, Costa-Nunes P, Vithayathil P, Pikaard CS. RNA polymerase V functions in Arabidopsis interphase heterochromatin organization independently of the 24-nt siRNA-directed DNA methylation pathway. Mol Plant. 2009;2:700–10. 56. Hyatt BA, Resnik ER, Johnson NS, Lohr JL, Cornfield DN. Lung specific developmental expression of the Xenopus laevis surfactant protein C and B genes. Gene Expr Patterns. 2007;7:8–14. 57. Jorgensen CB. Amphibian respiration and olfaction and their relationships: from Robert Townson (1794) to the present. Biol Rev Camb Philos Soc. 2000;75:297–345. 58. Nussbaum RA, Wilkinson M, New Genus A. Of lungless tetrapod: a radically divergent caecilian (Amphibia: Gymnophiona). Proc R Soc B Biol Sci. 1995; 261:331–5. 59. Huang L, Li J, Anboukaria H, Luo Z, Zhao M, Wu H. Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Rep. 2016;6:24069. 60. Chen Z, Lai TC, Jan YH, Lin FM, Wang WC, Xiao H, Wang YT, Sun W, Cui X, Li YS, et al. Hypoxia-responsive miRNAs target argonaute 1 to promote angiogenesis. J Clin Invest. 2013;123:1057–67. 61. Wu D, Talbot CC Jr, Liu Q, Jing ZC, Damico RL, Tuder R, Barnes KC, Hassoun PM, Gao L. Identifying microRNAs targeting Wnt/beta-catenin pathway in end-stage idiopathic pulmonary arterial hypertension. J Mol Med (Berl). 2016;94:875–85. 62. Zhang X, Yuan X, Zhu W, Qian H, Xu W. SALL4: an emerging cancer biomarker and target. Cancer Lett. 2015;357:55–62. 63. Zhu L, Huang F, Deng G, Nie W, Huang W, Xu H, Zheng S, Yi Z, Wan T. Knockdown of Sall4 inhibits intrahepatic cholangiocarcinoma cell migration and invasion in ICC-9810 cells. Onco Targets Ther. 2016;9: 5297–305. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Genomics Springer Journals

Comparative expression analysis identifies the respiratory transition-related miRNAs and their target genes in tissues of metamorphosing Chinese giant salamander (Andrias davidianus)

Free
15 pages
Loading next page...
 
/lp/springer_journal/comparative-expression-analysis-identifies-the-respiratory-transition-VH04s10jh7
Publisher
BioMed Central
Copyright
Copyright © 2018 by The Author(s).
Subject
Life Sciences; Life Sciences, general; Microarrays; Proteomics; Animal Genetics and Genomics; Microbial Genetics and Genomics; Plant Genetics and Genomics
eISSN
1471-2164
D.O.I.
10.1186/s12864-018-4662-5
Publisher site
See Article on Publisher Site

Abstract

Background: Chinese giant salamander (Andrias davidianus) undergoes a metamorphosis from aquatic larvae to terrestrial adults, with concomitant transfer of respiration from gills to lungs prior to metamorphosis. These two tissues, as well as skin, were sampled to identify the differentially expressed miRNAs. Results: High-coverage reference transcriptome was generated from combined gill, lung and skin tissues of metamorphosing juveniles, and lung tissue of adults: 86,282 unigenes with total length of approximately 77,275,634 bp and N50 of 1732 bp were obtained. Among these, 13,246 unigenes were assigned to 288 pathways. To determine the possible involvement of miRNAs in the respiratory transition, small RNA libraries were sequenced; 282 miRNAs were identified, 65 among which were known and 217 novel. Based on the hierarchical clustering analysis, the twelve studied samples were classified into three major clusters using differentially expressed miRNAs. We have validated ten differentially expressed miRNAs and some of their related target genes using qPCR. These results largely corroborated the results of transcriptomic and miRNA analyses. Finally, an miRNA-gene-network was constructed. Among them, two miRNAs with target genes related to oxygen sensing were differentially expressed between gill and lung tissues. Three miRNAs were differentially expressed between the lungs of larvae and lungs of adults. Conclusions: This study provides the first large-scale miRNA expression profile overview during the respiration transition from gills to lungs in Chinese giant salamander. Five differentially expressed miRNAs and their target genes were identified among skin, gill and lung tissues. These results suggest that miRNA profiles in respiratory tissues play an important role in the regulation of respiratory transition. Keywords: Deep sequencing, Respiratory system, miRNA, Andrias davidianus * Correspondence: 362205379@qq.com; yuanxh@ffrc.cn Shengyan Su, Yuheng Wang and Huiwei Wang contributed equally to this work. Department of Animal Husbandry & Veterinary Medicine, Jiangsu Polytechnic College of Agriculture and Forestry, Zhenjiang 212400, People’s Republic of China Key Laboratory of Genetic Breeding and Aquaculture Biology of Freshwater Fishes, Ministry of Agriculture; Freshwater Fisheries Research Center, Chinese Academy of Fishery Sciences, Wuxi 214081, People’s Republic of China Full list of author information is available at the end of the article © The Author(s). 2018 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. Su et al. BMC Genomics (2018) 19:406 Page 2 of 15 Background pattern of the expression of these target mRNAs, micro- The endangered Chinese giant salamander (Andrias RNAs (miRNAs) play an important role by specifically davidianus, Cryptobranchidae family), endemic to binding to the 3’UTR of mRNA and (usually) promoting China, is the world’s largest amphibian. It has been mRNA degradation [15]. miRNAs are endogenous ~ 22 called a “living fossil”, and it is referred to locally as nucleotide-long non-coding RNAs, which regulate a wawayu (baby fish) because its vocalization is somewhat range of essential cellular and biological processes by resemblant of a baby’s cry [1]. The Chinese government targeting mRNA for cleavage or translational repression has declared this species a class II protected species [2]. [16–23]. There are indications that miRNAs modulate Its ability to switch from aquatic to aerial respiration pulmonary development and maintain lung homeostasis and the associated morpho-functional adjustments [3, 4] [24, 25]. Misexpression of some miRNAs can contribute further contribute to its status of a species of great inter- towards neurodevelopmental disorders, which can also est to biologists. Studies of this species might provide affect the lung ventilation [26]. key information about body adaptions (anatomical, With two key developmental stages, Chinese giant physiological and molecular) that enable organisms to salamander is a typical amphibian, undergoing a transi- transition from water to land. tion from aquatic larvae to terrestrial adults, with con- Piiper (1982) divided the transition from water to land comitant transfer of respiration from gills to lungs prior into three sub-stages dependent on different organs: fish to the metamorphosis [27]. In order to unveil the role of gills, amphibian gills, and lungs and skin. The features of miRNAs in the adjustment of respiratory system to the respiration in transitional states showed that the adjust- transition from water to land, we have studied the ex- ment for gas exchange style in different breathing pression of miRNAs and their corresponding target models is related to the metabolic rate [5]. Burggren and genes in metamorphosing larvae of the Chinese giant Doyle (1986) reported that the transition from gill res- salamander. Real time quantitative polymerase chain re- piration to lung respiration during the development of action (qPCR) was used to study expression in gill, lung bullfrog larvae is accompanied by a decline in the regu- and skin tissues. We also assembled a high-coverage ref- lation of branchial ventilation frequency reflex and at- erence transcriptome, and determined the putative regu- tendant increasing dependence on breathing air [6]. latory network of miRNAs associated with the studied Based on the neurorespiratory pattern of gill and lung morpho-functional adjustments. This study contributes ventilation, physiologists found that spinal nerve II is to the understanding of the role of miRNA-mediated crucial for lung ventilation [7]. However, these studies regulatory networks in the transition from water to land merely identified the differences among different respir- of Chinese giant salamander. ation styles. In transitional forms, depending on both gill and lung ventilation, gill and lung breathing can be re- Methods placed or complemented by skin and or buccopharyn- Sample collection and RNA isolation geal breathing [8]. This adaptation is necessary for the Healthy Chinese giant salamander’s larvae (5 months-old) adaptive changes of respiratory properties of blood and and adults (20 months-old) were obtained from the ventilatory and circulatory flow rates [8]. Nitric oxide Jurong Chinese Giant Salamander Breeding Center in synthase (NOS)/nitric oxide (NO) levels in gills and Zhenjiang, Jiangsu, China. Specimens were anaesthetized lungs of Protopterus annectens during aestivation and using 2-phenoxyethanol, euthanized using tricaine metha- arousal changed in a tissue-specific manner in parallel nesulfonate (Sigma, USA), and promptly dissected. Three with organ readjustment [9]. However, despite numerous different tissues were collected from three larvae: skin efforts, signal transduction mechanisms that underlie (skin5m), gills (gill5m), and lungs (lungs5m); and only the cyclic branchial and pulmonary remodeling remain lungs from three adults (lungs20m) (Fig.1;Additional file 1: incompletely explained. Assunta et al. (2016) studied the Table S1). This ensured triplicate samples for four differ- water-land transition from the perspective of gene ex- ent tissue groups (skin5m, gill5m, lungs5m, lungs20m). pression patterns; they used lungfish transcriptome se- RNA was extracted from all samples using Trizol Reagent quence (which is a time-effective and economical (Invitrogen, USA) according to the manufacturer’sinstruc- method [10–13]) and phylogenetic relationship analysis tions. Purity of RNA was checked by NanoPhotometer® [14]. They identified 226 concatenated protein sequences spectrophotometer (IMPLEN, CA, USA). The final con- in all vertebrates and 59, 951 informative amino acid po- centration was determined using Qubit® RNA Assay Kit in sitions in both lungfish/tetrapod and the coelacanth/ Qubit®2.0 Fluorometer (Life Technologies, CA, and USA). lungfish + tetrapod nodes. They also identified several The integrity of RNA was assessed using the RNA Nano proteins in the tetrapod lung which were not found in 6000 Assay Kit and Agilent Bioanalyzer 2100 system the actinopterygian fish (pulmonary surfactants A, C, D (Agilent Technologies, CA, and USA). The handling of SFTPA, SFTPC, and SFTPD). Regarding the regulation animals was conducted in accordance with the guidelines Su et al. BMC Genomics (2018) 19:406 Page 3 of 15 Fig. 1 Chinese giant salamander respiratory system transition from gills, lungs and skin to respiration mainly using lungs (a) A transitioning Chinese giant salamander with gills and lungs; (b) Chinese giant salamander with lungs and without gills for the care and use of animals for scientific purposes set cBot Cluster Generation System using TruSeq PE Cluster by the Institutional Animal Care and Use Committee of the Kit v3-cBot-HS (Illumina) according to the manufacturer’s Freshwater Fisheries Research Center, Chinese Academy of instructions. Finally, the libraries were sequenced using Fishery Sciences, Wuxi, China. paired-end reads on an Illumina Hiseq 4000 platform. Library preparation and transcriptome sequencing De novo transcriptome assembly and prediction of Sequencing libraries were constructed using 3 μgof coding sequence RNA per sample using NEBNext®Ultra™ RNA Library In order to obtain high-quality clean reads, the raw data Prep Kit for Illumina® (NEB, USA), following the manu- were trimmed to remove reads containing adapter, reads facturer’s instructions. Index codes were added to attri- containing ploy-N and low-quality reads. Q20, Q30, GC- bute sequences to each sample. Briefly: mRNA was content and sequence duplication level of the clean data purified from total RNA using poly-T oligo-attached were calculated. Clean data were then assembled de magnetic beads and fragmented using divalent cations novo as described before [28]: the sequenced ‘left’ and under elevated temperature in NEBNext First Strand ‘right’ files (read1 and read2 files respectively) from all Synthesis Reaction Buffer (5X). First strand cDNA was samples were pooled into left.fq and right.fq files, re- transcribed using random hexamer primer and M-MuLV spectively. Transcriptome assembly was finalized using Reverse Transcriptase (RNase H). Second strand cDNA these two large files, and the longest fragments per gene synthesis was then performed using DNA polymerase I were designated as unigenes. All these operations were and RNase H. Remaining overhangs were converted into conducted using Trinity [29], with all parameters set to blunt ends using exonuclease/polymerase. After adenyla- default. tion of 3′ ends of DNA fragments, NEBNext adaptors with hairpin loop structure were ligated to prepare for Gene function annotation, SSR detection and primer hybridization. Library fragments were purified with design AMPure XP system (Beckman Coulter, Beverly, USA) to These unigenes were then queried against public databases select cDNA fragments with the length of 150~ 250 bp. for functional annotation: NCBI non-redundant protein In order to generate the products, we used 3 μl USER database using BLASTX searches [30], Protein Families Enzyme (NEW ENGLAND BioLabs Inc., USA) with (Pfam) [31], orthologous groups of genes (EggNOG) [32] size-selected, adaptor-ligated cDNA at 37 °C for 15 min, and Swiss-Prot database [33]. The cut-off of value for anno- followed by 5 min at 95 °C. PCR was performed with tation in NR, Swissport, KEGG databases was set E to 1-e5 Phusion High-Fidelity DNA polymerase, universal PCR using blastx, whereas for the Pfam database the value was primers and index (X) primer. Ultimately, we purified set E to 10 using hmmscan. In addition, the predicted CDS PCR products with AMPure XP system and assessed the were classified into EggNOG categories using HMMER library quality on the Agilent Bioanalyzer 2100 system. version 3.1 with an E-value cutoff of 1E-05 [34]. Kyoto We performed clustering of the index-coded samples on a Encyclopedia of Genes and Genomes database (KEGG) Su et al. BMC Genomics (2018) 19:406 Page 4 of 15 [35] was used to associate unigenes with metabolic path- the Agilent Bioanalyzer 2100 system. We performed the ways. Gene Ontology database (GO) [36] was used to clas- clustering of index-coded samples on a cBot Cluster sify the unigenes functionally into biological process, Generation System using TruSeq PE Cluster Kit v4- molecular functions and cellular components categories. cBot-HS (Illumina) according to the manufacturer′s in- GO term analysis was conducted by blast2go using the Nr structions. After cluster generation, library preparations annotation information. TransDecoder was used to predict were sequenced on an Illumina HiSeq X Ten platform the sequence of the Unigene coding region and its corre- and single-end reads were generated. sponding amino acid sequence, which was recommended by Trinity and Cuffinks. MIcroSAtellite identification tool Small RNA annotation, new miRNAs, and prediction of (MISA) [37] was used to detect SSRs in the transcriptome their target genes with the following parameters: definition (unit size - min. For small RNA sequencing reads, all sequences with repeats) = 1–10 2–63–54–55–56–5, and interruption adapter or poly-N contaminants and low-quality reads (max. Difference between two SSRs) = 100. Six types of SSRs were flittered. After that, sequences < 15 and > 35 nt were investigated:mono-,di-,tri-, tetra-, penta-, and hexa- were removed. Q20, Q30, GC-content and sequence du- nucleotide repeats. Primers for SSRs were designed using plication level of the clean data were calculated as well. Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/). Using Bowtie software, the clean reads were aligned with Rfam, Silva, GtRNAdb and Repbase databases to filter Quantification of gene expression levels the transfer RNAs (tRNA), small nuclear RNAs (snRNA) The read count for each gene was obtained by mapping , small nucleolar RNAs (snoRNA), ribosomal RNAs clean reads to assembled transcriptome using RSEM (rRNA) and other ncRNAs, as well as repeats. By com- [38]. Expression levels of unigenes were calculated as paring the remaining reads with known miRNAs from fragments per kb per million reads (FPKM). Differen- the miRBase, known miRNAs were detected and new tially expressed genes (DEGs) were identified using miRNAs predicted. Non-conserved miRNAs were pre- DESeq R package (1.10.1) [39], employing a model based dicted using miRdeep2 [41] with the following run pa- on the negative binomial distribution. Benjamini and rameters: “-g -1 -b 0”. Secondary structure prediction of Hochberg’s approach was used to adjust the resulting P new miRNAs was performed using Randfold tools [42] value, which was determined by false discovery rate with the parameters: s = simple mononucleotide shuf- (FDR) and fold change (FC) [40]. In the pairwise analysis fling, and number of randomizations = 99. As we did not of the four groups (12 samples with 3 replicates) de- have the whole genome sequence at disposal, RNA-seq scribed before, any gene with an adjusted P < 0.05 and data were used to predict new microRNAs. FC ≥ 2 was considered to be a differentially expressed gene (DEG). miRNA expression levels and functional annotation of their target genes Library preparation and sRNA sequencing To estimate the miRNA expression levels for each sam- A total amount of 1.5 μg RNA per sample was used as ple, sRNAs were mapped back onto the precursor se- input material for the RNA sample preparations. Se- quence [43]. From the mapping results we know the quencing libraries were generated using NEBNext®Ultra™ readcount, which represents the number of reads that small RNA Sample Library Prep Kit for Illumina® (NEB, are mapped to a certain miRNA. Then the expression of USA) following manufacturer’s recommendations, and miRNA in each sample was normalized as transcripts index codes were added to attribute sequences to each per million (TPM): TPM = readcount× 1,000,000/ sample. First of all, 3’ SR Adaptor for Illumina was Mapped Reads. Potential target genes of differentially mixed with RNA and nuclease-free water. This was expressed miRNAs were identified using miRnada v3.3a followed by incubation for 2 min at 70 °C in a preheated (run parameters: -sc 150.0 -en − 30 -scale 4.0 -go − 2.0 thermal cycler. Tube was then transferred to ice, and 3’ -ge − 8.0) and RNA hybrid v2.1.1 (run parameters: -m Ligation Reaction Buffer (2X) and 3’ Ligation Enzyme 50,000 -d 1.9,0.28 -b 1 -e − 30) programs. Intersection of Mix. The mixture was incubated for 1 h at 25 °C in a the target set was regarded as the target gene of a micro- thermal cycler. To prevent adaptor-dimer formation, the RNA. After these step, intersection of predicted targets SR RT Primer hybridizes to the excess of 3′ SR Adaptor were verified by target accessibility (PITA) [44].Target that has remained free after the 3′ ligation reaction, and genes of miRNAs were subjected to BLASTX searches transforms the single-stranded DNA adaptor into a (E-value threshold of 1.00e-5) and annotated against the double-stranded DNA molecule. After ligating the 5′ SR NCBI’s Nr and Pfam databases, as well as the EggNOG Adaptor, reverse transcription and PCR amplification and Swiss-Prot databases again. Target genes were also were performed. At last, PCR products were purified subjected to GO analysis, and classified into biological (AMPure XP system) and library quality was assessed on process, cellular component and molecular function Su et al. BMC Genomics (2018) 19:406 Page 5 of 15 categories. Metabolic pathways were analyzed using the assemble all the trimmed reads. About 77.83% of clean KEGG database. reads were mapped to the newly assembled transcrip- tome. Finally, a total of 86,282 unigenes with total length Validation of differentially expressed miRNAs and their of approximately 77,275,634 bp and N50 of 1732 bp corresponding target genes were obtained (Table 1). These values are higher than Expression levels of the following randomly selected the number of unigenes and N50 in previously reported miRNAs were determined using predesigned TaqMan Chinese giant salamander transcriptomes [47, 48]. MicroRNA assay (Life Technologies, Foster City, CA, Among these, 30,060 unigenes (34.84%) were < 300 bp in USA): aca-miR-142-5p, bta-miR-142-5p, gga-miR-34c-3p, length, 33,789 (39.16%) were > 301 and < 1000 bp, and aca-let-7a-5p, aca-miR-203-3p, aca-miR-203-5p, aja-miR- 22,433 (26.00%) were > 1000 bp (Additional file 2). The 142, dre-miR-203a-5p, efu-miR-223 and ipu-miR-142. overall GC ratio of the unigenes was 49.19%. Most of cDNA was synthesized from 100 ng of total RNA using the CDS were 200~ 300 bp in length (19,941, 31.01%), TaqMan® MicroRNA Reverse Transcription kit (Life Tech- followed by 100~ 200 bp (15,984, 24%), and 300~ 400 bp nologies). qPCR was performed using 7900HT Fast Real- (7652, 11.9%) (as predicted by Trans-Decoder). No CDS Time PCR System (Life Technologies), with primers listed with length < 100 were found, but 1130 (1.76%) of CDS in Additional file 1: Tables S2 and S3. Relative expression were > 3000 bp (Additional file 1: Table S5). A total of of miRNAs was determined using the 2-ΔΔCT method 13,199 SSRs were identified from 22,433 unigenes (> 1 kb) [45], with mamm U6 small nuclear 1 (U6), which belongs using MISA Perl script, with 816 (3.64%) unigenes con- to the class of small nuclear RNAs, used as an en- taining more than one SSR. Among the different types of dogenous reference [46]. No significant differences in SSRs, mono-nucleotide repeats were the most abundant the U6 expression were found among the four groups (10,332), followed by di-nucleotide (1233), tri-nucleotide (Additional file 1: Tables S4). All reactions were per- (687), tetra-nucleotide (88) and penta-nucleotide repeats formed in triplicate [37], data expressed as mean ± SE, (4) (Additional file 1: Table S6). Frequency distribution of and subjected to one-way analysis of variance these putative cSSRs was also calculated (Additional file 3). (ANOVA) using SAS8.0 program. On average, one cSSR locus was found for about every 5.9 kb of a unigene sequence. Results BLASTX and BLASTn searches against GO, KEGG, Assembly and annotation of Andrias davidianus eggnog and Nr databases were conducted to validate transcriptome and annotate the assembled unigenes; this produced To obtain a reference transcriptome for the study of 6390 (21.27%), 9236 (30.74%), 13,246 (44.09%), 16,180 microRNA regulation patterns during the transition (53.85%), 27,181 (90.48%), 28,247 (94.02%) unigene from gill respiration to lung respiration, samples from matches respectively (E-value ≤1.0E-5) (Additional file 1: individuals which only had lung tissue and from those Table S7). A set of 30,044 unigenes were annotated in possessing both gill and lung tissues at the same time these public databases via Trans-Decoder annotation. were used for RNA-Seq analysis. Their mix was used to Functional annotation of all unigenes against protein da- obtain the transcriptome sequence. A total of tabases showed that a total of 18,435 (61.36%) unigenes 41,878,756 bp of clean data were obtained after quality exhibit significant similarity to known proteins in the filtering (Table 1). As there are no available assembled Pfam database, and a total of 13,341 (44.40%) in the and annotated genomic sequences of the Chinese giant Swissprot database. Four species with most BLASTX salamander that could be used for transcriptome assem- hits against the Nr protein database accounted for bly, Trinity de novo assembler program was used to 42.86% of the identified unigenes: Xenopus tropicalis (14.16%), Chrysemys picta (13.89%), Chelonia mydas (7.76%) and Latimeria chalumnae (7.05%). These re- Table 1 Assembly quality statistics of the Chinese giant salamander sults are comparable to those reported in the larvae transcriptome of tiger salamander and the amphibian species usu- Parameter Value ally used for the respiratory transition studies - Xen- Raw reads 41,878,756 opus tropicalis [49]. The remaining 57.14% were High quality bases 12,482,332,994 distributed as follows: Pelodiscus annectens (5.38%), High quality reads % (≥Q30) 89.39% Xenopus laevis (3.37%), Anolis carolinensis (3.22%), Oncorhynchus mykiss (3.16%), Alligator mississippiensis Unigene number 86,282 (2.90%), Alligator sinensis (2.82%), and other species Unigene total length 77,275,634 (36.28%) (Additional file 4). In total, 9236 (10.70%) uni- Average length (bp) 895.62 genes were allocated into three main GO categories Unigene N50 length (bp) 1732 (Additional file 5): within the cellular component category, Su et al. BMC Genomics (2018) 19:406 Page 6 of 15 cell (5172) and cell part (5171) terms were most might have effect on the miRNA target recognition [50]. prominent; within the molecular function category, To identify the conserved miRNAs in A. davidianus,small GO terms were predominantly associated with bind- RNA sequences were compared with the currently avail- ing activity (5022) and catalytic activity (3622); the able mature miRNAs in the miRBase. The miRNAs that biological process category was dominated by cellular perfectly matched with conserved miRNA sequences or process (5471) and single-organism process (5039). had stem-loop precursors were identified as conserved These results are almost identical to Huang et al.’sre- miRNAs, which indicates their conserved function across port [48] and similar to be found by Eo et al. [49]. In various species. Finally, a total of 65 small RNAs in A. the the latter study [49] these terms were associated davidianus were identified as orthologs of known miR- with transcripts differentially expressed between gills NAs in other species (Additional file 1: Table S10). For and lungs. conserved miRNA, most of them ranged from 21 to 23 nt The unigenes were also queried (by BLASTX) in in length; sorted by abundance: 22 nt > 21 nt > 23 nt. KEGG database to further explore their biological func- These results are fully consistent with Yang et al.’s results tions and interactions. A total of 13,246 unigenes were in Chinese giant salamander [51]. Abundance of known assigned to 288 pathways, among which ribosome (345), miRNA family members ranged from 1 to 13. The most MAPK signaling pathway (311) and focal adhesion (303) abundant were let-1 and mir-142 (13 and 12 members re- were the largest three groups. In addition, 189 unigenes spectively). Among the remaining, two families had 9 were associated with oxidative phosphorylation, 61 with members, and four families had > 3 members and < 9 glutathione metabolism, 46 with metabolism of xenobi- members. To uncover novel miRNAs sequences, all non- otics by cytochrome P450 and 45 with drug metabolism- annotated small RNAs were searched against the newly- cytochrome P450. assembled transcriptome. We identified 46 non-conserved miRNAs belonging to 32 miRNA families, five of which Small RNA sequencing data and annotation (families) contained more than one member. Length dis- To assess the involvement of miRNAs in the respiratory tribution of non-conserved miRNAs was identical to that transition, small RNA libraries were constructed for se- of all miRNAs. Thus, length distribution of miRNAs quencing. A total of 183,422,095 clean reads were gener- might provide a window into the role of miRNAs compo- ated from 12 samples after removing low quality reads, nents in morpho-functional adjustments, because length adapter sequences, and sequences with lengths < 15 nt variation is related to enzymatic modifications, such as and > 35 nt (Additional file 1: Table S8). Q30 percentage RNA editing and exonuclease actives [52]. per individual ranged from 98.83% to 99.13%. Reads from all libraries were annotated using GenBank and Differential expression of miRNAs in gills, lungs and skin Rfam databases (Additional file 1: Table S9). By annota- DESeq program was used to determine whether miRNA tion, Rrnas, tRNAs, snRNAs and snoRNAs were sepa- profiles were different between salamanders depending rated from each other. tRNAs and rRNAs were the most predominantly on gills and skin for respiration and those abundant small RNA types. Length of mature miRNAs depending on lungs for respiration (Table 2). Between was mostly between 20 and 22 nt (27.65%), which were skin and gills in larvae, we found 77 differentially typical dicer-processed miRNA products (Additional file 6). expressed miRNAs, 46 of which were up-regulated and The most abundant small RNAs in all of the libraries were 31 were down-regulated. Among them, four miRNAs 21 nt in length, and a significant bias towards U was ob- were unique to skin and two to gills in larvae. All of served in the first nucleotides of all putative miRNA se- these miRNAs belong to novel miRNAs. Between skin quences across all libraries (Additional file 7). U also had and lungs in larvae, we identified 69 differentially high frequency (51%) in positions 9 and 12. Moreover, A expressed miRNAs, 37 of which were up-regulated and + U were dominant at the start position. These features 32 down–regulated. Among them, only one miRNA, Table 2 miRNA count among different tissues Type DEG Number Up-regulated Down-regulated skin_5 m vs gills_5m 77 46 31 skin_5 m vs lungs_5m 69 37 32 skin_5 m vs lungs_20m 82 55 27 gills_5 m vs lungs_5m 32 16 16 gills_5 m vs lungs_20m 45 30 15 lungs_5 m vs lungs_20m 3 1 2 5 m and 20 m indicate the age of sampled specimens Su et al. BMC Genomics (2018) 19:406 Page 7 of 15 named unconservative_c58354.graph_c0_367043, was identified by intersection of miRnada and RNA hybrid unique to skin in larvae, and seven miRNAs were unique prediction programs. The target intersection set was to gills in larvae. Between gills in larvae and lungs in verified by PITA program [44]. Among the 15,739 genes adults, we identified 82 differentially expressed miRNAs, predicted as possible targets, 8358 (53.10%) were suc- 55 up-regulated and 27 down–regulated. Among these, cessfully annotated against the NR, Swiss-Prot, GO, only one miRNA, named unconservative_c63346.graph_ KEGG and Pfam databases (Additional file 1: Table S11). c0_552043, was unique to skin in larvae and three miR- A total of 352, 291, 1016, 117 and 862 differentially NAs were unique to gills in larvae. Between gills and expressed miRNA target genes were found in pairwise lungs in larvae, we identified 32 differentially expressed comparisons of skin5 m vs gill5m, skin5 m vs lungs5m, miRNAs, half of which were up-regulated and half skin5 m vs lungs20m, gill5 m vs lungs5 m and gill5 m vs down–regulated. Here, six miRNAs were unique to lungs20m, respectively, but no differentially expressed lungs in larvae. Between gills in larvae and lungs in genes were found in the lungs5 m vs lungs20m compari- adults, we identified 45 differentially expressed miRNAs, son. Expression levels of ten target genes related to the 30 of which were up-regulated and 15 down–regulated. validated 12 miRNAs were also determined (Fig. 4). Among them, only one miRNA, named unconservative_ Among them, four genes had higher expression in skin, c60271.graph_c0_422386, was unique to gills in larvae while two genes had higher expression in gills (both and four miRNAs were unique to lungs in adults. Be- were compared to the lung tissue). One gene, semialde- tween the adult and larva lungs, three miRNAs were dif- hyde dehydrogenase, had higher expression in both gills ferentially expressed including one up-regulated and two and skin, while only centrosomal protein of 57 kDa iso- down-regulated. Results of all pairwise comparisons in- form X1 exhibited higher expression level in lungs. dicate that aca-miR-203-5p, dre-miR-203a-5p, ipu-miR- The top-scoring BLASTX hits against the Nr protein 203a and rno-miR-203a-5p were expressed only in the database indicated that the top four species were Chrys- lungs of larvae. Distribution of significant changes de- emys picta (17.12%), Xenopus (Silurana) (13.82%), Chelo- tected is shown as a volcano plot in the Additional file 8. nia mydas (8.22%) and Latimeria chalumnae (7.69%), According to the hierarchical clustering analysis, the 12 which accounted for nearly half of the miRNA targets studied samples were classified into three major clusters, identified against the Nr database. The remaining which corresponds to the three studied tissues (Fig. 2). miRNA targets were distributed among Pelodiscus Therefore, miRNA expression patterns indicate that sinensis (4.81%), Xenopus laevis (4.55%), Anolis caroli- miRNAs are involved in the regulation of adaptation of nensis (3.61%), Alligator mississippiensis (3.51%), A. respiratory tissues during the transition from aquatic lar- sinensis (3.20%), Python bivittatus (1.68%), and other vae to terrestrial adults. species (31.83%). Among the large number of differentially expressed miRNAs, we selected ten for further validation. These Construction of the miRNA-gene-network were chosen based on two criteria, their strong differen- The ultimate goal of our study was to test the hypothesis tial expression among the above four tissues (> 2-fold in- that miRNAs and their target genes are associated with crease) and their strong correlation with ventilation the transition from aquatic to aerial respiration. Thus, features, either strong positive or strong negative: a) no we identified five candidate miRNAs expressed in the expression in some tissues, while expression was ob- specific expression patterns in gill, skin and lung. In served in other tissue(s) in differential expression dataset order to identify putative functional modules, we (using DESeq); b) FDRs were lower than 0.05 for all constructed a miRNA-gene network (Fig. 5). Among pairwise comparisons. Based on these criteria, aca-miR- thefivemiRNAsinthe network, threeweredown- 142-5p, bta-miR-142-5p, gga-miR-34c-3p, aca-let-7a-5p, regulated during the respiratory transition: unconserva- aca-miR-203-3p, aca-miR-203-5p, aja-miR-142, dre-miR- tive_c62709.graph_c3_520410, unconservative_c41959. 203a-5p, efu-miR-223 and ipu-miR-142 were selected for graph_c0_140263 and unconservative_c60463.graph_ further validation experiments (Fig. 3). These results c1_429066; whereas one exhibited increased expression largely corroborated the results of transcriptomic during the respiratory transition: unconservative_ analysis. c60430.graph_c0_428023. The miRNAs exhibiting the widest and narrowest degrees in the miRNA-Gene Identification of target genes of differentially expressed network were unconservative_c60430.graph_c0_428023 miRNAs and their functional annotation and unconservative_c57205.graph_c0_339408, which The ultimate goal of our study was to test the hypothesis had 1361 and 3 target genes respectively. Regarding the that miRNAs and their target genes are associated with total differentially expressed genes, most of them were the transition from aquatic to aerial respiration. Poten- targeted by two and more miRNAs in the network, tial targets of differentially expressed miRNAs were whereas 18 genes were regulated by two miRNAs. Su et al. BMC Genomics (2018) 19:406 Page 8 of 15 Fig. 2 (See legend on next page.) Su et al. BMC Genomics (2018) 19:406 Page 9 of 15 (See figure on previous page.) Fig. 2 Hierarchical cluster analysis of differentially expressed miRNAs in the respiration tissue transfer process. Each row represents a different miRNA, and each column represents a different sample. Clustering based on parameter similarity is shown on the left (miRNAs) and above (samples) the figure. Increasing green intensity denotes decreasing gene expression, and increasing red intensity denotes increasing gene expression, with color scale shown on the right. The heatmap was constructed using expression profiling of each sample, with the data calculated from log10 (TPM + 1) GO categorization of these network targets indicated a 140263 and unconservative_c60463.graph_c1_429066 significant enrichment in the cellular process, single- were highly expressed in lungs of larvae. Unconserva- organism process, metabolic process, biological regulation, tive_c60430.graph_c0_428023 subsequently up-regulated cell, cell part, organelle, binding and catalytic activity the expressions of secretory phospholipase A2 receptor (Additional file 9). These targets were also aligned to the isoform X2 gene, which is functionally associated with KEGG database: a total of 203 differentially expressional oxidative stress-induced premature senescence. Two miRNA target genes were involved in 185 pathways, with down-regulated genes were pulmonary surfactant- ‘focal adhesion’ and ‘MAPK signaling pathway’ being the associated protein B precursor and pulmonary most abundant terms (Additional file 10). surfactant-associated protein B-like isoform X1, which are necessary for vertebrate lung development [56]. Discussion These results are consistent with Eo et al.’s report [49], Our understanding of the mechanisms involved in the and illustrate that miRNAs contribute to the observed transition from gill respiration to lung respiration re- morpho-functional adjustments. In amphibians, at least mains incomplete. In previous reports, candidate pro- 30% of oxygen exchange and most of the carbon dioxide teins [14], nitric oxide synthase (NOS)/nitric oxide [9], elimination occurs through the skin [57]. The skin of and conserved molecular networks [53] were reported to lungless salamanders even serves as their only respira- be involved in this complex process. Herein, we focused tory organ [58]. Among the five genes associated with on the expression profiles of miRNAs. In order to asso- respiration identified in Chinese giant salamander [27], ciate the miRNAs with their target genes, a reference NADH-ubiquinone oxidoreductase chain 5, ATPase sub- transcriptome of the Chinese giant salamander (A. davi- unit 6 and cytochrome b were differentially expressed be- dianus) was generated. As a result, along with already tween skin and other tissues in the present study. Huang assembled transcriptomes of the larval-stage tiger sala- et al. [59] found that one of the most highly expressed mander [49] and Chinese salamander [54], scientists genes related to the respiration in the skin of seven dif- aiming to study the expression of genes in the transition ferent anurans was NADH-ubiquinone oxidoreductase from aquatic to aerial respiration now have transcrip- (GO term: “catalytic activity”), which is consistent with tomes of three salamander species at disposal. our result. Therefore, all of these results illustrate that In the present study, miRNA expression profiles were gills, lungs and skin play an important role in the detected in three different tissues. Characteristics of the morpho-functional adjustments associated with respir- miRNAs detected, in particular 21 nt as the most com- ation transition via oxidative stress, lung development mon length and a significant bias towards U in the first and catalytic activity mechanisms. nucleotide, are considered to be typical for miRNAs Secondly, biological functions of miRNAs and their [55]. Based on these identified miRNAs, two steps were target genes involved in respiration regulation were clas- taken to better understand the association between miR- sified. Five known miRNAs were related to hypoxia- NAs and their target genes during the transition from induced respiration regulation. The let-7 miRNA family aquatic to aerial respiration. members are induced by hypoxia [60], and increased Firstly, we focused on the specific expression of miR- level of let-7a-5p targets major pulmonary arterial NAs in different tissues and different functions of these hypertension pathways [61]. Higher expression level of tissues during the respiration transition process. Expres- let-7a-5p, and its target gene sal-like protein 1 (SALL1), sion profiles of gill and skin tissues could be distin- in gills in comparison to lungs illustrated that gill respir- guished from the expression profile of lung tissue using ation might be a reflection hypoxia induced by the 115 miRNAs. Among these, we validated and confirmed change of environment: from oxygen-abundant terres- the differential expression of 12 miRNAs. Among them, trial to oxygen-scarce aquatic environment. Some stud- three new miRNAs differentially expressed between ies have reported that sal-like protein 1 plays essential lungs of larvae and adults show specific expression pat- roles in maintaining self-renewal and pluripotency of terns: unconservative_c60430.graph_c0_428023 was embryonic stem cells (ESCs) [62] and promotes the highly expressed in lungs of adults compared to lungs of intrahepatic cholangiocarcinoma cell migration [63]. larvae, while both unconservative_c41959.graph_c0_ These results indicate that cell differentiation and Su et al. BMC Genomics (2018) 19:406 Page 10 of 15 Fig. 3 Validation of differentially expressed miRNAs in the respiration tissue transfer process. Different letters above bars indicate that two means are significantly different (P < 0.05) migration may be one possible way of promoting adjust- induced apoptosis and fibrosis of cardiomyocytes. For ments in the respiration process by sal-like protein 1. miR-142, Lu et al. [65] found that hypoxia-induced cell Other three miRNAs have been associated with the in- proliferation and invasion could be inhibited by the hibition of hypoxia-induced adaptive changes: Wang overexpression of miR-142 or hypoxia inducible factor1 et al. [64] reported that miR-142-3-P inhibits hypoxia- (HIF-1α). Down-regulation of pulmonary miR-223 and Su et al. BMC Genomics (2018) 19:406 Page 11 of 15 Fig. 4 Validation of differentially expressed mRNAs in the respiration tissue transfer process. Different letters above bars indicate that two means are significantly different (P < 0.05). Note: c82457.graph_c0: sal-like protein 1; c54346.graph_c0: Semialdehyde dehydrogenase, NAD binding domain; c63687.graph_c2: sorting nexin-33; c53027.graph_c0: Ribosomal protein L11; c41045.graph_c0: Carbon-nitrogen hydrolase; c40974.graph_c0: Iron/zinc purple acid phosphatase-like protein C; c59945.graph_c0: glypican-5-like isoform X3; c60274.graph_c0: reverse transcriptase; c71119.graph_c0: centrosomal protein of 57 kDa isoform X1; c76183.graph_c0: trimethyllysine dioxygenase, mitochondrial isoform X3 upregulation of insulin-like growth factor 1 receptor ex- metabolism, as well as ribosomal structure and biogen- pression can be initiated by pulmonary hypoxia in the esis, this illustrates that metabolic ways were coordi- right heart failure [66]. In our study, the three miRNAs nated with adaptive changes in the skin respiratory mentioned above exhibited higher expression in lungs patterns. For two miRNAs, unconservative_c57205. than in skin and gills, which indicates that pulmonary graph_c0_339408 and unconservative_c62709.graph_c3_ breathing is susceptible to adaptive air respiration by 520410, we did not observe a significant difference in inhibiting the hypoxia-induced gene expression changes. expression among the studied tissues, but one of their However, miR-142 and miR-223 targets, Ribosomal pro- targeted genes exhibited a significant difference in tein L11 and Iron/zinc purple acid phosphatase-like pro- expression: trimethyllysine dioxygenase gene (validated tein C respectively, had higher expression in skin. As by qPCR) had higher expression in skin compared to they are involved in carbohydrate transport and other tissues. As its mitochondrial isoform X3 is known Su et al. BMC Genomics (2018) 19:406 Page 12 of 15 Fig. 5 miRNA-Gene-Network for the respiratory transition in Chinese giant salamander. The miRNA-Gene-Network was built using gene expression data and predicted interactions in the TargetScan miRNA database. The squares represent miRNAs and the circles represent their target genes, where red indicates up-regulation, and blue - down-regulation, purple – no significant regulation. Group1 is skin5 m and Group2 is gill5m, Group3 is lungs5 m and Group4 is lungs20m. There are four pairwise group comparisons (skin5 m vs gill5m, gill5 m vs lungs5m, gill5 m vs lungs20m and lungs5 m vs lungs20m) for differential miRNA expression regulation pattern analysis. So different proportions of different colors in squares and cycles correspond to the differential expression levels among the four groups (pairwise comparisons); for example, if half of the square is filled with red, it means that this miRNA was up-regulated in two of the four groups to be down-regulated in response to hypoxia (10% O2) components, qualitative and quantitative (immuno- [67], this result illustrates that these two miRNAs facili- fluorescence microscopy and Western blotting) pat- tate the respiratory system transition from water to land terns of hypoxia inducible factor1 (HIF1) expression by regulating the oxygen-related genes in lungs and skin. correlated inversely to that of nitric oxide synthase in Eo et al. [49] used larval tiger salamander and pyrose- both gills and lungs [9]. Synaptotagmin-like 4, associ- quencing to identify genes differentially expressed be- ated with the GO terms “phospholipid binding” and tween gill and lung tissues. They also found that “neurexin family protein binding”, is involved in intra- ‘pulmonary-associated protein isoform cra_a’ was exclu- cellular trafficking, secretion, and vesicular transport. sively expressed in the lungs. Among these differentially This gene indicates that lipid metabolism and nervous expressed miRNA-targeted genes, Carbon-nitrogen system may be related with the respiration transition. hydrolase (annotation in pFam) is associated with the Zabelinskii et al. [68] studied the lipid and fatty acid “nitrogen compound metabolic process” GO term. This composition of gills and lungs in 18 fish species and gene was involved in signal molecule changes in the gills seven mammalian species, and found different molar and lungs of African lungfish during the maintenance ratios between phospholipids and cholesterol: 2:1 in and arousal phases of aestivation [9]. They also reported fish and 3:1 in mammals. In vitro, cranial nerve (CN) localization and expression of nitric oxide synthase roots and spinal nerve (SN) roots had different fre- (NOS), Protein Kinase B, phospho-Akt (Akt), heat shock quency, and their amplitude bursts patterns corre- protein 90 (Hsp-90), hypoxia inducible factor α (HIF-1α) sponded to gill and lung ventilatory burst patterns, in a tissue-specific manner in parallel with organ re- respectively [69]. So, lipid metabolism and neural ac- adjustment in the gills and lungs of P. annectens during tivity also play a key role in the process of morpho- aestivation and arousal [9]. Among these molecular functional adjustments. Su et al. BMC Genomics (2018) 19:406 Page 13 of 15 Conclusions Additional file 9: Figure S8f. Go classification of differentially expressed This study provides the first large-scale report of miRNA targets for the miRNA-gene-network. Differentially expressed miRNA targets were allocated into three main GO categories: biological miRNA expression profiles and their predicted target process, cellular component and molecular function. Left vertical axis genes during the respiration transition from gills to represents the percentage of genes, while right vertical axis indicates lungs in Chinese giant salamander. A high-coverage the number of genes. (PNG 1185 kb) reference transcriptome was also generated and as- Additional file 10: Figure S9f. KEGG pathway analysis of differentially expressed miRNA targets for the miRNA-gene-network. The horizontal sembled. Two miRNAs with target genes related to axis represents the number of genes, and the percentage of this number oxygen sensing were differentially expressed between in the total number of genes in a given pathway. On the vertical axis are gill and lung tissues. Three miRNAs were differen- the names of KEGG pathways. (PNG 84 kb) tially expressed between lungs of larvae and lungs of adults. Differential expressions of several miRNAs were validated and confirmed; e.g., let-7a-5p was tar- Abbreviations geting major pulmonary arterial hypertension path- Akt: Phospho-Akt; CN: Cranial nerve; DEGs: Differentially expressed genes; EggNOG: Orthologous groups of genes; ESCs: Embryonic stem cells; FC: Fold ways. Their target gene annotation indicates that the change; FDR: False discovery rate; FPKM: Fragments per kb per million reads; observed miRNA regulation pattern is a response to gill5m: Larvae (5 months-old) gill; GO: Gene Ontology database; HIF-1α:Hypoxia the change in the levels of oxygen in the environ- inducible factor1; Hsp-90: Heat shock protein 90; KEGG: Kyoto Encyclopedia of Genes and Genomes database; lungs20m: Adults (20 months-old) lungs; ment. These results suggest that expression profiles of lungs5m: Larvae (5 months-old) lungs; miRNAs: MicroRNAs; MISA: MIcroSAtellite miRNAs in three respiratory tissues are correlated identification tool; NO: Nitric oxide; NOS: Nitric oxide synthase; PITA: Target with oxygen responses during the transition from accessibility; qPCR: Real time quantitative polymerase chain reaction; RNase H: M-MuLV Reverse Transcriptase; rRNA: Ribosomal RNAs; SALL1: Sal-like protein water to land. 1; SFTPA, SFTPC, and SFTPD: Pulmonary surfactants A, C, D; skin5m: Larvae (5 months-old) skin; SN: Spinal nerve; snoRNA: Small nucleolar RNAs; snRNA: Small nuclear RNAs; SSR: Microsatellite; TPM: Transcripts per million; tRNA: Transfer RNAs; U6: mamm U6 small nuclear 1 Additional files Acknowledgments We thank the students and staff of the Aquatic Genetic Laboratory, Freshwater Additional file 1: Supplementary Tables S1 to S11. (DOC 122 kb) Fisheries Research Center of Chinese Academy of Fishery Sciences for their kind Additional file 2: Figure S1. Size distribution of the assembled assistance in the study. The authors are grateful to Dr. Du for valuable technical unigenes. Horizontal axis gives different size intervals of the assembled support. unigenes, vertical axis gives the number of unigenes located in the specific size interval. (PNG 44 kb) Funding Additional file 3: Figure S2. Frequency distribution of the putative cSSRs This study was supported by grants from Central Public-interest Scientific observed. C is the number of SSRs present in compound formation; C* is Institution Basal Research Fund, Chinese Academy of Fishery Sciences the number of sequences containing more than one SSR; p1 – p5 represent (CAFS) (2016RC-LX03), New Aquaculture Project projects of Jiangsu Province the numbers of mono-, di-, tri-, tetra- and penta-nucleotides respectively. (Y2013–51) and the 2016 annual Jiangsu University “Qing Lan Project” young (PNG 21 kb) academic leaders training project (Jiangsu teachers [2016]15), the key Additional file 4: Figure S3. Distribution of the top BLASTX hits for technology programs of Zhenjiang (NY2016004) and the Key Projects in unigenes in the Nr databases. (PNG 64 kb) the National Science & Technology Pillar Program during the twelfth Additional file 5: Figure S4. GO classification of Chinese giant Five-Year Plan Period (2012BAD26B02), the National Nonprofit Institute salamander unigenes. 9236 (10.70%) unigenes were allocated into three Research Grant of CATAS-TCGRI (2013JBFM14). The funders had no role main GO categories: biological process, cellular component and in the design of the study, collection, analysis and interpretation of data, molecular function. Some of these unigenes were annotated with and in writing the manuscript. multiple GO terms. (PNG 109 kb) Additional file 6: Figure S5. Length distribution of novel miRNAs. Availability of data and materials Horizontal axis gives the length of novel miRNAs, vertical axis gives the The datasets generated during and/or analysed during the current study number of novel miRNAs. (PNG 110 kb) (assembled transcriptome and miRNA sequences) are available in NCBI’s Additional file 7: Figure S6. miRNA nucleotide bias. A: Nucleotide Transcriptome Shotgun Assembly (TSA) database (ftp://ftp.ncbi.nlm.nih.gov/ distribution at the first position of miRNA. Each row represents a genbank/tsa) under the accession numbers: PRJNA401072 and PRJNA398145 miRNA of a different length, and each column represents the (SRR5944935, SRR5944934, SRR5944939, SRR5944938, SRR5944937, distribution of bases in the first nucleotide position. B: Overall SRR5944936, SRR5944931, SRR5944930, SRR5944933, SRR5944932, nucleotide distribution. Each row represents a miRNA of a different SRR5944929, SRR5944928). length, and each column represents the overall distribution of the bases. (PNG 31 kb) Authors’ contributions Additional file 8: Figure S7. Volcano plots of differentially expressed SS, CJ, XJ and XP designed the experiments, supervised the study and miRNA for different respiration-related tissues. Red dots represent significantly modified the manuscript. SS, WY, WH, ZY, HW and HC performed the up-regulated genes and green dots significantly down-regulated genes in research, analyzed the data and wrote the manuscript. YX participated in the pairwise tissue comparisons (blue dots are genes with unchanged expression). design of the study, lab work, data analysis, data visualisation, molecular lab The x-axis represents the log2-transformed gene expression. The y-axis is the p work; SS carried out sequence alignments and contributed custom–made value (−log2) adjusted by Benjamini-Hochberg correction. All pairwise tissue software and codes; CJ and HW participated in sample collection; XJ comparisons revealed differentially expressed miRNAs. A: skin5 m vs gill5m; B: conceived and coordinated the study; XP co-wrote the manuscript, and skin5 m vs lungs5m; C: skin5 m vs lungs20m; D: gill5 m vs lungs5m; E: gill5 m reviewed the manuscript. All authors contributed comments to the final vs lungs20m; F: lungs5 m vs lungs20m. (PNG 742 kb) version of the manuscript, read and approved it. Su et al. BMC Genomics (2018) 19:406 Page 14 of 15 Ethics approval and consent to participate 15. Vila-Casadesus M, Gironella M, Lozano JJ. MiRComb: an R package to Study design and all experimental procedures involving animals were analyse miRNA-mRNA interactions. Examples across five digestive cancers. reviewed, approved and supervised by the Institutional Animal Care and Use PLoS One. 2016;11:e0151127. Committee of the Freshwater Fisheries Research Center of Chinese Academy 16. Xu Y, Chu L, Jin Q, Wang Y, Chen X, Zhao H, Xue Z. Transcriptome-wide of Fishery Sciences. Various measures were taken to minimize the suffering identification of miRNAs and their targets from Typha angustifolia by RNA- of the animals. Seq and their response to cadmium stress. PLoS One. 2015;10:e0125462. 17. Liu X, Zhu L, Liao S, Xu Z, Zhou Y. The porcine microRNA transcriptome Competing interests response to transmissible gastroenteritis virus infection. PLoS One. 2015;10: The authors declare that they have no competing interests. e0120377. 18. Coolen M, Bally-Cuif L. MicroRNAs in brain development and physiology. Curr Opin Neurobiol. 2009;19:461–70. Publisher’sNote 19. Ivey KN, Srivastava D. MicroRNAs as regulators of differentiation and cell fate Springer Nature remains neutral with regard to jurisdictional claims in decisions. Cell Stem Cell. 2010;7:36–41. published maps and institutional affiliations. 20. Bailey SG, Sanchez-Elsner T, Stephanou A, Cragg MS, Townsend PA. Regulating the genome surveillance system: miRNAs and the p53 super Author details family. Apoptosis. 2010;15:541–52. Key Laboratory of Genetic Breeding and Aquaculture Biology of Freshwater 21. Dykxhoorn DM. MicroRNAs and metastasis: little RNAs go a long way. Fishes, Ministry of Agriculture; Freshwater Fisheries Research Center, Chinese Cancer Res. 2010;70:6401–6. Academy of Fishery Sciences, Wuxi 214081, People’s Republic of China. 22. Skalsky RL, Cullen BR. Viruses, microRNAs, and host interactions. Annu Rev Department of Animal Husbandry & Veterinary Medicine, Jiangsu Microbiol. 2010;64:123–41. Polytechnic College of Agriculture and Forestry, Zhenjiang 212400, People’s 23. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. Republic of China. Wuxi Fisheries College, Nanjing Agricultural University, 2004;116:281–97. Wuxi 214081, People’s Republic of China. 24. Nana-Sinkam SP, Karsies T, Riscili B, Ezzie M, Piper M. Lung microRNA: from development to disease. Expert Rev Respir Med. 2009;3:373–85. Received: 4 January 2018 Accepted: 12 April 2018 25. Rupani H, Sanchez-Elsner T, Howarth P. MicroRNAs and respiratory diseases. Eur Respir J. 2013;41:695–705. 26. Nowak JS, Michlewski G. miRNAs in development and pathogenesis of the References nervous system. Biochem Soc Trans. 2013;41:815–20. 1. Wang X-M, Zhang K-J, Wang Z-H, Ding Y-Z, Wu W, Huang S. The decline of 27. Geng X, Wei H, Shang H, Zhou M, Chen B, Zhang F, Zang X, Li P, Sun J, Che the Chinese giant salamander Andrias davidianus and implications for its J, et al. Proteomic analysis of the skin of Chinese giant salamander (Andrias conservation. Oryx. 2004;38:197–202. davidianus). J Proteome. 2015;119:196–208. 2. Zhang L, Jiang W, Wang Q-J, Zhao H, Zhang H-X, Marcec RM, Willard ST, 28. Chen L, Fan J, Hu L, Hu Z, Xie Y, Zhang Y, Lou Y, Nevo E, Fu J. A Kouba AJ. Reintroduction and post-release survival of a living fossil: the transcriptomic analysis of bermudagrass (Cynodon dactylon) provides novel Chinese Giant salamander. PLoS One. 2016;11:e0156715. insights into the basis of low temperature tolerance. BMC Plant Biol. 2015; 3. Simons RS, Bennett WO, Brainerd EL. Mechanics of lung ventilation in a 15:216. post-metamorphic salamander, Ambystoma Tigrinum. J Exp Biol. 2000; 29. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis 203:1081–92. X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly 4. Ijspeert AJ, Crespi A, Ryczko D, Cabelguen JM. From swimming to from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29: walking with a salamander robot driven by a spinal cord model. Sci. 644–52. 2007;315:1416–20. 30. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ, 5. Piiper J. Respiratory gas exchange at lungs, gills and tissues: mechanisms Gapped BLAST. PSI-BLAST: a new generation of protein database search and adjustments. J Exp Biol. 1982;100:5–22. programs. Nucleic Acids Res. 1997;25:3389–402. 6. Burggren W, Doyle M. Ontogeny of regulation of gill and lung ventilation in 31. Finn RD, Coggill P, Eberhardt RY, Eddy SR. The Pfam protein families the bullfrog, Rana catesbeiana. Respir Physiol. 1986;66:279–91. database: towards a more sustainable future. Nucleic Acids Res. 2016;44: 7. Gdovin MJ, Torgerson CS, Remmers JE. Neurorespiratory pattern of gill and D279–85. lung ventilation in the decerebrate spontaneously breathing tadpole. Respir 32. Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walter MC, Rattei Physiol. 1998;113:135–46. T, Mende DR, Sunagawa S, Kuhn M, et al. eggNOG 4.5: a hierarchical 8. Lenfant C, Johansen K. Gas exchange in gill, skin, and lung breathing. Respir orthology framework with improved functional annotations for eukaryotic, Physiol. 1972;14:211–8. prokaryotic and viral sequences. Nucleic Acids Res. 2016;44:D286–93. 9. Garofalo F, Amelio D, Icardo JM, Chew SF, Tota B, Cerra MC, Ip YK. Signal 33. Bairoch A, Apweiler R. The SWISS-PROT protein sequence database and its molecule changes in the gills and lungs of the African lungfish Protopterus supplement TrEMBL in 2000. Nucleic Acids Res. 2000;28:45–8. annectens, during the maintenance and arousal phases of aestivation. Nitric 34. Finn RD, Clements J, Arndt W, Miller BL, Wheeler TJ, Schreiber F, Oxide. 2015;44:71–80. Bateman A, Eddy SR. HMMER web server: 2015 update. Nucleic Acids 10. Hu XG, Liu H, Jin Y, Sun YQ, Li Y, Zhao W, El-Kassaby YA, Wang XR, Mao JF. Res. 2015;43:W30–8. De novo transcriptome assembly and characterization for the widespread 35. Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto and stress-tolerant conifer Platycladus orientalis. PLoS One. 2016;11: encyclopedia of genes and genomes. Nucleic Acids Res. 1999;27:29–34. e0148985. 36. Chen L, Chu C, Lu J, Kong X, Huang T, Cai YD. Gene ontology and KEGG 11. Xiong H, Li Q, Liu S, Wang F, Xiong Z, Chen J, Chen H, Yang Y, Tan X, Luo pathway enrichment analysis of a drug target-based classification system. Q, et al. Integrated microRNA and mRNA transcriptome sequencing reveals PLoS One. 2015;10:e0126492. the potential roles of miRNAs in stage I endometrioid endometrial 37. Thiel T, Michalek W, Varshney RK, Graner A. Exploiting EST databases for the carcinoma. PLoS One. 2014;9:e110163. development and characterization of gene-derived SSR-markers in barley 12. Vidal EA, Moyano TC, Krouk G, Katari MS, Tanurdzic M, McCombie WR, (Hordeum vulgare L.). Theor Appl Genet. 2003;106:411–22. Coruzzi GM, Gutierrez RA. Integrated RNA-seq and sRNA-seq analysis 38. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data identifies novel nitrate-responsive genes in Arabidopsis thaliana roots. BMC with or without a reference genome. BMC Bioinf. 2011;12:323. Genomics. 2013;14:701. 39. Anders S, Huber W. Differential expression analysis for sequence count data. 13. Wang Y, Brahmakshatriya V, Lupiani B, Reddy SM, Soibam B, Benham AL, Gunaratne P, Liu HC, Trakooljul N, Ing N, et al. Integrated analysis of Genome Biol. 2010;11:R106. microRNA expression and mRNA transcriptome in lungs of avian influenza 40. Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery virus infected broilers. BMC Genomics. 2012;13:278. rate in behavior genetics research. Behav Brain Res. 2001;125:279–84. 14. Biscotti MA, Gerdol M, Canapa A, Forconi M, Olmo E, Pallavicini A, Barucca 41. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 M, Schartl M. The lungfish transcriptome: a glimpse into molecular accurately identifies known and hundreds of novel microRNA genes in evolution events at the transition from water to land. Sci Rep. 2016;6:21571. seven animal clades. Nucleic Acids Res. 2012;40:37–52. Su et al. BMC Genomics (2018) 19:406 Page 15 of 15 42. Bonnet E, Wuyts J, Rouze P, Van de Peer Y. Evidence that microRNA 64. Wang Y, Ouyang M, Wang Q, Jian Z. MicroRNA-142-3p inhibits hypoxia/ precursors, unlike other non-coding RNAs, have lower folding free energies reoxygenationinduced apoptosis and fibrosis of cardiomyocytes by than random sequences. Bioinformatics. 2004;20:2911–7. targeting high mobility group box 1. Int J Mol Med. 2016;38:1377–86. 43. Lin YL, Ma LT, Lee YR, Lin SS, Wang SY, Chang TT, Shaw JF, Li WH, Chu FH. 65. Lu Y, Ji N, Wei W. MiR-142 modulates human pancreatic cancer proliferation MicroRNA-like small RNAs prediction in the development of Antrodia and invasion by targeting hypoxia-inducible factor 1 (HIF-1alpha) in the cinnamomea. PLoS One. 2015;10:e0123245. tumor microenvironments. Biol Open. 2017;6:252–9. 66. Shi L, Kojonazarov B, Elgheznawy A, Popp R, Dahal BK, Bohm M, Pullamsetti 44. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. The role of site SS, Ghofrani HA, Godecke A, Jungmann A, et al. miR-223-IGF-IR signalling in accessibility in microRNA target recognition. Nat Genet. 2007;39: hypoxia- and load-induced right-ventricular failure: a novel therapeutic 1278–84. approach. Cardiovasc Res. 2016;111:184–93. 45. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using 67. Ganfornina MD, Perez-Garcia MT, Gutierrez G, Miguel-Velado E, Lopez-Lopez real-time quantitative PCR and the 2(−Delta Delta C (T)) method. Methods. JR, Marin A, Sanchez D, Gonzalez C. Comparative gene expression profile of 2001;25:402–8. mouse carotid body and adrenal medulla under physiological hypoxia. 46. Wang X, Li M, Wang Z, Han S, Tang X, Ge Y, Zhou L, Zhou C, Yuan Q, Yang J Physiol. 2005;566:491–503. M. Silencing of long noncoding RNA MALAT1 by miR-101 and miR-217 68. Zabelinskii SA, Brovtsyna NB, Chebotareva MA, Gorbunova OB, Krivchenko inhibits proliferation, migration, and invasion of esophageal squamous cell AI. Comparative investigation of lipid and fatty acid composition of fish gills carcinoma cells. J Biol Chem. 2015;290:3925–35. and mammalian lungs. A model of the membrane lipid component areas. 47. Li F, Wang L, Lan Q, Yang H, Li Y, Liu X, Yang Z. RNA-Seq analysis and gene Comp Biochem Physiol B Biochem Mol Biol. 1995;111:127–40. discovery of Andrias davidianus using Illumina short read sequencing. PLoS 69. Torgerson CS, Gdovin MJ, Remmers JE. Fictive gill and lung ventilation in One. 2015;10:e0123730. the pre- and postmetamorphic tadpole brain stem. J Neurophysiol. 1998;80: 48. Huang Y, Xiong JL, Gao XC, Sun XH. Transcriptome analysis of the Chinese 2015–22. giant salamander (Andrias davidianus) using RNA-sequencing. Genom Data. 2017;14:126–31. 49. Eo SH, Doyle JM, Hale MC, Marra NJ, Ruhl JD, DeWoody JA. Comparative transcriptomics and gene expression in larval tiger salamander (Ambystoma tigrinum) gill and lung tissues as revealed by pyrosequencing. Gene. 2012; 492:329–38. 50. Ji Z, Wang G, Xie Z, Zhang C, Wang J. Identification and characterization of microRNA in the dairy goat (Capra hircus) mammary gland by Solexa deep- sequencing technology. Mol Biol Rep. 2012;39:9361–71. 51. Huang Y, Yang YB, Gao XC, Ren HT, Sun XH. Identification and characterization of the Chinese giant salamander (Andrias davidianus) miRNAs by deep sequencing and predication of their targets. 3 Biotech. 2017;7:235. 52. Li M, Xia Y, Gu Y, Zhang K, Lang Q, Chen L, Guan J, Luo Z, Chen H, Li Y, et al. MicroRNAome of porcine pre- and postnatal development. PLoS One. 2010;5:e11541. 53. Cushing L, Jiang Z, Kuang P, Lu J. The roles of microRNAs and protein components of the microRNA pathway in lung development and diseases. Am J Respir Cell Mol Biol. 2015;52:397–408. 54. Che R, Sun Y, Wang R, Xu T. Transcriptomic analysis of endangered Chinese salamander: identification of immune, sex and reproduction-related genes and genetic markers. PLoS One. 2014;9:e87940. 55. Pontes O, Costa-Nunes P, Vithayathil P, Pikaard CS. RNA polymerase V functions in Arabidopsis interphase heterochromatin organization independently of the 24-nt siRNA-directed DNA methylation pathway. Mol Plant. 2009;2:700–10. 56. Hyatt BA, Resnik ER, Johnson NS, Lohr JL, Cornfield DN. Lung specific developmental expression of the Xenopus laevis surfactant protein C and B genes. Gene Expr Patterns. 2007;7:8–14. 57. Jorgensen CB. Amphibian respiration and olfaction and their relationships: from Robert Townson (1794) to the present. Biol Rev Camb Philos Soc. 2000;75:297–345. 58. Nussbaum RA, Wilkinson M, New Genus A. Of lungless tetrapod: a radically divergent caecilian (Amphibia: Gymnophiona). Proc R Soc B Biol Sci. 1995; 261:331–5. 59. Huang L, Li J, Anboukaria H, Luo Z, Zhao M, Wu H. Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Rep. 2016;6:24069. 60. Chen Z, Lai TC, Jan YH, Lin FM, Wang WC, Xiao H, Wang YT, Sun W, Cui X, Li YS, et al. Hypoxia-responsive miRNAs target argonaute 1 to promote angiogenesis. J Clin Invest. 2013;123:1057–67. 61. Wu D, Talbot CC Jr, Liu Q, Jing ZC, Damico RL, Tuder R, Barnes KC, Hassoun PM, Gao L. Identifying microRNAs targeting Wnt/beta-catenin pathway in end-stage idiopathic pulmonary arterial hypertension. J Mol Med (Berl). 2016;94:875–85. 62. Zhang X, Yuan X, Zhu W, Qian H, Xu W. SALL4: an emerging cancer biomarker and target. Cancer Lett. 2015;357:55–62. 63. Zhu L, Huang F, Deng G, Nie W, Huang W, Xu H, Zheng S, Yi Z, Wan T. Knockdown of Sall4 inhibits intrahepatic cholangiocarcinoma cell migration and invasion in ICC-9810 cells. Onco Targets Ther. 2016;9: 5297–305.

Journal

BMC GenomicsSpringer Journals

Published: May 29, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off