Evaluation of key miRNAs during early pregnancy in Kazakh horse using RNA sequencing

Ethics statement

All experimental procedures involving animals were approved by the Animal Care and Use Committee of Xinjiang Agricultural University, Urumqi, Xinjiang, China (animal protocol number: 2017008).

Animals and sample collection

All experiments were performed on eight healthy Kazakh horses. Eight lactating horses, parity four, were selected from Fuyun County. Lactating horses were allocated to the H group (milk yield >  7 kg) or the L group (milk yield <  3 kg) according to their milk yield (Table S1). A 300 mL sample of milk was collected from each horse, quickly frozen on liquid nitrogen, and brought back to the laboratory where the samples were kept at −80 °C until testing.

Milk fat extraction

A total of 50 mL (2700 g) samples of milk were centrifuged at 10 °C for 10 min; the upper layer of milk fat was removed and trizol was added to a 2:8 ratio. The mixture was thoroughly shaken and then frozen at −80 °C.

RNA extraction

Total RNA extraction was performed using the miRNeasy Mini Kit (217004, QIAGEN, Germany) according to the standard operating procedures provided by the manufacturer. The Qubit^®2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and the Nanodrop One spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA) were used to assess the concentration and quality of RNA. The integrity of the total RNA was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., USA) and samples with RNA integrity number (RIN) values above 7.0 were used for sequencing.

Small RNA library construction and sequencing

A Qiagen™ Small RNA Library Kit (Qiagen) was used to prepare small RNA-seq libraries according to the manufacturer’s protocol. The Agilent 2100 Bioanalyzer and an Agilent High Sensitivity DNA chip (Agilent Technologies) were used to assess the miRNA-seq libraries. The final product was a library of approximately 173–181 bp. The sequence tags were obtained via Illumina Nova Seq 6000 by the Shanghai Sinomics Corporation.

Small RNA-seq data analysis

MiRDeep software was used to predict the known and novel miRNAs (Friedländer et al., 2008). The reads were trimmed with Cutadapt v 1.9.1 (Martin, 2011). Reads shorter than 15 bp or longer than 40 bp were discarded. Quality control of the reads was performed using seqtk (1.0-r82-dirty) and R (3.4.2). The ends of the filtered reads were further trimmed to remove low-quality bases (Phred score ≤ 20). The clean reads were mapped to the horse genome (EquCab 2.0) using bowtie (1.0.0). Gene abundance was expressed as counts of exon model per million mapped reads (CPM). R (3.4.2) was used with a self-written script and other software was used with default parameters.

Identification of differentially expressed miRNAs (DEMs) and prediction analysis of DEMs

The DESeq R package (v. 1.8.3) was used to identify DEMs based on miRNA expression values by pairwise comparison to understand the functional role of DEMs in horse milk of the two horse groups (Wang et al., 2010). The cut-off criteria for DEMs were —log2(Fold change) — ≥ 1 and P-value ≤ 0.05. The target genes of these miRNAs were predicted using the miRanda (Betel et al., 2008) analytical tools. MiRNA target gene searching was performed using the MultimiR package, and the search consisted of predicted and validated target genes.

Enrichment analyses for DEM targets

GO annotations of differentially expressed miRNA target genes were performed using the gene ontology (GO) database (http://geneontology.org/); pathway analysis of differentially expressed genes was performed using the KEGG (Kyoto Encyclopedia of genes and genomes) database. Both the GO annotation and the KEGG pathway analysis had a P value of less than 0.05 as a significant enrichment criterion.

Validation of miRNA expression by real-time PCR (RT-PCR)

Libraries of cDNA were created using the HiScript II Q RT SuperMix. A total of 500 ng RNA was used per reaction, following manufacturer’s instructions. We used 10 µL 2 ×miRNA RT Reaction Buffer, 2 µL miRNA RT Mix with nuclease-free water added to bring the final volume to 20 µL. An additional aliquot was made with water in place of the RNA for a negative control. Tubes were held at 42 °C for 60 min followed by 3 min at 95 °C, then cooled to 4 °C and held at this temperature until use.

Reactions were immediately prepared for qPCR using the miRcute Plus miRNA qPCR Kit (Qiagen). Each reaction consisted of 5 μL 2 ×miRcute Plus miRNA Premix, 1 µL Amplification Primer (2 µM), 1 µL Reverse Primer (2 µM), 2 µL nuclease-free H₂O, and 1 µL cDNA. All reactions were performed in duplicate. The RNA-free reverse transcription reaction and a cDNA-free PCR reaction were used as negative controls.

Thermal cycling took place on a 7900 real-time qPCR system (ABI, USA), and consisted of an initial incubation at 50 °C for two min, 95 °C for 10 min, 95 °C for 15 s followed by 40 cycles at 60 °C for 1 min, 95 °C for 15 s, 60 °C for 15 s and 95 °C for 15 s. Porcine U6 snRNA (horse) and eca-miR-26a were used as an endogenous reference. Data were normalized by ΔCt = Ct (target) −Ct (mean of sample). The primer sequence of 5 miRNAs are shown in Table S2.

Overview of miRNA sequencing

A total of 40.3 million and 28.1 million raw reads were obtained from the milk of the H group and L group, respectively. 28,055,021 and 20,692,008 clean reads from the H and L horse milk libraries, respectively, were selected for further analysis after the reads were assessed for quality (Table 1). The QC results are shown in Table S3. Most of these miRNAs ranged from 20-24 nt, and the read-length distributions of the H group and L groups are shown in Fig. 1.

We identified 341 miRNAs in the H group and 333 miRNAs in the L group. The results of the analysis demonstrated that there were 83 significantly differentially expressed miRNAs, which were divided into two categories: known miRNAs comprising 32 miRNAs (27 up-regulated, 5 down-regulated) and novel miRNAs comprising 51 miRNAs (9 up-regulated, 42 down-regulated) (Table S4).

The effect of milk on the miRNA was also evidenced by a cluster analysis and a principle component analysis (PCA) (Fig. 2). Samples were clustered according to H and L groups (Fig. 2A). The first principal component (PC1) accounted for 41% of the total variance, and the second components (PC2) accounted for 19% of the total variance (Fig. 2B).

Target gene prediction and gene functional annotation

Target prediction is an important way to determine the functions of miRNAs. We predicted a total of 2,415 target genes among the differentially expressed miRNAs. GO gene function analysis showed the enrichment of some terms related to breast development, nervous system development, and kinase activity (Fig. 3). Results of the KEGG analysis showed that these target genes were involved mainly in the insulin signaling pathway, insulin secretion, breast cancer, oocyte division, and the thyroid hormone signaling pathway (Fig. 4).

Validation of differentially expressed miRNAs

To validate the sequencing data, five known miRNAs were selected for quantitative analysis. qPCR detection assays were used to confirm the expression of differentially expressed miRNAs in the H and L groups. The 5 miRNAs were verified using real-time PCR. The expression of miR-143, miR-145, miR-199a-3p and miR-221 were up-regulated, whereas miR-486-5p was down-regulated, which was in agreement with the RNA-sequencing results (Fig. 5).