Interplay Between MicroRNAs and Targeted Genes in Cellular Homeostasis of Adult Zebrafish (Danio rerio).

BACKGROUND: Cellular homeostasis is regulated by the intricate interplay between a plethora of signaling pathways and "energetic sensors" in organs. In order to maintain energy balance, induction or repression of metabolic pathways must be regulated and act in concert with the energetic demands of the cell at a given point in time. A new class of small noncoding RNAs, the microRNAs (miRNAs), has added yet further complexity to the control of metabolic homeostasis.
OBJECTIVE: Understanding the damages induced by toxins in the liver and the intestine as well as the interplay between the miRNome and transcriptome first requires baseline characterization in these tissues in healthy animals under cellular homeostasis.
METHODS: The liver (main site for detoxification) and the gut (primary exposure routes for contaminant exposure) were dissected out (wildtype fish), total and small RNA extracted, mRNA and miRNA libraries constructed and subjected to high throughput sequencing. Differential Expression (DE) analysis was performed comparing liver with gut and an "miRNA matrix" that integrates the miRNA-seq and mRNA-seq data was constructed.
RESULTS: Both the miRNome and transcriptome of the liver and gut tissues were characterized and putative novel miRNAs were identified. Exploration of the "miRNA matrix" regulatory network revealed that miRNAs uniquely expressed in the liver or gut tissue regulated fundamental cellular processes important for both organs, and that commonly expressed miRNAs in both tissues regulated biological processes that were specific to either the liver or the gut.
CONCLUSION: The result of our analyses revealed new insights into microRNA function in these tissues.

INTRODUCTION

Cell homeostasis is regulated by the intricate interplay between a plethora of signaling pathways and “energetic sensors” in tissues and organs. In order to maintain energy balance, induction or repression of metabolic pathways must be regulated so they act in concert with the energetic demands of the cell at a given point in time. A series of well-established signaling cascades, and conversely, the cells metabolic profile guides signal transduction events by incorporating metabolic sensors into the wide range of signaling molecules. A new group of small noncoding RNAs, the microRNAs (miRNAs), has added further complexity to the control of metabolic homeostasis [1]. MiRNAs are small non-coding RNAs that play important roles in gene regulation in animals and plants by pairing, with semi- or full-complementarity, to messenger RNAs (mRNAs) of protein-coding genes to guide their post-transcriptional repression; miRNAs bind to the 3’ untranslated region (UTR) of target mRNAs in combination with the RNA-inducing silencing complex (RISC) and stimulate either mRNA degradation or translational inhibition [2]. Computational models suggest that each miRNA can affect many genes for silencing [3], and a given gene can be the target of more than one miRNA, thereby generating highly intricate miRNAs networks with pervasive downstream regulatory effects. MiRNA networks regulate protein-coding gene expression and thereby a wide spectrum of biological functions fundamental for disease progression and adverse health outcomes [4-13].

The zebrafish (Danio rerio) is an excellent model for toxicology and disease research. The majority of human protein-coding human genes (70%) have a counterpart in zebrafish and this increases when just human disease genes are considered (84%) [14]. Its high fertility, relatively low maintenance costs, and ease of genome editing make it an appealing alternative to rodent models. Additionally, miRNA functions can be explored in zebrafish with a variety of transient assays that may be more cumbersome to implement in mice, including the insertion of double-stranded RNA molecules (dsRNAs), antisense oligonucleotides and fluorescently labeled reporters into fertilized eggs [15].

Compared to other fish models, miRNAs are best characterized in zebrafish [16]. Approximately 415 miRNAs specific to Danio rerio (dre-miRNAs) have been identified [17] and their role in early development and disease has been characterized (Table 4.2 in Freeman et al.’s review [16]). However, this field of research is new and the functions of several dre-miRNAs, as well as their tissue-specific expression profile, in the context of studies on endocrine disruption remain to be defined.

In this study, we focused on two important zebrafish tissues, the liver and the intestine. The gut represents one of the primary exposure routes for contaminants and the liver is the main site for detoxification [18]. Understanding how diseases or exposure to contaminants affects the interplay between the miRNome and transcriptome in the liver and intestine first requires baseline characterization in these tissues in healthy animals under cellular homeostasis. In order to achieve this objective, the liver and gut were dissected (wildtype fish), total and small RNA were harvested, mRNA and miRNA libraries generated and subjected to Illumina High Throughput Sequencing (HTS). We carried out Differential Expression (DE) analysis and compared liver miRNA expression to intestine using established bioinformatics pipelines [19-21]. Using this approach, known and putative novel miRNAs were identified. Finally, we constructed a “miRNA matrix” that integrates the mRNA-seq and miRNA-seq datasets and facilitates the exploration of the regulatory networks between the miRNome and the transcriptome and the impact they have on biological processes in these tissues.

MATERIALS AND METHODS

