Defining age- and lactocrine-sensitive elements of the neonatal porcine uterine microRNA-mRNA interactome.

Factors delivered to offspring in colostrum within 2 days of birth support neonatal porcine uterine development. The uterine mRNA transcriptome is affected by age and nursing during this period. Whether uterine microRNA (miRNA) expression is affected similarly is unknown. Objectives were to (1) determine effects of age and nursing on porcine uterine miRNA expression between birth and postnatal day (PND) 2 using miRNA sequencing (miRNAseq) and; (2) define affected miRNA–mRNA interactions and associated biological processes using integrated target prediction analysis. At birth (PND 0), gilts were euthanized, nursed ad libitum, or gavage-fed milk replacer for 48 h. Uteri were collected at birth or 50 h postnatal. MicroRNAseq data were validated using quantitative real-time PCR. Targets were predicted using an established mRNA database generated from the same tissues. For PND 2 versus PND 0 comparisons, 31 differentially expressed (DE) miRNAs were identified for nursed, and 42 DE miRNAs were identified for replacer-fed gilts. Six DE miRNAs were identified for nursed versus replacer-fed gilts on PND 2. Target prediction for inversely correlated DE miRNA–mRNA pairings indicated 20 miRNAs targeting 251 mRNAs in nursed, versus 29 miRNAs targeting 585 mRNAs in replacer-fed gilts for PND 2 versus PND 0 comparisons, and 5 miRNAs targeting 81 mRNAs for nursed versus replacer-fed gilts on PND 2. Biological processes predicted to be affected by age and nursing included cell-to-cell signaling, cell morphology, and tissue morphology. Results indicate novel age- and lactocrine-sensitive miRNA–mRNA relationships associated with porcine neonatal uterine development between birth and PND 2

Introduction

One of the defining characteristics of mammals is lactation. Nursing provides a conduit for delivery of both nutrients and milk-borne bioactive factors (MbFs) from mother to offspring in colostrum (first milk) via a lactocrine mechanism. Colostral MbFs can affect neonatal development as proposed in the lactocrine hypothesis [1, 2]. In pigs and other mammals, female reproductive tract development, initiated prenatally, continues postnatally [3]. Data for the pig show that disruption of normal lactocrine signaling from birth (postnatal day = PND 0), by substitution of a porcine milk replacer for colostrum, altered patterns of uterine gene expression by PND 2, and inhibited uterine endometrial gland development by PND 14 [4]. The observation that gilts consuming minimal amounts of colostrum on their day of birth, as reflected by low immunoglobulin immunocrit ratio [5], displayed reduced live litter size over four parities as adults [6, 7] provided strong support for the lactocrine hypothesis for maternal programming of uterine development and function.

Recently, RNA sequencing (mRNAseq) of the neonatal porcine uterus revealed both age- and lactocrine-sensitive changes in gene expression [8]. More than 3000 genes were differentially expressed in uteri at PND 2 as compared to PND 0. Lactocrine effects were also pronounced on PND 2, when more than 800 genes were differentially expressed by uterine tissues in nursed as compared to replacer-fed gilts. Mechanisms responsible for such age- and lactocrine-sensitive transcriptional changes in the neonatal uterus are unknown.

One mechanism by which transcriptomic changes can occur is through regulation by microRNAs (miRNAs). MicroRNAs are short noncoding RNAs (18–25 nucleotides in length) that regulate post-transcriptional gene expression through translational repression and/or mRNA destabilization and degradation [9]. Micro-RNAs target mRNAs, and a single miRNA can have multiple mRNA targets [10]. Moreover, disruption of miRNA processing affects biological processes governing development and function of the uterus [11]. In pigs, miRNAs were identified in adult endometrium [12-14] and in placental tissues [12, 15, 16] during early pregnancy. Endometrial [17] and placental [18] miRNA–mRNA interactions during pregnancy were also characterized. Little is known about miRNA expression in the developing uterus, or the roles of miRNAs in regulation of uterine gene expression during the perinatal period. Ideally, such studies should integrate miRNA and mRNA expression profiles generated from the same tissues. With these observations in mind, objectives of the present study were to (1) determine effects of age and nursing on the porcine uterine miRNA transcriptome between birth and PND 2 using miRNA sequencing (miRNAseq); and (2) define uterine miRNA–mRNA interactions and associated age- and lactocrine-sensitive biological processes in silico using integrated target prediction analysis.

Materials and methods

Animals and experimental design

Gilts (Sus scrofa domesticus) were born and raised from an established herd of crossbred (Duroc, Hampshire, Yorkshire and Landrace genetics) pigs at the Swine Unit of the New Jersey Agricultural Experiment Station, Rutgers University. All procedures involving animals were reviewed and approved by the Rutgers Institutional Animal Care and Use Committee and conducted in accordance with the Guide for the Care and Use of Agricultural Animals in Agriculture Research and Teaching [19]. Consideration was given to ensure that sows nursed litters of similar size and that treatments were balanced for potential effects of litter (n = 8).

At birth, gilts (n = 12) were assigned randomly to be either (1) sacrificed on PND 0, prior to nursing (n = 4); (2) nursed ad libitum from birth through 48 h of age (PND 2N, n = 4); or (3) gavage-fed a nutritionally complete, commercial pig milk replacer (30 mL/kg BW/2 h; Advance Liqui-Wean MSC Specialty Nutrition; Carpentersville, IL, USA) from birth through 48 h of age (PND 2R, n = 4) [20]. Gilts were euthanized and uterine tissues were collected on either PND 0 or 50 h of age. Uteri were trimmed of associated tissues and uterine wet weights (mg) were recorded. Uterine tissue samples were immersed in RNAlater (Life Technologies, Carlsbad, CA, USA) and stored at –80°C until total RNA was extracted.

