miRNA expression profiles of peripheral white blood cells from beef heifers with varying reproductive potential

Reproductive performance is the most critical factor affecting production efficiency in the cow-calf industry. Heifers with low reproductive efficiency may fail to become pregnant during the breeding season or maintain a pregnancy. The cause of reproductive failure often remains unknown, and the non-pregnant heifers are not identified until several weeks after the breeding season. Therefore, improving heifer fertility utilizing genomic information has become increasingly important. One approach is using microRNAs (miRNA) in the maternal blood that play an important role in regulating the target genes underlying pregnancy success and thereby in selecting reproductively efficient heifers. Therefore, the current study hypothesized that miRNA expression profiles from peripheral white blood cells (PWBC) at weaning could predict the future reproductive outcome of beef heifers. To this end, we measured the miRNA profiles using small RNA-sequencing in Angus-Simmental crossbred heifers sampled at weaning and retrospectively classified as fertile (FH, n = 7) or subfertile (SFH, n = 7). In addition to differentially expressed miRNAs (DEMIs), their target genes were predicted from TargetScan. The PWBC gene expression from the same heifers were retrieved and co-expression networks were constructed between DEMIs and their target genes. We identified 16 differentially expressed miRNAs between the groups (p-value ≤0.05 and absolute (log2 fold change ≥0.05)). Interestingly, based on a strong negative correlation identified from miRNA-gene network analysis with PCIT (partial correlation and information theory), we identified miRNA-target genes in the SFH group. Additionally, TargetScan predictions and differential expression analysis identified bta-miR-1839 with ESR1, bta-miR-92b with KLF4 and KAT2B, bta-miR-2419-5p with LILRA4, bta-miR-1260b with UBE2E1, SKAP2 and CLEC4D, and bta-let-7a-5p with GATM, MXD1 as miRNA-gene targets. The miRNA-target gene pairs in the FH group are over-represented for MAPK, ErbB, HIF-1, FoxO, p53, mTOR, T-cell receptor, insulin and GnRH signaling pathways, while those in the SFH group include cell cycle, p53 signaling pathway and apoptosis. Some miRNAs, miRNA-target genes and regulated pathways identified in this study have a potential role in fertility; other targets are identified as novel and need to be validated in a bigger cohort that could help to predict the future reproductive outcomes of beef heifers.

1 Introduction

Reproductive traits are essential for sustainable food production. Low reproductive capacity in beef heifers is inferred as a failure to become pregnant during the breeding season or maintain a pregnancy to calving (Lamb, 2013). The possible causes for reproductive failure in heifers may be attributed to poor reproductive and nutritional management, diseases, and genetics (Houghton et al., 1990; BonDurant, 2007; Larson and White, 2016; Shao et al., 2021). To overcome this problem, producers often develop more heifers than required for replacements and perform additional inseminations or hormonal treatments that consequently alter the current and subsequent cow-calf production ratio (Shao et al., 2021). In addition to the cost of artificial insemination (AI) and treatments, losses are further increased due to the rearing of the animal from birth to breeding age and then culling and replacing the animals with poor reproductive performance. Management practices to select heifers with high reproductive potential contribute to increasing efficiency; however, some heifers still fail to conceive. Considering this, the beef production sector is particularly interested in identifying biomarkers to predict reproductive efficiency.

In recent years, developments in molecular biology have provided new insights into the potential candidates and biomarkers for fertility and predicting reproductive outcomes. For instance, Binelli et al. (2015) identified transcriptome signatures in uterine biopsies from pregnant cows at 6 days post-AI in beef cattle. Forde et al. (2012) identified increased gene expression in the endometrium of pregnant cattle. Machine learning approaches have provided opportunities to unravel genomic signatures for fertility from omics data. These approaches are powerful in processing and modeling omics data with vast and diverse volumes (Li et al., 2022). Among the studies, Rabaglino and Kadarmideen, (2020) identified 50 genes from endometrial transcriptome as predictors of uterine receptivity to embryo transfer in cattle. Likewise, Diniz et al. (2022) reported nine potential candidates using a multi-tiered approach of machine learning and gene co-expression network on the transcriptome profile of uterine luminal epithelial cells. Although these studies provided insights into the different expression levels of potential candidates regulating fertility or reproductive outcome, the downstream molecular targets regulating the expression of these genes were not investigated. They were also unable to predict reproductive potential at weaning as they were conducted at or near the time of insemination.

Among the factors regulating gene expression, a class of RNA molecules, microRNAs (miRNAs), are known to regulate gene expression. Such regulation is mainly based on binding with messenger RNA (mRNA) targets and destabilizing them, thereby repressing protein production and translational silencing (Cannell et al., 2008). A plethora of studies identified miRNAs as biomarkers for pregnancy outcomes in humans (Barchitta et al., 2017; Ali et al., 2021; Vashukova et al., 2021; Xu et al., 2021) and cattle (Ioannidis and Donadeu, 2016; De Bem et al., 2017; Pohler et al., 2017; Gebremedhn et al., 2018; Lim H.-J. et al., 2021). Ioannidis and Donadeu reported an increased level of circulating miR-26a on days 16–24 of pregnancy in cattle (Ioannidis and Donadeu, 2016). Similarly, levels of circulating miR-221 and miR320a were increased in weeks 8, 12 and 16 of pregnancy in cattle (Lim H.-J. et al., 2021).

Collectively, these findings provide critical insights into the biological mechanisms determining different reproductive outcomes. However, most of the studies have focused on time periods in proximity to the breeding season. Considering the producer has invested time and resources in developing the replacement heifers, using these potential biomarkers at breeding is too late. Thus, it would be beneficial to trace back the potential biomarkers through weaning—when the selection of replacement animals usually takes place. Therefore, our main objective was to profile the miRNA levels in the peripheral white blood cells (PWBC) of beef heifers at weaning that could potentially be used to predict heifers with a high reproductive outcome following AI. We hypothesized that the miRNAs are differentially expressed, the genes regulated by these miRNAs are co-expressed, and the miRNA-gene networks are rewired in beef heifers at an early development stage, contributing to a varying reproductive outcome at maturity.

2 Materials and methods

2.1 Animal handling, sample collection, phenotyping, and classification

All procedures involving animals were approved by the Institutional Animal Care and Use Committee (IACUC) at Auburn University (IACUC protocol numbers 2015-2786 and 2019-3591). Cross-bred heifers (Angus-Simmental, n = 75) utilized in this study originated from and were housed at an outlying Alabama agricultural experiment station research and extension center (Auburn University). The heifers were weaned ∼238 days after birth. Phenotypic data, such as weaning age and weight, were recorded for each heifer. Immediately after weaning, 10 mL of blood was drawn into ethylenediaminetetraacetic acid (EDTA) vacutainers (Becton, Dickinson and Company, Franklin, NJ) from the jugular vein. The vacutainers were inverted 8–10 times and immersed in ice until immediate processing. In the laboratory, blood was processed to isolate peripheral white blood cells (PWBC), as described elsewhere (Banerjee et al., 2023). The PWBCs were stored at −80°C until further processing.

