Changes in circulating microRNA levels can be identified as early as day 8 of pregnancy in cattle

Poor reproductive performance remains a major issue in the dairy industry, with low conception rates having a significant impact on milk production through extended calving intervals. A major limiting factor is the lack of reliable methods for early pregnancy diagnosis. Identification of animals within a herd that fail to conceive within 3 weeks after insemination would allow early re-insemination and shorten calving intervals. In a previous study, we found an increase in plasma miR-26a levels in Day 16-pregnant relative to non-pregnant heifers, however changes in miRNA levels that early during pregnancy were very small which likely prevented the identification of robust biomarkers. In this study, we extended our analyses to a wider interval during pregnancy (Days 8 to 60, n = 11 heifers) with the rationale that this may facilitate the identification of additional early pregnancy miRNA biomarkers. Using small RNA sequencing we identified a total of 77 miRNAs that were differentially expressed on Day 60 relative to Day 0 of pregnancy. We selected 14 miRNAs for validation by RT-qPCR and confirmed significant differences in the expression of let-7f, let-7c, miR-30c, miR-101, miR-26a, miR-205 and miR-143 between Days 0 and 60. RT-qPCR profiling throughout Days 0, 8, 16 and 60 of pregnancy showed a distinct increase in circulating levels of miR-26a (3.1-fold, P = 0.046) as early as Day 8 of pregnancy. In summary, in contrast to earlier stages of pregnancy (≤ Day 24), marked differences in the levels of multiple miRNAs can be detected in circulation by Day 60 in cattle. Retrospective analyses showed miR-26a levels to be increased in circulation as early as Day 8, sooner than previously reported in any species, suggesting a biological role for this miRNA in the very early events of pregnancy.


Introduction
Based on government estimates for 2015 [1] the UK dairy herd comprises 1.8 million cattle which produce 14.3 billion litres of milk annually (7,944 litres per animal per annum). This generates approximately £3.7 billion per annum, which is of significant importance for the UK economy overall. However, for the past 50 years the dairy industry has suffered from reduced profitability, despite a continuous net increase in milk production (which has been estimated at 1% per cow annually) [2]. In negative correlation with milk production, fertility decreased during the same period, partially as a consequence of the introduction of North American Cells were obtained from whole blood samples (100 μL) that were centrifuged at 1,900 g for 10 min at 4˚C. The resulting supernatant was discarded and the cell pellet was re-suspended to a volume of 250 μL. The fresh blood cell solution was then used for RNA extraction as described below.
Additional body tissue samples were collected from 3-4 animals at a local abattoir and were placed in ice-cold PBS until arrival at the laboratory. Tissues were then dissected, labelled, weighed and snap-frozen using dry ice. All samples were stored at -80˚C until further use.

RNA extraction
RNA was extracted from 700 μL of plasma using TRIzol LS (Life Technologies, USA) following the manufacturer's protocol. Following homogenisation with TRIzol LS, 3.5 μL of exogenous cel-miR-39-3p (5.6 x 10 8 copies per sample; Qiagen, Germany) were spiked into the homogenate. Prior to RNA precipitation, 40 μg of glycogen (Sigma-Aldrich, USA) were added to each sample to facilitate precipitation of RNA in the presence of 0.1 volume of 5 M ammonium acetate salt (Sigma-Aldrich) as recommended by the manufacturer. RNA pellets were re-suspended in 20 μL of RNase-free water (Qiagen) and used immediately in downstream protocols or frozen at -80˚C until further use.
Blood cell suspensions (250 μL) were homogenised using 3 volumes of TRIzol LS and extracted as per manufacturer instructions. RNA pellets were re-suspended in 10 μL of RNasefree water and frozen at -80˚C until further use.
Other tissues (50 mg) were thawed in 1 mL of TRIzol LS and disrupted with Lysing Matrix D (MP Biomedicals, UK) and a FastPrep FP120 Tissue Disruptor (Thermo Electron, USA) using 2 consequent 30-second runs at a speed setting of 5. Following homogenisation, samples were centrifuged at 10,000 g for 5 min at 4˚C and the resulting supernatant was aspirated. The remaining steps for RNA extraction were carried out as described in the manufacturer's protocol and the resulting RNA pellets were re-suspended in 40 μL of RNase free water. RNA content and quality in blood-cell and tissue samples were monitored using the Nanodrop ND-1000 Spectrophotmeter (average RNA yield = 630 ng and A 260/280 = 1.9, Thermo Fischer Scientific, USA) and RNU6-2 expression levels.

