Transcriptional landscapes of Axolotl (Ambystoma mexicanum)

Developmental Biology xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Developmental Biology journal homepage: www.elsevier.com/locate/developmentalbiology

Original research article

Transcriptional landscapes of Axolotl (Ambystoma mexicanum) Juan Caballero-Péreza,b,1, Annie Espinal-Centenoa,1, Francisco Falcona, Luis F. García-Ortegaa, Everardo Curiel-Quesadab, Andrés Cruz-Hernándezc, Laszlo Bakod, Xuemei Chene, Octavio Martínezf, Mario Alberto Arteaga-Vázquezg, Luis Herrera-Estrellah, Alfredo ⁎ Cruz-Ramíreza, a Molecular and Developmental Complexity Group, Unidad de Genómica Avanzada-Laboratorio Nacional de Genómica para la Biodiversidad (LANGEBIO), CINVESTAV, Irapuato, Guanajuato, Mexico b Department of Biochemistry, Escuela Nacional de Ciencias Biológicas, Instituto Politécnico Nacional, Mexico City, Mexico c Universidad Autónoma de Querétaro, Santiago de Querétaro, Querétaro, Mexico d Umeå Plant Science Centre, Department of Plant Physiology, Umeå University, Sweden e Department of Botany and Plant Sciences, Institute of Integrative Genome Biology, University of California, Riverside, USA f Computational Biology Group, Unidad de Genómica Avanzada-Laboratorio Nacional de Genómica para la Biodiversidad (LANGEBIO), CINVESTAV, Irapuato, Guanajuato, Mexico g Instituto de Biotecnología y Ecología Aplicada (INBIOTECA), Universidad Veracruzana, Veracruz, Mexico h Metabolic Engineering Group, Unidad de Genómica Avanzada-Laboratorio Nacional de Genómica para la Biodiversidad (LANGEBIO), CINVESTAV, Irapuato, Guanajuato, Mexico

A R T I C L E I N F O

A BS T RAC T

Keywords: Coding and non-coding RNAseq Axolotl Regenerative tissues MicroRNAs

The axolotl (Ambystoma mexicanum) is the vertebrate model system with the highest regeneration capacity. Experimental tools established over the past 100 years have been fundamental to start unraveling the cellular and molecular basis of tissue and limb regeneration. In the absence of a reference genome for the Axolotl, transcriptomic analysis become fundamental to understand the genetic basis of regeneration. Here we present one of the most diverse transcriptomic data sets for Axolotl by profiling coding and noncoding RNAs from diverse tissues. We reconstructed a population of 115,906 putative protein coding mRNAs as full ORFs (including isoforms). We also identified 352 conserved miRNAs and 297 novel putative mature miRNAs. Systematic enrichment analysis of gene expression allowed us to identify tissue-specific protein-coding transcripts. We also found putative novel and conserved microRNAs which potentially target mRNAs which are reported as important disease candidates in heart and liver.

1. Introduction Among caudata amphibians, also known as salamanders, the Ambystoma genus comprises 32 closely related species ranging from southern Canada to Central Mexico. In Mexico, endemic species of this genus are found in the so called Transmexican Volcanic Belt, one of the most biodiverse regions of Mexico, which also displays high endemism (Ochoa and Villela, 2006). Among the 16 mexican Ambystoma species, the best studied is the named Axolotl (Ambystoma mexicanum). The Axolotl was established as a model organism over 100 years ago and constitutes an important vertebrate model for studying regeneration and tissue repair due to its capacity to regenerate several body structures such as tail, limb, central nervous system, and tissues of

⁎

1

the eye and heart (Roy and Gatien, 2008; Kragl et al., 2009; Voss et al., 2009; Khattak et al., 2014). In addition to its regeneration capacity, the Axolotl has also been focus of attention given its ability to retain the gross morphology of larvae once they become reproductive adults, a phenomenon known as neoteny or paedomorphosis (Voss et al., 2009). This means that the Axolotl does not suffer metamorphosis to become a terrestrial salamanders as other members of the genus do (Parra‐Olea et al., 2007). However, induced metamorphosis using thyroid hormone has also been implemented to unravel the molecular basis of metamorphosis (Huggins et al., 2012). The Axolotl is one of amphibians for which the most comprehensive scientific tools have been developed (Roy and Gatien, 2008; Voss et al., 2009) but unfortunately it is an endangered species for which more individuals can be in research

Corresponding author. E-mail address: [email protected] (A. Cruz-Ramírez). These authors contributed equally to this work and should be considered as joint first authors.

http://dx.doi.org/10.1016/j.ydbio.2017.08.022 Received 26 April 2017; Received in revised form 12 August 2017; Accepted 17 August 2017 0012-1606/ © 2017 Published by Elsevier Inc.

Please cite this article as: Caballero-Perez, J., Developmental Biology (2017), http://dx.doi.org/10.1016/j.ydbio.2017.08.022

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

gene regulation (Rand et al., 2004). The canonical pathway for the biosynthesis of miRNAs starts with the transcription of a long RNA molecule known as the primary microRNA (pri-miRNA) by the RNA polymerase II (Lee et al., 2004). Then, the protein complex known as the microprocessor, formed by the proteins Drosha and DiGeorge Critical Region 8 (DGCR8), recognizes the pri-miRNA, and then cleaves it generating a ~70 nt RNA hairpin called the pre-miRNA (Kim and Nam, 2006; Denzler and Stoffel, 2015). Pre-miRNAs are exported into the cytoplasmic space of the cell by the karyopherin Exportin-5 (Gwizdek et al., 2003; Yi et al., 2003), where the RNA-induced silencing complex (RISC) recognizes the molecule. RISC's functions include cleaving the pre-miRNA into a RNA duplex by the associated protein Dicer and loading said duplex into the protein AGO to perform gene regulation (Rand et al., 2004; Gregory et al., 2005). MicroRNAs are distributed across the animal kingdom, and it has been suggested that once a novel microRNA emerges, it is rarely lost in a lineage (Wheeler et al., 2009). Under this assumption, it is expected to find several microRNAs conserved in Ambystoma mexicanum. Since last decade, the number of publications addressing microRNA has increased in an overwhelming manner, going from the first publications (Lagos-Quintana et al., 2001), to the more than 26,000 papers that exist to date (Casey et al., 2015). This may be due to the quantity of protein coding genes regulated by these molecules, that has been inferred exceeds the 60% (Friedman et al., 2009), but also how these molecules seem to be connected to different biological phenomena, such as cancer and development (Alvarez-Garcia and Miska, 2005; Xie et al., 2013). In this work, we present a highly diverse transcriptomic data set for Axolotl by profiling both mRNAs and microRNAs from tissues and limbs with proven or potential ability to regenerate such as heart, liver, legs, gills and tail by using RNA-seq and small-RNA-seq technologies. We also focused on the enriched mRNAs and miRNAs in heart and liver.

Fig. 1. Ambystoma mexicanum. Wild type individuals of Axololt from Xochimilco in diverse developmental stages: (A) Late embryo. (B) Larvae. (C), (D) Adult. Bars in all pictures = 1 cm.