For breeding, 72 heifers were selected based on ideal body condition scores (5–6) and reproductive tract scores (≥4). The breeding protocol, estrus synchronization and fixed-time artificial insemination (AI) have been described previously (Banerjee et al., 2023). In brief, approximately 45 days before breeding by AI, pre-breeding examinations for each heifer were performed to evaluate the pubertal status. All heifers underwent an estrus synchronization and fixed-time AI protocol (7-Day CO-Synch + CIDR) (Dickinson et al., 2018). Fourteen days following AI, all heifers were exposed to fertile bulls for a 60-day natural service breeding season to ensure adequate opportunities for conception to occur. Pregnancy of the heifers was evaluated on day 75 post-AI by transrectal palpation and ultrasound. The presence and gestation age or absence of conceptus was confirmed 75 days following AI with ultrasound and transrectal palpation and was used to classify heifers as fertile to those that were pregnant to AI (FH, n = 35), pregnant to natural breeding (P-NB, n = 26) or subfertile that were not pregnant (SFH, n = 11). Heifers that were pregnant to AI (FH, n = 7) and subfertile (SFH, n = 7) and with comparable birth age, weaning age and body weight in both groups were selected for this study. A schematic representation of the study design and analysis steps is given in Figure 1.

FIGURE 1

FIGURE 1. Schematic representation of the study design and analysis steps.

2.2 RNA extraction from PWBC

Based on the pregnancy classification of heifers, total RNA was extracted from the 14 PWBC samples (FH and SFH) that were collected at the time of weaning. The RNA was extracted using the Trizol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s protocol. Additional RNA purification and DNase digestion steps using an RNA clean and concentrator kit (Zymo Research, Irvine, CA, United States) were included in the protocol. RNA was quantified on a Qubit fluorometer using a Qubit RNA broad-range assay kit (Life Technologies, Thermo Fisher Scientific Inc., MA, United States). The RNA integrity was assessed using the Agilent Bioanalyzer and the Agilent RNA 6000 Nano kit (Agilent, Santa Clara, CA, United States). The quality of small RNAs was determined using the Agilent 2100 Bioanalyzer Small RNA kit (Agilent, Santa Clara, CA, United States). The samples with average RNA integrity number (RIN) values >6.8 were further processed for small RNA library construction.

2.3 Library preparation and sequencing

The total RNA of each sample was diluted with RNase-free water to obtain a final concentration of 1 μg as a starting material. This diluted total RNA was used to prepare libraries using the protocols by NEXTflex small RNA-Seq kit v3 (Perkin Elmer). Following the protocol from the kit, the 5′ and 3’ adapters were ligated to the RNA fragments, which were then reverse-transcribed and amplified (18 cycles) to generate cDNA libraries. Each cDNA library was prepared using a different barcode primer for the ease of being multiplexed for sequencing. Libraries were cleaned using NEXTflex Cleanup beads (gel-free protocol), and the size distribution of the final library was assessed by Agilent Bioanalyzer high-sensitivity DNA assay (Agilent, Santa Clara, CA, United States). The quality-checked libraries were pooled and sequenced in the NextSeq 500 using the single-end 50 bp chemistry at Discovery Life Sciences (Hudson Alpha Institute of Biotechnology, Huntsville, AL, United States).

2.4 Data processing and miRNA expression profiling