Small RNA sequencing
Small RNA sequencing libraries were prepared from plasma samples using the TruSeq Small RNA Library Preparation Kit (Illumina, USA) and submitted to 36-base single-end sequencing on the HiSeq 2000 Sequencing System (Illumina). Raw sequencing data (publicly available on the GEO database accession GSE83509, [27]) were analysed using the sRNAbench 1.0 tool [28,29]. Briefly, the software was run in genome mode using default settings, with the bovine genome (bosTau4) and miRBase 21 (accessed on 30 April 2015) [30,31] as reference in order to identify bovine miRNAs (bta-miR) and human miRNA homologues (hsa-miR). A single nucleotide mismatch was allowed when mapping raw reads to known miRNA sequences and reads without sequencing adaptor, with undetermined bases and reads below 15 nucleotides (nt) in length being automatically removed from the analysis. For additional details about the integrated analysis steps in sRNAbench please refer to the software manual.
Normalised reads (reads per million mapped, RPMM) obtained from sRNAbench were filtered before being passed to edgeR 3.10.2 for differential expression analysis in R language 3.2.1 [32,33]. Prior to analysis, miRNAs which were detected with less than 25 RPMM in more than 8 samples per experimental group were excluded to reduce technical noise. Differential expression analysis was performed using edgeR's paired mode on 197 miRNAs (S1 File) using GLMfit. The statistical significance threshold was set to false discovery rate (FDR) = 0.05.

RT-qPCR
For plasma samples, cDNA was generated using the miScript II RT Kit (Qiagen) with 2 μL of RNA extract used in each 10 μL reaction. For blood cell and tissue samples, 500 ng of RNA were used in each 10 μL reaction. A Whatman-Biometra Thermocycler (Biometra, USA) was used for the reverse transcription reaction using the conditions recommended by the manufacturer.
The cDNA template was diluted 40-fold and added to 10 μL qPCR reactions that were prepared in 96-well format using the miScript SYBR Green PCR Kit (Qiagen). Template amplification was carried out using commercial species-specific miRNA primers (Qiagen) in an Agilent MX3000P qPCR system (Agilent Technologies, USA), following the cycling conditions suggested by the manufacturer. No-reverse transcriptase and no-template controls were used for each plate. Raw fluorescence data were collected using MxPro software (Agilent Technologies). The amplification efficiency of different primers ranged between 81% and 115%, with R 2 > 0.91. Expression levels were determined using freshly-made standard curves and data were analysed further using Microsoft Excel (Microsoft Corporation, USA). For plasma samples, expression levels were normalised to the mean expression of cel-miR-39-3p and miR-128 (miR-128 was the miRNA with the least variation across samples with coefficient of variation (CV) = 15%). For blood cells and tissues, expression data were normalised to RNU6-2. Samples with distinctly low RNU6-2 levels were excluded.
Statistical analyses were carried out using GraphPad Prism 7 (GraphPad Software, USA). Differences in miRNA expression between Day 0 and Day 60 were assessed using paired ttests. For data involving more than 2 time points one-way rmANOVA was used to determine the effect of Day, followed by Dunnett's test to compare each time point with the Day 0 control. In all cases, normality was tested using the Shapiro-Wilk normality test, outliers were identified with the ROUT test and then removed, and data were log 2 (x+1) transformed to meet the tests' normality criteria. Statistical significance thresholds for all tests were set to P = 0.05.

