Identification and Characterization of MicroRNAs in the Goat (Capra hircus) Rumen during Embryonic Development

The rumen is an important digestive organ in ruminants. Numerous regulatory factors including microRNAs (miRNAs) are involved in embryonic organ development. In the present study, miRNAs expressed in the rumens of goats (Capra hircus) and their potential roles in the pathways involved in rumen development were identified using high-throughput sequencing. Histological morphology revealed a distinct difference in each layer of rumen during the period from embryonic day 60 (E60) to embryonic day 135 (E135). We determined the expression profiles of miRNAs in the goat rumen, and identified 423 known miRNAs and 559 potentially novel miRNAs in the E60 and E135 embryonic rumen, respectively. Bioinformatics analysis annotated the 42 differentially expressed miRNAs and the top 10 most highly expressed miRNAs of the two libraries to 48 and 38 gene ontology categories, as well as to 168 and 71 Kyoto Encyclopedia of Genes and Genomes pathways, respectively. The expression patterns of eight randomly selected miRNAs were validated by stem-loop quantitative reverse transcription PCR, suggesting that the sequencing data were reliable. We profiled the genome-wide expression of rumen-expressed miRNAs at different prenatal stages of rumen tissues, revealing that a subset of miRNAs might play important roles in the formation of the rumen layers. Taken together, these findings will aid the investigation of dominant rumen-related miRNA sets and help understand the genetic control of rumen development in goats.


INTRODUCTION
The rumen, a specialist digestive organ of ruminants, plays a vital role in digestion metabolism. Its wall can be divided into four layers: the mucosa, submucosa, muscularis, and tunica adventitia. In goats (Capra hircus), rumen papillae appear to be separated by clear boundaries at embryonic day 64 (E64); by embryonic day 150 (E150), rumen papillae are fully developed and show signs of surface keratinization (Franco et al., 1992). To date, studies have focused on the effect of dietary composition and nutrition on rumen development (Coverdale et al., 2004;Gupta et al., 2016;Suarez-Mena et al., 2016). Naeem et al. (2012) also evaluated the expression of genes involved in rumen development, while several microRNAs (miRNAs) associated with development of the digestive tract have been reported (Hino et al., 2008;McKenna et al., 2010;Nguyen et al., 2010;Dalmasso et al., 2011;Kim et al., 2011;Ye et al., 2011;Chen et al., 2013). However, the role of miRNAs in the goat rumen during embryonic development has not been described.
MicroRNAs were originally identified as negative regulators of their target genes that binding to the sequences of 3 -untranslated region (3 -UTR) (Lim et al., 2005). They may also induce gene expression by targeting the promoter sequences of target genes (Place et al., 2008). Although miRNAs are encoded by only a small number of genes, they may regulate the expression of many mRNAs Lewis et al., 2005). Increasing evidence has suggested that miRNAs are implicated in embryonic development , cell differentiation (Ivey and Srivastava, 2010), proliferation (Brennecke et al., 2003), apoptosis (Cimmino et al., 2005), and tumorigenesis (Yamanaka et al., 2012;Han et al., 2015).
Little genetic information is available about the development of the rumen in goats, and differences between miRNA profiles during the formation of rumen layers are also unknown. Therefore, in the present study, we used RNA sequencing technology to understand the complex regulation of rumen development at gestational day 60 (before rumen papillae formation) and gestational day 135 (after rumen papillae formation). Systematical analyses were carried out of known miRNAs, potentially novel miRNAs, differentially expressed miRNAs, and their target genes. We also observed the morphology of rumen tissues at different development stages (E60 and E135) to provide an insight into the mechanisms of miRNA regulation of rumen development.

Ethics Statement
All animal experiments were performed in strict accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and given prior approval by the Institutional Animal Care and Use Committee at the College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China under permit No. DKY-B20110807. All animals were sacrificed humanely to minimize suffering.

Rumen Collection and Total RNA Isolation
Pregnant does were provided by the Station of the Jianyang Da'er Goat Breeding Center (Sichuan, China). Six pregnant does at two embryonic stages (E60 and E135) were anesthetized by Zoletil (4 mg/kg body weight) for cesarean section. Rumen tissues were collected from the embryos (n = 3 per stage). The six rumen tissues were immediately frozen in liquid nitrogen, and then stored at −80 • C until required for RNA isolation. Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The purity, concentration, and quality of the RNA were determined by a Nanodrop 2000 spectrophotometer (Thermo Scientific Nanodrop, Wilmington, DE, United States), Qubit 3.0 Fluorometer (Life Technologies, Camarillo, CA, United States), and Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States), respectively.

