Analysis of Stage-Specific Gene Expression Profiles in the Uterine Endometrium during Pregnancy in Pigs

The uterine endometrium plays a critical role in regulating the estrous cycle and the establishment and maintenance of pregnancy in mammalian species. Many studies have investigated the expression and function of genes in the uterine endometrium, but the global expression pattern of genes and relationships among genes differentially expressed in the uterine endometrium during gestation in pigs remain unclear. Thus, this study investigated global gene expression profiles using microarray in pigs. Diverse transcriptome analyses including clustering, network, and differentially expressed gene (DEG) analyses were performed to detect endometrial gene expression changes during the different gestation stages. In total, 6,991 genes were found to be differentially expressed by comparing genes expressed on day (D) 12 of pregnancy with those on D15, D30, D60, D90 and D114 of pregnancy, and clustering analysis of detected DEGs distinguished 8 clusters. Furthermore, several pregnancy-related hub genes such as ALPPL2, RANBP17, NF1B, SPP1, and CST6 were discovered through network analysis. Finally, detected hub genes were technically validated by quantitative RT-PCR. These results suggest the complex network characteristics involved in uterine endometrial gene expression during pregnancy and indicate that diverse patterns of stage-specific gene expression and network connections may play a critical role in endometrial remodeling and in placental and fetal development to establish and maintenance of pregnancy in pigs.


Introduction
The uterus plays a critical role in the control of the estrous cycle and pregnancy in pigs. During pregnancy, the uterus communicates with the conceptus (embryo/fetus and associated extraembryonic membranes) to establish and maintain pregnancy and it undergoes dramatic functional and morphological changes [1]. Synchronization of the developing embryo and appropriate endometrial remodeling is essential to a successful pregnancy, while failure leads to embryonic mortality. Accordingly, knowledge regarding the pattern of expression of uterine endometrial genes and their function in endometrial remodeling, embryo development and placentation during pregnancy is an important aspect of management of successful pregnancies. Thus, numerous studies have focused on elucidating the expression and function of uterine genes during pregnancy.
In pigs, maternal recognition of pregnancy, the extension of the functional lifespan of the corpora lutea, occurs on day (D) 12 of pregnancy, followed by embryo implantation, which continues until D18 of pregnancy [2,3]. Once implantation is accomplished, the uterus participates in the formation of a maternal component of a true epitheliochorial type placenta during mid-to-late pregnancy, which transports nutrients to the developing embryos and exchanges gases. During this period placentation is completed, and rapid organogenesis of the developing fetuses occurs [4,5]. In the late stage of pregnancy near term, the uterus experiences drastic changes in the expression of many genes and gene product function resulting from altered placental hormone secretion in order to prepare for parturition. For example, progesterone levels decrease, whereas estrogen levels increase [6]. These hormonal changes promote uterine contractility through their effects on myometrial contractile proteins, gap junction formation and increasing responsiveness of the uterus to oxytocin and prostaglandin (PG) F 2α [7].
To investigate the pattern of uterine endometrial gene expression and gene function during pregnancy, a one-by-one approach was used before the genomic era. Although this approach has significantly helped us to understand gene expression and function in the uterus, a genome-wide approach using a microarray technique allows us to more efficiently investigate global gene expression in the uterus during different stages of pregnancy or in different physiological or pathological conditions. Many studies have applied microarray analysis to investigate expression of uterine endometrial genes during various stages of development or pregnancy and under pathological conditions. For example, gene expression profiles in the uterine endometrium during the implantation period have been analyzed using a genome-wide microarray technique in humans, mice and cows [15][16][17]. In pigs, microarray-based experiments have also been carried out to analyze differentially expressed genes in the endometrium during the implantation stage due to early exposure to estrogen [18] or during the estrous cycle [19]. However, there have been no previous reports on the global patterns of expression of genes and the relationships among genes affecting expression and function in the uterine endometrium during pregnancy. Therefore, in this study we analyzed global gene expression profiles using a microarray technique and performed three types of analyses: (1) detection of differentially expressed genes (DEGs) by comparing genes expressed in the uterine endometrium on D12 of pregnancy with those on D15, D30, D60, D90 and D114 of pregnancy; (2) clustering analysis to group those DEGs based on expression patterns during pregnancy; and (3) network analysis to find hub genes that are correlated with expression of other genes. From these analyses, we sought to identify novel uterine marker genes at different pregnancy stages, and to provide valuable insight into the molecular mechanisms regulating dynamic expression of uterine endometrial genes during pregnancy.