Zebrafish Care

Male zebrafish were maintained at 26-29°C with a 14:10 h light-dark cycle. Water pH ranged from 7.0 to 7.6. Sponge filters provided aeration and filtration. Fish were fed twice daily with flaked food (Tetra, Melle, Germany). Fish were acclimated prior to extracting the tissues, i.e. liver and intestine for seven days. Procedures strictly followed The University of California San Diego, IACUC guidelines, AUP S09418. Fish were treated and sacrificed humanely. Liver and gut tissues were dissected out and instantly frozen in liquid nitrogen and then stored at -70°C.

RNA Extraction and Sequencing Library Preparation

RNA extraction from zebrafish liver and gut samples was performed using TRIzol reagent (Invitrogen/Thermo Fisher Scientific, Waltham, MA, USA) and further purification of the RNA was carried out with RNeasy Mini (Qiagen, Valencia, CA, USA). All RNA was subjected to on-column digestion of DNA to yield pure RNA (free from DNA contamination). The concentrations were determined by absorbance readings (OD) at 260 nm using an ND-1000 (Nanodrop, Wilmington, DE, USA). RNA was tested for structural integrity with the 6000 Nano LabChip assay from Agilent, (Santa Clara, CA, USA). Only RNA samples with RIN scores >7.0 were used for RNA-seq. For the zebrafish RNA-seq experiments, 10 individual fish tissues were pooled together. The TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA) was utilized to construct mRNA-seq libraries with 100-200 ng of total input RNA. The miRNA-seq libraries were constructed using the TruSeq Small RNA Prep kit and 1 μg input RNA.

High Throughput Sequencing (HTS)

HTS was performed using an Illumina HiSeq2000. Sequencing depth for RNA-seq was ~20 to 30 million reads per library and ~5 million reads for miRNA-seq. A single end 50 cycle sequencing approach was used. Data were subjected to standard Illumina best practices (>80% = Q30). The miRNA-seq and RNA-seq datasets have been submitted to the NCBI Gene Expression Omnibus, with accession designation number GSE108437 providing access to all of the data in this manuscript. By convention Danio rerio, miRNAs are labeled “dre-miRNAs”, i.e., dre-miR-133a.

Gene Level Analyses

Secondary analyses of the RNA-seq and miRNA-seq data were carried out on an OnRamp Bioinformatics Genomics Research Platform as previously described by us (OnRamp Bioinformatics, San Diego, CA, USA) [19]. The Analysis Engine utilizes an automated RNA-seq workflow to process data (Fig. , blue), including (1) FastQC [19]; (2) CutAdapt [22]; (3) STAR [23]; (4) HTSeq [24]; and (5) DESeq2 [25]. The Platform also provides the Tuxedo Suite [26] that includes TopHat2 [27-29], Cufflinks/Cuffmerge to assemble transcripts and merge assemblies, and Cuffdiff to perform differential expression analysis that delivers Fragments Per Kilobase of transcript per Million mapped reads (FPKM). Transcript count data from DE analyses of the liver and intestine samples were ranked according to their adjusted p-value (or q), which is the lowest False Discovery Rate (FDR), at which a transcript is significant. FDR was determined using the Benjamini-Hochberg multiple testing with an FDR cutoff at q ≤ 0.01 on the RNA-seq dataset.

The Comprehensive Analysis Pipeline for miRNA sequencing data (CAP-miR-Seq) was used for pre-processing, alignment, mature/precursor/novel miRNA detection and quantification, and data presentation (Fig. , green) [30]. The miRNA-seq data was aligned to Zv9 zebrafish genome using miRDeep [30, 31]. DE analysis of liver and intestine miRNAs was performed with EdgeR/miRDiff [32, 33] with an FDR cutoff of q ≤ 0.01.

Heatmaps of DE transcripts were generated using the heatmap.2 from gplots R-package [34]. Venn diagrams were generated using VENNY 2.1 online tool [35].

System Level Analyses

DE zebrafish transcripts were further analyzed with (1) Gene Ontology enRIchment anaLysis and visuaLizAtion

(GOrilla) [36, 37] and (2) REduce & VIsualize Gene Ontology (REViGO) [38].

We also utilized the orthologies available from Ensembl to attach a human gene ID to a given zebrafish one [39]. This “human centric” dataset was analyzed using (1) ToppFun [40], and (2) iPathwayGuide by Advaita Bioinformatics [41], which utilize the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [42], GO terms [43], miRNAs from miRBase [44] and TARGETSCAN [45], and diseases from KEGG.

Construction of miRNAs-mRNAs Matrix