Uterine RNA isolation and analysis

Total RNA (including miRNA and mRNA) was isolated from 50–60 mg of whole uterine tissue from one uterine horn/animal using the miRNeasy Mini Kit (Qiagen Inc., Valencia, CA, USA) following manufacturer's protocol. RNA quantity was determined using a Qubit 2.0 Fluorometer (Invitrogen; Carlsbad, CA, USA) and RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Applied Biosystems; Carlsbad, CA, USA). Samples with an RNA integrity number ≥8.0 were used for library preparation for miRNAseq.

Preparation of microRNA libraries

MicroRNAseq was performed at the Genomic Services Laboratory, HudsonAlpha Institute for Biotechnology (Huntsville, AL, USA). Total RNA (500 ng) from each uterine sample was used for RNA library preparation using the NEBNext Small RNA Library Prep Set for Illumina (New England Biolabs Inc., Ipswich, MA, USA) according to the manufacturer's protocol. Briefly, adapters were ligated to total RNA, multiplex primers were hybridized, and reverse transcription was accomplished using SuperScript III RT (Life Technologies, Grand Island, NY, USA) for 1 h at 50°C. Bar codes with uniquely indexed primers were attached to each cDNA library and amplified through six PCR cycles. Following PCR amplification, purification was done using the QIAquick PCR purification kit (Qiagen Inc., Valencia, CA, USA). Size selection of the libraries was then performed using a 3% dye-free agarose gel on the Pipin prep instrument (Sage Science, Beverly, MA, USA). Post size-selected miRNA library concentration was assessed using a Qubit 2.0 Fluorometer and library quality was determined using a DNA High Sense chip on an Agilent 2100 Bioanalyzer. Further, library quantification was performed using the quantitative real-time polymerase chain reaction (qPCR)-based KAPA Biosystem Library Quantification kit (Kapa Biosystems Inc., Woburn, MA, USA). Individual sample libraries were diluted to a final concentration of 1.25 nM and equimolar amounts of each sample were pooled prior to sequencing.

MicroRNA sequencing and data analysis

MicroRNAseq was performed using an Illumina HiSeq 2500 instrument (Illumina Inc., San Diego, CA, USA) at 50 bp single end condition, generating approximately 15 million reads per sample. Quality control checks on raw sequence data from each sample were performed using FastQC (Babraham Bioinformatics, London, UK). Raw reads were imported on a commercial data analysis platform Avadis NGS (Strand Scientifics, CA, USA). Adapter trimming was done to remove ligated adapters from the 3' end of the sequenced reads with only one mismatch allowed; poorly aligned 3' ends were also trimmed. Sequences shorter than 15 nucleotides in length were excluded from further analysis. Trimmed reads with low qualities (base quality score less than 30, alignment score less than 95, mapping quality less than 40) were removed. Filtered reads were used to extract and count miRNAs which were annotated with miRBase release 18 database [21-25]. Reads were grouped according to their respective identifiers followed by quantification of miRNA abundance [26]. Differentially expressed miRNAs, based on fold change (≥±2.0), were identified with respect to effects of neonatal age (PND 2N versus PND 0; PND 2R versus PND 0) and nursing (PND 2N versus PND 2R). Probability values for each differentially expressed miRNA were estimated by z-score calculations using a false discovery rate of 0.05. Data were subjected to principal component analysis and hierarchical clustering created with Spearman correlation coefficient. Relative fold change in miRNA abundance was illustrated using volcano plots generated for each comparison using R Programming (GNU General Public License; www.r-project.org).

Quantitative real-time polymerase chain reaction

TaqMan Advanced miRNA assays (Life Technologies) were used for qPCR validation of the miRNAseq data. The same RNA used to generate cDNA libraries for miRNAseq was also used to validate results by qPCR. Uterine RNA from individual animals was pooled to create PND 0, PND 2N, and PND 2R samples. For qPCR validation, miRNAs were selected at random from the population of miRNAs that increased or decreased in at least one of the three conditions. Reverse transcription of total RNA samples (5 ng) was done using gene-specific primers and the TaqMan Advanced MicroRNA cDNA Synthesis kit. Primer pairs specific to each miRNA and TaqMan Fast Advanced PCR Master Mix were used for amplification per manufacturer's recommendation. Primers were evaluated for quality by amplifying serial dilutions of the cDNA template. Control qPCR reactions included substitution of water in place of primers and template to ensure specific amplification in all assays. Dissociation curves for primer sets were evaluated to ensure that no amplicon-dependent amplification occurred.

Geometric means of qPCR CT (cycle threshold) values for miRNAs with high, medium, and low reads (ssc-miR-21, ssc-miR-101, and ssc-miR-127, respectively) were used for normalization as described by Maalouf [27]. Expression of reference miRNAs did not change significantly between samples from PND 0, PND 2N, or PND 2R gilts, as determined using data from miRNAseq and qPCR analyses. Data generated by qPCR were analyzed using the ΔΔCT method as described by Applied Biosystems (ABI User Bulletin 2, 2001). Pearson correlation coefficients were determined to compare miRNA expression fold-change results obtained by miRNAseq and qPCR.

