The Republic of Turkey's Ministry of Agriculture and Forestry's Institute of Bahri-Dagdaş International Agricultural Research supplied information on the fertility phenotypes of adult rams. Pregnancy rates from natural mating were used to determine the fertility of mature rams (n = 66) at least 3–4 years old during the breeding seasons of 2017–2018–2019. The ewes' estrus was detected using teaser rams, who were not permitted to mate, by covering the prepuce area. The teaser rams were presented to the ewes early in the morning for about 30 min. The estrus was considered to be ewes seeking, standing for teasing, and allowing mount attempts by teaser rams. A handler selected estrous ewes and brought them into an enclosure along with a randomly selected single ram for natural mating. Throughout the breeding season, estrus detection was maintained, and ewes were accepted to mate with the chosen ram at random. Ewes were regarded as pregnant if they did not return to estrus within 35 days of mating. In addition, the number of pregnant and non-pregnant ewes for each ram was confirmed by matching the mating and lambing dates according to the duration of pregnancies. The rams' fertility scores were rated based on their conception rates. The rams were split into two groups based on their fertility levels: the high-fertility (HF) group (n = 31; 94.5 ± 2.8%) and the low-fertility (LF) group (n = 25; 83.1 ± 5.73%). The average pregnancy rate was 89.4 ± 7.2% (n = 66), and the rams were divided into these groups based on their fertility levels. We conducted an a priori power analysis using G^*Power3 (V3.1.9.4) to test the differences between two independent group means using a two-tailed test, an effect size (d = 2.52), and an alpha level of .05. The result showed that a total sample of eight animals with two equal-sized groups of n = 4 was required to achieve a power of .80 for RNA profiling. However, we excluded the rams that did not have 50 mating and were not used, so we profiled four animals for each group out of a total of 56. During the breeding season, each ram served at least 50 ewes in both groups.

The Bahri-Dagdaş Research Center Ethical Committee, Turkey (Number: 22.12.2016/58), approved the animal procedures. We trained the rams to obtain sperm using an artificial vagina (AV) that enabled them to mount on teaser ewes during estrus. Rams were permitted to ejaculate into the AV upon mounting. The first three collections were discarded prior to the collection of research samples, which was followed by semen collection and processing for use in research. Fertility rates that were 1 standard deviation above or below the mean were termed outliers. Four rams with the highest fertility (pregnancy rate; % 99.2 ± 1.6) and four rams with the lowest fertility (pregnancy rate; % 73.6 ± 4.4) were selected for lncRNA sequencing with high confidence. Subsequently, we collected about 2 × 10⁹/ml spermatozoa per ejaculate, and then the aliquots from each sample were adjusted to a final concentration of 10⁷/ml in straws and frozen at −80°C until lncRNA analysis.

Prior to RNA isolation, we purified sperm by filtering the semen samples with a 500-mesh sieve to eliminate cell debris. Then, samples were treated with a somatic cell lysis solution (0.3% Triton X-100 and 0.1% SDS in DEPC-treated H₂O) for 30 min on ice to eradicate somatic cells, followed by microscopic analysis of non-sperm cell contamination. We isolated total RNA from the purified ram sperm (n = 4, for each group) samples using the SanPrep column microRNA miniprep kit (Bio Basic Inc, Canada) with slight modification using the manufacturer's protocols. We added 800 μl of a guanidine–thiocyanate lysis buffer enriched in 20 mM DL-dithiothreitol onto the pellet, and then sperm cells were homogenized by passing the samples through a 26-G needle syringe 20–25 times. After other contaminants were thoroughly removed and total RNA was attached to the membrane, an on-column DNase digestion was carried out to remove any traces of DNA contamination. We evaluated the concentration and integrity of the total RNA samples using a NanoDrop (Colibri Microvolume Spectrometer, Titertek-Berthold, Germany) and a 2100-Bioanalyzer with the RNA 2100 Nano Chip (Applied Biosystems, Carlsbad, CA, USA), respectively.