Raw sequence demultiplexed reads obtained from Hudson Alpha were initially assessed for sequencing quality using FastQC v0.11.9 (Andrew, 2010) and MultiQC v1.12 (Ewels et al., 2016). The reads were evaluated for quality based on average read length, adaptor content, per-sequence GC content, and sequence quality scores. After that, the 3′ adapter sequence was trimmed using Cutadapt (Martin, 2011) with the following parameters: -a TGG​AAT​TCT​CGG​GTG​CCA​AGG -minimum length 23. The reads were further trimmed using Cutadapt (Martin, 2011) to remove four bases from either side of each read following the small RNA trimming instructions (cutadapt -u 4 -u -4) (recommended by NEXTflex small RNA-Seq kit v3 (Perkin Elmer)). Trimmed Fastq files were checked for quality control with FastQC v0.11.9 and MultiQC v1.12. To profile both novel and known miRNA expression in the samples from the cleaned sequence data, the trimmed reads were processed using the miRDeep2 analysis workflow (Friedländer et al., 2012). Sequences were aligned to Ensembl’s ARS-UCD 1.2 Bos taurus reference genome (https://useast.ensembl.org/Bos_taurus/Info/Index) using the mapper.pl module in miRDeep2 and were further aligned with B. taurus precursor and matured miRNAs extracted from miRBase v22.1 (Kozomara et al., 2019). Only the alignments with zero mismatches in the seed region and those greater than 18 nt of a read mapped to the genome were retained, producing read counts for each sample.

2.5 Differential expression analysis

Initially, the mature miRNAs with zero counts were filtered out. Thereafter, the read counts were transformed to counts per million (CPM) using edgeR v3.28.1 (Robinson et al., 2010). For quality control, raw counts with CPM <1 in 50% of the samples were filtered out. To identify the differentially expressed miRNAs, the raw counts were analyzed in DESeq2 v1.26.0 (Love et al., 2014). The differentially expressed miRNAs (DEMIs) identified with a p-value ≤0.05 and absolute (log2 fold change) ≥ 0.5 were considered significant. The DEMIs were classified as up or downregulated based on the sign of log2 fold change in the SFH group. The up and downregulated miRNAs were visualized using a volcano plot constructed using the R-package EnhancedVolcano v1.4.0 (Blighe et al., 2018).

2.6 Target gene prediction

The target genes for DEMIs were predicted using TargetScan v8.0, which predicts the targets of miRNAs by searching the conserved 8mer, 7mer and 6mer sites that match the seed region of each miRNA (Lewis et al., 2005). The search for mRNA targets was specific to B. taurus. As a feature of TargetScan, the predictions were ranked based on targeting efficacy estimated using cumulative weighted context++ scores (Agarwal et al., 2015). The genes with high confidence (≤−0.4) in the cumulative weighted context++ score were selected as target genes for each miRNA (Agarwal et al., 2015).

2.7 mRNA-miRNA co-expression and network analysis

To identify the correlation pattern between the DEMI-gene targets identified by TargetScan, we also generated PWBC mRNA expression profiles from the same heifers as reported elsewhere (Banerjee et al., 2023). In brief, the mRNA profile of the same set of samples from FH (n = 7) and SFH (n = 7) was sequenced to generate paired-end 100 bp reads. After a quality check using FastQC v0.11.9 and MultiQC v1.12, the raw sequences were mapped using STAR aligner v2.7.5 (Dobin et al., 2013) to Ensemble’s ARS-UCD1.2 B. taurus genome as the reference. The raw counts per gene (obtained using the -quantMode function in STAR) were transformed to CPM using edgeR v3.28.1. The genes identified from the PWBC mRNA expression profiles were filtered to retain only those identified as miRNA targets based on TargetScan and were named ‘PWBC expressed TargetScan genes’. To determine the miRNA-gene correlation, we used the partial correlation and information theory (PCIT) approach. The CPM normalized ‘PWBC expressed TargetScan genes’ and DEMIs were used as an input for PCIT. The PCIT reports the significantly correlated pairs after comparing the possible trios of genes (Reverter and Chan, 2008), which allows us to build co-expression networks.

The networks were constructed using the correlations identified by PCIT between DEMIs and ‘PWBC-expressed TargetScan genes’ for each FH and SFH group in Cytoscape v3.8.2 (Shannon et al., 2003). The network for each group was analyzed using the Network Analyzer tool in Cytoscape v3.8.2 (Assenov et al., 2008). Next, to determine the differentially connected genes and miRNAs in each group, the network connectivity (K) measure for each network was standardized by taking a ratio of gene connectivity (degree) and maximum connectivity in each network. The differential connectivity (DK) measure was calculated as $\mathit{D}{\mathit{K}}_{\mathit{i}}=\mathit{K}\mathit{S}\mathit{F}{\mathit{H}}_{\left(\mathit{i}\right)}-\mathit{K}\mathit{F}{\mathit{H}}_{\left(\mathit{i}\right)}$. The DK values were transformed to a z-score and ±1.96 SD was considered significant (p ≤ 0.05). The connectivity gain or loss was evaluated for each gene or miRNA in the SFH group as the reference. The central reference network comparing both groups was constructed using DyNet (Cytoscape plug-in) (Goenawan et al., 2016). DyNet allows for identifying and visualizing nodes and edge rewiring in response to cellular signals and highlights the most rewired nodes based on the central reference network constructed with FH and SFH groups (Goenawan et al., 2016).

2.8 Minimum free energy for miRNA-target prediction

To cross-validate the specificity of the target genes identified from TargetScan and correlate to the DEMIs, we used the RNAhybrid (Rehmsmeier et al., 2004). To this end, we retrieved the mature miRNA sequence from the miRBase database (Kozomara et al., 2019) and the cDNA sequence of the genes from BioMart (Durinck et al., 2005). The minimum free energy (mfe) was predicted for the miRNA-mRNA sequence; several miRNAs predicted to have a low mfe indicate a relatively high affinity for miRNA-mRNA complexes. The cutoff for mfe was set to less than −20 kcal/mol for evaluating miRNA-mRNA affinity (Yen et al., 2019; de Lima et al., 2021). Average mfe was calculated for the genes with multiple transcripts.

2.9 Functional enrichment analysis

The pathway over-representation analysis for the genes correlated with miRNAs identified by PCIT was carried out using ClueGO v2.5.8 (Cytoscape plug-in) (Bindea et al., 2009). The redundant terms were clustered based on a kappa score of 0.4, and the B. taurus annotation was used as the background in ClueGO. The significant pathways were selected based on the corrected p-value with Bonferroni step-down feature ≤0.05.

3 Results

3.1 Data generation and mapping statistics

We used small RNA-Seq to identify the difference in the expression of miRNAs from PWBC between two groups, FH and SFH. The sequencing from all the samples yielded an average of 6.9 million reads. On average, 78.13% of reads were uniquely mapped to the B. taurus reference genome (Supplementary Table S1), yielding a dataset of 1,207 miRNAs in 14 samples. Following filtering to remove the zero counts, 614 miRNAs were retained. After the CPM filter approach, 341 miRNAs in 14 samples were used for further analysis.

3.2 Differentially expressed miRNAs (DEMIs)

We identified 16 DEMIs between the FH and SFH group out of 341 miRNAs [(p-value ≤0.05 and absolute (log2 fold change ≥0.5)] (Figure 2; Supplementary Table S2). The top 5 upregulated miRNAs in the SFH group with the highest fold change difference included bta-miR-677, bta-miR-1839, bta-miR-1434-3p, bta-miR-2332 and bta-miR-140, while the downregulated miRNAs were bta-miR-450b, bta-miR-2419-5p, bta-miR-92b, bta-miR-574 and bta-miR-2478.

FIGURE 2

FIGURE 2. Volcano plot of differentially expressed miRNAs between FH and SFH groups. Each dot corresponds to a miRNA. The log2fold change is represented in the x-axis represents, while the negative log (base 10) of the p-value in the y-axis. The horizontal dashed line represents the threshold with a p-value cutoff <0.05, while the vertical bars represent the absolute log2fold change >0.5. The 16 DEMIs are labeled in the plot. The left panel (0 to −3 of log2fold change) is downregulated, while the right panel (0–4 of log2fold change) is upregulated miRNAs.

3.3 Target gene prediction and overlapping with genes expressed in PWBC

The predicted gene targets for 16 DEMIs were retrieved from TargetScan. The genes with high confidence (≤−0.4) in the cumulative weighted context++ score were selected as target genes for each miRNA. We identified 1,365 genes targeted by the 16 DEMIs (Figure 3), out of which 741 were expressed in PWBC (Supplementary Table S3). Out of the 741 genes, nine were previously reported as differentially expressed in PWBC from the same set of FH and SFH heifers (Supplementary Table S3) (Banerjee et al., 2023). The nine genes and their miRNA targets are as follows: bta-miR-92b targeted KLF4 and KAT2B, bta-miR-2419-5p targeted LILRA4, bta-miR-1260b targeted UBE2E1, SKAP2 and CLEC4D, bta-let-7a-5p targeted GATM, MXD1 and bta-miR-1839 targeted ESR1.

FIGURE 3

FIGURE 3. Number of genes targeted by 16 DEMIs (ranked in descending order).

3.4 miRNA-mRNA co-expression and network analysis

We used the PCIT algorithm to build two networks (each from FH and SFH groups) to identify the coordinated expression patterns and potential regulatory role between genes and miRNAs. For the network, the correlation of CPM normalized PWBC-expressed TargetScan genes (n = 741) (see Methods) and DEMI (n = 16), corresponding to 7 samples in FH and 6 in the SFH group, was used. One sample was removed from the SFH group due to the poor quality of the mRNA sequence. We identified 33,014 and 29,280 correlations in FH and SFH groups, respectively, consisting of miRNA-gene, miRNA-miRNA, and gene-gene correlations (Supplementary Table S4). To retrieve the biologically meaningful miRNAs and reduce data dimensionality, we only retrieved the miRNA-miRNA and miRNA-gene pairs (Supplementary Table S5). We identified 926 and 689 significantly correlated pairs (r > |0.6|, p ≤ 0.05) in the FH and SFH groups (corresponding to 546 and 409 unique genes) (Supplementary Table S5). The network for each group was visualized using Cytoscape (Figures 4A, B). Of the 16 DEMIs, 15 miRNAs were differentially connected with only bta-miR-1260b not having differential connections. Overall, the SFH group exhibited a loss of miRNA-mRNA network connectivity (Supplementary Table S6). Next, we identified the correlation pattern of the miRNA-gene pair in our study that was identified from TargetScan (Supplementary Table S7). Of the 741 predicted target genes with the corresponding miRNA, 48 and 39 genes were correlated with 12 and 11 DEMIs in the FH and SFH groups, respectively (Supplementary Table S7; Figures 5A, B). Out of these pairs, we identified 30 and 27 negative correlations among the miRNA-gene pairs in FH and SFH, respectively (Supplementary Table S7). Some correlated pairs included bta-miR-2478 with SMIM7, bta-miR-574 with TRADD, HDHD2 and HSF1, and bta-miR-2332 with MAX in the FH group; and bta-miR-1260b with ERAS, bta-miR-2332 with SEC63 and HIF1A in the SFH group, exhibited a high negative correlation (r = 0.9). All the miRNA-gene pairs exhibited a mfe <−20 kcal/mol, except bta-miR-450b with SH3YL1 (−18.9 kcal/mol) and with B2M (−19.7 kcal/mol).

FIGURE 4

FIGURE 4. Regulatory networks of co-expressed genes with 16 DEMIs in (A) FH and (B) SFH groups. Nodes are genes significantly correlated with miRNAs, while edges are positive or negative interactions (correlations) between the miRNA and target genes. The blue diamond represents the DEMIs; green strokes represent positive correlations, while red strokes represent negative correlations.

FIGURE 5

FIGURE 5. Central reference network constructed using DyNet. (A) Network comparison based on the rewired node in the FH and SFH group. The network comprises of 650 nodes and 1547 edges. The blue diamond represents the DEMIs. Unique nodes are shown in red (FH) and green (SFH). Shared nodes are shown in white. (B) Central reference network showing the miRNA-gene correlated pair as identified from TargetScan (Supplementary Table S7). Unique nodes are shown in red (FH) and green (SFH). Shared nodes are shown in white. For ease of visualization, each miRNA-target gene pair is labeled with the same color.

3.5 Pathway analysis

To translate the genes targeted by the miRNAs from PCIT into biologically meaningful processes, we performed pathway enrichment analysis. The top significant (Bonferroni corrected Group p-value ≤0.05) KEGG pathways identified in the FH group include MAPK, ErbB, HIF-1, FoxO, p53, mTOR, T-cell receptor, insulin and GnRH signaling, apoptosis and pathways regulating pluripotency of stem cells (Supplementary Table S8). The top significant KEGG pathways of the SFH group include cell cycle, p53 signaling pathway and apoptosis (Supplementary Table S9). Additionally, the biological process in the SFH group includes a hormone-mediated signaling pathway, cellular response to stress, negative regulation of growth, inflammatory response to antigenic stimulus and signal transduction in response to DNA damage (Supplementary Table S9).

4 Discussion

Evidence suggests that miRNAs play an essential regulatory role in several biological processes involving cell proliferation, cell death, epigenetic changes, and apoptosis (Bartel, 2004; Cupp et al., 2009). All these processes have the potential to promote phenotypic variation within the population. Among the traits, phenotypes associated with reproductive outcome and fertility are of major interest to cattle producers, mainly because of the high production loss from a heifer of low reproductive potential. Specifically in reproductive biology, miRNAs have been demonstrated to be important regulators of embryonic development in humans and livestock, such as predicting pregnancy outcomes (Ioannidis and Donadeu, 2016; Hitit et al., 2022), embryo viability (Pohler et al., 2017) and implantation (Kose et al., 2022), and endometrial receptivity (Zang et al., 2021). Moreover, miRNAs regulate the gene targets by repressing or stimulating their expression (Orang et al., 2014). Thus, the present study was undertaken to identify the difference in the miRNA expression and the respective gene targets in the PWBC of beef heifers at weaning and their potential as predictors of future reproductive potential in beef heifers.

We identified five miRNAs (bta-miR-92b, bta-miR-2419-5p, bta-miR-1260b, bta-miR-1839 and bta-let-7a-5p) targeting differentially expressed genes from PWBC of the same beef heifers. Among them, bta-miR-92b is the most significant miRNA downregulated in the SFH group. Expression of miR-92b was significantly upregulated in the porcine placenta on day 90 of gestation (Su et al., 2010). In a study reported in cattle, miR-92b (downregulated expression) was associated with endometritis, an inflammatory response in the endometrium that causes reproductive disorder (Jiang et al., 2021). The authors further reported that over-expression of miR-92b significantly suppressed the activation of the toll-like receptor in lipopolysaccharide-mediated bovine endometrial cells, thereby reducing pro-inflammatory cytokines and inhibiting cell apoptosis (Jiang et al., 2021). Interestingly, in a separate study, miR-92 has been identified to regulate KLF4 (Hamik and Jain, 2012). In female fertility, KLF4 mediates the anti-proliferative effects of progesterone during the G0/G1 arrest in endometrial epithelial cells in humans (Shimizu et al., 2010). In our study, KLF4 was identified as one of the potential targets of bta-miR-92b from TargetScan and associated with pathways involved with post-embryonic development, cell-cell adhesion, and immune system development in the FH group. Another target gene associated with bta-miR-92b in our study was KAT2B. With the pathway analysis, we identified KAT2B underlying negative regulation of serine/threonine kinase activity and phosphorylation pathways in the FH group. Interestingly, KLF4 and KAT2B have been identified as downregulated in the PWBC of non-pregnant heifers compared to pregnant heifers (Banerjee et al., 2023). Furthermore, from the PCIT network analysis, we identified bta-miR-92b as more connected with genes/miRNA targets in fertile heifers (68 targets) than in subfertile heifers (31 targets). Among the common genes from TargetScan and PCIT, we found that the bta-miR-92b was correlated with the BTG2, REV3L, FOXN2 and NSMF genes in the FH group and DUSP10 in the SFH group. DUSP10 is associated with negative regulation of growth, phosphorylation, and transferase and kinase activity. DUSP10 plays a significant role in the innate and adaptive immune response by regulating the mitogen-activated protein (MAP) kinase modulated by kinase and phosphatase (Zhang et al., 2004). Despite these findings, further mechanisms of bta-miR-92b and its gene targets need to be explored for fertility and reproductive outcome in bovines.

Bta-miR-1260b was significantly downregulated in the SFH group. In previous studies, miR-1260 has been associated with infertility in humans (Butler et al., 2021). From TargetScan, we identified UBE2E1, SKAP2 and CLEC4D as predicted genes targeted by bta-miR-1260b. These genes were differentially expressed in the PWBC of beef heifers at weaning (Banerjee et al., 2023). SKAP2 and CLEC4D were associated with regulating the immune response pathway in the FH group. SKAP2 is a substrate of Src family kinase and regulates cellular processes, including proliferation, adhesion, migration and stress response (He et al., 2017). In a previous study, SKAP2 was detected in all developmental stages of mouse oocytes and depletion of SKAP2 caused the failure of spindle migration, polar body extrusion and cytokinesis defects (He et al., 2017). In a separate study in sheep, SKAP2 was differentially expressed in the granulosa cells from superstimulated lamb and ewe follicles and was associated with cellular growth, proliferation, and migration (Wu et al., 2016). CLEC4D was found to play important roles in immunity and homeostasis (Graham et al., 2012). Besides SKAP2 and CLEC4D, bta-miR-1260b had a high correlation (r > |0.9|) with ERAS and NSG1 in the SFH group, identified from the network analysis. Furthermore, ERAS was identified with a mfe of −30.5 kcal/mol confirming as a predictive target of bta-miR-1260b. In bovines, lipopolysaccharide treatment of bovine endometrial epithelial cells caused differential methylation of the NSG1 gene associated with inflammation and endometrial function (Jhamat et al., 2020). These studies suggest the role of the genes targeted by bta-miR-1260b with inflammation and immune response. Inflammation and immune response are connected at many levels with fertility and reproductive outcome in cattle (Fair, 2015); however, the detailed mechanisms still need to be explored.

Bta-miR-1839 was upregulated but had less miRNA-gene network connections (64 targets) in the SFH group compared to the FH group (69 targets). In a microarray-based mice study, an increase in the expression of miR-1839 was observed during the implantation periods (day 5) in the luminal epithelium and endometrium (Li et al., 2015). We identified ESR1 (estrogen receptor 1) targeted by the bta-miR-1839 through TargetScan. Interestingly, ESR1 was downregulated in the non-pregnant beef heifers compared to the pregnant group (Banerjee et al., 2023). In humans, ESR1 gene expression and protein abundance were disrupted in the endometrium of patients with severe preeclampsia (Garrido-Gomez et al., 2021). Likewise, the role of ESR1 with the cAMP signaling pathway was identified as critical in decidualization (Kaya Okur et al., 2016). Moreover, an association of genetic variants in ESR1 has been identified with recurrent pregnancy loss in women (Pan et al., 2014; Bahia et al., 2020). Previous studies reported ESR1 mediating the biological activity and proliferative effects of estrogen on the reproductive tissues, including ovarian follicular cells (Hewitt and Korach, 2003). Moreover, ERα (ESR1) knockout female mice were sterile and preovulatory follicles did not ovulate upon superovulation treatment (Hewitt and Korach, 2003). In our study, ESR1 was associated with the activation of immune response and pattern recognition receptor signaling pathway, further supported by previous human studies (Kovats, 2015).

The SFH group heifers exhibited bta-miR-2419-5p as significantly downregulated and with low miRNA-gene connectivity (29 targets) in our study and associated with LILRA4 identified by TargetScan. In a study on super stimulatory response in cattle, plasma miR-2419-5p was identified as downregulated in unstimulated low ovarian response heifers compared to high response heifers (Gad et al., 2020). LILRA4 is a member of leukocyte immunoglobulin-like receptors (LILR) that regulate innate and adaptive immune functions (Hogan et al., 2012). LILRA4 was significantly downregulated in the subfertile beef heifers at weaning compared to the heifers conceived to AI during the first breeding season (Banerjee et al., 2023). This was supported by the findings in a study on blood transcriptome in Holstein cows where LILRA4 was downregulated in the miscarriage cow group compared to the pregnant group (Zhao et al., 2019).

The importance of miRNA function is not only for regulating adaptive and innate immune response but also for cellular proliferation, trophoblast invasion and cellular differentiation (Hayder et al., 2018), which are important for a successful pregnancy. We identified bta-let-7a-5p as downregulated in the SFH group heifers. Expression of let-7a-5p in the placenta during the first trimester has been involved with low cytotrophoblast proliferation (Smith et al., 2021). In a similar context in humans, let-7a was downregulated in the blood plasma of females with recurrent pregnancy loss compared to controls (Jairajpuri et al., 2021). In a transcriptome study with human plasma and placenta, let-7a-5p was upregulated at 11–23 weeks of gestation compared to the 6–10 weeks group (Smith et al., 2021). For further insight, we identified bta-let-7a-5p with more miRNA-gene connectivity in the SFH group (33 targets) than in the FH group (29 targets). Furthermore, GATM and MXD1 genes were regulated by bta-let-7a-5p identified by TargetScan. GATM, a gene that encodes L-arginine and catalyzes the rate-limiting step in the synthesis of creatinine, was associated with the amino acid biosynthesis pathway in our study. This is further supported by a previous study where GATM is expressed during development and is imprinted in mouse placenta and yolk-sac (Sandell et al., 2003).

Apart from these miRNA-gene pairs, we identified bta-miR-574 with a high negative correlation (>0.85) with FKBP4 and N4BP1 in the SFH group and with TRADD, HDHD2 and HSF1 in the FH group. Bta-miR-574 exhibited a loss of gene connectivity in the SFH group. Interestingly, miR-574 was identified to be upregulated in women with preeclampsia and involved with endothelial dysfunction (Lip et al., 2018). Furthermore, miR-574 in porcine cumulus cells has been shown to suppress oocyte maturation (Pan et al., 2018). In our study, bta-miR-574 was correlated with the FKBP4 gene that was associated with response to lipid and steroid hormone and female pregnancy pathways. In humans, FKBP4 mRNA expression was decreased in the endometrium with endometriosis compared to controls (Yang et al., 2012). In fertile heifers, bta-miR-574 was correlated with TRADD, which was over-represented for apoptosis pathways. Furthermore, bta-miR-574 was correlated with HSF1, which was over-represented for nucleic acid transport. In previous studies, TRADD has been reported as regulating perinatal development and adulthood survival in mice (Dowling et al., 2019), while HSF1 was essential for reproductive success in pre-implantation embryos (Abane and Mezger, 2010). MiRNA bta-miR-2332 was negatively correlated with MAX in the FH group and with HIF1A, SEC63 and COX7B in the SFH group. HIF1A, in response to gonadotropins, activates steroidogenesis and cell proliferation in granulosa cells critical for ovulation (Lim M. et al., 2021), thereby exhibiting a potential role in fertility. We identified bta-miR-2478 downregulated with a loss of gene connectivity in the SFH group. Supporting our results, previous studies found that the expression of bta-miR-2478 was downregulated in bovine cumulus cells that did not cleave compared to ones that cleaved to form a blastocyst, thereby serving as an indicator of oocyte quality (Uhde et al., 2017). These miRNAs and miRNA-gene targets open up further avenues for research to be explored in relation to fertility and reproductive outcome in heifers.

From the network connectivity, we identified bta-miR-92b, bta-miR-574, bta-miR-1839, bta-miR-450b, bta-miR-2419-5p, bta-miR-874, bta-miR-2478 and bta-miR-1434-3p with a loss in network connectivity with the target genes and miRNAs, while bta-miR-2332, bta-miR-677, bta-let-7b, bta-let-7a-5p, bta-miR-1306, bta-miR-140 and bta-miR-6119-3p with a gain in network connections in the subfertile heifer group. This suggests that the loss of connectivity in heifers with low reproductive outcome is due to the rewiring of the major regulators. A potential explanation for the heifers not becoming pregnant could be attributed to the difference in the connectivity between the miRNAs and their target genes.

Some miRNAs and their targeted genes were previously reported to be associated with fertility and pregnancy outcomes in cattle and other species; herein, we reported novel miRNA-gene pairs that warrant further investigation. Combining the gene prediction of miRNAs from TargetScan and the genes differentially expressed in our previous study, we identified bta-miR-92b targeted KLF4, KAT2B, bta-miR-1839 targeted ESR1, bta-let-7a-5p targeted GATM and MXD1, bta-miR-1260b targeted UBE2E1, SKAP2 and bta-miR-2419-5p targeted LILRA4. Additionally, with the PCIT network analysis approach, we identified genes with high negative correlation targeted by the miRNAs, such as DUSP10 correlated with bta-miR-92b, ERAS with bta-miR-1260b, FKBP4 with bta-miR-574 and HIF1A with bta-miR-2332 in the SFH group. Most of the previously reported studies revealed a difference in the miRNA and target gene levels during pregnancy; however, there is a possibility that these miRNAs and genes have been dysregulated for a long time without any visible symptoms until identified with pregnancy complications. Therefore, our study provides insights into the differential expression of miRNAs and genes in beef heifers as early as weaning that might play a role in predicting future reproductive performance. Thorough knowledge of the interactions will likely be needed to improve our understanding of female fertility issues and potentially develop therapeutic targets. One of the limitations of this study was the sample size, i.e., 14 heifers (7 in each group). Therefore, confirming these targets in a bigger cohort and following-up with additional time points including post-weaning, before AI and during pregnancy in heifers is required to further validate these targets and determine their potential benefit to the beef cattle industry.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc&equals;GSE225854,GSE225854.

Ethics statement

The animal study was reviewed and approved by the the Institutional Animal Care and Use Committee (IACUC) at Auburn University (IACUC protocol numbers 2015-2786 and 2019-3591).

Author contributions

Conceptualization, PB and PD; methodology, PB, WD, SR, and PD; software, PB; Formal analysis, PB and WD; writing—original draft preparation, PB; writing—review and editing, PB, WD, PD, and SR; funding acquisition, PD. All authors have read and agreed to the published version of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Alabama Cattlemen’s Association, the Foundation for Food and Agriculture Research—grant no. FF-NIA19-0000000048, and by the Alabama Agricultural Experiment Station. The research was also supported by USDA—Agricultural Research Service—USDA-58-6010-1-005.

Acknowledgments

This project was financially supported by the Agricultural Research Service, U.S. Department of Agriculture, under Agreement No. 58-6010-1-005. This work used the resources of the High-Performance Computer Center—Easley Cluster from Auburn University.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1174145/full#supplementary-material

References

Abane, R., and Mezger, V. (2010). Roles of heat shock factors in gametogenesis and development. FEBS J. 277, 4150–4172. doi:10.1111/j.1742-4658.2010.07830.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Agarwal, V., Bell, G. W., Nam, J.-W., and Bartel, D. P. (2015). Predicting effective microRNA target sites in mammalian mRNAs. Elife 4, e05005. doi:10.7554/eLife.05005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ali, A., Hadlich, F., Abbas, M. W., Iqbal, M. A., Tesfaye, D., Bouma, G. J., et al. (2021). MicroRNA–mRNA networks in pregnancy complications: A comprehensive downstream analysis of potential biomarkers. Int. J. Mol. Sci. 22, 2313. doi:10.3390/ijms22052313

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrew, S. (2010). FastQC: A quality control tool for high throughput sequence data. Available at: https://www.bioinformatics.babraham.ac.uk/index.html.

Google Scholar

Assenov, Y., Ramírez, F., Schelhorn, S.-E., Lengauer, T., and Albrecht, M. (2008). Computing topological parameters of biological networks. Bioinformatics 24, 282–284. doi:10.1093/bioinformatics/btm554

PubMed Abstract | CrossRef Full Text | Google Scholar

Bahia, W., Soltani, I., Haddad, A., Soua, A., Radhouani, A., Mahdhi, A., et al. (2020). Association of genetic variants in estrogen receptor (ESR)1 and ESR2 with susceptibility to recurrent pregnancy loss in Tunisian women: A case control study. Gene 736, 144406. doi:10.1016/j.gene.2020.144406

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, P., Diniz, W. J. S., Hollingsworth, R., Rodning, S. P., and Dyce, P. W. (2023). mRNA signatures in peripheral white blood cells predict reproductive potential in beef heifers at weaning. Genes (Basel) 14, 498. doi:10.3390/genes14020498

PubMed Abstract | CrossRef Full Text | Google Scholar

Barchitta, M., Maugeri, A., Quattrocchi, A., Agrifoglio, O., and Agodi, A. (2017). The role of miRNAs as biomarkers for pregnancy outcomes: A comprehensive review. Int. J. Genomics 2017, 8067972. doi:10.1155/2017/8067972

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2004). MicroRNAs: Genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi:10.1016/S0092-8674(04)00045-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25, 1091–1093. doi:10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

Binelli, M., Scolari, S. C., Pugliesi, G., Van Hoeck, V., Gonella-Diaza, A. M., Andrade, S. C. S., et al. (2015). The transcriptome signature of the receptive bovine uterus determined at early gestation. PLoS One 10, e0122874. doi:10.1371/journal.pone.0122874

PubMed Abstract | CrossRef Full Text | Google Scholar

Blighe, K., Rana, S., and Lewis, M. (2018). EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling. Available at: https://github.com/kevinblighe/EnhancedVolcano.

Google Scholar

BonDurant, R. H. (2007). Selected diseases and conditions associated with bovine conceptus loss in the first trimester. Theriogenology 68, 461–473. doi:10.1016/j.theriogenology.2007.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Butler, A. E., Cunningham, T. K., Ramachandran, V., Diboun, I., Halama, A., Sathyapalan, T., et al. (2021). Association of microRNAs with embryo development and fertilization in women undergoing subfertility treatments: A pilot study. Front. Reprod. Heal. 3, 719326. doi:10.3389/frph.2021.719326

CrossRef Full Text | Google Scholar

Cannell, I. G., Kong, Y. W., and Bushell, M. (2008). How do microRNAs regulate gene expression? Biochem. Soc. Trans. 36, 1224–1231. doi:10.1042/BST0361224

PubMed Abstract | CrossRef Full Text | Google Scholar

Cupp, A. S., Matthews, J., Huff-Lonergan, E., Spurlock, D. M., and McLean, D. (2009). Cell biology symposium: The role of microRNA in cell function. J. Anim. Sci. 87, E19–E20. doi:10.2527/jas.2008-1552

PubMed Abstract | CrossRef Full Text | Google Scholar

De Bem, T. H. C., da Silveira, J. C., Sampaio, R. V., Sangalli, J. R., Oliveira, M. L. F., Ferreira, R. M., et al. (2017). Low levels of exosomal-miRNAs in maternal blood are associated with early pregnancy loss in cloned cattle. Sci. Rep. 7, 14319. doi:10.1038/s41598-017-14616-1

PubMed Abstract | CrossRef Full Text | Google Scholar

de Lima, A. O., Afonso, J., Edson, J., Marcellin, E., Palfreyman, R., Porto-Neto, L. R., et al. (2021). Network analyses predict small RNAs that might modulate gene expression in the testis and epididymis of Bos indicus bulls. Front. Genet. 12, 610116. doi:10.3389/fgene.2021.610116

PubMed Abstract | CrossRef Full Text | Google Scholar

Dickinson, S. E., Griffin, B. A., Elmore, M. F., Kriese-Anderson, L., Elmore, J. B., Dyce, P. W., et al. (2018). Transcriptome profiles in peripheral white blood cells at the time of artificial insemination discriminate beef heifers with different fertility potential. BMC Genomics 19, 129. doi:10.1186/s12864-018-4505-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Diniz, W. J. S., Banerjee, P., Rodning, S. P., and Dyce, P. W. (2022). Machine learning-based Co-expression network analysis unravels potential fertility-related genes in beef cows. Animals 12, 2715. doi:10.3390/ani12192715

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). Star: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi:10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Dowling, J. P., Alsabbagh, M., Del Casale, C., Liu, Z.-G., and Zhang, J. (2019). TRADD regulates perinatal development and adulthood survival in mice lacking RIPK1 and RIPK3. Nat. Commun. 10, 705. doi:10.1038/s41467-019-08584-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Durinck, S., Moreau, Y., Kasprzyk, A., Davis, S., De Moor, B., Brazma, A., et al. (2005). BioMart and bioconductor: A powerful link between biological databases and microarray data analysis. Bioinformatics 21, 3439–3440. doi:10.1093/bioinformatics/bti525