Integrated target prediction analysis

Human orthologs of mRNAs were identified using NCBI BLAST and a custom computer program written in Python (www.python.org) using NCBI modules within Biopython [28], as described elsewhere [8]. Human orthologs of miRNAs were identified using miRBase release 21 database [21-25]. Relationships between differentially expressed miRNAs and their respective differentially expressed mRNA targets were determined using Qiagen's Ingenuity Pathway Analysis MicroRNA Target Filter (IPA, Qiagen Redwood City, CA; www.qiagen.com/ingenuity). Differentially expressed uterine mRNA data, generated from the same tissues, were reported previously [8] and can be found in the GEO repository under series accession number GSE72388. Gene enrichment and functional annotation analyses were conducted using the Database for Annotation, Visualization, and Integrated Discovery (DAVID 6.7; http://david.abcc.ncifcrf.gov/) [29, 30] and IPA to identify enriched biological functions among differentially expressed gene transcripts, including miRNAs and mRNAs.

Results

Effects of age and nursing on neonatal uterine microRNA expression

Uterine miRNA data for PND 0, PND 2N, and PND 2R treatments grouped independently according to principal component analysis and Spearman correlation (Figure 1A and B). By examining age (PND 2N versus PND 0; PND 2R versus PND 0) and nursing (PND 2N versus PND 2R), differentially expressed miRNAs (≥two-fold, P < 0.05) were determined in neonatal porcine uteri. Expression analyses are illustrated as volcano plots (Figure 2A–C). Results indicated 31 differentially expressed miRNAs on PND 2N when compared to PND 0 and all of these miRNAs decreased with age (Figure 2A). When gilts were fed milk replacer for 2 days from birth, 42 mi-RNAs were differentially expressed on PND 2R as compared to PND 0 (Figure 2B). Of these differentially expressed miRNAs, one was increased and 41 were decreased on PND 2R compared to PND 0 (Figure 2B). A total of six miRNAs were differentially expressed in the uteri of nursed (PND 2N) as compared to replacer-fed (PND 2R) gilts (Figure 2C). Of these differentially expressed miRNAs, three were increased and three were decreased on PND 2N as compared to PND 2R (Figure 2C). Differentially expressed miRNAs in each comparison and their respective fold changes are presented in Supplementary Tables S1–S3. Quantitative real-time PCR results for six miRNAs used to validate miRNAseq data are shown in Figure 2D–F. Positive correlations between miRNAseq and qPCR results were identified (r = 0.76, P < 0.01).

Figure 1.

(A) Principal component analysis (PCA) plot and (B) Spearman correlation heat map from miRNAseq analysis. (A) PCA plot of pooled neonatal porcine uterine samples showing the clustering of expressed miRNAs for PND 0 (green), PND 2N (blue), and PND 2R (purple). (B) Heat map of pairwise correlations between pooled samples based on the Spearman correlation coefficients. Light green represents highest correlations and bright red represents lowest correlations.

Figure 2.

Volcano plots (left) demonstrating differences in uterine expression of miRNAs between (A) PND 2N versus PND 0, (B) PND 2R versus PND 0, and (C) PND 2N versus PND 2R. Upregulated miRNAs are indicated in red (fold change ≥ 2) and downregulated miRNAs are denoted in green. Black indicates miRNAs that did not change between groups. The horizontal line indicates P = 0.05. For each comparison, the total number of differentially expressed miRNAs is given (Δ value). The number of miRNAs that were downregulated (green) or upregulated (red) are shown. Results of qPCR validation (right) for six miRNAs identified by miRNAseq for (D) PND 2N versus PND 0; (E) PND 2R versus PND 0; (F) PND 2N versus PND 2R. White bars indicate miRNAseq fold change; black bars indicate qPCR fold change. An overall positive correlation between miRNAseq and qPCR results was identified (r = 0.76, P < 0.01).

To examine the overlap of differential miRNA expression, a Venn diagram was created to categorize expression domains through a three-way comparison (Figure 3). Overall, 62 miRNA species were differentially expressed (Figure 3). There was a single differentially expressed miRNA species (miR-184) that was common to each group (Figure 3). Comparison of differentially expressed miRNAs in PND 2N versus PND 0 revealed 16 uniquely expressed mi-RNAs, 13 of which were shared with PND 2R versus PND 0 and one of which was shared with PND 2N versus PND 2R (miR-296-5p; Figure 3). For PND 2R versus PND 0, 27 differentially expressed miRNAs were uniquely expressed and one was common with PND 2N versus PND 2R (miR-345-5p; Figure 3). A total of three differentially expressed miRNAs were uniquely expressed between nursed and replacer groups on PND 2 (Figure 3). Altogether, common and unique differentially expressed uterine miRNAs were identified between treatments.

Figure 3.

Venn diagram illustrating unique and overlapping differential miRNA expression domains for the three group comparisons [PND 2N versus PND 0 (green); PND 2R versus PND 0 (purple); PND 2N versus PND 2R (blue)]. Values indicate the number of differentially expressed miRNAs associated with unique and overlapping domains. Differentially expressed miRNAs determined to be unique to each comparison are listed in descending order for absolute fold change (AFC > 2-fold; red = upregulated, green = downregulated). Superscript 1 indicates miRNAs with same name found at different genomic loci.

Integrated target prediction analyses