laboratories around the globe than in the wild, since during the last couple of decades its population reached extinction levels due to the diverse threatening factors in its natural habitat at Xochimilco in Mexico City (Fig. 1; CONABIO, 2011). In the past decades, efforts to obtain global genetic and genomic data have been implemented, mainly using organisms that have lived for long periods in captivity, a feature that increases endogamy and deleterious effects. The great majority of expression data sets come from transcriptional expression profiles, obtained through DNA microarray, ESTs and RNA-seq, exclusively in function of tissue and limb regeneration (Habermann et al., 2004; Putta et al., 2004; Monaghan et al., 2007, 2009, 2012; Knapp et al., 2013; Stewart et al., 2013). There have also been reports for data obtained from genomic fragments, in some cases by BACs generation and sequencing (Voss et al., 2001; Smith, 2006; Smith et al., 2009). However, the full sequencing of the Axolotl genome has not been accomplished so far. The size of this genome is approximately 10 times bigger than that of the human (estimated in 32 Gbp), making its full sequencing an expensive and time consuming approach that only few laboratories or consortia can afford. One alternative to study the genome is focusing only on the genes being expressed in a specific time and tissue using massive sequencing of RNA molecules (RNA-seq). Since gene expression and the resulting coding and non-coding RNAs are the determinants of the cellular programs and phenotypes, we decided to generate and characterize transcriptomes as a starting point to unravel the molecular basis of tissue and organ homeostasis and they will be useful to compare with other RNAseq studies that focus on regeneration and metamorphosis developmental processes. In addition, this data can also be helpful in studies of biodiversity, evolution and ultimately the implementation of conservation initiatives. RNA-seq is a technological approach where messenger RNA is purified from samples, fragmented and sequenced by next-generation sequencers (NGS). The fragments obtained represent the transcriptome at base resolution, providing a quantifiable measure of the expression of the genes at the transcript level (Mortazavi et al., 2008). It also can be as specific as the quantification at the isoform level for genes with variable exons (Wang et al., 2009). This new technology has been quickly adapted for non-model organisms or, as in the case of Axolotl, when there is no reference genome available. Computational tools have been developed to reconstruct the transcriptome from the fragments, assemblers as Trinity (Grabherr et al., 2011) are commonly used to perform this task. RNA-seq allows to observe the presence and expression level of coding and non-coding genes. Among non-coding genes, microRNAs are small RNA molecules composed of ~ 22 nucleotides. Like other small RNAs, miRNAs associate with the Argonaute (AGO) family of proteins to carry out

2. Results 2.1. RNA-seq library generation and sequencing High quality RNA samples from heart, liver, gills, lung, front-leg and rear-leg were isolated from an adult A. mexicanum female organism. Since Axolotl from the wild are rarely found and are in the top list of endangered endemic organisms, we decided to sacrifice only one individual for sequencing and another two organisms were used for validation. RNAs were sequenced with Illumina HiSeq-2000 platform in a format of 2×100 bp. The six libraries produced 13549557, 10704899, 9519524, 13549693, 13549457 and 13549455 paired-reads for gills, heart, liver, rear leg, front leg and tail, respectively. Raw data was deposited in NCBI Short Read Archive (SRA) under accession SRP093628. 2.2. Transcriptome reconstruction Each tissue sample was assembled individually using Trinity v2.0.4 (Grabherr et al., 2011). Twelve additional samples were also assembled separately combining all samples from Stewart et al. (2013), (SRA accession numbers: SRR389053 to SRR389064). Combined assembly was not possible because those samples are 84 bp single-end reads but were added because that increases our gene coverage. We found that individual assembly of tissues produces larger and probably most complete forms of transcripts. We combined all assembled transcriptomes using CD-HIT-EST program keeping the best gene model produced by all combined sets. This produced 483,060 sequences with an average size of 719.5 bp (Fig. 2B). We searched for partial and complete ORFs in the sequences reconstructed, we keep only ORFs with a minimal size of 100 aa, detecting a 26.52% of potential protein-coding sequences which are encoded by 115,906 full 2

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 2. Transcriptome landscape reconstructed from RNA-seq data sets. (A) Proportions of coding, non-coding and ribosomal transcripts. (B) Transcript length distribution, the size is in base pairs. (C) Proportions of coding and non-coding transcripts per tissue. (D) Principal Component for first two components of the six samples.

transcripts in liver (2325), heart (1742), gills (847), tail (616), rear (432) and front leg (388) at > 1 cpm in the tissue but < 1 cpm in the rest of the tissues (Tables 1 and 2, Fig. 3). Due to their potential for biomedical research, we focused in liver and heart enriched transcripts, using more astringent criteria for transcript selection: a) whose expression was 4-fold over the maximal expression in the other tissues and b) showing a minimal expression of 3 cpm. Using these criteria, we identified 2137 liver-enriched and 928 heart-enriched transcripts. The top 10 enriched transcripts for these tissues are listed in Tables 1 and 2.

ORFs and 109,186 partial ORFs, therefore our reconstruction can be still improved to obtain the real full length transcripts from Axolotl genes and to determine single gene transcripts and isoforms (Fig. 2A). After quantification and normalization, we detected expression of 245, 429 transcripts in all tissues with a minimal expression of 0.1 cpm. Then we filtered the expression to be > 1 (cpm) to consider a transcript present in the tissue, and found 55664, 64959, 70333, 66542, 65348 and 67219 coding transcripts and 52861, 77073, 88412, 75778, 82338 and 84547 non-coding RNAs in liver, tail, rear leg, front leg, heart and gills, respectively (Fig. 2C and Supp. Fig. 1). The reconstructed transcriptome contains many transcripts that might be fragmented, but in proportion they are in the expected ratio of protein coding sequences as observed when compared to the human genome (Supp. Fig. 2). A Principal Component Analysis (PCA) for transcript expression in tissues showed that the populations of transcripts expressed in tail and legs are quite similar, while those of heart and liver are diverse from the other samples, indicating that there are specific transcripts related to the specialized functions of these two organs (Fig. 2D). We also found that a total of 15,003 transcripts are expressed in the 6 tissues (Fig. 3). There is a second group of transcripts (1683), which expression is shared among the legs, tail and gills (Fig. 3).

2.4. Ortholog mRNAs among Axolotl and Human The enriched genes detected were aligned with human protein sequences, keeping only the best match with an E-value < 1e-6. This gave us a subset of 1227 and 466 axolotl transcripts orthologs of 697 and 253 human proteins for liver and heart tissues, respectively. After Gene Ontology enrichment analysis (see Methods), we found that heart-related diseases in human are significantly high in our gene sets for heart-enriched genes in axolotl. For example, cardiomyopathy (DOID:0050700, DOID:0060036, DOID:12930) has a q-value < 10−8, and other categories are also significant such as heart disease (DOID:114, q < 10−6), hypertrophic cardiomyopathy (DOID:11984, q < 10−5), congestive heart failure (DOID:6000, q < 10−4), restrictive cardiomyopathy (DOID:397, q < 10−3) or atrial heart septal defect (DOID:1882, q < 10−3) (Fig. 4A, Supp. Table 1). In a similar analysis, we observed in liver-related diseases in human for our gene sets for liver-enriched genes in axolotl with categories such as inherited metabolic disorder (DOID:655, q < 10–13), disease of metabolism (DOID:0014667, q < 10–13), acquired metabolic disease

2.3. Tissue-enriched mRNAs In order to identify tissue-enriched transcripts, libraries were aligned to the reconstructed transcriptome and quantified their expression pattern among the sequenced samples. Uniquely mapped reads were counted to produce an expression matrix for all transcripts in the 6 samples. From these expression tables, we identified enriched 3

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 3. Expression patterns of reconstructed transcripts. A transcript is considered expressed in the tissue if its expression is higher than 3 cpm, with this criterion, the largest proportion of sequences are co-expressed in all tissues. Table 1 Top 10 Heart-enriched mRNAs. Enrichment is computed as the ratio of the tissue expression over the maximal expression in other tissues. SeqID

Enrichment

AveExpr (CPM)

Human Gene Symbol

Human Gene Description

Amex_g262475 Amex_g256644 Amex_g281412 Amex_g216739 Amex_g280879 Amex_g263235 Amex_g446849 Amex_g254785 Amex_g254847 Amex_g217600

> 3000 > 3000 > 3000 > 3000 2032.04 960.97 923.30 539.56 514.61 378.14

4.68 4.33 3.92 3.26 809.04 104.35 126.24 135.16 988.30 34.56

KCNG4 JPH3 LDB3 MYH8 TNNI3 NPPA TNNI3 HSPB7 TNNT2 MYH7