PubMed Abstract | CrossRef Full Text | Google Scholar

Ewels, P., Magnusson, M., Lundin, S., and Käller, M. (2016). MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048. doi:10.1093/bioinformatics/btw354

PubMed Abstract | CrossRef Full Text | Google Scholar

Fair, T. (2015). The contribution of the maternal immune system to the establishment of pregnancy in cattle. Front. Immunol. 6, 7. doi:10.3389/fimmu.2015.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Forde, N., Duffy, G. B., McGettigan, P. A., Browne, J. A., Mehta, J. P., Kelly, A. K., et al. (2012). Evidence for an early endometrial response to pregnancy in cattle: Both dependent upon and independent of interferon tau. Physiol. Genomics 44, 799–810. doi:10.1152/physiolgenomics.00067.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi:10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Gad, A., Sánchez, J. M., Browne, J. A., Nemcova, L., Laurincik, J., Prochazka, R., et al. (2020). Plasma extracellular vesicle miRNAs as potential biomarkers of superstimulatory response in cattle. Sci. Rep. 10, 19130. doi:10.1038/s41598-020-76152-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Garrido-Gomez, T., Castillo-Marco, N., Clemente-Ciscar, M., Cordero, T., Muñoz-Blat, I., Amadoz, A., et al. (2021). Disrupted PGR-B and ESR1 signaling underlies defective decidualization linked to severe preeclampsia. Elife 10, e70753. doi:10.7554/eLife.70753