Based on the previous results, it was hypothesized that a portion of the uterine transcriptomic changes observed due to age and nursing were regulated by differentially expressed miRNAs. Therefore, integrated miRNA–mRNA analyses were conducted, via Qiagen's Ingenuity Pathway Analysis (IPA) miRNA Target Filter, to explore differentially expressed uterine miRNAs and mRNA targets between treatments. Figure 4 summarizes the workflow and output for IPA data integration and bioinformatics analyses. Table 1A–C lists the differentially expressed miRNAs and top five differentially expressed mRNA targets. Supplementary Tables S4–S6 list all differentially expressed mRNA targets for each comparison.

Figure 4.

High-throughput sequencing workflow and output. Uterine tissues were obtained from gilts (1) at birth on PND 0, before nursing; (2) after nursing ad libitum from birth through 48 h of age (PND 2N); or (3) after gavage feeding a commercial pig milk replacer for 48 h (PND 2R). Total uterine RNA was isolated and both small RNAs and mRNAs were sequenced as described in section Materials and Methods. Asterisk indicates mRNAseq data for uterine mRNAs as described in Rahman et al [8]. The numbers of differentially expressed miRNAs and mRNAs for each comparison are presented. Integrated target prediction analysis was conducted using IPA as described in section Materials and Methods. The numbers of differentially expressed miRNAs targeting differentially expressed mRNAs are presented for each comparison. Individual miRNAs targeting mRNAs are listed (bottom; absolute fold change > 2; red = upregulated, green = downregulated).

Table 1.

Integrated target prediction analysis for differentially expressed uterine miRNAs and mRNA targets between (A) PND 2N versus PND 0, (B) PND 2R versus PND 0, and (C) PND 2N versus PND 2R as identified by IPA.

(A) PND 2N versus PND 0
miRNA	No. of targeted mRNAs2	Top 5 differentially expressed target mRNAs
let-7d-5p	30	AGXT2, PTPRO, SLC10A2, DAPK2, DUSP9
miR-106a	26	ERICH3, PTPRO, ZBTB38, DAPK2, ESR1
miR-1343	26	C14orf105, PLEKHB1, LRRN2, SRRM4, PTPRO
miR-135	22	GABRG1, CXCL10, STMN4, FAM25A, HNRNPA3
miR-136	17	ATP5F1, SRRM4, FAM229B, SLC7A3, TMEM232
miR-146a-5p	19	CCL8, CAMP, S100A12, SEPT14, CNTF
miR-184	10	PLEKHB1, BST2, CX3CR1, F5, C21orf62
miR-193a-5p	16	C14orf105, NRXN1, THEMIS2, CPXM2, FXYD3
miR-205	11	WDR77, C17orf97, ZBTB38, AFF3, CXorf21
miR-206	32	CXCL11, CCL2, FXYD3, SRRM4, BRI3BP
miR-210	7	CAPN9, SF3B3, PCYT1B, CEND1, ATXN10
miR-221-3p	14	BBOX1, GABRG1, SEPT14, NRXN1, CXCL11
miR-296-5p	42	CXCL10, HRASLS5, SERPINA1, FXYD3, SRRM4
miR-339-5p	29	GABRG1, NRXN1, RBM3, UNC45B, PTPRO
miR-374a-5p	13	CCL8, GABRG1, CCL2, ZBTB38, AFF3
miR-451	5	CXorf21, BATF, NR5A2, KIAA0101, FAM159B
miR-455-3p	20	XCL1, MS4A2, CAMK2N1, ANKRD34C, CBLN2
miR-486	11	BUB1B, CAMK2N1, AFF3, VTCN1, INMT
miR-504	26	UNC45B, ZBTB38, SAA4, FEZ1, EPHA8
miR-708-5p	28	C14orf105, PLEKHB1, SKA1, COL17A1, SRRM4
(B) PND 2R versus PND 0
miRNA	No. of targeted mRNAs 2	Top 5 differentially expressed target mRNAs
miR-106a	35	VCL, DAPK2, PTH, SCRT2, PTPRD
miR-10a-5p	32	ATP5F1, RBM3, OPALIN, SKA1, UNC45B
miR-129a-5p	19	AFF3, CALM1, AUTS2, SF3B3, CAMK2N1
miR-129b	30	SNTN, CFAP61, VNN2, BPI, KLK15
miR-135	39	KALRN, OPALIN, STMN4, SLC9A4, PTPRD
miR-149	61	VCL, FAM216B, CAPN8, TDRD3, MURC
miR-181b	34	KALRN, TRIM64/TRIM64B, THBS4, CCL8, AGT
miR-184	17	FOXA2, BST2, DLX3, SCRT2, C21orf62
miR-185	54	KALRN, UNC45B, LRRC38, GPX6, SOX5
miR-193a-5p	30	VCL, KLK4, DIO2, THEMIS2, C14orf105
miR-196a	18	CALM1, SCRT2, PDE11A, PABPC1L2A, INMT
miR-206	38	CWC15, SOX5, CXCL11, CCL2, CALM1
miR-221-5p	22	BBOX1, CXCL11, ACTC1, SEPT14, RBP2
miR-30b-3p	70	TRIM64/TRIM64B, VCL, AGT, SNTN, PATE3
miR-331-5p	15	SNTN, HR, INMT, FAM19A3, GEN1
miR-335	19	KALRN, NPAS4, GBP4, NOS1, LSM1
miR-345-5p	103	CRISPLD2, PIEZO2, GRIK3, HSPA12A, EZR
miR-34c	59	CAV3, WDR77, HOXA13, VCL, DAPK2
miR-361-3p	72	GGT1, TMEM242, FAM92B, ACTC1, NPAS4
miR-362	14	TBXAS1, ANGPTL7, PTPRD, F5, ANKRD34B
miR-451	11	MEGF6, CRELD2, GDAP1, CXorf21, FAM159B
miR-486	18	AFF3, BUB1B, CAMK2N1, INMT, UBTF
miR-504	41	SAA4, KALRN, SAA2-SAA4, UNC45B, FAM196A
miR-505	23	TDRD3, MPL, ASCL4, LZTFL1, NOS1
miR-615	14	BPIFB4, CFAP61, LRRN2, ITSN1, SESN2
miR-628	22	FAM229B, SKA1, TMEM242, CALM1, CRYBA4
miR-671-5p	34	WDR77, OPALIN, VNN2, PTPRD, SF3B3
miR-769-3p	19	VCL, GPX6, DIO2, NPAS4, SCRT2
miR-9843-3p	32	MYH7, SOCS2, EED, PTPRD, C17orf80
(C) PND 2N versus PND 2R
miRNA	No. of targeted mRNAs 2	Top 5 differentially expressed target mRNAs
miR-184	4	FOXA2, LIPG, NOS1, SIRPA
miR-296-5p	49	RAB11FIP5, HEYL, GPR37L1, ADCY2, OXCT1
miR-345-5p	20	GALNTL6, OXCT1, GRIK3, PCDH10, HSPA12A
miR-490-3p	6	KALRN, CLEC18A/CLEC18C, SP5, SLC5A3, MAGI1
miR-582	7	KALRN, PBK, TTN, AGTPBP1, YTHDF3