Potassium channel, voltage gated modifier subfamily G, member 4 Junctophilin 3 LIM domain binding 3 Myosin, heavy chain 8, skeletal muscle, perinatal Troponin I type 3 (cardiac) Natriuretic peptide A Troponin I type 3 (cardiac) Heat shock protein family B (small) member 7 Troponin T type 2 (cardiac) Myosin, heavy chain 7, cardiac muscle, beta

Table 2 Top 10 Liver-enriched Am-mRNAs. Enrichment is computed as the ratio of the tissue expression over the maximal expression in other tissues. SeqID

Enrichment

AverageExpr (CPM)

Human Gene Symbol

Human Gene Description

Amex_g2308 Amex_g258591 Amex_g2477 Amex_g259376 Amex_g281047 Amex_g281060 Amex_g269346 Amex_g272336 Amex_g1395 Amex_g765

> 3000 > 3000 > 3000 > 3000 > 3000 > 3000 > 3000 > 3000 > 3000 > 3000

428.02 216.18 212.59 135.85 103.96 86.42 77.90 69.47 67.88 62.79

A2M CRIM1 F9 RETN LGALS2 LEAP2 ABCB11 INHBB RDH11 HPN

alpha-2-macroglobulin cysteine rich transmembrane BMP regulator 1 (chordin-like) coagulation factor IX resistin lectin, galactoside-binding, soluble, 2 liver expressed antimicrobial peptide 2 ATP binding cassette subfamily B member 11 inhibin beta B retinol dehydrogenase 11 hepsin

(DOID:0060158, q < 10−9), lipid metabolism disorder (DOID:3146, q < 10−9), fatty liver disease (DOID:9452, q < 10−8), familial hyperlipidemia (DOID:1168, q < 10−8) or hepatobiliary disease (DOID:3118, q10−8). (Fig. 4B, Supp Table 1).

closely resemble the predicted patterns, indicating that our RNA-seq analyses generated reliable tissue-enriched patterns.

2.5. Validation of mRNA expression patterns

RNA samples from four tissues (liver, tail, front-leg and rear-leg) were processed and sequenced with Illumina Hiseq. 2000 with the sRNA protocol (as described in methods). The four libraries produced 26298923, 42955452, 21521750 and 16505286 reads for the front-leg, rear-leg, liver and tail, respectively. Raw data was deposited in NCBI Short Read Archive (SRA) under accession SRP093628.

2.6. Small-RNA library generation and sequencing

In order to validate our transcriptome, we selected a group of 10 genes based on their characteristic expression pattern in the 6 tissues sequenced. The main criteria to select them was that they were specific of, or highly enriched in, a given tissue as predicted by the in-silico analyses (Fig. 5, panels A and C). The expression pattern of each of the selected transcripts was determined by qPCR assays using biological duplicates (Fig. 5, panel B and D). Our experiments demonstrated that although the qPCR-based analyses for some transcripts showed slight differences with the predicted behavior from or transcriptomes along the 6 tissues, the tissue-enriched expression patterns of all of them

2.7. MicroRNA expression landscape Our analysis of small-RNA libraries showed a behavior of sequence length distribution with a good quality and homogeneity across samples. Such sequences have a preferential length of ∼ 22 nt, probably 4

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 4. Statistically significant human-disease ontologies identified by term-enrichment of human-ortholog genes. Tissue-enriched genes from Axolotl (A) heart, (B) liver.

expression in front and rear leg). However, one of the candidates for validation, miR454b, behaved completely different in the qPCRobtained expression pattern when compared to that obtained by the in-silico analyses of the miRNA-Seq data. This might be the result of additional members of this miRNA family. Besides the discrepancies previously described, most the miRNAs showed similar behavior in both qPCR and in silico analyses, suggesting that the RNA-seq analyses are reliable.

belonging to miRNAs (Fig. 6A). Using the dataset of the total of sRNAs sequences we searched for conserved miRNAs by using miRBase v.21 database (Kozomara and Griffiths-Jones, 2014). This analysis revealed that from a population of several dozens of millions of sequences obtained, we identified a total of 312 (front leg), 301 (Rear Leg), 278 (Liver) and 270 (Tail) sequences that mapped the miRBase dataset and passed a threshold of 20 reads (see methods and Suppl. Table 6). We also found a total of 522 putative pre-miRNAs which possibly harbour mature miRNAs. Across all samples 448 corresponding to putative novel miRNAs and 74 conserved miRNAs (Suppl. Table 5), while the number of putative mature miRNAs identified as conserved miRNAs was 352, shared across samples. The lower number in mature miRNAs when compared with the precursors might be expected since miRNAs with the same sequence can be matured from diverse premiRNAs (Volk and Shomron, 2011). Using normalized count per tissue, we identify the tissue-enriched microRNAs using a fold difference of 2X to be considered enriched relative to the maximal expression in the other tissues in cpm. We identified 35 liver-enriched miRNAs and 17 miRNAs enriched in Tail, while the front and rear Leg did not show any evident enriched miRNAs (Table 3).

2.9. Putative novel MicroRNAs in Axolotl In order to identify novel putative miRNAs which were not previously observed in miRBase, we established the following criteria: putative novel miRNAs most have a suitable precursor (Friedländer et al., 2012), such that small-RNAs would align to the either side of the hairpin and not the loop, and that it has a lower folding free energy than a random RNA sequence (Bonnet et al., 2004) (Fig. 8). By using a miRDeep threshold of 5, and by filtering pre-miRNAs containing conserved miRNAs (54 different mature miRNAs distributed in 74 different precursors), we found 297 different mature miRNAs, which mapped to 448 pre-miRNAs.

