The Dynamics of microRNA Transcriptome in Bovine Corpus Luteum during Its Formation, Function, and Regression

The formation, function, and subsequent regression of the ovarian corpus luteum (CL) are dynamic processes that enable ovary cyclical activity. Studies in whole ovary tissue have found microRNAs (miRNAs) to by critical for ovary function. However, relatively little is known about the role of miRNAs in the bovine CL. Utilizing small RNA next-generation sequencing we profiled miRNA transcriptome in bovine CL during the entire physiological estrous cycle, by sampling the CL on days: d 1–2, d 3–4, and d 5–7 (early CL, eCL), d 8–12 (mid CL, mCL), d 13–16 (late CL, lCL), and d > 18 (regressed CL, rCL). We characterized patterns of miRNAs abundance and identified 42 miRNAs that were consistent significantly different expressed (DE) in the eCL relative to their expression at each of the analyzed stages (mCL, lCL, and rCL). Out of these, bta-miR-210-3p, −2898, −96, −7-5p, −183-5p, −182, and −202 showed drastic up-regulation with a fold-change of ≥2.0 and adjusted P < 0.01 in the eCL, while bta-miR-146a was downregulated at lCL and rCL vs. the eCL. Another 24, 11, and 21 miRNAs were significantly DE only between individual comparisons, eCL vs. the mCL, lCL, and rCL, respectively. Irrespective of cycle stage two miRNAs, bta-miR-21-5p and bta-miR-143 were identified as the most abundant miRNAs species and show opposing expression abundance. Whilst bta-miR-21-5p peaked in number of reads in the eCL and was significantly downregulated in the mCL and lCL, bta-miR-143 reached its peak in the rCL and is significantly downregulated in the eCL. MiRNAs with significant DE in at least one cycle stage (CL class) were further grouped into eight distinct clusters by the self-organizing tree algorithm (SOTA). Half of the clusters contain miRNAs with low-expression, whilst the other half contain miRNAs with high-expression levels during eCL. Prediction analysis for significantly DE miRNAs resulted in target genes involved with CL formation, functionalization and CL regression. This study is the most comprehensive profiling of miRNA transcriptome in bovine CL covering the entire estrous cycle and provides a compact database for further functional validation and biomarker identification relevant for CL viability and fertility.

The formation, function, and subsequent regression of the ovarian corpus luteum (CL) are dynamic processes that enable ovary cyclical activity. Studies in whole ovary tissue have found microRNAs (miRNAs) to by critical for ovary function. However, relatively little is known about the role of miRNAs in the bovine CL. Utilizing small RNA next-generation sequencing we profiled miRNA transcriptome in bovine CL during the entire physiological estrous cycle, by sampling the CL on days: d 1-2, d 3-4, and d 5-7 (early CL, eCL), d 8-12 (mid CL, mCL), d 13-16 (late CL, lCL), and d > 18 (regressed CL, rCL). We characterized patterns of miRNAs abundance and identified 42 miRNAs that were consistent significantly different expressed (DE) in the eCL relative to their expression at each of the analyzed stages (mCL, lCL, and rCL). Out of these, bta-miR-210-3p, −2898, −96, −7-5p, −183-5p, −182, and −202 showed drastic up-regulation with a fold-change of ≥2.0 and adjusted P < 0.01 in the eCL, while bta-miR-146a was downregulated at lCL and rCL vs. the eCL. Another 24, 11, and 21 miRNAs were significantly DE only between individual comparisons, eCL vs. the mCL, lCL, and rCL, respectively. Irrespective of cycle stage two miRNAs, bta-miR-21-5p and bta-miR-143 were identified as the most abundant miRNAs species and show opposing expression abundance. Whilst bta-miR-21-5p peaked in number of reads in the eCL and was significantly downregulated in the mCL and lCL, bta-miR-143 reached its peak in the rCL and is significantly downregulated in the eCL. MiRNAs with significant DE in at least one cycle stage (CL class) were further grouped into eight distinct clusters by the self-organizing tree algorithm (SOTA). Half of the clusters contain miRNAs with low-expression, whilst the other half contain miRNAs with high-expression levels during eCL. Prediction analysis for significantly DE miRNAs resulted in target genes involved with CL