Owing to the fact that an individual miRNA can elicit pervasive effects by targeting multiple transcripts for silencing and that any one transcript can be targeted by multiple miRNAs, our goal was to construct a “miRNAs-mRNAs matrix” to fully comprehend the impact that the identified miRNAs have on the respective liver and intestinal transcriptomes, and to characterize the pathways and processes regulated by this miRNA-mRNA matrix. First, for each miRNA expressed in the liver and the gut tissues (count ≥ 10) [44], predicted targets were identified with TargetScanFish 6.2 (using a context score ≥ -0.05) [46-48]. Next, we generated an “adjacency matrix” of miRNAs and predicted target genes present in the RNA-seq dataset from the Tuxedo Suite (conservative cutoff set at FPKM ≥ 5) [49].

From this matrix, several output figures are generated: 1) miRNA impact table with the sum of the predicted genes found within the list of DE genes (2nd column: sum of the predicted target genes that each miRNA of interest has in the FPKM ≥ 5 gene list), as well as the percentage of all predicted targets (3rd column: (sum / total # of target genes all miRNAs of interest have in the FPKM ≥ 5 gene list) x 100) and the percentage of all DE genes (4th column: (Sum / all DE genes in the FPKM ≥ 5 list, including genes that are not targeted by miRNAs of interest) x 100) that this sum represents, 2) a barplot illustrating percentages from the 3rd and 4th columns, 3) a dendrogram showing the clustering of miRNAs determined based on the number of common genes they target for silencing, and 4) a sample identity heatmap.

RESULTS

Overview of mRNA- and miRNA-sequencing

In this study, two mRNA and two miRNA libraries from the liver and the gut were sequenced to obtain the complete set of mRNAs and miRNAs transcribed under cellular homeostasis (i.e. in the absence of exposure to toxicants) in these two tissues. The depth of sequencing averaged to 2,743,528 (± 257,530) and 5,734,461 (± 534,447) reads in the liver and gut small RNA libraries respectively. Approximately the same percentage of miRNA was identified in both tissues as shown in the pie charts (Supplementary Fig. ): 93% and 97% of the liver and gut RNA quantifications are respectively labeled as miRNA. The depth of sequencing for the total RNA libraries averaged to 25,245,624 (±5,125,531) and 28,651,653 (±2,951,565) reads for the liver and gut respectively. The length distribution based on total abundance was very similar in the small RNA libraries for liver1/liver2 and gut1/gut2, ranging from 21 to 24 nucleotides (nt) (Supplementary Fig. ). Our data show that the most dominant sequence reads were 22 nt small RNAs, which is consistent with small RNA lengths in animals [50].

Characterization of miRNome & Transcriptome in Liver and Gut Tissues

In order to investigate previously identified and uncover novel miRNAs in the zebrafish homeostatic liver and gut tissues, the mappable sequences obtained from HTS were aligned to zebrafish miRNA sequences (miRDeep, genome Zv9). In deep sequencing experiments, the read counts of miRNA in libraries can be used as an index to estimate their relative expression abundance [50]. In the liver, 224 known miRNAs were identified (count ≥ 10). Under the same criteria, 233 miRNAs were identified in the gut tissue. When comparing liver to gut miRNomes, we determined that 215 miRNAs were commonly expressed in both tissues (Fig. ), and that 9 and 18 miRNAs were uniquely expressed in the liver and the gut, respectively.

Based on the Tuxedo Suite output, we obtained FPKM values for each gene in the RNA-seq dataset; FPKM reflects gene relative expression abundance. We identified 5,988 and 5,817 genes in the liver and gut transcriptomes, respectively. When comparing liver to gut transcriptomes, we determined that 4,438 genes were commonly expressed in both tissues, while 1,550 and 1,379 genes were uniquely expressed in the liver and the gut, respectively (Fig. ).

ToppFun Functional Enrichment Analysis of the Transcriptomes

We utilized Ensembl homology relationships to attach a human gene ID to a given zebrafish one in order to permit a systems level analysis. We exploited the “Transcriptome, Ontology, Phenotype, Proteome and Pharmacome annotations based gene list FUNctional enrichment analysis” (ToppFun) tool and better-characterized Gene Ontology (GO) terms available in human compared to zebrafish. We entered 3 gene lists in ToppFun called “common”, “unique liver” and “unique gut” based on the comparison performed previously (Fig. ). ToppFun provided output files with GO terms for Biological Process (BP), Molecular Function (MF) and Cellular Component (CC); from these files, the GO term IDs and p-values were entered in REViGO, which simplifies GO term analysis via a clustering algorithm based on semantic similarity measures (Fig. ; Supplementary Fig. ). ToppFun also provided a list of enriched Pathways (PW) for common, unique to the liver and unique to gut datasets (Supplementary Tables ).

In the gene set common to liver and gut tissues, the most enriched BP are related to cellular metabolism including cellular catabolism, cellular amide/peptide/organic acid metabolism and oxidation-reduction process (Fig. ). Other enriched BP terms include organonitrogen compound biosynthesis, translational initiation, protein localization to endoplasmic reticulum, interspecies interaction between organisms and ribonucleoprotein complex biogenesis. RNA binding, oxidoreductase activity, structural constituent of ribosome and cofactor binding are amongst the most enriched MF (Supplementary Fig. ) followed by NADH dehydrogenase (ubiquitone) activity, oxidoreductase activity, acting on NAD(P)H, coenzyme binding and electron carrier activity. The most pertinent CC terms are related to mitochondria (Supplementary Fig. ; mitochondrial part, mitochondrion and respiratory chain), to the envelope (organelle envelope), and to the endoplasmic reticulum (ribosome, cytosolic ribosome and the ribonucleoprotein complex).

In the gene set unique to the liver, BP terms associated with RNA processing and metabolism are most enriched (tRNA and ncRNA metabolism, cleavage involved in rRNA processing, nucleic acid phosphodiester bond hydrolysis) as well as mitochondrion organization (Fig. ). SNAP receptor activity and methyltransferase activity are the two most enriched MFs (Supplementary Fig. ). Mitochondrial part and mitochondrion dominate the CC terms, followed by SNARE and envelope related terms (Supplementary Fig. ).

In the gene set unique to gut tissue, terms related to the metabolic process are the most enriched BP terms, including lipid/sphingolipid/cellular lipid/carbohydrate metabolism, organic substance and regulation of catabolism (Fig. ). Other significant BP terms are related to the cellular component organization (positive regulation of cellular component and endomembrane system organization) and to the immune system process (antigen processing and presentation, defense response). Protein binding (enzyme, lipid, cytoskeletal, macromolecular complex and ATP binding) and transporter activity (symporter, amide/organic anion transmembrane transporter activity) represent the most enriched MF in the dataset unique to the gut (Supplementary Fig. ). In terms of CC, the most significant components are related to intestine-specific cell projections, such as brush border, cluster of actin-based cell projections (microvillus) and cell projection membrane, followed by membrane region (plasma membrane region), vacuole (lytic vacuole) and apical part of the cell (Supplementary Fig. ).

Network of miRNAs and Target Genes

To take advantage of the fact that in this study, both RNA-seq and miRNA-seq were performed on the same sample, we decided to build the networks of miRNAs and their predicted target genes in order to characterize the impact that miRNAs have on the liver and gut transcriptomes. In the initial step of our analysis, predicted target genes for each miRNA of interest were identified using TargetScanFish (genes with a context score ≥ -0.05 were selected). Once these “miRNA-target genes” lists were established, we then determined whether the predicted target genes listed were present in the liver and gut transcriptomes (FPKM ≥ 5) that we obtained by RNA-seq. We decided to create several networks based on 2 scenarios described in Fig. (.

Fig. (4)

Scenarios to analyze networks of miRNAs and their target genes. A) In the 1st scenario, we analyzed the impact of the 9 miRNAs uniquely expressed in the liver onto the liver transcriptome (5988 genes, blue box) and the impact of the 18 miRNAs uniquely expressed in the gut onto the gut transcriptome (5817 genes, red box). B) In the 2nd scenario, we created a network between the 215 miRNAs commonly expressed in both liver and gut tissues and the 1,550 genes uniquely expressed in the liver (green box) as well as another network between the 215 miRNAs commonly expressed and the 1,379 genes uniquely expressed in the gut (purple box). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)