2.8. Validation of MicroRNA expression patterns 2.10. Target analyses of Putative novel miRNAs results In order to validate the expression patterns of miRNAs, which were predicted by the in-silico analyses, we selected a group of 6 miRNAs based on the expression patterns in the 4 tissues sequenced. Similarly to our analyses of mRNAs, the main criteria to select them was that they were specific of, or highly enriched in, one of the 4 organs/limbs. Under such criteria, we selected 4 miRNAs that, based on the in-silico analyses, showed to be highly expressed in liver: mir-122, miR-233, miR425-5p and miR3600, and two that showed a low expression also in this organ: miR-429 and miR-454b (Fig. 7 left graphics). Our qPCR analyses (Fig. 7 right graphics) showed that the expression patterns of mir-122, miR-233, miR425-5p and miR3600 were highly like those predicted by our in-silico analyses (Fig. 7 left graphics). In the case of miR-429 we observed similarity among the expression pattern obtained by qPCR and that of the predicted by the bioinformatic analyses, since only the predicted pattern for 3 of the 4 organ/tissues was similar among both experiments (a lower expression in liver, and a higher

We also searched for potential targets of the 297 putative novel miRNAs with a bioinformatic screen using MiRanda (Enright et al., 2003). We considered a miRNA targeting a sequence if the miRNA has a strong 5′ seed, at least to bind sites in the sequences and a maximal Energy of −20 kcal/mol. We found that all novel miRNAs have at least one potential target in our transcriptome and an average target number of 213 transcripts (Supp. Table 2). Interestingly, the targets of 176 novel miRNAs seem to be related to development in embryonic states or tissues, and those of 144 miRNAs are apparently related to transcription with DNA-template, a typical function for microRNAs. We also observed that 28 miRNAs target mRNAs that are reported as involved in controlling neuron development (Supp. Table 3). These data are quite important, since they open a new pathway to follow in future functional studies that focus on the role of these miRNAs in diverse developmental processes, including regeneration. 5

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 5. Predicted differential expression patterns of transcripts and validation by qPCR. In total, the expression pattern of ten transcripts was validated. (A), (C) graphics for each transcript showing the expression levels as obtained from the RNA-seq data of six different tissues (heart, liver, gills, tail, front leg and rear leg). (B),(D) graphics of each transcript validated for two biological replicates by qPCR are expressed as relative quantification [2 ^ (- Δ Cт)], considering as the endogenous control GAPDH transcripts, data are the mean of two biological replicates, each one with three technical replicates.

annotate 225,092 transcripts which show partial or complete open reading frames (ORFs) from which 115,906 are complete ORFs which are putative protein-coding transcripts that matched with the human and xenopus gene set. Our work renders more information when compared with a previous RNA-seq approach which focused in

3. Discussion We have obtained and analyzed a multi-tissue transcriptional landscape from Ambystoma mexicanum by deep sequencing of mRNAs and miRNAs expressed in six different tissues. We could 6

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 6. MiRNA expression in 4 different tissues. (A) Sequence length distribution after trimming. (B) Principal Components projection of first two components. (C) total miRNAs identified by tissue.

factor in cancer progression and metastasis (Zeng and Tang, 2014). In humans SULT1A1 (Sulfotransferase 1A1) transcriptional expression is also liver-enriched, while its up-regulation in tissues where it is normally absent causes cellular alterations. For example, induced expression of this gene has been associated with human breast cancer (Yao-Borengasser et al., 2014). We also found liver-enriched transcripts that are not as highly expressed as the previously described. This is the case of the transcript of the human homolog gene SERPIND1 (Fig. 5 A-B), which codes for Heparin cofactor II, a very important coagulation factor which proper function and expression is elemental to avoid thrombosis and that has been pointed as key for the relationship among coagulation disorders in patients with chronic liver disease (reviewed in (Peck-Radosavljevic, 2007)). A similar story can be told for axolotl heart-enriched transcripts and its human homologs. In the top heart-enriched transcripts is AmLBD3, its human homolog is called LBD3 or ZASP, that encodes a LIM domain binding 3 protein which function is to stabilize the basic units of muscles, the sarcomere. Mutations in this gene cause myofibrillar myopathy, an autosomal dominant muscular dystrophy in humans (Selcen and Engel, 2005). The axolotl homolog of the human MYH8 gene is also highly expressed in heart, mutation of this gene is associated with the hereditary syndrome Carney complex, a heart– hand syndrome that involves cardiac myxomas. Moreover, mutations in this myosin are associated with cardiac tumorigenesis (Veugelers et al., 2004). Although not so highly expressed, another heart-specific gene is the axolotl homolog of the human GATA4 transcription factor (Fig. 5 A-B). Mutations in GATA4 cause defects in cardiomyocyte proliferation during embryogenesis, a condition that leads to the development of congenital heart disorders in humans (Misra et al., 2012). These findings are quite relevant, since they point to the use of the axolotl as an experimental model not only for regeneration studies, but also as a model to study factors associated to heart and liver

profiling gene expression during axolotl limb regeneration (Stewart et al., 2013), where a total of 10,242 protein-coding transcripts were identified when mapped using the human gene set. Recently, in Bryant et al., (2017), the authors reported a new transcriptome for Axolotl using other tissues and organs. We compared our reconstructed transcriptome and found that our assembled contigs are more robust (average size: 719.5 bp vs 511.6 bp, median size: 354 bp vs 288 bp, total contigs: 480,881 vs 1,554,055) covering 84% of the sequences (Supp. Fig. 3). Our study profiles a diverse transcriptomic landscape and of Ambystoma mexicanum protein-coding transcripts and microRNAs along internal organs and limbs. Some transcriptional profiles are observed to be similar in expression, specifically in heart, legs and tail, which share also tissue composition, while the profile of the liver is more specialized (Supp. Fig. 4). We defined tissue-enriched expression of both mRNAs and miRNAs, which agree with the specialized function of organs such as heart and liver. Also, when compared to human multi-tissue data sets (Melé et al., 2015), we found that the axolotl shares 950 protein-coding transcripts which expression is also specific in human liver and heart. Importantly, among the shared tissuespecific transcripts, 125 are involved in key developmental and physiological processes, while the perturbation of 344 of them in humans has been associated with disorders and congenital diseases. For example, we found among the most liver-enriched transcripts homologs of the human genes A2M, CRIM1 and SULT1A1. Indeed, AmA2M, the most liver-enriched transcript, codes for an alpha 2 macroglobulin that in humans has been characterized as a biomarker to predict liver fibrosis in hepatitis C patients (Ho et al., 2010). The Human homolog of AmCRIM1 encodes a Cysteine-rich motor neuron1 protein, an antagonist of bone morphogenetic proteins (BMPs), a key player in placental development, organogenesis, angiogenesis and kidney disease and it has been recently proposed as a potential risk 7

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Table 3 Tissue-enriched miRNAs. Enrichment is computed as the ratio of the tissue expression over the maximal expression in other tissues. Tissue expression (cpm) miRNA Liver Enriched and Highly Specific* amx-miR-142-3p_1 amx-miR-2970-5p amx-miR-425-5p_1 amx-miR-317 amx-miR-142-3p_2 amx-miR-34c amx-miR-3600 amx-miR-425-5p_2 amx-miR-122 amx-miR-1260 amx-miR-22-3p amx-miR-449a-5p amx-miR-11-3p amx-miR-8-5p amx-miR-216b amx-miR-217 amx-miR-22b-3p amx-miR-216a amx-miR-148a amx-miR-191 amx-miR-103b amx-miR-192-5p amx-miR-215-5p amx-miR-142b amx-miR-150–5p amx-miR-34c-3p amx-miR-126-3p amx-miR-148 amx-miR-194-5p amx-miR-22 amx-miR-223 amx-miR-122-3p amx-miR-3591 amx-miR-30a-5p amx-miR-192 amx-miR-142–5p Tail enriched and Highly specific* amx-miR-124 amx-miR-9-5p amx-miR-375_1 amx-miR-375_2 amx-miR-92 amx-miR-182 amx-miR-10c-5p amx-miR-9-7-3p amx-miR-10b-5p_1 amx-miR-10a amx-miR-375_3 amx-miR-124b-5p amx-miR-10b-5p_2 amx-miR-375_4 amx-miR-135a-1-3p amx-miR-96-5p amx-miR-219-3p

Liver

Front Leg

Tail

Rear Leg

Enrichment

54.74 1389.86 25.40 17.82 63.16 80.28 75787.61 1393.09 4394.57 254.17 425.11 41.54 45.19 58.24 *22.60 225.12 66.95 *364.20 121,205.10 41,216.47 2181.85 84,611.69 6669.60 38.31 1845.29 *33.68 17,819.80 13.61 707.07 10.95 1567.12 *203.22 267.78 53,720.04 57.26 2832.64

13.88 318.36 5.01 5.29 21.53 2.93 11312.91 168.10 16.05 22.19 36.26 7.93 14.35 7.18 0.00 3.21 10.86 2.74 25,132.24 9544.25 789.61 478.71 48.73 9.82 117.20 0.28 4579.56 3.78 14.07 1.13 232.51 0.28 1.51 13,567.84 0.00 919.75

7.85 511.44 16.13 12.64 3.05 1.74 49600.44 931.32 3.92 191.85 269.02 7.41 34.88 52.32 0.00 1.74 44.91 2.62 47,639.69 36,396.21 1594.06 474.82 31.83 6.10 831.48 0.00 7154.54 6.54 10.90 6.98 625.24 0.00 0.44 16,660.50 0.00 817.09

4.17 305.40 4.17 4.25 5.82 0.47 14301.22 256.95 7.55 41.37 40.11 4.17 8.34 9.36 0.08 1.18 14.87 1.81 22,885.24 12,466.77 621.10 261.67 17.85 5.58 322.31 0.39 3938.36 3.22 4.56 1.89 151.95 0.16 0.63 10,547.14 0.08 448.62

3.94 2.72 1.57 1.41 2.93 27.42 1.53 1.50 273.72 1.32 1.58 5.24 1.30 1.11 287.12 70.11 1.49 132.98 2.54 1.13 1.37 176.75 136.87 3.90 2.22 85.64 2.49 2.08 50.25 1.57 2.51 717.34 177.22 3.22 727.60 3.08

1.54 4.49 2.39 489.53 57.54 210.10 0.98 0.28 1.40 5905.83 995.35 0.14 1.54 16.98 0.28 0.00 0.00

5.67 25.40 1.04 262.73 16.15 2491.31 2.17 6.23 5.57 13,030.67 340.36 0.57 2.93 8.59 0.28 5.10 0.00

63.22 1164.59 34.88 9205.55 1784.17 53,101.63 44.04 205.80 87.64 298,700.07 13,702.59 *10.46 65.84 218.01 *11.34 69.76 *11.77

3.46 14.08 1.89 426.76 19.11 2814.67 2.44 3.15 5.90 13,834.98 510.36 0.24 2.52 7.79 0.31 3.30 0.00

11.16 45.84 14.62 18.80 31.01 18.87 18.06 33.02 14.86 21.59 13.77 18.47 22.49 12.84 36.03 13.68 –

From the total miRNAs, 82 showed a broad expression pattern among the 4 samples, while 53 showed a tissue-enriched expression pattern (Table 3). Tissue-specificity of axolotl miRNAs opens new paths to explore the function of these molecules in the context of a regenerative model system when miRNAs have homologs in humans which are also expressed in a tissue-specific manner. This is the case of Am-miR-122 that is highly expressed in both axolotl and human liver (Fig. 7, reviewed in (Cruz-Santos et al., 2016)). Indeed, this miR-122 has been used to reprogram hMSCs (human mesenchymal stem cells) to hepatocyte-like cells (reviewed in (Bandiera et al., 2015)). Among other liver-enriched miRNAs found in our study are Am-miR-223, AmmiR-425-5p and Am-miR-3600 (Fig. 7). miR-223 down-regulation is

diseases in the context of an organism that is able to regenerate. In conclusion, these novel findings can relaunch the axolotl as a model system for biomedical research with a broader spectrum than just regeneration. We also deep sequenced small RNAs (sRNAs) of 20–38 bp in size from heart, liver, rear leg and front leg. We then focused in the analysis of those with a length in which the majority could be putative miRNAs and identified a total of 352 miRNAs which have been previously annotated in miRBase, suggesting conservation in at least another species besides A. mexicanum. We also identified 297 mature miRNA sequences that did not map to miRBase, which could be considered as putative novel miRNAs. 8

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Fig. 8. Novel miRNAs predicted in Axolotl tissues. (A) Expression profiles of the most expressed novel mature miRNAs in different tissues (B) Secondary structure predicted for three of the novel miRNA, the mature sequence is highlighted in red.

associated with human hepatocellular carcinoma (Wong et al., 2008), while down-regulation of miR-425-5p expression inhibits proliferation, invasion and migration of gastric carcinoma mouse cells (Zhang et al., 2015). These findings position miR-425-5p as potential metastasisassociated gene and a target for gene therapy. The 297 novel putative microRNAs showed characteristics of known microRNAs in terms of structure and putative targets, however, the validation of this novel microRNAs is required and we plan the respective experiments in the future. When we compare our microRNA detection approach with other previous reports, we found that only 34 microRNAs are shared with those reported in King and Yin (2016), and 151 are shared with Gearhart et al. (2015). The differences among datasets can be explained by the fact that the methodology and experiment design are not similar. Such approaches during limb regeneration (King and Yin, 2016) and focusing on tail regeneration (Gearhart et al., 2015) are quite different from ours. Our approach focuses on profiling miRNAs from tissues/organs in a context where there is no developmental program, such as regeneration, driving differential gene expression as a consequence of cell reprogramming events. Moreover, ours is the first study that focuses on identifying also the pre-miRNAs sequences using non-coding transcripts from RNA-seq data from the same tissues where miRNAs were profiled. Besides our findings on the conserved miRNAs, by using Friedländer et al. algorithm (Friedländer et al., 2012), we predicted that the 297 non-previously annotated putative miRNAs based on: a) having a pre-miRNA with hairpin structure that has been proven nonrandom based on the software RandFold (Bonnet et al., 2004), b) the mature sequences align to a side of the hairpin and not the loop c) have sufficient counts in our miRNA-seq analyses to pass miRDeep2 score cut-off. The existence of a high amount of non-coding RNAs in our transcriptome landscape and the finding of putative novel miRNAs correlates well with the predicted genome size of A. mexicanum (approximately 32 Gbp). In fact, it has been previously reported that

