MicroRNAs (miRNA) are 18-24 nucleotide (nt) long eukaryotic noncoding sequences that act as post-transcriptional regulators of gene and gene networks through either mRNA degradation or translational repression. These small RNAs are generated from their source genes mostly through bi-directional transcription by RNA polymerase-II. This generates a primary miRNA (pri-miRNA) which is then processed by the RNAse-III-like enzyme, Drosha, in the nucleus. The pri-miRNA is then processed into a ~70-nt long hairpin-shaped precursor miRNA (pre-miRNA), which is transported out of the nucleus by the Exportin-5 protein. Once inside the cytoplasm, the pre-miRNA undergoes a round of Dicer-mediated cleavage to generate the duplex miRNA. The more thermodynamically stable mature miRNA strand is bound by an Argonaute (AGO) protein and forms the functional ribonucleoprotein complex known as the RISC (RNA-induced silencing complex). The RISC is then guided by the mature miRNA to a complimentary or partially complimentary mRNA target site to carry out gene silencing either through mRNA cleavage or translational repression [1][2][3]. The less stable complement to the mature miRNA will usually undergo degradation and be lost [4][5][6].
miRNAs in animals are known to be most effective in targeting gene transcripts when there is full complementarity between the 3′ untranslated region (UTR) of the target mRNA with a 7-8 nt stretch in the 5' end of miRNA (i.e., the seed sequence) [7][8][9][10]. Nearly 60% of the protein coding genes are subject to miRNA regulation [11]. In animals, https://doi.org/10.1016/j.ab.2020.113781 Received 12 February 2020; Received in revised form 11 May 2020; Accepted 19 May 2020 miRNAs are known to be involved in regulation of several critical biological processes-including the time and pattern of development, apoptosis, control of the cell cycle, disease pathogenesis and host response, and stress response [12][13][14][15].
Birds are the closest extant relatives to crocodilians in the clade Archosauria, which also encompasses non-avian dinosaurs and pterosaurs [16,17]. Crocodilian genomes have been evolving very slowly over the past several million years when compared to birds [18]. Understanding the evolution, regulation and adaptive capabilities of crocodilian genomes and their genetic diversity can therefore provide information on how slowly evolving genomes manage to stay viable in the face of competition from other taxa and changing environmental conditions. Reptilian genomic studies are also necessary for better understanding of the patterns of genomic evolution across amniotes (mammals, birds and non-avian reptiles).
miRNAs provide us with the opportunity to investigate the potential causes behind the slow evolution of crocodilians. At present, only two studies, both done in American alligators (Alligator mississippiensis), have examined small RNA dynamics in crocodilians. Rice et al. [19] targeted a related crocodilian, Alligator mississippiensis, but only a single tissue. In an earlier study, some conserved miRNAs, also in Alligator mississippiensis, were annotated by mapping small RNA reads to miRNAs from other species [20]. For this study, we used the most recent annotated genome assembly for the saltwater crocodile, Crocodylus porosus [21], to identify and characterize 644 miRNAs using Ilumina's® deep sequencing HiSeq 2500 platform. The miRNAome includes known and novel miRNAs. Tissue specific patterns of expression were also observed, providing the first comprehensive analysis of crocodilian miRNAs.
In addition to examining the general catalog of these small RNAs, we were interested in the genomic context of miRNAs. One feature that reflects the low rate of genomic change in crocodilians is the lack of recent transposable element (TE) accumulation in crocodilian genomes [18,22]. Transposable elements are genetic features that can copy and, in some cases, physically move from one location to another in a genome. Such activity can result in the modification of gene expression through multiple processes including disruption of coding and noncoding regions, promotion of non-homologous recombination, exaptation, and epigenetic regulation. All of these processes could impact the evolution of diversity, both taxonomic and genomic [23][24][25][26][27][28][29][30][31]. TEs are broadly classified into two classes, Class I retrotransposons and Class II DNA transposons [32]. An increasing body of research has suggested that TEs are major players in the evolution of gene regulation [33][34][35][36][37] and miRNAs have been known to be either derived from or involved in the silencing of TEs [38][39][40][41][42][43][44][45]. Because better understanding of regulatory networks and TEs can help us assess how gene expression phenotypes evolve and the TE landscape is maintained in the genome, we chose to investigate the connection between these two genomic features.
Ten hatchling saltwater crocodiles (of which five were phenotypically normal and five were runts) were euthanized for immediate tissue collection (University of Sydney Animal Ethics protocol N00/5-2012/ 3/5729) using methods previously described [46,47]. All tissues (Table 1) were flash frozen in liquid nitrogen. Tissues included were belly skin, brain, cloaca, heart, jaw skin, liver, small intestine, stomach, spleen, testis, and tongue. Normal vs. runt status of subjects was determined based on size differences. The five runts were, on average, 359 days old with a total length of 386.2 ± SD 23.6 mm compared to the normal, healthy crocodiles, which were 373 days old and with an average length of 740.2 ± 73.5 mm. These animals were denoted R1-R5 and N1-N5 respectively.
Total RNA was isolated using Trizol (Life Technologies Inc.) and the Direct-Zol RNA Miniprep Kit (Zymo Research Inc.). Total RNA was quantified using the Qubit RNA assay and Qubit 2.0 fluorometer (Life Technologies Inc.). RNA integrity was confirmed using the Total RNA Standard Sensitivity assay reagents and Fragment Analyzer system (Agilent Inc.). Only samples having an RNA Integrity Number (RIN) greater than 7 or samples with clear 18s and 28s ribosomal RNA peaks as seen with the Bioanalyzer were used in further studies, the only exception being one sample each from the brain and the heart where a RIN of < 7 was considered acceptable. Note that some of the brain and heart samples had RINs > 7 and these were compared to RIN < 7 samples later in the procedure.
One microgram of total RNA was used as input material in building each NEBNext Small RNA Library Prep for Illumina (NEB Inc.). To prepare samples for sequencing, small non-coding RNA species underwent 3′ RNA ligation and reverse transcriptase primer hybridization followed by 5' RNA ligation and first strand cDNA generation using Superscript III reverse transcriptase (Life Technologies Inc.). The resulting cDNA which contained a portion of the adapter sequences for Illumina sequencing, underwent PCR enrichment (94 °C for 30 s, 12 cycles [94 °C for 15sec, 62 °C for 30 s, 70 °C for 15 s], 70 °C for 5 min) using a universal primer and a multiplex specific primer. The final product was small RNA libraries with sample-specific 6-nt indexes that allow for multiplex sequencing on the Illumina platform. Following PCR, samples were purified using a 1.8x SPRI bead purification (Beckman Coulter Inc) and the final sequencing construct was size selected on the Pippin Blue system (Sage Science Inc.) using a 3% w/v gel cassette and internal marker for collecting the range of fragments from 135 to 170 base pairs (bp).
Library validation was performed using the NGS Hi Sensitivity assay reagents and the Fragment Analyzer system followed by quantitation using the Qubit HS DNA assay and a qPCR kit for Illumina (Kapa Biosystems Inc). Libraries were diluted using the Qubit or qPCR information and loaded onto a single HiSeq2500 (Illumina Inc.) flow cell. Single-end 75 bp reads were generated using the pooled libraries. Pooled sequence data was demultiplexed into individual files with bcl2fastq (v2) representing the sixteen unique libraries. Using a series of programs available in the FASTX-toolkit (v 0.0.14; http://hannonlab. cshl.edu/fastx_toolkit/), sequences were quality filtered. In brief, the Fastx_clipper program was used to remove adapter sequence (5′-AGA TCG GAA GAG CAC ACG TCT GAA CTC CAG TCA C-3′) from each read as well as culling reads with ambiguous nucleotides ("N"). Reads which lacked adapter sequence or were less than 10-nt after clipping were not included in downstream analyses. These reads were then piped into the Fastq_quality_filter program and any reads with Phred quality values less than 20 across more than half of the read were removed. Small RNA reads have been deposited in the Gene Expression Omnibus archive and can be accessed through GEO accession number GSE150461.
The miRNAome in C. porosus was identified and characterized in three phases: Identification of conserved miRNAs, identification of novel crocodilian miRNAs, and characterization of expression profiles. After initial library preparation, sequencing and filtering, miRDeep2 (version 2_0_0_8) [48] was used to identify and annotate both known and novel miRNA. Because of the relatively limited current knowledge of crocodilian miRNA, with no information available for the saltwater crocodile, two birds-the chicken and zebra finch (the only two birds with data in miRBase v21), were chosen as taxa for comparison. The mature and hairpin precursor sequences for the above avian miRNAs, as downloaded from miRBase (v21; www.mirbase.org), were used as a training set after using a custom UNIX script (https://github.com/ davidaray/cPorosus_miRNAs) to select chicken and zebra finch mature miRNAs and corresponding hairpin precursor sequences. This script was also used to process all remaining steps in our pipeline unless otherwise indicated. All mature and hairpin precursor sequences were concatenated for both species in two separate files and denoted AVIAN_ MATURE and AVIAN_HAIRPIN, respectively, to be used as training sets in the miRDeep2 analysis. The mapper.pl script (part of the miRDeep2 suite) was then used to map the deep sequencing reads to a reference genome. Default options for mapper.pl were used except reads less than 18 bases were removed (-l 18), identical reads were collapsed (-m), and sequences with non-canonical bases were removed (-j). miRDeep2 predicted miRNAs were scored based on the ability of the flanking genomic sequence to form a hairpin structure typical of miRNAs. Throughout the rest of this work, any miRNA identified in the saltwater crocodile but previously known and catalogued in miRbase were considered as "known", while those miRNAs that were not found in miRbase were designated as "novel" in crocodilians.
The miRDeep2 core algorithm assigns a distinct "score" to each mapped small RNA read based on its conformity to several miRNA biogenesis hallmarks. Based on a previous study [48], any miRNA with a miRDeep2 score of < 2 was removed from further analyses. After prediction by the miRdeep2 algorithm, several steps were taken to remove low confidence predictions from the dataset and subsequent analyses. First, all predicted known and novel miRNA hairpin precursor sequences were filtered to remove matches to tRNA and/or rRNA datasets. Any matches to either of these two datasets irrespective of evalue or bitscore, were removed from the predicted entries. The eukaryotic genomic tRNA database (http://gtrnadb.ucsc.edu) and the rRNA dataset (http://www.arb-silva.de/) were used. BLASTn in conjunction with our UNIX script was used to select out unique predicted miRNA entries that were produced from the miRDeep2.pl run but did not match to the rRNA and/or tRNA datasets.
There were some cases where multiple hairpin precursor sequences were identified from overlapping regions, often offset by only a few bases. We removed any overlapping miRNA predictions. The procedure was accomplished in two steps using Bedtools (v 2.26.0) [49]. First, the 'sort' and 'merge' functions from Bedtools were used to identify overlaps between the predicted miRNA precursor sequences. Only one instance was found where a novel and a known predicted miRNA candidate precursor overlapped and that pair was removed from both lists. For all remaining precursor candidates, the script re-moveOverlappingMirnas.pl was used to retain the highest scoring partner from each overlapping pair and remove the partners with lower scores. The score here refers to a variable representing the difference between start and end coordinate of the miRNA precursor. Essentially, any shorter precursor(s) completely overlapping and "contained within" another precursor (hence lower score) was discarded. The longer precursor (essentially the highest scoring partner) was retained.
The miRDeep2 core algorithm assigns a distinct "score" to each mapped small RNA read based on its conformity to several miRNA biogenesis hallmarks. Based on a previous study [48], a miRDeep2 score of 2 was chosen as a cutoff value to remove any candidate miRNA that had a miRDeep2 score of < 2.
The miRDeep2 algorithm utilizes the RNAfold (v 2.0) program to predict secondary structures of putative miRNA precursors, where the core algorithm assigns a "yes" flag to the "randfold" value for potential true miRNA candidate precursors. At this step, another custom bash script was used to select those miRNA candidates that were assigned a "yes" on the significant randfold value category by the miRDeep2.pl script. These selected miRNA candidates were then piped into the output file of the miRDeep2score filter from the previous step and used to generate a preliminary list of high-quality filtered miRNA candidates.
We removed any candidate mature miRNAs that mapped fewer than 10 times on the hairpin precursor sequence. The output file was ultimately piped into genHQseq.pl script to generate a final list of highquality mature miRNA and their corresponding hairpin precursor files in FASTA format.
Final manual filtering steps were included as follows. Any candidate miRNA with ten or more single nucleotide repeats were removed irrespective of prediction by the miRDeep2 program. Ten such candidate miRNAs were removed and none of these candidate miRNAs were from the previously 'known' category. The remaining candidates were scanned for duplicates (generated due to biogenesis from multiple loci in the genome or being retained as a result of existing in both chicken and zebra finch) and only one unique candidate sequence, based on its highest normalized expression was retained. All retained predicted miRNA candidates had at least a 70% estimated probability of being a true miRNA and a final total of 644 unique high-quality candidate miRNAs along with their hairpin precursors were identified and annotated in C. porosus. This included 132 known and 532 novel miRNA candidates. A general overview representing key steps of the miRNA prediction and annotation procedure is provided in Fig. 1.
Quality filtered small RNA reads from the individual libraries were mapped to the curated precursor and mature miRNA sequences. Once mapped, the quantifier.pl script from the miRDeep2 package was used to calculate both raw counts and normalized counts for each of the distinct mature miRNA candidates. Expression scores were weighted based on the number of small RNA read mapping positions. We applied a serial nomenclature to the predicted novel miRNAs in C. porosus. For known miRNAs, the miRNA family names from the avian datasets as well as other species were retained with the addition of the prefix "cPor". Two separate heatmaps were generated expressing normalized read counts of the novel and known miRNA candidates using R v3.3 and output from quantifier.pl script. Separate normalized expression profile for several distinct miRNA candidates were generated due to their unique expression profile in a certain tissue type (library) or a group of tissues (libraries). All novel miRNAs were deposited in miRBase.
miRNA target identification was performed by using two independent target prediction tools -miRanda (version 3.3a) and RNAhybrid (version 2.1.2). For miRanda, the option "-strict" was enforced so that the miRNA seed region had to show complete complementarity to a 3′UTR region. Naturally, this option prevented detection of target sites containing gaps or non-canonical base pairing. A free energy of -20 kcal/mol and an alignment threshold score of 170 were chosen as cutoffs in option parameters for miRanda. These stricter values were used to reduce the possibility of false positives when compared to previous analyses [50][51][52][53][54][55]. A FASTA file containing all the 3′ UTR regions of C. porosus genes were provided as the subject for miRNA queries against the newly annotated salt water crocodile genome [21]. Bedtools' "getfasta" command was used to identify the 3′UTR regions of genes and pipe it into the miRanda command line script.
For RNAhybrid, RNAcalibrate (a part of the RNAhybrid package) was run on the mature miRNA sequence FASTA file and the same 3′UTR sequence file as above. The free energy cut off was kept identical to miRanda (-20 kcal/mol) and the p value was set at ≤ 0.01. In an effort to reduce false positives for RNAcalibrate and RNAhybrid, the seed region folding option of the mature miRNA (-f 2, 8) was again enforced.
The TE profile in C. porosus was determined using Repeatmasker (version 4.0.6) [56]. The '-species Crocodylus' option was invoked. rmblastn version 2.2.27+ and RM database and RepBase version 20150807 was used. The coordinates of the mature miRNA sequences, corresponding star and hairpin precursor sequences were all obtained from the miRDeep2 run. The Bedtools "intersect" package was used to identify overlaps with coordinates generated by RepeatMasker.
We used MapMi (version 1.5.9-b31) to determine phylogenetic distributions of the identified miRNAs relative to other vertebrates. Mature miRNAs were mapped to the American alligator (A. mississippiensis; GCA_000281125.4), the Indian gharial (G. gangeticus; JRWT00000000), human (Homo sapiens, GCA_000001405.25), platypus (Ornithorhynchus anatinus; GCF_000002275.2), African clawed frog (Xenopus tropicalis; GCA_000004195.1), chicken (Gallus gallus; GCA_000002315.3), Anna's Hummingbird (Calypte anna; GCA_000699085.1), zebra finch (Taeniopygia guttata; ABQF01000000), green anole (Anolis carolinensis; GCA_000090745.1), painted turtle (Chrysemys picta; AHGY00000000.2) and the saltwater crocodile (C. porosus), by allowing 1 mismatch between the miRNAs and each genome. The genome assemblies were accessed from NCBI Genbank and www.ensembl.org (Ensembl Release 90). MapMi results were processed to retain only unique miRNAs in each species and numbers of shared and lineage-specific loci were plotted on a tree depicting established relationships [57]. Because the core algorithm of MapMi differs from that of miRDeep2, the total number of miRNAs mapped was reduced compared to the 644 reported with miRDeep2. For example, the MapMi algorithm assigns a score based on several biogenesis criteria to a predicted miRNA in a species. Based on matches with miRBase, the default cutoff for the MapMi score is 35 and removal of all miRNAs below this score is recommended by the MapMi program. These factors resulted in 29 fewer miRNAs being mapped through MapMi.
644 unique miRNAs, along with their corresponding star and hairpin precursor sequences were identified. A majority of the small RNA reads fell in the size class of 21-nt and 22-nt. In animals, miRNAs predominate at these two size classes [58][59][60]. The total number of collapsed raw reads for each size class of small RNA reads from all 16 libraries combined is shown in Fig. 2. Although all small RNA size classes do not qualify for miRNA, our dataset of 644 miRNAs had size classes ranging between 18 through 22-nt. A 5′-T (or U) is a hallmark of eukaryotic miRNA biogenesis [61] and this bias was observed in the C. porosus data (Supplemental Fig. 1, Supplementary Material online).
As expected, the mature miRNAs, being the functional and stable component, were associated with higher read counts compared to the star sequences [62,63]. Out of the 644 miRNAs identified, 139 were known miRNAs from birds or other species as catalogued in miRbase (ver 21). Many highly conserved miRNAs like the let-7 family, miR-10, miR-1, miR-30, miR-7, miR-16 were identified in C. porosus suggesting that the additional, novel miRNAs are likely valid. A representation of some known and novel miRNAs in crocodilians and their hairpin precursors are shown as identified by miRDeep2 in this study (Fig. 3). A list of all mature miRNAs, both known and novel, in crocodilians identified in this study is provided in Supplemental file S1, Supplemental Material online. Several miRNAs with varying numbers of loci of origin were found. One deeply conserved miRNA (let7a/j -5p), had four distinct loci, while several others had 3 loci (cPor-miR-N143, cPor-miR-133-5p, cPor-miR-135a/b-5p, cPor-miR-199-5p to name a few).
miRNAs are often expressed in a tissue specific manner. This can be attributed to signaling, compartmentalization, or pattern formation during development. It can also occur in response to localized stress or defense mechanisms and other factors [64][65][66][67][68]. The expression profile of the known and novel miRNAs in C. porosus varied, with the widely conserved miRNAs (and some novel crocodilian miRNAs) displaying very high expression in all tissue types and some novel crocodilian miRNAs peaking only in certain tissue types. There were also several miRNA candidates with very low expression in all tissue types. Despite this low expression, these candidate miRNAs passed all criteria of miRNA biogenesis, possessed all hallmarks of a true miRNA (as processed through the miRDeep2 pipeline) and were hence retained in the final dataset. Heatmaps were generated representing normalized expression profile for the top 20 known and novel crocodilian miRNAs as found in the different small RNA libraries in the saltwater crocodile (Fig. 4a andb).
We identified miRNA targets by using the most recent annotated C. porosus genome as reference [21] and enforcing strict seed complementarity at the 3' end in miRanda [52]. Previous studies indicate that this improves the binding score and provides more efficient target prediction [69][70][71]. To be conservative, we used miRanda with an alignment match score threshold of 170 and free energy threshold of -20 kcal/mol. One or multiple targets were identified for 583 of the 644 unique miRNAs. We attribute the reduction in number of identified miRNAs with targets from the total number of identified miRNAs to the steep alignment threshold score cutoff of 170, which is high compared to earlier studies that used a threshold score of 155-160 [52,53]. When the cutoff score was lowered to 160 (still high considering previous publications) while keeping the free energy value the same, miRanda identified targets for 640 out of the 644 miRNAs. We also performed miRNA target identification using RNAhybrid [72][73][74][75][76].
A majority of the TE-associated miRNA precursors are associated with Class II TEs, aka DNA transposons. A total of 283 miRNA-TE overlaps were found, of which nearly 78.3% were DNA transposons (subclass I + subclass II) overlaps, while 21.6% of the overlaps arose from Class I TEs, i.e. retrotransposons. The remaining 0.1% were derived from unidentified TEs.
Many of the TE-associated miRNA precursors appear to arise from these DNA transposons as they had the same orientation as the TE. Others were complimentary to the TE, possibly indicating the TE to be a target of miRNA-mediated regulation. However, this hypothesis is not testable using data for this study. Not all TE-derived miRNA precursors were fully overlapped by their respective TEs. Therefore, additional analyses were performed to identify four separate categories of miRNAprecursor and TE overlap situations. We expected the most prevalent category of these four types would be when the entire miRNA precursor (including the mature and star sequences) would overlap the TE. This was indeed the case, as ~92.5% of the miRNA-TE hits were found to be complete overlaps of the precursor and TE (Fig. 5). A list of all miRNA precursors that completely overlap with a TE is provided in Supplementary file S2, Supplemental Materials Online.
To investigate the taxonomic distributions of C. porosus miRNAs, we selected ten species representing all four classes of vertebrates. miRNAs were mapped to the genomes of each species with MapMi [77]. Given established relationships [57], Fig. 6 shows the relationship of miRNA presence/absence within each of the selected vertebrates.
This work represents the first de novo analysis of miRNAs in the saltwater crocodile. Small RNAs were isolated, sequenced and analyzed from eleven different tissues in as many as four different individuals per tissue. A total of 16 different small RNA libraries were prepared. More than 187 million reads were generated with each library represented by an average of ~11.7 million reads. Stringent quality control measures were taken to reduce low quality reads, and only those, where more than 50% of the reads scored higher than Phred 20 were accepted for downstream analysis. After filtering, the number of high-quality reads per library varied from just under a quarter million to more than 2.25 million.
For all 644 distinct miRNAs identified in this study, miRDeep2 identified mature and star sequences as well as miRNA precursors and precursor coordinates. A score calculated based on the miRdeep2 algorithm was assigned to each miRNA candidate. The higher the score, usually the better the predicted miRNA candidate. Examples of a mature miRNA and star sequence represented on their hairpin precursor as generated by miRDeep2 are shown in Supplemental Fig. 2, Supplementary Material online. The terms "gga" and "tgu" refer to miRNAs being already identified and present in chicken (Gallus gallus) and zebra finch (Taeniopygia guttata), respectively. A list of all mature miRNA consensus sequences, their corresponding consensus star sequences and precursor sequences generated by the miRDeep2 program is provided in Supplementary file S3, Supplementary Material online.
Several tissue specific miRNA expression profiles could be identified and corroborated from previous studies. One such example is miR-21, which is known to target several apoptotic genes in human spleen [12]. Given that the spleen is the site for expression of various apoptotic genes like Bcl-2, Bak-1, and Caspase-3, it was expected that the crocodilian spleen would also show high expression of miR-21 in C. porosus and this was observed as well (Fig. 4a). In murine tissues (Mus musculus), miR-21 is detected at much higher levels in liver, lung, kidney as compared to the brain [78]. Though kidney and lung tissues were absent in our small RNA library data, miR-21 expression data in C. porosus brain and liver tissues exhibited a similar pattern. miR-133 has been associated with mammalian cardiac muscle and regulation of cardiac Fig. 3. Representation of five selected hairpin precursors of three known (A-C) and two novel crocodilian miRNAs (D-E) identified in the saltwater crocodile. For the known miRNAs (A-C), both the mature miRNA and star sequences are represented in red. For the novel crocodilian miRNAs (D-E), the mature miRNA is in red and the star sequence is in purple, while the loop is golden. hypertrophy as well as functioning as an inhibitor of cardiomyocyte proliferation [79]. Both miR-1 and miR-133 are bicistronic in the mammalian genome and can be thus transcribed together. They are encoded by loci that are duplicated bicistronic and hence they possess identical miRNA mature sequence. Both of these miRNAs are expressed in the heart [80]. We found that both miR-1 and miR-133 achieved peak expression in heart tissue in C. porosus.
Experimental validation of these numerous miRNAs was beyond the scope of this work. However, certain miRNAs like miR-1, miR-20, miR-21, miR-27, miR-30, miR-122, miR-126, miR-132, miR-199, miR-214, miR-338 and miR-363 identified in this study were found in the adult chicken [81] and chicken embryo [82,83] with small RNA blot validation. As expected, the patterns of tissue specific expression of these miRNAs are similar in both chicken and the saltwater crocodile. While some of these miRNAs are highly conserved, some are not and that certainly points to an interesting aspect of miRNA expression and evolution in these two relatively closely related species. Our expression profiles clearly demonstrate highly tissue specific as well as conserved expression patterns (Fig. 7) suggesting that our analysis accurately identifies expression profiles for miRNAs.
With this assurance, we investigated the expression patterns of selected novel crocodilian miRNAs that exhibited tissue-specific patterns. cPor-miR-N1 peaked in expression in the gastro-intestinal regions (small intestine and stomach tissues) of the crocodile while two of the other novel candidates -cPor-miR-N16 and cPor-miR-N17 peaked in the spleen tissue. Some others like cPor-miR-N170 had highest expression in the tongue. Interestingly, cPor-miR-N341 -one of the novel miRNA candidates, had expression in all testes tissues, but it's expression was essentially zero to very low in all other tissues. Specific functions of these miRNAs are unknown but these patterns suggest regulatory processes related to these tissues. A comprehensive list of the normalized expression profiles of all identified miRNAs in this study as determined by the quantifier.pl script of the miRDeep2 program has been provided in Supplemental file S4, Supplemental Material online. Lastly, although not the focus of this work-it was observed that the total normalized expression of miRNAs in the heart tissue type was higher in the normal hatchlings as compared to the runts. The opposite was true for the miRNA expression pattern found in the brain between the two types of hatchlings. In future work, it will be interesting to follow up on a comparative miRNA expression profiles of normal and runt hatchlings.
Although miRNA targeting in animals primarily takes place in the 3′ UTR region of targeted genes, reports have shown the 5′ UTR of genes and open reading frames of some highly repeat rich elements to be miRNA targets as well [84][85][86]. Binding and targeting are most effective if there is complete complementarity between the target with the seed region, the 2-8th position on the 5' end of the mature miRNA [1,[87][88][89][90][91][92]. It has been noted previously that target prediction tools which require Watson-Crick base pairing with the query miRNA seed region, perform better in predicting true targets as compared to those programs that do not have this requirement [93,94].
As mentioned above, we used miRanda and RNAhybrid to predict the miRNA targets. RNAhybrid mainly works based on free energy of folding of miRNA around the target and a calibration step improves overall target estimation for the miRNAs. For RNAhybrid, we used the same value for free energy cutoff as that of miRanda but found that the miRanda results were more comprehensive. Thus, the miRNA-target data shown here only corresponds to those identified by miRanda. The four miRNA candidates that failed to be predicted by miRanda with a lowered complementary match cutoff score of 160 included miR-32-5p, miR-126-5p, miR-N512 and miR-208a. All of these miRNA candidates had zero normalized counts in multiple libraries and extremely low total normalized count overall. However, having satisfied all criteria of being a true miRNA by the miRDeep2 pipeline and possessing all hallmarks of a true miRNA biogenesis, these candidates were retained. It is possible, that with a larger target dataset, these miRNAs might find a suitable target protein.
Because of the significant role miRNAs play in posttranscriptional gene regulation and the increasing focus on TEs as a major player in the same, we investigated the relationship between the two. Our results indicate that nearly half (44%) of the miRNAs we identified overlapped in some way with TE loci and 94% of those arose completely from within the TE. DNA transposons are the mobile element fraction most highly associated with miRNAs. A mixed distribution of both complementary and sense orientated population of miRNA (with respect to the TEs) suggested that miRNAs could target the TEs from which they originated. This suggests a potential mechanism for suppressing the expression of the associated TEs. This is interesting in that TEs have not accumulated at significant rates in crocodilians over the past several million years [18].
Our results are in good agreement with previous studies that demonstrate the role of TEs in the evolution of functional miRNAs [44,95]. Specifically, 62% of miRNAs in the vespertilionid bat Eptesicus fucus have been known to be associated with TEs while that percentage changes to 23.9% and 16.2% in dog and horse respectively [96]. In the same vein, miRNA association has been found to be more abundant with retrotransposons in other mammals, including humans. In addition, miRNAs have been also known to be primarily associated with DNA transposons in the zebrafish and the western clawed frog [97]. This further relates TEs to the evolution of phenotypic and phylogenetic diversity. Several hypotheses exist suggesting that TEs play a significant role in mechanisms of diversification [98][99][100]. There is also some evidence for the temporal correlation between bursts of TE accumulation (and therefore miRNA deposition) with speciation [101][102][103][104][105][106]. Thus, we believe that the TE-derived miRNAs described here are good candidates for further investigation with regards to identifying lineagespecific functional differences in gene expression across Crocodylia.
Finally, when analyzing the origin of C. porosus miRNAs, we observed that as taxonomic distance increases, the number of shared miRNAs decreases. The vast majority of miRNAs recovered from C. porosus were identifiable in other crocodilians. As would be expected given, the more recent divergence of gharials from true crocodiles (~55 mya vs. ~85 mya for the alligator divergence) [107], more miRNAs were shared by G. gangeticus and C. porosus, when compared to the two available alligator genomes (A. mississippiensis and A. sinensis) (Fig. 8).
We identified both known and novel miRNAs in the saltwater crocodile, many of which are orthologous to several other vertebrates. A weakness of this work was our inability to validate miRNAs by identifying associations with RISC proteins or demonstrating definitive impacts on mRNAs and/or phenotype. Future work including target validation through 5' Rapid amplification of cDNA ends (RACE) as well as approaches like degradome assays [108,109] might provide target validation of some of these miRNAs identified here. However, our conservative scoring strategies and the fact that many of the miRNAs we identified were orthologous with other known and validated miRNAs suggests that our pipelines were effective.
The data here provide an intriguing start to analyses of gene regulation in this ecologically important clade that is, nonetheless, rarely studied from a molecular perspective. Additional analyses, especially those that can be compared with the only other extant member of Archosauria, birds, could provide insight into the ancient and extinct relatives of both groups, the ever-interesting dinosaurs as well as into the evolution of gene regulation in non-mammalian vertebrates. For example, similar and/or identical miRNA target sites can reveal potentially conserved target sites among vertebrates. More direct applications arise from an increased understanding of gene regulation as it may relate to the development of commercially important features, especially hide, which is highly valued.