Animals and tissue preparation
All experimental procedures involving animals were conducted in accordance with the Guide for Care and Use of Research Animals in Teaching and Research and were approved by the Institutional Animal Care and Use Committee of Yonsei University (Approval No. YWC-P120). Eighteen Landrace and Yorkshire crossbred sexually mature gilts of similar age (8-10 months) were artificially inseminated at the onset of estrus (D0) and 24 h later with fresh boar semen and assigned to D12, D15, D30, D60, D90, or D114 pregnancy groups (n = 3 pigs per day). The reproductive tracts of gilts were obtained immediately after they were slaughtered at a local slaughterhouse and uterine endometrial tissues were harvested immediately. Pregnancy was confirmed by the presence of apparently normal filamentous conceptuses in uterine flushings on D12 and D15 and the presence of embryos and placenta at later dates. Uterine flushings were obtained by introducing and recovering 50 ml of phosphate buffered saline (PBS) (pH 7.4) after hysterectomy (25 ml/uterine horn). Endometrial tissues were dissected from the myometrium and placental tissues, collected from four different areas of the middle portion of the uterine horn, and mixed for analysis to reduce the heterogeneity of gene expression. Endometrial tissues on D30 to D114 were collected from where placental tissues were removed. Endometrial tissues were snap-frozen in liquid nitrogen and stored at -80°C for RNA extraction.

Gene expression profiling
Total RNA extraction. Total RNA was extracted from endometrial tissues using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's recommendations. The quality of total RNA was evaluated using NanoDrop (NanoDrop Technologies, Wilmington, DE) and Experion (BioRad, Hercules, CA).
Target preparation, microarray hybridization and scanning. The GeneChip Porcine Genome Array (Affymetrix, Santa Clara, CA) used in this study contained 24,123 probe sets outlining 19,675 transcripts, which represented 11,265 genes [20]. Five micrograms of total RNA from porcine endometria were used for labeling. Probe synthesis from total RNA samples, hybridization, detection, and scanning were performed according to standard protocols from Affymetrix at Seoulin Bioscience Molecular Biology Center (Seoul, Korea). Briefly, cDNA was synthesized from total RNA using the One-Cycle cDNA Synthesis Kit (Affymetrix). Single stranded cDNA was synthesized using Superscript II reverse transcriptase and T7-oligo (dT) primers at 42°C for 1 h. Double stranded cDNA was obtained by using DNA ligase, DNA polymerase I and RNase H at 16°C for 2 h, and followed by T4 DNA polymerase at 16°C for 5 min for gap filling. After clean up with the Sample Cleanup Module (Affymetrix), double-strand cDNA was used for in vitro transcription (IVT). cDNA was transcribed using the GeneChip IVT Labeling Kit (Affymetrix) in the presence of biotin-labeled CTP and UTP. After clean up with the Sample Cleanup Module (Affymetrix), 10-15 μg of labeled cRNA was fragmented from 35 to 200 bp by fragmentation buffer (Affymetrix). Fragmented cRNA was hybridized to the porcine genome microarray chips (Affymetrix) at 45°C for 16 h according to the Affymetrix standard protocol. After hybridization, the arrays were washed in a GeneChip Fluidics Station 450 with a non-stringent wash buffer at 25°C followed by a stringent wash buffer at 50°C. After washing, the arrays were stained with a streptavidin-phycoerythrin complex. After staining, intensities were determined with the GeneChip scanner 3000 (Affymetrix) controlled by Gene Chip Operating System (GCOS) Affymetrix software.

