Integration of miRNA and mRNA Co-Expression Reveals Potential Regulatory Roles of miRNAs in Developmental and Immunological Processes in Calf Ileum during Early Growth

This study aimed to investigate the potential regulatory roles of miRNAs in calf ileum developmental transition from the pre- to the post-weaning period. For this purpose, ileum tissues were collected from eight calves at the pre-weaning period and another eight calves at the post-weaning period and miRNA expression characterized by miRNA sequencing, followed by functional analyses. A total of 388 miRNAs, including 81 novel miRNAs, were identified. A total of 220 miRNAs were differentially expressed (DE) between the two periods. The potential functions of DE miRNAs in ileum development were supported by significant enrichment of their target genes in gene ontology terms related to metabolic processes and transcription factor activities or pathways related to metabolism (peroxisomes), vitamin digestion and absorption, lipid and protein metabolism, as well as intracellular signaling. Integration of DE miRNAs and DE mRNAs revealed several DE miRNA-mRNA pairs with crucial roles in ileum development (bta-miR-374a—FBXO18, bta-miR-374a—GTPBP3, bta-miR-374a—GNB2) and immune function (bta-miR-15b—IKBKB). This is the first integrated miRNA-mRNA analysis exploring the potential roles of miRNAs in calf ileum growth and development during early life.


Introduction
MicroRNAs (miRNAs) are small (~22 nucleotides) endogenous RNA molecules that regulate gene expression post-transcriptionally by targeting principally the 3 untranslated region (3 UTR) of genes and to some extend the 5 UTR, introns and coding region of mRNAs [1]. They play key roles in a wide range of biological processes [2]. In bovine, miRNAs have been shown to play important roles in embryonic development [3][4][5], mammary gland [6,7] and adipose tissue [8] functions and also in the regulation of production traits such as milk yield [6,9,10], milk quality [11] and diseases like mastitis [12,13] and Johne's disease [14]. The importance of miRNAs in gut development and disease (mostly inflammatory bowel disease) has been extensively studied in humans [15][16][17][18]. For instance, several miRNAs have been reported to be relevant for different aspects of inflammatory bowel disease (IBD), including miRNAs important for intestinal fibrosis (miR-29, miR-200b, miR-21, miR-192),

Identification of Known miRNA and Novel miRNA Discovery
The identification of known miRNAs was performed with miRBase v21 http://www.mirbase. org/) (Kozomara and Griffiths-Jones, 2014), while novel miRNA discovery was achieved with miRDeep2 v2.0.0.8 (https://github.com/rajewsky-lab/mirdeep2) [27]. MiRDeep2 was designed to detect miRNAs from deep sequence reads using a probabilistic algorithm based on the miRNA biogenesis model. The core and quantifier modules of miRDeep2 were applied to discover novel miRNAs in the pooled dataset of all the libraries while the quantifier module was used to profile the detected miRNAs in each library. MiRDeep2 score higher than five was used as cuff point for the identification of novel miRNAs. Subsequently, a threshold of 10 counts per million and present in ≥2 libraries was applied to remove lowly expressed miRNAs. MiRNAs meeting these criteria were further used in downstream analyses including differential expression (DE) analysis.