The goal of the 1st scenario was to determine the regulatory function of miRNAs that are uniquely expressed in an organ onto the corresponding transcriptome. Therefore, we analyzed the impact of the 9 miRNAs uniquely expressed in the liver on the liver transcriptome (5,988 genes) and the impact of the 18 miRNAs uniquely expressed in the gut on the gut transcriptome (5,817 genes) (Fig. ). Two networks were generated and a GO analysis was then performed on each network using ToppFun and REViGO.

The goal of the 2nd scenario was to examine how the environment in which a miRNA is expressed affects its function by taking commonly expressed miRNAs and analyzing their specific function in the liver or the gut by only considering uniquely expressed genes in these organs (Fig. ). Therefore, we created a network between the 215 miRNAs commonly expressed in both liver and gut tissues and 1,550 genes uniquely expressed in the liver as well as another network between the 215 miRNAs commonly expressed and the 1,379 genes uniquely expressed in the gut, and performed GO analyses on both networks.

In the 1st scenario, we determined that 1) only 5 miRNAs out of 9 uniquely expressed in the liver had predicted targets in the liver transcriptome (Fig. and 5B), 2) only 2,052 genes out of 5,988 genes in the liver transcriptome were predicted targets for these 5 miRNAs, 3) dre-miR-144-3p had the most predicted target genes in the liver transcriptome (1,065 genes) followed by dre-miR-460-5p (793 genes), -730 (534 genes), -190b (420 genes) and -733 (101 genes), and 4) dre-miR-144-3p and -460-5p target similar genes as shown in the dendrogram (Fig. ). In the gut, we determined that 1) only 7 miRNAs out of 18 uniquely expressed had predicted targets in the gut transcriptome (Fig. and 6B), 2) only 2,062 genes out of 5,817 genes in the gut transcriptome were predicted targets for these 7 miRNAs, 3) dre-miR-137-3p had the most predicted target genes in the gut transcriptome (676 genes) followed by dre-miR-429b (673 genes), -1788-5p (603 genes), -736 (481 genes), -7a (402 genes), -196c (221 genes) and -196d (221 genes), 4) dre-miR-137-3p and -1788-5p target similar genes as shown in the dendrogram (Fig. ), and 5) dre-miR-196c and -196d share the same predicted targets.