INTRODUCTION
In mammalian reproduction the corpus lutem (CL) is a pivotal gland, primarily because of its role in regulating the length of estrous cycle. The CL development is completed rapidly after the egg ovulation from the follicular remnants and its structural features remain stable, it instead secretes the hormone progesterone (P4), at levels required for ovarian cyclical activity and maintenance of pregnancy when it occurs (Reynolds and Redmer, 1999;Davis and Rueda, 2002;Stocco et al., 2007;Rekawiecki et al., 2008). As a highly dynamic tissue, exhibiting periodical growth, functionalization and regression, the CL is an ideal example for the study of diverse cellular processes that occur in it: cell proliferation and apoptosis, immune system activation, down regulation of functional proteins, and cellular matrix remodeling events (Davis et al., 2003;Schams and Berisha, 2004;Sugino and Okuda, 2007;Yoshioka et al., 2013).
Over the past few years many studies have helped elucidate the regulatory events underlying bovine CL development and their function/s at the transcription level. These studies have highlighted and further reinforced the importance of several classes of growth factors as key regulators of CL formation (angiogenesis), function (Berisha et al., 2000b(Berisha et al., , 2002, and CL regression (Reynolds et al., 2000;Neuvians et al., 2004;Berisha et al., 2008Berisha et al., , 2010Berisha et al., , 2016Berisha et al., , 2017Miyamoto et al., 2010;Nishimura and Okuda, 2010). Although the well-orchestrated transcriptional changes in bovine ovary and CL are described in considerable details (Berisha et al., 2000a(Berisha et al., ,b, 2002(Berisha et al., , 2008(Berisha et al., , 2010(Berisha et al., , 2016(Berisha et al., , 2017Neuvians et al., 2003Neuvians et al., , 2004Kliem et al., 2007;Shimizu et al., 2007;Sarkar et al., 2011), very little in regard to the post-transcriptional regulation has been addressed. In particular the role of post-transcriptional regulators such as miRNAs that act locally and in a cell-type specific manner is not clear yet (Wienholds and Plasterk, 2005). MiRNAs are short non-coding and single-stranded nucleic acids that regulate the translation of mRNA transcripts during tissue development and differentiation (Lee et al., 1993;Bartel, 2004;Kloosterman and Plasterk, 2006;Stefani and Slack, 2008). They regulate the abundance of gene protein products by inhibiting translation or by degrading the targeted mRNA itself (Shivdasani, 2006). Gene expression in the CL as an endocrine gland is subjected to dynamic regulations that are likely influenced by miRNAs from development through luteolysis. Indeed, there is experimental evidence supporting the involvement of miRNAs in the formation of CL. Using a Dicer1 (an RNase III enzyme required for processing miRNAs) hypomorphic mouse line (Dicer d/d ), Otsuka et al. (2008) demonstrated that the females form this line were infertile. The observed infertility was shown to be due to the altered CL development resulting mainly from an impaired vascular network. The mechanism for this phenotype was linked to the up-regulation of the anti-angiogenic transcription factors: tissue inhibitor of metalloproteinase 1 (TIMP1) and platelet factor 4 (PF4). Both transcription factors were shown to be predicted targets of the miRNAs: 17-5p and let-7b, and indeed their expression were restored to baseline by the injection of these two miRNAs. In bovine the preovulatory luteinizing hormone (LH) surge is the key event leading to major cellular phenotype changes in follicles and changes in CL transcriptome (Sood et al., 2006;Christenson et al., 2013). The ongoing cellular alterations suggest there might be stage specific populations of miRNAs in the developing, fully functional and then regressing CL, which could be induced by the secretion of other hormones e.g., P4 .
The published data involving miRNAs in the regulation of CL function in large farm animals are rather limited and conducted on either whole ovary tissue (Tripurani et al., 2010;Sontakke et al., 2014) or in the CL from a specific stage (Maalouf et al., 2014(Maalouf et al., , 2016, consequently a complete and detailed picture is missing. Therefore the objective of this study was to analyze the miRNA transcriptome dynamics at different timely defined CL classes (days). Specifically, we sequenced miRNAs transcriptome during early CL formation (C1: d 1-2, C3: d 3-4, and C5: d 5-7), mid CL function (C8: d 8-12), late CL function (C13: d 13-16) and during the CL regression (C18: d > 18). We identified 42 miRNAs that were consistently significant DE across the analyzed CL classes. Further 24, 11, and 21 miRs showed significant DE only between individual comparisons eCL vs. mCL, lCL, and rCL, respectively.
Results from the present study will create knowledge on miRNA transcriptome in the field of animal reproduction and provide a fertile ground for any further functional validation studies or biomarker identification beyond reproduction.

Animals and Tissue Collection
The corpora lutea were collected within 20 min post the routine slaughter from German Fleckvieh Cattle at local abattoirs. The isolated luteal tissue was immediately frozen in liquid nitrogen and stored at −80 • C until used in experiments.
The stage of the estrous cycle was determined by macroscopic examination as previously described (Berisha et al., 2016(Berisha et al., , 2017. In brief: macroscopic evaluation of the whole ovaries by assessing their color, consistency and of the uterus by assessing its color, consistency and the mucus were considered for the determination of CL stages. According to that the CL were assigned into the following classes (days): early CL (eCL, C1: d 1-2, C3: d 3-4, C5: d 5-7), mid CL (mCL, C8: d 8-12), late CL (lCL, C13: d 13-16) and regressed CL (rCL, C18: d >18). A total of 18 corpora lutea, each from individual animals, were included in the study. Three corpora lutea were used for each analyzed CL classes (C1, C3, C5, C8, C13, and C18) and subsequent small RNA-Seq.

RNA Extraction and MicroRNA Quantification
For the extraction of total RNA and subsequent purification of miRNA, slices from frozen CL were cut and processed using the miRNeasy Mini Kit (Qiagen, Germany) according to the manufacture's protocol. RNA was quantified using a Nanodrop spectrophotometer (Thermo Fisher Scientific, Germany), purity was determined based on the optical density ratio and mRNA quality was assessed on a Bioanalyzer 2100 (Agilent Technologies, Germany) using the RNA 6000 Nano Assay. Only CLs with a RIN value higher than 8 were considered (Mueller and Schroeder, 2004;Fleige and Pfaffl, 2006). In order to determine the percentage of small RNA relative to total RNA, samples were further analyzed using the Small RNA Kit (Agilent Technologies, Germany), which allows the quantification of small RNA ranging in size between 6 and 150 nts.

Preparation of Library for Small RNA Sequencing
Library preparation was performed using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England BioLabs Inc., USA). A total of 200 ng of purified RNA was used as starting material. This was converted into barcoded cDNA according to Protocol E7330 (New England BioLabs Inc., USA) and sequenced on the Illumina HiSeq2500 platform (Illumina Inc., USA). Briefly the steps are: 3 ′ -adaper ligation, followed by reverse transcription primer hybridization (to prevent adaptor-dimer formation), 5 ′ -adapter ligation, a reverse transcription reaction and PCR enrichment (Spornraft et al., 2014). The PCR products were then purified with the MinElute PCR Purification Kit (Qiagen, Germany) and their length checked on the Bioanalyzer 2100 (Agilent, Germany) using the DNA 1000 Kit (Agilent, Germany). The cDNA constructs were size selected by loading the PCR products on a 4% agarose gel. Bands corresponding to products of 140-160 bp in length were cut (Supplementary Figure 1) and processed by the MinElute Gel Extraction Kit (Qiagen, Germany). The size of the cDNA library was validated using the High Sensitivity DNA Kit (Agilent Technologies, Germany) (Supplementary Figure 1). Library constructs were then sequenced on the HiSeq2500 Illumina platform.

Sequencing Quality and Mapping
The sequencing data were processed using in-house tools (Spornraft et al., 2014(Spornraft et al., , 2015Buschmann et al., 2016;Schanzenbach et al., 2017). The trimming program, Btrim (Kong, 2011), was used to trim the adaptor sequence from the 3 ′ ends and all reads without detectable adaptor sequences were removed from the dataset. Next the length distribution and the base call accuracy, represented by the phred quality scores (Q score), were calculated using the quality control software FastQC (Babraham Bioinformatics, UK, Version 0.10.1). In order to prevent bias and false-positive mapping by degenerated RNA or other small RNA material, all reads shorter than16 nts were removed from the dataset. The remaining reads were filtered for bovine rRNA, tRNA, snRNA, and snoRNA sequences by mapping them to the RNAcentral database (RNAcentral Consortium, 2015). The filtered dataset was aligned to the most recent miRBase database for mature bovine miRNAs (release 21) (Kozomara and Griffiths-Jones, 2014) by using the BowTie short read aligner (Langmead et al., 2009). Default parameters and settings were used except for the "best" alignment algorithm, which was adapted by allowing only a single mismatch across the whole sequence. All aligned reads were subsequently sorted and indexed by SAMtools (Li et al., 2009), the read counts were then generated by retrieving the sum of hits per miRNA sequence.
RNA pools lacking an RT enzyme from different CL samples as well as an additional water control sample were included in each qPCR run. The expression data were normalized using the Hs_RNU6-2_1 gene and measurements were undertaken following the Pfaffl method (Pfaffl, 2001).

Statistical Analysis
Statistical analysis on read counts, representing the frequency of the sequenced miRNAs, was performed using the DESeq2 R script (v.3.1.2). This methodology is based on a negative binomial distribution, linking the mean and variance by local regression (Love et al., 2014) and determines the adjusted P-value (FDR; Benjamini and Hochberg adjustment) allowing for across sample comparison. A read count noise-cut-off >50 was applied for all miRNAs identified as DE. Distribution of significantly DE miRNAs across the estrous cycle was assessed by the web application Bio Venn (Hulsen et al., 2008). SOTA: Self-organizing Tree Algorithm and Unsupervised Hierarchical Clustering (HCL) The SOTA and HCL were used to search for correlating miRNA expression patterns (Dopazo and Carazo, 1997;Herrero et al., 2001). All miRNAs with average read counts >50 and at least one cycle stage of significant DE were included in the search. The mean centered expression was calculated as follows: mean read values of replicate CL samples (n = 3) for each cycle stage was assessed and log2 transformed. In the next step, the mean value of all samples was subtracted from the cycle stage specific value. HCL and SOTA were performed using the MeV_4 Software (v. 4.9.0) (Saeed et al., 2003) with default settings. For SOTA a limit of seven cycles was applied allowing for a total of eight generated clusters.

Micro RNA Target Identification and Pathway Analysis
Identification of putative miRNAs targets for DE miRNAs common to all conditions (# 42miRNAs) was performed using the miRTarBase release 6.0 (Chou et al., 2016) and TargetScan release 7.0 (http://www.targetscan.org). We selected miRNAtarget interactions (MTIs) that were experimentally validated in human and rodents by RT-PCR, reporter assay or Western blot as detailed in http://mirtarbase.mbc.nctu.edu.tw/ and did an in silico prediction of their targets within the bovine genome using the TargetScan (v7.0).
The open access DIANA miRPath v.30 software (Vlachos et al., 2015) was used to generate pathways enriched for genes targeted by miRNAs identified as DE between individual comparisons: eCL vs. mCL, lCL, and rCL. The latest version of DIANA microT-CDS (v5.0) with a microT threshold of 0.8 and a false discovery rate (FDR) correction for the threshold p < 0.05 were applied. To exclude any pathways containing only a few targeted gene nodes the modified Fisher's exact test was applied. DIANA-miRPath database and the miRTarBase consist of validation studies that have been mainly performed on human or rodents and less in bovine, however miRNAs are highly functionally conserved between species and their nomenclature is based on seed sequence similarities (Friedman et al., 2009).

Small RNA Sequence Quality Assessment
Each sample, representing one animal, was analyzed for the presence of high quality RNA reads. The assessment of the quality of sequenced samples was performed using the fastqc algorithm. Both, per base and per sequence, quality scores were considered. In addition, the mean sequence quality (Phred Score) for all 18 sequenced CL samples was evaluated. Overall, 98% of all reads showed a quality (Q) score of 37 (Supplementary Figure 2). It is known that base calls with a Q-score higher than 30, contain one incorrect base call in 1,000 and therefore have an error probability of 0.001 (Ewing and Green, 1998). According to the HiSeqIlumina platform, if more than 90% of reads have an average Q-score of 30 or greater, they are considered reads of high quality (Illumina, 2011). The total sequence length showed a clear pattern in its distribution with a peak at 22 nts (Supplementary Figure 2), which corresponds to the known length of miRNA and the mean GC content per sequence was around 47% (Supplementary Figure 2), ensuring thus non-biased annotation of the sequencing data, due to the ubiquitously genome presence of GC regions.

Evaluation of Small RNA Sequencing Data
The mean number of sequence reads across the 18 small RNA Seq analyses were 15,244,516 ± 4,070,073 (SD) (Supplementary Table 1). In total 9,211,818 ± 2,048,812 (SD) sequences passed the trimming threshold, whilst 6,031,779 ± 2,039,033 (SD) did not and were excluded from any further analysis. The sequences that passed trimming were further mapped to the RNAcentral database to remove all small RNAs, which were not microRNA in nature (RNAcentral Consortium, 2015). The resulting number of filtered sequences for each sample was 3,683,423 ± 653,417 (SD), while a mean of 5,730,324 ± 1,438,850 (SD) sequences did not map to RNAcentral. These sequences were annotated by mapping them to the bovine mirbase (version 21, Kozomara and Griffiths-Jones, 2014). This allowed a mean 898,679 ± 289,266 (SD) sequences across the CL samples to be mapped onto annotated miRNAs, making up 5.1 ± 2.9% (SD) of the total read counts (Supplementary Table 1).
The overall numbers of read were well-correlated between eCL and the three other stages ( Figure 1A). Yet, by the acquisition of overall data behavior as can be seen in Figure 1B, we found that the eCL was clearly separated from mCL, lCL, and rCL. Performing principal components of scores and loadings (SL-plot), two distinct and significantly different clusters were identified, one corresponding to the eCL class ( Figure 1B, loadings, top right) and the other to the mCL, lCL, and rCL classes ( Figure 1B, loadings bottom right). The two clusters were separated by principal component 1 (PC1) which explained 99.8% of the total data variance and clearly differentiated eCL from mCL, lCL, and rCL. PC1 is determined by the highly abundant and significantly DE miRs: miR-21-5p ( Figure 1B, scores upper quadrant and Figure 1C) and miR-143 ( Figure 1B, scores lower quadrant and Figure 1F). The PC2 is mainly determined by less abundant, hence significant DE miRs: miR-146a (

Identification of miRNAs in the CL at Different Cycle Stages
Over all CL samples we found 439 miRNAs with more than one read in average, of them 205 account for 99.7% of the total reads. Furthermore, 8.65% of the annotated miRNAs had 1,000-250,000 reads, 17.2% had between 50 and 1,000 reads and 74.15% show fewer than 50 reads (Supplementary Table 2).
The expression level of the most abundant miRNAs was determined by their number of read counts, which reflects how frequently they were sequenced. The top 20 most abundantly expressed miRNAs for each cycle stage are shown in Figures 1C-F. Interestingly, we found that two miRNAs, bta-miR-21-5p and bta-miR-143, show constantly greater expression and account for more than 40% of the total number of reads in the CL samples across the cycle. A decrease in number of read counts, and hence abundance, was seen for bta-miR-21-5p from early to mid, late and the subsequent regression whilst the opposite was observed in the case of bta-miR-143 (Figures 1C-F). Several other miRNAs were observed that vary in expression across the CL cycle. For example, bta-miR-202 was highly expressed during the early CL, however its number of reads decreases in the subsequent CL classes and in the rCL is not among the top 20 most abundantly expressed miRNAs ( Figure 1D). Conversely, bta-miR-30a-5p is highly abundant in mCL, lCL, and rCL, but not in the eCL (Figures 1B-D). Notable, the list of top 20 most abundantly expressed miRNAs shows similar expression pattern Frontiers in Genetics | www.frontiersin.org FIGURE 1 | Acquisition of read correlation, general data behavior and miRNAs abundance. (A) Correlation of NGS normalized reads between the comparisons eCL vs. mCL, lCL, and rCL. (B) Scores and loadings plot (SL-plot) showing the specific clustering of miRNAs from six different CL classes. The normalized reads are clustered using DESeq2. Each CL class consists of three CL from individual animals (e.g., C1, C1.1, andC1.2 belong to the CL class C1). (C-F) Top 20 most abundant miRNAs at each CL class. Black bars at either mCL, lCL, or rCL represent miRNAs abundant at these CL classes and absent in the eCL. The orange bars at eCL represent miRNAs abundant at this CL class and absent in the subsequent CL classes.
in the eCL, mCL, and lCL, while greater discrepancy is observed between eCL and rCL (Figures 1A,D, bars with either black or orange color).

Significantly DE miRNAs in the CL during Estrous Cycle
Using the DESeq2 algorithm (Anders and Huber, 2010) we identified differentially expressed miRNAs at the comparisons between stages, namely eCL to mCL, lCL, and to rCL. In the comparison between eCL and mCL samples we identified 92 DE miRNAs (Supplementary Table 3A). Of these, 42 miRNAs were also found to be DE in two other CL classes, lCL and rCL, vs. the eCL (Figures 2A-C, Table 1).
In contrast, 24 miRNAs were DE regulated only between individual comparison of eCL and mCL, while a total of 19 miRNAs were common to comparison of eCL to mCL and eCL to lCL. Other 7miRNAs were common to the comparisons of eCL to mCL and eCL to rCL (Supplementary Tables 3A, 4A). The complete list of DE miRNAs in eCL vs. the other cycle stages as well as all common or unique miRNAs to a given comparison can be accessed as supplementary information (Supplementary  Tables 3A-C, 4A,B).
Overall, 21 miRNAs out of the 42 miRNAs significantly DE across all comparisons were seen to be upregulated in the mid, late and regressed CL stages relative to their expression in the eCL (Figures 2B,C ). Yet, seven of them: bta-miR-210, −2898, −96, −7-5p, −183-5p, −182, and −202 show a fold change (FC) upregulation ≥2.0 with an adjusted P < 0.01 in at least one of the three CL classes (mCL, lCL, and rCL) vs. the eCL. In contrast, only one miRNA bta-miR-146a was downregulated with a FC A total of 42 miRNAs were common to all comparisons. The read counts for each miRNA are given as the base mean of three CLs replicates. The mean read count ratio between either early CL (eCL) and mid CL (mCL), or late CL (lCL) and regressed CL (rCL), is shown as the fold change (FC). MiRNAs are listed in the order of their adjusted P-value (Padj). Down-regulated miRs: FC < 1. Up-regulated miRs: FC >1. The fold changes and adjusted P-values were calculated using DESeq2. ‡ FDR, Benjamini and Hochberg correction.
≥ 2.0 and adjusted P < 0.01 in the lCL and rCL vs. the eCL. MiRNAs with FC ≥ 2.0 and adjusted P < 0.01 are indicated with asterisk ( Figures 2B,C  The RT-qPCR fold change regulation was calculated using the 2 −ddCt . Significance between the CL classes was assessed by Student's t-test on the normalized Ct values (dCt-Value). Significant regulation (p< 0.05) are indicated in bold. Detailed information for individual miRNA level of significance regulation is provided in the Supplementary Figure 3. NGS normalized reads are significantly regulated in all three comparisons and are detailed on Table 1.
As the cycle progresses the number of miRNAs found to be significantly DE was slightly decreased and 84 miRNAs were regulated between eCL early and lCL (Supplementary Table 3B). Of these 11 are unique as beeing regulated only in eCL vs. lCL and not in other comparisons (Figure 2A, Supplementary Table  4B). We next compared eCL to rCL and identified a total of 82 DE miRNAs (Supplementary Table 3C
For six miRNAs (bta-miR-210, −96, −182, −183, −146a, and −202) out of the eight selected miRNAs the qPCR results were consistent with those from the NGS and all of them are upregulated in the eCL vs. the three other stages. Statistical significance was obtained and indicated in Tables 2, 3, Supplementary Figure 3.
Two miRNAs (bta-miR-2898 and −7) however, show only a trend in the regulation comparing qPCR and NGS results. This is not unusually, as different studies with human tissues have shown only partial agreement between expression data generated by different platforms (Lee et al., 2014). The trends in regulation of miR-2898 and miR-7 and corresponding significances are provided in Supplementary Figure 3.
Bta-miR-21-5p and -143, two highly abundant miRNAs counting for the separation of eCL from the other CL classes, did also exhibit the same trend of regulation comparing the qPCR and NGS results. The regulation was statistically significant between eCL and mCL for bta-miR-21-5p and between eCL and all comparisons in the case of bta-miR-143 (Tables 2, 3).
Although significant regulation could be confirmed for eight out of the ten selected miRNAs, the qPCR data show lower fold-change regulation than the fold-change regulation assessed by NGS. This is in agreement with results by different studies profiling miRNAs in tissue and plasma Donadeu, 2016, 2017).

SOTA and Unsupervised Hierarchical Clustering (HCL) of miRNAs
Eight different clusters consisting of miRNAs with similar expression profiles were generated using SOTA ( Figure 3A) and the greatest distance is calculated between the cycle stage C1 and C18. Each cluster consists of 10 to 36 miRNAs, with an average cluster diversity of 0.87 ± 0.02 ( Figure 3B). According to HCL (Figure 3C, cycle time-point individually considered, see Materials and Methods section) the cycle stage C1 and C3 consist of miRNAs with similar expression pattern and the same is observed for cycle stage C8 and C13.
SOTA clusters are clearly split into two separate groups, each with an equal number of clusters. Furthermore, half of the clusters contain miRNAs with low-expression, whilst the other half contains miRNA with high-expression levels at the beginning of CL development. Again, whilst no significant changes were found during the very early CL (cycle stage: C1 and C3), for each cluster, except for cluster 1, a change in expression can be seen between samples from the cycle stage C5 (days 5-7) and those from C8 (days 8-12). This represents the time at which the CL reaches its functional state (Reynolds and Redmer, 1999). The clusters labeled 1, 2, 3, and 8 are especially interesting as the miRNAs in these clusters are primarily expressed differentially ( Figure 3A) during the functional stage (cluster 1: up regulation between days 5 and 12, cluster 2 and 3: up regulation between days 8 and 12, cluster 8: down regulation between days 8 and 18 (Figures 3A,B). In contrast, miRNAs in other clusters are sequentially up or down regulated during the CL life span (cluster 4: up regulation; cluster 5, 6, 7: down regulation (Figures 4A,B). Interestingly the two clusters were observed to generally exhibit opposite expression patterns, for instance cluster 2 vs. cluster 8 and cluster 4 vs. cluster 5. The list with all individual clusters and their corresponding miRNAs can be accessed in the supplementary data (Supplementary Table 5).

"In-Silico" Functional Annotation and Pathways Enrichment Analysis
The search on the miRTarBase for functional annotation of the DE miRNAs resulted in the identification of numerous angiogenesis promoting genes, such as vascular endothelial growth factor A (VEGFA), fibroblast growth factors 2 (FGF2), fibroblast growth factors 9 (FGF9), hypoxia inducible factor A (HIFA), angiopoietin 1 (ANGPT1), and angiopoietin 2 (ANGPT2). reads in the eCL. The fold changes and adjusted P-values were calculated using DESeq2 (**adjusted P < 0.01). Error bars are shown as log fold-change standard error (lfcSE). Note: only miRNAs with FC > 2.0 are marked with asterisk; adjusted P-values for miRNAs with FC < 2 are presented in Table 1.
Using the DIANA tool miRPath v3.0 (Vlachos et al., 2015) with a threshold of P < 0.05 we searched for pathway enrichment for genes targeted by miRNA species that were unique to and DE only between individual comparisons of eCL and other three CL classes. The main molecular pathways all being those related to eCL are listed on Figure 4A (eCL vs. mCL, lCL, and rCL). The top pathways which were targeted by 3 or more unique miRNAs in the assessments between: eCL and mCL were the fatty acid biosynthesis, transforming growth factor-B (TGFB) and Wntsignaling pathway; between: eCL and lCL, the glutamatergic synapse, axon guidance, thyroid hormone signaling, oxytocin and prolactin synthesis and between: eCL and rCL, the ECMreceptor interaction, the phosphatidylinositol 3 ′ -kinase (PI3K)-Akt signaling and mitogen-activated protein kinase (MAPK)pathway. MiRNAs that contributed to gene enrichment for the top 5 KEGG are given separately for each of the comparison in Figure 4B, whereas all significantly enriched pathways and their respective number of genes and number of miRNAs can be found in Supplementary Table 6.

DISCUSSION
The series of sequencing data from this study provide the most comprehensive profile of miRNAs transcriptome in bovine CL during the entire physiological estrous cycle. We used small RNA Seq (RNA Seq) to identify miRNAs expressed in the CL and found 439 miRNAs with at least one read, a number that is comparable to what other recent studies in bovine reproductive tissue have reported (Hossain et al., 2009;Maalouf et al., 2014). The identified miRNAs in this study were further annotated and grouped in two distinct clusters by principal component analysis of scores and loadings (SL-plot). The assessment of the overall data behavior (SL-plot, Figure 1B) revealed two distinct clusters, the one corresponding to the CL samples from d 1 to 7 (eCL) (SL-plot, Figure 1B, loadings: C1-C5.2) and the other to CL samples from d 8 to 12 (mCL), d 13-16 (lCL), and d > 18 (rCL). Furthermore the observed clustering of samples explained 99.8% of the variance in our data set, which we argue primarily to be composed by the different expression patterns of bta-miR-21-5p, −143, −202, −182, −96, −181, −10a, and −146a (SL-plot, Figure 1B, Scores). Irrespective of the cycle stage we found a number of miRNAs to be highly expressed in the CL throughout its lifespan (Figures 1C-F). In total, contributed to the top twenty most abundantly sequenced miRNAs in the bovine CL. Each of these miRNAs was also significantly DE in at least one of the individual comparisons between the eCL and three other  Table 3). MiRNA sequences, which align with and hence correspond to bta-miR-21-5p and bta-miR-143-3p account for >40% of the overall reads (Supplementary  Table 2). Moreover, they show opposing and significantly DE at specific CL classes, namely bta-miR-21-5p is significantly upregulated in eCL relative to mCL and lCL (FC = 1.078 and adjusted P = 6.84E-04; FC = 1.119 and adjusted P = 2.71E-04, Supplementary Table 3, eCL vs. mCL, eCL vs. lCL), while bta-miR-143 is significantly downregulated in the eCL vs. the rCL (FC = −1.15 and adjusted P = 3.04E-06, Supplementary Table  3

eCL vs. rCL).
A network of forty-two miRNAs (Figures 2B,C) was identified as consistently and significantly DE in the eCL class relative to their expression at the mCL, lCL, and rCL. Of them five miRNAs (bta-miR-10b, -miR-26a, -miR-30a-5p, -miR-450a, and -miR-202) were as well-classified among the top 20 most abundantly expressed miRNAs. The remaining miRNAs showed sequentially up-or downregulation in eCL vs. mCL, lCL, and rCL classes (Table 1). However seven miRNAs (bta-miR-210-3p, −2898, −96, −7-5p, −183-5p, −182, and −202) show drastic upregulation with a FC ≥ 2.0 between the eCL and the three other CL classes, while bta-miR-146a is the only miRNA found to be differentially downregulated at the same extent in the eCL vs. lCL and rCL (Figures 2B,C, miRNAs with FC ≥ 2.0 and adjusted P < 0.01 are indicated with asterisks). Furthermore all miRNAs were validated by RT-qPCR and the trend of regulation was significant for five of them, while bta-miR-2898 and bta-miR-7 show only a trend in regulation (Table 3, Supplementary Figure 3).
The identification of bta-miR-183-5p, -96, and -182 with increased expression in the eCL vs. the other three cycle stages (Figures 2B,C) is consistent with results from previous studies in bovine (Sontakke et al., 2014;Gebremedhn et al., 2016). All three miRNAs are shown to regulate bovine granulosa cell proliferation by targeting FOXO1 (Gebremedhn et al., 2016). Although, miR-96 was experimentally validated to also target FOXO1 both in human (Guttilla and White, 2009) and bovine (Gebremedhn et al., 2016), its role in bovine CL (luteal cells) has been partially explained. Using human luteinized granulosa cells (hLGC) and bovine CL from day 1 to 4 of the cycle, (Mohammed et al., 2017) identified miR-96 as the key regulator of luteal cell survival by modulating the effect of FOXO1 on steroidogenesis (P4). However, in contrast to hLGC they could not show an inhibitory effect for miR-96 on P4 production in bovine CL. Emphasizing thus the limitation in miR-target interaction validation when using a cross-species approach.
Our sequencing results provide the expression pattern for all three mRNAs: miR-96, -183-5p, and -182 in bovine CL also throughout the cycle and not only for the very early days of the cycle (d 1-4). Interestingly, we found all three miRNAs to peak in their number of reads in the eCL followed by a significant down-regulation in mCL and further steadily upregulation in lCL and rCL (Table 1). This observation suggests miR-96, -183-5p, and -182 to be involved not only during very early CL development in the process of follicle-luteal transition as suggested by Gebremedhn et al. (2016), but also in the steps thereafter when the CL becomes functional (mCL and lCL) and then regresses (rCL). Indeed, a sustained P4 production is required throughout the functional stage of the CL (eCL, mCL, and lCL) and in contrast a decline in P4 production signals a regressing CL (Davis and Rueda, 2002). Since these three miRNAs belong to the same cluster and exhibit similar expression dynamic in bovine CL stages analyzed in this study (Table 1) they presumable regulate other target genes than FOXO1 coordinately to ensure P4 synthesis in the subsequent sates.
Further, we identified bta-miR-202 not only as highly abundant in eCL, mCL, and lCL (Figures 1C-E) but also drastically upregulated in eCL vs. the rCL ( Table 1, FC = 2.3 and adjusted P = 3.88E-10; Table 3, RT-qPCR FC = 0.2 and P = 0.0015). Sontakke et al. (2014), provide the first evidence for the involvement of bta-miR-202 in bovine reproductive tissue and show that bta-miR-202 regulates the viability of mural follicles. Notably, they did a body-wide screening that revealed bta-miR-202 to be expressed exclusively in the gonads. Studies in human have linked mir-202 with different types of cancer acting thus as tumor suppressor or oncogene (reviewed by Lee and Dutta, 2009). In line with the physiology of regressing CL, which is characterized by an enhanced cell death (Stocco et al., 2007) and in contrast to eCL where cell death most be prevented (Sugino and Okuda, 2007), bta-miR-202 peaks in its number of reads in the eCL and the reads fall drastically in the regressing CL ( Table 1, DESeq2 normalized reads in eCL = 31 ′ 881; reads in rCL = 5 ′ 617). The role of bta-miR-202 in the bovine CL is not known, however in the eCL development bta-miR-202 might be acting in similar manner as in cancers: to block cell proliferation, invasion, or apoptosis.
The expression data generated by sequencing and further confirmed by RT-qPCR (Tables 2, 3 and Supplementary Figure 3) highlights for the first time the dynamic regulation of miRNAs over the time of the estrous cycle in bovine CL. Furthermore our results are consistent with findings from profiling studies of circulating miRNAs during the bovine estrous cycle and pregnancy Donadeu, 2016, 2017). Considering that CL's functionality is critical for a physiological cycle and a successful pregnancy, the identification of circulating miRNAs originating from the CL is of particular importance. Recently Ioannidis and Donadeu (2017) identified changes in circulating miRNAs as early as day 8 of the pregnancy in cattle and provide important candidate biomarker for early pregnancy detection. Bta-miR-26a and bta-miR-143, among others were identified as potential biomarkers. Intriguing, bta-miR-143 was shown to be the only miR found in circulation but not originating from the blood, thus letting us hypotheses it could be of endocrine origin and is secreted from the CL. Future studies profiling at the same time circulating miRNA and miRNAs from the CL will facilitated biomarker identification important for early CL function and other CL functional stages.
With regard to the most abundantly expressed miRNAs bta-miR-21-5p is of particular interest since many studies have demonstrated both its up-regulation in response to postovulatory LH surge (Fiedler et al., 2008), as well as its regulatory effect on different apoptosis genes to promote cell survival (Asangani et al., 2008;Carletti et al., 2010;Buscaglia and Li, 2011). Survival of the ovarian granulosa cells is decisive for their transition into luteal cells and along with the process of follicle basal membrane cells degradation it enables angiogenesis and early CL development (Reynolds and Redmer, 1999;Skarzynski et al., 2013). The cellular regulation of both of these events corresponds to the highest expression level of bta-miR-21-5p, specifically in the early CL ( Figure 1C, DESeq2 normalized reads = 320 × 10 3 ). This observation and published data from mice studies (Fiedler et al., 2008;Carletti et al., 2010) let us to hypothesize bta-miR-21-5p had an indirect effect on bovine eCL development. Initially, it is the inflammatory like response elicited after the LH (Richards et al., 2008) surge that is shown to increase miR-21 expression from pre-miR-21 via the IL6 and transcription factor STAT3 (Löffler et al., 2007;Han et al., 2012) and following its abundance that facilitates angiogenesis via the suppression of the reversion-inducing, cysteine-rich protein with kazal motifs (RECK) gene, which is shown to negatively regulate transcription of MMP1 and MMP9 (Takagi et al., 2009;Han et al., 2012). MMPs activity is required for cellular degradation of basal membrane prior to angiogenesis (Oh et al., 2001). Both of these MMPs are shown by our group (Kliem et al., 2007) to be highly expressed in bovine in the eCL (d 1-7) and mid . Accordingly, we propose that the increase in STAT3 mediated bta-miR-21-5p expression leads to downregulation of apoptosis genes (Buscaglia and Li, 2011) and the RECK mRNA levels, allowing thus granulose cell transition in luteal cells and MMPs to be active for the degradation of basal membrane cells and propagation of new blood vessel form existing ones (angiogenesis). Furthermore, we identified bta-miR-146a as drastically downregulated in eCL vs. the mCL and lCL. Bta-miR-146a has just recently been shown to act as a negative feedback regulator of bovine inflammation and innate immunity through downregulation of the TLR4/TRAF6/NF-κB pathway (Wang et al., 2017) and we propose bta-miR-146a supports the activity of bta-miR-21-5p by downregulating the post-ovulation inflammatory mediators such as IL6 and IL8 in the early CL in bovine.
Bta-miR-7-5p and bta-miR-2898 are also identified as drastically upregulated in eCL vs. the mCL and lCL, rCL respectively (Figures 2B,C). Studies implicating bta-miR-7-5p in bovine are limited to its presence and identification in immunerelated tissues (Coutinho et al., 2007) and to date no functional validation studies are performed in bovine or other species. Whereas, recent studies show bta-miR-2898 expression to be upregulated in mammary gland tissues of mastitis-infected cows (Wang et al., 2014). This implicates bta-miR-2898 in the immune response relevant to mastitis, however in the miRTarBase experimentally validated interactions for this miRNA are not provided. Consequently its role in bovine CL needs further investigation.
CL formation has similarities to the processes involved with wound healing and tumorogenesis, both of which that are enabled by extensive angiogenesis (Niswender et al., 2000). Interestingly, in mice miR-21 has also been shown to induce tumor angiogenesis by targeting PTEN and activating AKT and ERK1/2 signaling pathways, which then in turn enhances VEGF and hypoxia-inducible factor 1-alpha (HIF1A) expression . Our group showed that in bovine CL, the highest expression of VEGF and VEGFR-2 is detected during days 1-7 (eCL) of the estrous cycle (Berisha et al., 2000b(Berisha et al., , 2017 and that is concurrent with increased bta-miR-21-5p expression we observed in the eCL (Figure 1C, d 1-7). Thus, similar to studies with mice  and because of its high abundance in the eCL, bta-miR-21-5p might activate AKT and ERK1/2 signaling pathways in bovine CL and promote angiogenesis through enhancement of VEGF expression at that time.
Around day 8 of the cycle, the hormonal environment becomes optimal and the CL produces P4 at levels sufficient to maintain the cellular phenotype at the mid CL stage (Wiltbank, 1994). However, shortly after this and in the absence of pregnancy the CL becomes responsive to prostaglandin F2 alpha (PGF2A) and starts regressing (Skarzynski et al., 2013). The regression of bovine CL is a complex process that starts with a decline in the synthesis of P4 (functional luteolysis), followed by tissue involution and cell death (Berisha et al., 2010). It involves multiple different factors, including IFNG, TNF, ANGPT2, IGF1, STAT3, SMAD2/4, TGFb1, IL8, and IL18 (Sugino and Okuda, 2007;Hou et al., 2008;Berisha et al., 2010). Most of these factors are experimentally validated gene targets for human miRNAs (Table 3), however very few are also in silico predicted for bovine miRNAs ( Table 3, miRNAs in bold).
The regression of CL is accompanied by extensive extracellular matrix remodeling and deposition of collagen type I (COL1A1) (Irving-Rodgers et al., 2006). Regression of the CL is likely to be facilitated by the second most abundantly expressed miRNA identified in this study bta-miR-143, which is significantly upregulated in rCL relative to the eCL ( Figure 1F, Supplementary  Table 3C) and COL1A1 is an experimentally validated target gene for has-miR-143 and successfully predicted for bta-miR-143 (Table 3). Further miR-143 targets genes related to P4 synthesis and cell death including COX2, IGF1R, IFNG, TNFA, and Bcl-2 (Table 3). These genes are relevant for CL regression, as the abundance of most of them is altered upon PGF2A administration. Some of them, specifically TNFA, act as intraovary factors and induce apoptotic death by modulating the Bcl-2 family of genes, and caspase 3 (Skarzynski et al., 2005;Kliem et al., 2009). Furthermore miR-143 has been shown to promote porcine granulose cell apoptosis and follicular atresia by suppression the follicle stimulating hormone receptor (FSHR) (Du et al., 2016).
Another interesting finding of this study follows from the self-organizing tree analysis (SOTA). Using SOTA we showed that differentially expressed miRNAs across the 21 day cycle demonstrate a specific pattern of expression. Based on a correlation of their expression, miRNAs were grouped in eight distinct clusters (Figures 3A,B). The clusters themselves were furthermore clearly split into two groups, one consisting of the miRNAs upregulated at the beginning of the cycle (Figure 3A, bottom 4 panel), and the other of miRNAs upregulated in the mid-, late-, and regressed CL stages (Figure 3A, top 4 panel). We theorize that this observation is due to "progesterone induced miRNAs expression." MiRNAs in cluster 7 (miR-29b, −31, −210, and let-7a) for instance are shown to be increased with progestin in human glandular epithelium derived from the uterus (Kuokkanen et al., 2010). However, despite this, to date no progestin response elements have been identified in the promoter region of miRNAs shown to be regulated by progestin. Therefore, the data should be viewed with this in mind and also by considering the other non-genomic means (Kuse et al., 2015), by which P4 can induce miRNA expression and may also be considered as determining factor for the cycle stage dependent miRNA expression. Notable examples of this include miRNAs having been shown to stably circulate in the blood in exosomes and therefore that they may also regulate gene expression in a paracrine manner (Zhang et al., 2015). There is much to be discovered about this exciting mechanism but it may represent one mode of action for miRNAs in the regulation of this dynamic gland.
The miRNAs we identified as differentially expressed only between individual comparisons of eCL to mCL are likely targeting genes involved in fatty acid biosynthesis, as well as TGFB-, Wnt-, and FoxO-signaling ( Figure 4A). Further, pathways enriched with genes targeted by miRNA DE between only between eCL and lCL included the glutamatergic synapse, axon guidance, thyroid hormone signaling pathway, oxytocin and prolactin signaling ( Figure 4B). Pathway enrichment for genes targeted by miRNAs DE only between eCL and rCL revealed the ECM-receptor interaction as the most enriched pathway, followed by lysine degradation, glycosphingolipid biosynthesis-lacto and neolacto series and the focal adhesion signaling pathway ( Figure 4C). Which individual miRNAs contributed to pathway enrichment is accessible as additional file to this paper.
Overall, the pathways we see to be enriched for genes targeted by miRNA DE only between individual compression are known to be involved in the early luteal CL development, formation and luteinisation of CL.

CONCLUSION
This study utilizes a small RNA sequencing technology to temporally profile miRNA transcriptome in the CL from cycling cows. It is the first study, to our knowledge, to analyze the dynamics of miRNAs expression at six different timely defined CL classes, throughout the entire 21 days of the bovine estrous cycle.
It enabled the identification of bta-miR-21-5p and bta-miR-143 as two species that were highly abundant regardless of CL developmental or functional stage. We demonstrated that 42 miRNAs are consistent significantly different expressed at mCL, lCL, and rCL compared to their expression in the eCL. Out of them bta-miR- 210-3p, -2898, -96, -7-5p, -183-5p, -182, and -202 show drastic up-regulation between the eCL and three subsequent CL classes, while bta-miR-146a is the only miRNA found to be differentially downregulated at the same extent in the eCL vs. lCL and rCL. These miRNAs have predicted targets involved in regulating critical angiogenic and antiangiogenic, as well as luteotropic and luteolytic genes.
In addition, we identified groups of miRNAs DE only at individual comparison's: 24 miRNAs in eCL vs. mCL only, 11 miRNAs in eCL vs. lCL only and 21 miRNAs in eCL vs. rCL only. Each of these groups of DE miRNAs was found to regulate distinct pathways involved in the progression of the CL from early development, to fully functional and regressing phase.
Using the self-organizing tree algorithm (SOTA) we grouped all DE miRNAs into eight clusters based on their expression profiles. These clusters were clearly split into two larger groups, consisting of miRNAs up regulated at either the beginning of the cycle or the mid-, late-and regressed stages of the CL cycle.
Considering the very limited access to healthy human CL tissue for research reasons, the similarities between early CL development and tumorogensis; and the similarity in ovarian physiology between cattle and human, this study is of particular importance and provides novel insights in the dynamics of miRNA transcriptome during the physiological estrous cycle.

AVAILABILITY OF DATA AND MATERIALS
The sequencings files are deposited at the ENA database and are available under the accession number: PRJEB20028.

AUTHOR CONTRIBUTIONS
RMG and BB: Conceived of the study and participated in its design; RMG: Wrote the manuscript. RMG and CIS: Performed the library preparation for Illumina HiSeq2500 sequencing platform; BK, RMG, and CIS: Conducted statistical analysis; BB, MWP, IR, and RT-C: Reviewed the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was supported by the German Research Foundation (DFG BE 3189/2-2) and the Ministry of Education, Science and Technology, Pristina, Kosovo (FPVSh -2-4051-2).