Each RNA sample was utilized to prepare 2 g of RNA for the RNA library, and ribosomal RNA was first eliminated using the Ribo-ZeroTM rRNA Removal Kit from Epicenter. Then, using the NEBNext^® UltraTM-Directional RNA Library Prep Kit from Illumina^® (NEB, USA) in compliance with the commercial kit protocol, sequencing libraries were developed utilizing the rRNA-depleted RNA. Reverse transcription was used to create the first strand of the cDNA following fragmentation with an average length of 200 bp. The Agilent Bioanalyzer 2100 equipment was used to evaluate library quality following product purification using the AMPure XP system. As a result, libraries underwent sequencing using the Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA, United States), and Novogene Corporation (Beijing, China) generated 150-bp paired-end reads.

Initially, rRNA, adapter sequences, empty reads, and low-quality reads were eliminated from the raw data. All trimmed reads were confirmed to satisfy the quality threshold (Q-score; Q20 and Q30) to ensure that there was no bias in the evaluation step toward approaches that favor maximum read. The Phred scale indicating the reliability of base-calling, with Q20 representing a base call accuracy of 99% (or a 1% chance of error) and Q30 representing a base call accuracy of 99.9% (or a 0.1% chance of error) was used as the quality score. The Ovis aries (v4.0) reference genome was indexed with Bowtie v2.0.6, and the processed paired-end reads were mapped to that genome with HISAT2 2.1.0 (37). Each sample's mapped reads were constructed using StringTie (v1.3.1) (38). Finally, Cuffcompare, a program included in Cufflinks (v2.1.1), was used to annotate the assembled transcripts.

To classify newly screened lncRNAs with respect to their positional relationship with known mRNAs, putative lncRNAs were identified. To minimize the false-positive rate (FDR), assembled transcripts were classified to retrieve putative lncRNAs, such as lincRNA, antisense lncRNA, intronic lncRNA, and sense-overlapping. (A) Transcripts with a single exon were eliminated. (B) Transcripts with fewer than 200 nucleotides were eliminated. (C) Using Cufflinks v2.1.1, the annotated lncRNAs in the database were used to exclude the transcripts that overlapped with the exon region of the database annotation. (D) All transcripts with modest levels of expression [FPKM thresholds were set for the categorization of transcript expression levels; genes with very low or no expression (FPKM < 0.5), and FPKM = 0.5 was chosen as the cutoff to filter out the average read coverage per transcript which was much higher than the other transcript] were omitted (FPKM of a single exon transcript). Using three methods, namely, Coding-Non-Coding-Index (CNCI) (39), Pfam-scan (40), and coded potential calculator (CPC) (41), all estimated transcripts with coding potential were filtered out, and a set of putative lncRNAs was compiled from those with noncoding potential. Using Cuffcompare, the distinct categories of lncRNAs were obtained.

The FPKM value was used to evaluate levels of mRNA and lncRNA expression. Cuffdiff (v2.1.1) was used to determine lncRNA FPKM values. Later, the statistically significant DE genes were quantified by a log₂-fold change higher or equal to 2 (P-value < 0.05) or P-adjust < 0.05 (applied correction for multiple testing to the P-values to FDR), using Ballgown (42).

DE lncRNAs were selected to determine target genes. Pearson's correlation was used to assess potential coexpression between lncRNAs and mRNAs. A Pearson correlation >0.7 and a P-value of 0.05 were used to determine a positive association between an lncRNA and an mRNA. LncRNAs can act as cis regulators by remodeling factors onto local chromatin. We described cis-modulated genes as protein-coding genes that were coexpressed with one dysregulated lncRNA and were within 30 kb upstream or downstream in genomic distance in the same allele. To participate in certain biological processes, key transcription factors (TFs) are regulated in a trans manner by unique lncRNAs. As a result, we matched these lncRNAs' coexpressed mRNAs to mRNAs that were regulatory targets of specific TFs to anticipate that these lncRNAs might be involved in pathways controlled by these TFs.