By projecting the identified zebrafish predicted target genes from the liver and the gut transcriptomes onto their human orthologs, we were able to perform a GO analysis using ToppFun (Fig. ). Based on the target genes of the 5 miRNAs expressed only in the liver, the top most enriched GO terms were 1) top 7 BP: cellular catabolism, oxidation-reduction process, cofactor/peptide/organic acid/cellular amide metabolism and organonitrogen compound biosynthesis (Fig. ), 2) top 5 MF: cofactor and RNA binding, oxidoreductase activity, coenzyme binding and oxidoreductase activity, acting on CH-OH group of donors (Supplementary Fig. ), and 3) top 5 CC: mitochondrion, mitochondrial part, endoplasmic reticulum membrane, envelope and organelle envelope (Supplementary Fig. ). Based on the target genes of the 7 miRNAs expressed only in the gut, the top most enriched GO terms were 1) top 7 BP: cellular catabolism, organic acid/ lipid/cellular lipid/monocarboxylic acid metabolism, oxidation-reduction process and regulation of catabolism (Fig. ), 2) top 6 MF: RNA/coenzyme/cofactor/enzyme binding, oxidoreductase activity and oxidoreductase activity, acting on CH-OH group of donors (Supplementary Fig. ), and 3) top 6 CC: mitochondrion, mitochondrial part, endoplasmic reticulum, organelle envelope, envelope and vacuole (Supplementary Fig. ). Interestingly, many of the GO terms identified here based on the unique matrix of miRNAs for the liver and the gut are similar, suggesting that these biological processes, molecular functions and cellular components are regulated by different miRNAs.

In the 2nd scenario, we determined that 1) only 164 miRNAs out of 215 commonly expressed in the liver and the gut had predicted targets in the unique liver transcriptome, 2) only 1,009 genes out of 1,550 genes uniquely expressed in the liver were predicted targets for these 164 miRNAs, and 3) dre-miR-2184 had the most predicted target genes in the unique liver transcriptome (290 genes) (Table ). In the gut, we determined that 1) only 164 miRNAs out of 215 commonly expressed in the liver and the gut had predicted targets in the unique gut transcriptome, 2) only 950 genes out of 1,379 genes uniquely expressed in the gut were predicted targets for these 164 miRNAs, and 3) dre-miR-2184 had the most predicted target genes in the unique gut transcriptome (329 genes) (Table ).

By taking the human orthologs of the identified predicted target genes from the unique liver and the unique gut transcriptomes, we carried out a GO analysis using ToppFun (Fig. ). Based on the target genes of the 164 miRNAs commonly expressed in the liver and the gut tissues, the top most enriched GO terms in the unique liver transcriptome were 1) top 5 BP (Fig. ): tRNA and ncRNA metabolism, mitochondrion organization, RNA processing and RNA (guanine-N7) methylation, 2) MF (Supplementary Fig. ): no MF was identified, and 3) top 9 CC (Supplementary Fig. ): mitochondrion, mitochondrial part and transmembrane space, nucleolus, envelope and organelle envelope lumen, transferase and catalytic complex. Based on the same 164 miRNAs commonly expressed in both tissues, the most enriched GO terms in the unique gut transcriptome were 1) top 7 BP (Fig. ): lipid/sphingolipid/cellular lipid metabolism, endomembrane system organization, positive regulation of cellular component organization, vesicle-mediated transport and regulation of intracellular signal 
transduction, 2) top 5 MF (Supplementary Fig. ): enzyme/protein complex/ATP/macromolecular complex binding, and amide transmembrane transporter activity, and 3) top 6 CC (Supplementary Fig. ): brush border, cluster of actin-based cell projections, membrane region, cell-cell junction, apical part of cell and vacuole. Interestingly, many of the top GO terms identified here based on the matrix of miRNAs are commonly expressed in the liver and the gut which is quite different, suggesting that this pool of miRNAs regulate different biological processes, molecular functions and cellular components depending on which tissue they are expressed from.

Differential Expression Analyses - Liver Versus Gut

