Genome-wide microRNA screening in Nile tilapia reveals pervasive isomiRs' transcription, sex-biased arm switching and increasing complexity of expression throughout development

Abstract

MicroRNAs (miRNAs) are key regulators of gene expression in multicellular organisms. The elucidation of miRNA function and evolution depends on the identification and characterization of miRNA repertoire of strategic organisms, as the fast-evolving cichlid fishes. Using RNA-seq and comparative genomics we carried out an in-depth report of miRNAs in Nile tilapia (Oreochromis niloticus), an emergent model organism to investigate evo-devo mechanisms. Five hundred known miRNAs and almost one hundred putative novel vertebrate miRNAs have been identified, many of which seem to be teleost-specific, cichlid-specific or tilapia-specific. Abundant miRNA isoforms (isomiRs) were identified with modifications in both 5p and 3p miRNA transcripts. Changes in arm usage (arm switching) of nine miRNAs were detected in early development, adult stage and even between male and female samples. We found an increasing complexity of miRNA expression during ontogenetic development, revealing a remarkable synchronism between the rate of new miRNAs recruitment and morphological changes. Overall, our results enlarge vertebrate miRNA collection and reveal a notable differential ratio of miRNA arms and isoforms influenced by sex and developmental life stage, providing a better picture of the evolutionary and spatiotemporal dynamics of miRNAs.

Figure 1

RNA-seq overview. (A) Reads obtained from RNA-seq. (B) Size distribution of “size selected” reads used for the analysis. (C) Read count of canonical and isomiR reads.

Figure 2

Patterns of miRNAs distribution. Known (A) and Novel (B) miRNAs detected (>5 reads) among developmental stages and adult tissues samples. *dpf = days post-fertilization; developmental stages, liver and eye correspond to pool samples.

Figure 3

Genome organization of miRNAs genes in Nile tilapia genome. Black: known miRNAs. Red: novel miRNAs described in current study. Green: novel miRNAs described previously by Huang et al. and Yan et al.; MiRNAs in brackets: clustered miRNAs.

Figure 4

Categorization and quantification of isomiRs sequences. Isoforms were categorized into canonical sequences; Shorter – isomiR sequences with reduced length when compared to canonical sequence; Longer – isomiR sequences with increasing nucleotide base; Seed shifted – isomiRs with any modification at 5′ end (nucleotide gain or loss) that changed seed region (2–7 nt region); and Non-template - IsomiRs with any nucleotide mismatch with the precursor sequence. Henceforth canonical sequences were defined as the most miRNA 5p and 3p expressed sequences of combined samples.

Figure 5

Individual profiles of miRNAs showing switches in arm counts at early and adult stages. X-axis represents 3p miRNA mapped raw read counts (log₂) and Y-axis represents 5p miRNA mapped raw read counts (log₂). Diagonal line represents same expression rates from 5p and 3p. Data series on the right side of the line shows greater proportion of 5p over 3p, while data series on the left side of the line shows a greater proportion of 3p over 5p.Threshould (dashed red line) represents Fold change ≥2 or ≤−2.

Figure 6

miRNA arm switching events detected and ranked in Nile tilapia. The miRNA 5p and 3p expression prevalence in each sample is represented by normalized fold change rates. Rates were normalized by the difference between the miRNA fold change from each sample and the mean of miRNA fold change of samples to collapse the series curves aiming at the improvement of data visualization. The 5p and 3p variations were ranked among samples organized from the greatest (prevalent 5p) to the smallest (prevalent 3p) fold rates per each pre-miRNA. The pre-miRNAs were also ranked from the greatest to the smallest distance between max and min fold change among samples (from dark to light blue series).

Figure 7