Genomic cluster and pathway analyses
Information on genomic miRNA clusters in cow was obtained from miRBase 21 (accessed 23/ 06/2016) [31]. Genomic clusters were numbered in order of decreasing size as shown on miR-Base, and the overlap with our sequencing gene list was identified manually.
Target and pathway analysis was carried out using DIANA miRPath 3.0 [34]. Experimentally validated miRNA targets were identified using DIANA-TarBase 7.0 and the Pathway Union option (an a posteriori method) was used for pathway analysis. The integrated Fischer's exact test followed by FDR adjustment were used for statistical analysis. Enriched pathways from the KEGG database were exported from the tool along with the corresponding FDR-corrected P values. A significance threshold of FDR = 0.05 was applied to the corrected P values.

Results and discussion
Small RNA sequencing of bovine plasma We sequenced a total of 22 plasma samples from Days 0 and 60 of pregnancy. On average, we generated 16.1 million raw reads per sample (Table 1). There was a distinct peak in read lengths between 20-23 nt, which matches the length of most mature miRNAs ( Fig 1A). Out of the total raw reads, 5.6 million (34.7%) corresponded to bovine and human mature miRNA sequences ( Table 1, Fig 1B). Across all samples, we identified up to 389 unique miRNAs present at more than 10 reads. The 15 most abundant miRNAs in bovine plasma are listed in Table 2. Eleven of these miRNAs were also among the most abundant in our previous cattle study [23].
From our comprehensive miRNA profiling dataset we sought to identify stable miRNAs which could be utilised for data normalisation in subsequent steps. The miRNA with the lowest variability across samples was miR-128 (CV = 15%) followed by miR-151-3p and miR-425-5p (Table 3). This result was confirmed by NormFinder, which listed miR-128 as the second least variable miRNA in the dataset after miR-181b, with stability values of 0.07 and 0.06, respectively (S1 File, Table 3). We choose miR-128 as an endogenous normaliser for subsequent RT-qPCR assays over miR-181b on the basis of its higher abundance (see 'Methods').  Table 4; S1 File). Fourteen miRNAs changed by more than 2-fold between groups and the largest difference was observed for let-7c, which increased 4.5-fold on Day 60 of pregnancy.