PubMed Abstract | CrossRef Full Text | Google Scholar

Gebremedhn, S., Salilew-Wondim, D., Hoelker, M., Held-Hoelker, E., Neuhoff, C., Tholen, E., et al. (2018). Exploring maternal serum microRNAs during early pregnancy in cattle. Theriogenology 121, 196–203. doi:10.1016/j.theriogenology.2018.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Goenawan, I. H., Bryan, K., and Lynn, D. J. (2016). DyNet: Visualization and analysis of dynamic molecular interaction networks. Bioinformatics 32, 2713–2715. doi:10.1093/bioinformatics/btw187

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, L. M., Gupta, V., Schafer, G., Reid, D. M., Kimberg, M., Dennehy, K. M., et al. (2012). The C-type lectin receptor CLECSF8 (CLEC4D) is expressed by myeloid cells and triggers cellular activation through syk kinase. J. Biol. Chem. 287, 25964–25974. doi:10.1074/jbc.M112.384164

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamik, A., and Jain, M. K. (2012). MiRrored regulation of KLF2 and KLF4. Arterioscler. Thromb. Vasc. Biol. 32, 839–840. doi:10.1161/ATVBAHA.112.245563

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayder, H., O’Brien, J., Nadeem, U., and Peng, C. (2018). MicroRNAs: Crucial regulators of placental development. Reproduction 155, R259–R271. doi:10.1530/REP-17-0603