We assessed DE mRNAs for gene ontology (GO) enrichment analysis using g:Profiler (43). GO terms possessing a corrected P-value < 0.05 were accepted as significantly enriched by DE genes. The statistical enrichment of lncRNA target genes was examined for KEGG pathway functional analysis based on the reactome pathway database by WebGestalt (WEB-based Gene SeT AnaLysis Toolkit) with P < 0.05 and FDR < 5.0% (44). The following parameters were specified for the enrichment analysis. A particular organism was labeled as Ovis aries (sheep). Sequential GO analyses (biological process; BP, cellular component; CC, and molecular function; MF) were conducted. The g:SCS method is used to compute multiple testing corrections for P-values based on GO and pathway enrichment analysis, with the user threshold set at 0.05.

Data were analyzed through SPSS software (version 22.0). Statistical plots were generated using GraphPad Prism 9 (GraphPad Software, USA). During the experimental process, four biological replicates were included with the measurement repeated twice. LF vs. HF groups were analyzed using an independent t-test. Significance was accepted at P-value = 0.05.

In the current study, we pooled two cDNA libraries out of the eight total RNAs (LF and HF, each = 4) obtained from low- and high-fertility ram sperm. The total number of raw readings obtained from all cDNA libraries was 91,108,174. Upon filtering the reads, a total of 89,252,823 clean reads were obtained. The Q20 (%) percentages were 96.06 and 95.39 for LF and HF, respectively. The Q30 (%) percentages were shown to be 90.96 and 89.50 for LF and HF, respectively. The percentages of GC content (%) were shown to be LF, 60.89, and HF, 62.37 (Table 1). Furthermore, after aligning the clean reads with the ovine reference genome through the TopHat2 algorithm, we discovered that the total mapped reads or fragments referring to all samples exceeded 75% and were mapped in the reference genome (Table 2).

The distribution intensity of the total mapped reads was split up and computed for both the plus and minus strands within each chromosome (Supplementary Table 1). The distributions of reads in the known RNA types for LF and HF are depicted in Table 3. The corresponding data related to these distributions can be found in Supplementary Table 2.

We applied the lncRNA filtering method to determine the novel lncRNA candidates in ram sperm. The specific filtering technique is outlined in Supplementary Table 3. Using this lncRNA filtering process, we identified a total of 14,352 lncRNAs, of which 11,209 were novel candidates, the data for which are supplied in Supplementary Table 4. We utilized three distinct software packages (CNCI, CPC, and Pfam) to estimate the potential protein-coding ability of the transcripts (Figure 1).

The primary characteristics of the lncRNAs' biotype distribution and length, exon intensity, and ORF length were examined as a consequence of the sequence analysis. Using Cufflinks, we divided lncRNAs into three groups: long intervening noncoding RNA (lincRNA), antisense lncRNA, and intronic lncRNA. Of these, 8.3%, 8.5%, 9.8%, 72.6%, and 0.8% of the attained lncRNAs were sense intronic, antisense, sense overlapping, lincRNA, and others, respectively (Figure 2A). Using this classification scheme, we revealed that the majority of the sperm-specific lncRNAs (72.6%) were lincRNAs. The exon lengths of the obtained lncRNAs were between 126 and 20,040 bp, and approximately 37% of the lncRNAs were demonstrated to be intense in 946–1,356 bp length (Figure 2B), with 1,542 bp being the median value. The exon length of the obtained mRNAs ranged between 75 and 26,726 bp, with around 10% of the mRNAs being intense at 835–1,025 bp (Figure 2B). Furthermore, lncRNAs have been shown to be more abundant in the first and fourth exons, whereas mRNAs were found to be more abundant in the first and seventh exons (Figure 2C). The length of the ORF of the lncRNAs was between 54 and 20,037 bp, with a median value of 1,461 bp (Figure 2D). For mRNAs, the ORF length was reported as 3–26,460, and the median value was 1,332.