For case where miRNA expression is downregulated and mRNA expression is upregulated or miRNA expression is upregulated and mRNA expression is downregulated.

Individual miRNAs can target overlapping target populations of mRNAs.

With respect to pigs nursed from birth (PND 2N versus PND 0), there were 31 differentially expressed miRNAs and 3,283 differentially expressed mRNA transcripts. Of the 31 differentially expressed miRNAs, 20 had 251 predicted mRNA targets (Figure 4; Table 1A). Relative abundance of all of these miRNAs decreased, and their respective mRNA targets increased, on PND 2 in nursed gilts (Figure 4; Table 1A).

For pigs fed replacer from birth (PND 2R versus PND 0), there were 42 differentially expressed miRNAs and 4,662 differentially expressed mRNA transcripts (Figure 4). Of these 42 miRNAs, 29 were predicted to target 585 mRNA transcripts (Figure 4). Relative abundance of all but one of these miRNAs decreased (Figure 4). The single miRNA that increased (miR-345-5p) on PND 2 in replacer-fed gilts was predicted to have 103/585 differentially expressed mRNA targets (Table 1B). The other 28 differentially expressed miRNAs and the number of predicted mRNA targets are presented in Table 1B.

For nursed versus replacer-fed gilts on PND 2 (PND 2N versus PND 2R), there were six differentially expressed miRNAs and 896 differentially expressed mRNA transcripts. Five of the six lactocrine-sensitive miRNAs were predicted to target 81 mRNAs (Figure 4). These five included miR-184, miR-296-5p, miR-345-5p, miR-490-3p, and miR-582 (Figure 4; Table 1C). The number of predicted mRNA targets and the top five differentially expressed mRNAs are presented in Table 1C.

Integrated functional annotation analyses

Enriched biological processes in neonatal porcine uteri associated with mRNAs targeted by miRNAs affected by age in nursed gilts (PND 2N versus PND 0) as identified by DAVID are shown in Table 2A. Functional annotation by DAVID analysis revealed terms including “defense response,” “inflammatory response,” “response to wounding,” “immune response,” and “cellular homeostasis” (Table 2A). Selected functional annotation categories identified by IPA for inversely correlated, differentially expressed transcripts (miRNAs and mRNAs) in neonatal porcine uteri between PND 2N versus PND 0 are shown in Figure 5A and Supplementary Table S7. MicroRNA–mRNA interactions were predicted to be involved with multiple biological processes including “cellular movement,” “cell-to-cell signaling and interaction,” “cellular function and maintenance,” and “tissue morphology” (Figure 5A).

Table 2.

Top 10 enriched biological processes in neonatal porcine uteri associated with miRNA–mRNA interactions between (A) PND 2N versus PND 0, (B) PND 2R versus PND 0, and (C) PND 2N versus PND 2R as identified by DAVID functional annotation analysis.

(A) PND 2N versus PND 0
Functional terms of overrepresented biological processesa	Enrichment scoreb
Defense response (24)	4.8
Inflammatory response (15)	3.7
Response to wounding (19)	3.3
Immune response (22)	3.2
Acute-phase response (5)	2.6
Locomotory behavior (11)	2.3
Cellular homeostasis (15)	2.2
Behavior (15)	2.2
Response to calcium ion (5)	2.1
Chemotaxis (8)	2.1
(B) PND 2R versus PND 0
Functional terms of overrepresented biological processes a	Enrichment score b
Vitamin metabolic process (10)	3.2
Cellular chemical homeostasis (25)	2.9
Ion homeostasis (26)	2.8
Cell–cell signaling (34)	2.7
Cellular ion homeostasis (24)	2.7
Fat-soluble vitamin metabolic process (6)	2.4
Chemical homeostasis (29)	2.4
Cellular metal ion homeostasis (15)	2.4
Cellular homeostasis (27)	2.4
Response to temperature stimulus (9)	2.3
(C) PND 2N versus PND 2R
Functional terms of overrepresented biological processes a	Enrichment score b
Cell differentiation in hindbrain (3)	2.9
Regulation of neurotransmitter levels (4)	2.5
Neurotransmitter metabolic process (3)	2.4
Neuron differentiation (7)	2
Cytosolic calcium ion homeostasis (4)	1.9
Transmission of nerve impulse (6)	1.8
Regulation of neuron differentiation (4)	1.8
Regulation of membrane potential (4)	1.7
Chemical homeostasis (7)	1.7
Cellular ion homeostasis (6)	1.7