PubMed Abstract | CrossRef Full Text | Google Scholar

He, S.-W., Xu, B.-H., Liu, Y., Wang, Y.-L., Chen, M.-H., Xu, L., et al. (2017). SKAP2 regulates Arp2/3 complex for actin-mediated asymmetric cytokinesis by interacting with WAVE2 in mouse oocytes. Cell Cycle 16, 2272–2281. doi:10.1080/15384101.2017.1380126

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewitt, S., and Korach, K. (2003). Oestrogen receptor knockout mice: Roles for oestrogen receptors alpha and beta in reproductive tissues. Reproduction 125, 143–149. doi:10.1530/rep.0.1250143

PubMed Abstract | CrossRef Full Text | Google Scholar

Hitit, M., Kose, M., Kaya, M. S., Kırbas, M., Dursun, S., Alak, I., et al. (2022). Circulating miRNAs in maternal plasma as potential biomarkers of early pregnancy in sheep. Front. Genet. 13, 929477. doi:10.3389/fgene.2022.929477

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, L., Bhuju, S., Jones, D. C., Laing, K., Trowsdale, J., Butcher, P., et al. (2012). Characterisation of bovine leukocyte ig-like receptors. PLoS One 7, e34291. doi:10.1371/journal.pone.0034291

PubMed Abstract | CrossRef Full Text | Google Scholar

