Abstract Simple Summary The objective of this study was to identify differentially expressed (DE) sperm and seminal plasma microRNAs (miRNAs) in high- and low-fertile Holstein bulls (four bulls per group), integrate miRNAs to their target genes, and categorize target genes based on predicted biological processes. Out of 84 bovine-specific, prioritized miRNAs analyzed by RT-PCR, 30 were differentially expressed in high-fertile sperm and seminal plasma compared to low-fertile sperm and seminal plasma, respectively (p ≤ 0.05, fold regulation ≥5 magnitudes). Interestingly, expression levels of DE-miRNAs in sperm and seminal plasma followed a similar pattern. Highly scored integrated genes of DE-miRNAs predicted various biological and molecular functions, cellular process, and pathways. Further in silico analysis revealed categorized genes may have a plausible association with pathways regulating sperm structure and function, fertilization, and embryo and placental development. In conclusion, highly DE-miRNAs in bovine sperm and seminal plasma could be used as a tool for predicting reproductive functions. Since the identified miRNA-mRNA interactions were mostly based on predictions from public databases, the causal regulations of miRNA-mRNA and the underlying mechanisms require further functional characterization in future studies. Abstract Recent advances in high-throughput in silico techniques portray experimental data as exemplified biological networks and help us understand the role of individual proteins, interactions, and their biological functions. The objective of this study was to identify differentially expressed (DE) sperm and seminal plasma microRNAs (miRNAs) in high- and low-fertile Holstein bulls (four bulls per group), integrate miRNAs to their target genes, and categorize the target genes based on biological process predictions. Out of 84 bovine-specific, prioritized miRNAs analyzed by RT-PCR, 30 were differentially expressed in high-fertile sperm and seminal plasma compared to low-fertile sperm and seminal plasma, respectively (p ≤ 0.05, fold regulation ≥ 5 magnitudes). The expression levels of DE-miRNAs in sperm and seminal plasma followed a similar pattern. Highly scored integrated genes of DE-miRNAs predicted various biological and molecular functions, cellular process, and pathways. Further, analysis of the categorized genes showed association with pathways regulating sperm structure and function, fertilization, and embryo and placental development. In conclusion, highly DE-miRNAs in bovine sperm and seminal plasma could be used as a tool for predicting reproductive functions. Since the identified miRNA-mRNA interactions were mostly based on predictions from public databases, the causal regulations of miRNA-mRNA and the underlying mechanisms require further functional characterization in future studies. Keywords: dairy bulls, fertility, sperm, seminal plasma, micro RNA, mRNA, bioinformatics, reproductive function 1. Introduction MicroRNAs (miRNAs) are small non-coding RNAs, expressed in various tissues [[28]1,[29]2,[30]3] and in biofluids [[31]3,[32]4,[33]5,[34]6] that regulate the target genes’ expression post-transcriptionally and translationally. Extracellular miRNAs are detected in various body fluids, including blood serum and plasma, urine, saliva, semen, and milk [[35]7,[36]8]. Although miRNAs are potential biomarkers for developmental stages and diseases, neither the origin of extracellular miRNAs nor their communication between cells of origin and the target cells, are well known. Such studies are hindered by the challenges of isolating and characterizing RNA. Recently, we characterized mature miRNA expression patterns in boar sperm and seminal plasma [[37]9,[38]10]. Seminal plasma is a multifaceted biofluid that provides a nutritious and protein-rich environment that is necessary for the proper functioning of sperm such as sperm maturation, metabolism, motility, modification of sperm membranes, capacitation, acrosome reaction, interaction with oviductal and uterine epithelium, and fertilization [[39]11,[40]12]. Components of seminal plasma are energy substrates (fructose, sorbitol, pyruvate, and glyceryl phosphocholine), organic compounds (citric acid, peptides, proteins, lipid, hormones, and cytokines) including essential (threonine, tryptophan, and lysine) and non-essential amino acids (arginine, glycine, serine, tyrosine, and phenylalanine), nitrogenous compounds (ammonia, urea, uric acid, and creatinine), ions (Na^+, K^+, Zn^2+, Ca^2+, Mg^2+, Cl^−, and PO[4]^3−), reducing substances (ascorbic acids and hypo-taurine), and reactive oxygen species scavenging enzymes (glutathione peroxidase 5, thioredoxin, glutathione-s-transferase M1, superoxide dismutase 1, and peroxiredoxin) [[41]13,[42]14,[43]15,[44]16]. MicroRNAs are highly stable in biofluids. Expression of miRNAs is detectable in sperm and seminal plasma of all animal species [[45]17,[46]18,[47]19,[48]20,[49]21,[50]22]. The miRNA expression is greater in sperm compared to seminal plasma and follows a similar expression trend [[51]9,[52]10]. The miRNAs may transmigrate between sperm and seminal plasma [[53]23,[54]24]. Bovine sperm and seminal plasma consist of or are associated with several miRNAs; however, functions of most miRNAs remain unknown, and several miRNAs are involved in fertilization [[55]25,[56]26]. Absence, presence, under-, or over-expression of specific genes regulated by miRNAs could alter sperm functions, reduce fertilizing ability, and thus lower the fertility of an ejaculate. The objective of this study was to elucidate differentially expressed (DE) miRNAs in sperm and seminal plasma of high- and-low fertile bulls, integrate miRNAs to their target genes, and categorize target genes to predict biological processes. The hypothesis was that sperm and seminal plasma miRNAs expression differ between high- and low-fertile bulls. Differentially expressed miRNAs of high-fertile bulls and their top-ranked integrated target genes will be involved in the biological process critical for regulating sperm function, fertilization, and embryo development. 2. Materials and Methods 2.1. Ethics Statement This study was performed in strict accordance with the ethics, standard operating procedure, and use of the animal cells and biofluids for research. The protocol was approved by the institutional animal care and use committee of the Washington State University. 2.2. Semen Sample Processing Fresh bull semen was collected during the spring season using an artificial vagina. Holstein bulls (n = 8; age 3.4 ± 0.11 year) with high (n = 4) and low (n = 4) fertility based on sire conception rates ([57]https://www.aipl.arsusda.gov/reference/arr-scr1.htm and [58]https://uscdcb.com) (accessed on 10 August and 7 December 2016) were selected [[59]Table S1; high-fertility bulls with SCR ≥ 4 (mean reliability: 91%); low fertility bulls with SCR ≤ −2 (mean reliability: 82%), 6 percentage points mean difference in sire conception rates between high vs. low fertile bulls in this study]. Sire conception rate (SCR) evaluation is released for active artificial insemination Holstein bulls that have a minimum of 300 breeding in the last 48 months (at least 100 breeding in the last 12 months) in at least 10 herds. Immediately after the collection, % progressive motility and % abnormal spermatozoa were determined using Society for Theriogenology standards. None of the semen samples contained either immature germ or somatic cells. Undiluted semen was centrifuged at 1000× g at 4 °C for 20 min. The sperm pellet was washed in PBS and centrifuged three times to ensure the removal of seminal plasma. Thereafter, the seminal plasma supernatant was centrifuged at 16,000× g at 4 °C to ensure the removal of residual sperm. Aliquots of the sperm pellets and seminal plasma supernatant were flash-frozen to −196 °C and stored at −80 °C. Prior to use, it was thawed at room temperature for 15 min. All of the procedures were performed for sperm and seminal plasma samples from each bull separately. 2.2.1. Individual Progressive Motility (%) Individual progressive motility was assessed in an ejaculate diluted with warmed semen extender. A drop of diluted sperm was placed on a slide, covered with a coverslip, and examined at 400× using phase contrast microscopy. The proportion of sperm that were moving progressively across the field of view was estimated. 2.2.2. Abnormal Sperm (%) The morphology of individual sperm was determined by examining an eosin-nigrosin stained semen smear under oil immersion (1000×). A 4 or 5 mm drop of warm stain near the end of a warm microscopic slide and a drop of semen near the stain were placed, mixed using a Pasteur pipette, and a smear was prepared by pulling a drop of stained semen slowly across the slide. Percentage abnormal sperm was determined by estimating morphology of at least 100 sperm. 2.3. RNA Isolation 2.3.1. Sperm RNA was isolated from sperm [[60]9,[61]10] using RNeasy plus Universal Mini Kit (Qiagen Inc., Valencia, CA, USA), following the manufacturer’s instructions. In brief, 750 μL QIAzol Lysis Reagent (Qiagen Sciences Inc., Germantown, MD, USA) was added to the sperm pellet (approximately 100 × 10^6 sperm), completely homogenized, and held at room temperature for 5 min to help dissociation of nucleoprotein complexes. Then, 100 μL genomic DNA (gDNA) eliminator solution was added, the mixture was shaken vigorously to eliminate gDNA, and then 150 μL of chloroform was added. Following vigorous shaking and 2 to 3 min incubation at room temperature, the mixture was centrifuged (12,000× g, 15 min, 4 °C), the upper aqueous phase was harvested, 1.5 volume of 100% ethanol was added, and the mixture was mixed thoroughly by vigorous pipetting. The sample was then layered on a RNeasy mini spin column (included in the RNeasy Plus Universal Mini Kit), centrifuged (≥8000× g, 15 s, room temperature) and bound RNA was washed by centrifugation (≥8000× g, 15 s, room temperature) using buffers RWT and RPE consecutively. The RNA then was eluted using 60 × L RNase-free water. The concentration was measured in a Thermo Scientific NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). The purity of RNA was determined by determining the ratio of absorbance at 260 and 280 nm and samples were stored at −80 °C. 2.3.2. Seminal Plasma Small RNAs were purified from seminal plasma [[62]9,[63]10] using a miRNeasy serum/plasma kit (Qiagen Inc., Valencia, CA, USA). The kit included phenol/guanidine-based lysis of samples and silica-membrane column-based isolation of small RNAs. The kit was designed to isolate cell-free small RNAs. QIAzol reagent (750 μL) was added to 150 μL of thawed samples, mixed thoroughly, incubated for 5 min at room temperature, and then 3.5 μL miRNeasy serum/plasma spike-in-controls (lyophilized C. elegans miR-39 mimic; 1.6 × 10^8 copies/μL) and 150 μL chloroform were added. The mixture was shaken vigorously for 15 s, incubated for 2 to 3 min at room temperature, and then centrifuged (12,000× g for 15 min at 4 °C). The upper aqueous phase was harvested, avoiding the interphase, mixed with ~1.5 volumes of 100% ethanol and mixed by pipetting. Approximately half of the mix was transferred into a RNeasy MinElute spin column in a 2 mL collection tube and centrifuged (12,000× g, 15 s, room temperature). The flow-through was discarded and this step was repeated with the rest of the sample using the same column. The RNA was bound to the membrane; this membrane with bound RNAs was washed sequentially through 700 μL of buffer RWT and 500 μL of buffer RPE (≥8000× g, 15 s, room temperature). Then, the bound RNA was washed through 500 μL of 80% ethanol by centrifugation (≥8000× g, 2 min, room temperature). The RNeasy MinElute column was dried by centrifugation (12,000× g, 5 min), and the flow-through and collection tube were discarded. Small RNAs were then eluted by placing the membrane in 24 μL RNase-free water and centrifuging (12,000× g, 1 min). 2.4. Complementary DNA Synthesis of Sperm and Seminal Plasma miRNAs Total RNA containing miRNA was used as starting material. Mature miRNA was reverse transcribed into cDNA using miScript II RT kit (Qiagen Inc., Valencia, CA, USA). In brief, template RNA was thawed on ice, and 10× miScript Nucleics mix, 5× miScript HiSpec buffer, and RNase-free water were thawed at room temperature. Reaction components for a 20 μL reaction were 4 μL of HiSpec buffer, 2 μL Nucleics mix, 2 μL of reverse transcriptase enzyme mix, and 12 μL of RNA template containing 250 ng RNA. Reverse-transcription reaction components were gently mixed, briefly centrifuged (2000× g, 10 s), and kept on ice. The mixture was incubated at 37 °C for 60 min and then at 95 °C for 5 min in a Thermocycler (Thermo Fisher Scientific, San Francisco, CA, USA). After incubation, the product was placed on ice. The product containing cDNA equivalent to 250 ng RNA was diluted with 90 μL nuclease-free water and stored at −20 °C prior to real-time PCR. 2.5. Sperm and Seminal Plasma Mature Mirna Profiling Using Real-Time PCR Real-time PCR profiling of sperm and seminal plasma mature miRNAs using miScript miRNA PCR arrays in combination with the miScript SYBR Green PCR Kit (Qiagen Sciences Inc., Germantown, MD, USA), which contained the miScript Universal reverse primer and QuantiTect SYBR Green PCR Master Mix, was performed. A bovine miRNome miScript miRNA PCR array 96-well plate 1 ([64]Table 1) was used in this study. The array plate consisted of specific primers to identify 84 highly prioritized bovine mature miRNAs from the most current miRNA genome database, as annotated in miRBase version 20 ([65]www.miRBase.org) (accessed on 19 August 2015) [[66]5]. A set of controls facilitated threshold cycle normalization, assessing reverse transcription performance, and assessing PCR performance. It should be noted that the 84 miRNAs specific to bovine species were selected based on the information available via experiments and computations in the literature and miRbase. In addition, genes involved in male reproductive functions and associated miRNAs were also taken into consideration. Table 1. Bovine miRBase profiler plate, consisting of primers for 84 target miRNAs and control genes. Layout 1 2 3 4 5 6 7 8 9 10 11 12 A bta-let-7f bta-miR-101 bta-miR-103 bta-miR-125a bta-miR-125b bta-miR-126-3p bta-miR-128 bta-miR-145 bta-miR-148a bta-miR-151-3p bta-miR-151-5p bta-miR-16b B bta-miR-181a bta-miR-18a bta-miR-18b bta-miR-199a-5p bta-miR-205 bta-miR-20a bta-miR-21-5p bta-miR-221 bta-miR-222 bta-miR-26a bta-miR-26b bta-miR-27a-3p C bta-miR-27b bta-miR-29a bta-miR-300-5p bta-miR-30d bta-miR-31 bta-miR-320a bta-miR-34b bta-miR-484 bta-miR-499 bta-miR-99a-5p bta-miR-7a-5p bta-let-7d D bta-let-7g bta-let-7i bta-miR-17-5p bta-miR-107 bta-miR-10a bta-miR-10b bta-miR-122 bta-miR-124b bta-miR-127 bta-miR-132 bta-miR-138 bta-miR-139 E bta-miR-140 bta-miR-142-3p bta-miR-142-5p bta-miR-148b bta-miR-150 bta-miR-15b bta-miR-17-3p bta-miR-17-5p bta-miR-181b bta-miR-181c bta-miR-186 bta-miR-191 F bta-miR-192 bta-miR-193a-3p bta-miR-193a-5p bta-miR-199a-3p bta-miR-199b bta-miR-200a bta-miR-200b bta-miR-200c bta-miR-20b bta-miR-210 bta-miR-21-3p bta-miR-214 G bta-miR-215 bta-miR-218 bta-miR-22-5p bta-miR-23a bta-miR-23b-3p bta-miR-24-3p bta-miR-25 bta-miR-29b bta-miR-29c bta-miR-30a-5p bta-miR-30c bta-miR-30e-5p H cel-miR39-3p cel-miR39-3p SNORD42B SNORD69 SNORD61 SNORD68 SNORD96A RNU6-6P miRTC miRTC PPC PPC [67]Open in a new tab Characterized 84 target miRNAs (plate well positions A1–G12) and controls (plate well positions H1–H12). The reaction mix was prepared with 1375 μL of 2× QuantiTect SYBR Green PCR master mix, 275 μL of 10× miScript universal primer, 1000 μL of RNase free water, and 100 μL of template cDNA for each 96-well plate. A 25 μL volume of the reaction mixture was added to each well and the template was amplified in StepOnePlus cycler (Applied Biosystems, Foster City, CA, USA). Cycling conditions consisted of an initial heating step at 95 °C for 15 min. Forty cycles included 15 s denaturation step at 94 °C, 30 s annealing step at 55 °C, and 30 s extension step at 70 °C. Dissociation curve analysis was performed to verify specificity and identity. 2.6. Bovine Mature miRNA PCR Array Analysis Eighty-four high-priority bovine mature miRNAs were selected from the miRNA genome database to elucidate DE-miRNAs in sperm and seminal plasma samples. Reverse transcription and positive controls were chosen to ensure the efficiency of the array, reagents, and instrument. Raw CT data (.xlsx file format) were uploaded to the data analysis center ([68]http://pcrdataanalysis.sabiosciences.com/pcr/arrayanalysis.php) (accessed on 10 May 2017 and 16 July 2018). Data quality control was examined to assess amplification reproducibility and reverse transcription efficiency, and to detect any other contamination in amplified samples [[69]5,[70]10]. The controls were cel-miR-39-3p, SNORD61, SNORD68, SNORD72, SNORD95, SNORD96A, RNU6-2, miRTC, and PPC. The CT values of samples were calibrated to the CT values of cel-miR-39-3p. Global CT mean of expressed miRNAs was chosen to normalize the target sperm and seminal plasma miRNAs. SNORDs and RNU6-2 served as internal normalizers. Two reverse transcription controls and two positive controls ensured the efficiency of the array, the reagents, and the instrument. The distribution of CT values and raw data average in both groups were reviewed. Average ΔCT, 2^−ΔCT, fold change, P-value, and fold regulation were calculated in the web-based program, and p-values were included in subsequent graphical analyses. Average CT values were converted into linear 2^−ΔCT values and p-values were calculated with a student’s t-test. 2.7. Bioinformatics Analysis 2.7.1. Conserved Nucleotide Sequences Nucleotide sequences of DE-miRNAs of bovine, human, and mice were retrieved from miRBase ([71]www.mirbase.org) (accessed on 10 December 2021). Sequences of miRNAs from humans and mice were compared with the sequences of bovine species to ensure the similarity [[72]27,[73]28]. 2.7.2. Identification of Target and Predicted Genes of Differentially Expressed miRNAs The target genes of DE-miRNAs were predicted using miRNet ([74]http://www.mirnet.ca/) (accessed on 16 December 2021) [[75]29]. This tool integrates data from different miRNA databases (TarBase, miRTarBase, and miRecords). The prediction analysis was performed for upregulated and downregulated DE-miRNAs separately. 2.7.3. Construction of Protein-Protein Interaction Network and Screening of Hub Genes The protein-protein interaction (PPI) network of DE-miRNAs’ predicted target genes was performed by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database ([76]http://stringdb.org/) (accessed on 16 December 2021) [[77]30]. Gene Ontology (GO) functional annotation for biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for predicted target genes of DE-miRNAs were performed. A p-value of <0.05 was regarded as statistically significant. The PPI interaction was exported to the Cytoscape software (version 3.9) and visualized [[78]31]. The hub genes were selected out as the top 20 nodes of the PPI network using the Maximal Clique Centrality (MCC) method [[79]32], which had a better performance on the precision of predicting top essential proteins. Further analysis was performed using ClueGO (Robertson and Sharkey, 2016) to integrate GO terms as well as KEGG pathways and create a functionally nested or organized GO/pathway term (k score = 3). This task analyzes one set of genes or compares two lists of genes and comprehensively visualizes functionally grouped genes [[80]33]. 2.7.4. Gene Ontology and Functional Annotation Analysis To understand the biological meaning behind DE-miRNAs and integrated genes, biological processes were investigated in the PANTHER (Protein ANalysis THrough Evolutionary Relationships) Classification System. ([81]www.geneontology.org/, [82]www.pantherdb.org/) (accessed on 22 December 2021). Further, the important GO terms and their pertinent cooccurring terms were analyzed using QuickGO ([83]https://www.ebi.ac.uk/QuickGO) (accessed on 22 December 2021). 2.7.5. Real-Time Polymerase Chain Reaction for Determining mRNA Expression of Hub Genes Genes such as dna methyltransferase 1 (DNMT1), forkhead box P3 (scurfin) (FOXP3), phosphatase and tensin homolog (PTEN) and Zinc Finger E-Box Binding Homeobox 1 (ZEB1) were selected from the group of hub genes to substantiate the mRNA expressions in sperm and seminal plasma samples. Total RNA extraction and complementary DNA synthesis were performed as previously described [[84]34,[85]35]. In brief, sperm and seminal plasma samples were used to extract RNA by TRizol (Invitrogen, Carlsbad, CA, USA). The RNA concentration and quality were determined using NanoDrop 1000 spectrophotometer and all RNA samples were treated with DNAse I (Invitrogen) to remove the DNA contaminant. Complementary DNA was synthesized using the iScript cDNA synthesis kit (Bio-Rad Laboratories Inc., Hercules, CA, USA) and stored at −20 °C. Specific primer pairs ([86]Table 2) for the hub genes were designed using primer-BLAST ([87]www.ncbi.nlm.nih.gov/tools/primer-blast/, accessed on 1 July 2022). Prior to real-time PCR, ethidium bromide-stained electrophoresis gel for the amplicon of the expected size was performed ([88]Figure S1). Real-time PCR was conducted using Fast SYBR Green Master Mix (Applied Biosystems, Foster City, CA, USA) as previously described [[89]25] following the manufacturer’s instruction. Endogenous control glyceraldehyde-3-phosphate dehydrogenase (GADPH) was used to normalize the threshold cycle (CT) values. Fold comparisons were made between the high- and low-fertile groups. Table 2. Forward and reverse primer sequence for quantitative real-time polymerase chain reaction amplification of mRNA for canine testis samples. Gene Forward Primer Reverse Primer Product Length Accession Number ZEB1 AAAGCAGCAGGGCGAGTTAT TATGGGGTTGGCACTTGGTG 181 [90]NM_001206590.1 DNMT1 TATCGGCTGTTCGGCAACAT GGCAGCCTCCTCCTTGATTT 153 [91]NM_182651.2 FOXP3 CAGCGGACACTCAACGAGAT AACTCATCCACGGTCCACAC 164 [92]XM_024987818.1 PTEN GCAGCTTCTGCCATCTCTCT ATGCTTTGAATCCAAAAACCTTACT 235 [93]NM_001319898.1 GADPH GTGAAGGTCGGAGTGAACGG ATTGATGGCGACGATGTCCA 93 [94]NM_001034034.2 [95]Open in a new tab DNMT1, dna methyltransferase 1; FOXP3, forkhead box P3 (scurfin); PTEN, phosphatase and tensin homolog; ZEB1, Zinc Finger E-Box Binding Homeobox 1; GADPH, glyceraldehyde-3-phosphate dehydrogenase. 2.7.6. Protein Immunoblots Proteins such as DNMT1, FOXP3, PTEN, and ZEB1, were selected to substantiate the protein expressions in sperm and seminal plasma samples. Western blots for sperm and seminal plasma samples were performed by methods described previously [[96]36]. In brief, protein extraction methods included: addition of protease and phosphatase inhibitor to sperm and seminal plasma samples accordingly; homogenization; lysate incubation at 4 °C for 45 min; centrifugation (at 12,000× g for 20 min); and determination of protein concentrations. Protein samples were then electrophoresed through 12% SDS-PAGE gel (Bio-Rad Laboratories, Philadelphia, PA, USA) and then transferred onto PVDF membrane (Bio-Rad Laboratories, Hercules, CA, USA). Sixty micrograms of protein lysate were used per lane. The membrane blots were incubated in 10% goat serum in PBS to block non-specific binding. After overnight incubation at 4 °C with primary antibodies (mouse monoclonal to DNMT1 (Catalog # MA5-16169), rabbit polyclonal to PTEN (Catalog # 600-401-859), and mouse polyclonal to ACBT (sc-47778) from Santa Cruz Biotechnology, Santa Cruz, CA, USA), membranes were washed in buffer containing 2% animal serum and 0.1% detergent. The membranes were then incubated in secondary antibodies (goat anti-mouse IgG-FITC for DNMTA1 and ACBT (sc-2010; Santa Cruz Biotechnology, Dallas, TX, USA) and goat anti-rabbit IgG-FITC for PTEN (sc-2012; Santa Cruz Biotechnology, Dallas, TX, USA) for 1 h at room temperature. The blots were then washed and scanned using the Pharos FX Plus system (Bio-Rad Laboratories, Hercules, CA, USA). FITC fluorophore was excited at 488 nm and read at the emission wavelength of 530 nm. All possible negative controls were included. 2.7.7. Statistical Analyses to Determine Differences in mRNA Expression The average mRNA CT values were converted into linear 2^−ΔCT values and differences in relative expressions between high- and low-fertile bulls were calculated with a student’s t-test. Data were analyzed with a statistical software program (SAS version 9.4 for Windows, SAS Institute, Cary, NC, USA). Correlation coefficients were estimated using PROC CORR to determine the association between miRNA and mRNA fold changes. The RT-PCR data were analyzed by t-test, using 2-DDCt values to ascertain statistical significance of any differences in mRNA expressions between high- and low-fertile bulls. For all statistical analyses, p ≤ 0.05 was considered significant. 3. Results The sire conception rate (SCR), progressive motility (%), and abnormal sperm (%) for high- and low-fertile Holstein bulls used in the study are given in [97]Supplementary File S1. On miRNA semiquantitative profiling, 32 miRNAs were differentially expressed (p ≤ 0.05; fold change magnitude cut-off at 5, high) in high compared with low-fertile sperm and seminal plasma ([98]Figure 1A). Of these, 20 miRNAs were upregulated (≥5), and 12 miRNAs were down regulated (≤−5) in the high-fertile sperm and seminal plasma relative to low-fertile sperm and seminal plasma (p < 0.05). When fold change cut-off at 2 was considered, there were 56 miRNAs (33 upregulated (≥2) and 23 down regulated (≤−2)) were differentially expressed in high compared with low-fertile sperm and seminal plasma (p ≤ 0.05). Further, when fold change cut-off at 10 (very high) was considered, there were 11 miRNAs (9 upregulated (≥10) and 2 down regulated (≤−10)) were differentially expressed in high compared with low-fertile sperm (p < 0.001); whereas 10 miRNAs (6 upregulated (≥10) and 4 down regulated (≤−10)) were differentially expressed in high compared with low-fertile seminal plasma ([99]Figure 1B; p < 0.001). Figure 1. [100]Figure 1 [101]Open in a new tab (A) Fold regulation of differentially expressed sperm and seminal plasma miRNAs in high-fertile compared to low-fertile Holstein bulls. Of 84 bovine-specific well-characterized miRNAs investigated, 20 were greater (p ≤ 0.05; fold ≥ 5) and 12 were lower (p ≤ 0.05; fold ≤ −5) in sperm and seminal plasma in high-fertile Holstein bulls; (B) Fold regulation of highly differentially expressed sperm and seminal plasma miRNAs in high-fertile compared to low-fertile Holstein bulls. Of 84 bovine-specific well-characterized miRNAs investigated, 11 miRNA [9 upregulated (≥10) and 2 down regulated (≤−10)) were differentially expressed in high compared with low-fertile sperm (p < 0.001); whereas 10 miRNA (6 upregulated (≥10) and 4 down regulated (≤−10)) were differentially expressed in high compared with low-fertile seminal plasma. Interestingly, the differential expression of miRNAs in high- vs. low-fertile sperm and seminal plasma followed a similar trend. Fold changes were slightly higher for all upregulated and slightly lower for all downregulated miRNAs in sperm compared with seminal plasma. Fold changes for all 84 miRNAs in sperm and seminal plasma for high-fertile bulls (relative to low-fertile bulls) are given in [102]Supplementary Files S1 and S2, respectively. The ratio of sperm: seminal plasma miRNA fold expression differences varied from 0.75 to 1.40 in high-fertile bulls. Since DE-miRNAs (upregulated and downregulated) had a similar pattern for sperm and seminal plasma, the miRNA-mRNA interaction analysis, construction of PPI network, and Gene Ontology and functional annotation analysis were performed once using differentially expressed miRNAs with a fold change cut-off at 5. Nucleotide sequence similarities for the DE-miRNAs for humans, mice, and cattle are presented in [103]Table S2. Bovine sequences were very similar to human and mouse nucleotide sequences. Therefore, human miRNA IDs were used to construct miRNA-mRNA interaction network and functional enrichment analysis. The miRNA-mRNA interaction analysis for the upregulated miRNAs resulted in 27881 target genes and 75 predicted genes ([104]Supplementary File S3). The PPI for the 75 predicted genes (73 nodes and 188 edges, PPI enrichment p < 1.0 × 10^−16 revealed 299 significantly enriched GO terms biological processes (False Recovery Rate, p < 0.05) and 35 significant (False Recovery Rate, p < 0.05) KEGG enrichment pathways ([105]Supplementary File S4). The miRNA-mRNA interaction analysis for the downregulated miRNAs resulted in 24380 target genes and 57 predicted genes ([106]Supplementary File S5). The PPI for the 57 predicted genes (56 nodes and 186 edges, PPI enrichment p < 1.0 × 10^−16) revealed 391 significantly enriched GO terms biological processes (False Recovery Rate, p < 0.05) and 50 significant (False Recovery Rate p < 0.05) KEGG enrichment pathways ([107]Supplementary File S6). The PPI networks for predicted genes for upregulated miRNAs and downregulated miRNAs were separately constructed ([108]Figure 2A,B) using the STRING database and the Cytoscape software. According to the degree (MCC method), the top 20 hub genes in the networks for upregulated miRNAs and downregulated miRNAs were screened and presented in [109]Figure 3A,B. To decipher functionally nested gene ontology and pathway annotation networks for the predicted genes of upregulated and downregulated miRNAs in higher fertile sperm and seminal plasma, ClueGo nested network analysis was performed, and the results are presented in [110]Figure 4 and [111]Figure 5. Figure 2. [112]Figure 2 [113]Open in a new tab STRING protein-protein interaction (PPI) network. (A) PPI network for the upregulated DE-miRNAs predicted 75 genes (73 nodes and 188 edges, PPI enrichment p < 1.0 × 10^−16; (B) PPI for the upregulated DE-miRNAs predicted 57 genes (56 nodes and 186 edges, PPI enrichment p < 1.0 × 10^−16); the color nodes represent proteins. The edges represent interactions. Figure 3. [114]Figure 3 [115]Open in a new tab Interaction of hub genes of DE-miRNAs in the PPI network. (A) PPI network of top genes for highly upregulated DE-miRNAs; (B) PPI network of the top genes for highly downregulated DE-miRNAs. DE-miRNAs, differentially expressed microRNAs; PPI, protein-protein interaction; the PPI among hub genes for upregulated DE-miRNAs were greater compared with hub genes for down-regulated DE-miRNAs. Color red to yellow denotes high to low degree of expression. Black lines indicate interactions between genes. Figure 4. [116]Figure 4 [117]Open in a new tab ClueGO analysis of up-regulated genes in sperm and seminal plasma from high-fertile bulls. (A) GO/pathway terms specific for upregulated genes. The bars represent the number of genes associated with the terms. The percentage of genes per term is shown as bar label; (B) overview chart with functional groups including specific terms for upregulated genes; (C) functionally grouped network with terms as nodes linked based on their kappa score level (≥0.4), where only the label of the most significant term per group is shown. The node size represents the term enrichment significance. Functionally related groups partially overlap. The color gradient shows the gene proportion of each cluster associated with the term. Single (*) or double (**) asterisk indicate significant enriched GO terms at the p < 0.05 and p < 0.01 statistical levels. Figure 5. [118]Figure 5 [119]Open in a new tab ClueGO analysis of down regulated genes in sperm and seminal plasma from high-fertile bulls. (A) GO/pathway terms specific for down regulated genes. The bars represent the number of genes associated with the terms. The percentage of genes per term is shown as bar label; (B) overview chart with functional groups including specific terms for down regulated genes; (C) functionally grouped network with terms as nodes linked based on their kappa score level (≥0.4), where only the label of the most significant term per group is shown. The node size represents the term enrichment significance. Functionally related groups partially overlap. The color gradient shows the gene proportion of each cluster associated with the term. Single (*) or double (**) asterisk indicate significant enriched GO terms at the p < 0.05 and p < 0.01 statistical levels. Differentially expressed miRNAs, associated hub genes, and their linked reproductive functions are given in [120]Table 3. Finally, GO terms such as fertilization, implantation, trophoblast formation, placental development, and in utero embryonic development (from the 20 hub gene predicted 640 terms) were used to retrieve (QuickGO) relevant co-occurring GO terms ([121]Supplementary File S7) and from the GO functional annotation terms, fertilization ([122]Figure 6A) and progeny development ([123]Figure 6B) were constructed. Table 3. Differentially expressed miRNAs, associated predicted hub genes and their linked reproductive functions. Upregulated miRNAs Hub Genes Reproductive Functions Species References