Regulating life or death: Potential role of microRNA in rescue of the corpus luteum

Molecular and Cellular Endocrinology 398 (2014) 78–88

Contents lists available at ScienceDirect

Molecular and Cellular Endocrinology j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / m c e

Regulating life or death: Potential role of microRNA in rescue of the corpus luteum Samar W. Maalouf a, Wan-Sheng Liu a, Istvan Albert b, Joy L. Pate a,* a b

Department of Animal Sciences, Center for Reproductive Biology and Health, The Pennsylvania State University, University Park, PA 16802, United States Department of Biochemistry and Molecular Biology, The Pennsylvania State University, University Park, PA 16802, United States

A R T I C L E

I N F O

Article history: Received 31 May 2014 Received in revised form 7 October 2014 Accepted 8 October 2014 Available online 14 October 2014 Keywords: Corpus luteum MicroRNA Cattle High-throughput sequencing Maternal recognition of pregnancy

A B S T R A C T

The role of miRNA in tissue biology has added a new level of understanding of gene regulation and function. The corpus luteum (CL) is a transitory endocrine gland; the dynamic nature of the CL makes it a candidate for regulation by miRNA. Rescue of the CL from luteolysis is essential for the maintenance of pregnancy in all eutherian mammals. Using next generation sequencing, we profiled miRNA expression in the bovine CL during maternal recognition of pregnancy. We identified 590 luteal miRNA, of which 544 were known and 46 were novel miRNAs. Fifteen (including 3 novel) miRNAs were differentially expressed between CL of pregnant vs. cyclic animals. Target analysis of the differentially expressed miRNA resulted in genes involved in regulating apoptosis and immune response, providing evidence that miRNAs regulate the intracellular pathways that lead to either luteal regression or survival. © 2014 Elsevier Ireland Ltd. All rights reserved.

1. Introduction The corpus luteum (CL) is a transitory endocrine gland that forms from the cellular remnants of the ovulated ovarian follicle and secretes progesterone. At the end of a normal estrous or menstrual cycle in which fertilization has not occurred, the CL will rapidly regress, allowing for follicular maturation and initiation of a new cycle. Maternal recognition of pregnancy (MRP) is the response of the mother to conceptus signals that result in rescue of the CL. Luteal rescue, and thus continued progesterone production, is a fundamental requirement for the establishment of pregnancy. In humans, the conceptus signal responsible for luteal rescue is human chorionic gonadotropin (hCG). In ruminants, MRP depends on conceptus production of interferon-tau (IFNT), and occurs around days 15– 17 after estrus in cows (Bazer et al., 2009). IFNT acts in the uterus to alter prostaglandin production, facilitating rescue of the CL from luteolysis, but there are also changes in the CL that render it less sensitive to the luteolytic effects of prostaglandin F2α (PGF2A; Pratt et al., 1977; Silvia and Niswender, 1986). It is now recognized that IFNT is released from the uterus into the maternal circulation (Bott et al., 2010), and interferon-stimulated genes are increased in blood cells (Yankey et al., 2001) and in the CL (Oliveira et al., 2008; Yang et al., 2010). Interferon alpha (IFNA), which has biological properties similar to IFNT, suppresses cytokine-induced prostaglandin

* Corresponding author. Department of Animal Science, 324 Henning Building, University Park, PA 16802, United States. Tel.: +1 814 863 0558; fax: +1 814 863 6042. E-mail address: [email protected] (J.L. Pate). http://dx.doi.org/10.1016/j.mce.2014.10.005 0303-7207/© 2014 Elsevier Ireland Ltd. All rights reserved.

production (Benyo and Pate, 1992), expression of class II major histocompatibility molecules (Pate, 1995), and cytokine-induced apoptosis in luteal cells in vitro (Petroff et al., 2001). Interferonstimulated gene 15 (ISG15) is induced in ovine luteal cells by IFNT (Romero et al., 2013). Together, these findings suggest a possible direct effect of IFNT in the CL during MRP. Little is known regarding the molecular events that occur in the CL to facilitate luteal rescue. We hypothesized that these events involve regulation by microRNA, and therefore tested if changes in miRNA expression occur in the CL during MRP. The emergence of microRNA as epigenetic regulators of gene expression has changed the paradigm of tissue regulation and function (Ambros, 2008; Berezikov, 2011). Only 2% or less of the transcribed genomic DNA is translated into protein (Djebali et al., 2012; Maeda et al., 2006). Hence, the majority of RNA transcripts are noncoding RNA. MicroRNAs are small noncoding RNAs that regulate signaling pathways by targeting functional genes and modulating their expression. The mature (21–22 nt) miRNA binds to a seed sequence in the 3′ untranslated region (UTR), but may also bind in the 5′UTR of target mRNA (Lytle et al., 2007; Orom et al., 2008), to downregulate protein expression either through mRNA degradation or inhibition of translation. However, some evidence exists that miRNA may also enhance expression of some target genes (Orom et al., 2008; Place et al., 2008). MicroRNAs are either transcribed as separate genes or from introns of protein-coding mRNA (reviewed in Baley and Li, 2012). Regulation of gene expression by miRNA appears to be a conserved evolutionary process that provides very precise control of biological pathways, yet adds another level of complexity to understanding biological organisms.

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

The role of miRNA in reproduction has been documented by several studies. Dicer knockout mice exhibit an array of reproductive defects, including dysfunctional folliculogenesis, abnormal oocyte maturation and ovulation, luteal insufficiency and infertility (Hong et al., 2008; Lei et al., 2010; Nagaraja et al., 2008; Otsuka et al., 2008). Steroid hormones can both regulate and be regulated by miRNA (Fiedler et al., 2008; Macias et al., 2009; Sirotkin et al., 2009). Using either cloning based techniques or next generation sequencing, a number of studies identified miRNA in normal ovarian tissue from several species, including cattle (Hossain et al., 2009; Huang et al., 2011; Miles et al., 2012; Tripurani et al., 2010). Several studies have demonstrated a role of specific miRNA in cells from male and female reproductive tissues, as reviewed by Baley and Li (2012), Donadeu et al. (2012), Hossain et al. (2012) and Nothnick (2012). However, most of these used either whole ovarian tissue or luteal tissue collected only during the estrous cycle. Comprehensive miRNA profiling of the corpus luteum at the critical point of MRP is lacking. Despite the essential role of the CL in establishment of pregnancy in all mammals, the exact mechanisms regulating luteal survival during MRP are unknown. Therefore, we undertook this study to investigate the expression of miRNA and perform prediction analyses of the biological pathways involved in luteal rescue that are likely targets of regulation by miRNA. 2. Materials and methods 2.1. Animals, synchronization protocol, and artificial insemination Eight Holstein heifers (age 18–24 months) were watched daily for signs of estrus. At least 5 days after detection of estrus (day 0 = day of estrus), heifers were injected intramuscularly with 25 mg of prostaglandin F2α (Estrumate, Intervet, Summit, NJ). At the ensuing estrus, half of the heifers were artificially inseminated with two inseminations 6–8 hr apart. The other half was left unbred. On day 17 after breeding, the heifers were slaughtered at a USDA approved abattoir. Pregnancy was confirmed by identification of a conceptus in uterine flushes and the presence of a functional corpus luteum was confirmed by measuring plasma progesterone by enzyme linked immunosorbent assay (ELISA). Animal housing and treatments were in compliance with the Institutional Animal Care and Use Committee (IACUC) protocols at the Pennsylvania State University, University Park, PA. 2.2. Tissue collection and RNA extraction Corpora lutea of 8 heifers were collected on day 17 of the estrous cycle (cyclic, C; n = 4) or day 17 of pregnancy (pregnant, P; n = 4). The corpora lutea were immediately quartered and frozen in liquid nitrogen for RNA extraction. Total RNAs were extracted using the mirVana miRNA isolation kit (Life Technologies, Grand Island, NY) according to the manufacturer’s protocol. RNA quality was examined using ExperionTM automated electrophoresis system (BioRad, Hercules, CA). Samples with RNA quality index (RQI) of ≥9 (Best RQI = 10) were selected and shipped on dry ice to Beijing Genomic Institute (BGI; Hong Kong, China) for small RNA transcriptome sequencing analysis. 2.3. Small RNA sequencing and quality control RNA quality and integrity were re-examined after shipping using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). Samples with RIN > 8 were processed according to BGI pipeline protocols (Fig. 1). Single-end sequencing with a read length of 50 nt was performed on an Illumina HiSeq 2000 sequencer. Briefly, 18–30 nt length small RNA underwent 3′-adapter ligation followed by 5′-adapter ligation to obtain a 62–75 nt long sequence. A reverse transcription

79