Fig. 7. Predicted differential expression patterns of mature miRNAs and validation by qPCR. To the left, Graphics for each miRNA showing the expression levels as obtained from the miRNA-seq data of four different tissues (liver, tail, front leg and rear leg). To the right, graphics of each miRNA validated for three biological replicates by qPCR are expressed as relative quantification using SNU6 as endogenous control. data are the mean of two biological replicates, each one with three technical replicates.

9

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

The final assembly was annotated first identifying their main ORF and then comparing the peptide sequence with the human proteins (GRCh38) and NCBI Non-Redundant database (downloaded Nov. 2015) with Blast+ v2.2.30 (Camacho et al., 2008) keeping the best match with E-value < 1e-6.

the big size of the genome is not due only to the existence or more repetitive elements but also to the size of the intronic and intergenic regions, which is 10 times longer than orthologous vertebrate introns (Smith et al., 2009). This suggests that these big introns and intergenic regions might code for a population of small and long ncRNAs that might be more numerous than those expressed in vertebrate species with smaller genomes. However, several future analyses need to be done to validate such hypothesis. Finally, our curated and analyzed data set, based in a tissue/organenriched design, is an approach that improves the use of axolotl as a model system, not only to understand diverse developmental processes such as regeneration and metamorphosis, but also to comprehend the molecular mechanisms that rely behind specific genetic-based human diseases.

4.4. Identification of tissue-enriched genes The reference transcriptome assembled was used to quantify the expression in each tissue and discover which genes are tissue-specific or tissue-enriched. We aligned the original fragments to this transcriptome with Bowtie v1.1.1 (Langmead et al., 2009) keeping only the reads uniquely mapped to one unigene. This is a more conservative strategy considering alternatives such as RSEM or eXpress protocols which tends to improve quantification in transcriptomes with large number of gene isoforms. However, our transcriptome is still partial and the detection of full isoforms is not complete, therefore we consider such methods could not improve the results. Then we count the number of reads per transcript and create an expression table for all tissues. The data was normalized and transformed to count per million (cpm) with edgeR (Robinson et al., 2010), the average expression for any unigene was 3 cpm. We searched for enriched genes using a 4-fold enrichment over the maximal expression in other tissues with a minimal expression of 3 cpm as thresholds.

4. Materials and methods 4.1. Animal dissection and tissue processing Adult female Ambystoma mexicanum specimens from Xochimilco, which parental individuals were collected from the wild and did not suffer endogamy due to a careful ex situ rearing and reproduction process, were obtained from the Unit of Environmental Management at the Centro de Investigaciones Acuáticas de Cuemanco, Universidad Autónoma Metropolitana of México City, Campus Xochimilco. Two individuals with the same size, 18 cm in length head-tail, were anesthetized in 0.01% benzocaine and dissected to cut or extract the heart, liver, rear legs, front legs, gills and tail. The organs and limbs of one individual were processed for RNA extraction and this RNA was used for the generation of libraries of mRNAs and sRNAs. A similar process was followed with the organs and limbs of 2 different female individuals with similar size and from the same source and used for cDNA synthesis and q-PCR validations. All animal experiments were performed according to the Mexican Official Norm (NOM-062-ZOO-1999) “Technical Specifications for the Care and Use of Laboratory Animals” based on the Guide for the Care and Use of Laboratory Animals “The Guide”, 2011, NRC, USA, with the Federal Register Number # BOO.02.08.01.01.0095/2014, awarded by the National Health Service, Food Safety and Quality (SENASICASAGARPA). The Institutional Animal Care and Use Committee (IACUC) from the CINVESTAV approved the project “Management and husbandry of Ambystoma Spp and experimental processing of tissue for functional analyses and genetic expression” ID animal use protocol number: 0209-16.