Houghton, P. L., Lemenager, R. P., Horstman, L. A., Hendrix, K. S., and Moss, G. E. (1990). Effects of body composition, pre- and postpartum energy level and early weaning on reproductive performance of beef cows and preweaning calf gain. J. Anim. Sci. 68, 1438–1446. doi:10.2527/1990.6851438x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ioannidis, J., and Donadeu, F. X. (2016). Circulating miRNA signatures of early pregnancy in cattle. BMC Genomics 17, 184. doi:10.1186/s12864-016-2529-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jairajpuri, D. S., Malalla, Z. H., Mahmood, N., Khan, F., and Almawi, W. Y. (2021). Differentially expressed circulating microRNAs associated with idiopathic recurrent pregnancy loss. Gene 768, 145334. doi:10.1016/j.gene.2020.145334

PubMed Abstract | CrossRef Full Text | Google Scholar

Jhamat, N., Niazi, A., Guo, Y., Chanrot, M., Ivanova, E., Kelsey, G., et al. (2020). LPS-treatment of bovine endometrial epithelial cells causes differential DNA methylation of genes associated with inflammation and endometrial function. BMC Genomics 21, 385. doi:10.1186/s12864-020-06777-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, K., Yang, J., Song, C., He, F., Yang, L., and Li, X. (2021). Enforced expression of miR-92b blunts E. coli lipopolysaccharide-mediated inflammatory injury by activating the PI3K/AKT/β-catenin pathway via targeting PTEN. Int. J. Biol. Sci. 17, 1289–1301. doi:10.7150/ijbs.56933

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaya Okur, H. S., Das, A., Taylor, R. N., Bagchi, I. C., and Bagchi, M. K. (2016). Roles of estrogen receptor-α and the coactivator MED1 during human endometrial decidualization. Mol. Endocrinol. 30, 302–313. doi:10.1210/me.2015-1274

PubMed Abstract | CrossRef Full Text | Google Scholar

Kose, M., Hitit, M., Kaya, M. S., Kırbas, M., Dursun, S., Alak, I., et al. (2022). Expression pattern of microRNAs in ovine endometrium during the peri-implantation. Theriogenology 191, 35–46. doi:10.1016/j.theriogenology.2022.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kovats, S. (2015). Estrogen receptors regulate innate immune cells and signaling pathways. Cell. Immunol. 294, 63–69. doi:10.1016/j.cellimm.2015.01.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozomara, A., Birgaoanu, M., and Griffiths-Jones, S. (2019). miRBase: from microRNA sequences to function. Nucleic Acids Res. 47, D155–D162. doi:10.1093/nar/gky1141

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamb, G. C. (2013). Criteria for selecting replacements at weaning, before breeding, and after breeding. Vet. Clin. North Am. Food Anim. Pract. 29, 567–578. doi:10.1016/j.cvfa.2013.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Larson, R. L., and White, B. J. (2016). Reproductive systems for north American beef cattle herds. Vet. Clin. North Am. Food Anim. Pract. 32, 249–266. doi:10.1016/j.cvfa.2016.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, B. P., Burge, C. B., and Bartel, D. P. (2005). Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are MicroRNA targets. Cell 120, 15–20. doi:10.1016/j.cell.2004.12.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Li, L., Xu, Y., and Yang, J. (2022). Machine learning meets omics: Applications and perspectives. Brief. Bioinform 23, bbab460. doi:10.1093/bib/bbab460

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Jia, J., Gou, J., Zhao, X., and Yi, T. (2015). MicroRNA-451 plays a role in murine embryo implantation through targeting Ankrd46, as implicated by a microarray-based analysis. Fertil. Steril. 103, 834–844.e4. doi:10.1016/j.fertnstert.2014.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, H.-J., Kim, H. J., Lee, J. H., Lim, D. H., Son, J. K., Kim, E.-T., et al. (2021a). Identification of plasma miRNA biomarkers for pregnancy detection in dairy cattle. J. Anim. Reprod. Biotechnol. 36, 35–44. doi:10.12750/JARB.36.1.35

