Transcriptomic Changes Associated with Pregnancy in a Marsupial, the Gray Short-Tailed Opossum Monodelphis domestica

Live birth has emerged as a reproductive strategy many times across vertebrate evolution; however, mammals account for the majority of viviparous vertebrates. Marsupials are a mammalian lineage that last shared a common ancestor with eutherians (placental mammals) over 148 million years ago. Marsupials are noted for giving birth to highly altricial young after a short gestation, and represent humans’ most distant viviparous mammalian relatives. Here we ask what insight can be gained into the evolution of viviparity in mammals specifically and vertebrates in general by analyzing the global uterine transcriptome in a marsupial. Transcriptome analyses were performed using NextGen sequencing of uterine RNA samples from the gray short-tailed opossum, Monodelphis domestica. Samples were collected from late stage pregnant, virgin, and non-pregnant experienced breeders. Three different algorithms were used to determine differential expression, and results were confirmed by quantitative PCR. Over 900 opossum gene transcripts were found to be significantly more abundant in the pregnant uterus than non-pregnant, and over 1400 less so. Most with increased abundance were genes related to metabolism, immune systems processes, and transport. This is the first study to characterize the transcriptomic differences between pregnant, non-pregnant breeders, and virgin marsupial uteruses and helps to establish a set of pregnancy-associated genes in the opossum. These observations allowed for comparative analyses of the differentially transcribed genes with other mammalian and non-mammalian viviparous species, revealing similarities in pregnancy related gene expression over 300 million years of amniote evolution.