Data analysis
Differentially expressed genes (DEGs). The quality of the array image was assessed as described in the Affymetrix GeneChip expression analysis manual. A robust multi-array averaging procedure (RMA) was implemented in the R statistical package [21]. Expression values were computed from the raw CEL files by applying the RMA model of probe-specific correction for perfect-match probes. The corrected probe values were then normalized via quantile normalization, and a median polish was applied to compute one expression measure from all probe values using the RMA package. Whole expressions were log 2 -transformed after normalized values were calculated by RMA and quantile normalization. The expression of individual genes on D12 to D114 was compared using the linear models for microarray data (LIMMA) analysis [22]. The Benjamini-Hochberg correction for false discovery rate (FDR) [23] was used for all probe-level normalized data. We defined genes as differentially expressed only if they met the criteria of FDR adjusted P-value < 0.05 in the unpaired Welch t-test. In addition, the DEGs were classified as either up-or down-regulated genes depending on fold change (+ or -, respectively) by calculating log 2 (gene expression level of later stage/gene expression level of earlier stage).
Clustering analysis. Soft clustering data were obtained using the Mfuzz package implemented in R [24,25]. The raw ratios for the time profiles of DEG were log 10 transformed and then normalized such that, for each profile, the mean was zero and the standard deviation was one. The transformed profiles were then clustered using the Mfuzz package. We used the fuzzy c-means (FCM) clustering algorithm, which is a part of the package. FCM clustering is a soft partitioning clustering method that requires two main parameters (c = number of clusters, m = fuzzification parameter) and uses Euclidean distance as the distance metric. FCM assigns to each profile a membership value in the range (0, 1) for each c cluster. The algorithm iteratively assigns the profile to the cluster with the nearest cluster center while minimizing an objective function. Parameter m plays an important role in deriving robust clusters that are not greatly influenced by noise and random artifacts in data. For our analysis, the optimal values of c and m were derived by the iterative refinement procedure as previously described [26]. The final clustering was done with parameters c = 8 and m = 1.25.
Functional annotation clustering analysis using the DAVID tool. For investigating gene-set level patterns, we employed the DAVID functional annotation tool [27]. Annotated porcine gene information for DEG probe identifications was obtained from NetAffx (http:// www.affymetrix.com/analysis/index.affx). Because of the limited number of annotated genes for the porcine genome in NetAffx resulting from the limited availability of the full-length porcine cDNA sequence, we used human gene symbols annotated for the Affymetrix porcine genome microarray probe identifications for the genes whose annotations were not available in NetAffx, as described in Tsai et al. [20]. Functional annotation clustering analysis using porcine and human gene symbols in the DAVID program was conducted at the highest stringency to identify the biological function of DEGs in the uterine endometrium during pregnancy. Annotation clusters having a low enrichment score (< 1.3) were filtered out because most of the included biological terms were not significant (P-value < 0.05).
Weighted gene co-expression network analysis (WGCNA) and module construction. Module detection and characterization were performed using customized R software functions. The absolute value of the Pearson correlation between expression profiles of all pairs of genes was determined. Then, the Pearson correlation measure was transformed into a connection strength measure by using a power function Edge ij ¼ jr 2 ij j b [28]. The connectivity measure for each gene is the sum of the connection strengths (correlation β) between the gene and all the other genes in the network. Gene expression networks, like virtually all types of biological networks, exhibit an approximate scale-free topology. The coefficient of determination between log p(k) and log(k) was used to determine how well a resulting network fit the scale-free topology for a range of values. The scale-free topology criterion [28] was used to determine the power. To group genes with coherent expression profiles into modules, we used average linkage hierarchical clustering, which uses the topological overlap measure as dissimilarity. The topological overlap of two nodes reflects their similarity in terms of the commonality of the nodes they connect to each other [29,30]. A height cutoff value of 0.94 was chosen and obtained 5 gene co-expression modules in this study.

Quantitative real-time RT-PCR
To technically validate the detected results from microarray analysis such as clustering, network, and DEG finding analysis in the uterine endometrium, real-time RT-PCR was performed using the Applied Biosystems GeneAmp 7300 Sequence Detection System (Applied Biosystems, Foster City, CA) using the SYBR Green method. Complementary DNA was synthesized from 4 μg total RNA isolated from different uterine endometrial tissues, and newly synthesized cDNAs (total volume of 21 μl) were diluted 1:4 with nuclease-free water and then used for PCR. The Power SYBR Green PCR Master Mix (Applied Biosystems) was used for PCR reactions. The final reaction volume of 20 μl included 2 μl of cDNA, 10 μl of 2X Master mix, 2 μl of each primer (100 nM), and 4 μl of dH 2 O. PCR cycle parameters were 50°C for 2 min and 95°C for 15 min followed by 40 cycles of 95°C for 30 s, 60°C for 30 s and 72°C for 30 s. Sequences of primer pairs for PCR reactions and product sizes are listed in Table 1. The results were reported as the relative expression to the level on D12 of pregnancy after normalization of the transcript amount to the two endogenous control genes, porcine ribosomal protein L7 (RPL7) and ubiquitin B (UBB), by the 2 -ΔΔCT method [31]. We performed the Pearson correlation test in order to technically validate the WCGNA results.

Identification of DEGs in different stages of pregnancy in the uterine endometrium
We analyzed gene expression profiles from microarray data consisting of 18 endometrial samples (three per each of the 6 selected days of pregnancy) to detect DEGs during pregnancy. For accurate estimation of DEGs, we performed probe and scale normalization using RMA and quantile normalization, respectively. Then, DEGs were detected by comparing genes expressed in the endometrium on D12 of pregnancy with those on D15, D30, D60, D90 and D114 of pregnancy using the LIMMA R package. As a result, we identified DEGs (FDR adjusted P < 0.05) as shown in Table 2. Interestingly, we observed that the number of DEGs during the late stage of pregnancy tended to increase when genes expressed on D12 of pregnancy were compared with genes expressed on D15, D30, D60, D90, and D114 of pregnancy: 256, 449, 1,620, 1,467 and 2,146 genes were significantly up-regulated, and 248, 636, 1,892, 1,947, and 1,884 genes were significantly down-regulated, respectively. Furthermore, we visualized the expression pattern of DEGs during pregnancy, as shown in Fig 1A and 1B, and rearranged all of the up-or down-regulated DEGs by day of pregnancy and their identity ( Table 2). Interestingly, many DEGs did not overlap with different stage groups, which suggests that detected DEGs are very stage-specific. In DEG analysis, we observed a total of 6,991 genes without overlapped genes in different stage groups. More detailed information regarding these DEGs is listed in S1 and S2 Tables.
Clustering analysis to determine the global expression pattern of genes during pregnancy Next, to distinguish the detected DEGs defined above in the endometrium throughout all stages of pregnancy, we conducted clustering analysis using FCM clustering implemented in the mFuzz R package. From the analysis, the optimal number of clusters with a similar expression pattern during pregnancy was calculated to be 8 (Fig 2). Each cluster included 523, 1,059, 831, 820, 1202, 737, 926 and 893 DEGs, respectively. More detailed clustering results with DEGs are shown in S3 Table. In cluster 1, expression of genes increased over the gestation  -ACT GGC TGT CTG TCT TGC AGT GAT T-3'  104  60°C   R: 5'-TGC TAT CAC GAA CCC AGT TAA GGT G-3'   CCL28  NM_001024695  F:5'-GAT GTG CCC CTT TAC TGT TCC TCT T-3  Functional annotation clustering analysis of the DEGs in each cluster to determine biological function during pregnancy Having determined the global gene expression pattern of DEGs in the uterine endometrium during pregnancy and having identified groups of genes showing a similar expression pattern, we next sought to determine the function of genes in each cluster. To answer this question, we used DAVID functional annotation clustering analysis with annotated gene symbols. We used the highest stringency cut-off criteria for reducing false positives, and obtained enriched functional annotation groups composed of significant gene ontology (GO) terms or biological features in 8 clusters (Table 3 and S4 Table).
In cluster 1 where gene expression increased over gestation and then declined at term, annotation groups, blood vessel development/morphogenesis, ion/cellular chemical homeostasis, and regulation-related terms were observed. In cluster 2 where gene expression decreased over gestation from early to term pregnancy, annotation groups included negative regulation of cellular biosynthetic/cellular metabolic processes and cellular protein catabolic /proteolysis processes. In cluster 5 where gene expression decreased on D60 to D90 of pregnancy, annotation groups included RNA processing/splicing, regulation of gene expression/macromolecule biosynthetic processes and cellular glucan/polysaccharide/energy reserve metabolic processes. In cluster 8 where gene expression was highest on D12 of pregnancy and remained low thereafter, annotation groups included blood vessel development/differentiation and neuronal projection development/cell morphogenesis.
Weighted gene co-expression network analysis (WGCNA) to determine groups of closely interrelated genes Next, further efforts were made to determine the interrelationship among genes expressed in the uterine endometrium during pregnancy. We applied WGCNA and identified groups (modules) of closely interrelated genes that had similar patterns of connection strengths to other genes, or high topological overlap [32]. Module detection and characterization were performed using customized R software functions, and gene co-expression modules were identified by   average linkage hierarchical clustering with a topological overlap matrix (TOM). As a result, 5 modules in different colors corresponding to branches of the dendrogram as visualized in the TOM plot were obtained by cutting the tree of the dendrogram at the height of 0.94 (Fig 3A). For each module, a heatmap was produced with rows corresponding to genes and columns corresponding to samples ordered by pregnancy stage (Fig 3B). Modules in yellow, turquoise, blue, green, and brown colors included 61, 816, 506, 39, and 246 genes, respectively. A list of genes in each module is provided in S5 Table. Interestingly, heatmaps for 5 modules showed a common pattern of gene expression during pregnancy; genes up-regulated on D12 and D15 of pregnancy were down-regulated after D30 of pregnancy, and genes down-regulated on D12 and D15 of pregnancy were up-regulated after D30 of pregnancy. We evaluated the intramodular connectivity of genes to each other to identify hub genes that could be markers for specific physiological functions in the uterus. In this network analysis study, hub genes were defined as highly connected genes (> 10 connected nodes). Twenty-five genes that had a connectivity of absolute scaled K value greater than 0.9 in each module were selected as hub genes (Table 4). In the yellow module, SOX13 and MRPL2 showed the highest connectivity to other genes positively and negatively, respectively. HNRNPA2B1 and REXO2 showed the highest connectivity to other genes positively and negatively, respectively, in the turquoise module. In the blue module, ALPPL2, RANBP17 and NF1B had the highest connectivity to other genes positively, and SPP1 and CST6 had the highest connectivity to other genes negatively. In the green module, GINS1, CDC2 and KPNA2 were highly connected to genes positively, and NOLA2 had the highest connectivity to other genes negatively. MYD88 showed the highest connectivity to other genes positively, and ATP1A1 showed the highest connectivity to other genes negatively in the brown module.

Real-time RT-PCR analysis to validate detected hub genes and correlated genes
To demonstrate the correlation of hub genes determined in WGCNA and confirm clustering analysis of DEGs, we conducted real-time RT-PCR for selected genes. We chose RANBP17 and SPP1 in the blue module in WGCNA because the blue module included a large number of genes, and RANBP17 and SPP1 showed the highest absolute connectivity to other genes in WGCNA. In addition, some genes correlated with hub genes were additionally selected for real-time RT-PCR. Pearson correlation was calculated between hub genes RANBP17 and SPP1 and the other genes, as shown in S6 and S7 Tables. Genes positively correlated with RANBP17 included SAL1, CCL28, GALT and NF1B3, and negatively correlated genes included SPP1, SLPI, MPZL2, UFBP and CST6 (S6 Table). Genes positively correlated with SPP1 included SLPI, MPZL2, UFBP, CST6 and SUCLA2, and negatively correlated genes included RANBP17, SAL1 and CCL28 (S7 Table).
As shown in Fig 4, a heatmap and correlation-based network plot revealed the relationships among the hub genes RANBP17 and SPP1 and their correlated genes SAL1, NF1B, SLPI, MPZL2, UFBP, CST6, CCL28, and SUCLA2. The general expression pattern of these genes in real-time RT-PCR matched the expression pattern of genes in each cluster where those genes were grouped, confirming the accuracy of the microarray analysis. In Fig 4A, expression of most validated genes dramatically changed after D15 of pregnancy except SAL1 and CCL28, the expression of which changed after D12 of pregnancy. We also performed network analysis in real-time RT-PCR in order to identify their relationship. The correlation-based network plot showed these relationships and the Pearson correlation test was used to determine significance. As shown in Fig 4B, significant correlation (FDR adjusted P-value < 0.05) and their direction (positive or negative) are visualized as a network plot. As a result, all genes were significantly correlated each other. As expected, RANBP17 and SPP1 were observed as hub genes in realtime RT-PCR. NFIB, SAL1, and CCL28 were highly positively correlated with RNABP17 and SPP1, UABP2, SLP1, MPZL2, and CTS6 were highly negatively correlated with RNABP17. These results confirmed the hub genes detected in WGCNA.

Discussion
The uterine endometrium plays an essential role in the establishment and maintenance of pregnancy and undergoes morphological and functional changes during pregnancy in mammals. To support pregnancy, expression of uterine endometrial genes alters dynamically depending  on the stage of gestation and is affected by hormones and cytokines of ovarian, conceptus, chorionic and other origins. Although many studies have investigated the expression and function of endometrial genes in pigs, the genes expressed in the uterine endometrium and the cellular and molecular function of those gene products during pregnancy are not completely understood. In this study, we analyzed the global expression patterns of genes in the endometrium from the early to the late stage of pregnancy using a microarray technique in pigs. As a result, we determined the following: 1) genes are expressed differentially during the different stages of pregnancy; 2) groups of genes show similar expression patterns during pregnancy. We also evaluated the functional annotation of grouped genes and the correlation of genes expressed in the endometrium during pregnancy and found that hub genes closely correlated with other genes. To our knowledge, this is the first report investigating the global expression patterns of uterine endometrial genes throughout gestation in pigs. Expression of endometrial genes is regulated by many factors and changes dynamically during pregnancy. In this study we identified DEGs by comparing expression levels of endometrial genes on D15, D30, D60, D90 and D114 of pregnancy with those on D12 of pregnancy, which is when maternal recognition of pregnancy occurs and embryo implantation begins in pigs [1]. Using clustering analysis we grouped those DEGs into 8 clusters based on the overall expression pattern throughout pregnancy, and the genes in each cluster and their functional annotation groups represented unique physiological functions during specific periods of pregnancy. Expression of DEGs in cluster 1 increased throughout gestation and declined at term. Genes following this pattern included ACSL1, ACSL4, CYP19A1, ICAM1, IGFBP1, IGFBP5, IL10, LIF, SLPI, SPP1 and UFBP. Some genes such as SLPI, SPP1 and UFBP have been previously investigated [33][34][35], and our results coincide with the expression pattern of those genes in previous studies. Because genes in this cluster were highly expressed during mid-to-late pregnancy when placental growth and fetal development reach the maximum levels [4,5], it is likely that genes in this cluster are involved in support of fetal/placental development and endometrial remodeling. Indeed, functional annotation clustering analysis showed that the gene function of this group includes blood vessel development/morphogenesis, ion/cellular chemical homeostasis, positive regulation of protein modification, and mesenchymal cell differentiation. Genes in cluster 5 showed a pattern of decreased expression on D60 to D90 during pregnancy and included ATP1B1, CTGF, ESR1, ITGB8, OAS1, OGT, and PGR. Functional annotation groups for the genes in this cluster include RNA processing/splicing, regulation of gene expression/ macromolecule biosynthetic process and cellular glucan/polysaccharide/energy reserve metabolic process. Genes in cluster 7 showed a pattern of gradually increasing expression from D12 to D90 of pregnancy and had the highest levels at term. These genes include ATP4B, CTSL2, ESD, FBP, MMP3, NLN, PLAT, and TIMP1. Functional annotation groups in this cluster include protein amino acid glycosylation/glycoprotein biosynthetic processes and modification-dependent macromolecule catabolic/protein catabolic processes, which suggests that remarkable endometrial remodeling occurs in the preparation for parturition at term. Genes in cluster 8 showed highest expression on D12 of pregnancy and decreased thereafter during pregnancy, and include ALCAM, CCL28, EGF, FGF7, HIF1A, IGF1, IL10RB, LPAR3, OCLN, SAL1 and SULT1E1. Functional annotation groups include blood vessel development/differentiation and neuronal projection development/cell morphogenesis. The highest levels of expression on D12 of pregnancy for the genes such as FGF7 [36], LPAR3 [37], SAL1 [38], and SULT1E1 [39,40] have already been reported, confirming the accuracy of microarray analysis in this study.
In WGCNA to determine groups of closely interrelated genes expressed in the uterine endometrium during pregnancy, we identified 5 groups (modules) of closely interrelated genes that had similar patterns of connection strengths to other genes, or high topological overlap. Since hub genes are known to play a critical role in the genetic interaction network, we assumed that these genes function in expression and regulation of other genes in the same module in the uterine endometrium during pregnancy and that they might be markers for predicting specific uterine gene function. Twenty-four hub genes were identified based on connectivity of absolute scaled K values of more than 0.9 in each module. Among them, SPP1, CST6, and RANBP17 were detected as hub genes in the blue module. Indeed, many studies have shown the critical roles of these genes in the uterine endometrium.
The importance of SPP1 in endometrial function during pregnancy is well understood in many species, including humans, rabbits, pigs, and ruminants [41]. SPP1, also called osteopontin and early T-cell activator factor 1, is an acidic glycosylated phosphoprotein found in all body fluids and in a variety of tissues. SPP1 acts on various biological processes such as cell migration, inflammation, activation of B and T cells, and bone formation [41]. In pigs, SPP1 is expressed in uterine endometrial luminal epithelial cells from early stage of pregnancy and is also expressed in glandular epithelial cells from mid-to late-stage of pregnancy [34]. SPP1 can undergo dramatic post-translational modifications such as phosphorylation, glycosylation, and cleavage, forming different sizes of variants including native 75-kDa, 45-and 25-kDa forms [41]. SPP1, an extracellular matrix (ECM) protein, binds multiple integrins through Arg-Gly-Asp-mediated or alternative integrin recognition sequences to effect cell-cell and cell-ECM adhesion [42]. In pigs, SPP1 binds specific integrins and promotes trophectoderm cell migration and attachment to luminal epithelial cells [43]. These suggest that SPP1 plays a critical role in cell-cell and cell-ECM adhesion at the maternal-conceptus interface, especially in pigs, which form a true non-invasive epitheliochorial placentation. Furthermore, several studies have shown that the SPP1 gene is associated with litter size and prenatal survival. The porcine SPP1 gene is included in quantitative trait loci (QTL) associated with reproductive performance such as ovulation rate, litter size, and prenatal survival on chromosome 8 [44]. A polymorphism has been detected in an intron region of the SPP1 gene [45,46], and is associated with birth body weight, growth rate, and carcass traits in the Landrace X Jeju (Korea) black pig F 2 population [46]. In addition, it has been shown that expression of SPP1 is higher in the uterine endometrium of the Meishan pig, a highly prolific pig breed, than in the hyperprolific Large White breed, suggesting the association of SPP1 with placental efficiency during pregnancy [47]. Therefore, the finding of SPP1 as a hub gene in the uterine endometrium during pregnancy in this study and the role of SPP1 in the uterine endometrium during pregnancy and association of the SPP1 gene with porcine litter size and prenatal survival indicates that SPP1 is a critical factor for regulating uterine-conceptus interactions and maintaining a successful pregnancy.
In this study, CST6 was also identified as a hub gene in the uterine endometrium during pregnancy. CST6 (cystatin 6 or cystatin E/M) is an inhibitor of cathepsins (CTSs), including CTSB, CTSL, and legumain (LGMN) [48], which are lysosomal cysteine proteases acting on degradation of extracellular matrix molecules and activation of intracellular pre-proteins [49]. Expression of CTSs and CSTs in the uterine endometrium has been shown in many species including rodents, humans, ruminants and pigs [50][51][52], and it has been suggested that CTSs and CSTs play important roles in endometrial remodeling in the uterine endometrium during the reproductive cycle and pregnancy [52]. In pigs, CTSB, CTSL1, and LGMN are expressed by endometrial luminal epithelial and chorionic epithelial cells at the maternal-conceptus interface during pregnancy and endometrial expression of CTSB and CTSL1 is induced by progesterone [51,53], suggesting that their action is required to remodel uterine endometrial and placental tissues and facilitate transplacental transport of nutrients. CST6 is expressed in the uterine endometrium during the estrous cycle and in pregnancy, and also in the chorionic epithelia of the placental membrane with increasing levels during late pregnancy [53], suggesting that cell type-specific expression and function of CST6 is critical in regulation of CTS action for appropriate maternal-fetal interactions.
The result of this study showed that RANBP17 is a hub gene in the endometrium during pregnancy in pigs. RANBP17 is a member of the importin-β superfamily of nuclear transporter receptors that are involved in the nucleocytoplasmic shuttling of many cargo proteins and interact directly with RanGTP to modulate the compartment-specific binding of their substrate [54]. Importin-βs play essential roles in diverse cellular processes such as gene expression, signal transduction, oncogenesis, cell division, and nuclear envelope assembly [55]. Expression of RANBP17 has been detected in many human and mouse tissues, including brain, heart, lung, pancreas, placenta, and testis, with the highest levels of expression in testis [54], and levels of RANBP17 mRNAs increases with ischemic and dilated cardiomyopathy [56]. Although the specific function of RANBP17 has not been well understood, it has been shown that RANBP17 interacts with a basic helix-loop-helix transcription factor, E12, and increases transcriptional activity of E12, suggesting that E12 may be a cargo protein for RANBP17 [57]. Even though expression and function of RANBP17 in the female reproductive tract have not been determined, our result suggests that RANBP17 may play a critical role in the shuttling and regulation of protein cargoes such as transcription factors in the uterine endometrium during pregnancy.
Using real-time RT-PCR analysis to determine endometrial expression of hub genes, RANBP17 and SPP1 and their correlated genes and to demonstrate the correlation of hub genes determined in WGCNA, we found that the genes positively or negatively correlated with hub genes showed the same or opposite patterns of expression to each other in the uterine endometrium during pregnancy. These results confirmed the correlation among the hub genes. Furthermore, general expression patterns of the hub genes and their correlated genes matched the patterns determined by clustering analysis of DEGs.
In conclusion, our study revealed the complex network characteristics involved in uterine endometrial gene expression during pregnancy stages. Our findings suggest that the diverse patterns of stage-specific gene expression and network connections observed during pregnancy play critical roles in endometrial remodeling and in placental and fetal development in pigs. Further study of the hub genes identified in this study would provide insight into the interrelationship of genes expressed in the uterine endometrium and their endometrial function during pregnancy in pigs.
Supporting Information S1