Hematoxylin and Eosin Staining
The fixed rumen samples were dehydrated and embedded in paraffin. Sections (5-µm) were cut using a microtome (LeicaRM2016, Leica Microsystems, Germany) and placed on glass slides. The sections were stained with hematoxylin and eosin (H&E) and viewed using an Olympus BX51 microscope (Olympus, Tokyo, Japan) equipped with a Nikon DS-Fi1 camera (Nikon, Tokyo, Japan).

Small RNA Library Construction and Sequencing
Approximately 1.5 µg of total RNA was used to generate small RNA libraries with TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, United States) according to the manufacturer's protocol. Briefly, small RNAs were ligated to 3 and 5 RNA adapters. The ligated small RNAs were reverse transcribed and amplified by PCR. The PCR products were then purified by denaturing polyacrylamide gel electrophoresis (PAGE). After examining the concentration and insert size by a Qubit 3.0 Fluorometer (Life Technologies, Camarillo, CA, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States), the libraries were subjected to 50 bp single-end reads sequencing on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, United States).

Sequencing Data Analysis and Identification of miRNAs
Raw reads in fastq format were filtered out the low-quality reads, reads containing poly-N (>10% unknown bases), and reads containing adapter using in-house scripts (fastq v1.1 1 ) written in Perl. Then reads were trimmed and cleaned by removing the sequences smaller than 18 nt or longer than 30 nt. In addition, Q30 of the clean data were also calculated. The clean reads were aligned with Silva, GtRNAdb, Rfam, and Repbase databases, respectively, to filter rRNA, tRNA, snRNA, snoRNA, other ncRNA, and repeat sequences using the Bowtie v.1.1.0 software (Langmead et al., 2009). The derived reads were used to detect known goat miRNAs by aligned to the goat miRNA precursors present in the miRBase v21 2 as described in our previous study (Guo et al., 2016a). Clear reads were also mapped to the goat reference genome (CHIR_1.0, NCBI) using the Bowtie aligner. The mapped reads were further denoted as known miRNAs under the following criteria: (i) aligned to the mature known miRNA precursors of other animal species in miRBase, (ii) raw count for known miRNAs was quantified and normalized using transcripts per million (TPM) mapped reads. Only those miRNAs, an overall expression level greater than one RPM in six libraries, were denoted as true known miRNAs. The unannotated reads which mismatched to any public databases or goat genome were further analyzed to identify novel miRNAs using Mirdeep2 v2.0.0.7 (Friedlander et al., 2012), with the default options (e.g., minimal free energy < −18 kcal/mol) as described in our previous study (Guo et al., 2016a).

Differential Expression Analysis of miRNAs
To assess the differentially expressed miRNAs, the DESeq2 package in the statistics software R (Version 3.0.1) was used with standard parameters (Love et al., 2014). In order to avoid false positive results, the two criteria: only |log 2 (Fold change)| ≥ 1 and the false-discovery rate (FDR) adjusted p-value ≤ 0.01 were used to determine the differentially expressed miRNAs.

Target Prediction and Functional Annotation of Target Genes
TargetScan (Lewis et al., 2003), RNAhybrid (Kruger and Rehmsmeier, 2006), and miRanda (Betel et al., 2010) prediction tools were used to predict the target genes of differentially expressed miRNAs. To filter out the false positives, the predicted target genes were retained only when they were identified by all the three programs.

Validation of Differentially Expressed miRNAs
To validate the sequencing data, eight randomly selected miRNAs (including six known and two novel miRNAs) were validated by stem-loop quantitative real-time reverse transcription PCR (qRT-PCR). Total RNA (5 ng) was reverse transcribed into cDNA using the Mir-X TM miRNA First-Strand Synthesis Kit (TaKaRa, Dalian, China) according to the user's manual. To detect the expression patterns of target genes of differentially expressed miRNAs, cDNA was synthesized using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China). The cDNA products were diluted 1:5 (v/v) with EASY Dilution (TaKaRa, Dalian, China) and stored at −20 • C. Primers used in this study are shown in Supplementary Table S8.
Quantitative real-time reverse transcription PCR was performed using a CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, Singapore) in a volume of 10 µL including 5 µL of SYBR Green Real Time PCR Master Mix (TaKaRa, Dalian, China), 0.2 µL each of forward and reverse primer, 3.8 µL of RNase-Free dH 2 O, and 0.8 µL of cDNA. U6 snRNA was used as the reference for miRNA and target gene candidates (Guo et al., 2016a). All reactions were performed in triplicate, and the relative expression level was normalized by 2 − C T method.