primer, complementary to 3′-adapter, was used for reverse transcription. After PCR amplification, the purified PCR products were sequenced by HiSeq 2000 sequencer (BGI protocol, Hafner et al., 2008). Using Illumina base-calling algorithm, each base was assigned a quality score ranging from 0 (low) to 41 (high). Low quality reads (>4 bases with a quality score < 10) were eliminated from total reads. The resulting clean reads were analyzed for length distribution to determine the composition of small RNAs in each of the 8 samples. Counts of reads for each length category were normalized for each sample as follows: (number of reads for each length category/total reads) × 100. Uniformity of sequencing was assessed through the comparison of the percentages of common and sample-specific sequences in all 8 samples, two samples at a time (i.e. total of 28 pairwise comparisons). Samples with high count of common tags and low count of specific tags indicated a good uniformity of sequencing (BGI protocol). The Illumina reads have been submitted to NCBI’s Sequence Read Archive (SRA) database with the study accession number SRP048796.

2.4. Small RNA annotation The clean reads were aligned to the bovine genome (UMD3.1) using short oligonucleotide alignment program (SOAP) to analyze their expression and distribution. Mapped reads were further aligned against Rfam (Gardner et al., 2009; Griffiths-Jones et al., 2003, 2005) to classify and identify small RNA according to the different families of noncoding RNA (ncRNA). Annotated miRNAs were aligned to bovine miRNA precursors and mature miRNA in miRBase (http://microrna.sanger.ac.uk/ release 19, August 2012; Griffiths-Jones, 2004; Griffiths-Jones et al., 2006, 2008; Kozomara and Griffiths-Jones, 2011, 2014). Mature miRNAs were counted to determine their expression in cyclic vs. pregnant groups.

2.5. Quantitative polymerase chain reaction (qPCR) TaqMan small RNA assays (Life Technologies) were used for quantitative PCR. Total RNA samples (10 ng) were reverse transcribed using gene specific primers and TaqMan® MicroRNA Reverse Transcription kit. Specific primer pairs and probes to each miRNA and TaqMan Universal PCR Master Mix II were used for amplification as recommended by the manufacturer. MicroRNA-98, miR-374b, and miR-24 were used as reference miRNA representing high, medium and low reads respectively. The expression of each did not change significantly between samples from cyclic and pregnant animals, as tested by small RNAseq and qPCR analyses. The geometric mean of the amplification of reference miRNA was used for normalization. Relative expression of miRNA was expressed as 2−ΔCT.

2.6. MicroRNA target prediction, functional annotation and pathway analysis of the predicted target genes Target prediction was performed using TargetScan (release 6.2: June 2012) for bovine species (http://www.targetscan.org), miRanda http://www.microrna.org/microrna/home.do, microT-CDS http:// www.microrna.gr/microT-CDS, TarBase (Vergoulis et al., 2012), and QIAGEN’s ingenuity Target Explorer (QIAGEN’s, https:// targetexplorer.ingenuity.com/) for human homologs. Functional analysis of targeted genes was performed using the data annotation, visualization and integrated discovery (DAVID) software (http:// david.abcc.ncifcrf.gov/home.jsp), PANTHER classification system 9.0 (http://www.pantherdb.org) selected for Bos taurus, and QIAGEN’s Ingenuity® Pathway Analysis (IPA®, QIAGEN Redwood City, http:// www.qiagen.com/ingenuity).

80

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

Total RNA isolated from 8 bovine corpora lutea on day 17 of the estrous cycle (N=4) or pregnancy (N=4)

Single-end small RNA sequencing (50 nt tags, ~14.8 million reads/sample)

Filter low quality tags Full length small RNA tags

Unannotated sRNA

Potential miRNA with seed edit

Prediction of novel miRNA

Identification of Exon/Intron and of different RNA Species other than miRNA (such as rRNA, tRNA, snRNA, snoRNA, piRNA)

Differential expression

Bovine miRNA Annotated to miRBase v. 19

Identification of Known miRNA

Target prediction

GO enrichment, KEGG pathway, and ingenuity pathway analysis Fig. 1. Workflow of small RNA Seq in bovine corpora lutea (CL). Total RNA were isolated from the CL of cyclic (n = 4) and pregnant (n = 4) heifers on day 17 of the estrous cycle. Small RNA libraries were prepared for single end sequencing. miRBase release 19 was used for miRNA annotation. Results of small RNASeq were analyzed for differential expression by DESeq. TargetScan was used to identify putative genes targeted by the differentially expressed known miRNA. Gene ontology and pathway analyses of human homolog of targeted genes were generated by DAVID, PANTHER, and IPA.

2.7. Novel miRNA prediction Novel miRNAs were predicted using Mireap (http://sourceforge.net/projects/mireap/) prediction software (BGI software). Predictions were based on the secondary structure, dicer cleavage site and the minimum free energy of the unannotated small RNA sequences that mapped to the bovine reference genome (UMD3.1). Predicted novel miRNA had to fulfill the following conditions: 1. Align to intron or antisense exon regions. 2. Form a hairpin secondary structure with a mature miRNA present in one arm of the hairpin precursor. 3. The mature miRNA and its complementary strand (miRNA*) have 2-nucleotide 3′-overhangs. 4. Hairpin precursors lack large internal loops or bulges. 5. The secondary structure of the hairpin is steady with free energy of hybridization lower than or equal to −18 kcal/mol. 6. The number of mature miRNA with predicted hairpin must be no fewer than 5 in the alignment result. The

expression of novel miRNA was obtained by adding the counts of miRNA with no more than 3 mismatches on the 3′ and 5′ ends and no mismatch in the middle from the alignment results. The novel miRNA sequences were submitted to miRBase and assigned official names. 2.8. Statistical analysis High-throughput sequencing assays such as miRNA-Seq provide quantitative readouts in the form of count data. To infer differential expression in such data an estimation of data variability throughout the dynamic range and a suitable error model are required. We applied a widely used methodology called DESeq (Anders and Huber, 2010) that is based on the negative binomial distribution, with variance and mean linked by local regression. Adjusted p values (P adj) were computed to correct for multiple

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

comparisons. To compare length distribution, number of total sRNA and sRNA mapped to the genome between the two treatment groups, Student’s t-test was performed to determine a significant difference.

81

A 100

Frequency (%)

3. Results and discussion 3.1. Small RNA sequencing results

Cyclic Pregnant

*

80 60 40

* *

20

The aim of this study was to profile miRNA expression in the bovine CL on day 17 of pregnancy compared to that in nonpregnant animals on the same day. Plasma progesterone was not different in cyclic (6.12 ± 1.77 ng/ml) vs. pregnant (6.66 ± 2.42 ng/ml) heifers, confirming that CL were not undergoing luteal regression at the time of collection. To determine the expression of miRNA in the CL during maternal recognition of pregnancy, we utilized a genomewide approach. Total RNA from corpora lutea of cyclic and pregnant heifers was analyzed by single-end small RNA sequencing. Fig. 1 summarizes the steps to select high quality RNA reads and subsequent bioinformatic analyses. The sequencing of each sample resulted, on average, in ~14.8 million reads, with 99.5% classified as high quality reads (Table 1). Contaminants such as number of reads with no 3′-adapter, no insert, 5′-adapter contaminants, sequences shorter than 18 nt, or reads with a poly A tail were eliminated and the resulting clean reads constituted an average of 99.22% of high quality reads (Table 1). Length distribution analysis of clean reads resulted in 77.9 ± 4.0% of small RNA falling in the 21–22 nt length range, typical for miRNA (Lee et al., 2002), and ~22% were in the 23–30 nt range (Fig. 2A). By comparing length distribution in samples from cyclic and pregnant animals, we observed a small but significant (p < 0.01) increase in the percent of 22 nt long reads and a small but significant (p < 0.01) decrease in the percent of 23 nt and 24–30 nt long reads in the pregnant group (Fig. 2A). This was an interesting observation, but the biological significance of such a difference is yet to be determined. To assess the uniformity of sequencing among the 8 different samples, the common and sample specific sequences were compared pairwise for all 8 samples. Fig. 2B depicts a Venn diagram comparing representative samples, C1 and P4. The average of all 28 pairwise comparisons resulted in 98.8% of common small RNA and 1.2% specific tags indicating a uniformity of sequencing among the different samples (BGI protocol).

0 21

B

22 23 Length (nt)

24-30

Total small RNA 1.71% 0.27% 98.02%

C1 specific Common to C1 & P4 P4 specific Fig. 2. Quality control of HiSeq raw data. (A) Length distribution of small RNA in cyclic (open bars) and pregnant (solid bars) groups. Each value represents the mean ± SD (n = 4). * Denotes significant difference between groups (p < 0.01). (B) Venn diagram comparing two representative samples (C1 and P4) for common and specific sequences of total small RNA. A high percent of common total reads vs. a low percent of specific reads indicates uniformity of sequencing between these two samples.

3.2. Identification of known miRNA in the bovine corpus luteum The percent of full-length small RNA tags (clean reads) that mapped to the bovine genome (UMD3.1) was 87% on average, which correlated with the completeness of the bovine genome sequence assembly in version UMD3.1. By comparing the percent of mapped reads in samples from cyclic (84.79 ± 0.60) to those from pregnant (88.96 ± 0.90) animals, a small but significant (p < 0.01) increase in

Table 1 Small RNA sequencing quality control. Type

C1

C2

C3

C4

P1

P2

P3

P4

Total readsa High qualityb (% of total reads) 3′ adapter nullc Insert nulld 5’ adapter contaminantse Shorter than 18 ntf Poly Ag Clean readsh (% of high quality reads)

16.52e+6 16.47e+6 (99.71) 10053 3179 102926 31859 170 16.32e+6 (99.10)

14.17e+6 14.13e+6 (99.70) 8253 2085 99696 19020 41 14.00e+6 (99.09)

14.43e+6 14.38e+6 (99.71) 10065 3262 91309 13175 37 14.27e+6 (99.18)

12.83e+6 12.78e+6 (99.62) 15992 1727 66170 11272 10 12.69e+6 (99.26)

14.84e+6 14.79e+6 (99.69) 15525 2104 61673 4097 15 14.71e+6 (99.44)

16.88e+6 16.83e+6 (99.69) 18752 1481 57127 2785 21 16.75e+6 (99.52)

13.20e+6 13.16e+6 (99.66) 14852 4271 81049 10890 45 13.05e+6 (99.16)

15.89e+6 15.84e+6 (99.66) 21165 1905 127276 3517 137 15.68e+6 (99.03)

a

Total reads: total sequenced reads. High quality: number of high quality reads with no N, no more than 4 bases whose quality score is less than 10 and no more than 6 bases whose quality score is less than 13. c 3′ adapter null: number of reads with no 3′ adapter. d Insert null: number of reads with no insert. e 5′ adapter contaminants: number of 5′ contaminants. f Shorter than 18 nt: number of tags shorter than 18 nt. Acceptable tags are between 18 and 30 nt. g Poly A: number of reads with a poly A tail. h Clean reads: number of clean reads after adapter and contaminants are removed. b

82

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

Table 2 Total small RNA mapped to the genome. Cyclic

Total sRNA Total sRNA mapped to the genome % sRNA mapped to the genome miRNA Unannotated miRNA Other small RNA a

p-valuea

Pregnant

Mean

SD

Mean

SD

1.43e+07 1.21e+07 84.79 82.45 15.36 2.19

1.50e+06 1.21e+06 0.60 1.62 0.63 1.11

1.50e+07 1.34e+07 88.96 87.66 11.13 1.21

1.57e+06 1.53e+06 0.90 1.2 0.95 0.32

0.53 0.24 2.54e-4 2.09e-3 3.02e-4 0.14

p-value depicts the difference between cyclic and pregnant groups (n = 4).

the percent of small RNA in samples from pregnant animals was observed (Table 2). Furthermore, reads that mapped to the genome were also mapped to Rfam in order to identify the different types of ncRNA. As a result, 85% of the annotated small RNA were identified as miRNA, 13% were unannotated sequences, and the remaining 2% mapped to introns, exons, repeats or other types of small RNA (Table 2). These results are in agreement with what is known of miRNA, which constitutes 52–87% of small RNA in the ovary (Ahn et al., 2010; Hossain et al., 2009; Huang et al., 2011; Mishima et al., 2008; Ro et al., 2007). Consistent with the length distribution and unannotated results, there was an increase in the percentage of annotated miRNA in the CL of pregnant animals compared to that of cyclic animals (Table 2). MicroRNAs were aligned against miRBase (version 19) and the total number of annotated miRNA expressed in the CL was 544 (Supplementary Appendix S1). Of these 544 miRNAs, 40% had only 1–4 reads, 50% had reads between 5 and 4500, and 10% had reads above 4500. The percentages were not different when we compared results from cyclic and pregnant groups. MicroRNAs that share the same seed sequence and function are grouped into families (Ding et al., 2011; Kozomara and Griffiths-Jones, 2011; Zou et al., 2014). Analysis of the 544 miRNAs resulted in 237 families of bovine miRNA, of which 13 miRNA families (miR-2284, miR-2300, miR-2319, miR2329, miR-2363, miR-2404, miR-2450, miR-2887, miR-2957, miR3432, miR-3604, miR-6526, miR-6536) were bovid specific. Twenty one miRNA families were not found in human, mouse or rat species but could be present in lower vertebrates or invertebrates. A complete list of all 237 families and the species in which they have been documented can be found in Supplementary Appendix S2. It is interesting to note that the number of total miRNA we identified in the CL is very close to the number of total unique miRNA identified in the bovine ovary (472) by Hossain et al. (2009) or in adult porcine gonads (673) by Li et al. (2011). MicroRNA-21, Let-7, miR-

320, miR-140, miR-199a-3p were among the most abundant miRNA in the bovine CL, consistent with the pig ovaries and testis (Li et al., 2011). Hossain et al. (2009) found let-7, miR-21 as well as miR23b, miR-24, miR-27a, miR-126, and miR-143 to be the most abundant miRNA in bovine ovaries. The abundance of these miRNAs may indicate their involvement in housekeeping regulation of the gonads. 3.3. Differential expression of known miRNA in the bovine corpus luteum Differential expression analysis using DEseq resulted in 9 miRNAs that were significantly different (P adj < 0.05) in the CL of cyclic vs. pregnant animals (Table 3). Of the 9 differentially expressed miRNAs, 8 (miR-2443, miR-345-3p, miR-1307, miR-2389, miR-222, miR-3655p, miR-6517 and miR-541) were downregulated and 1 (miR302b) was upregulated in the pregnant group (Table 3). In addition, 3 miRNAs (miR-1-1, miR-10b, and miR-122) showed a strong tendency (p < 0.07) to be upregulated in the pregnant group (Table 3). Quantitative PCR was performed on selected miRNA, representing low, medium and high reads (Fig. 3). MicroRNA-302b, which had the lowest number of reads (<6), was amplified by qPCR in samples from pregnant, but not cyclic, animals. This pattern of expression confirmed the results obtained by next generation sequencing. MicroRNA-1-1, which had reads above 3000, was detected by qPCR and its pattern of expression followed the same pattern observed by small RNAseq analysis, albeit with somewhat lower fold change between groups. Analysis of two miRNA with medium number of reads, miR-222 and miR-1307, by qPCR did not result in confirmation of the treatment effect observed in the NGS results. Both of these miRNAs were downregulated in early pregnancy based on NGS data, but qPCR did not reveal a difference in these miRNAs due to treatment. Although many researchers choose to ‘confirm’ the results

Table 3 Differentially expressed miRNA in the CL on day 17 of the estrous cycle or pregnancy. ID

Base meanaA (cyclic; C)

Base meanaB pregnant; P)

Fold changeb (P/C)

P valc

P adjd

miR-1-1 miR-10b miR-122 miR-222 miR-302b miR-345-3p miR-365-5p miR-541 miR-1307 miR-2389 miR-2443 miR-6517

3992.96 3959.297 164.95 3364.70 0.49 509.24 1238.03 500.29 3343.79 13.14 265.89 74.99

8511.07 6403.14 265.88 1846.40 6.44 117.69 689.97 230.26 844.05 1.54 61.02 37.14

2.13 1.62 1.61 0.55 13.28 0.23 0.56 0.46 0.25 0.11 0.23 0.50

0.00105 0.00111 0.00151 3.16E-05 0.000664 3.94E-09 0.000334 0.000695 2.41E-06 1.13E-05 1.76E-11 0.000395

0.0547 0.0547 0.0686 0.00344 0.0420 1.07E-06 0.0264 0.0420 0.000438 0.00154 9.57E-09 0.0264

a b c d

Base mean is the number of reads divided by the size factor (normalization constant) of the sample. Fold change is the ratio of expression in pregnant/cyclic. P val = p-value. P adj = p-value adjusted for multiple comparisons.

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

NGS

0.05

0.05

NGS (Reads x 102)

0.10

80

60

60

40

40

20

20

40

30

30

20

20

10

10

0

Cyclic

Pregnant

0

NGS (Reads x 102 )

40

Cyclic

Pregnant

miR-222

0

40

30

30

20

20

10

10

0

Cyclic

Pregnant

qPCR (2-ΔCT)

40 miR-1307

Pregnant

qPCR (2-ΔCT x 10-2)

0

Cyclic

100

80

0.00

0.00

miR-1-1

qPCR (2-ΔCT x 10-2)

0.10

NGS (Reads x 102 )

100

0.15

miR-302b

qPCR (2-ΔCT)

NGS (Reads x 102)

0.15

qPCR

0

Fig. 3. qPCR analysis of selected differentially expressed miRNA. Bars depict the base mean values generated for next generation sequencing (NGS) data by DESeq (solid bars) or expression relative to reference miRNA by qPCR analysis (2−ΔCt; open bars) comparing miRNA in CL from pregnant vs. cyclic animals (n = 4 per group).

83

Ma et al. (2011) have shown, both by microarray and quantitative PCR analyses, that 7 miRNAs (miR-378, miR-22-5p, miR-147, miR-877, miR-24, miR-107, and miR-148b) were upregulated in nonregressed compared to regressed bovine CL. Although all 7 were expressed in our dataset, none of these miRNAs were differentially expressed between the cyclic and pregnant groups. The stage of the CL used by Ma et al. (2011) was determined based on blood serum progesterone concentration and macroscopic observation of the ovary, which precludes exact definition of luteal stage, and did not include a comparison of CL from pregnancy, as did the current study. Further, the current study did not include regressing CL. Therefore, it is not surprising that miRNA expression detected in this study was different from that of Ma et al. The most abundant miRNA identified in bovine ovaries or follicular and luteal tissues (Let-7, miR143-3p, miR-23b, miR-24, miR-27a, miR-652, miR-21, miR-126; reviewed by Donadeu et al., 2012) were also abundant in our dataset, with the exception of miR-652 (<70 reads), but were not differentially expressed in the CL on day 17 of the estrous cycle or pregnancy, suggesting that these miRNAs are involved in housekeeping functions rather than life or death decisions of the CL. Thirty three percent of pig gonadal miRNAs were differentially expressed in ovaries compared to testes (Li et al., 2011), including miR-222, which was also differentially expressed in our NGS data. However, according to Hossain et al. (2009), miR-222 was present in cumulus cells but not the ovarian cortex or corpus luteum. In summary, there is general agreement in miRNAs that were expressed in the CL (this study) with previous reports, but those miRNAs that were differentially expressed in CL of cyclic vs. pregnant animals represent a novel set of miRNA that may be regulators of luteal rescue. 3.4. Target prediction

of NGS using qPCR, as is very appropriately done following microarrays, it is erroneous to assume that the precision of qPCR (typically used to detect relative differences) is equal to, or superior to, NGS (a method for absolute quantification). First, the normalization methods for NGS and qPCR are completely different. In qPCR, the gene of interest is normalized to a reference gene that is not different between treatments, but still may exhibit some variability among individual samples. In DEseq of the NGS data, the assumption is that no miRNAs are changed by treatment, and are normalized to the geometric mean of all miRNA reads (Dillies et al., 2012). The difference in normalization methods is one reason that NGS and qPCR do not always produce the same result. Secondly, the high fidelity polymerase enzyme used in small RNA qPCR assays may fail to amplify transcripts if a single mismatch occurs at the 3′ end. Further, primers used in qPCR may not hybridize with all isomiRs that are detected by NGS (Morin et al., 2008). Therefore, comparison of NGS and qPCR data must be done in light of these considerations. MicroRNA-302b exhibited the greatest fold change in response to pregnancy (13-fold), but had the lowest number of reads. Although it is tempting to focus on this miRNA that demonstrated the greatest fold change relative to the experimental group, the biological significance of a very low number of miRNA copies is questionable (Ambros, 2008; Bartel, 2004). Therefore, it remains to be determined if the increase in miR-302b is biologically significant, and if miR-302b counts would further increase after day 17 of pregnancy. Alternatively, miR-1-1 increased by just 2 fold in the pregnant group, but the number of miR-1-1 reads was greater than any other differentially expressed miRNA (Table 3), suggesting that miR-1-1 could have considerable impact on translational regulation in the CL. The other 8 miRNA decreased by 1.8–9 fold, with miR2389 decreasing to the greatest extent. However, miR-2389 also had the lowest number of reads in that group (Table 3).

Appendix S7: Supplementary Table S1 summarizes the number of targets identified for each miRNA family in each of the different algorithms available. The identification of valid miRNA targets is a major challenge. Several algorithms, as well as several methods for identifying miRNA targets, have been utilized (Berezikov, 2011; Ruby et al., 2007). TargetScan and miRanda predict miRNA targets based on seed complementarity. MicroRNAs that share the same seed sequence are grouped in one family in TargetScan as depicted in Appendix S7: Supplementary Table S1. This family nomenclature is different from the family referred to in miRBase. Other algorithms such as microT-CDS rely on the thermodynamic stability of miRNA–mRNA base pairing. TarBase and IPA predictions are based on experimentally validated targets. Since the precision of target prediction by any available algorithm is low (Zheng et al., 2013), targets predicted by more than 2 algorithms would be more likely bonafide targets. However, not all the differentially expressed miRNAs were recognized in each algorithm. MicroRNA-6517, which was identified only for bovine species in miRBase, had no predicted targets in any of the algorithms utilized. MicroRNA-10b, miR-122, miR365-5p, miR-2443, and miR-2389 were not recognized by the Diana tool database (microT-CDS and TarBase). No targets were predicted for miR-1-1 and miR-365-5p in miRanda. Furthermore, miR2443, miR-345-3p, miR-365-5p, and miR-6517 were not found in the IPA database, and therefore had no predicted targets. TargetScan was the only algorithm that predicted targets for all differentially expressed miRNAs, except for miR-6517. Therefore, the subsequent gene ontology and pathway analysis in this study were conducted on targets predicted by TargetScan. MicroRNA-6517 was omitted from subsequent target analyses. The genes targeted by miRNA are cell context dependent. One miRNA may target several genes, and one gene may have binding sequences for several miRNAs. The multiplicity of miRNA binding sites in the UTR of target genes may play a synergistic or additive function in the regulation of such targets (Bartel, 2009; Berezikov,

84

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

2011; Brennecke et al., 2005; Farh et al., 2005; Grimson et al., 2007). Hence, the role of miRNA in gene expression would most likely be that of a fine-tuning process rather than an ON/OFF switch. One gene may be targeted by up and downregulated miRNA at the same time in order to attain the optimum concentration required for a specific function (Herranz and Cohen, 2010). To gain insight into the potential function of differentially expressed miRNA during MRP, we analyzed the combined 2431 predicted targets (human homologs, Supplementary Appendix S3) of the 9 differentially expressed miRNAs, as well as the additional 3 miRNAs that showed a strong tendency for upregulation on day 17 of pregnancy. Of these targets, 3 genes (neuronal PAS domain protein 3 (NPAS3), AP2 associated kinase 1 (AAK1), and mitogen-activated protein kinase kinase kinase 2 (MAP3K2) were common targets of 5 differentially expressed miRNAs, 18 were common targets of 4 miRNAs, 68 were common to 3 miRNAs, and 410 genes were commonly targeted by 2 miRNAs. The remaining 1932 genes were specific targets for one of the differentially expressed miRNAs. Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide (YWHAZ), calcium/calmodulindependent kinase kinase 2 beta (CAMKK2), guanine nucleotidebinding protein (G protein) beta polypeptide 1 (GNB1) and regulatory G protein signaling 2 (RGS2) have been suggested to play a role in the acquisition of luteal sensitivity to PGF2A (Goravanahally et al., 2009; Juengel et al., 1998). In the CL that is insensitive to PGF2A, YWHAZ expression is high and CAMKK2 expression is low. In the CL that is sensitive to PGF2A, the reverse expression pattern is observed. Moreover, YWHAZ codes for a protein kinase C inhibitor that inhibits CAMKK2 activity (Davare et al., 2004), leading to increased resistance of the CL to luteolysis. Our target analysis showed YWHAZ and CAMKK2 to be predicted targets of miR-1 and miR345-3p, respectively. Earlier studies have shown the CL of pregnancy to be less sensitive to PGF2A than the CL of the estrous cycle (Pratt et al., 1977; Silvia and Niswender, 1986). Based on their known roles in acquisition of luteal sensitivity to PGF2A, we would expect YWHAZ to be upregulated and CAMKK2 to be downregulated in pregnancy. However, any conclusion about their expression pattern in the rescued CL requires further analysis by RNASeq and/or Western blot analyses. In a study by Romero et al. (2013), serpine 1 and thrombospondin 1 (THBS1) were upregulated in the ovine CL during luteolysis, and PGF2A stimulated THBS1 in bovine CL (Zalman et al., 2012) and luteinized granulosa cells (Farberov and Meidan, 2014). THBS1 promotes apoptosis in endothelial cells by activation of the CD36/ p59Fyn/caspase-3/p38 MAPK cascade and upregulation of mediators of apoptosis such as cytokines (TNF, FasL), and Bax (Jimenez et al., 2000; Nor et al., 2000; Rege et al., 2009). In bovine luteal endothelial cells, THSB1 inhibited cell migration and induced apoptosis (Zalman et al., 2012), and knockdown of THBS1 in luteal endothelial cells or luteinized granulosa cells reduced concentrations of cleaved caspase 3 (Farberov and Meidan, 2014). Moreover, THSB1 is a predicted target for miR-1-1 and miR-222. In our study, miR222 was downregulated while miR-1-1 was upregulated on day 17 of pregnancy, but the ratio of normalized reads of miR-1-1 to miR222 was about 3.5 fold greater in the CL of pregnancy than in the CL of the cycle, suggesting that THSB1 may be a target of miRNA involved in rescue of the CL. Pentraxin 3 (PTX3), interleukin 6 (IL6), vascular endothelial growth factor (VEGF) and luteinizing hormone receptor (LHR) mRNA were greater in the CL of pregnant compared to CL of nonpregnant ewes on Day 15 (Romero et al., 2013). The current list of predicted targets included neuronal pentraxin 1 (NPTX1), NPTX receptor (NPTXR), VEGFA, IL6 signal transducer (IL6ST), and IL6 receptor (IL6R). VEGFA is predicted to be a target of miR-1-1 and miR302b, while NPTX1, NPTXR, IL6R and IL6ST are predicted to be targeted by miR-365-5p, miR-541, miR-222 and miR-302b, respec-

GO: Biological Process 5 14 64 70 15 6

18

18 43

47

7

22 cellular component organization or biogenesis (GO:0071840) cellular process (GO:0009987) localization (GO:0051179) apoptotic process (GO:0006915) reproduction (GO:0000003) biological regulation (GO:0065007) response to stimulus (GO:0050896) developmental process (GO:0032502) multicellular organismal process (GO:0032501) biological adhesion (GO:0022610) metabolic process (GO:0008152) immune system process (GO:0002376) Fig. 4. Gene ontology of biological pathways. The 142 genes enriched for regulation of apoptosis and cell death were classified by PANTHER according to their gene ontology. The numbers on the pie chart represent the number of genes sharing the same gene ontology nomenclature.

tively. Although further studies are required to determine expression of these genes in the bovine CL during maternal recognition of pregnancy, their differential expression in the ovine CL supports our target prediction of these genes as regulated by miRNA in the CL. 3.5. Functional analysis of target genes Using the functional annotation cluster (FAC) tool in DAVID (Huang et al., 2009a, 2009b) at the highest stringency, the 2431 genes were clustered into 82 enriched GO terms, of which 9 had an enrichment score ≥ 1.3 (equivalent to a p-value of 0.05; Appendix S7: Supplementary Table S2). Of the enriched clusters, the cluster with the highest enrichment score (ES = 4.79) and the lowest false discovery rate (FDR < 0.01) included 142 of the 2431 predicted genes that were enriched for apoptosis and cell death (Cluster #1; Supplementary Appendix S4). Other clusters with significant enrichment scores (>1.3) were involved in the regulation of protein kinase activity (ES = 3.95), nuclear import (ES = 2.59), protein-lysine-Nmethyltransferase and histone methyltranserase activities (ES = 1.91), response to glucose (ES =1.62), and other metabolic pathways (Appendix S7: Supplementary Table S2). All these pathways fall nicely into our prediction that the differentially expressed miRNA would target cellular processes important for cell death and survival. Moreover, we explored the gene ontology biological pathways in which the 142 genes enriched in cluster 1 participate, using Panther classification system (Huaiyu et al., 2012, 2013). Fig. 4 shows that the 142 apoptosis enriched genes were further clustered into several GO terms, of which 18 genes were involved in the apoptotic process (GO: 0006915), 7 genes were involved in reproduction (GO:0000003), and 14 genes were involved in the immune system process (GO: 0002376). Fig. 5 shows the child GO terms of apoptotic process (Fig. 5A), immune system process (Fig. 5B) and immune response (GO: 0042116; Fig. 5C).

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

A.

Apoptotic process (GO:0006915)

B.

Immune System Process (GO:0002376)

2

4

6

4

Induction of apoptosis (GO:0006917)

Macrophage activation (GO:0042116)

Immune response (GO:0006955) Negative regulation of apoptotic process (GO:0043066) C. Immune Response (GO:0006955)

1

2

1

response to interferon-gamma (GO:0034341) natural killer cell activation (GO:0030101)

85

Further analysis of the 89 common targets of 3 or more miRNA by IPA revealed that the major pathways associated with these predicted target genes are immune response pathways (p < 0.001; Fig. 6). Luteal function and regression involve immune cell-mediated homeostatic and destabilizing mechanisms, and there is considerable evidence from our laboratory and others that cytokines and immune cells facilitate both the inhibition in progesterone production and initiation of apoptosis in luteal cells (rev. Walusimbi and Pate, 2013). The ability of PGF2A to effectively cause luteolysis is dependent on upregulation of immune response genes in the CL (Atli et al., 2012). Further, there is now evidence that the embryonic signal, IFNT, escapes the uterus and can activate interferon-stimulated genes in the CL (Bott et al., 2010). Thus, we have hypothesized that rescue of the CL during MRP requires regulation of the functions of the resident immune cells in the CL. Subtle changes in resident T cells were described during MRP (Poole and Pate, 2012). The results of this pathway analysis further support the concept that resident immune cells are key players in the fate decision that besets the CL around the time of MRP, and further implicates miRNAs as regulators of immunelike responses in the CL. Table 4 summarizes the genes in the six top signaling pathways (p < 0.0001; B-cell receptor signaling, PKCΘ signaling in T lymphocytes, RANK signaling in osteoclasts, NGF signaling, calcium signaling, and T cell receptor signaling). In addition to roles in immune response pathways, a number of these genes are important regulators of steroidogenesis and PGF2A action in the CL (Arvisais et al., 2010; Carambula et al., 2003; Wiltbank et al., 1991).

B cell mediated immunity (GO:0019724)

3.7. Novel miRNA identified in bovine corpora lutea Fig. 5. Child gene ontology. (A) Child GO of the parent apoptotic process, (B) child GO of the parent immune system process, and (C) child GO of the parent immune response. The numbers on the pie chart represent the number of genes sharing the same gene ontology.

3.6. Functional analysis of predicted genes targeted by 3 or more miRNAs Effects of miRNA on a given target may depend on the simultaneous expression of other miRNAs targeting the same gene (rev. in Berezikov, 2011). Hence, we further analyzed the predicted genes that were targets of at least 3 differentially expressed miRNAs (89 genes; Supplementary Appendix S3) using the DAVID clustering function. By applying the highest stringency condition as in Section 3.5, target genes were clustered into 15 functional clusters, of which 5 had an enrichment score of ≥1.3 (Appendix S7: Supplementary Table S3). Negative regulation of cellular catabolic process and positive regulation of transcription and gene expression were among the highly enriched GO terms.

Unannotated small RNAs were analyzed by Mireap prediction software to predict novel miRNA. A total of 226 sequences were identified by BGI. Forty-six of these sequences were assigned official nomenclature by miRBase. Of these novel miRNA, 25 were not found in any other species, the rest were variations of already identified miRNA or novel mature sequences of previously known stem loop sequences. Analysis for differential expression by DESeq resulted in 3 differentially expressed novel miRNAs (Table 5, and Supplementary Appendix S5). Bos taurus-miR-10178 was only present in the CL of cyclic animals, whereas bta-miR-10163 was only present in the CL of pregnant animals. Bos taurus-miR-10170 was downregulated in the CL of pregnant animals compared to the CL of cyclic animals by 6.7 fold. In order to gain insight into the biological importance of these novel miRNAs, we utilized Diana microT3.0 target prediction software (Maragkakis et al., 2009a, 2009b), to generate predicted targets based on the sequence of each miRNA. We set the threshold to the least stringent conditions (threshold score = 7.3 with maximum stringency score being 19) in order to include targets from all 5 novel miRNAs. A total of 484 genes

Table 4 List of genes associated with the top 6 canonical pathways by IPA analysis. Symbol

Entrez gene name

Location

Type(s)

Entrez gene ID for human

CALM1 (includes others) CAMKK1 CBL CREB1 CREB5 HDAC4 MAP3K2 MAP3K13 MAP3K14 NFAT5 NFATC3 PIK3R3 POU2F1

Calmodulin 1 (phosphorylase kinase, delta)

Cytoplasm

Other

801|805|808

Calcium/calmodulin-dependent protein kinase kinase 1, alpha Cbl proto-oncogene, E3 ubiquitin protein ligase cAMP responsive element binding protein 1 cAMP responsive element binding protein 5 Histone deacetylase 4 Mitogen-activated protein kinase kinase kinase 2 Mitogen-activated protein kinase kinase kinase 13 Mitogen-activated protein kinase kinase kinase 14 Nuclear factor of activated T-cells 5, tonicity-responsive Nuclear factor of activated T-cells, cytoplasmic, calcineurin-dependent 3 Phosphoinositide-3-kinase, regulatory subunit 3 (gamma) POU class 2 homeobox 1

Cytoplasm Nucleus Nucleus Nucleus Nucleus Cytoplasm Cytoplasm Cytoplasm Nucleus Nucleus Cytoplasm Nucleus

Kinase Transcription regulator Transcription regulator Transcription regulator Transcription regulator Kinase Kinase Kinase Transcription regulator Transcription regulator Kinase Transcription regulator

84254 867 1385 9586 9759 10746 9175 9020 10725 4775 8503 5451

86

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

Ratio= # of molecules in the data set #molecules in pathway -log(p-value)= strength of association of data set molecules with a pathway

B Cell Receptor Signaling PKCΘ Signaling in T Lymphocytes RANK Signaling in Osteoclasts NGF signaling Calcium signaling T Cell receptor Signaling Glucocorticoid Receptor Signaling Neurotrophin/TRK Signaling PI3K Signaling in B lymphocytes GNRH Signaling

FLT3 Signaling in Hematopoeitic Progenitor Cells Role of Macrophages, Fibroblasts and Endothelial cells in Rheumatoid Arthritis Role of Osteoblasts, Osteoclasts and Chondrocytes in Rheumatoid Arthritis Cardiac Hypertrophy Signaling FGF Signaling April Mediated Signaling Phospholipase C Signaling B Cell Activating Factor Signaling SAPK/JNK Signaling

Fig. 6. Functional analysis by IPA: target genes common to ≥3 miRNA were analyzed by IPA. The top 19 pathways are plotted. The lower the p-value the more significant is the correlation of that function with target genes.

were targeted by differentially expressed novel miRNA (Supplementary Appendix S6). One gene, pyridoxal kinase (PDXK), was targeted by miR-10163 and miR-10178. PDXK is a conserved gene in several species and is important in the synthesis of the active coenzyme pyridoxal-5′-phosphate (PLP) from vitamin B6 precursor (Yang et al., 1998). Other targets included cytochrome P450 XXI (CYP21A2) involved in steroid hormone biosynthesis (Ryan and Engel, 1957), retinoic acid receptor RXR-alpha (RXRA), a member of the steroid hormone receptor superfamily of transcription factors involved in cell proliferation and apoptosis (Liu et al., 2000), adenylate cyclase type 1 (ADCY1), calcium channel (CACNA1C) and sodium channel subunit beta-2 precursor (SCN2B), all known to be involved in steroidogenesis (Niswender et al., 2000). Target analysis of all 484 targets of the 3 differentially expressed novel miRNA by DAVID, under conditions of highest

stringency as in Section 3.5, resulted in 43 clusters, of which regulation of virus-host interaction was cluster #1 (ES = 1.53). The other clusters included processes such as negative regulation of transcription, regulation of apoptosis, metal ion binding, regulation of fatty acid metabolic process, and others. Therefore, these results show that novel miRNA also targets apoptosis and immune response signaling in the CL during maternal recognition of pregnancy. 4. Conclusion Rescue of the CL is required for the establishment of pregnancy in all eutherian mammals. The role of miRNA in the regulation of luteal function is not well described, and no studies have explored miRNA expression during the critical period of maternal recognition of pregnancy. In this study we profiled miRNA

Table 5 Novel microRNAs that are differentially expressed in the CL on day 17 of the estrous cycle or pregnancy. ID

Sequence

Base meanaA (cyclic; C)

Base meanaB (pregnant; P)

P valb

P adjc

bta-miR-10163 bta-miR-10178 bta-miR-10170

TGACAAGCATGGTTCTCTCTCT CCTCCTGGGACCTGAGGCAGCA ACCGGCGCAGCGGCTGACCTGG

0.00 8.47 16.19

16.45 0.00 2.43

5.62E-05 4.84E-04 4.48E-03

1.46E-03 8.40E-03 4.66E-02

a b c

Base mean is the number of reads divided by the size factor (normalization constant) of the sample. P val = p-value P adj = p-value adjusted for multiple comparisons.

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

expression in corpora lutea on day 17 of the estrous cycle or pregnancy, to identify miRNA with a putative role in the rescue of the CL during maternal recognition of pregnancy. Using next generation sequencing, we identified 544 total miRNAs expressed in the CL, of which 9 were differentially expressed and 3 miRNAs expressed a hightendency for differential expression in CL of cyclic vs. pregnant animals. Eight miRNAs (miR-2443, miR-345-3p, miR1307, miR-2389, miR-222, mmiR-365-5p, miR-6517, and miR541) were downregulated, 1 miRNA (miR-302b) was significantly (p = 0.042) upregulated, and 3 miRNAs (miR-1-1 and miR-10b (p = 0.054) and miR-122 (p = 0.069)) had a strong tendency to be upregulated in the CL of pregnancy. Target analysis of these miRNAs resulted in 2431 predicted targets, of which 68 were targeted by 3 or more differentially expressed miRNA. Gene ontology and pathway analysis of these predicted common targets showed enrichment for apoptosis and regulation of immune responses. Both of these pathways are known to be involved in luteolysis and survival of the CL. In addition, 46 novel miRNAs were detected in the CL during maternal recognition of pregnancy. Statistical analysis of the expression of the 46 novel miRNA resulted in 3 novel miRNAs that were differentially expressed in the CL of cyclic vs. pregnant animals. Target analysis of these sequences resulted in several genes involved in steroidogenesis, apoptosis and regulation of immune responses. The similar clustering of predicted target genes of both known and novel miRNAs provides strong support that the miRNAs differentially expressed during MRP are important in the regulation of luteal function and survival. Future studies to confirm changes in predicted targets and demonstrate functional effects on the CL will shed new light on the regulation of this dynamic tissue. Acknowledgements The authors express gratitude to Troy L. Ott, Manasi Manohar Kamat, and Sreelakshmi Vasudevan for facilitation of animal management and tissue collection. The authors also thank Ti-Cheng Chang and Zhihong Liu for advice in bioinformatics programs. This project was supported by Agriculture and Food Research Initiative Competitive Grant no. 2012-67015-30212 from the USDA National Institute of Food and Agriculture. Appendix: Supplementary material Supplementary data to this article can be found online at doi:10.1016/j.mce.2014.10.005. References Ahn, H.W., Morin, R.D., Zhao, H., Harris, R.A., Coarfa, C., Chen, Z.J., et al., 2010. MicroRNA transcriptome in the newborn mouse ovaries determined by massive parallel sequencing. Mol. Hum. Reprod. 16, 463–471. doi:10.1093/molehr/gap017. Ambros, V., 2008. The evolution of our thinking about microrNAs. Nat. Med. 14 (10), 1036–1040. doi:10.1038/nm1008-1036. Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biol. 11, R106. doi:10.1186/gb-2010-11-10-r106. Arvisais, E., Hou, X., Wyatt, T.A., Shirasuna, K., Bollwein, H., Miyamoto, A., et al., 2010. Prostaglandin F2 alpha represses IGF-I-stimulated IRS1/phosphatidylinositol-3kinase/AKT signaling in the corpus luteum: role of ERK and p70 ribosomal S6 kinase. Mol. Endocrinol. 24, 632–643. doi:10.1210/me.2009-0312. Atli, M.O., Bender, R.W., Mehta, V., Bastos, M.R., Luo, W., Vezina, C.M., et al., 2012. Patterns of gene expression in the bovine corpus luteum following repeated intrauterine infusions of low doses of prostaglandin F2alpha. Biol. Reprod. 86, 130. doi:10.1095/biolreprod.111.094870. Baley, J., Li, J., 2012. MicroRNAs and ovarian function. J. Ovarian Res. 5, 8. doi:10.1186/ 1757-2215-5-8. Bartel, D.P., 2004. MicroRNAs: genomics, biogenesis, mechanism and function. Cell 116 (2), 281–297. doi:10.1016/S0092-8674(04)00045-5. Bartel, D.P., 2009. MicroRNAs: target recognition and regulatory functions. Cell 136, 215–233. doi:10.1016/j.cell.2009.01.002. Bazer, F.W., Spencer, T.E., Johnson, G.A., Burghardt, R.C., 2009. Comparative aspects of implantation. Reproduction 138, 195–209. doi:10.1530/REP-09-0158.

87

Benyo, D.F., Pate, J.L., 1992. Tumor necrosis factor-alpha alters bovine luteal cell synthetic capacity and viability. Endocrinology 130 (2), 854–860. doi:10.1210/ endo.130.2.1733731. Berezikov, E., 2011. Evolution of miRNA diversity and regulation in animals. Nat. Rev. Genet. 12 (12), 846–860. doi:10.1038/nrg3079. Bott, R.C., Ashley, R.L., Henkes, L.E., Antoniazzi, A.Q., Bruemmer, J.E., Niswender, G.D., et al., 2010. Uterine vein infusion of interferon tau (IFNT) extends luteal life span in ewes. Biol. Reprod. 82, 725–735. doi:10.1095/biolreprod.109 .079467. Brennecke, J., Stark, A., Russell, R.B., Cohen, S.M., 2005. Principles of microRNA-target recognition. PLoS Biol. 3, e85. doi:10.1371/journal.pbio.0030085. Carambula, S.F., Pru, J.K., Lynch, M.P., Matikainen, T., Goncalves, P.B., Flavell, R.A., et al., 2003. Prostaglandin F2 alpha- and FAS-activating antibody-induced regression of the corpus luteum involves caspase-8 and is defective in caspase-3 deficient mice. Reprod. Biol. Endocrinol. 1, 15. doi:10.1186/1477-7827 -1-15. Davare, M.A., Saneyoshi, T., Guire, E.S., Nygaard, S.C., Soderling, T.R., 2004. Inhibition of calcium/calmodulin dependent protein kinase kinase by protein 14-3-3. J. Biol. Chem. 279, 52191–52199. Dillies, M.A., Rau, A., Aubert, J., Hennequet-Antier, C., Jeanmougin, M., Servant, N., et al., 2012. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief. Bioinform. 14 (6), 671–683. doi:10.1093/bib/bbs046. Ding, J., Zhou, S., Guan, J., 2011. MiRFam: an effective automatic miRNA classification method based on n-grams and a multiclass SVM. BMC Bioinformatics 12, 216. doi:10.1186/1471-2105-12-216. Djebali, S., Davis, C.A., Merkel, A., Dobin, A., Lassmann, T., et al., 2012. Landscape of transcription in human cells. Nature 489 (7414), 101–108. doi:10.1038/ nature11233. Donadeu, F.X., Schauer, S.N., Sontakke, S.D., 2012. Involvement of miRNAs in ovarian follicular and luteal development. J. Endocrinol. 215, 323–334. doi:10.1530/JOE12-0252. Farberov, S., Meidan, R., 2014. Functions and transcriptional regulation of thrombospondins and their interrelationship with fibroblast growth factor-2 in bovine luteal cells. Biol. Reprod. 91, 58. doi:10.1095/biolreprod .114.121020. Farh, K.K., Grimson, A., Jan, C., Lewis, B.P., Johnston, W.K., Lim, L.P., et al., 2005. The widespread impact of mammalian microRNAs on mRNA repression and evolution. Science 310, 1817–1821. doi:10.1126/science.1121158. Fiedler, S.D., Carletti, M.Z., Hong, X., 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. Gardner, P.P., Duab, J., Tate, J.G., Nawrocki, E.P., Kolbe, D.L., et al., 2009. Rfam: updates to the RNA families database. Nucleic Acids Res. 37, D136. doi:10.1093/nar/ gkn766. Goravanahally, M.P., Salem, M., Yao, J., Inskeep, E.K., Flores, J.A., 2009. Differential gene expression in the bovine corpus luteum during transition from early phase to midphase and its potential role in acquisition of luteolytic sensitivity to prostaglandin F2 alpha. Biol. Reprod. 80, 980–988. doi:10.1095/biolreprod. 108.069518. Griffiths-Jones, S., 2004. The microRNA registry. Nucleic Acids Res. 32, D109–D111. doi:10.1093/nar/gkh023. Griffiths-Jones, S., Bateman, A., Marshall, M., Khanna, A., Eddy, S.R., 2003. Rfam: an RNA family database. Nucleic Acids Res. 31, 439–441. doi:10.1093/nar/ gkg006. Griffiths-Jones, S., Moxon, S., Marshall, M., Khanna, A., Eddy, S.R., Bateman, A., 2005. Rfam: annotating noncoding RNAs in complete genomes. Nucleic Acids Res. 33, D121–D124. doi:10.1093/nar/gki081. Griffiths-Jones, S., Grocock, R.J., van Dongen, S., Bateman, A., Enright, A.J., 2006. MiRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 34, D140–D144. doi:10.1093/nar/gkj112. Griffiths-Jones, S., Saini, H.K., van Dongen, S., Enright, A.J., 2008. MiRBase: tools for microRNA genomics. Nucleic Acids Res. 36, D154–D158. doi:10.1093/nar/gkm952. Grimson, A., Farh, K.K., Johnson, W.K., Garrett-Engele, P., Lim, L.P., Bartel, D.P., 2007. MicroRNA targeting specificity in mammals: determinants beyond seed pairing. Mol. Cell 27, 91–105. doi:10.1016/j.molcel.2007.06.017. Hafner, M., Landgraf, P., Ludwig, J., Rice, A., Ojo, T., Lin, C., et al., 2008. Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods 44, 3–12. doi:10.1016/j.ymeth.2007.09.009. Herranz, H., Cohen, S.M., 2010. MicroRNAs and gene regulatory networks: managing the impact of noise in biological systems. Genes Dev. 24, 1339–1344. doi:10.1101/ gad.1937010. Hong, X., Luense, L.J., McGinnis, L.K., Nothnick, W.B., Christenson, L.K., 2008. Dicer1 is essential for female fertility and normal development of the female reproductive system. Endocrinology 149, 6207–6212. doi:10.1210/en .2008-0294. Hossain, M.M., Ghanem, N., Hoelker, M., Rings, R., 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. Hossain, M.M., Sohel, M.M.H., Schellander, K., Tesfaye, D., 2012. Characterization and importance of microrNAs in mammalian gonadal functions. Cell Tissue Res. 349, 679–690. doi:10.1007/s00441-012-1469-6. Huaiyu, M., Anushya, M., Thomas, P.D., 2012. PANTHER in 2013: modeling the evolution of gene function, and other gene attributes, in the context of phylogenetic trees. Nucleic Acids Res. 41, D377–D386. doi:10.1093/nar/ gks1118.

88

S.W. Maalouf et al./Molecular and Cellular Endocrinology 398 (2014) 78–88

Huaiyu, M., Anushya, M., Casagrande, J.T., Thomas, P.D., 2013. Large-scale gene function analysis with the PANTHER classification system. Nat. Protoc. 8, 1551–1566. doi:10.1038/nprot.2013.092. Huang, D.W., Sherman, B.T., Lempicki, R.A., 2009a. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi:10.1038/nprot.2008.211. Huang, D.W., Sherman, B.T., Lempicki, R.A., 2009b. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 3, 1–13. doi:10.1093/nar/gkn923. Huang, J., Ju, Z., Li, Q., Hou, Q., Wang, C., Li, J., et al., 2011. Solexa sequencing of novel and differentially expressed microRNAs in testicular and ovarian tissues in Holstein cattle. Int. J. Biol. Sci. 7 (7), 1016–1026. doi:10.7150/ijbs.7. 1016. Jimenez, B., Volpert, O.V., Crawford, S.E., Febbraio, M., Silverstein, R.L., Bouck, N., 2000. Signals leading to apoptosis-dependent inhibition of neovascularization by thrombospondin-1. Nat. Med. 6, 41–48. Juengel, J.L., Melner, M.H., Clapper, J.A., Turzillo, A.M., Moss, G.E., Nett, T.M., et al., 1998. Steady-state concentrations of mRNA encoding two inhibitors of protein kinase C in ovine luteal tissue. J. Reprod. Fertil. 113, 299–305. doi:10.1530/ jrf.0.1130299. Kozomara, A., Griffiths-Jones, S., 2011. MiRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 39, D152–D157. doi:10.1093/nar/ gkq1027. Kozomara, A., Griffiths-Jones, S., 2014. MiRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi:10.1093/nar/gkt1181. Lee, Y., Jeon, K., Lee, J.T., Kim, S., Kim, V.N., 2002. MicrorNA maturation: stepwise processing and subcellular localization. EMBO J. 21, 4663–4770. doi:10.1093/ emboj/cdf476. Lei, L., Jin, S., Gonzalez, G., Behringer, R.R., Woodruff, T.K., 2010. The regulatory role of Dicer in folliculogenesis in mice. Mol. Cell. Endocrinol. 315, 63–73. doi:10.1016/j.mce.2009.09.021. 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, B., Lee, H.-Y., Weinzimer, S.A., Powell, D.R., Clifford, J.L., Kurie, J.M., et al., 2000. Direct functional interactions between IGFBP-3 and RXR-a regulate transcriptional signaling and apoptosis. J. Biol. Chem. doi:10.1074/jbc. M002547200. Lytle, J.R., Yario, T.A., Steitz, J.A., 2007. Target mRNAs are repressed as efficiently by microRNA-binding sites in the 5’ UTR as in the 3’ UTR. Proc. Natl. Acad. Sci. U.S.A. 104, 9667–9672. doi:10.1073/pnas.0703820104. Ma, T., Jiang, H., Gao, Y., Zhao, Y., Dai, L., Xiong, Q., et al., 2011. Microarray analysis of differentially expressed microRNAs in non-regressed and regressed bovine corpus luteum tissue; microRNA-378 may suppress luteal cell apoptosis by targeting the interferon gamma receptor 1 gene. J. Appl. Genet. 52, 481–486. doi:10.1007/s13353-011-0055-z. Macias, S., Michlewski, G., Caceres, J.F., 2009. Hormonal regulation of microRNA biogenesis. Mol. Cell 36, 172–173. doi:10.1016/j.molcel.2009.10.006. Maeda, N., Kasukawa, T., Oyama, R., Gough, J., Firth, M., Engström, P.G., et al., 2006. Transcript annotation in FANTOM3: mouse gene catalog based on physical cDNAs. PLoS Genet. 2 (4), e62. doi:10.1371/journal.pgen.0020062. Maragkakis, M., Alexiou, P., Papadopoulos, G.L., Reczko, M., Dalamagas, T., Giannopoulos, G., et al., 2009a. Accurate microRNA target prediction correlates with protein repression levels. BMC Bioinformatics 10, 295. doi:10.1186/14712105-10-295. Maragkakis, M., Reczko, M., Simossis, V.A., Alexiou, P., Papadopoulos, G.L., Dalamagas, T., et al., 2009b. DIANA-microT web server: elucidating microRNA functions through target prediction. Nucleic Acids Res. 37, W273–W276. doi:10.1093/nar/ gkp292. Miles, J.R., McDaneld, T.G., Wiedmann, R.T., Cushman, R.A., Echternkamp, S.E., Vallet, J.L., et al., 2012. MicroRNA expression prolife in bovine cumulus-oocyte complexes: possible role of let-7 and miR-106a in the development of bovine oocytes. Anim. Reprod. Sci. 130, 16–26. doi:10.1016/j.anireprosci. 2011.12.021. Mishima, T., Takizawa, T., Luo, S.-S., Ishibashi, O., Kawahigashi, Y., Mizuguchi, Y., et al., 2008. MicroRNA (miRNA) cloning analysis reveals sex differences in miRNA expression profiles between adult mouse testis and ovary. Reproduction 136, 811–822. doi:10.1530/REP-08-0349. Morin, R.D., O’Connor, M.D., Griffith, M., Kuchenbauer, F., Delaney, A., Prabhu, A.-L., et al., 2008. Application of massively parallel sequencing to microRNA profiling and discovery in human embryonic stem cells. Genome Res. 18 (4), 610–621. doi:10.1101/gr.7179508. Nagaraja, A.K., Andreu-Vieyra, C., Franco, H.L., Ma, L., Chen, R., Han, D.Y., et al., 2008. Deletion of Dicer in somatic cells of the female reproductive tract causes sterility. Mol. Endocrinol. 22, 2336–2352. doi:10.1210/me.2008-0142. Niswender, G.D., Juengel, J.L., Silva, P.J., Rollyson, M.K., McIntush, E.W., 2000. Mechanisms controlling the function and life span of the corpus luteum. Physiol. Rev. 80, 1–29. Nor, J.E., Mitra, R.S., Sutorik, M.M., Mooney, D.J., Castle, V.P., Polverini, P.J., 2000. Thrombospondin-1 induces endothelial cell apoptosis and inhibits angiogenesis by activating the caspase death pathway. J. Vasc. Res. 37, 209–218.

Nothnick, W.B., 2012. The role of micro-RNAs in the female reproductive tract. Reproduction 143 (5), 559–576. doi:10.1530/REP-11-0240. Oliveira, J.F., Henkes, L.E., Ashley, R.L., Purcell, S.H., Smirnova, N.P., Veeramachaneni, D.N., et al., 2008. Expression of interferon (IFN)- stimulated genes in extrauterine tissues during early pregnancy in sheep is the consequence of endocrine IFN-tau release from the uterine vein. Endocrinology 149, 1252–1259. doi:10.1210/ en.2007-0863. Orom, U.A., Nielsen, F.C., Lund, A.H., 2008. MicroRNA-10a binds the 5’UTR of ribosomal protein mRNAs and enhances their translation. Mol. Cell 30, 460–471. doi:10.1016/j.molcel.2008.05.001. Otsuka, M., Zheng, M., Hayashi, M., Lee, J.D., Yoshimo, 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. Pate, J.L., 1995. Involvement of immune cells in regulation of ovarian function. J. Reprod. Fertil. Suppl. 49, 365–377. Petroff, M.G., Petroff, B.K., Pate, J.L., 2001. Mechanisms of cytokine-induced death of cultured bovine luteal cells. Reproduction 121 (5), 753–760. doi:10.1530/ rep.0.1210753. Place, R.F., Li, L.C., Pookot, D., Noonan, E.J., Dahiya, R., 2008. MicroRNA-373 induces expression of genes with complementary promoter sequences. Proc. Natl. Acad. Sci. U.S.A. 105, 1608–1613. doi:10.1073/pnas.0707594105. Poole, D.H., Pate, J.L., 2012. Luteal microenvironment directs resident T lymphocyte function in cows. Biol. Reprod. 86, 29. doi:10.1095/biolreprod.111.092296. Pratt, B.R., Butcher, R.L., Inskeep, E.K., 1977. Antiluteolytic effect of the conceptus and of PGE2 in ewes. J. Anim. Sci. 45, 784–791. Rege, T.A., Stewart, J., Jr., Dranka, B., Benveniste, E.N., Silverstein, R.L., Gladson, C.L., 2009. Thrombospondin-1-induced apoptosis of brain microvascular endothelial cells can be mediated by TNF-R1. J. Cell. Physiol. 218, 94–103. Ro, S., Song, R., Park, C., Zheng, H., Sanders, K.M., Yan, W., 2007. Cloning and expression profiling of small RNAs expressed in the mouse ovary. RNA 13, 2366–2380. doi:10.1261/rna.754207. Romero, J.J., Antoniazzi, A.Q., Smirnova, N.P., Webb, B.T., Yu, F., Davis, J.S., et al., 2013. Pregnancy-associated genes contribute to antiluteolytic mechanisms in ovine corpus luteum. Physiol. Genomics 45 (22), 1095–1108. doi:10.1152/physiolgenomics.00082.2013. Ruby, J.G., Stark, A., Johnston, W.K., Kellis, M., Bartel, D.P., Lai, E.C., 2007. Evolution, biogenesis, expression, and target predictions of a substantially expanded set of Drosophila microRNAs. Genome Res. 17, 1850–1864. doi:10.1101/gr. 6597907. Ryan, K.J., Engel, L.L., 1957. Hydroxylation of steroids at carbon 21. J. Biol. Chem. 225, 103–114. Silvia, W.J., Niswender, G.D., 1986. Maintenance of the corpus luteum of early pregnancy in the ewe. IV. Changes in luteal sensitivity to prostaglandin F2 alpha throughout early pregnancy. J. Anim. Sci. 63, 1201–1207. Sirotkin, A.V., Ovcharenko, D., Grossmann, R., Laukova, M., Mlyncek, M., 2009. Identification of microRNAs controlling human ovarian cell steroidogenesis via a genome-scale screen. J. Cell. Physiol. 219, 415–420. Tripurani, S.K., Xiao, C., Salem, M., 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. Vergoulis, T., Vlachos, I., Alexiou, P., Georgakilas, G., Maragkakis, M., Recsko, M., et al., 2012. Tarbase 6.0: capturing the exponential growth of miRNA targets with experimental support. Nucleic Acids Res. 40 (D1), D222–D229. doi:10.1093/nar/ gkr1161. Walusimbi, S.S., Pate, J.L., 2013. Physiology and endocrinology symposium: role of immune cells in the corpus luteum. J. Anim. Sci. 91, 1650–1659. doi:10.2527/ jas.2012-6179. Wiltbank, M.C., Diskin, M.G., Niswender, G.D., 1991. Differential actions of second messenger systems in the corpus luteum. J. Reprod. Fertil. Suppl. 43, 65–75. Yang, L., Wang, X.L., Wan, P.C., Zhang, L.Y., Wu, Y., Tang, D.W., et al., 2010. Up-regulation of expression of interferon-stimulated gene 15 in the bovine corpus luteum during early pregnancy. J. Dairy Sci. 93, 1000–1011. doi:10.3168/jds.2009-2529. Yang, Y., Tsui, H.-C.T., Man, T.-K., Winkler, M.E., 1998. Identification and function of the pdxy gene which encodes a novel pyridoxal kinase involved in the salvage pathway of pyridoxal 5’-phosphate biosynthesis in Escherichia coli K-12. J. Bacteriol. 180, 1814–1821. Yankey, S.J., Hicks, B.A., Carnahan, K.G., Assiri, A.M., Sinor, S.J., Kodali, K., et al., 2001. Expression of the antiviral protein Mx in peripheral blood mononuclear cells of pregnant and bred, non-pregnant ewes. J. Endocrinol. 170 (2), R7–R11. doi:10.1677/joe.0.170R007. Zalman, Y., Klipper, E., Farberov, S., Mondal, M., Wee, G., Folger, J.K., et al., 2012. Regulation of angiogenesis-related prostaglandin F2alpha-induced genes in the bovine corpus luteum. Biol. Reprod. 86, 92. Zheng, H., Fu, R., Wang, J.T., Liu, Q., Chen, H., Jiang, S.W., 2013. Advances in the techniques for the prediction of microRNA targets. Int. J. Mol. Sci. 14, 8179–8187. doi:10.3390/ijms14048179. Zou, Q., Mao, Y., Hu, L., Wu, Y., Ji, Z., 2014. miRClassify: an advanced webserver for miRNA family classification and annotation. Comput. Biol. Med. 45, 157–160. doi:10.1016/j.compbiomed.2013.12.007.