Predicted Target Genes of miRNAs
In order to investigate the functions of the most highly expressed miRNAs and differently expressed (DE) miRNAs, we firstly predicted their target mRNAs. Perl scripts (targetscan_60.pl and targetscan_61_context_scores.pl) (http://targetscan.org) were used to predict target mRNAs and to calculate their context scores, respectively. Predicted target mRNAs with context + scores above 95th percentile were further used [9,10,30]. The predicted target mRNAs were then filtered against the mRNA transcriptome obtained from ileum tissues of the same animals. Only predicted target genes that were expressed in the mRNA transcriptome of the ileum tissues of the animals [23] were retained for further analysis.

miRNA-mRNA Co-Expression Analysis and Target Gene Enrichment
For miRNA-mRNA co-expression, the Pearson correlation coefficient between target mRNAs (retained above) and DE miRNAs were calculated. A miRNA-mRNA pair was considered co-expressed if it had a negative and significant correlation value at FDR < 0.05. The mRNAs significantly correlated with miRNAs were then used for downstream target gene ontology and KEGG pathways enrichment using ClueGO (http://apps.cytoscape.org/apps/cluego) [31]. For ClueGO analysis, a hypergeometric test was used for enrichment analyses and Benjamini-Hochberg [29] correction was used for multiple testing correction (FDR < 0.05). Since KEGG pathways enrichment relied on the human database (due to lack of information in bovine), we used a less stringent threshold (uncorrected p-value < 0.05) to declare if a pathway was significantly enriched. Interactions between miRNAs and mRNAs were visualized with Cytoscape (http://www.cytoscape.org/) [32].

Real-Time Quantitative PCR
The method of real-time quantitative PCR was used to validate the expression of four DE (bta-miR-142-5p, miR-146a, miR-24-3p and miR-374b) and two non-DE (bta-miR-486-5p and miR-193b) miRNAs. The same total RNA used in miRNA-sequencing was used. Total RNA was reverse transcribed with Universal cDNA Synthesis Kit II from Exiqon (Exiqon Inc., Woburn, MA, USA), following the manufacturer's instructions. ExiLENT SYBR ® Green Master Mix Kit (Exiqon, Woburn, MA, USA) and the miRCURY LNA™ Assay (Exiqon, Woburn, MA, USA) specific for each miRNA listed above were used to perform Quantitative qPCR on a StepOne Plus System (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions. Bta-miR-126-3p was used as endogenous control to assess the expression level of miRNAs using the comparative Ct (∆∆Ct) method. Bta-miR-126-3p was selected as an endogenous control based on its consistent expression throughout all the analyzed samples on D33 and D96.
Novel miRNAs were considered to have a minimum MiRDeep2 score of five, as shown in Table S3. After removing lowly expressed reads, a total of 307 known and 81 novel miRNAs satisfying the conditions of having at least 10 read counts per million and present in a minimum of two libraries were used for DE analysis (Table S4a,b).

Figure 1.
Volcano plot depicting miRNA differential expression results. Each dot represents a miRNA. Green and red dots represent miRNAs significantly differentially expressed at FDR < 0.05 and FDR < 1 × 10 11 , respectively. Black dots represent miRNAs that were not differentially expressed. Differential expression analysis was accomplished with DeSeq2. The DE miRNAs (220 miRNAs) were predicted to target 11,691 mRNAs (Table S6b). Using mRNA transcriptome data of the same samples, 1560 mRNAs out of the predicted 11,691 mRNAs, were significantly and negatively correlated with their targeting miRNAs (Table S6c). Bta-miR-2285f had the highest number of target genes (172), while AGO2 gene was the most popular target for DE miRNAs (targeted by 25 DE miRNAs) ( Table S6c). Other common target genes for DE miRNAs were SLC25A46, KCTD13 and PAXIP1, each targeted by 9 DE miRNAs ( Table S6c). The GO enrichment analyses of the 1560 target genes (significantly and negatively correlated with miRNA) indicated that 158, 26 and 28 of them were significantly enriched in BP-, CC-and MF-GO terms, respectively (Table S7a-c). The most enriched BP-, CC-and MF-GO terms were cellular macromolecule metabolic process (FDR = 9.38 × 10 10 ), intracellular (FDR = 3.37 × 10 19 ) and organic cyclic compound binding (FDR = 1.19 × 10 4 ), respectively ( Table 5, Figures 2-4). Moreover, 16 KEGG pathways were significantly enriched for the target genes (1560) of 220 DE miRNAs, and peroxisome (p = 0.004) and Hedgehog signaling pathways (p = 0.006) were the most significantly enriched (Table S7d, Figure 5). Moreover, among the 1560 target genes negatively correlated with miRNAs, 278 were also significantly DE between D33 and D96 in our previous study (Table S8) [23]. The 278 genes were the targets for 64 DE miRNAs. SOX4 was the most common target, since it was targeted by 6 different miRNAs (bta-miR-191, bta-miR-30e-5p, bta-miR-15-11508, bta-miR-2285f, bta-miR-92b and bta-miR-2285q). Meanwhile, bta-miR-2285f and bta-miR-874 had the highest number of target genes (37 and 28, respectively) ( Figure 6).  Figure 2. The ClueGO results for biological processes gene ontology terms enrichment for target genes (mRNAs) of differentially expressed miRNAs and relationships between them. The nodes (round shape) represent gene ontology terms, node color represents the level of significance as indicated in the legend, while node size reflects the number of genes enriched in each gene ontology term.

Figure 2.
The ClueGO results for biological processes gene ontology terms enrichment for target genes (mRNAs) of differentially expressed miRNAs and relationships between them. The nodes (round shape) represent gene ontology terms, node color represents the level of significance as indicated in the legend, while node size reflects the number of genes enriched in each gene ontology term. The nodes (round shape) represent gene ontology terms, node color represents the level of significance as indicated in the legend, while node size reflects the number of genes enriched in each gene ontology term.

Real Time Quantitative PCR Validation
The RNA-Seq expression of 6 miRNAs was validated by qPCR. Two of them (bta-miR-486-5p and bta-miR-193b) were non DE, while four of them (bta-miR-142-5p, miR-146a, miR-24-3p and miR-374b) were DE between D33 and D96 by RNA-Seq. Observed fold changes for DE miRNAs between both methods were similar, except for the non-DE miRNAs, where an opposite trend was observed (Figure 7).

Real Time Quantitative PCR Validation
The RNA-Seq expression of 6 miRNAs was validated by qPCR. Two of them (bta-miR-486-5p and bta-miR-193b) were non DE, while four of them (bta-miR-142-5p, miR-146a, miR-24-3p and miR-374b) were DE between D33 and D96 by RNA-Seq. Observed fold changes for DE miRNAs between both methods were similar, except for the non-DE miRNAs, where an opposite trend was observed (Figure 7).

Discussion
Physiologically, major metabolic changes that take place in calf gastrointestinal tract following the transition from liquid to solid food are accompanied by rapid changes in gene expression controlled by the signal-mediated coordination of transcriptional and post-transcriptional mechanisms [33]. Previously, we reported that ~20% of bovine genes were significantly DE between the pre-weaning (day 33) and the post-weaning (day 96) period, and enrichment analysis revealed the importance of DE genes in biological processes necessary for the switch in nutrition and developmental stage from the pre-weaning to the post-weaning period [23]. While it is well known that miRNAs are important for the regulation of these processes, little is known about how they participate in the regulation of ileum functions from the pre-weaning to the post-weaning period.
Highly expressed and DE miRNAs identified in this study suggest potential roles in ileum developmental processes during the transition from the pre-weaning to the post-weaning period. However, some highly expressed miRNAs such as bta-miR-21-5p, bta-miR-26a, bta-miR-148a and bta-let-7a-5p (Table 1) were also highly expressed in other tissues such as milk fat [10], milk whey/somatic cells [25], and mammary gland epithelial cells [13]. Interestingly, bta-miR-143, the most highly expressed miRNA, was also among the most abundant miRNAs in bovine testis [34] and reported as the most highly expressed miRNA in ileum tissue of calves at the pre-and post-weaning periods [35]. It was suggested that bta-miR-143 might regulate key genes involved in differentiation of connective tissue cells, the major components of the gut; hence, its high abundance might be important for the regulation of rapid development and growth of the gastrointestinal tract during early life. Indeed, enrichment analysis of target genes of the top 20 abundant miRNAs indicated enrichment in many biological processes and molecular function GO terms (Tables 2 and S5) involved in developmental processes, therefore supporting important roles for highly abundant miRNAs, including bta-miR-143, in these processes. Further supporting evidence was derived from enriched KEGG pathways crucial for cellular development processes such as FoxO signaling pathway, cell cycle, MAPK signaling pathway, and p53 signaling pathway (Table 3). Moreover, we also detected 81 novel miRNAs in this study that were lowly expressed. The most abundant novel miRNA, bta-miR-22_24033, accounted for only 0.3% of the total read

Discussion
Physiologically, major metabolic changes that take place in calf gastrointestinal tract following the transition from liquid to solid food are accompanied by rapid changes in gene expression controlled by the signal-mediated coordination of transcriptional and post-transcriptional mechanisms [33]. Previously, we reported that~20% of bovine genes were significantly DE between the pre-weaning (day 33) and the post-weaning (day 96) period, and enrichment analysis revealed the importance of DE genes in biological processes necessary for the switch in nutrition and developmental stage from the pre-weaning to the post-weaning period [23]. While it is well known that miRNAs are important for the regulation of these processes, little is known about how they participate in the regulation of ileum functions from the pre-weaning to the post-weaning period.
Highly expressed and DE miRNAs identified in this study suggest potential roles in ileum developmental processes during the transition from the pre-weaning to the post-weaning period. However, some highly expressed miRNAs such as bta-miR-21-5p, bta-miR-26a, bta-miR-148a and bta-let-7a-5p (Table 1) were also highly expressed in other tissues such as milk fat [10], milk whey/somatic cells [25], and mammary gland epithelial cells [13]. Interestingly, bta-miR-143, the most highly expressed miRNA, was also among the most abundant miRNAs in bovine testis [34] and reported as the most highly expressed miRNA in ileum tissue of calves at the pre-and post-weaning periods [35]. It was suggested that bta-miR-143 might regulate key genes involved in differentiation of connective tissue cells, the major components of the gut; hence, its high abundance might be important for the regulation of rapid development and growth of the gastrointestinal tract during early life. Indeed, enrichment analysis of target genes of the top 20 abundant miRNAs indicated enrichment in many biological processes and molecular function GO terms (Tables 2 and S5) involved in developmental processes, therefore supporting important roles for highly abundant miRNAs, including bta-miR-143, in these processes. Further supporting evidence was derived from enriched KEGG pathways crucial for cellular development processes such as FoxO signaling pathway, cell cycle, MAPK signaling pathway, and p53 signaling pathway (Table 3). Moreover, we also detected 81 novel miRNAs in this study that were lowly expressed. The most abundant novel miRNA, bta-miR-22_24033, accounted for only 0.3% of the total read counts. Forty-four novel miRNAs were significantly DE between D33 and D96, therefore suggesting roles in the regulation of ileum gene expression during the early period of growth. Nevertheless, novel miRNAs identified in this study will enrich the bovine miRNome as well as enhance knowledge of the potential roles of miRNAs in calf GIT. However, further functional validations to clarify the roles of identified miRNAs in the development of calf gut during the early period of growth are needed.
Bta-miR-374a was the most significantly DE miRNA in this study (Table 4). Bta-miR-374a was found to be differentially expressed between lactating and non-lactating cows [36]. Bta-miR-374a potentially targeted 36 different genes ( Table S6b) and some of them might be important for ileum functions such as EIF2AK4, FBXO18, GTPBP3 and GNB2, etc. EIF2AK4 is an important transcription factor in host response to infection with pathogenic bacteria associated with Crohn's disease [37]. Moreover, FBXO18, GTPBP3 and GNB2 have been reported to be significantly DE in calf gastrointestinal tract between the pre-and post-weaning periods [23]. FBXO18 encodes a member of the F-box protein family with function in phosphorylation-dependent ubiquitination [38]. FBXO18 is implicated in the regulation of stress-induced apoptosis processes and homologous recombination in familial and sporadic breast cancer [39]. Cells deficient in FBXO18 were unable to activate the cytotoxic-stress-induced cascade, resulting in increased cell survival [40]. GTPBP3 is an important gene for mitochondrial functions and a mutation in this gene resulted in defective mitochondrial energy production through oxidative phosphorylation [41]. GNB2 is important for neuronal apoptosis and was induced by lidocaine in the rat [42]. Nevertheless, the functions of these genes (FBXO18, GTPBP3 and GNB2) in the ileum are unknown. Bta-miR-15b, the second most significantly up-regulated miRNA, belongs to the miR-15b family cluster. This miRNA cluster can target cell cycle proteins and the anti-apoptotic Bcl-2 gene to control cell proliferation and apoptosis [43]. Bta-miR-15b might also play a role in mastitis disease development in cows [44]. Moreover, Liang et al. [35] also reported that bta-miR-15b was significantly DE between 0-day-old and 7-day-old calves and its expression correlated with bacterial population, thus suggesting roles in the regulation of gut development, immune, and digestive functions. Furthermore, bta-miR-15b potentially targeted IKBKB (Figure 6), a gene known as an essential molecule for NF-κB signaling pathway [45] with important roles in both innate and acquired immunity [46].
Interestingly, bta-miR-26a, one of the most significantly DE miRNA (Table 4), was also one of the most highly expressed miRNA ( Table 1). The human homologue of this miRNA plays important roles in Crohn's disease [15]. In cattle, this miRNA regulates the expression of PCK1 gene, which is important for semen quality and longevity of Holstein bulls [47]. Bta-miR-26a potentially targeted BID gene and their expressions were significantly correlated in this study ( Figure 6). BID is a pro-apoptotic member of the Bcl-2 protein family with roles in the regulation of apoptosis processes [46]. Therefore, bta-miR-26a might have roles in ileum development via targeting the BID gene. Bta-miR-455-5p was the most down-regulated miRNA between D33 and D96. Bta-miR-455-5p was reported to be important for the function of granulosa cells of subordinate and dominant follicles during the early luteal phase of the bovine estrous cycle [48]. In humans, this miRNA homologue down regulated RAB18 gene in gastric cancer [49]. In fact, bta-miR-455-5p also potentially targeted RAB18 gene in this study ( Figure 6).
Among the top 20 DE miRNAs in this study, bta-miR-142-3p, bta-miR-142-5p, bta-miR-191, bta-miR-146a, bta-miR-210 and bta-miR-424-5p were also found to be DE in calf ileum at the pre-weaning and weaning periods [35]. Some of these miRNAs have been reported to play important roles in immune functions. For example, bta-miR-146a inhibited the mRNA and protein expression levels of TRAF6 gene and acted as a negative feedback regulator of bovine inflammation and innate immunity through down regulation of the TLR4/TRAF6/NF-κB pathway in bovine mammary epithelial cells [50]. Furthermore, bta-miR-142-5p was important for bovine alveolar macrophage response to Mycobacterium bovis infection [51].
As expected, the target genes of DE miRNAs were enriched in important biological process GO-terms related to metabolic processes (such as cellular macromolecule metabolic process, macromolecule metabolic process and regulation of metabolic process) ( Table 5 and Figure 2) and molecular function GO-terms related to metabolism of the macromolecule compound (such as organic cyclic compound binding, heterocyclic compound binding, nucleic acid binding and protein binding) ( Table 5 and Figure 4), thus suggesting roles in the regulation of these processes. Interestingly, the target genes of DE miRNAs were also enriched in GO-terms like transcription factor activity and sequence-specific DNA binding (FDR = 3.63 × 10 −3 ) (Table 5), thus suggesting their importance in transcription factor activity. The interaction between miRNAs and transcription factors to regulate gene expression in biological processes is well documented [52,53]. MiRNAs might inhibit transcription factor activities by either directly inhibiting the expression of their encoding genes or by inhibiting other gene(s) that have impact on their activities [53,54]. Previously, we observed that some transcription factors might play important roles in mediating miRNA regulatory functions in cow milk yield and milk component traits [9]. In humans, several miRNAs have been reported to participate in the regulation of intestinal transcription factors; miR-196b inhibited the GATA6 intestinal transcription factor to control intestinal cell homeostasis and tumorigenesis in colon cancer patients [55], while miR-30 family controlled intestinal epithelial cell proliferation and differentiation by targeting SOX9 (transcription factor) and other genes in ubiquitin ligase pathway [56]. In fact, as mentioned above, the most DE miRNA in this study (bta-miR-374) also potentially targeted a transcription factor (EIF2AK4) known to be important for human Crohn's disease [37]. The most important pathway enriched for DE miRNAs target genes was Hedgehog signaling pathway (p = 0.003, Table S7, Figure 5). Hedgehog signaling pathway is important for cell growth, survival and fate, as well as for normal embryonic development [57,58]. This pathway also has multiple patterning functions during mammalian gut development [59]; therefore, it may be important for ileum functions during the early part of life. Another important pathway enriched for target genes of DE miRNAs was peroxisomes pathway (p = 0.006, Table S7, Figure 5). The peroxisomes pathway is crucial for metabolic processes such as fatty acid oxidation, biosynthesis of ether lipids, and free radical detoxification [60]. Since one of the main functions of the ileum is to absorb bile salts, one of the products of fatty acid oxidation, enrichment of the peroxisome pathway supports its role in normal ileum function. Another important role of the ileum is vitamin absorption and the enriched pathway, vitamin digestion and absorption (Table S7 and Figure 5), might reflect the changes in gene expression for different vitamin requirements between the pre-and post-weaning periods. Other enriched pathways also reflect the importance of miRNAs in the regulation of genes involved in lipid metabolism (phosphatidylinositol signaling system and sphingolipid signaling pathway), protein metabolism (valine, leucine and isoleucine biosynthesis and protein processing in endoplasmic reticulum) and intracellular signaling (cGMP-PKG signaling, FoxO signaling and neurotrophin signaling pathway) in the development of the ileum during the early part of life.

Conclusions
This is the first integrated miRNA-mRNA analysis characterizing the function of miRNAs in calf ileum during early life. Eighty-one novel miRNAs were identified that will enrich the bovine miRNome repertoire and contribute to the understanding of regulatory processes in calf ileum. This study highlighted potential roles of bta-miR-143, bta-miR-192, bta-miR-26a and bta-miR-21-5p in growth and developmental processes during the transition from the pre-weaning to the post-weaning period. This study also suggested roles for DE miRNAs in metabolic processes, metabolism of the macromolecule compound, transcription factor activities, as well as involvement in pathways related to metabolism (peroxisomes), vitamin digestion and absorption, lipid and protein metabolism, and intracellular signaling (Hedgehog signaling, GMP-PKG signaling, FoxO signaling, neurotrophin signaling pathway). Moreover, several DE miRNAs-DE mRNAs pairs such as bta-miR-374a-FBXO18, bta-miR-374a-GTPBP3 and bta-miR-374a-GNB2 with potential roles in tissue development, and bta-miR-15b-IKBKB with vital roles in immune functions were revealed. This study, therefore, provided insights on miRNA expression and their potential functions in calf ileum development during early life, which might facilitate identification of miRNA biomarkers for growth, nutritional and disease challenges during the pre-and post-weaning periods.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4409/7/9/134/s1. Table S1: Mapping statistics of miRNA sequencing reads; Table S2: Length distribution (nt) of miRNA reads; Table  S3: Novel miRNAs identified by miRDeep2; Table S4: Novel and known miRNAs and their read count summary; Table S5; The 20 most highly expressed miRNAs and their predicted target genes (a); enriched gene ontology terms (b) and enriched KEGG pathways (c); Table S6; Differentially expressed miRNAs between day 33 (pre-weaning) and day 96 (post-weaning) (a); predicted target genes of DE miRNAs (b) and differentially expressed genes which are both negatively correlated and predicted targets of miRNAs (c); Table S7; Gene ontology and KEGG pathways enriched for target genes of differential expressed miRNAs; Table S8; Differentially expressed genes that were targets of DE miRNAs and also negatively correlated.