Values within parentheses indicate the number of annotated mRNAs targeted by miRNAs that are involved with the corresponding functional term.

Enrichment scores were calculated by taking the geometric mean of the P-values associated with the differentially expressed transcripts involved in the corresponding annotation cluster (in –log10 scale).

Figure 5.

Selected functional annotation categories for differentially expressed transcripts (miRNAs and mRNAs) associated with miRNA–mRNA interactions between (A) PND 2N versus PND 0; (B) PND 2R versus PND 0; and (C) PND 2N versus PND 2R as identified by IPA. Values within parentheses indicate the number of annotated mRNAs targeted by miRNAs that are involved in the corresponding functional term. Enrichment scores were calculated by taking the geometric mean of the P-values associated with differentially expressed transcripts (in –log10 scale). Dashed vertical line indicates P = 0.05.

Enriched biological processes in neonatal porcine uteri associated with mRNAs targeted by miRNAs affected by age in replacer-fed gilts (PND 2R versus PND 0) as identified by DAVID are shown in Table 2B. Functional annotation by DAVID analysis revealed terms including “cell-cell signaling,” “cellular ion homeostasis,” and “cellular homeostasis” (Table 2B). Selected functional annotation categories identified by IPA for inversely correlated, differentially expressed transcripts (miRNAs and mRNAs) in neonatal porcine uteri between PND 2R versus PND 0 are shown in Figure 5B and Supplementary Table S8. MicroRNA–mRNA interactions were predicted to be involved with multiple biological processes similar to the PND 2N versus PND 0 comparison, including: “cell morphology,” “cellular growth and proliferation,” “tissue morphology,” and “cell-to-cell signaling and interaction” (Figure 5B). However, the miRNAs involved and the mRNA targets were distinct (Supplementary Table S8).

Enriched biological processes in neonatal porcine uteri associated with mRNAs targeted by miRNAs affected by nursing (PND 2N versus PND 2R) as identified by DAVID are shown in Table 2C. Functional annotation by DAVID analysis identified categories of interest, including “cellular ion homeostasis,” “chemical homeostasis,” and “cellular component morphogenesis” (Table 2C). Further investigation of lactocrine-sensitive miRNAs and their mRNA targets by comparison of PND 2N versus PND 2R uterine transcripts by IPA identified “cell-to-cell signaling and interaction,” “organ, organismal, tissue, and cellular development,” and “cell, organ, and tissue morphology” (Figure 5C; Supplementary Table S9).

For the three miRNAs unique to the PND 2N versus PND 2R comparison, two (miR-490-3p and miR-582) were predicted to target 12 mRNAs. Biological processes associated with these miRNA–mRNA relationships, identified by DAVID analysis, included “organ morphogenesis” and “muscle tissue morphogenesis” (data not shown). Using IPA, functional annotation categories previously observed, such as “cell, organ, and tissue morphology,” “organ, organismal, tissue, and cellular development,” as well as “reproductive system development and function” were identified (Supplementary Table S10).

Predicted microRNA–mRNA uterine interactome networks affected by age and nursing

Ingenuity Pathway Analysis enables illustration of miRNA–mRNA interactome networks. Shown here are examples for “cell-to-cell signaling and interaction” networks illustrating the size, directionality (positive or negative fold change), and overlap for each comparison (Figure 6). Additional examples for the “tissue morphology” network are provided in Supplemental Figure S1. In each case, upregulated transcripts are shown in red and downregulated transcripts are green.

Figure 6.

Integrated uterine miRNA–mRNA interactomes illustrating effects of (A) PND 2N versus PND 0; (B) PND 2R versus PND 0; and (C) PND 2N versus PND 2R on cell-to-cell signaling and interaction. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. Inset provides an opportunity to study detailed relationships within the interactome network. IPA legend key (left): (1) cytokine/growth factor; (2) enzyme; (3) G-protein coupled receptor; (4) growth factor; (5) ion channel; (6) kinase; (7) peptidase; (8) phosphatase; (9) transcription regulator; (10) transmembrane receptor; (11) transporter; (12) other; and (13) mature miRNA.

Interactome networks presented here depict important dynamic relationships associated with effects of age and nursing on the uterine transcriptome. The skeletons of interactome networks are identical, illustrating overlap in response domains (Figure 6). These relationships are seen readily in inset details (Figure 6A–C) and can be studied in Supplemental Figures S2–S4. Note, the number of affected elements in “cell-to-cell signaling and interaction” networks (denoted by red and green) differs between experimental groups (Figure 6A–C), as does the direction (up or down) and degree (color intensity) of change for specific elements of each interactome. For example, forkhead box A2 (FOXA2), a predicted miRNA target, is unchanged between PND 0 and PND 2N (Figure 6A inset), upregulated on PND 2 compared to PND 0 in replacer-fed gilts (Figure 6B inset), and downregulated in nursed as compared to replacer-fed gilts on PND 2 (Figure 6C inset). Similar complex relationships were identified for other interactome networks (Supplemental Figures S1, S5–S7, and other data not shown).