Histological Observation of the Embryonic Goat Rumen
To observe the differences between E60 and E135 goat rumens, H&E-stained sections were prepared. As shown in Figure 1, each layer of E60 rumen was thinner than that those of the E135 rumen. At E135, rumen papillae were more developed than at E60, appearing as evaginations of the stratum basale toward the stratum granulosum. Moreover, cells in the lamina propria were more closely packed, the diameters of muscle fibers were increased, and the spaces between muscle fibers were wider. Finally, more numerous blood vessels were found in the tunica muscularis.

Overview of Sequencing Data
To obtain a transcriptome of the caprine rumen at different developmental stages, two libraries were constructed from goats at E60 and E135. An average of 18,810,062 and 21,162,680 raw reads was obtained from E60 and E135, respectively ( Table 1).
After removing low-quality sequences and adapters, an average of 15,500,932 and 19,252,745 clean reads, respectively, was retained ( Table 1). The rRNA, scRNA, snRNA, snoRNA, tRNA, and repeat sequences were filtered from clean reads to obtain an average of 13,630,047 and 17,583,220 reads, respectively, which were mapped to the goat genome (CHIR_1.0, NCBI) using the miRDeep2 program (Supplementary Table S1). The reads between 18 and 30 nt were BLASTed against the miRBase, which identified 423 known miRNAs and 559 novel miRNAs (Supplementary Table S2).

Expression Analysis of miRNAs in the Embryonic Goat Rumen
As shown in Figure 2A, most miRNAs were 22 nt in length, followed by 21 and 23 nt, indicating good sequencing quality. Expression density distributions were analyzed according to the normalized reads. As shown in Figure 2B, the expression levels were highly consistent among the six tested samples and between the different embryonic periods (Supplementary Figure S1). Differential expression analysis found that 42 miRNAs were differentially expressed, of which 22 miRNAs were up-regulated and 20 were down-regulated: 27 of the miRNAs were known while 15 were novel (Table 2). However, sample size is an important point for the determination of DE miRNAs, less abundant miRNAs might be hidden in a small sample size. We also focused on miRNAs with high TPM in each library. The reads of the top 10 most highly expressed miRNAs accounted for 50.6 and 63.9% of the total average reads, respectively, in each library. A total of 11 miRNAs are listed in the top 10 of the two libraries, of which 9 (except miR-7-5p and miR-145-5p) are found in both libraries ( Table 3).

GO Enrichment Analysis and KEGG Pathway Analysis of Target Genes
Target gene prediction found that 223 genes were predicted to be targeted by 18 differentially expressed miRNAs, while 65 genes were predicted to be targeted by 6 of the top 10 most highly expressed miRNAs (Supplementary Tables S3, S4).
Supplementary Table S5 shows the annotation of target genes for all miRNAs.
Gene ontology enrichment analysis of these putative target genes revealed that 175 and 54 target genes of differentially expressed and the top 10 most highly expressed miRNAs, respectively, were annotated. As shown in Figure 3, the most enriched GO categories of differentially expressed miRNAs were consistent with those of the top 10 most highly expressed miRNAs. The analysis of cellular component showed that most genes were involved in cell and cell part categories. The analysis of molecular function showed that most genes were clustered into binding and catalytic activities. Moreover, the analysis of biological process showed that most genes were related to single organism and cellular processes.
Kyoto Encyclopedia of Genes and Genomes pathway annotation showed that 168 and 71 biological functions were annotated by target genes of differentially expressed and the top 10 most highly expressed miRNAs, respectively (Supplementary  Tables S6, S7). The enriched pathways for target genes of differentially expressed miRNAs mainly included the Ras signaling pathway, the phosphoinositide 3-kinase (PI3K)-Akt signaling pathway, focal adhesion, and regulation of the actin cytoskeleton. KEGG pathway analysis based on predicted targets revealed that the top 10 most highly expressed miRNAs were involved in the mitogen-activated protein kinase (MAPK) signaling pathway, the Ras signaling pathway, the PI3K-Akt signaling pathway, and olfactory transduction.

Quantitative Real-Time-PCR Validation of the Differentially Expressed miRNAs
Eight miRNAs, including six known and two novel miRNAs, were randomly selected to validate the reliability of the sequencing data. As shown in Figure 4, relative expression changes of miRNAs were represented using log 2 (E135/E60). All selected miRNAs showed expression changes that were consistent with the deep sequencing results, suggesting that the sequencing data were credible.

