Rreze M. Gecaj1
Corina I. Schanzenbach2,3
Benedikt Kirchner2,4Michael W. Pfaffl2Irmgard Riedmaier5,6
Ry Y. Tweedie-Cullen7
Bajram Berisha1,2*- 1Department of Animal Husbandry and Biotechnology, Faculty of Agriculture and Veterinary, University of Pristina, Pristina, Albania
- 2Department of Animal Physiology & Immunology, Weihenstephan, Technical University of Munich, Munich, Germany
- 3Department of Environmental Systems Science, ETH Zürich, Zurich, Switzerland
- 4Dr. von Hauner Children's Hospital, Ludwig Maximillian Universität, München, Germany
- 5Department of Physiology, TUM School of Life Sciences Weihenstephan, Technical University of Munich, Freising, Germany
- 6Eurofins Medigenomix Forensik GmbH, Ebersberg, Germany
- 7Department of Health Science and Technology, ETH Zürich, Zurich, Switzerland
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.
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, 2002, 2004), and CL regression (Reynolds et al., 2000; Neuvians et al., 2004; Berisha et al., 2008, 2010, 2016, 2017; Miyamoto 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,b, 2002, 2004, 2008, 2010, 2016, 2017; Neuvians et al., 2003, 2004; Kliem 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 (Dicerd/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 (Yao et al., 2010).
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; Li et al., 2011; Sontakke et al., 2014) or in the CL from a specific stage (Maalouf et al., 2014, 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.
Materials and Methods
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, 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, 2015; Buschmann 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.
Reverse Transcription and Real-Time PCR
The miScript PCR assays from Qiagen were used for the validation of the NGS results. The regulated miRNAs were quantified from 100 ng of the total RNA (n = 18 animals) and reversed transcribed using the 5x miScript HighSpec buffer. The RT-PCR was carried out in a Rotor—Disc 100 on the Rotor-Gene Q, Real-Time Detection System. For validation the following miScript assays were used: Bt_mir-210_1, Bt_mir-96_1, Bt_mir-146a_1, Bt_mir-182_1, Bt_mir-183_1, Bt_mir-2898_1, Bt_mir-7_1, Bt_mir-202_1, Bt_mir-21-5p_1, and Bt_mir-143_1.
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 miRNA-target 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).
Results
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.
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.
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 (Table 1, significantly downregulated in mCL, lCL, and rCL vs. the eCL,) and miR-96 (Table 1, significantly upregulated in mCL, lCL, and rCL vs. the eCL). Nevertheless the cumulative effect of other DE miRNAs (Table 1) also counts for the cluster separation of our data.
Table 1. Differentially expressed bovine miRNAs throughout the CL lifespan.
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 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).
Figure 2. Distribution of differentially expressed (DE) miRNAs across all comparisons and log2 fold-changes of miRNAs consistently DE in all three CL classes vs. the eCL. (A) Venn diagram showing the overall number of DE miRNAs and their relation to each of the analyzed CL classes being: early CL (eCL) vs. mid CL (mCL), late CL (lCL), and regressed CL (rCL), respectively. (B) log2 fold-changes of miRNAs with <1,000 reads in the eCL. (C) log2 fold-changes of miRNAs with >1,000 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.
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, Table 1, FC >1). The remaining 21 miRNAs were downregulated for the same comparison (Figures 2B,C, Table 1, FC < 1). Yet, seven of them: bta-miR-210, −2898, −96, −7-5p, −183-5p, −182, and −202 show a fold change (FC) up-regulation ≥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 ≥ 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, Table 1). Other significantly regulated miRNAs with a FC below 2.0 can be found in the Table 1. MiRNAs upregulated in eCL vs. the three other stages are in bold only with a positive FC, while miRNAs downregulated in the eCL vs. all three other stages are indicated with a negative FC.
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), of which 21 are unique to this comparison (Figure 2A, Supplementary Table 4B). Details of all regulated miRNAs with their respective log2FC and their adjusted P-values across all comparisons can be accessed in the supplementary data (Supplementary Tables 3A–C), while the distribution of common and unique miRNAs DE between the eCL and other CL classes is provided as Supplementary Tables 4A,B.
RT-qPCR and miRNA Validation
Eight miRNAs (bta-miR-210, -96, -7, -182, -183, -2898, -146a, and -202) that were found to be DE regulated across the different CL classes as well as two miRNAs (bta-miR-21-5p and bta-miR-143) identified as the most abundantly expressed miRNA species are chosen for RT-qPCR validation.
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.
Table 2. Fold-change regulation of DE miRNAs during the estrous cycle obtained by RT-qPCR and NGS.
Table 3. Selected candidate genes related to CL development, function, and regression that are targeted by DE miRNAs common to all: eCL vs. mCL, lCL, and rCL.
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 (Ioannidis and 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.
Figure 3. Clusters of mean centered Illumina read counts (log2 scale) with different miRNA expression profiles during early CL 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) Expression graphs of all clusters, blue: down regulation, red: up regulation, pink line: mean expression. (B) Detailed unsupervised hierarchical clustering (HCL) dendrogram generated by Pearson Corelation Coeficient. Rows indicate single DE miRNAs and columns show CL classes across the oestrous cycle. (C) Dendrogram including hierarchical clustering of the eight clusters and mean expression values for each cycle stage. SOTA and HCL are performed using MeV.
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).
Figure 4. KEGG pathways enriched with genes targeted by more than 3 DE miRNAs from the comparison of individual CL classes. (A) Pathways enriched with genes targeted by differentially regulated miRNAs in either: eCL vs. mCL, eCL vs. lCL or eCL vs. rCL. Top 15 most significantly enriched pathways have been plotted. Lower p-values indicate more significance in terms of the biological function targeted by miRNAs in the pathway. (B) List of miRNAs that contributed to the 5 most enriched pathways between the comparisons: eCL vs. mCL, eCL vs. lCL, and eCL vs. rCL.
“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). Furthermore, genes involved with remodeling of the extracellular matrix (ECM) were also amongst the top predicted targets of the DE miRNAs. Examples include thrombospondin 1 (THBS1), tissue inhibitor of metalloproteinase 1 (TIMP-1), matrix-degrading proteinase 1 (MMP-1), and matrix-degrading proteinase 9 (MMP-9). Multiple cell cycle related genes were also seen as targets: phosphatase and tensin homolog (PTEN), cyclin-dependent kinase inhibitor 1B (CDKN1B), cyclin-dependent kinase inhibitor 1A (CDKN1A), cyclin D2 (CCND2), and cell division protein kinase 6 (CDK6) as well as genes involved in the CL regression (see Table 3).
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 Wnt-signaling pathway; between: eCL and lCL, the glutamatergic synapse, axon guidance, thyroid hormone signaling, oxytocin and prolactin synthesis and between: eCL and rCL, the ECM-receptor 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, 22 miRNAs (bta-miR-21-5p, -miR-143, -miR-10b, -let-7i, -miR-202, -miR-148a, -let-7f, -miR-3600, -miR-99a55p, -let-7a-5p, -miR-27b, -miT-100, -let7g, -miR-26a, -miR-378, -miR-30d, -miR-125b, -450a, -miR-30e-5p, -let-7b, -miR-199a-3p, and -miR-26c) 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 CL classes (Supplementary 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 up-regulation 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 up-regulation 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 (Ioannidis and 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 × 103). 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 CL (d 8–12). 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 immune-related 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 (Liu et al., 2011). 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, 2017 unpublished data) 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 (Liu et al., 2011) 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.
In contrast, the elevated expression of HIF1A, which has been documented by our group to be the highest in bovine CL during d 1–7 (Berisha et al., 2017) may in turn regulate the expression of other miRNAs and likely supports the activity of bta-miR-21-5p. In agreement, bta-miR-210-3p one of the drastically upregulated miRNA in the eCL vs. mCL (FC = 2.4, adjusted P = 2.26E-16, Figure 2B, Table 1), was shown to be induced under hypoxic conditions (Ivan and Huang, 2014). Hypoxic conditioning is documented in the early CL of monkeys (Tesone et al., 2005) and hypoxia is critical for establishing vascularization in bovine CL (Nishimura and Okuda, 2010). Furthermore, HIF1A is an experimentally validated target gene for human micro-210-3p (hsa-mir-210-3p) (Table 3), but not for bta-miR-210-3p (http://www.targetscan.org). Nevertheless, HIF1A is an experimentally validated target for has-miR-17-5p and -20a (http://mirtarbase.mbc.nctu.edu.tw/) and is also successfully in silico predicted for bta-miR-17-5p and -20a (http://www.targetscan.org). Notable, bta-miR-17-5p and bta-miR-20a are DE miRNAs identified in our study and show a FC < 2.0 (Table 1).
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 intra-ovary 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).
Conflict of Interest Statement
The author IR was employed by the company Eurofins Medigenomix Forensik GmbH. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
The authors thank Mrs. Angela Sachsenhauser at the Animal Physiology and Immunology, TUM for her valuable support during the sample preparation.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2017.00213/full#supplementary-material
References
Anders, S., and Huber, W. (2010). DESeq: differential expression analysis for sequence count data. Genome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106
Asangani, I. A., Rasheed, S. A., Nikolova, D. A., Leupold, J. H., Colburn, N. H., Post, S., et al. (2008). MicroRNA-21 (miR-21) post-transcriptionally downregulates tumor suppressor Pdcd4 and stimulates invasion, intravasation and metastasis in colorectal cancer. Oncogene 27, 2128–2136. doi: 10.1038/sj.onc.1210856
Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi: 10.1016/S0092-8674(04)00045-5
Berisha, B., Meyer, H. H. D., and Schams, D. (2010). Effect of prostaglandin F2 alpha on local luteotropic and angiogenic factors during induced functional luteolysis in the bovine corpus luteum. Biol. Reprod. 82, 940–947. doi: 10.1095/biolreprod.109.076752
Berisha, B., Schams, D., Rodler, D., Sinowatz, F., and Pfafl, M. W. (2016). Expression and localization of members of the thrombospondin family during final follicle maturation and corpus luteum formation and function in the bovine ovary. J. Reprod. Dev. 62, 501–510. doi: 10.1262/jrd.2016-056
Berisha, B., Schams, D., Kosmann, M., Amselgruber, W., and Einspanier, R. (2000a). Expression and localisation of vascular endothelial growth factor and basic fibroblast growth factor during the final growth of bovine ovarian follicles. J. Endocrinol. 167, 371–382. doi: 10.1677/joe.0.1670371
Berisha, B., Schams, D., Kosmann, M., Amselgruber, W., and Einspanier, R. (2000b). Expression and tissue concentration of vascular endothelial growth factor, its receptors, and localization in the bovine corpus luteum during estrous cycle and pregnancy. Biol. Reprod. 63, 1106–1114. doi: 10.1095/biolreprod63.4.1106
Berisha, B., Schams, D., and Miyamoto, A. (2002). The expression of angiotensin and endothelin system members in bovine corpus luteum during estrous cycle and pregnancy. Endocrine 19, 305–312. doi: 10.1385/ENDO:19:3:305
Berisha, B., Schams, D., Rodler, D., Sinowatz, F., and Pfaffl, M. W. (2017). Expression pattern of HIF1alpha and vasohibins during follicle maturation and corpus luteum function in the bovine ovary. Reprod. Domest. Anim. 52, 130–139. doi: 10.1111/rda.12867
Berisha, B., Sinowatz, F., and Schams, D. (2004). Expression and localization of fibroblast growth factor (FGF) family mambers during the final growth of bovine ovarian follicles. Mol. Reprod. Dev. 67, 162–171. doi: 10.1002/mrd.10386
Berisha, B., Steffl, M., Welter, H., Kliem, H., Meyer, H. H. D., Schams, D., et al. (2008). Effect of the luteinising hormone surge on regulation of vascular endothelial growth factor and extracellular matrix-degrading proteinases and their inhibitors in bovine follicles. Reprof. Fertil. Dev. 20, 258–268. doi: 10.1071/RD07125
Buscaglia, L. E. B., and Li, Y. (2011). Apoptosis and the target genes of microRNA-21. Chin. J. Cancer 30, 371–380. doi: 10.5732/cjc.30.0371
Buschmann, D., Haberberger, A., Kirchner, B., Spornraft, M., Riedmaier, I., Schelling, G., et al. (2016). Toward reliable biomarker signatures in the age of liquid biopsies - How to standardize the small RNA-Seq workflow. Nucleic Acids Res. 44, 5995–6018. doi: 10.1093/nar/gkw545
Carletti, M. Z., Fiedler, S. D., and Christenson, L. K. (2010). MicroRNA 21 blocks apoptosis in mouse periovulatory granulosa cells. Biol. Reprod. 83, 286–295. doi: 10.1095/biolreprod.109.081448
Chou, C.-H., Chang, N.-W., Shrestha, S., Hsu, S.-D., Lin, Y.-L., Lee, W.-H., et al. (2016). miRTarBase: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 44, D239–D247. doi: 10.1093/nar/gkv1258
Christenson, L. K., Gunewardena, S., Hong, X., Spitschak, M., Baufeld, A., and Vanselow, J. (2013). Preovulatory LH surge effects on follicular theca and granulosa transcriptomes. Mol. Endocrinol. 27, 1153–1171. doi: 10.1210/me.2013-1093
Coutinho, L. L., Matukumalli, L. K., Sonstegard, T. S., Van Tassell, C. P., Gasbarre, L. C., Capuco, A. V., et al. (2007). Discovery and profiling of bovine microRNAs from immune-related and embryonic tissues. Physiol. Genomics 29, 35–43. doi: 10.1152/physiolgenomics.00081.2006
Davis, J. S., and Rueda, B. R. (2002). The corpus luteum: an ovarian structure with maternal instincts and suicidal tendencies. Front. Biosci. 7:d1949–d1978. doi: 10.2741/davis1
Davis, J. S., Rueda, B. R., and Spanel-Borowski, K. (2003). Microvascular endothelial cells of the corpus luteum. Reprod. Biol. Endocrinol. 1:89. doi: 10.1186/1477-7827-1-89
Dopazo, J., and Carazo, J. M. (1997). Phylogenetic reconstruction using an unsupervised growing neural network that adopts the topology of a phylogenetic tree. J. Mol. Evol. 44, 226–233. doi: 10.1007/PL00006139
Du, X., Zhang, L., Li, X., Pan, Z., Liu, H., and Li, Q. (2016). TGF-β signaling controls FSHR signaling-reduced ovarian granulosa cell apoptosis through the SMAD4/miR-143 axis. Cell Death Dis. 7:e2476. doi: 10.1038/cddis.2016.379
Ewing, B., and Green, P. (1998). Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Res. 8, 186–194. doi: 10.1101/gr.8.3.186
Fiedler, S. D., Carletti, M. Z., Hong, X., and Christenson, L. K. (2008). Hormonal regulation of MicroRNA expression in periovulatory mouse mural granulosa cells. Biol. Reprod. 79, 1030–1037. doi: 10.1095/biolreprod.108.069690
Fleige, S., and Pfaffl, M. W. (2006). RNA integrity and the effect on the RT-PCR performance. Mol. Aspects Med. 27, 126–139. doi: 10.1016/j.mam.2005.12.003
Friedman, R. C., Farh, K. K. H., Burge, C. B., and Bartel, D. P. (2009). Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 19, 92–105. doi: 10.1101/gr.082701.108
Gebremedhn, S., Salilew-Wondim, D., Hoelker, M., Rings, F., Neuhoff, C., Tholen, E., et al. (2016). MicroRNA-183-96-182 cluster regulates bovine granulosa cell proliferation and cell cycle transition by coordinately targeting FOXO1. Biol. Reprod. 94:127. doi: 10.1095/biolreprod.115.137539
Guttilla, I. K., and White, B. A. (2009). Coordinate regulation of FOXO1 by miR-27a, mir-96, and miR-182 in breast cancer cells. J. Biol. Chem. 284, 23204–23216. doi: 10.1074/jbc.M109.031427
Han, L., Yue, X., Zhou, X., Lan, F. M., You, G., Zhang, W., et al. (2012). MicroRNA-21 Expression is regulated by b-catenin/STAT3 pathway and promotes glioma cell invasion by direct targeting RECK. CNS Neurosci. Ther. 18, 573–583. doi: 10.1111/j.1755-5949.2012.00344.x
Herrero, J., Valencia, A., and Dopazo, J. (2001). A hierarchical unsupervised growing neural network for clustering gene expression patterns. Bioinformatics 17, 126–136. doi: 10.1093/bioinformatics/17.2.126
Hossain, M. M., Ghanem, N., Hoelker, M., Rings, F., Phatsara, C., Tholen, E., et al. (2009). Identification and characterization of miRNAs expressed in the bovine ovary. BMC Genomics 10:443. doi: 10.1186/1471-2164-10-443
Hou, X. Y., Arvisais, E. W., Jiang, C., Chen, D. B., Roy, S. K., Pate, J. L., et al. (2008). Prostaglandin F2 alpha stimulates the expression and secretion of transforming growth factor B1 via induction of the early growth response 1 gene (EGR1) in the bovine corpus luteum. Mol. Endocrinol. 22, 403–414. doi: 10.1210/me.2007-0272
Hulsen, T., de Vlieg, J., and Alkema, W. (2008). BioVenn-a web application for the comparison and visualization of biological lists using area-proportional Venn Diagrams. BMC Genomics 9:488. doi: 10.1186/1471-2164-9-488
Illumina (2011). Quality Scores for Next-Generation Sequencing. Technical note, Quality scores for Next Generation Sequencing. 1–2.
Ioannidis, J., and Donadeu, F. X. (2016). Circulating microRNA profiles during the Bovine Oestrous Cycle. PLoS ONE 11:e0158160. doi: 10.1371/journal.pone.0158160
Ioannidis, J., and Donadeu, F. X. (2017). Changes in circulating microRNA levels can be identified as early as day 8 of pregnancy in cattle. PLoS ONE 12:e0174892. doi: 10.1371/journal.pone.0174892
Irving-Rodgers, H. F., Friden, B. E., Morris, S. E., Mason, H. D., Brannstrom, M., Sekiguchi, K., et al. (2006). Extracellular matrix of the human cyclic corpus luteum. Mol. Hum. Reprod. 12, 525–534. doi: 10.1093/molehr/gal060
Ivan, M., and Huang, X. (2014). miR-210:fine-tuning the hypoxic response. Adv. Exp. Med. Biol. 772, 205–227. doi: 10.1007/978-1-4614-5915-6_10
Kliem, H., Berisha, B., Meyer, H. H. D., and Schams, D. (2009). Regulatory changes of apoptotic factors in the bovine corpus luteum after induced luteolysis. Mol. Reprod. Dev. 76, 220–230. doi: 10.1002/mrd.20946
Kliem, H., Welter, H., Kraetzl, W. D., Steffl, M., Meyer, H. H. D., Schams, D., et al. (2007). Expression and localisation of extracellular matrix degrading proteases and their inhibitor during the oestrous cycle and after induced luteolysis in the bovine corpus luteum. Reproduction 134, 535–547. doi: 10.1530/REP-06-0172
Kloosterman, W. P., and Plasterk, R. H. A. (2006). The diverse functions of microRNAs in animal development and disease. Dev. Cell 11, 441–450. doi: 10.1016/j.devcel.2006.09.009
Kong, Y. (2011). Btrim: a fast, lightweight adapter and quality trimming program for next-generation sequencing technologies. Genomics 98, 152–153. doi: 10.1016/j.ygeno.2011.05.009
Kozomara, A., and Griffiths-Jones, S. (2014). MiRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi: 10.1093/nar/gkt1181
Kuokkanen, S., Chen, B., Ojalvo, L., Benard, L., Santoro, N., and Pollard, J. W. (2010). Genomic profiling of microRNAs and messenger RNAs reveals hormonal regulation in microRNA expression in human endometrium. Biol. Reprod. 82, 791–801. doi: 10.1095/biolreprod.109.081059
Kuse, M., Sakumoto, R., and Okuda, K. (2015). Genomic and non-genomic effects of progesterone on prostaglandin (PG) F2alpha and PGE2 production in the bovine endometrium. Reprod. Fertil. Dev. 28, 1588–1597. doi: 10.1071/RD14490
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
Lee, I.-H., Hong, X., Mathur Sh, C., Sharma, M., Rastogi, A., Sharma, P., et al. (2014). A detailed analysis of next generation sequencing reads of microRNA expression in Barrett Esophagus: absolute versus relative quantification. BMC Res. Notes 7:212. doi: 10.1186/1756-0500-7-212
Lee, R. C., Feinbaum, R. L., and Ambros, V. (1993). The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell 75, 843–854.
Lee, Y. S., and Dutta, A. (2009). MicroRNAs in cancer. Ann. Rev. Pathol. 4, 199–227. doi: 10.1146/annurev.pathol.4.110807.092222
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Li, M., Liu, Y., Wang, T., Guan, J., Luo, Z., Chen, H., et al. (2011). Repertoire of porcine MicroRNAs in adult ovary and testis by deep sequencing. Int. J. Biol. Sci. 7, 1045–1055. doi: 10.7150/ijbs.7.1045
Liu, L. Z., Li, C., Chen, Q., Jing, Y., Carpenter, R., Jiang, Y., et al. (2011). Mir-21 induced angiogenesis through AKT and ERK activation and HIF-1a expression. PLoS ONE 6:e19139. doi: 10.1371/journal.pone.0019139
Löffler, D., Brocke-Heidrich, K., Pfeifer, G., Stocsits, C., Hackermüller, J., Kretzschmar, A. K., et al. (2007). Interleukin-6-dependent survival of multiple myeloma cells involves the Stat3-mediated induction of microRNA-21 through a highly conserved enhancer. Blood 110, 1330–1333. doi: 10.1182/blood-2007-03-081133
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 1–34. doi: 10.1186/s13059-014-0550-8
Maalouf, S. W., Liu, W. S., Albert, I., and Pate, J. L. (2014). Regulating life or death: potential role of microRNA in rescue of the corpus luteum. Mol. Cell. Endocrinol. 398, 78–88. doi: 10.1016/j.mce.2014.10.005
Maalouf, S. W., Smith, C. L., and Pate, J. L. (2016). Changes in microRNA expression during maturation of the bovine corpus luteum: regulation of luteal cell proliferation and function by microRNA-34a. Biol. Reprod. 94, 71. doi: 10.1095/biolreprod.115.135053
Miyamoto, A., Shirasuna, K., Shimizu, T., Bollwein, H., and Schams, D. (2010). Regulation of corpus luteum development and maintenance: specific roles of angiogenesis and action of prostaglandin F2alpha. Soc. Reprod. Fertil. Suppl. 67, 289–304.
Mohammed, B. T., Sontakke, S. D., Ioannidis, J., Duncan, W. C., and Donadeu, F. X. (2017). The adequate corpus luteum: miR-96 promotes luteal cell survival and progesterone production. J. Clin. Endocrinol. Metab. 102, 2188–2198. doi: 10.1210/jc.2017-00259
Mueller, O., and Schroeder, A. (2004). RNA Integrity Number (RIN) – standardization of RNA quality control application. Nano 1–8.
Neuvians, T. P., Berisha, B., and Schams, D. (2004). Vascular Endothelial Growth Factor (VEGF) and Fibroblast Growth Factor (FGF) Expression during Induced Luteolysis in the Bovine Corpus Luteum. Mol. Reprod. Dev. 67, 389–395. doi: 10.1002/mrd.20032
Neuvians, T. P., Pfaffl, M. W., Berisha, B., and Schams, D. (2003). The mRNA expression of the members of the IGF-system in bovine corpus luteum during induced luteolysis. Domest. Anim. Endocrinol. 25, 359–372. doi: 10.1016/j.domaniend.2003.08.002
Nishimura, R., and Okuda, K. (2010). Hypoxia is important for establishing vascularization during corpus luteum formation in cattle. J. Reprod. Dev. 56:110–116. doi: 10.1262/jrd.09-162E
Niswender, G. D., Juengel, J. L., Silva, P. J., Rollyson, M. K., and McIntush, E. W. (2000). Mechanisms controlling the function and life span of the corpus luteum. Physiol. Rev. 80, 1–29.
Oh, J., Takahashi, R., Kondo, S., Mizoguchi, A., Adachi, E., Sasahara, R. M., et al. (2001). The membrane-anchored MMP inhibitor RECK is a key regulator of extracellular matrix integrity and angiogenesis. Cell 107, 789–800. doi: 10.1016/S0092-8674(01)00597-9
Otsuka, M., Zheng, M., Hayashi, M., Lee, J. D., Yoshino, O., Lin, S., et al. (2008). Impaired microRNA processing causes corpus luteum insufficiency and infertility in mice. J. Clin. Invest. 118, 1944–1954. doi: 10.1172/JCI33680
Pfaffl, M. W. (2001). A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 29:e45. doi: 10.1093/nar/29.9.e45
Rekawiecki, R., Kowalik, M. K., Slonina, D., and Kotwica, J. (2008). Regulation of progesterone synthesis and action in bovine corpus luteum. J. Physiol. Pharmacol. 59 (Suppl. 9), 75–89.
Reynolds, L. P., and Redmer, D. A. (1999). Growth and development of the corpus luteum. J. Reprod. Fertil. Suppl. 54, 181–191.
Reynolds, L. P., Grazul-Bilska, A. T., and Redmer, D. A. (2000). Angiogenesis in the corpus luteum. Endocrine 12, 1–9. doi: 10.1385/ENDO:12:1:1
Richards, J. S., Liu, Z., and Shimada, M. (2008). Immune-like mechanisms in ovulation. Trends Endocrinol. Metab. 19, 191–196. doi: 10.1016/j.tem.2008.03.001
RNAcentral Consortium (2015). RNAcentral: an international database of ncRNA sequences. Nucleic Acids Res. 43, D123–D129. doi: 10.1093/nar/gku991
Saeed, A. I., Sharov, V., White, J., Li, J., Liang, W., Bhagabati, N., et al. (2003). TM4: a free, open-source system for microarray data management and analysis. Biotechniques 34, 374–378.
Sarkar, M., Schilffarth, S., Schams, D., Meyer, H. H. D., and Berisha, B. (2011). The expression of thrombopoietin and its receptor during different physiological stages in the bovine ovary. Reprod. Domest. Anim. 46, 757–762. doi: 10.1111/j.1439-0531.2010.01736.x
Schams, D., and Berisha, B. (2004). Regulation of corpus luteum function in cattle–an overview. Reprod. Domest. Anim. 39, 241–251. doi: 10.1111/j.1439-0531.2004.00509.x
Schanzenbach, C. I., Kirchner, B., Ulbrich, S. E., and Pfaffl, M. W. (2017). Can milk cell or skim milk miRNAs be used as biomarker for early pregnancy detection in cattle? PLoS ONE 12:e0172220. doi: 10.1371/journal.pone.0172220
Shimizu, T., Berisha, B., Schams, D., and Miyamoto, A. (2007). Expression of angiopoietin (ANPT)-1, ANPT-2 and their receptors in dominant follicles during periovulatory period in GnRH-treated cow. Reprod. Domest. Anim. 42, 221–224. doi: 10.1111/j.1439-0531.2006.00748.x
Shivdasani, R. A. (2006). MicroRNAs: regulators of gene expression and cell differentiation. Blood 108, 3646–3653. doi: 10.1182/blood-2006-01-030015
Skarzynski, D. J., Jaroszewski, J., and Okuda, K. (2005). Role of tumor necrosis factor-alpha and nitric oxide in luteolysis in cattle. Domest. Anim. Endocrinol. 29, 340–346. doi: 10.1016/j.domaniend.2005.02.005
Skarzynski, D., Piotrowska-Tomala, K., Lukasik, K., Galvão, A., Farberov, S., Zalman, Y., et al. (2013). Growth and regression in bovine corpora lutea: regulation by local survival and death pathways. Reprod. Domest. Anim. 48(Suppl. 1), 25–37. doi: 10.1111/rda.12203
Sontakke, S. D., Mohammed, B. T., McNeilly, A. S., and Donadeu, F. X. (2014). Characterization of microRNAs differentially expressed during bovine follicle development. Reproduction 148, 271–283. doi: 10.1530/REP-14-0140
Sood, P., Krek, A., Zavolan, M., Macino, G., and Rajewsky, N. (2006). Cell-type-specific signatures of microRNAs on target mRNA expression. Proc. Natl. Acad. Sci. U.S.A. 103, 2746–2751. doi: 10.1073/pnas.0511045103
Spornraft, M., Kirchner, B., Haase, B., Benes, V., Pfaffl, M. W., and Riedmaier, I. (2014). Optimization of extraction of circulating RNAs from plasma - Enabling small RNA sequencing. PLoS ONE 9:e107259. doi: 10.1371/journal.pone.0107259
Spornraft, M., Kirchner, B., Pfaffl, M. W., and Riedmaier, I. (2015). Comparison of the miRNome and piRNome of bovine blood and plasma by small RNA sequencing. Biotechnol. Lett. 37, 1165–1176. doi: 10.1007/s10529-015-1788-2
Stefani, G., and Slack, F. J. (2008). Small non-coding RNAs in animal development. Nat. Rev. Mol. Cell Biol. 9, 219–230. doi: 10.1038/nrm2347
Stocco, C., Telleria, C., and Gibori, G. (2007). The molecular control of corpus luteum formation, function, and regression. Endocr. Rev. 28, 117–149. doi: 10.1210/er.2006-0022
Sugino, N., and Okuda, K. (2007). Species-related differences in the mechanism of apoptosis during structural luteolysis. J. Reprod. Dev. 53, 977–986. doi: 10.1262/jrd.19047
Takagi, S., Simizu, S., and Osada, H. (2009). Reck negatively regulates matrix metalloproteinase-9 transcription. Cancer Res. 69, 1502–1508. doi: 10.1158/0008-5472.CAN-08-2635
Tesone, M., Stouffer, R. L., Borman, S. M., Hennebold, J. D., and Molskness, T. A. (2005). Vascular endothelial growth factor (VEGF) production by the monkey corpus luteum during the menstrual cycle: isoform-selective messenger RNA expression in vivo and hypoxia-regulated protein secretion in vitro. Biol. Reprod. 73, 927–934. doi: 10.1095/biolreprod.105.039875
Tripurani, S. K., Xiao, C., Salem, M., and Yao, J. (2010). Cloning and analysis of fetal ovary microRNAs in cattle. Anim. Reprod. Sci. 120, 16–22. doi: 10.1016/j.anireprosci.2010.03.001
Vlachos, I. S., Zagganas, K., Paraskevopoulou, M. D., Georgakilas, G., Karagkouni, D., Vergoulis, T., et al. (2015). DIANA-miRPath v3.0: deciphering microRNA function with experimental support. Nucleic Acids Res. 43, W460–W466. doi: 10.1093/nar/gkv403
Wang, X. G., Huang, J. M., Feng, M. Y., Ju, Z. H., Wang, C. F., Yang, G. W., et al. (2014). Regulatory mutation in the A2M gene are involved in the mastitis susceptibility in diary cow. Anim. Genet. 45, 28–37. doi: 10.1111/age.12099
Wang, X. P., Luoren, Z. M., Zan, L. S., Li, F., and Li, N. (2017). Bovine miR-146a regulates inflammatory cytokines of bovine mammary epithelial cells via targeting the TRAF6 gene. J. Dairy Sci. 100, 7648–7658. doi: 10.3168/jds.2017-12630
Wienholds, E., and Plasterk, R. H. (2005). MicroRNA function in animal development. FEBS Lett. 579, 5911–5922. doi: 10.1016/j.febslet.2005.07.070
Wiltbank, M. C. (1994). Cell types and hormonal mechanisms associated with mid-cycle corpus luteum function. J. Anim. Sci. 72, 1873–1883. doi: 10.2527/1994.7271873x
Yao, N., Yang, B.-Q., Liu, Y., Tan, X.-Y., Lu, C.-L., Yuan, X.-H., et al. (2010). Follicle-stimulating hormone regulation of microRNA expression on progesterone production in cultured rat granulosa cells. Endocrine 38, 158–166. doi: 10.1007/s12020-010-9345-1
Yoshioka, S., Abe, H., Sakumoto, R., and Okuda, K. (2013). Proliferation of luteal steroidogenic cells in cattle. PLoS ONE 8:e84186. doi: 10.1371/journal.pone.0084186
Keywords: microRNA, small RNA Seq, bovine, corpus luteum, estrous cycle, SOTA
Citation: Gecaj RM, Schanzenbach CI, Kirchner B, Pfaffl MW, Riedmaier I, Tweedie-Cullen RY and Berisha B (2017) The Dynamics of microRNA Transcriptome in Bovine Corpus Luteum during Its Formation, Function, and Regression. Front. Genet. 8:213. doi: 10.3389/fgene.2017.00213
Received: 27 September 2017; Accepted: 30 November 2017;
Published: 15 December 2017.
Edited by:
Emanuele Buratti, International Centre for Genetic Engineering and Biotechnology, ItalyReviewed by:
Maurizio Romano, University of Trieste, ItalyMauricio Fernando Budini, Universidad de Chile, Chile
Copyright © 2017 Gecaj, Schanzenbach, Kirchner, Pfaffl, Riedmaier, Tweedie-Cullen and Berisha. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Bajram Berisha, bajram.berisha@uni-pr.edu