CrossRef Full Text | Google Scholar

Lim, M., Thompson, J. G., and Dunning, K. R. (2021b). Hypoxia and reproductive health: Hypoxia and ovarian function: Follicle development, ovulation, oocyte maturation. Reproduction 161, F33–F40. doi:10.1530/REP-20-0509

PubMed Abstract | CrossRef Full Text | Google Scholar

Lip, S., Boekschoten, M., van Pampus, M., Scherjon, S., Plösch, T., and Faas, M. (2018). 103. Dysregulated circulating microRNAs in preeclampsia: The role of miR-574-5p and miR-1972 in endothelial dysfunction. Pregnancy Hypertens. 13, S22. doi:10.1016/j.preghy.2018.08.068

CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17, 10–12. doi:10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

Orang, V. A., Reza, S., and Kazemzadeh-Bavili, M. (2014). Mechanisms of miRNA-mediated gene regulation from common downregulation to mRNA-specific upregulation. Int. J. Genomics 2014, 1–15. doi:10.1155/2014/970607

CrossRef Full Text | Google Scholar

Pan, B., Toms, D., and Li, J. (2018). MicroRNA-574 suppresses oocyte maturation via targeting hyaluronan synthase 2 in porcine cumulus cells. Am. J. Physiol. Physiol. 314, C268–C277. doi:10.1152/ajpcell.00065.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, H., Suo, P., Liu, C., Wang, J., Zhou, S., Ma, X., et al. (2014). The ESR1 gene in unexplained recurrent spontaneous abortion. Syst. Biol. Reprod. Med. 60, 161–164. doi:10.3109/19396368.2013.877540

PubMed Abstract | CrossRef Full Text | Google Scholar

Pohler, K. G., Green, J. A., Moley, L. A., Gunewardena, S., Hung, W., Payton, R. R., et al. (2017). Circulating microRNA as candidates for early embryonic viability in cattle. Mol. Reprod. Dev. 84, 731–743. doi:10.1002/mrd.22856

PubMed Abstract | CrossRef Full Text | Google Scholar

Rabaglino, M. B., and Kadarmideen, H. N. (2020). Machine learning approach to integrated endometrial transcriptomic datasets reveals biomarkers predicting uterine receptivity in cattle at seven days after estrous. Sci. Rep. 10, 16981. doi:10.1038/s41598-020-72988-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rehmsmeier, M., Steffen, P., Hochsmann, M., and Giegerich, R. (2004). Fast and effective prediction of microRNA/target duplexes. RNA 10, 1507–1517. doi:10.1261/rna.5248604

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverter, A., and Chan, E. K. F. (2008). Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics 24, 2491–2497. doi:10.1093/bioinformatics/btn482

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandell, L. L., Guan, X.-J., Ingram, R., and Tilghman, S. M. (2003). Gatm, a creatine synthesis enzyme, is imprinted in mouse placenta. Proc. Natl. Acad. Sci. 100, 4622–4627. doi:10.1073/pnas.0230424100

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Owen, O., Baliga, N. S., Wang, J. T., Daniel Ramage, N. A., et al. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 6. doi:10.1101/gr.1239303.metabolite

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, B., Sun, H., Ahmad, M. J., Ghanem, N., Abdel-Shafy, H., Du, C., et al. (2021). Genetic features of reproductive traits in bovine and buffalo: Lessons from bovine to buffalo. Front. Genet. 12, 617128. doi:10.3389/fgene.2021.617128

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimizu, Y., Takeuchi, T., Mita, S., Notsu, T., Mizuguchi, K., and Kyo, S. (2010). Krüppel-like factor 4 mediates anti-proliferative effects of progesterone with G₀/G₁ arrest in human endometrial epithelial cells. J. Endocrinol. Invest. 33, 745–750. doi:10.1007/BF03346681

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M. D., Pillman, K., Jankovic-Karasoulos, T., McAninch, D., Wan, Q., Bogias, K. J., et al. (2021). Large-scale transcriptome-wide profiling of microRNAs in human placenta and maternal plasma at early to mid gestation. RNA Biol. 18, 507–520. doi:10.1080/15476286.2021.1963105

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, L., Zhao, S., Zhu, M., and Yu, M. (2010). Differential expression of microRNAs in porcine placentas on Days 30 and 90 of gestation. Reprod. Fertil. Dev. 22, 1175–1182. doi:10.1071/RD10046

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhde, K., van Tol, H., Stout, T., and Roelen, B. (2017). MicroRNA expression in bovine cumulus cells in relation to oocyte quality. Non-Coding RNA 3, 12. doi:10.3390/ncrna3010012

PubMed Abstract | CrossRef Full Text | Google Scholar

Vashukova, E. S., Kozyulina, P. Y., Illarionov, R. A., Yurkina, N. O., Pachuliia, O. V., Butenko, M. G., et al. (2021). High-Throughput sequencing of circulating MicroRNAs in plasma and serum during pregnancy progression. Life 11, 1055. doi:10.3390/life11101055

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Lin, J., Li, X., Han, B., Wang, L., Liu, M., et al. (2016). Transcriptome profile of one-month-old lambs’ granulosa cells after superstimulation. Asian-Australasian J. Anim. Sci. 30, 20–33. doi:10.5713/ajas.15.0999

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, P., Ma, Y., Wu, H., and Wang, Y.-L. (2021). Placenta-derived MicroRNAs in the pathophysiology of human pregnancy. Front. Cell Dev. Biol. 9, 646326. doi:10.3389/fcell.2021.646326

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Zhou, Y., Edelshain, B., Schatz, F., Lockwood, C. J., and Taylor, H. S. (2012). FKBP4 is regulated by HOXA10 during decidualization and in endometriosis. Reproduction 143, 531–538. doi:10.1530/REP-11-0438

PubMed Abstract | CrossRef Full Text | Google Scholar

Yen, P.-S., Chen, C.-H., Sreenu, V., Kohl, A., and Failloux, A.-B. (2019). Assessing the potential interactions between cellular miRNA and arboviral genomic RNA in the yellow fever mosquito, Aedes aegypti. Viruses 11, 540. doi:10.3390/v11060540

PubMed Abstract | CrossRef Full Text | Google Scholar

Zang, X., Zhou, C., Wang, W., Gan, J., Li, Y., Liu, D., et al. (2021). Differential MicroRNA expression involved in endometrial receptivity of goats. Biomolecules 11, 472. doi:10.3390/biom11030472

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Blattman, J. N., Kennedy, N. J., Duong, J., Nguyen, T., Wang, Y., et al. (2004). Regulation of innate and adaptive immune responses by MAP kinase phosphatase 5. Nature 430, 793–797. doi:10.1038/nature02764

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, G., Li, Y., Kang, X., Huang, L., Li, P., Zhou, J., et al. (2019). The study of blood transcriptome profiles in Holstein cows with miscarriage during peri-implantation. Asian-Australasian J. Anim. Sci. 32, 38–48. doi:10.5713/ajas.17.0793

PubMed Abstract | CrossRef Full Text | Google Scholar