DISCUSSION
Histological analysis revealed that the goat E60 rumen is composed of four layers, and that the rumen mainly undergoes Frontiers in Genetics | www.frontiersin.org a process of tissue growth during the period from E60 to E135. Rudimentary pillars have previously been observed in the E46 rumen, whose wall lacked a complete structure (Garcia et al., 2012). This indicates that samples earlier than E46 should be used to identify miRNAs with important roles in rumen development.
Embryonic rumen development involves the formation and differentiation of the rumen wall, and differentiation of the epithelial layer. The molecular mechanism for this complex biological process has not been reported yet. However, numerous studies have focused on the development and regulatory mechanisms of intestinal epithelial cells, and Wnt and Notch signaling pathways were identified as determining factors for cell fate (Nakamura et al., 2007). In this study, differentially expressed miRNAs such as miR-182 (Kim et al., 2012), miR-23a (Wang et al., 2012), andmiR-195-3p (Chen et al., 2011) were reported to mainly take part in the regulation of cell proliferation, differentiation, and apoptosis. Moreover, miR-34a was shown to promote renal senescence and regulate neural stem cell differentiation (Aranha et al., 2011;Bai et al., 2011), while miR-592 regulates neuronal apoptosis through p75 NTR (Irmady et al., 2014). These findings suggest that differentially expressed miRNAs are involved in cell differentiation and apoptosis, mitochondrial metabolism and neuronal development of the goat rumen. Moreover, the top 10 most highly expressed miRNAs, such as miR-26a-5p (Guo et al., 2016b) and miR-27b-3p (Tao et al., 2015), were also implicated in cell proliferation. Interestingly, let-7i-5p, let-7f-5p, and miR-21-5p are highly expressed in the porcine ovary and testis, suggesting that they are associated with secretory function (Li et al., 2011). miR-143 and miR-145 are co-expressed miRNAs that are not expressed in colonic epithelial tissue (Kent et al., 2014). Deep sequencing results revealed that they were highly expressed in the goat rumen, although the expression of miR-145 was significantly lower than that of miR-143 in both E60 and E135, which may be caused by post-transcriptional regulation. The expression of the conserved miRNAs in goat rumen revealed that the top 10 abundant miRNAs might act important roles in the   development or physiology of rumen layers. These miRNAs could be studied as the preferred candidate miRNAs for regulation analysis of rumen development in ruminants. However, it is also crucial to point that DE miRNAs associated with rumen development could be limited by the small sample size in our study. Therefore, the evaluation of DE miRNAs in goat rumen with larger sample size might produce more reliable and compatible results. Functional enrichment analysis revealed that the differentially expressed and top 10 most highly expressed miRNAs shared similar functions (Supplementary Table S5 and Figure 3). GO annotation analysis revealed that a large proportion of the genes of miRNA targets in goat rumen were involved in the (extra)cellular structures and cytoskeleton, which had binding and catalytic functions. Consistently, no specifically expressed miRNAs were found between E60 and E135. The most enriched pathways, including MAPK, Ras, and PI3K-Akt signaling pathways, are implicated in gastric cell proliferation, intestinal epithelial cell growth, and intestinal epithelial cell proliferation, respectively (Osaki et al., 2011;Waseem et al., 2014;Depeille et al., 2015). In addition, KEGG analysis also showed that approximately 24% pathways of the top 10 miRNAs were committed to internal secretion, cancer, leukemia, and virus infection diseases (Supplementary Table S7). These results indicated that some miRNAs might be involved in cell proliferation, apoptosis, and phagocytosis. Further research is required to demonstrate the relationship between miRNAs and target mRNAs transcripts and their pathways in the development of the goat rumen.

CONCLUSION
This study identified miRNAs of the E60 and E135 goat rumen and evaluated expression of these miRNAs at the two time points. A total of 423 known miRNAs and 559 novel miRNAs were identified. Analysis of 42 differentially expressed and the 11 top most highly expressed miRNAs revealed that most of their target genes are involved in cellular activities such as proliferation, differentiation, and apoptosis, which is consistent with histological observations. To identify miRNAs that play essential roles in early rumen development, samples before E46 should be examined. However, our current results provide useful information for further research, which will contribute to the elucidation of the molecular mechanisms involved in rumen development.

AUTHOR CONTRIBUTIONS
TZ conceived and designed the experiments. JH and PX performed the experiments. JH and SZ analyzed the data. LW, JG, LL, and HZ participated in sample collection. HZ, LN, and TZ contributed reagents/materials. JH and TZ wrote the manuscript.

FUNDING
This work was partially supported by the grants from National Natural Science Foundation of China (31501936) and the Chinese Domestic Animal Germplasm Resources Infrastructure.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2017.00163/full#supplementary-material FIGURE S1 | Box plot of the expression levels among the six tested rumen samples.