4.5. Human and Axolotl disease relationship approach The enriched gene list from AmRNAs was converted to the equivalent human ortholog genes searching for ortholog genes in a best-match with a translated search with Blast+ (blastx) to the human proteome obtained from Ensembl r86 for GRCh38, keeping only best hit with E-value < 0.000001. The human ortholog genes were analyzed with clusterProfiler (Yu et al., 2012) and DOSE (Yu et al., 2015) packages in R using a q-value threshold of 0.05 for disease-enrichment analysis. 4.6. sRNA library generation and sequencing Total RNAs extracted from liver, rear-leg, front-leg and gills for library generation and sequencing, was processed to obtain small RNAs in the desired size range of 18–34 nt, which were recovered by size fraction on a 15% Urea Polyacrylamide gel. Library construction was performed using the Illumina Small RNA Sample Preparation kit per manufacturer's instructions. In brief, small RNA fragments were ligated with 3′ and 5′ adapters sequentially, followed by a reverse transcription reaction and a 15-cycle PCR reaction. PCR products with the expected size were excised and eluted from a 6% Urea Polyacrylamide gel, and pooled for sequencing on an Illumina HiSeq. 2000. The different samples were decoded during the bioinformatics analysis stage based on the barcodes in the primers.

4.2. RNA-seq library generation and sequencing Total RNA was extracted using TRIZOL reagent (Invitrogen). High quality RNAs were used for library preparation after evaluated using bioanalyzer: 5067-1511 Agilent RNA 6000 Nano Kit. All libraries were generated using TruSeq protocols and sequenced in one lane using the 2 × 100 Illumina HiSeq. 2000 platform. Libraries and sequencing were generated at the LANGEBIO genomics services.

4.7. qPCR for mRNA quantification Total RNA was extracted using Trizol Reagent (Invitrogen). Total RNA samples of the different organs from two adult female axolotls (front leg, rear leg, tail, gill, liver and heart) was reverse-transcribed to first-strand cDNA with reverse transcriptase was done according to the manufacturer's protocol High Capacity RNA to cDNA kit (Applied Biosystems), we used 5 μg of total RNA per 20-μL reaction. The expression levels of specific mRNAs were determined by qPCR using SYBR GREEN PCR master mix in 20-μL reaction, the primer sequences used in these experiments are enlist in Supp. Table 4. To quantify the relative expression level we did two biological replicates with three independent qPCR reactions for each data point and we included the efficiency in the calculations, all data points were normalized to GAPDH expression levels.

4.3. Transcript assembly and annotation We used Trinity v2.0.6 (Grabherr et al., 2011) to reconstruct the transcriptomes individually, given the high computer resources required to perform an assembly using all samples. Each sample was assembled using default parameters but including the “–trimmomatic” step to remove sequencing adapters and low quality fragments based on their base quality. Each set (six from our tissues and one extra from the combined limb regeneration time course) was assembled and the resulting unigenes were combined in one single transcriptome using CD-HIT-EST v4.6 (Limin et al., 2012) using global alignment with 0.9 identity to build clusters and define the best representative consensus sequence. 10

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

other individuals). The total RNA from front leg, rear leg, tail and liver was used for miRNA expression quantification. The first-strand cDNA synthesis and qRT-PCR were done according to manufacturer's protocol NCode kit (Invitrogen). A Step-One Real-Time PCR system (Applied Biosystems) was used, the primer sequences used in these experiments are listed in Supp. Table 4. To quantify the relative expression levels we did 3 biological replicates with 3 independent qPCR reactions for each data point and we included the efficiency in the calculations, all data points were normalized to SNU6 expression levels.

4.8. MicroRNA identification and annotation Adapters were removed using Cutadapt v1.5 (Martin, 2011) under the parameters “-e 0.05 -O 10 -m 15”, then low quality reads were removed using Trimmomatic v0.33 (Bolger et al., 2014) with parameters “SLIDINGWINDOW:4:22 MINLEN:15”. miRBase (Kozomara and Griffiths-Jones, 2014) version 21 was used to identify conserved miRNAs. Redundancy in the database was suppressed by removing sequences that by taking 2 nucleotides at most were identical to another one under the same conditions. The identification of miRNA was performed keeping the best hits with BLAST (with the parameters “-reward 5 -penalty −4 -gapopen 25 -gapextend 10 -word size 7”). Only those miRNAs with a perfect match in at least 18 nt were considered as conserved miRNAs. Lastly, miRNAs were collapsed after mapping with BLAST by removing 1 nucleotide in the 3′ hang, and checking if it was identical to others. Finally, we set a final cut-off of 20 reads on at least one sample to consider the miRNA as present in the samples.

Declarations Ethics approval and consent to participate Not applicable. Consent for publication

4.9. Pre-microRNA identification Not applicable. Non-coding sequences of the transcriptome were used as template, using the software mirDeep2 v2.0.0.7 (Friedländer et al., 2012) to find conserved and novel precursors. Mappings of the sRNA against the non-coding sequences of the transcriptome was performed with Bowtie (Langmead et al., 2009) (parameters: “-n 0 -l 15”). All the reads from the sRNA samples were put into one file to obtain global results in one run and to maximize the possibility of finding all pre-miRNA. All the precursors with a mirDeep2 score < 5 or coming from isoforms from the same transcript were discarded.

Availability of data and material The datasets generated during the current study are available in the NCBI SRA repository, under accession SRP093628. mRNA expression patterns per tissue can be searched at http://www.cruzomics.net/ ambystoma-mexicanum/ Competing interests

4.10. MicroRNA target prediction

The authors declare that they have no competing interests

The novel putative miRNAs were used to identify potential target sequences aligning the mature miRNA to the reconstructed transcriptome. We used MiRanda (Enright et al., 2003) to predict potential targets using the following parameters “-strict -E -20” and filter the hits with a strong 5′ seed and a maximal free energy of −20 kcal/mol. An additional filter keeps only miRNA:Targets with more than one binding site to increase our prediction confidence. Potential targets were mapped to the human proteome by using the best translated hit with Blastx (Camacho et al., 2008), keeping only the best hit with an Evalue < 1e-6 and analyzed with R: Bioconductor packages clusterProfiler (Yu et al., 2012) selecting significant Gene Ontologies categories with q-value < 0.05.

Funding This work was supported by The Swedish International Research links grant 2014-9040-114152-32, the UC MEXUS Collaborative Grant 2011-UC MEXUS-19941-44-OAC7, Consejo Nacional de Ciencia y Tecnología grants: Fronteras de la Ciencia Grants CONACYT FOINS301 and Ciencia Básica I0017-CB-2015-01-000000000252126. Authors' contributions Conceived and designed the experiments: ACR, LHE, JCP, LFGO, MAAV, AEC. Performed the experiments: ACR, AEC, MAAV, XC, ACH. Analyzed the data: JCP, LFGO, ACR, OMV, AEC, XC, MAAV, FFC and LB. Contributed reagents/materials/analysis tools: LHE, XC, MAAV. Wrote the paper: JCP and ACR.

4.11. Differential expression analyses of microRNAs We used the method proposed by García-Ortega and Martínez (2015) to estimate overdispersion in our RNA-seq experiments without biological replicates. This approach is based in the estimation of over dispersion in a set of housekeeping genes, selected by their low specificity (Martínez and Reyes-Valdés, 2008) and logarithmic ratio of the maximum likelihood statistics (Stekel et al., 2000). This approximation allows us to correctly determine at least a subset of the differential expressed mature miRNA's and pre-miRNA's. Differential expression analyses were carried out by using the R package edgeR v3.8.6 (Langmead et al., 2009), using the common overdispersion values obtained as proposed by Robinson and Smyth (2008) for differential expression analysis for RNA-seq experiments without biological replicates.