The differential expression of previously characterized mature miRNAs in the liver compared to gut tissue, performed with EdgeR/miRdiff (genome Zv9, Supplementary Fig. ), revealed 108 DE miRNAs (q ≤ 0.01 and absolute logFC ≥ 1), amongst which, 29 are more highly expressed and 79 are less expressed in the liver as compared to gut tissue (Supplementary Table ).

The differential expression of total mRNAs in the two tissues was assessed with DESeq2 (genome GRCz10, Supplementary Fig. ) and revealed 6,248 DE zebrafish genes (q ≤ 0.01 and absolute log2FC ≥ 1), amongst which, 3,160 are more highly expressed and 3,088 are less expressed in the liver as compared to gut tissue. Out of the 6,248 zebrafish DE genes, 4,473 have a human ortholog that we entered into iPathwayGuide. In summary, 25 pathways were found to be significantly impacted (Table ). In addition, 279 GO terms and 47 miRNAs were found to be significantly enriched. No diseases were represented in this analysis. Most of the genes in the “Ribosome, Oxidative phosphorylation, Carbon metabolism, Glycolysis /Gluconeogenesis, Fatty acid degradation and Complement & coagulation cascades” pathways are upregulated and most of the genes in the “Vitamin and Protein digestion & absorption” are downregulated in the liver tissue as compared to intestinal tissue (Supplementary Tables - and Supplementary Figs. -) bar plots and KEGG for each pathway).

Identification of Predicted Novel miRNAs

The CAP-miR-Seq pipeline also allows the discovery of the predicted novel miRNAs. Following conservative criteria (count ≥ 10 and detected in both liver and both gut samples), we uncovered 21 novel miRNAs in the liver and 28 novel miRNAs in the gut. Our data show that the expression level of novel miRNAs is considerably lower on average compared to the expression level of already known miRNAs (Fig. ), as reflected in the average (ave) and the maximum (max) count obtained in each dataset. A comparison revealed that 8 novel miRNAs are commonly expressed in the liver and the gut tissues while 13 and 20 are uniquely expressed in the liver and the gut. respectively (Fig. ). It is important to note that on chromosome 1, there are 4 isoforms of the same miRNA with genomic coordinates, a couple of base pairs apart from nucleotide 10170922-10170926; this miRNA is commonly expressed in both tissues.

DISCUSSION

To our knowledge, most miRNome studies on zebrafish have focused on embryonic stage [51-56] while a few have examined miRNAs in the adult brain [57] and liver [58]. Her et al. [58] demonstrated that gankyrin, a chaperone during the assembly of the 26S proteasome, contributes to the development of liver steatosis via dysregulation of hepatic miRNAs miR-16, miR-27b, miR-122, and miR-126, and elevated gene expression linked to lipid metabolism and apoptosis. Consequently, the goals of our study were to (1) fully characterize the liver and gut miRNomes in healthy animals under homeostasis, without exposure to endocrine disrupting chemicals (2) identify predicted novel miRNAs, and (3) examine the regulatory networks that exist between the miRNome and transcriptome in each tissue to better understand the impact that miRNAs have on biological processes, pathways and molecular functions.

Here we showed that: (1) 224 and 233 miRNAs are expressed in the liver and the gut tissues, respectively. Contrary to Plasterk et al. [59] who concluded that many miRNAs are tissue-specific and exhibit striking organ-specific expression, we determined that 215 miRNAs are commonly expressed in both tissues, and that only 9 and 18 miRNAs are uniquely expressed in the liver and the gut respectively. (2) 5,988 and 5,817 mRNAs are expressed in the liver and the gut, respectively, 4,438 genes are commonly expressed in both tissues, 1,550 and 1,379 mRNAs are expressed uniquely in the liver and the gut, respectively, and (3) 4,438 commonly expressed genes regulate basic cellular processes and biological pathways (i.e. cellular metabolism, oxidation-reduction process and mitochondrial energy production), while the unique transcriptomes regulate processes and pathways that are phenotypical diverse: 1,550 genes unique to the liver regulate RNA processing and metabolism, mitochondrial organization and methyltransferase activity while 1,379 genes unique to the gut regulate lipid and sphingolipid metabolism, immune system process, protein binding and transporter activity molecular functions as well as intestine-specific cell projections components such as brush border and microvillus.

Our laboratory focuses on the study of endocrine disruption using zebrafish as a model organism for toxicology research. A profusion of both organic and inorganic contaminants has been introduced into the environment over the past century. Many are manufacturing chemicals and possess structural homology to steroid hormones to be able to bind to hormone receptors or enzymes that control the levels of steroid hormones [60, 61] and thus disrupt normal endocrine physiology in aquatic species, humans, and wildlife [62-64]. The prelude to examining the impact of CEC on biological systems such as the zebrafish is to first have a clear understanding of their normal homeostatic state. Here we investigated the interplay between miRNAs and their predicted target mRNAs in cellular homeostasis in the liver and the gut of adult male zebrafish and performed a systems level and GO analyses to capture the key biological processes and pathways regulated by this network. Future work will examine the effects of endocrine disruptors on this network.