Introduction
The transition from oviparity to viviparity has occurred independently multiple times in vertebrate evolution. Indeed, most vertebrates including cartilaginous fishes, bony fishes, amphibians, squamate reptiles, and mammals have lineages that transitioned to viviparity from oviparous ancestors [1]. In many viviparous vertebrates the transition away from oviparity involved the evolution of the placenta, which serves as an interface between the maternal system and the conceptus. The placenta is a specialized organ that facilitates the transfer of gases, nutrients, and waste as well as providing a physical barrier to protect the fetus from pathogens [2].
Based on phylogenetic relationships between metatherian (marsupial) and eutherian ("placental") mammals, viviparity appears to have evolved once early in the therian lineage [3][4][5]. There are three general categories of placental invasiveness in therians. The least invasive is epitheliochorial placentation, which is a superficial apposition of fetal and maternal membranes, where as the most is the hemochorial placenta where fetal trophoblast invades the uterine endometrium and remodels maternal capillaries [6,7]. The third major type, an intermediary endotheliochorial placenta, also invades into the endometrium but does not directly contact maternal circulation [7]. Eutherians vary from species to species in placental invasiveness [6]. Primates and rodents typically have hemochorial placentation, whereas ruminants like cows and sheep use epitheliochorial placentas [8,9]. Most branches of eutherians also have species that use endotheliochorial placentation, such as canines and bats [6,9].
Prior to the widespread use of molecular phylogenetics, most theories of placental evolution were based on assuming a progression from less invasive to more invasive [10]. The current view is that the ancestral eutherian placenta was either hemochorial or endotheliochorial [9,[11][12][13][14][15]. The less invasive epitheliochorial placenta in eutherian species is likely an evolutionary innovation to give the mother greater control over nutrient and gas exchange [14,15]. There is a correlation between forming an epitheliochorial placenta and giving birth to more precocial offspring in eutherians [14]. Garratt and colleagues [15] hypothesized that this may be due to maternal and fetal interests deviating relatively quickly after birth in many eutherian species with epitheliochorial placentation.
The marsupials are the sister group to eutherians that give birth to relatively altricial young [16]. In contrast to the eutherian placentation trend, marsupials tend to form epitheliochoriallike placentas despite their altricial offspring. Moreover most marsupials use the yolk sac, as opposed to the allantois, to contact the chorion and form the trophoblast that is apposed to maternal membranes [16,17]. A notable exception is within the family Peramelidae (bandicoots and echymiperas) where an invasive chorioallantoic placenta forms in the last few days before parturition [16,18].
The morphological changes of the placenta throughout pregnancy have been well characterized in several marsupial species such as Sminthopsis crassicaudata, Macropus eugenii, and M. domestica [19][20][21][22][23]. M. domestica forms an invasive endotheliochorial yolk sac placenta in the last 60 hours prior to birth and, characteristic for marsupials, a permeable shell coat separates the fetal and maternal membranes for the majority of gestation [21]. M. domestica is a member of the Didelphidae family, considered to be among the least derived marsupial lineage and likely represents the ancestral reproductive characteristics of marsupials [17,21]. Therefore M. domestica is potentially an important model for discovering other basal characteristics of reproduction in marsupials and perhaps even therians in general.
More recent analyses of gene expression in reproductive tissues during opossum pregnancy have made important contributions to our understanding of the evolution of the placenta in eutherian mammals [24,25]. Kin and colleagues described the distribution of important eutherian transcription factors (TFs) in M. domestica endometrium, determining the ancestral cell type of eutherian endometrial stromal fibroblasts (ESFs) that undergo decidualization [24]. Lynch and colleagues recently revealed that endometrially-expressed transposable elements (TEs), ancient in the eutherian lineage, have played an important role in the evolution of decidual tissues; a tissue type lacking in marsupials [25]. This has demonstrated key adaptations in endometrial tissue that may have enabled the successful extended gestation periods seen in eutherians.
Needed for a greater understanding of the mammalian adaptation to viviparity are analyses of those genes differentially transcribed between pregnant and non-pregnant tissues within a non-eutherian mammal. This may also facilitate the discovery of genes crucial to marsupial pregnancy that may or may not be shared with eutherians or, indeed, other viviparous vertebrates. Presented here is just such a comparison with an analysis of the pregnant and non-pregnant uterine transcriptomes of M. domestica.

Results
To investigate changes in gene expression patterns associated with opossum pregnancy, transcriptome analyses were performed on uterine tissues from late stage pregnant (P), virgin (V), and non-pregnant (N) experienced breeders. Three biological replicates were used for each condition. Paired-end reads were generated on the Illumina HiSeq 2000 platform and aligned to the M. domestica genome. The number of mapped reads per sample ranged from 17.1 to 25.9 million at the high end, with an average of 20.4 million (S1 Table). The numbers of aligned reads per sample were more than sufficient for differential expression analysis with treatment groups containing three biological replicates [26].
Comparing the global transcriptomes individually, the nine samples segregated into pregnant (P) and non-pregnant groups (N and V) based on Jensen-Shannon (JS) distances (Fig 1). The P transcriptomes grouped together and were more similar to each other than they were to any other samples. Of the non-pregnant individuals, all three N grouped together and two of the three V individuals grouped together. Virgin opossum V1 did not group with the other virgin animals, but her uterine transcriptome was more similar to the other non-pregnant animals than to the pregnant animals (Fig 1). For this reason analyses both including and excluding V1 were performed when testing the P group against the V and N groups. No other replicate in any group was dissimilar enough to warrant being separated for further analyses. Three replicates per treatment group appeared to be adequate since they grouped together well overall indicating consistency between biological replicate transcriptomes (Fig 1). The P group was tested for differentially transcribed genes compared to N, all three V, and V2 and V3 alone (V2 & V3) (Fig 2).
The three sample groups, V, N, and P, were compared in a pair-wise manner and independently assessed for differential transcription using Cuffdiff, DESeq, and edgeR. Since these algorithms assess differential expression using different parameters they identified partially, but not fully, overlapping sets of genes (Fig 2). Genes that were identified as being differentially transcribed by at least two of the three programs were used for this study (red triangles in Fig 2).
Transcripts of 2202 genes were significantly more abundant in the P uterus compared to N (Fig 2A), whereas 2910 genes were significantly less ( Fig 2B). When the P samples were compared to all three V only 986 gene transcripts increased in abundance (Fig 2C). However when V1 was excluded this number was 2061 (Fig 2E), more similar to the comparison with N ( Fig  2A). There were 1798 and 2972 decreased gene transcripts in the same comparisons, respectively (Fig 2D and 2F). Therefore, including the outlier V1 suppressed the number of differentially genes. Its exclusion made the V group more similar to N in numbers of differentially transcribed genes.

Validation of differential gene transcript abundance
To validate the Illumina data, we developed quantitative PCR (qPCR) assays for a subset of genes that were differentially represented in the different groups. These were membrane cofactor protein (CD46), MAC-inhibitory protein (CD59), interleukin-6 (IL6), lysosomal-associated membrane protein 1 (LAMP1), and lumican (LUM). The direction of log2-fold change of gene expression in P compared to N samples in the Illumina datasets was recapitulated by the qPCR (Fig 3). The same was found in a P vs V comparison. These genes were not significantly differentially expressed in V compared to the N group in either Illumina or by qPCR (S1 Fig).

Pregnancy-associated genes in the opossum
To determine gene sets that were explicitly differentially transcribed during pregnancy, the intersections of differentially transcribed genes in the P set compared to both non-pregnant sets (N and V) were identified. Such genes should represent a minimal, high confidence "pregnancy-associated" set for the opossum. In the intersections of all three comparisons 932 genes were significantly more abundant and 1482 genes were less abundant in the P group (Fig 4). These increased and decreased differentially transcribed pregnancy-associated genes are listed in S2 and S3 Tables. Gene Ontology (GO) analyses of the pregnancy defining genes revealed 18 Overlap of genes called as significant in differential expression programs. Venn diagrams of the numbers of gene transcripts found to be significantly more abundant (A, C, and E) or less abundant (B, D, and F) when comparing the dataset from the P group to N (A and B), V (C and D) and V2 & V3 only (E and F). Genes were tested for differential expression according to Cuffdiff (pink), DESeq (yellow), and edgeR (blue). Red triangles denote the cases where at least two out of three programs were in agreement.
GO terms were overrepresented and 6 GO terms were underrepresented (Table 1). In the pregnancy-associated genes with decreased transcription during pregnancy, 25 GO terms were overrepresented and two were underrepresented ( Table 2). The most overrepresented GO term in the increased in pregnant gene set was ion transport (p = 9.86E-06) ( Table 1); the most overrepresented in the decreased in pregnant gene set was ectoderm development (p = 2.23E-09) ( Table 2).

Differential transcript abundance between pre-vs. post-estrus
There is evidence that marsupial uterine differentiation is not completed until the first estrus occurs [16,27]. The data sets generated for this study allow for investigation of transcriptomic differences between the N and V groups that might be linked to this differentiation. The V and N samples were transcriptomically similar to each other when compared to the P group ( Fig  1). In order to determine the genes that might be linked to uterine differentiation, the preestrus V dataset was compared to the post-estrus N dataset. The intersections of N vs V and N vs V2 & V3 were examined and 8 genes were significantly more abundant and 4 were less abundant in V relative to N, respectively (Table 3). Gene transcription levels in P vs N according to differential expression programs compared to qPCR data. Log2-fold changes in the transcription of membrane cofactor protein (CD46), MAC-inhibitory protein (CD59), interleukin 6 (IL-6), lysosomal-associated membrane protein 1 (LAMP1), and lumican (LUM) in the pregnant group compared to the non-pregnant past breeder group. Log2-fold changes were according to Cuffdiff (pink bars), DESeq (yellow bars), edgeR (blue bars), and qPCR using the Vandesompele method (black bars). Red line indicates the threshold of log2-fold change needed for significance according to the Vandesompele method of relative quantification of qPCR data. * Adjusted p < 0.05 for differential expression programs or >2 log2-fold change for qPCR.

Conservation of genes regulated during pregnancy
Transcriptomic changes in M. domestica uterus during pregnancy were compared to previously published uterine transcriptome studies in pig (Sus scrofa), skink (Chalcides ocellatus), and seahorse (Hippocampus abdominalis) [28][29][30]. For this meta-analysis, "top-100" gene sets were identified. The top-100 significantly increased genes in opossum pregnancy are the 100 genes with the greatest transcript abundance among the 932 increased in P (Fig 4A). These were then compared to the published analyzed differentially expressed gene sets for the pig, skink, and seahorse [28][29][30]. Fifty-five of the top 100 gene transcripts most abundant during pregnancy for the opossum were identified as also increased in pregnant skink (Table 4). Only 11 and 1 of the selected top 100 opossum genes were identified as increased in pregnant pig and seahorse, respectively (Table 4). Of the 1482 genes identified as being decreased in pregnant (Fig 4B), the top 100 is that subset that had the greatest abundance of transcripts in the opossum N samples. Forty-three of these genes were also identified as significantly decreased in skink pregnancy (Table 5). Only 8 of these were identified in the pig pregnancy and none in the seahorse analyses (Table 5).

Discussion
Pregnancy involves substantial tissue remodeling and physiological changes. This study was aimed at investigating the changes that occur at the level of uterine gene transcription in a model marsupial species. Marsupials are a distinct lineage of mammals noteworthy for their short gestation times, relative immaturity at birth, and minimal placentation [16]. Viviparity likely has a single common origin in marsupials and eutherians, which last shared a common ancestor at least 148 million years (MY) ago, making marsupials an important comparative contrast to pregnancy in eutherians [3][4][5].
The comparison between pregnant and non-pregnant samples was chosen to identify the most differentially transcribed uterine genes during pregnancy. All the P and N samples used in this study came from opossums that had successfully given birth to at least one litter prior to Pregnancy-associated gene sets. Venn diagrams describing the intersection of genes that were considered differentially transcribed in comparisons of the P group to the N, V, and V2 & V3 control groups. The genes that were differentially transcribed in the P group compared to all three control groups (red circles) were considered to be differentially transcribed at a high level of confidence for this study.
this study. Based on what genes were being transcribed and the transcript number, the P transcriptomes were distinctly different from those of non-pregnant animals, whether N or V (Fig 1).
Sets of pregnancy-associated genes were identified as those at the intersection of differentially transcribed in the P versus all non-pregnant animals (N and V). Identifying the overlap in differentially transcribed genes in P vs N, P vs V, and P vs V2 & V3 generated gene sets with of high confidence for being significantly increased or decreased during opossum pregnancy. It is noteworthy that there were more genes with decreased transcription in the opossum uterus than increased transcription during pregnancy (n = 550; the difference between Fig 4A and 4B).
The top 100 most differentially transcribed opossum genes in the P and N groups were compared to results from pregnant uterine transcriptome studies in the pig, skink, and seahorse [28][29][30]. The pattern of differentially expression during pregnancy in skink had the most in common with the results in M. domestica. This similarity is unlikely due to shared structural similarities between the opossum and skink since the skink (C. ocellatus) has an epitheliochorial chorioallantoic placenta [28,31,32]. Rather the similarity in transcriptomes is likely due to the skink study focusing on a pregnancy time point more similar to that used here in the opossum. The opossum P samples were from uteruses preceding parturition by only 12-24 hours and there was likely to be late-pregnancy and parturition-specific gene transcription represented in the datasets. In contrast there was less overlap with the transcriptome in the pig. This may due to the pig study focusing primarily on early pregnancy stages and also the pig's chorioallantoic epitheliochorial placenta type and the source of embryonic nutrition, which primarily histotrophic rather than hemotrophic [29,33,34]. Likewise, there was very little overlap with the seahorse, which may be due to evolutionary distance or, again, the pregnancy time points examined.
Since there was substantial overlap in genes differentially transcribed in both the opossum and the skink, we focused on this comparison further. Keratin 18 (KRT18) is has been previously shown to be essential to maintaining pregnancy in mice [35] and KRT18 is increased during pregnancy in both opossum and skink (Table 4). KRT18 is the only keratin gene that is increased in pregnant opossum (S2 Table). This is in contrast to eutherian pregnancy, such as in pigs and humans, where in addition to KRT18, KRT8 and KRT19 are also up-regulated [29,36]. Other notable genes transcripts encoding structural molecules that are more abundant consistently in opossum and skink pregnancy are fibronectin (FN1) and fibulin (FBLN1) ( Table 4). Both of these have been shown to be important to maintaining pregnancy in eutherian species such as mice [37,38] and this conservation across a long evolutionary history is consistent with such proteins being critical to placental formation. In contrast serpin peptidase inhibitor clade member 2 (SERPINE2), an up-regulated tissue remodeling gene in human, mouse, pig, and even seahorse pregnancy [29,30,39,40], was among the most decreased gene transcripts during opossum pregnancy (Table 5). SERPINE2 is also not differentially regulated in skinks [28]. Therefore, although SERPINE2 transcription in pregnancy is shared between a teleost and eutherians, this does not appear to be a critically conserved gene. Not surprisingly, an area of broad conservation was hormonal regulation. For example, progesterone receptor (PGR) had fewer transcripts during pregnancy in both opossum and skink (Table 5), and down-regulation of PGR in endometrial cells is observed in eutherian species at implantation [41][42][43][44]. Kin and colleagues used immunohistochemistry (IHC) to demonstrate that term pregnant M. domestica uterine epithelium and uterine glands express PGR weakly in comparison to non-pregnant uterine tissues [24]. The transcriptome results described here are consistent with that finding (S4 Table). Oxytocin receptor (OXTR) is more abundant in opossum at terminal pregnancy and in human pregnancy OXTR expression increases over time peaking prior to parturition and playing a role in inducing contractions [45][46][47][48]. Interestingly, Yamashita and Kitano [49] hypothesized that the evolution of OXTR associated with labor contractions in eutherians occurred after the marsupial-eutherian split. However, OXTR was a top 100 most abundant gene in pregnant opossum raising the possibility of its role in parturition being more ancient than the marsupial-eutherian split.
Prostaglandins have been shown to be important to normal birth mechanisms in other marsupials such as the tammar wallaby [50,51]. Injection of exogenous prostaglandin F2 (PGF2) have also been shown to induce birthing in marsupials [50,52]. The pregnant samples used in this study were taken from opossums that were within 24 hours of giving birth and it is probable that gene transcript values in these pregnant transcriptomes were influenced by birth mechanisms at work. Genes involved in prostaglandin synthesis and detection such as prostaglandin E synthase (PTGES), prostaglandin-endoperoxide synthase 2 (PTGS2), and prostaglandin E receptor 4 subtype EP4 (PTGER4) were all included in the increased in pregnancy-associated set (S2 Table).
Eight solute carrier proteins (SLC39A9, SLC41A2, SLC16A5, SLC16A1, SLCO2A1, SLC39A8, SLC34A2L, SLC2A1) are included in the top 100 most abundant genes in opossum during pregnancy (Table 4). This is consistent with ion transport being the most highly overrepresented GO term in the datasets. Two of these, SLC39A9 and SLC41A2 are zinc and magnesium ion  transporters that are also up-regulated during skink pregnancy [28,53,54]. Transporters of cholesterol and chylomicrons like ATP-binding cassette sub-family A member 2 (ABCA2) and apolipoprotein genes APOA1, APOA4, and APOB are also in the top 100 most abundant gene transcripts in opossum pregnancy (Table 4) which is consistent with what is known about cholesterol transport across eutherian placentas [55]. Two different non-pregnant control groups were used in this study due to concern over previous observations on sexual maturation in M. domestica [27]. Briefly, M. domestica females do not appear to reach reproductive maturity until they enter estrus and mate for the first time, a   process dependent on male pheromones. Therefore, both virgin (V) females (pre-first estrus) as well as non-pregnant (N) experienced past breeders (post-first estrus) were used. In addition to providing controls for the pregnant state, this also allowed for a direct comparison between presumably immature pre-estrus and more mature post-estrus uterine tissues. In the analyses comparing the V group to P and N, only 12 genes were uniquely differentially transcribed in common (Table 3). Interestingly three of the four genes less abundant (Ig J chain, Ig kappa chain, IgG heavy chain constant domain) in virgin uterine tissue were genes associated with antibody expression (Table 3). This may indicate a need for humoral immunity in uterine tissues only after copulation and could point to a mechanism for controlling sexually transmitted infections in opossums. The morphological changes described by Stonerook and Harder [27], which is primarily organ weight, do not appear to be substantially reflected in changes in gene regulation in the uterus. Actin alpha cardiac muscle 1 (ACTC1) and keratin 13 (KRT13) are two of eight significantly more abundant genes in the virgin samples but no other major structural genes are differentially transcribed between the non-pregnant groups ( Table 3). The V group had the most variation between biological replicates than any other group as illustrated by V1 being an outlier (Fig 1). Nonetheless, V1 was transcriptomically more similar to the non-pregnant samples than to the pregnant ones. V1's outlier status may have been due to her being 4.8 months and old at the time of collection, whereas the other two virgin animals (V2 and V3) were 10 months old. M. domestica individuals are generally considered to be sexually mature at 5-6 months old [56], but V1's uterus was likely in a different stage of development from the older virgin animals. Aside from V1, all biological replicates were more similar to each other than to samples from other treatment groups (Fig 1). Non-pregnant individuals (groups N and V) were more similar to each other than they were to the pregnant individuals (Fig 1). This was an expected result and provided justification for retaining V1 in the analyses of differential transcription. There were more differentially transcribed genes shared between the P vs N and P vs V2 & V3 comparisons than between the P vs V and P vs V2 & V3 comparisons (Fig 4). This confirmed that V1 was skewing the differential transcription results in the P vs V comparison.
Kin and colleagues identified TFs spatially in M. domestica endometrium tissue sections in non-pregnant, pre-implantation pregnancy, and term pregnancy using IHC [24]. They concluded that the stromal mesenchymal cells seen in M. domestica endometrium are homologous to eutherian ESFs based on the expression of TFs and certain cytoskeletal proteins, and no marsupial has definitive decidualization. Two of the transcription factors (TFs) examined by Kin and colleagues were found in our pregnancy-associated gene sets. These were PGR, which was less abundant in pregnant samples, and CEBPB which was more abundant in pregnant samples (S4 Table). All samples examined in this study did have gene transcripts from the TFs evaluated in the Kin et al study (S4 Table). This is not surprising since Kin and colleagues observed transcription factors in a spatial manner by comparing IHC-stained endometrium of pregnant and non-pregnant M. domestica, whereas our study used whole sections of tissue without separating by cell type. Indeed, a limitation of the data sets used in this study is that they contain multiple cell and tissue types and therefore a specific transcript cannot be ascribed to a particular cell or tissue. Therefore in future studies on the genes identified as pregnancy-associated here, a spatial analysis of gene expression would be informative.
In conclusion, the analysis of the opossum late pregnancy uterine transcriptome enabled the discovery of a set of genes that can be described as being pregnancy-associated in the opossum. Curiously, morphological changes associated with sexual maturation in the opossum uterus reported previously are not obviously reflected in the transcriptome of this species. These analyses also revealed significant conservation of gene regulation among distantly related species. Reptiles and mammals are separated by over 300 MY of evolution and the opossum and skink do not share a common viviparous ancestor [57,58]. While the similar nature of the gene regulation is likely indicative of convergence on common mechanisms in the evolution of live birth, the transcriptomic similarities in these studies could be explained by the similar pregnancy stages examined here and in the skink.

Animals and tissue collection
This study was approved under protocol numbers 13-100920-MCC and 15-200334-B-MC from the University of New Mexico Institutional Animal Care and Use Committee. Animals were from a captive-bred research colony housed at the University of New Mexico Department of Biology Animal Research Facility. Founders were obtained from M. domestica Populations 1 and 2 at the Southwest Foundation for Biomedical Research, San Antonio, TX [59]. All animals were fed ad libitum a diet of Labdiet Short Tailed Opossum #2 pellets (90%, Lab Supply TX) and dried mealworms (10%, Lab Supply TX). They were given drinking water ad libitum in sterile glass sipper bottles. All animals were housed individually in standard rat caging and bedding. When not breeding adult females and males were caged separately. All animals used in tissue collection were euthanized by inhaled isoflurane overdose (>5%) until breathing stopped for one minute.
Uterine tissues were collected from virgin (V), pregnant (P), and non-pregnant (N) experienced breeders. Non-pregnant past breeders had at least one successful pregnancy prior to harvesting. Uterine horns were excised and separated from surrounding tissues. In the case of pregnant uteri, the uterine horns were opened laterally and the intact embryos and amnions removed. Endometrium was separated from myometrium by teasing apart the tissue layers in shallow petri dishes filled with RNALater buffer (Ambion). In pregnant samples invasive fetal chorion and yolk sac membranes could not be fully separated from maternal endometrium. RNA extractions were either performed immediately or tissues were stored in RNALater overnight at 4°C. For longer-term storage, excess RNALater was removed and tissues were stored at -80°C until use.

RNA extraction
Total RNA was isolated by homogenizing tissue in TRIzol (Ambion) using a sterile glass homogenizer followed by chloroform extraction and centrifugation at 4°C for 15 min. Protein and DNA contamination was removed using the PureLink RNA Mini (Ambion) and On-Column DNase Treatment (Invitrogen) kits following manufacturers' recommended protocols. RNA was stored in RNase-free water at -80°C. RNA quality was assessed using a 2100 Bioanalyzer (Agilent) and concentration was determined using a Qubit 2.0 Fluorometer (Life Technologies).

cDNA library synthesis and Illumina sequencing
All cDNA library synthesis for high-throughput RNA-seq and Illumina sequencing procedures were performed at the National Center for Genome Resources (NCGR) in Santa Fe, New Mexico, USA. A Sciclone G3 Automated Liquid Handling Workstation (Caliper Life Sciences) with a Multi TEC Control for heating and cooling steps (INHECO) was used for the majority of liquid handling. Poly-A RNA was isolated using magnetic RNA Purification Beads (Dynabeads: Invitrogen). RNA was fragmented by high temperature (95°C for 8 min). SuperScript II enzyme (Intvitrogen), First Strand Master Mix (Illumina), and random hexamer primers (Illumina) were added to RNA to conduct first strand cDNA synthesis by reverse transcription. Second Strand Master Mix (Illumina) was used to conduct second strand cDNA synthesis, resulting in double-stranded cDNA. Double-strand cDNA was isolated using Agencourt AMPure XP beads (Beckman Coulter). Blunt end cDNA was prepared by removing 3' overhangs and filling in 5' overhangs using End Repair Mix (Illumina). 3' ends were then adendylated using A-Tailing Mix (Illumina). Universal and barcoded TruSeq Adapters (Illumina) were ligated to cDNA ends. Polymerase Chain Reaction (PCR) was conducted with KAPA HiFi HotStart ReadyMix (KAPA Biosystems) and PCR Primer Cocktail (Illumina) to selectively amplify the cDNA fragments that had TruSeq adapters ligated to both ends. PCR products were purified using Agencourt AMPure XP beads. PCR product quality was assessed using a 2100 Bioanalyzer and cDNA quantity was measured by Nanodrop ND-1000 (Thermo Scientific). Samples were normalized to 10nM equimolar concentrations and pooled prior to flow cell injection. All Illumina sequencing was performed on an Illumina HiSeq 2000 instrument (Illumina) to generate 50bp paired-end reads. Read data was de-multiplexed (segregated by library index) using Illumina CASAVA v1.8.2 to produce FASTQ files. Thereafter, NCGR's custom contaminant filtering pipeline removed anomalous sequences (i.e. Illumina PhiX control, library adapters, primer dimmers, and library indexes not part of the experiment). Quality assessment of read data was performed using FASTQC.

Data access
All high-throughput sequence data sets generated for this study have been deposited at NCBI Sequence Read Archive (SRA) under the accession numbers SRR2969483, SRR2969536, SRR2970443, SRR2972728, SRR2972729, SRR2972792, SRR2972837, SRR2972840, and SRR2972848.
Illumina read quality control, alignment, and quantification All alignment and differential expression bioinformatics programs were run through the online Lumenogix platform (api.lumenogix.com, [60]). Reads were aligned to the annotated M. domestica genome (assembly: monDom5, annotation: Ensembl release 76) using Tophat2 [61,62]. For a summary of Illumina read counts for this study see S1 Table. Read-count based differential expression was performed by Cufflinks, which normalizes by transcript length and apportions multi-mapped reads, and HTSeq which counts at the gene level and counts only uniquely mapped reads. Differential expression analysis was performed in three different algorithms: Cuffdiff, DESeq, and edgeR. Filtering was performed using log2-fold changes and adjusted p-value 0.05.

Quantitative PCR
All RNA used for qPCR originated from the same RNA samples used for Illumina sequencing. Prior to cDNA synthesis, all RNA used in cDNA synthesis was treated with DNase using the TURBO DNA-free Kit (Ambion) according to manufacturers' recommended protocols. 500ng of total RNA was used to make cDNA libraries by reverse transcriptase PCR (RT-PCR) using the SuperScript III First Strand Synthesis kit (Invitrogen) according to the manufacturer's instructions. All cDNA libraries were made in triple replicates and then pooled to reduce bias generated during reverse transcription. Transcription levels of specific genes were assessed by qPCR using ABsolute Blue QPCR SYBR Green ROX Mix (Thermo Scientific) according to manufacturer's instructions. Primers used in qPCR were designed for the M. domestica genome and according to manufacturer's recommendations for qPCR primer properties. Primer sequences and properties are in S2 Table. All qPCR reactions were performed in triplicate repeats on an ABI Prism 7000 Sequence Detection System (Applied Biosystems) using the default cycling parameters with a dissociation step. qPCR data were assessed by the Vandesompele method [63,64] using tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein zeta polypeptide (YWHAZ) and TATA box-binding protein (TBP) as reference genes. These reference genes were chosen based on relatively consistent transcription across Illumina data samples (S5 Table), as well as literature recommending them as reference genes for placental gene expression studies [65,66]. All pregnant and virgin qPCR gene transcription levels described here were relative to the non-pregnant past breeder group, and greater than 2-fold change was considered significant. Genes were chosen for qPCR confirmation based on being differentially transcribed in the Illumina data set, having >20 reads per sample in the Illumina data set, and novel qPCR primers passing efficiency testing. All primers were generated using the online tool Primer 3 (http://primer3.ut. ee/, [67]) with M. domestica gene sequences from the MonDom5 genome assembly. Primer sets used for qPCR and their associated properties are listed in S5 Table. Differential transcription analysis Three separate differential expression analysis tools, Cuffdiff version 2.2.1, [68,69], DESeq version 1.2.1 [70], and edgeR version 2.0.5 [71], were used to assess relative gene transcription between the P, N and V groups. All Venn diagrams were created using the R package VennDiagram [72]. The R package cummeRbund [73] was used to produce a dendrogram describing Jensen-Shannon (JS) distances between sample transcriptomes. The dendrogram was modified for readability in Microsoft PowerPoint. Since DESeq and edgeR use read counts and Cuffdiff uses Fragments Per Kilobase of transcript per Million mapped reads (FPKM) to calculate differential expression, the log2-fold change and associated adjusted p-value of genes were reported in cases where differential expression analysis results of specific genes were compared. The adjusted p-value threshold was 0.05 for all used differential expression programs.

Gene Ontology analysis
Gene Ontology (GO) analysis was performed using online PANTHER tools (http://www. pantherdb.org/ [74,75]. Specific gene sets were tested for statistical overrepresentation or underrepresentation of PANTHER GO-Slim Biological Process terms. Gene sets were compared to the whole uterine transcriptome gene set (all genes with expression of 1 FPKM in at least one uterine sample transcriptome) and GO terms with p < 0.05 (Bonferroni corrected for multiple testing) were considered significant. The "Unclassified" GO category was ignored in all analyses.
Supporting Information S1 Fig. Gene transcription levels in P vs V according to differential expression programs compared to qPCR data. Log2-fold changes in the expression of membrane cofactor protein (CD46), MAC-inhibitory protein (CD59), interleukin 6 (IL-6), lysosomal-associated membrane protein 1 (LAMP1), and lumican (LUM) in the virgin group compared to the non-pregnant past breeder group. Log2-fold changes were according to Cuffdiff (pink bars), DESeq (yellow bars), edgeR (blue bars), and qPCR using the Vandesompele method (black bars). Red line indicates the threshold of log2-fold change needed for significance according to the Vandesompele method of relative quantification of qPCR data. None significant. (TIF) S1