Integration of lncRNA and mRNA Transcriptome Analyses Reveals Genes and Pathways Potentially Involved in Calf Intestinal Growth and Development during the Early Weeks of Life

A better understanding of the factors that regulate growth and immune response of the gastrointestinal tract (GIT) of calves will promote informed management practices in calf rearing. This study aimed to explore genomics (messenger RNA (mRNA)) and epigenomics (long non-coding RNA (lncRNA)) mechanisms regulating the development of the rumen and ileum in calves. Thirty-two calves (≈5-days-old) were reared for 96 days following standard procedures. Sixteen calves were humanely euthanized on experiment day 33 (D33) (pre-weaning) and another 16 on D96 (post-weaning) for collection of ileum and rumen tissues. RNA from tissues was subjected to next generation sequencing and 3310 and 4217 mRNAs were differentially expressed (DE) between D33 and D96 in ileum and rumen tissues, respectively. Gene ontology and pathways enrichment of DE genes confirmed their roles in developmental processes, immunity and lipid metabolism. A total of 1568 (63 known and 1505 novel) and 4243 (88 known and 4155 novel) lncRNAs were detected in ileum and rumen tissues, respectively. Cis target gene analysis identified BMPR1A, an important gene for a GIT disease (juvenile polyposis syndrome) in humans, as a candidate cis target gene for lncRNAs in both tissues. LncRNA cis target gene enrichment suggested that lncRNAs might regulate growth and development in both tissues as well as posttranscriptional gene silencing by RNA or microRNA processing in rumen, or disease resistance mechanisms in ileum. This study provides a catalog of bovine lncRNAs and set a baseline for exploring their functions in calf GIT development.


Introduction
From birth, calves undergo major physiological and dietary changes including adaptation to extrauterine life (first week of life), pre-ruminant stage (3 to 5 months of age) and weaning [1]. After birth, sufficient colostrum ingestion and microbial colonization is crucial for adequate development of the gastro-intestinal tract (GIT) and mucosal immune system [1,2]. These developmental transitions are accompanied by rapid changes in gene expression controlled by signal-mediated coordination of transcriptional and post-transcriptional mechanisms [3]. The transcriptional mechanisms include the activities of non-coding RNA (ncRNA) and transcription factors as well as epigenetic modifications. While only about 2% of the mammalian genome is transcribed as proteins, about 75-90% is transcribed as ncRNA, the vast majority being long non-coding RNA (lncRNA) [4].

RNA Isolation
Total RNA from rumen (n = 32) and ileum (n = 16) tissues (30 mg/sample) was extracted using miRNeasy Kit (Qiagen Inc., Toronto, ON, Canada) following manufacturer's recommendations. Briefly, tissues were homogenized in 700 µL TRIzol Reagent (Life Technologies, Burlington, ON, Canada) with a Polytron homogenizer (Polytron PT 10-35 GT, Kinematica AG, Luzern, Switzerland) using a 7 mm probe for 10 s at 12,000 rpm, repeated three times with incubation on ice between repetitions. Following 5 min incubation at room temperature, 140 µL chloroform was added to the mixture and vortexed vigorously for 20 s. The mix was centrifuged (15 min at 12,000× g at 4 • C) and 1.5 volumes ethanol (100%) was added to the aqueous phase. Samples were then transferred to the column and washed steps performed according to kit's recommendations. Elution was done twice using 30 µL nuclease-free water each time. Total RNA (10 µg) was subjected to DNase treatment with Turbo DNA-free Kit (Ambion Inc. Foster City, CA, USA) to remove contaminating genomic DNA. Nanodrop ND-1000 (NanoDrop Technologies, Wilmington, DE, USA) was used to determine the concentration of RNA (before and after DNase treatment) and the integrity was assessed with Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) using the RNA 6000 Nano Labchip Kit (Agilent Technologies) after DNase treatment. All samples had a RNA integrity number (RIN) ≥8.0.