Lessons Learned from the Matrix of miRNAs-mRNAs

We constructed a matrix of interacting miRNAs-mRNAs and performed GO analysis of the targeted mRNAs. This allowed us to explore two interesting possibilities. In the 1st scenario, many of the top most enriched GO terms identified from the unique miRNAs expressed in the liver and the gut were nearly the same in both tissues, including cellular metabolism and catabolism, oxidation-reduction process, RNA and cofactor binding, oxidoreductase activity, mitochondrial part, endoplasmic reticulum and envelope, suggesting that they are tightly regulated by different miRNAs. Previous studies suggested that well-conserved miRNAs are responsible for the control of fundamental cellular processes [59, 65]. Additionally, many of the conserved miRNAs have a salient organ-specific expression signature in zebrafish [52]. Given that dre-miR-144-3p and -miR-460-5p (Fig. ) as well as dre-miR-137-3p and -miR-7a (Fig. ) are well conserved (Supplementary Data: Conservation analysis) and have numerous predicted target genes found in the liver and gut transcriptomes, it is consistent with the studies cited above that they regulate similar basic biological processes and pathways in these two tissues.

In the 2nd scenario, many of the top GO terms identified from the commonly expressed miRNAs were quite different, suggesting that this pool of 164 miRNAs regulates different biology depending on the tissue, and more specifically which cell type, they are expressed from. A great example of this multi-faceted personality of miRNA is miR-133. Although miR-133 is expressed mainly in cardiomyocytes where it regulates hypertrophic genes such as RhoA, Cdc42, NELFA, WHSC2 and IP3RII [66, 67] as well as proliferative and apoptotic genes including CCND2 and SRF [68, 69], it is also expressed in cardiac fibroblasts and its function in regulating fibrosis, by targeting fibrotic genes such as Col1A1 and CTGF [70, 71], has been shown by in vitro and in vivo studies [68, 72].

Discussion DE Analysis: Liver Versus Intestine

The DE analysis revealed 25 enriched pathways, most of which contained upregulated genes (in the liver as compared to gut) such as “Ribosome, Oxidative phosphorylation, Carbon metabolism, Glycolysis/Gluconeogenesis, Fatty acid degradation and Complement & coagulation cascades” pathways, highlighting that the liver is an extremely active organ that conducts a wide variety of functions including vascular functions (under normal physiological conditions, the liver comprises 10-15% of the total body blood volume, it functions as a supply center when the peripheral circulation is encumbered and releases blood into the circulation when a blood loss occurs, and synthesizes approximately 50% of all circulating lymph), immunological functions (the Kupffer cells act as macrophages and establish a component of the body’s phagocytic system), metabolic functions (carbohydrate metabolism such as glycogenesis/glycogenolysis/ gluconeogenesis, protein metabolism including transamination/deamination/urea synthesis, lipid metabolism such as lipogenesis/fatty acid oxidation/lipoprotein synthesis/phos- pholipid/cholesterol synthesis), secretory functions (the liver secretes bile acids required for the digestion and absorption of fat and fat-soluble vitamins from the small intestine) and excretory & detoxification functions [18, 73]. In fact, after the skin, the liver is the largest organ in the body, comprising up to 2.5% of lean body mass [73].

Several enriched pathways contained downregulated genes in the liver as compared to gut, including “Vitamin and Protein digestion & absorption”, two pathways that are clearly more relevant in the gut tissue since the intestine is responsible for the breakdown (digestion) and absorption of various foods and liquids needed to sustain life. The large intestine absorbs vitamins from food as well as vitamins made by intestinal bacteria such as Vitamin K (important for the process of blood clotting), B1, B2, B6 and B12 and Biotin [74].

The presence of Parkinson’s and Alzheimer’s among the significant pathways reported in Table is due to pathway cross-talk. Pathway cross-talk is a phenomenon through which a pathway influences another either directly through common genes or through inter-pathway signals or linkage. The cross-talk phenomenon is well known and has been reported in two landmark papers [75, 76]. In particular, Donato et al. [75] also reported the very same pathways, i.e. Parkinson’s and Alzheimer appearing as significant in an experiment completely unrelated to these diseases. They investigated the phenomenon in great detail and identified the cause as being the inclusion of mitochondrial genes in the Alzheimer’s, Parkinson’s and Huntington’s disease.

Discussion Novel miRNAs

Several novel miRNAs were identified in our analysis. We uncovered 21 novel miRNAs in the liver and 28 novel miRNAs in the gut, amongst which, 8 are commonly expressed in both tissues while 13 and 20 are uniquely expressed in the liver and the gut respectively. The expression level of novel miRNAs is much lower on average than the expression level of already known miRNAs, which is consistent with the previous reports that concluded that novel miRNAs are often species-specific and weakly expressed [50, 77]. Further examination is needed to clearly define the functions of these novel miRNAs in a tissue-specific context and at different developmental stages.

