The aims of this study were to 1) identify the earliest transcriptional response of the bovine endometrium to the presence of the conceptus (using RNAseq), 2) investigate if these genes are regulated by interferon tau (IFNT) in vivo, and 3) determine if they are predictive of the pregnancy status of postpartum dairy cows. RNAseq identified 459 differentially expressed genes (DEGs) between pregnant and cyclic endometria on day 16. Quantitative real-time PCR analysis of selected genes revealed PARP12, ZNFX1, HERC6, IFI16, RNF213, and DDX58 expression increased in pregnant compared with cyclic endometria on day 16 and were directly upregulated by intrauterine infusion of IFNT in vivo for 2 h (P < 0.05). On day 13 following estrous endometrial expression of nine genes increased [ARHGAP1, MGC127874, LIMS2, TBC1D1, FBXL7, C25H16orf71, LOC507810, ZSWIM4, and one novel gene (ENSBTAT00000050193)] and seven genes decreased (SERBP1, SRGAP2, AL7A1, TBK1, F2RL2, MGC128929, and WBSCR17; P < 0.05) in pregnant compared with cyclic heifers. Of these DEGs, significant differences in expression between pregnant and cyclic endometria were maintained on day 16 for F2RL2, LIMS2, LOC507810, MGC127874, TBC1D1, WBSCR17, and ZSWIM4 (P < 0.05) both their expression was not directly regulated by IFNT in vivo. Analysis of the expression of selected interferon-stimulated genes in blood samples from postpartum dairy cows revealed a significant increase (P < 0.05) in expression of ZXFX1, PARP12, SAMD9, and HERC6 on day 18 following artificial insemination in cows subsequently confirmed pregnant compared with cyclic controls. In conclusion, RNAseq identified a number of novel pregnancy-associated genes in the endometrium of cattle during early pregnancy that are not regulated by IFNT in vivo. In addition, a number of genes that are directly regulated by short term exposure to IFNT in vivo are differentially expressed on day 18 following estrus detection in the blood of postpartum dairy cows depending on their pregnancy status.
- RNA sequencing
- gene expression
- estrus cycle
- interferon-stimulated gene
in order for successful pregnancy to be established and maintained in cattle, the corpus luteum (CL) must be maintained beyond the time when luteolysis would normally occur. It has been well established that the conceptus-derived pregnancy recognition signal in both cattle and sheep is interferon tau (IFNT) and in cattle, pregnancy recognition must occur by day 16 to ensure luteolysis does not occur (6, 36). IFNT exerts its effects both locally on the uterus as well as in a paracrine fashion on extra uterine tissues. A number of studies have demonstrated that IFNT alters the transcriptome of the endometrium during early pregnancy to increase expression of genes and/or transporters that enhance conceptus elongation and to establish uterine receptivity to implantation (reviewed by Ref. 46). These studies have shown that IFNT during the peri-implantation period of pregnancy in cattle induces the expression of large numbers of interferon-stimulated genes (ISGs; Refs. 3, 14, 19, 26, 31, 52). Indeed, the transcriptional response of the endometrium to the developing conceptus can act as a sensor of the developmental competency of the conceptus (5, 32). However, the majority of these studies have focused on the peri-implantation period of pregnancy during or after the period of pregnancy recognition in cattle (3, 14, 26, 31, 52).
Recent studies by our group and others sought to identify the differences between the endometrial transcriptome of pregnant and cyclic heifers at distinct developmental stages of embryo development. According to the Affymetrix microarray platform the first detectable differences between the endometria of pregnant and cyclic heifers occurred on day 15 (4) and 16 (14). Despite this, interaction between the conceptus and the maternal environment is required during the preimplantation stages of pregnancy in cattle (45), and it is known that the bovine embryo expresses mRNA for IFNT as early as the blastocyst stage of development (27, 29), indicating that the embryo may elicit localized or subtle responses from the endometrium prior to pregnancy recognition. In addition, transcriptomic profiling of bovine embryos at the initiation of conceptus elongation (11) as well as trophoblast cells cultured in vitro (12, 50) identified the expression of genes that code for secreted proteins. It is therefore reasonable to hypothesize that these conceptus-derived factors may stimulate more understated, localized gene expression changes in the endometrium during early pregnancy that are independent of IFNT.
Next-generation sequencing methods such as RNA sequencing (RNAseq) provide a far more precise measure of transcript abundance than other methods developed to date; as it is unbiased (no prior knowledge required), allowing complete coverage (i.e., very high sensitivity) and allowing true genome-wide discovery (33, 34). Data from human studies comparing the use of RNAseq and Affymetrix microarrays were highly correlated; however, the findings suggested that RNAseq was the more sensitive of the two methods as it identified 30% more differentially expressed genes (DEGs) than the microarray platform, as validated by quantitative real-time PCR (qPCR) (33). It has also been shown that data obtained using RNAseq are highly reproducible, and this holds true for both biological and technical replicates (34).
We hypothesized that a more sensitive transcriptomic platform such as RNAseq would identify subtle transcriptional responses of the endometrium to the conceptus prior to pregnancy recognition in cattle on day 16. The specific objectives of this study were 1) to determine if RNAseq could identify any changes in the endometrium of pregnant compared with cyclic heifers at earlier stages of embryo development than previously reported, 2) to determine what key biological pathways are affected in the endometrium during these stages of early pregnancy, 3) to identify if the expression of these genes is directly regulated by IFNT in vivo and 4) to determine if the expression of interferon-stimulated genes (ISGs) in blood are predictive of the pregnancy status of dairy cows.
MATERIALS AND METHODS
All experimental procedures involving animals were licensed by the Department of Health and Children, Ireland, in accordance with the Cruelty to Animals Act (Ireland 1876) and the European Community Directive 86/609/EC and were sanctioned by the Animal Research Ethics Committee of University College Dublin. Unless otherwise stated all chemicals were sourced from Sigma (Dublin, Ireland).
The aim of this experiment was to determine whether RNAseq technology can identify differences in endometrial gene expression between pregnant and cyclic heifers prior to maternal recognition of pregnancy in cattle.
Experimental design and tissue collection.
Tissue collection was performed as previously described as part of a larger experimental design (14). In brief, the estrus cycles of cross-bred beef heifers (predominantly Charolais and Limousin cross) were synchronized using a controlled internal drug release (CIDR) device (1.94 g progesterone; InterAg, Hamilton, New Zealand) placed into the vagina for 8 days with intramuscular administration of a prostaglandin F2α analog (equivalent to 0.5 mg of PG, Estrumate; Shering-Plough Animal Health, Hertfordshire, UK) 3 days prior to CIDR removal. Only heifers detected in standing estrous were used (n = 54). Thirty-eight heifers were inseminated, and 16 were left as noninseminated cyclic controls. Heifers were then randomly assigned to be slaughtered on day 13 or 16 of early pregnancy/estrous cycle corresponding to the initiation of conceptus elongation (day 13) and maternal recognition of pregnancy (day 16) in pregnant heifers, respectively. All reproductive tracts were processed ∼30 min postslaughter, whereby the uterine horn ipsilateral to the CL was flushed with phosphate-buffered saline containing 5% fetal calf serum. In inseminated heifers, the presence of the conceptus was determined by examination under a stereomicroscope. Tissue was collected only from heifers from which a conceptus at the appropriate stage of development was recovered.
After being flushed, the uterine horn ipsilateral to the CL was opened longitudinally, and 300 mg strips of endometrium (predominantly intercaruncular tissue) were immersed in 1:5 wt/vol of RNAlater. This tissue was transported back to the laboratory on ice, stored at 4°C for 24 h, removed from RNAlater, and stored at −80°C prior to RNA extraction. RNA extraction was performed using the TRIzol reagent (Invitrogen, Carlsbad, CA) as per manufacturer's instructions and RNA cleanup and on column DNase digestion performed (Qiagen, Crawley, West Sussex, UK). RNA quality and quantity were determined with the Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA), and only samples with an RNA integrity number (RIN) of >8.5 were processed further (n = 5 samples per treatment per time point).
RNA sequencing tissue processing.
Unless otherwise stated, all materials were obtained from Illumina (http://www.illumina.com). Analysis of endometrial gene expression was carried out on the Illumina GA2 sequencer using the standard Illumina protocol for sequencing cDNA samples. We diluted 10 μg of total RNA to a volume of 50 μl, and this was incubated at 65°C for 5 min to disrupt the secondary structures. Following washing with binding buffer, total RNA was added to the Dynal oligo(dT) beads and rotated at room temperature for 5 min. The mRNA fraction was eluted in 10 mM Tris·HCl by incubation at 80°C for 2 min and subsequently at 65°C for 5 min; this was repeated twice to ensure maximal mRNA yield from the sample. mRNA was then fragmented by addition of 1 μl of fragmentation buffer (Ambion, Carlsbad, CA) to 9 μl of mRNA and incubated at 70°C for 5 min. Then, 1 μl of stop buffer was added, and the sample immediately placed on ice. To elute the fragmented mRNA, 3 M NaOAc (pH 5.2), glycogen (5 μg/μl), and 100% EtOH were added, and samples were incubated at −80°C for 30 min, spun at 4°C for 25 min at 14,000 rpm; the pellet washed with 70% EtOH and then air-dried. The fragmented mRNA was then resuspended in RNase DNase-free H2O. First-strand synthesis was performed using random hexamers (3 μg) and Superscript II (Invitrogen) by incubating the fragmented mRNA in a thermocycler (G-Storm, Surrey, UK) at 25°C for 10 min, 42°C for 50 min, 70°C for 15 min, and placed on ice. This was then used for second-strand cDNA synthesis using RNaseH (2 U) and DNA Pol I (50 U), which was incubated at 16°C in a thermocycler for 2.5 h and purified using the QIAquick PCR spin column as per manufacturer's instructions (Qiagen).
Library preparation and cluster generation.
The ends of the double-stranded cDNA fragments were blunted by incubation at 20°C for 30 min with T4 DNA Polymerase (15 U), Escherichia coli DNA Polymerase I Klenow Fragment (5 U), and T4 PNK (50 U) to remove the 3′ overhangs and fill in the 5′ overhangs. The DNA was then purified using the QIAquick PCR spin column. A single “A” base was then added to the 3′ end of the fragmented cDNA by incubation at 37°C for 30 min with 15 U of Klenow 3′-5′ exo-enzyme, and the DNA was further purified using the QIAquick MinElute column (Qiagen). Ligation of the proprietary Illumina genomic adapters to the DNA fragments was performed by incubation for 15 min at room temperature with T4 DNA Ligase in a 10:1 molar ratio of adapter to genomic DNA insert. The DNA was purified with the QIAquick MinElute column and gel purification performed (2% agarose 400 ng/ml ethidium bromide) to remove all unligated adapters. A gel band of 200 bp was excised and purified with the QIAquick gel extraction kit. Finally, PCR enrichment of the purified adapter ligated DNA was performed using Phusion polymerase (New England Biolabs, Ipswich, MA) and PCR primers under the following cycle conditions: 98°C for 30 s, 15 cycles of 98°C for 10 s, 65°C for 30 s 72°C for 30 s, and 72°C for 5 min. The PCR product was then purified with a QIAquick column, and the quality and quantity of the purified product determined using the Agilent Bioanalyzer DNA 1000 chip and the Qubit (Invitrogen) prior to loading the flow cell. Cluster generation on a single-read flow cell was carried out on the Illumina cluster generation station according to the manufacturer's instructions and each flow cell underwent 36 cycles of sequencing on an Illumina GA2 sequencer.
RNAseq data analysis.
The RNAseq samples were processed through the standard software pipeline for the Genome Analyzer (http://bioinfo.cgrb.oregonstate.edu/docs/solexa/SCS2_01_IPAR1_01_Release_Notes.pdf). The CASAVA module from Illumina software was used to process RNAseq data. Data were aligned against the BosTau4 genome, and a pseudo chromosome containing potential splice junction sequences was generated. These sequences contained 32 bases from the end of one exon and 32 bases from the start of another exon. At least a 4 bp overlap on the shorter side of the splice junction was required, and all potential splice variants for single transcripts were generated. The ensGene table from the UCSC genome browser (http://hgdownload.cse.ucsc.edu/goldenPath/bosTau4/database/ensGene.txt.gz: Oct 2007 BosTau4) was used to provide exon location information to the CASAVA module. The moderated negative binomial test from the edgeR Bioconductor library (41) was used to generate the lists of differentially expressed RNAseq transcripts which were displayed as transcripts per million. An false discovery rate (FDR)-adjusted P value of <0.05 was used as the cut-off for determining significance. The comparative analysis was restricted to the 26957 protein coding transcripts in version 52 of Ensembl (http://www.ensembl.org).
Quantitative real-time PCR.
Quantitative real-time PCR (qRT-PCR) was used to validate the findings of the RNAseq data. In addition to the samples used for RNAseq analysis, independent samples in each treatment group were used for validation giving n = 8 per treatment (pregnant vs. cyclic) per time point (day 13 vs. 16). In brief, 5 μg of total RNA was converted to complementary DNA (cDNA) using Superscript III (Invitrogen) and random hexamers as per manufacturer's instructions. All primers were designed using Primer-BLAST software (Supplementary Table S1 and http://www.ncbi.nlm.nih.gov/tools/primer-blast/).1 Each qRT-PCR reaction was carried out on the 7500 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA) with 50 ng of cDNA, at a concentration of 300 nM, 10 μl Sybrgreen mastermix (Applied Biosystems) in a final reaction volume of 20 μl. Cycling conditions were as follows: 2 min at 50°C, 10 min at 95°C, and 40 cycles of 95°C for 15 s, and 60°C for 1 min and were carried out with the inclusion of a dissociation curve to ensure specificity of amplification. A standard curve was included for each gene of interest as well as for the normalizer gene to obtain primer efficiencies. Analysis of transcript abundance was determined by the standard curve method by comparing the gene of interest to the reference gene (RPL19). Significance differences between treatments were determined using Student's t-test with significance set as P < 0.05.
Gene ontology analysis.
Analysis of the gene ontology (GO) terms was performed using the GOstats package of Bioconductor (13). The RNAseq reads were first filtered as outlined in the GOstats vignette or help-page in Bioconductor. This provided a “gene universe” that represented the set of expressed genes in all experimental conditions i.e., pregnant and cyclic endometria on both day 13 and day 16. Data were filtered to reduce the number of “false positives” resulting from the analysis, i.e., those GO terms that were marked as statistically significant when in truth they are not. For each list of identified DEGs from a given comparison e.g., pregnant versus cyclic heifers on day 16, a conditional hypergeometric statistical test was performed using a cut-off P value of 0.01. This generated the overrepresented GO nodes, i.e., those associated with the probe list more than would be expected by chance based on the gene universe.
Ingenuity pathway and interaction network analysis.
Ingenuity pathway analysis (IPA; http://www.ingenuity.com) was carried out for each list of DEGs to determine the functions and/or pathways associated with each of the gene lists and to identify which molecules from a given list of DEGs interacted with one another. To identify the pathways that were associated with each list of DEGs a P value for a given function was calculated by considering the number of functional analysis molecules that participate in that function and the total number of molecules that are known to be associated with that function/pathway in the Ingenuity Knowledge Base. Functions and pathways with a P value <0.05 were considered significant for a given list of DEGs. The ratio indicates the number of genes for a particular pathway in a gene list compared with the total number of genes in a given pathway.
Interaction networks for each list of DEGs were generated by identifying genes that serve as molecules of interest that interact with other molecules in the Ingenuity Knowledge Base. These serve as “seed” molecules to generate a network with a limit of 35 molecules for each interaction network generated. The networks generated were scored based on the number of network-eligible molecules they contained. The network score was based on the hypergeometric distribution and was calculated using a right-tailed Fisher's exact test. The higher the score for a given interaction network, the lower the probability of finding the observed number of network eligible molecules in a given network by chance.
The aim of this experiment was to confirm in vivo if the DEGs identified in the endometrium in experiment 1 were specifically induced by conceptus-derived IFNT and to localize their expression to either the caruncular or intercaruncular region of the endometrium.
Infusion of IFNT to the uterine endometrium in vivo.
On day 14 of a synchronized estrous cycle nonbred heifers received a single intrauterine infusion (25 ml/horn) of 0.9% NaCl solution containing recombinant ovine IFNT (200 μg/ml) or vehicle as previously described (43). Following 2 h of infusion, heifers were slaughtered (n = 5 per treatment), and the uterine horn where infusion took place was opened longitudinally. The caruncular and intercaruncular regions were carefully dissected from one another, immediately snap-frozen in liquid nitrogen, and stored at −80°C until being processed for qRT-PCR analysis as described for experiment 1.
The aim of this experiment was to confirm if ISGs are detectable in blood during early pregnancy in postpartum dairy cows and may be predictive of the pregnancy status of a cow.
Animal model and gene expression analysis.
Dairy cows from a commercial herd (Lyons Research Farm) ∼80 days postpartum had their estrous cycles synchronized as described for experiment 1. All cows were observed for estrus, and only those observed were utilized further (estrus = day 0). Cows were then randomly assigned to either cyclic controls (n = 4) or were inseminated (n = 16). We collected 3 ml of whole blood by jugular vein puncture into Tempus Blood RNA Tubes (Applied Biosystems) on day 13, 16, 18, and 25 postestrus detection of the synchronized cycle and stored at 4°C for 3 h. All inseminated cows were scanned via transrectal ultrasonogrophy on day 35 (Aloka SSD-900; ALOKA, Zug, Switzerland) and further subdivided into those that were inseminated and pregnant and inseminated and nonpregnant.
Following storage of blood samples at 4°C, RNA was extracted from whole blood as per manufacturer's instructions (Applied Biosystems). Following DNase treatment, quality and quantity of RNA were determined using the Agilent Bioanalyzer (Agilent Technologies) and the NanoDrop 1000 (Thermo Fischer Scientific, Wilmington, DE). All samples had RINs of >8.0. First-strand cDNA synthesis was performed from 1 μg of total RNA using the high-capacity cDNA reverse transcription kit as per manufacturer's instructions (Applied Biosystems). qRT-PCR analysis was carried out a described for experiments 1 and 2 on genes known to be directly regulated by IFNT from experiment 2 as well as previously published results (14).
Data were checked for normality and homogeneity of variance using histograms, qplots, and formal statistical tests in the UNIVARIATE procedure (version 9.1.3; SAS Institute, Cary, NC). Data that were not normally distributed were transformed by raising the variable to the power of lambda. The appropriate lambda value was obtained by conducting a Box-Cox transformation analysis using the TRANSREG procedure of SAS. The transformed data were used to calculate P values. The corresponding least squares means and SE of the nontransformed data are presented in the results for clarity. A repeated-measures mixed-model ANOVA (MIXED procedure of SAS) was conducted to determine the effect of pregnancy status and time on the relative expression level for each gene of interest. Fixed effects included pregnancy status (pregnant or cyclic), time, and their interaction. The interaction term, if not statistically significant (P > 0.10), was subsequently excluded from the final model. Animal within status was included as a random effect in the model, with the most appropriate covariance structure between records within cow determined by minimizing the Akaike information criterion (AIC). Models were run under compound symmetry, unstructured, autoregressive, or Toeplitz variance-covariance structures. The model with the least AIC value was selected. Differences between treatments were determined by F-tests using type III sums of squares. The PDIFF command incorporating the Tukey test was applied to evaluate pairwise comparisons between treatment means.
Overall Expression of Transcripts in Pregnant and Cyclic Endometria Detected Using RNAseq and Comparison With Affymetrix Microarray Analysis
In total, RNAseq detected an average of 5.5 million reads in endometria from pregnant and cyclic heifers with a range of 3.2 million–13.1 million. Approximately 85% of these reads aligned to the bovine genome and 21,547 out of 28,054 Ensembl transcripts identified with at least 1 read in 1 sample. Only RNAseq (and not microarray analysis) detected differences between pregnant and cyclic heifers on day 13. Comparing the expression of the same samples (n = 5 pregnant and n = 5 cyclic) using the Affymetrix microarray platform revealed only 58 differentially expressed Ensembl transcripts between pregnant and cyclic endometria on day 16 in contrast to the 459 identified in the same sample set by RNAseq.
Differences in Endometrial Gene Expression Between Pregnant and Cyclic Heifers Identified With RNAseq
Analysis of data derived from day 13 tissue identified 16 DEGs between pregnant and cyclic heifers (Fig. 1). The expression of nine of these genes increased in pregnant compared with cyclic endometria [P < 0.05: ARHGAP1, MGC127874, LIMS2, TBC1D1, FBXL7, C25H16orf71, LOC507810, ZSWIM4, and one novel gene (ENSBTAT00000050193)], while the expression of seven others decreased (P < 0.05: SERBP1, SRGAP2, AL7A1, TBK1, F2RL2, MGC128929, and WBSCR17; Supplementary Table S2). Given the small number of DEGs identified, no overrepresented gene ontologies, IPA pathways or interaction networks were associated with these genes.
In contrast to day 13, comparison of endometria from pregnant and cyclic heifers on day 16 identified 459 DEGs of which 283 increased in expression and 176 decreased in pregnant heifers compared with their cyclic counterparts (Supplementary Table S3). Using the GOStats vignette of Bioconductor, we categorized the overrepresented GO terms for these DEGs into the biological processes of cellular macromolecule metabolic process, protein metabolic process, and cellular protein metabolic process, among others (Table 1). IPA was used to categorize these 459 DEGs into pathways and interaction networks. The top five pathways associated with DEGs on day 16 of pregnancy were involved in 1) interferon signaling (10 genes from a total of 36 in the pathway: 10/36), 2) role of pattern recognition receptors in recognition of bacteria and viruses (11/89), 3) activation of IRF by cytosolic pattern recognition receptors (10/74), 4) retinoic acid-mediated apoptosis signaling (8/70), and 5) role of PKR in interferon induction and antiviral response (6/46). Interaction networks, which identify direct and indirect interaction between DEGs, were categorized into the biological processes of 1) antimicrobial response, inflammatory response, and infection mechanism, 2) infection mechanism, antimicrobial response, and inflammatory response, 3) inflammatory response, antimicrobial response, and infection mechanism, 4) cell morphology, cellular compromise, hepatic system development, and function, as well as 5) connective tissue disorders, immunological disease, and inflammatory disease (Table 1).
qRT-PCR Validation of RNAseq Data
A number of the genes identified as differentially expressed on day 16 have been previously validated by qRT-PCR (14). However, we chose to further validate a panel of genes (18 in total) with the largest fold change difference between pregnant and cyclic heifers, as well as some genes that have not been previously characterized as altered as early as day 16 of pregnancy in the bovine endometrium (Fig. 2A). Of the 10 genes that decreased in expression in pregnant compared with cyclic endometria on day 16, five (SLC34A2, CA8, PLVAP, ADAMTS13, and ALAD; P < 0.05) were validated by qRT-PCR. Although expression of the remaining five genes was lower in pregnant compared with cyclic heifers, this was not statistically significant (LDB1, MBTPS1, MATN2, PER1, and TRAF7; P > 0.05). Of the eight genes whose expression significantly increased in pregnant heifers, all eight validated (PARP12, ZNFX1, HERC6, IFIT1, IFI16, AGRN, RNF213, and DDX58; P < 0.05) following qRT-PCR analysis. Furthermore, infusion of IFNT into the uterine horn of cyclic heifers for 2 h directly upregulated the expression of HERC6, ZNFX1, PARP12, IFI16, RNF213, and DDX58 in the caruncular region and HERC6, IFI16, RNF213, and DDX58 in the intercaruncular region of the endometrium (Fig. 2, B and C, respectively).
Validation of Novel Conceptus-induced Genes in the Endometrium
Given that the 16 DEGs were identified between pregnant and cyclic heifers on day 13 we chose to validate their expression by qRT-PCR (Fig. 3A). The expression of MGC127874, LIMS2, TBC1D1, LOC507810, and ZSWIM4 increased significantly (P < 0.05) in pregnant endometria, while the expression of both MGC128929 and TBK1 approached significance (P = 0.08 and P = 0.07, respectively). Of the seven genes identified as decreased in pregnant compared with cyclic heifers, only WBSCR17 and F2RL2 validated by qRT-PCR (P < 0.05).
Of the nine DEGs on day 13 that were validated by qRT-PCR, to further confirm that the expression of these genes was altered due to the presence of a conceptus, we measured their expression in pregnant and cyclic endometria on day 16. Expression of both WBSCR17 and F2RL2 was significantly decreased in pregnant compared with cyclic heifers on day 16 (P < 0.05, Fig. 3B). The expression of LIMS2, TBC1D1, LOC507810, ZSWIM4, and MGC127874 increased significantly in day 16 pregnant versus cyclic endometria (P < 0.05). Although the expression of MGC128929 and TBK1 did increase in pregnant compared with cyclic heifers, the difference was not statistically significant (P = 0.07 and P = 0.09, respectively; Fig. 3B).
Following on from qRT-PCR validation of genes whose expression increased in pregnant compared with cyclic endometria as early as day 13 and day 16, we aimed to identify whether these genes were directly regulated by IFNT in vivo. There was no effect of IFNT infusion on the expression of LIMS2, TBC1D1, ZSWIM4, MGC127874, MGC128929, or TBK1 in either the caruncular or intercaruncular regions of the endometrium (P > 0.1), confirming these genes are not regulated by conceptus derived IFNT in vivo. However, the expression of LOC507810 increased significantly, while MGC127874 approached significance (P = 0.07) in the caruncular region following IFNT infusion.
Pregnancy Status Effects on ISGs in Blood
Of the 16 cows that were inseminated, six were confirmed as pregnant on day 35 following estrus detection. There was no significant effect of pregnancy status or any day × status interaction on the expression of MGC127874, LOC507810, BST2, EIF4E, IFI16, IFI44, ISG20, or UCRP (P > 0.05). These data were then merged, and the effect of day was determined. There was no significant day effect on the expression of MGC127874, LOC507810, or IFI16 (P > 0.05, Table 2); however, in all three groups BST2, IFI44, and ISG20 expression increased significantly on day 18 compared with day 16 (P < 0.05). In addition, EIF4E expression was higher (P < 0.05) on day 18 and 25 compared with early days of the estrous cycle/pregnancy, while UCRP expression significantly increased on day 25 in all three groups.
The expression of five genes (HERC6, PARP12, SAMD9, RNF213, and ZNFX1) was temporally altered and a significant day × pregnancy status interaction observed. HERC6 expression increased on day 18 in all three groups; however, on day 25, HERC6 expression was significantly lower in controls compared with either of their artificially inseminated (AI) counterparts (P < 0.05, Fig. 4A).
The expression of both PARP12 and RNF213 increased on day 18, and there was significantly higher expression of both these genes in the AI-confirmed pregnant compared with the cyclic control and AI nonpregnant groups (P < 0.05; Fig. 4, B and C). On day 25, however, the expression of both PARP12 and RNF213 was only significantly different between AI groups and controls.
Both SAMD9 and ZNFX1 also displayed peak expression on day 18 of the estrous cycle/pregnancy, which reduced significantly on day 25 for SAMD9 expression only (Fig. 4, D and E). The expression of both these genes was higher in cyclic controls on day 13 compared with both AI groups (P < 0.05), but by day 18 expression of SAMD9 and ZNFX1 was significantly higher in the AI confirmed pregnant group compared with controls: an expression pattern that was maintained only for ZNFX1 on day 25 where expression of ZNFX1 was higher in both AI groups compared with cyclic controls (P < 0.05).
In this study we hypothesize that subtle changes in gene expression in the endometrium are detectable prior to maternal recognition of pregnancy on day 16. Using next-generation sequencing we have identified 1) novel differentially expressed genes on day 16 of pregnancy, some of which are directly regulated by IFNT in vivo, 2) the networks and pathways associated with genes that are altered in the endometrium during conceptus elongation, and 3) the expression of ZSWIM4, TBC1D1, LIMS2, TBK1, LOC507810, MGC127874, and MGC128929 as being increased in pregnant compared with cyclic endometria as early as day 13 of pregnancy and 4) the expression of only LOC507810 as being directly regulated in the caruncular region of the endometrium by short-term exposure to IFNT in vivo.
Previous studies from our group (14) using a microarray approach demonstrated that the first detectable differences in endometrial gene expression between pregnant and cyclic heifers were observed on day 16, coincident with pregnancy recognition (6, 36). This previous microarray analysis involved the use of a larger cohort of heifers, including those that had been exposed to a high progesterone environment (14), resulting in advanced conceptus elongation and identification of 764 DEGs between pregnant and cyclic heifers on day 16. However, analysis of the microarray data from the same cohort of heifers utilized for RNAseq analysis identified only 58 DEGs (Supplementary Table S4). This indicates that the advanced conceptus elongation contributed significantly to the pregnancy recognition transcriptomic signature observed in the endometrium of pregnant compared with cyclic heifers and is consistent with published data at later stages of pregnancy (2, 31, 52). The data presented in this study demonstrate the increased sensitivity of RNAseq given it identified up to eight times more DEGs in pregnant compared with cyclic endometria on day 16 in heifers where no manipulation of the progesterone status of the animal was performed.
Of the 402 DEGs identified on day 16 between pregnant and cyclic endometria by RNAseq, a number had previously been shown to be regulated by IFNT in vivo (14). Nonetheless, we chose to determine whether or not IFNT regulated the expression of seven additional genes (14). IFI16 (interferon-inducible protein 16) is a member of the family of p-200 proteins originally identified in human lymphoid cell lines. Its expression can be induced by interferon gamma (IFNG), and it functions to negatively regulate cell growth (reviewed by Ref. 1). The expression of IFI16 has been reported previously as increased in pregnant endometria of cattle (day 15, 17, 18, and 20; Refs. 3, 4, 26, 31, 52), and recent studies have demonstrated that infusion of IFNA into the uterine lumen for 3 days increases IFI16 expression. However, the current study using RNAseq determined that the expression of IFI16 increased on day 16 of pregnancy due to the direct induction by IFNT in vivo, i.e., following infusion of IFNT for only 2 h in vivo. Similarly, increases in DDX58 expression have been identified in endometria of pregnant compared with cyclic heifers (26, 31, 52), as well as regulation by IFNA in vivo in cattle (4). DDX58 is a member of the DExH box family encodes for a protein that contains two caspase recruitment domains (CARDs) at the NH2 terminus and an RNA helicase domain at the COOH terminus. In human endometrial cells its expression is regulated by IFNG (55) and addition of recombinant ovine IFNT to 2fTGH (parental) cells derived from a human fibrosarcoma increased DDX58 expression in a STAT1-dependent manner. Both the spatial and temporal expression of this gene in the ovine endometria during early pregnancy has been shown to be directly regulated by IFNT (44). RNAseq has identified that the presence of the conceptus on day 16 induces expression of this gene in the bovine endometrium, and this expression is directly regulated by conceptus-derived IFNT. These data demonstrate that the expression of these genes in the endometrium are a result of the type 1 IFN signaling cascade induced during the process of pregnancy recognition in cattle.
Interestingly, the additional four genes directly regulated by IFNT in vivo on day 16 of pregnancy are all involved or potentially involved in the process of ubiquitination. Ubiquitination, the cellular process whereby proteins are degraded, requires three enzymes: the E1, which is ubiquitin-activating enzyme; the E2 proteins, which are the ubiquitin conjugase enzymes; and finally the E3 proteins, known as the ubiquitin ligase enzymes. These E3 enzymes guarantee specific substrate recognition, and this large family of proteins comprises two main subtypes, the HECT and RING containing proteins. Notably this process of ubiquitination is utilized by the innate immune system to elicit type 1 IFN responses in cells as well as mediating other processes.
PARP12 [poly(ADP-ribose) polymerase (PARP) superfamily 12], is one of three members of the PARP superfamily that are involved in ADP-ribosylation. This gene has previously been shown to be upregulated in the bovine endometrium during the peri-implantation period (26, 52) and localized to the superficial and deep uterine glands (26). It is thought that PARP12 is involved in protein modification and cell cycle/growth phase regulation. PARP12 was shown by Bauersachs et al. (3) to contain a WWE domain that is normally encoded by genes whose proteins are involved in ubiquitination. HERC6, a low molecular weight member of the HERC family, contains a HECT domain and an RCC1-like domain (RLD). This RLD domain is a guanine nucleotide exchange factor for Ran and participates in transport of proteins from the nucleus to the cytoplasm and is involved in mitotic spindle formation. HERC6 has been previously detected in human placenta and uterine microvascular endothelial and smooth muscle cells; however, its expression is comparably low compared with other members of the HERC family (23). Unlike human HERC6, mouse HERC6 has been shown to be an ISG15 E3 ligase, meaning that its activity is ligase dependent and it is involved in the process of ISG15-ylation (51). According to the limited literature available for RNF213, this protein is a member of the RING family of E3 ligase enzymes and is a zinc-binding domain that is involved in protein-protein interactions and displays E3 ubiquitin ligase activity (49). ZNFX1 (zinc finger, NFX1-type containing 1) expression is increased during viral infection in the brain of Atlantic cod as well as other known ISGs (ISG15, RSAD2, B2M; Ref. 40), indicating it functions in the antiviral process of type 1 interferon signaling. The current study found both ZNFX1 and RNF213 to be directly regulated by IFNT in caruncular regions in vivo. Collectively, the increase in expression of these genes (PARP12, HERC6, RNF213, and ZNFX1) during early pregnancy and their direct regulation by IFNT in vivo may point toward a role for genes in the ubiquitination/ISG15-ylation process in the endometrium during early pregnancy, where it is known that restricted expression of genes/proteins, in particular those involved in immune response, is required to establish and maintain uterine receptivity to implantation (46).
The data presented involving differences between pregnant and cyclic endometria on day 16, as well as previous reports in the literature on early pregnancy in cattle (14, 31), demonstrate that the predominant effector of change in the pregnant endometrial transcriptome is due to type 1 IFN signaling, which is predominantly a consequence of IFNT secretion by the conceptus. It is known that the embryo expresses mRNA for IFNT as early as the blastocyst stage of development (27, 29), and it is probable that prior to pregnancy recognition, the embryo may elicit local effects on the endometrium that have yet to be identified, and the increased sensitivity achievable with RNAseq may allow such subtle changes to be detected. The fact that a number of DEGs were identified on day 13 by RNAseq is consistent with previous reports comparing different transcriptomic and proteomic platforms (54); additional novel genes and pathways were identified, possibly due to the reduced redundancy and bias, as well as increased sensitivity of RNAseq (33). Validation of these novel DEGs identified on day 13 has shown some discrepancies between the RNAseq and qRT-PCR platforms. This occurred for genes where the fold change differences between pregnant and cyclic endometria were quite small and has previously been shown to occur in genes with low expression values (33, 34). Despite these inconsistencies, the expression of LIMS2, TBC1D1, ZSWIM4, MGC127874, LOC507810, MGC128929, and TBK1 was validated by qRT-PCR analysis on day 13, and evidence that these differences were due to the presence of the conceptus was provided by maintained increase in their expression on day 16 in pregnant compared with cyclic heifers. However, none of these genes were directly regulated by short-term exposure to IFNT in vivo, and a survey of data from Bauersachs et al. (4) confirms that longer-term infusion of a type I interferon (for 3 days in vivo) does not affect the expression of these genes.
The gene LIMS2 (LIM and senescent cell antigen-like domains 2) encodes for a protein that consists of five LIM domains, each containing two zinc-finger domains in tandem that are separated by a two-amino acid residue hydrophobic linker and have been identified in all eukaryotes whose genomes have been extensively studied (25). LIMS2 competes with LIMS1 for binding with integrin-linked protein kinase; however, little is known about this protein in cattle. The murine homologs Pinch1 and Pinch2 are detectable in mouse uterine homogenate (9), and Pinch2 protein is localized to focal adhesions (47). Pinch2 null mice are fertile and develop normally, whereas Pinch1 are embryonic lethal (28), and cells that do not express Pinch1 show a reduced capacity for cellular adhesion when tested with a large number of extracellular matrix substrates. Both Pinch 1 and 2 are thought to stabilize integrin-linked protein kinase in focal adhesions in particular those containing the integrin subunits-β1 and -β3 in mice (reviewed by Ref. 53). Therefore, it is possible that increased LIMS2 expression as early as day 13 in pregnant endometria may act in establishing uterine receptivity, including binding to the integrin subunit-β1 expressed in the luminal epithelium as well as the stroma later in pregnancy (30) to establish the formation of focal adhesions.
TBC1D1 [TBC1 (tre-2/USP6, BUB2, cdc16) domain family, member 1] prevents the stimulation of a member of the facilitative glucose transporter family, SLC2A4, which is known to actively transport glucose (reviewed by Ref. 10). Previous studies carried out in sheep have shown that a number of facultative and sodium-dependent glucose transporters are expressed in the endometrium and are regulated by progesterone and/or IFNT (16). Both SLC2A4 mRNA and protein localized to the luminal and glandular epithelium as well as the stroma and its expression seems to be increased by the coordinate actions of progesterone and IFNT but not IFNT alone (16). Moreover, total recoverable glucose is increased in the lumen of pregnant compared with cyclic ewes (17). It seems counterintuitive that TBC1D1 expression is increased in pregnant endometria on day 13 and day 16, preventing availability of SLC2A4 to transport glucose into the uterine lumen. However, other members of the glucose transporter family are expressed in the endometrium of cattle, some of which are altered by progesterone and pregnancy (15) and increased TBC1D1 expression may be involved in the tight regulation of availability of glucose transporters to regulated glucose availability in the uterine lumen to the developing conceptus.
ZSWIM4 (zinc finger, SWIM-type containing 4) is thought to be involved in zinc and metal ion binding (42). No data are available for ZSWIM4 in the literature; however, a murine homolog zswim2 is contained in the MEX protein and is required for MEX E3 ubiquitin ligase activity in the testes (35). This is of interest as if other members of the ZSWIM family such as ZSWIM4 are involved in the ubiquitination pathway; it is possible that ZSWIM4 plays a role in the ubiquitination pathway in the bovine endometrium as, in the present study by day 13 of pregnancy, ZSWIM4 is differentially expressed. It has been proposed that ubiquitin cross-reactive protein and other proteins involved in the ubiquitination pathway are induced by IFNT in the pregnant endometrium to tag endometrial cytosolic proteins for degradation to enhance uterine receptivity to implantation (24, 38, 39). Further evidence that ZSWIM4 may be an early mediator of this protein degradation process is its increased expression is maintained to day 16, although further investigation into its regulation is required.
The transcript MGC127874 (p60MONOX brain-derived rescue factor) is protein coding, and interaction studies have shown that it interacts with both the prion protein PrP and alpha synuclein (SNCA), which is a component of the Parkinson's and Alzheimer's pathways. The LOC507810 gene has also been shown to be protein coding, but so far nothing is known of its function. What is clear is these genes seem to be directly regulated by IFNT infusion for 2 h in vivo and may mediate the early IFNT response in the endometrium. While IFNT does affect the expression of genes to the greatest extent in the endometrium during early pregnancy, it is known that a number of genes are altered in the transition of the embryo from a blastocyst to an elongated filamentous conceptus (reviewed by Ref. 7). It is possible that additional conceptus-derived proteins, e.g., prostaglandins, may alter the expression of some of these genes identified by RNAseq on day 13, but this requires further analysis.
Given that RNAseq detected a larger number of DEGs directly regulated by IFNT in vivo, analysis of the expression of these genes in blood may be predictive of the pregnancy status of postpartum dairy cows. It is known that IFNT acts on extrauterine tissues (8, 37), and reports from the literature have demonstrated the potential of measuring ISGs in circulation as being predictive of the pregnancy or nonpregnancy status of cattle (18, 20, 21, 48). Analysis of a panel of ISGs demonstrated to be directly regulated by short-term exposure to IFNT in vivo revealed not all ISGs are predictive of the pregnancy status of postpartum dairy cows. Nonetheless, three genes identified by RNAseq as differentially expressed on day 16 in the endometrium (ZNFX1, SAMD9, RNF213) were found to be significantly altered in the blood of postpartum dairy cows on day 18 following estrous detection. More importantly, the expression of these genes was highest and significantly different in pregnant compared with inseminated cows that were subsequently identified as nonpregnant on day 35 following ultrasound examination. However, the fact that a number of ISGs were not predictive of the pregnancy status of the animal, while surprising, may be explained by the fact that IFNT and IFNA act through the same signaling mechanisms. Indeed, incidents of bovine viral diarrhea in cattle have been shown to induce the expression of classical type 1 IFN genes in peripheral blood mononuclear cells (22). We propose therefore that measuring the expression of these genes in dairy cows on day 18 following insemination as a predictor of pregnancy success may be limited due to propensity of these genes for underlying subclinical infections, which may lead to false positives.
In summary, using RNAseq technology, we have identified DEGs in the endometrium of pregnant compared with cyclic heifers on day 16, the earliest stage of pregnancy these transcripts have yet been identified. The expression of a number of these genes (HERC6, ZNFX1, PARP12, IFI16, RNF213, and DDX58) is directly regulated by IFNT in the bovine endometrium in vivo, and expression of ZNFX1, SAMD9, and RNF213 in blood may be predictive of the pregnancy status of dairy cows as early as day 18 following insemination. RNAseq has also identified and qRT-PCR validated a number of genes that are altered between pregnant and cyclic endometria as early as day 13, and these differences are maintained to day 16. The expression of the majority of these genes is not directly regulated by IFNT in vivo and may therefore be regulated by, as yet unknown, conceptus-derived factors. However, the regulation of MGC127874 by IFNT may provide evidence of the earliest IFNT response of the endometrium to early pregnancy in cattle.
This work was supported by Science Foundation Ireland under grant numbers 06/INI/B62 and 07/SRC/B1156 (the opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Science Foundation Ireland). The authors also acknowledge Science Foundation Ireland for the provision of computational facilities and support as well as the ANR-08-GENM-037_MEETAC grant for financial support.
No conflicts of interest, financial or otherwise, are declared by the author(s).
Author contributions: N.F., N.M.-A., O.S., M.A.C., T.F., J.F.R., P.L., and A.C.O.E. conception and design of research; N.F., G.B.D., J.A.B., O.S., and T.F. performed experiments; N.F., P.A.M., J.P.M., A.K.K., and B.J.L. analyzed data; N.F. and P.A.M. interpreted results of experiments; N.F. prepared figures; N.F. drafted manuscript; N.F., G.B.D., P.A.M., J.A.B., J.P.M., A.K.K., N.M.-A., O.S., B.J.L., M.A.C., T.F., J.F.R., P.L., and A.C.O.E. edited and revised manuscript; N.F., G.B.D., P.A.M., J.A.B., J.P.M., A.K.K., N.M.-A., O.S., B.J.L., M.A.C., T.F., J.F.R., P.L., and A.C.O.E. approved final version of manuscript.
The authors acknowledge the help of all present and previous graduate students, postdoctoral scientists, and technical staff at University College Dublin and Institut National de la Recherche Agronomique (INRA) for assistance in sample collection. We thank Gilles Charpigny, Pierrette Reinaud, and Christophe Richard (INRA, Jouy-en-Josas, France) for the kind gift of recombinant ovine IFNT and generation of the in vivo IFNT model.
↵1 The online version of this article contains supplemental material.
- Copyright © 2012 the American Physiological Society