Library Preparation and RNA-Sequencing
Total RNA (2 µg) was subjected to ribosomal RNA (rRNA) depletion using Ribo-Zero Gold rRNA Removal Kit (Illumina Inc., San Diego, CA, USA) following company recommendations. Generation of libraries for sequencing was made using NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Whitby, ON, Canada) with NEBNext Multiplex Oligos for Illumina ® (New England Biolabs) for barcoding the multiplexed samples. Library concentration was measured using the Quant-iT PicoGreen double-stranded DNA (dsDNA) Assay Kit (Life Technologies). Fragment size was estimated using High Sensitivity DNA Analysis Kit (Agilent) with Agilent 2100 Bioanalyzer (Agilent) and quantified by real time quantitative polymerase chain reaction (qPCR) using the Kapa Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems, Wilmington, MA, USA).
Six libraries were pooled in equimolar amounts and paired-end sequenced (2 × 126 bp) in one lane on a High Throughput Model flow cell on an Illumina HiSeq 2500 system by The Centre for Applied Genomics, The Hospital for Sick Children, Toronto (http://www.tcag.ca/).

Sequence Data Processing, Alignment and Identification of Genes
Demultiplexed sequence files in fastq format were processed using a pipeline developed by McGill University and Genome Quebec Innovation Centre (MUGQIC, http://gqinnovationcenter.com/). Briefly, adaptor sequences were removed with Trimmomatic software v0.32 [37] which was set to keep reads longer than 32 bp with a minimum phred score of 30. Alignment of reads to the bovine genome (UMD3.1, v84) was accomplished with STAR v2.5.1b [38]. Following alignment, Picard tools [39] and RNA-SeQC [40] were used to generate quality control read metrics such as percentage of mapped reads and genomic localization (intronic, intergenic, etc.). Uniquely mapped and properly paired reads were used in transcript construction with Cufflinks (v2.2.1) [41] for mRNA identification. Constructed transcripts were compared with Ensembl bovine gene annotation (release 84) to identify expressed mRNAs using Cuffcompare [41].
For lncRNA identification, alignment files or properly mapped reads from all the samples of the same tissue were sorted and merged using Samtools v1.3.1 [42] prior to the transcript construction step with Cufflinks [41] to generate a unique set of all transcripts per tissue. To identify expressed lncRNAs, only transcripts >200 bp were kept and compared with Ensembl bovine gene annotation (release 84) to remove annotated transcripts (mRNA) and transcripts overlapping with other noncoding RNA species (transfer RNA (tRNA), rRNA, small nuclear RNA (snRNA), small nucleolar (snoRNA) and microRNA (miRNA)) using Cuffcompare [41]. Specifically, transcripts with class code "i" (predicted transcript fall entirely within a reference intron), "u" (predicted transcript is intergenic in comparison with known reference transcripts) and "x" (exon of predicted transcript overlaps reference but lies on the opposite strand) when compared against Ensembl-ncRNA databases were retained. Retained transcripts were then assessed for their coding potential using Coding-Non-Coding Index (CNCI) [43] and Coding-Potential Assessment Tool (CPAT) [44] and intersecting transcripts across CNCI score < 0 and CPAT score < 0.5 were blasted against the Swiss-prot database to further filter transcripts with the ability to code for a protein (evalue < 1 × 10 −5 ) using usearch [45] and also transcripts that possessed an open reading frame with the ability to code for a peptide of 100 or more amino acids.
The retained transcripts were compared with known bovine lncRNA annotation from NONCODE2016 database [46] using Cuffcompare. Transcripts with class codes "=" (predicted transcript has exactly the same introns as the reference transcript), "c" (predicted transcript is contained within the reference transcript) and "j" (predicted transcript is a potential novel isoform that shares at least one splice junction with a reference transcript) were classified as known bovine lncRNAs while the rest were classified as novel lncRNAs. The expression of lncRNAs (known and novel) and mRNAs were quantified in each sample using HTSeq-count (v0.6.1p1) [47] with default settings (-s reverse).

Differential mRNA and lncRNA Expression Analyses
Differential gene expression (DE) was accomplished with DeSeq2 (v1.14.1) [48], an R package [49] that implements a negative binomial model in DE analysis. Those genes (mRNA or lncRNA) with DESeq2 normalized counts >5 and in at least 10% of libraries were considered truly expressed and were used in DE and pathway analyses. The differences in gene expression (mRNA or lncRNA) between D33 and D96 were considered significant at Benjamini and Hochberg [50] corrected p-value < 0.05.

Gene Ontology and Pathways Enrichment of Differentially Expressed mRNAs
Differentially expressed genes were enriched for gene ontology (GO) and Kyoto Enciclopedia of Genes and Genomes (KEGG) pathways using ClueGO [51]. In this enrichment analysis, a hypergeometric test was used. For GO enrichment, the Bovine GO database and Benjamini and Hochberg (BH) [50] correction for multiple testing controlled p-values and GO with corrected p < 0.05 was considered significantly enriched. Meanwhile, since there are no KEGG databases for bovine species in ClueGO, the human KEGG database was used as background and pathways with uncorrected p < 0.05 were considered significantly enriched.

Identification of lncRNA cis Target Genes and cis Target Gene Enrichment
Since lncRNAs can cis regulate mRNAs, we performed co-expression analysis between mRNAs and lncRNAs in 50 Kb flanking regions of identified lncRNAs. Firstly, we computed Pearson's correlation coefficients between each pair of lncRNA-mRNA. The mRNAs having significant correlations with lncRNAs at p.BH < 0.05 were considered potential cis target genes for those lncRNAs. Potential cis target genes of lncRNAs were subjected to enrichment analysis with ClueGO using the same procedure as for DE mRNAs. The interaction between DE lncRNAs and their cis target genes were visualized using Cystoscope v3.0 [52].

Real Time Quantitative Polymerase Chain Reaction
Real time quantitative PCR was performed to verify the expression levels of genes DE by the method of RNA-Seq. Four potential housekeeping genes (PPIB, ATP5B, GPI and PGK1) were randomly selected from the list of non-DE genes by RNA-Seq in both tissues. Two DE genes (IFI6 and OAS1X) were randomly selected from the list of DE genes common to both tissues while three DE genes in ileum tissue (OAS2, Mx1 and UBA7) and three DE genes in rumen tissue (ACTA2, CA3 and HERC6) were randomly selected from the list of DE genes unique to these tissues. Moreover, two lncRNAs were selected from rumen DE lncRNAs (rXLOC_042149 and rXLOC_022071; rXLOC denotes lncRNAs expressed in the rumen). Transcript-specific primers for the two lncRNAs and 12 mRNAs were designed using Integrated DNA Technologies Assay tool for real time quantitative PCR (Integrated DNA Technologies Inc., Skokie, IL, USA) (Supplementary Table S1). Reverse transcription was performed with the SuperScript III Reverse Transcriptase (Life Technologies Inc., Burlington, ON, Canada), using aliquots (1 µg) of the same total RNA used in RNA-Seq. The complementary DNA (cDNA) samples were diluted to 20 ng/µL. Real-time PCR reaction mix was composed of 5 µL Power SYBR Green PCR Master Mix (Life Technologies Inc.), 3 µL cDNA and 0.1 U AmpErase Uracil N-Glycosylase (UNG) (Life Technologies). Forward and reverse primer concentrations in the mix are presented in Supplementary Table S1. StepOne Plus Real-Time PCR System (Life Technologies) was used to perform qPCR reactions. The thermal cycling conditions were composed of a step for UNG treatment at 25 • C for 5 min followed by an initial denaturation/activation step at 95 • C for 10 min, 45 cycles at 95 • C for 30 s, 60 • C for 60 s. The experiments were carried out in triplicate for each data point. The relative quantification of gene expression was determined using the 2−∆∆Ct method [53]. Normalization with GPI and ATP5B as housekeeping genes was done prior to statistical analysis for ileum and rumen samples, respectively. The housekeeping genes (GPI and ATP5B) were tested and found to be the most stable out of the four potentially stably expressed genes in their respective tissues using NormFinder. Statistical analysis (Student's t-test with homogeneous variances) was done using MIXED procedure and SAS software (Statistical Analysis System, release 9.4, 2002-2012, SAS Institute Inc., Cary, NC, USA).

RNA-Sequencing and Identification of Expressed mRNA and lncRNA Genes in the Gastrointestinal Tract of Calves
Next-generation RNA sequencing of 48 libraries from ileum and rumen tissues of calves generated a total of 2.5 billion reads (Supplementary Table S2). Following adaptor trimming and size selection, 2.4 billion reads (96%) with length >32 bp and having a phred score >30 were further processed. Out of this number, 91.84% were successfully mapped to the bovine genome (UMD.3.1). Out of successfully mapped reads, 86.17% uniquely mapped reads (Supplementary Table S2) were used in transcript construction.
A total of 24,616 expressed mRNAs were retained after removing transcripts that were not annotated to Ensemble bovine gene annotation release 84. After further filtering of genes with <5 normalized read counts and present in less than 10% of libraries, 15,905 and 15,628 genes in the ileum and rumen tissues, respectively were retained for further analyses (Supplementary Table S3).
After filtering and coding potential evaluation steps, lncRNA transcripts with ≥5 normalized counts and present in ≥10% of libraries resulting in 1568 (63 known and 1505 novel) and 4243 (88 known and 4155 novel) lncRNAs were considered expressed in ileum and rumen tissues, respectively (Supplementary Table S4a-d) and used in further analyses. Among them, 51 lncRNAs were common to both tissues (Supplementary Table S4e).

Differentially Expressed mRNAs and Enrichment Analyses
To identify differentially expressed genes, raw read counts of retained transcripts (mRNA or lncRNA) were imported into DESeq2 [48]. DESeq2 normalizes read counts by calculating a size factor for each sample to correct for library size and RNA composition bias. A total of 3310 and 4217 mRNAs were DE between D33 and D96 in Ileum and rumen tissues, respectively (Figures 1a and 2a Table S6c). The top 20 most significant DE mRNAs are shown in Table 1 and Figures 1a and 2a. Embigin (EMB) and radical S-adenosyl methionine domain containing 2 (RSAD2) were the most significant DE mRNAs for rumen and ileum tissues, respectively (Table 1).  A total of 459 biological process (BP) and 52 molecular function (MF) GO terms were significantly enriched for rumen DE mRNAs and the most enriched terms were sensory perception A total of 459 biological process (BP) and 52 molecular function (MF) GO terms were significantly enriched for rumen DE mRNAs and the most enriched terms were sensory perception (p = 1.3 × 10 −41 ) and G-protein coupled receptor activity (p = 2.4 × 10 −55 ), respectively (Figure 1b and Supplementary  Table S7a). Moreover, rumen DE mRNAs were also significantly enriched for 36 KEGG pathways and the most significantly enriched pathway was Ribosome (p = 1.00 × 10 −05 ) (Figure 1c and Figure S1 and Table S7b).
A total of 590 BP-GO and 110 MF-GO terms were significantly enriched for ileum DE mRNAs and the most enriched terms were neurological system process (p = 2.2 × 10 −43 ) and signaling receptor activity (p = 5 × 10 −50 ), respectively ( Figure 2b and Table S7c). Moreover, ileum DE mRNAs were significantly enriched for 86 KEGG pathways and the most significantly enriched pathway was Ubiquitin mediated proteolysis (p = 23.0 × 10 −12 ) (Figure 2c and Table S7d).
KEGG pathways and the most significantly enriched pathway was Ribosome (p = 1.00 × 10 ) (Figure 1c and Figure S1 and Table S7b).
A total of 590 BP-GO and 110 MF-GO terms were significantly enriched for ileum DE mRNAs and the most enriched terms were neurological system process (p = 2.2 × 10 −43 ) and signaling receptor activity (p = 5 × 10 −50 ), respectively (Figure 2b and Table S7c). Moreover, ileum DE mRNAs were significantly enriched for 86 KEGG pathways and the most significantly enriched pathway was Ubiquitin mediated proteolysis (p = 23.0 × 10 −12 ) (Figure 2c and Table S7d).

Identification and Enrichment Analyses of cis Target Genes of lncRNAs
To identify potential cis target genes of lncRNAs, a matrix of Pearson's correlation was computed using the expression data of 15,628 mRNAs and 4243 lncRNAs and between 15,905 mRNA and 1567 lncRNAs for rumen and ileum tissues, respectively. A total of 632 mRNAs (within 50 Kb surrounding regions of lncRNAs) were significantly correlated (p < 0.05) with 884 lncRNAs and were considered their potential cis targets in rumen tissues (Supplementary Table S8a). Eight genes (WWOX, TTLL7, STAG2, BMPR1A, NEDD4L, KIAAO319L, COA1 and SUDS3) were the potential cis targets of 27, 21,

Identification and Enrichment Analyses of cis Target Genes of lncRNAs
To identify potential cis target genes of lncRNAs, a matrix of Pearson's correlation was computed using the expression data of 15,628 mRNAs and 4243 lncRNAs and between 15,905 mRNA and 1567 lncRNAs for rumen and ileum tissues, respectively. A total of 632 mRNAs (within 50 Kb surrounding regions of lncRNAs) were significantly correlated (p < 0.05) with 884 lncRNAs and were considered their potential cis targets in rumen tissues (Supplementary Table S8a). Eight genes (WWOX, TTLL7, STAG2, BMPR1A, NEDD4L, KIAAO319L, COA1 and SUDS3) were the potential cis targets of 27,21,20,19,18,17,16,16 and 15 lncRNAs in rumen tissue, respectively while rXLOC_044298 and rXLOC_036836 potentially cis targeted 13 and 6 mRNAs, respectively (Supplementary Table S8a). GO enrichment analysis indicated that the 632 cis target genes of rumen lncRNAs were significantly enriched for 129 BP-GO and 13 MF-GO terms (Supplementary Table S9a, Figure 3) and 22 KEGG pathways (Supplementary Table S9b and Figure 4).

Differentially Expressed lncRNAs in Ileum and Rumen Tissues
A total of 14 and 525 lncRNAs were significantly DE between D33 and D96 in ileum and rumen tissues, respectively (Table S10a,b). The top significant DE lncRNAs are shown in Table 2. rXLOC_027852 and iXLOC_002882 were the most significant DE lncRNAs in rumen and ileum tissues, respectively. For rumen, 57 out of 525 DE lncRNAs had potential cis roles by correlating significantly with 77 mRNAs (Supplementary Table S10c). Most of the DE lncRNAs had only one potential cis target mRNA while rXLOC_025037 and rXLOC_006791 had the most potential cis target mRNAs (Figure 7). For ileum, 7 out of 14 DE lncRNAs had potential cis regulatory roles, targeting (significant correlation) five different mRNAs (Table S10d) and iXLOC_008931 potentially cis targeted three different mRNAs (Figure 7). Furthermore, one mRNA was the potential cis target of three lncRNAs (iXLOC_009320, iXLOC_009324 and iXLOC009321) in ileum tissue (Figure 7). Some of the potential cis target genes of DE lncRNAs were also significantly DE between D33 and D96 (Table 3 and Figure 7).

Real Time Quantitative PCR Confirmation of RNA-Seq Results
RNA-Seq results were mostly confirmed by qPCR analysis. For ileum samples, observed trend for all nine genes (four housekeeping and five DE genes) was similar to RNA-Seq results (Figure 8a). For rumen samples, the qPCR result was similar to RNA-Seq for four DE mRNAs, two lncRNAs and two housekeeping genes. However, one mRNA only showed a tendency towards significance (HERC6: p < 0.1) by the method of qPCR while one lncRNA (rXLOC_022071) was significantly DE by RNA-Seq but not by qPCR, although the fold change by both methods was the same (4.1-fold change, Supplementary Table S8b). Furthermore, two potential housekeeping genes PGK1 (p = 0.08) and PPIB (p = 0.06) in rumen tissue tended to be DE by the method of qPCR.

Real Time Quantitative PCR Confirmation of RNA-Seq Results
RNA-Seq results were mostly confirmed by qPCR analysis. For ileum samples, observed trend for all nine genes (four housekeeping and five DE genes) was similar to RNA-Seq results (Figure 8a). For rumen samples, the qPCR result was similar to RNA-Seq for four DE mRNAs, two lncRNAs and two housekeeping genes. However, one mRNA only showed a tendency towards significance (HERC6: p < 0.1) by the method of qPCR while one lncRNA (rXLOC_022071) was significantly DE by RNA-Seq but not by qPCR, although the fold change by both methods was the same (4.1-fold change, Supplementary Table S8b). Furthermore, two potential housekeeping genes PGK1 (p = 0.08) and PPIB (p = 0.06) in rumen tissue tended to be DE by the method of qPCR.

Transcriptome mRNA Transition from Pre-to Post-Weaning Period in Rumen and Ileum Tissues
The rumen and ileum are important organs for the digestion and absorption of nutrients in ruminants. The molecular mechanisms and key drivers underlying developmental changes in rumen and ileum in the early part of life have been revealed in several omics studies such as transcriptome analysis [54,55] and microbiome analysis [56]. In this study, we observed that 20.8% (33,107 genes)

Transcriptome mRNA Transition from Pre-to Post-Weaning Period in Rumen and Ileum Tissues
The rumen and ileum are important organs for the digestion and absorption of nutrients in ruminants. The molecular mechanisms and key drivers underlying developmental changes in rumen and ileum in the early part of life have been revealed in several omics studies such as transcriptome analysis [54,55] and microbiome analysis [56]. In this study, we observed that 20.8% (33,107 genes) and 27.0% (4217) bovine genes were significantly DE between the pre-weaning (D33) and post-weaning (D96) periods in ileum and rumen tissues, respectively (Figures 1 and 2, Supplementary Table S6). High variation of gene expression in rumen and ileum tissues between the pre-and post-weaning periods might reflect changes in anatomical size increases, metabolic changes as well as adaptation to different diets. Calves were fed with different diets (milk replacer and solid feed (starter diet and hay)) which require the activities/interactions of different groups of enzymes and pathways to digest and absorb them. As a consequence, distinct groups of genes and pathways are required in each period.
In rumen, several high DE genes are important for immune functions such as lymphocyte antigen 6 family member D (LY6D), mucin 1 (MUC1), embigin (EMB) and EYA transcriptional coactivator and phosphatase 2 (EYA2). LY6D protein initiates the first stage of B cell development [57] while MUC1 protein plays important roles in the first line defense against invading pathogens [58]. Meanwhile, EMB and EYA2 have important roles in development processes [59,60]. As expected, DE genes were enriched in many biological processes related to regulation of biological and cellular processes, which reflect the major developmental changes in the physiology and functions of the rumen from pre-to post-weaning stages ( Figure 1). Notably, top BP-and MF-GO terms are involved in G-protein-coupled receptors, which are important components of cell signaling and also vital processes required for the normal growth and development of cells [61]. Interestingly, top enriched MF-GO terms (e.g., poly(A) RNA binding, RNA binding, transcription factor binding and nucleic acid binding, etc.) for DE rumen genes involved in RNA processing reflects the many transcriptional activities that take place during growth and development. The most important pathway enriched for rumen DE genes was ribosomes, an important pathway for the regulation of protein biosynthesis. Other interesting pathways have roles in cellular signaling such as thyroid hormone, B cell receptor and Notch signaling pathways. The regulation of thyroid hormone is crucial for animals to respond to changes in nutrient availability and requirements and, to homeorhetic changes during different physiological stages [62]. Therefore, it is not surprising that this pathway was the most significantly enriched pathway in this study. Meanwhile, the enrichment of B cell receptor signaling pathway might reflect changes in immune activities in the rumen during the early period of growth, while the enrichment of Notch signaling pathways might reflect changes in cellular activities in rumen from the pre-weaning to the post-weaning stages.
In ileum, radical S-adenosyl methionine domain containing 2 (RSAD2) gene was the most significant DE gene between D33 and D96 (p.BH) = 6.05 × 10 −11 ). This gene plays a key role in the innate immune response system and has a wide variety of antiviral activities [63]. Another most significant DE gene, 2 -5 -oligoadenylate synthetase 2 (OAS2) (p.BH = 2.34 × 10 −08 ), functions as a regulator in the innate immune response process [64,65]. Similar to rumen DE genes, many ileum DE genes were enriched in BP-GO terms involved in cellular signaling and system development while MF-GO terms are involved in RNA processing and cell cycle regulation (Figure 2, Table S7c). Notably, several BP-GO terms (negative regulation of type I interferon production, interferon-beta production and positive regulation by symbiont of host I-κB kinase/nuclear factor-κB (NF-κB) cascade) and KEGG pathways (T cell receptor signaling pathway, B cell receptor signaling pathway and tumor necrosis factor (TNF) signaling pathway) related to immune activities were significantly enriched for ileum DE genes. Previously, Malmuthuge et al. [66] examined the expression of immune related genes (Toll-like receptors (TLR), peptidoglycan recognition protein 1 (PGLYRP1) and antimicrobial peptides (β-defensin)) in GIT tissues from pre-(three weeks of age) and post-weaned (six weeks of age) calves and reported a down regulation of most TLR genes after weaning, coinciding with a drastic change in the microbiome, as well as the morphological and functional changes in the GIT. These results may reflect a developmental shift toward that of a mature ruminant GIT to prevent unnecessary inflammatory responses to a shifting intestinal microbiome [66]. However, we did not observe significant differential expression of TLR genes between the pre-and post-weaning periods in this study.
In fact, the most important pathway enriched for DE genes in ileum was ubiquitin mediated proteolysis. Protein ubiquitination mainly functions as a signal for 26S proteasome dependent protein degradation and therefore plays an important role in eukaryotic cellular processes. Another interesting enriched pathway was insulin signaling pathway which was not enriched for rumen DE genes. The insulin signaling pathway is crucial for carbohydrate metabolism. Therefore, the expression of genes related to this process was changed more in the Ileum than in the rumen. Using functional prediction from microbiome data, Meale et al. [67] reported that genes related to carbohydrate metabolism decreased post-weaning suggesting a shift in nutrient metabolism away from the lower GIT to ruminal fermentation. Some other significantly enriched pathways were related to growth control and disease (e.g., mechanistic target of rapamycin (mTOR) signaling pathway), regulation of cell cycle (e.g., cell cycle pathway) and cell matrix adhesions (e.g., focal adhesion pathway).

lncRNAs in Rumen and Ileum Tissues and Predicted Functions
In this study, we identified 1567 (62 known and 1505 novel) and 4243 (88 known and 4155 novel) lncRNAs in ileum and rumen tissues, respectively (Supplementary Table S4). The differences in the number of identified lncRNAs in the rumen and ileum tissues might reflect their different roles in each tissue since lncRNAs are known to be tissue specific [4,20,21]. Similar to lncRNA studies in other bovine tissues [21,[29][30][31]33] (the majority of identified lncRNAs were 200-999 bp long (44.70%) and composed of one transcript (95.85%). Moreover, in the same trend with other studies [4,21,28], the majority of identified lncRNA transcripts are located in the intergenic regions.
It is well documented that lncRNAs may act in cis or trans to regulate the activities of genes [15,[68][69][70][71][72]. Due to current lack of information, the functions of bovine lncRNAs can be inferred either by their physical position to protein coding genes such as their potential cis target genes (genes neighboring lncRNAs) or by their co-expression relationship with mRNAs (guilt by association methods) [73][74][75][76]) or determined experimentally. To infer the function of lncRNAs, we identified their cis target genes (genes located within 50 Kb flanking regions of lncRNAs), especially those having significant co-expression correlation (p.BH < 0.05) with lncRNAs. The combination of two filtering options (distance and strength of correlation) strengthened the confidence of potential cis target gene prediction for these lncRNAs but might lead to the loss of information since there is no clear definition of the distance from lncRNAs that should qualify a gene as its cis target. Moreover, it should be noted that, the co-expression (guilt by association) analysis might not completely reflect the causal relationship between lncRNAs and mRNAs.
Nevertheless, we identified 632 potential cis target genes for 884 lncRNAs (Supplementary Table S8a) in the rumen. A most interesting potential cis target gene was WWOX, predicted as potential cis target gene to 27 different lncRNAs. This gene encodes a member of the short-chain dehydrogenases/reductases (SDR) protein family and is known as a tumor suppressor [77]. Furthermore, WWOX has been shown to be regulated by lncRNA rP11-190D6.2 [78]. Since WWOX was not significantly DE between pre-and post-weaning (p.BH = 0.13) periods, we suspect that it performs similar functions in both periods. Another interesting gene was bone morphogenetic protein receptor type 1A (BMPR1A), predicted as cis target of 13 lncRNAs (rXLOC_042928, rXLOC_042929, rXLOC_042930, rXLOC_042932, rXLOC_042934, rXLOC_042935, rXLOC_042936, rXLOC_042938, rXLOC_042939, rXLOC_042940, rXLOC_042942, rXLOC_042945 and rXLOC_042946) is important in Juvenile Polyposis Syndrome (JPS) in humans [79]. JPS is characterized by predisposition to hamartomatous polyps in the GIT [80]. Given the important function of BMPR1A gene in the GIT health of humans, further studies are required to establish it roles in the bovine GIT. rXLOC_044298 and rXLOC_036836 potentially cis targeted 13 and 6 mRNAs, respectively (Supplementary Table S8a) but none of them was significantly DE between the two periods. Therefore, further functional studies to characterise their roles are needed. rXLOC_036836 (NONBTAG014380.2) is a known lncRNAs located on bovine chromosome 5 (from 104,262,200 to 104,262,521 bp) but no report on it function in the rumen exist. Enrichment results indicated that a high number of potential cis target genes are involved in posttranscriptional gene silencing by RNA or microRNA processing. Although the interaction between lncRNAs and miRNAs is well documented [81,82], no information on such interaction is available in the rumen. The enrichment of potential cis target genes of lncRNAs suggest potential roles of lncRNAs in the regulation of cell cycle, metabolic process, cellular metabolic processes and cell macromolecule catabolic process (Figure 3) which are essential processes for normal development. Moreover, rumen lncRNAs play potential roles in the regulation of pathways involving energy sources (oxidative phosphorylation pathway), cell signaling (p53 signalling or ErbB signaling) and immunity (T cell receptor signaling, Fc-γ-mediated phagocytosis and NF-κB signaling) (Figure 4). The importance of these pathways in the development and functions of the rumen have been reported by several authors [83][84][85]. Perhaps, the most interesting results for rumen lncRNAs is the fact that several DE lncRNAs also potentially cis targeted mRNAs that were DE between pre-and post-weaning periods. A notable rumen DE lncRNA is rXLOC_025037, which also potentially cis targeted four DE genes (WDR97, MAF1, EXOSC4 and GPAA1) (Figure 7). Glycosylphosphatidylinositol anchor attachment 1 (GPAA1) gene serve an important function in linking proteins to the cell surface membrane while exosome component 4 (EXOSC4) is related to cyclin-dependent kinases (CDK)-mediated phosphorylation and removal of cell division cycle 6 (CDC6) and deadenylation-dependent mRNA decay [86], however it is not clear how these genes function in the rumen.
In ileum, 610 lncRNAs potentially cis-targeted 439 mRNAs (Supplementary Table S8b) and BMPR1A and ST3 β-galactoside α-2,3-sialyltransferase 1 (ST3GAL1) were amongst genes potentially cis-targeted by the highest number of lncRNAs. The function of BMPR1A has been described above, while ST3GAL1 is an important gene in cancer [87] but little is known about its functions in the ileum. Microtubule cytoskeleton organization and endomembrane system organization were the most significantly enriched GO terms for potential cis target genes of ileum lncRNAs ( Figure 5, Supplementary Table S9c). These GO terms are essential to many cellular processes. Immune related GO terms such as modulation by virus of host morphology, physiology, lymphatic endothelial cell differentiation and regulation of viral release from host cell were also significantly enriched ( Figure 5, Supplementary Table S9c) suggesting that ileum lncRNAs might have roles in the regulation of immunity pathways. Moreover, we also observed that inflammatory mediator regulation of transient receptor potential (TRP) channels pathway was significantly enriched for ileum lncRNAs potential cis target genes (ADCY7, MAP2K6, MAPK8 and PLCG1). The mammalian TRP superfamily comprises 28 TRP cation channels [88] involved in the regulation of gastrointestinal motility and absorption [89]. Similar to rumen, ileum lncRNAs might also regulate genes in oxidative phosphorylation pathway and lipid metabolism related pathways (glycosphingolipid biosynthesis and phosphatidylinositol signaling system). Interestingly, three DE lncRNAs also potentially cis regulated the G-patch domain-containing 2 (GPATCH2) gene with potential roles in spermatogenesis [90] and tumor growth during breast cancer [91]. However, the function of this gene in ileum is still unknown. Another notable gene is cysteine-and glycine-rich protein 1 (CSRP1), a member of the cysteine-rich protein (CSRP) family involved in regulatory processes important for development and cellular differentiation. Variants of CSRP1 have been associated with growth and carcass traits in cattle [92].

Conclusions
To the best of our knowledge, this is the first study to characterize the lncRNA expression in calf's ileum and rumen tissues. This study has shown the potential molecular mechanisms and pathways involved in the early development of calf GIT. The enrichment of DE mRNAs confirmed the major drivers for the switch from pre-to post-weaning involved in the regulation of growth, development and immune function. Relevant genes for these processes such as LY6D, MUC1, EMB and EYA2 in rumen and RSAD2 and OAS2 in ileum tissues were identified. This study also characterized lncRNAs that may be involved in developmental processes of calves GIT and identified potential lncRNAs and the potential cis target genes that might be important for the regulation of the development of the GIT. Therefore, this study has set a baseline for the further exploration of the roles of lncRNAs in calf's GIT.

Supplementary Materials:
The following are available online at www.mdpi.com/2073-4425/9/3/142/s1. Figure S1: KEGG pathways enriched for differentially expressed mRNAs in rumen tissue, Figure S2: KEGG pathways enriched for differentially expressed mRNAs in ileum tissue, Table S1: Genes and primer sequences used in qPCR validation of RNA-Sequencing data, Table S2: RNA-Seq read mapping statistics, Table S3: Genes (mRNA) in (a) rumen and (b) ileum tissues retained for further analyses, Table S4: Identified lncRNA genes in ileum and rumen tissues of calves. (a) Ileum known lncRNAs; (b) Ileum novel lncRNAs; (c) Rumen known lncRNAs; (d) Rumen novel lncRNAs; (e) LncRNAs common to ileum and rumen tissues, Table S5 Table S8: Potential cis target genes of lncRNAs in (a) rumen and (b) ileum tissues and (c) differentially expressed lncRNAs and their differentially expressed mRNAs, Table S9: (a) GO terms enriched for rumen lncRNA potential cis target genes; (b) KEGG pathways enriched for rumen lncRNA potential cis target genes; (c) GO terms enriched for ileum lncRNA potential cis target genes; (d) KEGG pathways enriched for ileum lncRNA potential cis target genes, Table S10: Differentially expressed lncRNAs between day 33 and day 96 in rumen (a) and ileum (b) and potential cis target mRNAs for DE lncRNAs in (c) rumen and (d) ileum tissues.

Availability of Sequence Data:
The sequence data has been submitted to Gene Expression Omnibus database (BioProject PRJNA400067) and it is available through this link: https://www.ncbi.nlm.nih.gov/bioproject/400067.