CONCLUSION

Using high throughput sequencing, we fully characterized the miRNomes and transcriptomes of the liver and the gut in adult male zebrafish, as well as uncovered novel miRNAs. Our miRNAs-mRNAs matrix revealed interesting characteristics for miRNAs that are either uniquely or commonly expressed in a given tissue, characteristics that may be associated with how well-conserved an miRNA is and its function. This study provides insightful data on these tissues at homeostasis and provides a platform for further investigation into the role of the miRNome network in zebrafish, a field that is still in its infancy.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE

The study was approved by the University of California San Diego, Institutional Animal Care and Use Committee; Animal Use Protocol (AUP) No. S09418.

HUMAN AND ANIMAL RIGHTS

No humans were involved in this study. All procedures strictly followed the University of California San Diego, Institutional Animal Care and Use Committee (IACUC) Guidelines.

CONSENT FOR PUBLICATION

Not applicable.

Table 1

Impact on the liver unique transcriptome of the 215 miRNAs commonly expressed in the liver and the gut (2nd scenario). Only the top 20 miRNAs are shown here based on the number of predicted genes found.

miRNA	Predicted Genes Found	% of Targets	% of Expressed Genes
dre-miR-2184	290	28.7	18.7
dre-miR-23b	284	28.1	18.3
dre-miR-23a-3p	284	28.1	18.3
dre-miR-101a	254	25.2	16.4
dre-miR-101b	253	25.1	16.3
dre-miR-141-3p	246	24.4	15.9
dre-miR-200a-3p	246	24.4	15.9
dre-miR-181c-5p	235	23.3	15.2
dre-miR-181b-5p	235	23.3	15.2
dre-miR-181a-5p	234	23.2	15.1
dre-miR-30e-5p	222	22.0	14.3
dre-miR-30c-5p	221	21.9	14.2
dre-miR-30b	220	21.8	14.2
dre-miR-30a-5p	217	21.5	14.0
dre-miR-30d	216	21.4	13.9
dre-miR-19b-3p	215	21.3	13.9
dre-miR-19d-3p	215	21.3	13.9
dre-miR-19a-3p	215	21.3	13.9
dre-miR-19c-3p	215	21.3	13.9
dre-miR-725-3p	209	20.7	13.5

Table 2

Impact on the gut unique transcriptome of the 215 miRNAs commonly expressed in the liver and the gut (2nd scenario). Only the top 20 miRNAs are shown here based on the number of predicted genes found.

miRNA	Predicted Genes Found	% of Targets	% of Expressed Genes
dre-miR-2184	329	34.6	23.8
dre-miR-101a	323	34.0	23.4
dre-miR-101b	322	33.9	23.3
dre-miR-181c-5p	319	33.6	23.1
dre-miR-181a-5p	319	33.6	23.1
dre-miR-181b-5p	319	33.6	23.1
dre-miR-23a-3p	302	31.8	21.9
dre-miR-23b	300	31.6	21.7
dre-miR-19d-3p	256	26.9	18.6
dre-miR-19c-3p	256	26.9	18.6
dre-miR-19b-3p	256	26.9	18.6
dre-miR-19a-3p	256	26.9	18.6
dre-miR-30e-5p	246	25.9	17.8
dre-miR-30b	243	25.6	17.6
dre-miR-30c-5p	243	25.6	17.6
dre-miR-725-3p	242	25.5	17.5
dre-miR-93	240	25.3	17.4
dre-miR-17a-5p	239	25.2	17.3
dre-miR-30d	238	25.1	17.2
dre-miR-20b-5p	238	25.1	17.2

Table 3

iPathwayGuide - Biological Pathways. Only the top 20 pathways are shown here based on q-value.

Biological Pathways	q-value
Ribosome	1.0E-24
Metabolic pathways	3.2E-21
Oxidative phosphorylation	1.3E-10
Carbon metabolism	3.7E-07
Glycolysis/Gluconeogenesis	2.3E-06
Systemic lupus erythematosus	3.5E-05
Fatty acid degradation	7.2E-05
Tryptophan metabolism	1.1E-04
Complement and coagulation cascades	1.2E-04
Parkinson’s disease	1.8E-04
Alzheimer’s disease	2.3E-04
Non-Alcoholic Fatty Liver Disease (NAFLD)	2.7E-04
Huntington’s disease	2.7E-04
Staphylococcus aureus infection	1.5E-03
Vitamin digestion and absorption	2.0E-03
Pyruvate metabolism	2.2E-03
Protein digestion and absorption	2.2E-03
Glycine, serine and threonine metabolism	2.4E-03
Lysine degradation	2.4E-03
Biosynthesis of amino acids	3.2E-03