Entropy analysis of miRNAs. Shannon’s entropy was used to quantify miRNA diversity for developmental stages through adulthood. We observed an increase of diversity at early developmental stages (2 to 5 dpf), where cell differentiation and organ formation occurs, followed by stabilization in late larval (10 dpf) and adult (6 month years old) stages. Developmental stage labels and images were obtained from Fujimura and Okada, 2007 under Creative Commons Attribution 4.0 International (CC BY 4.0 - authors “Koji Fujimura and Norihiro Okada”, title “Fig. 1. Overview of the developmental stages of Nile tilapia O. niloticus.”, date “14 March 2007”, url “http://onlinelibrary.wiley.com/doi/10.1111/j.1440-169X.2007.00926.X/full”, no endorsement, only esthetical adaptations: image vectorized, background removed, shadow included and image rotated. Adult Oreochromis niloticus image from https://www.usgs.gov/ under the Public domain.

Figure 8

Detection of novel Nile tilapia miRNAs in diverse vertebrate species by miRDeep2 approach and ortholog/paralog annotation. (A) Vertebrate evolutionary tree with the number of novel miRNAs (inside circles) detected in each species through genome-wide matching searches by miRDeep2 software. Phylogenetic tree is a handmade tree derived by merging published trees by Brawand et al. and Amemiya et al.. Homo sapiens, Ornithorhynchus anatinus, Petromyzon marinus, Danio rerio, Salmo salar, Gallus gallus and Monodelphis domestica images were obtained from https://commons.wikimedia.org/wiki/Main_Page under the Public Domain; Macaca mulata and Mus musculus images from https://pixabay.com under the Public Domain; Pundamilia nyererei, Haplochromis burtoni and Neolamprologus brichardi adapted by permission from Macmillan Publishers Ltd: Nature (Brawand et al.), copyright 2014, https://www.nature.com/; Xenopus tropicalis and Xenopus laevis images from http://journals.plos.org/plosbiology/ under Creative Commons Attribution 4.0 International (CC BY 4.0 - authors “Patrick Narbonne, David E. Simpson and John B. Gurdon”, title “Figure S1. Characterization of the hybrids formed by the cross-fertilization of X. laevis eggs with X. tropicalis sperm.”, date “15 November, 2011”, url “https://doi.org/10.1371/journal.pbio.1001197.s001”, no endorsement, only esthetical adaptations: image vectorized, background removed, shadow included and image rotated); Metriaclima zebra or Maylandia zebra and Anolis carolinensis images from http://www.publicdomainpictures.net under CC0 1.0 Universal: CC0 1.0 Public Domain Dedication; Oreochromis niloticus image from https://www.usgs.gov/ under the Public domain. (B) Comparative features of 42 novel miRNAs among species, including sense orientation, paralogs abundance, and arms detection. Other 55 novel miRNAs were tilapia-specific, and therefore, excluded from this comparative analysis. (C) Orthologs and paralogs relationship cardinalities considering common anscestral and the “young” gene copies from novel Nile tilapia miRNA genes. Orthologs were infered from reconciled gene and species trees obtained with RAxML and TreeBest using predicted data and homologs miRNA data from miRBase (see methodology).

Figure 9

Tracing the evolution of oni-miRNA-n780. (A) Syntheny block shows a high conservation among teleost genes. The origin of oni-miRNA-n780 paralogs is associated with extra rounds of genome duplication between in teleosts, since the lambrey presents only one copy of the miRNA cluster 29d/a and cichlids show three copies. In amphibians, reptiles and mammals, oni-miR-n780 paralog could have been lost, while the cluster 29b/a present in mammals or 29d/a present in amphibians might have emerged from a same ancestor paralogs found in lampreys. Syntheny blocks were obtained from Genomicus and altered with homology data from Ensembl about oni-mir-n780 paralogs from other vertebrate species. (B) A manual check using the gene annotation in cichlids and the homology level among genes shows a high conservation of the synthetic block among the cichlid species. However, among M. zebra and N. brichardi there was an inversion of the genomic block. Each color represents a different gene and its transparency level represents the homology level (high homology/low transparency to low homology/high transparency). The oni-mir-n780 is represented by the red symbol in the center of the image. (C) A reconciled gene and species tree of oni-miR-n780 using RAxML and TreeBest to reconstructing the historical evolutive events among species. We found that the mir-29b from Apis mellifera is the closest common miRNA ancestor of oni-miR-n780. (D) Alignment of novel mir-n780 precursor sequences with homologs found in lamprey, amphibians and other teleost fish. Sequence similarity of oni-mir-n780 within all species is around 68.8%, although inside cichlid clade it has 100% of identity. The mature miRNA 3p is highly conserved, showing 95.2% of similarity among vertebrats and 100% of similarity when compared among cichlid fishes, B. belcheri, X. laevis, X. tropicalis and S. salar; by contrast the mature miRNA 5p form is more variable presenting 80.8% of similarity within vertebrates, but still highly conserved in the cichlid genomes where they show of 100% of identity. Another features include a G-to-A mutation in the seed region of the mature miRNA 5p sequences exclusively detected in cichlids, which might refer to a possible specific functional regulatory role of oni-miR-n780 in the genome of cichlids. Alignment data was retrieved from Geneious.