Identification of differentially expressed miRNAs on Day 60 of pregnancy
Given the identification of miRNA clusters associated with pregnancy in other species [35], we next sought to determine whether the miRNAs differentially expressed between Days 0 and 60 of pregnancy belonged to specific genomic clusters. A total of 82 miRNA clusters, each containing 2 to 47 miRNAs, are registered for cow in miRBase 21. Using this as reference, we found that 31 out of the 77 differentially expressed miRNAs belonged to a cluster, being present together with at least another miRNA from the same cluster (Table 5). In total, 13 miRNA clusters were represented among differentially expressed miRNAs. Among these were the miR-379~656 cluster (known as C14MC in humans), represented by 5 miRNAs, and the miR-17~92 and miR-106b~25 clusters, represented by 5 and 2 miRNAs, respectively ( Table 5). The C14MC cluster is known to be expressed in the placenta and is involved in foetal development; circulating levels of C14MC miRNAs increase during human pregnancy [36,19,37,38]. The miR-17~92 and miR-106b~25 clusters are paralogues, with a demonstrated role in vascular endothelial growth factor (VEGF)-mediated angiogenesis through targeting of phosphatase and tensin homolog (PTEN), and a possible role in placental immune tolerance during the first trimester of pregnancy; both clusters are highly expressed in the human placenta [36,[39][40][41]. Finally, several members of the let-7 family known to be expressed in the human placenta  [36] were present across multiple genomic clusters (Table 5). Overall, these results suggest functionally conserved roles of a subset of pregnancy-associated miRNA clusters across species.
To gain insight on the roles of the differentially expressed miRNAs during pregnancy we carried out pathway enrichment analysis, using the DIANA web-tool [34], separately using experimentally validated targets of all the up-regulated and down-regulated miRNAs on Day 60 of pregnancy (39,555 and 17,439 targets, respectively). There was significant overlap between enriched pathways targeted by up-regulated and down-regulated miRNAs (S2 File), and these common pathways included fatty acid synthesis and metabolism, extracellular cell membrane (ECM)-receptor interaction, adherens junction, Hippo signalling and the thyroid hormone, p53 and cell cycle pathways. Up-regulated miRNAs in pregnancy also targeted the forkhead O (FOXO) and transforming growth factor beta (TGF-β) signalling pathways. These results suggest roles of miRNAs in regulating cell metabolism as well as proliferation/survival/ migration processes involved in tissue growth during pregnancy.

Validation of sequencing data using RT-qPCR
Previous studies with human tissues have shown only partial agreement in miRNA expression data generated using different profiling platforms [42]. Thus, we set out to validate our results from sequencing using RT-qPCR. Fourteen differentially expressed miRNAs were selected for validation (Table 6) based on both their relative abundance (RPMM) and fold-change between Days 0 and 60 (Table 4, S1 File). For 11 out of the 14 miRNAs, RT-qPCR results agreed with those from sequencing ( Table 6, Fig 2A), although statistical significance was obtained only for 7 miRNAs, namely let-7c, -7f, miR-101, -30c, -26a, -143 and -205. Six miRNAs were expressed at significantly higher levels (up to 6-fold, P < 0.05) on Day 60 compared to Day 0, while miR-205 was significantly lower on Day 60 (1.5-fold, P < 0.05, Fig 2B, Table 6). Our subsequent analyses focused on these 7 miRNAs.
Potential roles for these miRNAs during pregnancy can be speculated based on previous literature and the results of gene ontology analyses above (S2 File). Four of the seven miRNAs, miR-143, miR-101, miR-30c and miR-26a, have been shown in previous studies to regulate insulin signalling and lipid metabolism and to be involved in diabetes [43][44][45][46], and thus they may play important roles regulating changes in maternal metabolism occurring early during pregnancy. In addition, miR-26a and let-7f have been shown to regulate macrophage-and Tcell-mediated immunity [47][48][49][50] which would be consistent with a role in establishing maternal immuno-tolerance, a key function during early pregnancy. Finally, miR-101, miR-205 and miR-26a are known regulators of angiogenesis [51][52][53][54], a key process during development of foeto-placental and maternal tissues during pregnancy. Tissue-wide expression of plasma miRNAs associated with pregnancy Next, to investigate the origin of these miRNAs, and at the same time gain further insight on their potential roles during early pregnancy, we profiled the expression of all 7 miRNAs across 14 different bovine tissues, including blood cells. A total of 4 outliers were identified and removed from each of let-7f, let-7c, miR-101 and miR-143 datasets. Six miRNAs (let-7c, let-7f, miR-101, miR-143, miR-30c and miR-26a) were expressed in many different body tissues including the uterus and the placenta, with no particular tissue being clearly enriched for any particular miRNA (Fig  3), overall consistent with data from human tissues [55,56]. Notably, relative expression levels for most miRNAs were lowest in both early (Day 70) and term placenta, indicating the developing conceptus may not be a significant source of these miRNAs in plasma. In addition, expression profiles suggested that, for all miRNAs except miR-143 and miR-205, blood cells are a significant source of plasma levels, consistent with roles in immuno-regulation as suggested above. Interestingly, mean levels of miR-205 were more than 8,000-fold higher in skin than in any other tissue (Fig 3). In mammals, miR-205 was reportedly expressed in a range of tissues including the thymus, mammary gland, early embryo, digestive tissue and skin [57][58][59]55], however, different studies highlighted the importance of miR-205 in regulating keratinocyte function [60,61], which is consistent with our results. Whether skin is indeed a significant contributor to circulating levels of miR-205 should be explored in future studies. A caveat associated with our tissue profiling approach is that all tissues analysed (other than placenta) were collected from non-pregnant animals, which would not account for changes in tissue miRNA expression that may only occur during pregnancy. However, since our goal was to determine whether any of the miRNAs analysed (Fig 3) was specifically expressed or highly enriched in a particular tissue, which we can reasonably assume will not depend on pregnancy status, then the use of tissues from non-pregnant heifers for profiling those miRNAs was justified. Profiling of circulating miRNAs throughout early pregnancy Finally, we sought to obtain further insight on these pregnancy-associated miRNAs by determining their temporal changes in circulation throughout early pregnancy. We hypothesised that some of the miRNAs identified on Day 60 may actually begin increasing at earlier stages of pregnancy (at or before Day 16), which could provide useful biomarkers of early pregnancy, possibly in the form of a miRNA panel. Accordingly, we profiled all 7 miRNAs by RT-qPCR on sequential plasma samples collected on Days 0, 8, 16 and 60 of pregnancy. We identified and removed one outlier value for each of miR-101, miR-26a and let-7f on Day 8, and miR-19a on Day 0. Analyses of the remaining data showed a significant effect of Day of pregnancy (P < 0.05) for all 7 miRNAs (Fig 4). Moreover, mean levels of many miRNAs did not change significantly until Day 60 (P < 0.05), with some miRNAs increasing gradually (non-significantly) throughout early pregnancy, e.g. miR-30c. Remarkably, however, mean levels of miR-26a were distinctly increased as early as Day 8 (3.1-fold; P = 0.046), resulting from a net increase in miR-26a levels (ranging from 1.3-to 11.7-fold) in 8 of 10 animals between Days 0 and 8. In addition, levels of miR-101 tended to increase (2.3-fold, P = 0.069) on Day 16 relative to Day 8 with no further elevation in mean levels of that miRNA up to Day 60 (Fig 4). Changes in miRNA levels in circulation this early during pregnancy (i.e., Day 8) have not been reported before and suggest a novel role for miR-26a, possibly through its immunomodulatory effects [48,49], very early during the establishment of pregnancy, a proposition that clearly deserves further investigation. Although miR-26 has also been shown to play different roles presumably unrelated to pregnancy [62][63][64], the robust (3-fold) increase in plasma miR-26 levels between Days 8 and 0 of pregnancy suggest that miR-26a could be potentially used as a diagnostic biomarker as early as the first week of pregnancy, which could significantly facilitate reproductive management of modern dairy herds. Lastly, because changes in miR-26a and miR-205 showed opposite trends during early pregnancy we calculated the ratio between these 2 miRNAs to determine whether this could provide an improved diagnostic biomarker; we found the miR-26a / 205 ratio to increase 7.5-fold (P < 0.001) between Days 0 and 8 of pregnancy, providing an even more robust, potential molecular indicator of very early pregnancy in cattle (Fig 4).

Conclusions
Our novel results using small RNA sequencing and RT-qPCR profiling show wide changes in circulating miRNA profiles occurring by Day 60 of pregnancy in cattle, in contrast with the  [23]. Specifically, we identify a subset of miRNAs (let-7f, let-7c, miR-30c, miR-101, miR-26a, miR-205 and miR-143), the levels of which increase distinctly in circulation (up to 6-fold) in Day 60 pregnant relative to non-pregnant (Day 0) cows, and which provide novel molecular candidates involved in the establishment of pregnancy in cattle. Significantly, an increase in plasma levels of miR-26a was identified as early as Day 8, much earlier than previously reported for any miRNA during pregnancy in any species, suggesting one of the biological role of this miRNA may be in the very initial stages of pregnancy, providing a potential unique diagnostic biomarker of early pregnancy. Overall, the subset of miRNAs identified in our study provide the basis to elucidate specific roles of miRNAs during early pregnancy and for more comprehensive biomarker identification studies in the future.