The expressions associated with transcripts of lncRNA, mRNA, novel lncRNAs, and TUCP were evaluated using Cuffdiff. The most common approach for predicting gene expression levels, FPKM, is based on the effects of gene length on sequence depth and fragments described in the RNA sequence. The FPKM values associated with samples were calculated (Supplementary Table 5). Upon comparing FPKM values between LF and HF, no significant change was observed (Figures 3A, B). Moreover, the expression levels related to novel lncRNAs were the highest among the transcripts. TUCP and novel mRNA transcripts were greater than those of lncRNA and mRNAs, while lncRNA transcripts were similar to those of mRNA (Figures 3C, D).

Cuffdiff was employed to identify differentially expressed lncRNAs, mRNAs, and TUCPs. As a result, in HF and LF ram sperm, 93 lncRNAs and 275 mRNAs were reported to be differentially expressed (DE). We discovered 49 lncRNAs that are significantly upregulated and 44 that are downregulated between LF and HF groups. In addition, we discovered that 141 mRNAs were upregulated, whereas 134 were downregulated. Volcano plots show the up and downregulation highlights (Figures 4A, B; Supplementary Table 6). We also examined the expression patterns of DE lncRNAs and mRNAs using hierarchical clustering analysis, which serves as an additional approach for shedding light on differentially expressed genes by grouping genes with comparable expression patterns. The transcript FPKMs were utilized for hierarchical clustering, with discrete colors indicating the direction of the expression level. The clustering of genes on the left was caused by similar expressions (fold change >2, P < 0.05) and LF and HF ram sperm, whereas the expression change from blue to red indicated them as increasingly increased (Figures 4C, D).

Following a study at the transcriptome level, 275 mRNAs were selected for enrichment analysis due to their differential expression levels (141 upregulated and 134 downregulated). The top eight GO annotations terms in biological process (BP) were significantly enriched in the DE mRNAs, namely, the biological process (GO:0008150), cellular process (GO:0009987), regulation of cellular process (GO:0050794), organelle organization (GO:0006996), flagellated sperm motility (GO:0030317), sperm motility (GO:0097722), and cilium movement involved in cell motility (GO:0060294). In addition, in cellular component, intracellular anatomical structure (GO:0005622), sperm flagellum (GO:0036126), and motile cilium (GO:0031514) were among the top GO terms and binding (GO:0005488), protein binding (GO:0005515), enzyme binding (GO:0019899), and protein C-terminus binding (GO:0008022) in molecular function (MF; Figure 5A). The detailed definition of the GO terms is presented in Supplementary Figure 1. According to the statistics of the pathway enrichment, the top 16 KEGG pathways, including complex I biogenesis, the citric acid (TCA) cycle, respiratory electron transport, respiratory electron transport, mitochondrial translation, and gene silencing by RNA, had higher concentrations of DE mRNAs. ABHD2, AK1, CABS1, ROPN1, SEPTIN2, SLIRP, and TEKT3 genes were commonly enriched for 9 + 2 motile cilium, sperm flagellum, and motile cilium. Moreover, CABS1, CCDC39, CFAP97D1, ROPN1, SLIRP, TEKT3, and TTC12 were commonly enriched in flagellated sperm motility, sperm motility, and cilium movement involved in cell motility (Supplementary Figure 1).

To better understand how lncRNAs act in both cis and trans manner, we predicted their targets. There were 471 mRNAs that were discovered as targets of 10 lncRNAs (five upregulated and five downregulated), using 30 kb as the cutoff (Supplementary Figure 2). According to the GO enrichment study results, 55 significant GO terms were found (corrected P-value 0.05). In MF, the top five GO keywords were nucleoside phosphate binding, protein binding, small-molecule binding, binding, and nucleotide binding (Supplementary Figure 3). In BP, negative regulation of the cellular metabolic process, system development, response to chemical stimuli, organonitrogen compound metabolic process, and multicellular organism development were among the top five GO terms and cytoplasm, cell junction, cytosol, nucleoplasm, and synapse in CC (Figure 5B).