Acknowledgements The authors acknowledge Profr. Fernando Arana-Magallón and Biol. Arturo Vergara-Iglesias from (CIBAC-UAM Xochimilco) for their support in animal management and UMA San Miguel Tecocomulco for reproduction and conservation efforts. Dr. Enrique Ibarra-Laclette (LANGEBIO-CINVESTAV) for preliminary data analyses (non-published), Biol. Fernando Hernández (LANGEBIO-CINVESTAV) for technical support, Mr. Jaime Calzada (LANGEBIO-CINVESTAV) for support in photographic work. ACR wish to dedicate this manuscript to Tepeapulco, Hidalgo; México, the Birthplace of Anthropologic studies in the American continent and the home of Ambystoma velasci.

4.12. qPCR for miRNA quantification Appendix A. Supporting information Total RNA was extracted using TRizol Reagent (Invitrogen). We used three biological replicates (we kept small RNA fraction of our original miRNA RNAseq experiment and included RNAs from two

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.ydbio.2017.08.022. 11

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

Limin, F., Beifang, N., Zhengwei, Z., Sitao, W., Weizhong, L., 2012. CD-HIT: accelerated for clustering the next generation sequencing data. Bioinformatics 28 (23), 3150–3152. http://dx.doi.org/10.1093/bioinformatics/bts565. Martin, M., 2011. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17 (1), 10–12. http://dx.doi.org/10.14806/ ej.17.1.200. Martínez, O., Reyes-Valdés, M.H., 2008. Defining diversity, specialization, and gene specificity in transcriptomes through information theory. Proc. Natl. Acad. Sci. USA 105 (28), 9709–9714. Melé, M., Ferreira, P.G., Reverter, F., DeLuca, D.S., Monlong, J., Sammeth, M., Young, T.R., Goldmann, J.M., Pervouchine, D.D., Sullivan, T.J., et al., 2015. The human transcriptome across tissues and individuals. Science 348, 660–665. http:// dx.doi.org/10.1126/science.aaa0355. Misra, C., Sachan, N., Rothrock, C.M., Koenig, S.N., Nichols, H.A., Guggilam, A., Lucchesi, P.A., Pu, W.T., Srivastava, D., Garg, V., 2012. Congenital heart diseasecausing Gata4 mutation displays functional deficits in vivo. PLoS Genet. 8 (5), e1002690. http://dx.doi.org/10.1371/journal.pgen.1002690. Monaghan, J.R., Walker, J.A., Page, R.B., Putta, S., Beachy, C.K., Voss, S.R., 2007. Early gene expression during natural spinal cord regeneration in the salamander Ambystoma mexicanum. J. Neurochem. 101, 27–40. Monaghan, J.R., Epp, L.G., Putta, S., Page, R.B., Walker, J.A., Beachy, C.K., Zhu, W., Pao, G.M., Verma, I.M., Hunter, T., et al., 2009. Microarray and cDNA sequence analysis of transcription during nerve-dependent limb regeneration. BMC Biol. 13 (7:1). http://dx.doi.org/10.1186/1741-7007-7-1. Monaghan, J.R., Athippozhy, A., Seifert, A.W., Putta, S., Stromberg, A.J., Maden, M., Gardiner, D.M., Voss, S.R., 2012. Gene expression patterns specific to the regenerating limb of the Mexican axolotl. Biol. Open 15, 937–948. http://dx.doi.org/ 10.1242/bio.20121594. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628. Ochoa, L.M.O., Villela, O.A.F., 2006. Áreas de diversidad y endemismo de la herpetofauna mexicana. UNAM, Mexico. Parra‐Olea, G., Recuero, E., Zamudio, K.R., 2007. Polymorphic microsatellite markers for Mexican salamanders of the genus Ambystoma. Mol. Ecol. Note 7 (5), 818–820. Peck-Radosavljevic, M., 2007. Coagulation disorders in chronic liver disease. Aliment. Pharmacol. Ther. 26 (s1), 21–28. Putta, S., Smith, J.J., Walker, J.A., Rondet, M., Weisrock, D.W., Monaghan, J., Samuels, A.K., Kump, K., King, D.C., Maness, N.J., et al., 2004. From biomedicine to natural history research: EST resources for ambystomatid salamanders. BMC Genom. 13, 54 . Rand, T.A., Ginalski, K., Grishin, N.V., Wang, X., 2004. Biochemical identification of Argonaute 2 as the sole protein required for RNA-induced silencing complex activity. Proc. Natl. Acad. Sci. USA 101 (40), 14385–14389. Robinson, M.D., Smyth, G.K., 2008. Small sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics 9, 321–332. Robinson, M.D., McCarthy, D.J., Smyth, G.K., 2010. EdgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. Roy, S., Gatien, S., 2008. Regeneration in axolotls: a model to aim for!. Exp. Gerontol. 43, 968–973. Selcen, D., Engel, A.G., 2005. Mutations in ZASP define a novel form of muscular dystrophy in humans. Ann. Neurol. 57 (2), 269–2676. Smith, J.J., Voss, S.R., 2006. Gene order data from a model amphibian (Ambystoma): new perspectives on vertebrate genome structure and evolution. BMC Genom. 29 (7), 219. Smith, J.J., Putta, S., Zhu, W., Pao, G.M., Verma, I.M., Hunter, T., Bryant, S.V., Gardiner, D.M., Harkins, T.T., Voss, S.R., 2009. Genic regions of a large salamander genome contain long introns and novel genes. BMC Genom. 13, 10–19. http://dx.doi.org/ 10.1186/1471-2164-10-19. Stekel, D.J., Git, Y., Falciani, F., 2000. The comparison of gene expression from multiple cDNA libraries. Genome Res. 10 (12), 2055–2061. Stewart, S., Alexander-Rascon, C., Tian, S., Nie, J., Barry, C., Chu, L.F., Ardalani, H., Wagner, R.J., Probasco, M., Bolin, J.M., et al., 2013. Comparative RNA-seq analysis in the unsequenced Axolotl: the oncogene burst highlights early gene expression in the blastema. PLoS Comput. Biol. 9, e1002936. http://dx.doi.org/10.1371/ journal.pcbi.1002936. Veugelers, M., Bressan, M., McDermott, D.A., Weremowicz, S., Morton, C.C., Mabry, C.C., Lefaivre, J.F., Zunamon, A., Destree, A., Zunamon, A., et al., 2004. Mutation of perinatal myosin heavy chain associated with a carney complex variant. N. Engl. J. Med. 351 (5), 460–469. Volk, N., Shomron, N., 2011. Versatility of MicroRNA biogenesis. PLoS One 6 (5), e19391. Voss, S.R., Smith, J.J., Gardiner, D.M., Parichy, D.M., 2001. Conserved vertebrate chromosome segments in the large salamander genome. Genetics 158, 735–746. Voss S.R., Epperlein H.H., Tanaka E.M. Ambystoma mexicanum, the Axolotl: A versatile amphibian model for regeneration, development, and evolution studies. In: Emerging Model Organisms: A laboratory manual- Vol. 2. Cold Spring Harbor Press, 2009, pp. 397–413. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10 (1), 57–63. Wheeler, B.M., Heimberg, A.M., Moy, V.N., Sperling, E.A., Holstein, T.W., Heber, S., Peterson, K.J., 2009. The deep evolution of metazoan microRNAs. Evol. Dev. 11 (1), 50–68. Wong, Q.W., Lung, R.W., Law, P.T., Lai, P.B., Chan, K.Y., Wong, K.F., 2008. N.MicroRNA-223 is commonly repressed in hepatocellular carcinoma and potentiates expression of Stathmin1. Gastroenterology 135 (1), 257–269.

References Alvarez-Garcia, I., Miska, E.A., 2005. MicroRNA functions in animal development and human disease. Development 132 (21), 4653–4662. Bandiera, S., Pfeffer, S., Baumert, T.F., Zeisel, M.B., 2015. MiR-122 - a key factor and therapeutic target in liver disease. J. Hepatol. 62 (2), 448–457. Bolger, A.M., Lohse, M., Usadel, B., 2014. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30 (15), 2114–2120. http://dx.doi.org/10.1093/ bioinformatics/btu170. Bonnet, E., Wuyts, J., Rouzé, P., Van, de Peer, Y., 2004. Evidence that microRNA precursors, unlike other non-coding RNAs, have lower folding free energies than random sequences. Bioinformatics 20 (17), 2911–2917. Bryant, D.M., Johnson, K., DiTommaso, T., Tickle, T., Couger, M.C., Payzin-Dogru, D., Lee, T.J., Leigh, N.D., Kuo, T., Davis, F.G., et al., 2017. Tissue-mapped axolotl de novo transcriptome enables identification of limb regeneration factors. Cell Rep. 18 (3), 762–776. Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., Madden, T.L., 2008. BLAST+: architecture and applications. BMC Bioinform. 10, 421. Casey, M.C., Kerin, M.J., Brown, J.A., Sweeney, K.J., 2015. Evolution of a research field a micro(RNA) example. Peer J. 3, e829. Cruz-Santos, M., Aragón-Raygoza, A., Espinal-Centeno, A., Arteaga-Vázquez, M., Bako, L., Cruz-Hernández, A., Cruz-Ramirez, A., 2016. The role of microRNAs in animal cell reprogramming. Stem Cells Dev. 25 (14), 1035–1049. http://dx.doi.org/ 10.1089/scd.2015.0359. Denzler, R., Stoffel, M., 2015. The long, the short, and the unstructured: a unifying model of miRNA biogenesis. Mol. Cell 60 (1), 4–6. Enright, A.J., John, B., Gaul, U., Tuschl, T., Sander, C., Marks, D.S., 2003. MicroRNA targets in Drosophila. Genome Biol. 5 (1), 1. Friedländer, M.R., Mackowiak, S.D., Li, N., Chen, W., Rajewsky, N., 2012. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40 (1), 37–52. Friedman, R.C., Farh, K.K., Burge, C.B., Bartel, D.P., 2009. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 19 (1), 92–105. García-Ortega, L.F., Martínez, O., 2015. How many genes are expressed in a transcriptome? Estimation and results for RNA-seq. PLoS One 10 (6), e0130262. Gearhart, M.D., Erickson, J.R., Walsh, A., Echeverri, K., 2015. Identification of conserved and novel MicroRNAs during tail regeneration in the mexican axolotl. Int. J. Mol. Sci. 16 (9), 22046–22061. http://dx.doi.org/10.3390/ijms160922046. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., et al., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29 (7), 644–652. http://dx.doi.org/10.1038/nbt.1883. Gregory, R.I., Chendrimada, T.P., Cooch, N., Shiekhattar, R., 2005. Human RISC couples microRNA biogenesis and posttranscriptional gene silencing. Cell 123 (4), 631–640. Gwizdek, C., Ossareh-Nazari, B., Brownawell, A.M., Doglio, A., Bertrand, E., Macara, I.G., Dargemont, C., 2003. Exportin-5 mediates nuclear export of minihelixcontaining RNAs. J. Biol. Chem. 278 (8), 5505–5508. Habermann, B., Bebin, A.G., Herklotz, S., Volkmer, M., Eckelt, K., Pehlke, K., Epperlein, H.H., Schackert, H.K., Wiebe, G., Tanaka, E.M., 2004. An Ambystoma mexicanum EST sequencing project: analysis of 17,352 expressed sequence tags from embryonic and regenerating blastema cDNA libraries. Genome Biol. 5, R67. Ho, A.S., Cheng, C.C., Lee, S.C., Liu, M.L., Lee, J.Y., Wang, W.M., Wang, C.C., 2010. Novel biomarkers predict liver fibrosis in hepatitis C patients: alpha 2 macroglobulin, vitamin D binding protein and apolipoprotein AI. J. Biomed. Sci. 17, 58. http:// dx.doi.org/10.1186/1423-0127-17-58. Huggins, P., Johnson, C., Schoergendorfer, A., Putta, S., Bathke, A.C., Stromberg, A.J., Voss, S.R., 2012. Identification of differentially expressed thyroid hormone responsive genes from the brain of the Mexican Axolotl (Ambystoma mexicanum) Comparative biochemistry and physiology. Toxicol. Pharmacol.: CBP 55 (1), 128–135. http://dx.doi.org/10.1016/j.cbpc.2011.03.006. Khattak, S., Murawala, P., Andreas, H., Kappert, V., Schuez, M., Sandoval-Guzmán, T., Crawford, K., 2014. Tanaka EM.optimized axolotl (Ambystoma mexicanum) husbandry, breeding, metamorphosis, transgenesis and tamoxifen-mediated recombination. Nat. Protoc. 9 (3), 529–540. http://dx.doi.org/10.1038/ nprot.2014.040. Kim, V.N., Nam, J.W., 2006. Genomics of microRNA. Trends Genet. 22 (3), 165–173. King, B.L., Yin, V.P., 2016. A conserved microRNA regulatory circuit is differentially controlled during limb/appendage regeneration. PLoS One 11 (6), e0157106. http:// dx.doi.org/10.1371/journal.pone.0157106. Knapp, D., Schulz, H., Rascon, C.A., Volkmer, M., Scholz, J., Nacu, E., Le, M., Novozhilov, S., Tazaki, A., Protze, S., et al., 2013. Comparative transcriptional profiling of the axolotl limb identifies a tripartite regeneration-specific gene program. PLoS One 1, e61352. http://dx.doi.org/10.1371/journal.pone.0061352. Kozomara, A., Griffiths-Jones, S., 2014. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42 (D1), D68–D73. Kragl, M., Knapp, D., Nacu, E., Khattak, S., Maden, M., Epperlein, H.H., Tanaka, E.M., 2009. Cells keep a memory of their tissue origin during axolotl limb regeneration. Nature 460, 60–65. Lagos-Quintana, M., Rauhut, R., Lendeckel, W., Tuschl, T., 2001. Identification of novel genes coding for small expressed RNAs. Science 294 (5543), 853–858. Langmead, B., Trapnell, C., Pop, M., Salzberg, S.L., 2009. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. Lee, Y.S., Nakahara, K., Pham, J.W., Kim, K., He, Z., Sontheimer, E.J., Carthew, R.W., 2004. Distinct roles for Drosophila Dicer-1 and Dicer-2 in the siRNA/miRNA silencing pathways. Cell 117 (1), 69–81.

12

Developmental Biology xxx (xxxx) xxx–xxx

J. Caballero-Pérez et al.

biological themes among gene clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. Yu, G., Wang, L., Yan, G., He, Q., 2015. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics 31 (4), 608–609. Zeng, H., Tang, L., 2014. CRIM1, the antagonist of BMPs, is a potential risk factor of cancer. Curr. Cancer Drug Targets 14 (7), 652–658. Zhang, Z., Li, Y., Fan, L., Zhao, Q., Tan, B., Li, Z., Zang, A., 2015. MicroRNA-425-5p is upregulated in human gastric cancer and contributes to invasion and metastasis in vitro and in vivo. Exp. Ther. Med. 9, 1617–1622.

Xie, B., Ding, Q., Han, H., Wu, D., 2013. miRCancer: a microRNA–cancer association database constructed by text mining on literature. Bioinformatics, Btt014. Yao-Borengasser, A., Rogers, L.J., Edavana, V.K., Penney, R.B., Yu, X., Dhakal, I.B., Williams, S., Kadlubar, S.A., 2014. Sulfotransferase 1A1 (SULT1A1) gene expression is regulated by members of the NFI transcription factors in human breast cancer cells. BMC Clin. Pathol. 14 (1), 1. http://dx.doi.org/10.1186/1472-6890-14-1. Yi, R., Qin, Y., Macara, I.G., Cullen, B.R., 2003. Exportin-5 mediates the nuclear export of pre-microRNAs and short hairpin RNAs. Genes Dev. 17 (24), 3011–3016. Yu, G., Wang, L., Han, Y., He, Q., 2012. ClusterProfiler: an R package for comparing

13