Discussion

Through lactocrine mechanisms, bioactive factors are delivered from mother to offspring as a consequence of nursing and support neonatal development [1, 2, 31–33]. Recently, both age and lactocrine effects on the neonatal porcine uterine transcriptome were defined [8]. Using the same tissues, present results illustrate similar effects on miRNA expression profiles associated with uterine development in nursed as compared to replacer-fed gilts between birth and PND 2. Further, integrated target prediction analyses conducted in silico extended previous mRNAseq findings [8] by revealing novel age- and lactocrine-sensitive miRNA–mRNA interactions and biological processes associated with porcine neonatal uterine development. Results suggest that miRNAs may contribute to regulation of uterine gene expression post-transcriptionally between birth and PND 2 and that imposition of a lactocrine-null state by replacer feeding during this period dysregulates the normal uterine developmental program at transcriptional and, potentially, post-transcriptional levels.

The first few days of neonatal life encompass a period of organizational transition for the porcine endometrium. Differentiation of uterine glandular epithelium (GE) from luminal epithelium (LE) is marked by onset of estrogen receptor (ESR1) expression in nascent GE, evident by 24 h postnatal [6]. Transition of the endometrium from a nonproliferative to a proliferative state is associated with progression of morphogenetic events supportive of uterine gland development by PND 3 [34-36]. The fact that uterine histogenesis proceeds normally prior to PND 60 in gilts ovariectomized at birth [34, 37] emphasizes the importance of extraovarian postnatal uterotrophic support. Nursing for 12 h from birth is necessary to support development of porcine uterine and cervical tissues to PND 2 [38], and imposition of a lactocrine-null state from birth by milk-replacer feeding altered patterns of endometrial cell proliferation and cell compartment-specific gene expression during this period [4]. Present data extend earlier observations [4, 8], indicating that nursing and maternally derived, lactocrine-active factors constitute one source of such uterotrophic support.

Neonatal uterine miRNAs expressed at PND 0, PND 2N, and PND 2R clustered independently, as indicated by principal component analysis, indicating distinct miRNA response domains for each condition. Characteristics of age- and lactocrine-sensitive response domains can be defined in terms of size, directionality (positive or negative fold change), and overlap. For age comparisons (PND 2N versus PND 0 and PND 2R versus PND 0), absolute domain size was smaller in nursed (31 differentially expressed miRNAs) as compared to replacer-fed gilts (42 differentially expressed miRNAs) by PND 2. Moreover, with the exception of a single miRNA in PND 2R gilts, relative uterine miRNA expression decreased from birth to PND 2 in both groups. Consequently, as supported by results of in silico target prediction analysis, imposition of a lactocrine-null condition in replacer-fed gilts resulted in a larger number of differentially expressed miRNAs with a larger number of predicted mRNA targets on PND 2. While overlap in differential uterine miRNA expression for age comparisons was substantial, the majority of differentially expressed miRNAs were unique to each domain (PND 2N versus PND 0 and PND 2R versus PND 0). These miRNA relationships and differences in the size of mRNA target pools (251 for nursed versus 585 for replacer-fed) suggest substantial divergence in transcriptomic signatures by PND 2 between these two conditions, with potential to effect very different organizational trajectories [4, 8]. The size of the response domain for uterine miRNA expression on PND 2 in nursed versus replacer-fed gilts (six differentially expressed miRNAs) was smaller than those defined for both age comparisons. This undoubtedly reflects the fact that the age component is absent from this comparison. While overlap was observed in all domains, three of six differentially expressed miRNAs were unique to the PND 2N versus PND 2R domain, indicating a lactocrine-specific response.

A few miRNAs identified here, including miR-296-5p, miR-345-5p, and miR-582, are expressed in adult uterine and placental tissues [12, 13, 15, 18, 39]. However, none were identified in neonatal uterine tissues prior to this report. Two miRNAs for which uterine expression was both unique and lactocrine sensitive on PND 2 included miR-582, which increased, and miR-345-5p, which decreased in nursed as compared to replacer-fed gilts. These miRNAs were also expressed differentially in adult human endometrium [39]. Specifically, expression of both miR-582-5p and miR-345 decreased in late proliferative as compared to mid-secretory endometrium [39]. Authors proposed that these miRNAs could function to suppress endometrial cell proliferation during the secretory phase of the menstrual cycle [39]. In the pig, imposition of a lactocrine-null state for 2 days from birth reduced epithelial cell proliferation by PND 2 in both LE and GE [4]. Thus, lactocrine effects on neonatal uterine miR-582 and miR-345-5p expression observed here may be affecting events associated with regulation of endometrial cell proliferation and nascent uterine gland development.

Uterine expression of miR-451, which decreased due to age in both nursed and replacer-fed gilts, increased in estrogen-treated, ovariectomized mice [40]. Additionally, wild-type mice that received endometrial fragments from miR-451-deficient mice displayed increased uterine expression of fibrinogen alpha chain precursor and fewer endometriotic lesions [41]. This was suggested to reflect effects on the biochemistry of cell adhesion. The extent to which miR-451 may be regulating similar organizationally important processes in the neonatal porcine uterus remains to be determined.

The fact that miRNA–mRNA interactions have functional consequences for mRNA stability and translational efficiencies [9, 10, 39] emphasizes the importance of defining such interactions. In the pig, miRNA–mRNA interactions are implicated in uteroplacental and pregnancy biology. Unique uterine miRNA signatures were defined during implantation [12–14, 42] and in pregnancy [15, 16]. Integrated analysis of miRNA–mRNA networks in placental tissues of Large White and Qingping sows revealed mRNAs and miRNAs associated with onset of labor, as well as a subset of genes that may play a role in regulation of gestation length [18]. Endometrial mRNAs and miRNAs in pregnant sows were also associated with extreme prolificacy phenotypes [17]. Present results now implicate miRNA–mRNA interactions in neonatal porcine uterine development.

Potential mRNA targets and related interactions were identified here using DAVID and IPA. Results generated using DAVID provided information regarding predicted mRNA targets of differentially expressed miRNAs [30]. Prediction of miRNA–mRNA interactions and their potential functions required use of IPA. This program incorporates miRNA and mRNA data, fold-change values, and experimentally validated miRNA–mRNA interactions from miRecords [43], TarBase [44], and predicted interactions from TargetScan [45]. Use of both DAVID and IPA enabled objective assessment of miRNA–mRNA interactions within this complex dataset.

Present results, taken together with data for the neonatal porcine uterine mRNA transcriptome [8], implicate many biological processes associated with uterine development between birth and PND 2. Among enriched biological processes associated with miRNA–mRNA interactions related to effects of age and nursing were those likely to be involved with cytodifferentiative and morphogenetic events affecting development of the neonatal uterine wall [46]. Biological processes including “cell-to-cell signaling and interaction,” “cellular assembly and organization,” “cellular function and maintenance,” and “cell morphology and tissue morphology” were determined to be age and lactocrine sensitive. Ingenuity Pathway Analyses of miRNA–mRNA interactions unique to nursing (PND 2N versus PND 2R) identified similar biological processes. Consistently, previous RNAseq analyses [8] identified similar processes for the mRNA transcriptome alone. Observations are consistent with data indicating that events associated with imposition of a lactocrine-null state for 2 days from birth set the stage for significant changes in endometrial development that are evident by PND 14 in the neonatal pig [4].

Careful evaluation of related miRNA–mRNA networks showed that, while some elements of these networks are shared, others are unique to nursed and lactocrine-null domains. Interactome networks show the size of the response domain, number of affected elements, direction (up or down) and degree of expression (color intensity) for each element of the domain of comparison. For example, the number of elements for “cell-to-cell signaling and interaction” is the same for each comparison, while the number of affected elements and direction of predicted effects are comparison specific. By way of illustration, expression of FOXA2, implicated as a mediator of murine uterine gland genesis [47] and a target for miR-184, is predicted to be unaffected by age in nursed gilts, upregulated in replacer-fed gilts between birth and PND 2, and downregulated on PND 2 in nursed as compared to replacer-fed gilts. These predicted effects on FOXA2 expression are likely to reflect complex miRNA–mRNA interactions affected by age and nursing. The miRNA–mRNA interactions identified here remain to be proven functionally.

Robust computational models used to predict miRNA targets notwithstanding, precise identification of targeted transcripts remains challenging [45]. While imperfect, the in silico approach to miRNA target prediction taken here benefits from the fact that expression profiling of both miRNA and mRNA populations was conducted on the same uterine tissues [8]. Complimentary miRNA and mRNA data generated through these efforts, used in concert with in silico target prediction analyses, will support studies designed to define functional miRNA–mRNA interactions.

Results of miRNAseq analyses identified both age- and lactocrine-sensitive miRNAs in the neonatal porcine uterus between birth and PND 2. Integrated target prediction analyses extended previous mRNAseq findings by revealing novel age- and lactocrine-sensitive miRNA–mRNA interactions predicted to regulate biological processes associated with porcine neonatal uterine development during this period. Results reinforce the importance of nursing from birth and lactocrine signaling on establishment of an optimal uterine developmental program [4, 8].

Supplementary data are available at online.

Supplemental Figure S1. Integrated uterine miRNA–mRNA interactomes illustrating effects of (A) PND 2N versus PND 0, (B) PND 2R versus PND 0, and (C) PND 2N versus PND 2R on tissue morphology. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. Inset provides an opportunity to study detailed relationships within the interactome network. IPA legend key (left): (1) cytokine/growth factor; (2) enzyme; (3) G-protein coupled receptor; (4) growth factor; (5) ion channel; (6) kinase; (7) peptidase; (8) phosphatase; (9) transcription regulator; (10) transmembrane receptor; (11) transporter; (12) other; and (13) mature miRNA.

Supplemental Figure S2. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2N versus PND 0 on cell-to-cell signaling and interaction. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Supplemental Figure S3. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2R versus PND 0 on cell-to-cell signaling and interaction. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Supplemental Figure S4. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2N versus PND 2R on cell-to-cell signaling and interaction. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Supplemental Figure S5. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2N versus PND 0 on tissue morphology. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Supplemental Figure S6. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2R versus PND 0 on tissue morphology. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Supplemental Figure S7. Enlarged integrated uterine miRNA–mRNA interactome illustrating effects of PND 2N versus PND 2R on tissue morphology. Red denotes increased and green denotes decreased transcript expression (P < 0.05). Color intensity indicates degree of change. See legend (bottom right) for symbol identification.

Click here for additional data file.