Identification of circular RNAs in porcine sperm and evaluation of their relation to sperm motility

Circular RNAs (circRNAs) are emerging as a novel class of noncoding RNAs which potential role as gene regulators is quickly gaining interest. circRNAs have been studied in different tissues and cell types across several animal species. However, a thorough characterization of the circRNAome in ejaculated sperm remains unexplored. In this study, we profiled the sperm circRNA catalogue using 40 porcine ejaculates. A complex population of 1,598 circRNAs was shared in at least 30 of the 40 samples. Generally speaking, the predicted circRNAs presented low abundances and were tissue-specific. Around 80% of the circRNAs identified in the boar sperm were reported as novel. Results from abundance correlation between circRNAs and miRNAs together with the prediction of microRNA (miRNA) target sites in circRNAs suggested that circRNAs may act as miRNA sponges. Moreover, we found significant correlations between the abundance of 148 exonic circRNAs and sperm motility parameters. Two of these correlations, involving ssc_circ_1458 and ssc_circ_1321, were confirmed by RT-qPCR using 36 additional samples with extreme and opposite sperm motility values. Our study provides a thorough characterization of circRNAs in sperm and suggests that circRNAs hold potential as noninvasive biomarkers for sperm quality and male fertility.


Results
Phenotypic data. The Table S1). We obtained an average of 40.7 M total RNA-seq sequencing reads per sample, 98.2% of which passed the quality control filters. After read mapping, the unmapped reads were used for circRNA prediction (Supplementary Table S1). For small RNA-seq, we obtained an average of 7.3 M reads per sample. The vast majority of these reads (99.2%) passed quality controls and 81.5% mapped to the porcine genome (Supplementary  Table S1). We identified 95 miRNAs from the list of 306 that are annotated in swine.
Only a small fraction of genes produced more than one circRNA. These were considered circRNA hotspots. The circRNAs from hotspot genes are often produced from different back-splicing events of one exon with several others 21 . We detected 12 genes with 5 or more circRNA isoforms (Table 2). Again, four of these genes, namely TESK2, SPATA19, PTK2 and SLC5A10 have been directly related to sperm function and fertility.
The boar sperm has a highly specific circRNAome. We compared our circRNA catalogue with equivalent available datasets from other studies in human (including different brain sections and cell lines) 22 , mice   Table 2. circRNA hotspot genes in swine sperm. 12 genes harbored 5 or more exonic circRNAs.
(several brain segments, cell types and embryonic stem cells) 22 and swine (lung, skeletal muscle, fat, heart, liver, spleen, kidney, ovary, testis and 5 brain sections) 21,23 . Twenty-four % and 11.3% of the boar sperm circRNAs had potential orthologs in human and mouse, respectively (Table 3). On the other hand, 20.3% of the porcine sperm circRNAs were also present in other porcine tissues 21,23 (Table 3). The porcine tissues showing higher overlap with sperm were testes (11.6%) and brain cortex (11.3%) ( Table 3). Comparing this in the opposite direction, 4.9% (6,636 circRNAs) of the circRNAs already annotated in any pig tissue were also present in our pig sperm list. This value was 12 times lower (0.4%) when evaluating the much bigger (90,067 circRNAs) human circRNA catalog (Table 3).
Sperm circRNAs do not follow an age-dependent pattern. We assessed whether sperm circR-NAs depicted an age-accumulating profile as it has been previously observed in rat testes 24 , and rat and mouse brains 24,25 . All the boars from our dataset were sexually mature with ages ranging between 9 and 54 months of age. Sexual maturity in boars is a process that starts at the age of 8 months and finalizes when animals reach 2 years. Thus, we divided the samples in those coming from boars approaching sexual maturity with ages below 2 years old and those produced by mature pigs with ages above 2 years old. There was no significant difference in the number of circRNAs identified (P-value: 0.68, Wilcoxon rank sum test) nor in their RNA abundance (P-value: 0.948) between the two groups. We repeated the analysis considering only extreme ages: young (N = 4; between 8.6 and 9.2 month old) and mature (N = 4; 29.9 to 54.6 month old). Again, there was no difference in the number of circRNAs identified (P-value: 0.89) or in their abundance (P-value: 0.2).

circRNA-miRNA interaction network.
To evaluate the potential role of exonic circRNAs as miRNA sponges, we built a circRNA:miRNA co-expression network with the RNA abundances of the 1,261 exonic circR-NAs (from total RNA-seq libraries) and the 95 miRNAs (from short RNA-seq libraries). The analysis identified 1,882 significant correlations that included 95 miRNAs and 458 exonic circRNAs. In parallel, we also performed an in silico prediction of miRNA target sites involving the miRNAs and exonic circRNAs identified in our data.
We assessed miRNA target sites in the exonic circRNA sequences based on sequence complementarity using the tool miRanda 26 . The analysis yielded 4,987 potential targets involving all the 95 miRNAs and 1,103 circRNAs. To reduce the proportion of false-positives, only the 81 interactions (from 34 miRNAs and 65 circRNAs, each from a different gene) that were identified in both approaches were used for network visualization (Fig. 2). Fifty-three of these 65 circRNAs were predicted to interact with only one miRNA (Fig. 2). Twelve circRNA were predicted to interact with two or more miRNAs. The circRNAs displaying the largest number of predicted interactions with miRNAs were ssc_circ_0954 and ssc_circ_1454, which connected with four miRNAs each (Fig. 2). On the other hand, miR-28-5p and miR-26a were predicted to be regulated by 10 and 9 different circRNAs, respectively. Similarly, miR140-3p and miR-423-5p were potentially targeted by five different circRNAs each (Fig. 2).
Correlation of circRNAs with sperm motility and circRNA validation. A total of 148 exonic cir-cRNAs from 142 genes showed significant correlations between their abundance and sperm motility traits (Supplementary Table S4). More in detail, 24, 83, 24 and 41 circRNAs correlated with the percentage of total motile spermatozoa, VCL, VAP and VSL, respectively (Supplementary Table S4).  Table 3. Concordance between the circRNAs catalogue of the boar sperm and the circRNA list in other species and pig tissues. Number of circRNAs after Sscrofa11.1 liftover from human, mice and swine (detailed by tissue).

Species/Tissues
Column 'Pig sperm' shows the proportion of swine sperm circRNAs that were identified in the given species and tissues. Column 'Other species/tissues' of the total circRNAs lists the proportion of circRNAs from that tissue or species that found a homolog in the boar sperm.
Scientific RepoRtS | (2020) 10:7985 | https://doi.org/10.1038/s41598-020-64711-z www.nature.com/scientificreports www.nature.com/scientificreports/ To confirm the existence of these circRNAs and their phenotypic correlations we undertook a Reverse Transcription quantitative PCR (RT-qPCR) and a Sanger Sequencing approach. First, we randomly selected two circRNAs with different RNA abundance levels to test whether we could confirm the bioinformatically predicted circRNAs. The chosen circRNAs were: ssc_circ_1141 (from PTGES3) with 75.0 CPM, and ssc_circ_0670 (from BAZ2B) with 10.6 CPM, both also detected in human 22 (hsa_circ_0008137 and hsa_circ_0002463, respectively) and in swine testes 23 . The PCR amplification of the two circRNAs resulted in a single electrophoretic band of the expected size ( Supplementary Fig. S1a) and Sanger Sequencing confirmed the back-splice junction (BSJ) ( Supplementary Fig. S2a,b).
Then, we used the same approach to confirm the existence of 8 circRNAs selected for their significant correlation with at least one motility trait (Supplementary Table S4). These circRNAs were further selected based on the fact that: (i) they were present at high abundances, or (ii) they showed significant correlation with at least 1 trait, or (iii) they had been previously identified in pig 23  . ssc_circ_0839 from PAIP2 did not amplify (data not shown) and ssc_circ_0118 from PDE10A (also identified in pig testes and in human as hsa_circ_0078638) displayed 2 amplification bands ( Supplementary Fig. S1b). These two circRNAs were discarded for further analysis. Thus, we confirmed the existence of these six circRNAs.
The RT-qPCR levels of the 6 circRNAs were measured in 36 animals presenting extreme and opposite values of sperm motility (N = 18 for each phenotypic distribution tail) from a dataset of 300 boar ejaculates with phenotype records. None of the 36 samples was included in the RNA-seq study. The two sample groups displayed significant phenotypic differences for all the traits: percentage of total motile spermatozoa (P-value: 2.65 ×10 −9 , Wilcoxon rank sum test), VCL (P-value: 2.20 ×10 −10 ), VSL (P-value: 3.22 ×10 −7 ) and VAP (P-value: 3.22 ×10 −7 ). The RT-qPCR assays presented efficiencies between 99.6% and 105.2%.
Two of the six circRNAs showed significant differences between the 2 sperm motility groups (Fig. 3b). These 2 circRNAs were ssc_circ_1458 from LRBA (P-value: 0.049) and ssc_circ_1321 from PAPOLA (P-value: 0.035). A third circRNA, ssc_circ_1132, from LIN7A, showed significant differences between the two motility groups (P-value: 0.008) (Fig. 3b) but in the opposite direction than expected according to the RNA-seq data and was consequently considered as not validated due to inconclusive results. The other 3 circRNAs, ssc_circ_1101 from KHDRBS3, ssc_circ_0437 from ULK4 and ssc_circ_1061 from ZNHIT6 did not present significant differences between the 2 groups (Fig. 3b).

Discussion
This study provides the first RNA-seq based circRNA repertoire of mature spermatozoa in a mammalian species. We have also assessed the sperm circRNA role as miRNA sponges and evaluated their associations with sperm motility as a parameter of semen quality. Our results provide further evidence and expand the relevance of spermatozoa RNAs in sperm biology and quality, as recently described in other studies 18 . circRNAs have been www.nature.com/scientificreports www.nature.com/scientificreports/ reported as promising potential biomarkers for human health including cancer, diabetes, cardiovascular diseases or the pathogenesis of pre-eclampsia 15 . Here, we provide novel data supporting the involvement of circRNAs on the male's reproductive function as reflected by their association with several sperm motility parameters.
We have identified nearly 1,600 boar sperm circRNAs that are robustly present in most samples (at least 30). The sperm circRNA repertoire showed similar circular genomic characteristics including the proportion of genomic overlap with protein coding and intergenic regions, the number of exons and the circRNA length ( Fig. 1), when compared to data previously reported in other tissues from swine 23 , human 17 or rat 24 . Nonetheless, the list of boar sperm circRNAs showed a modest overlap with these other datasets (Table 3). Not surprisingly, the largest overlap was with porcine testes (11.6%), probably, due to the fact that the male gonads are mainly composed of cells from the spermatogenic lineage including spermatozoa. The reduced concordance between swine sperm and other tissues and species (Supplementary Table S5) may be partly influenced by technical differences involving the processing of samples and the data analysis between the studies. Moreover, opposite scenarios were identified in swine and human. In swine, there was high sperm-specificity, with nearly 80% of the circRNAs being present only in sperm (Table 3). In contrast, human sperm gave a positive signal for the majority (72%) of the 13,617 circRNAs from circBase 22 that were queried using microarray technology 18 . A biological explanation for this difference between swine and human cannot be rule out but it is unlikely. The mRNA and miRNA components and payload of the human and porcine sperm have been previously interrogated and are highly ressemblant 20 . For this reason, we hypothesize that the differences observed are merely technical and related to the processing of the samples and the data.
We investigated the functional relevance of sperm circRNAs under the hypothesis that their function is associated to the host gene. Four genes (ATP6V0A2, PPA2, PAIP2 and PAXIP1) harboring the 15 most abundant exonic circRNAs (Table 1)   PPA2 is located in the mitochondrial membrane and might be involved in the production of ATP, control of molecular processes linked to the launching of sperm capacitation and sperm motility 28 . PAIP2 is an inhibitor of translation and is linked to sperm maturation and male fertility 29 . Finally, PAXIP1 is critical for genome stability and chromatin condensation and is associated with developmental arrest of spermatocytes, testicular atrophy and infertility in knockout mice 30 . We also identified 12 hotspot genes producing five or more circRNAs each ( Table 2). Some of these hotspot genes were related to sperm function and fertility. TESK2 is a protein kinase mostly expressed in round spermatids and is predicted to play a role in early stages of spermatogenesis 31 . PTK2, is essential for embryo development 32 . SPATA19 is critical for sperm mitochondrial function in relation to sperm motility and fertilization ability 33 .
The ontology enrichment analysis of the genes harboring circRNAs pointed towards epigenetic related functions, which are essential for chromatin condensation in sperm and for the reprogramming of gene expression upon egg fertilization and during embryo development (Supplementary Table S3). Gene enrichment analysis also signaled towards spermatogenesis and developmental processes, also involving the embryo development related genes, DHX36, IPMK, RICTOR, CDC73 and ANGPT1 (Supplementary Table S3). These functions are in line with previous works studying sperm mRNAs 9,20 , thereby providing further basis for the regulatory role that circRNAs might have on their cognate linear mRNAs.
A previous study in rat testes identified a dynamic circRNA age-dependent pattern of expression and suggested a relation between their abundance and function with the male's sexual maturity and spermatogenesis 24 . For this reason, we sought to investigate whether, like in rat testes, circRNAs also accumulate through age in the porcine sperm. Our data suggested that there is no association between the mature sperm circRNAs and the boar's age. This could be explained by different scenarios. First, while the study on rat testes compared animals with broad age differences (between 2 and 104 weeks) 24 , all our samples were from adult pigs with less dissimilar ages (between 9 and 54 months) taking into account that swine has a longer life cycle. Second, the study in rats assessed the whole testis payload, which includes cells from the spermatogenic lineage and also, leukocytes, endothelial, Sertoli and Leydig cells. Our study in pigs exclusively queried mature-selected ejaculated spermatozoa. Cells with high proliferation rates 34 such as spermatogonial stem cells (SSC) have been suggested to accumulate less circRNAs due to passive thinning out during their continuous self-renewal and proliferation until they become spermatozoa. Spermatozoa originates from SCC through a highly constant and orchestrated process and upon maturation, they enter a silent state with no apparent transcriptional activity. This would imply that in spermatozoa, circRNAs do not accumulate through age. Interestingly, our results indicate that circRNAs in mature sperm are stable and thus, provide further support for the potential of circRNAs as noninvasive biomarkers for male reproductive conditions.
To shed light into the functional relevance of circRNAs as miRNA sponges we built an interaction network (Fig. 2). We combined RNA-seq benchwork data of circRNA and miRNA abundances and in silico prediction of miRNA target sites to increase the reliability of the circRNA-miRNA relationship predictions. This network of 81 interactions contained some circRNA genes and miRNAs previously related to semen quality and male fertility. Remarkably, 2 circRNAs, ssc_circ_0954 and ssc_circ_1454 were linked to 4 different miRNAs each. The first, ssc_circ_0954 arises from DCDC2C, a gene identified in the human sperm flagellum end-piece with a suggested role on microtubule dynamics by acting as a depolymerization/polymerization balancing system 35 . This circRNA was predicted to interact, among others, with miR-361-3p which was found to be dysregulated in subfertile men 36 , miR-423-5p, which was altered in oligozoospermic men 37 and miR-28-5p, a miRNA that was dysregulated in normozoospermic infertile individuals 38 . The other cirRNA, ssc_circ_1454, is transcribed from MTHFD2L, a mitochondrial isozyme from the folate cycle metabolic pathway, a vitamin that has been also related to semen quality and fertility in men 39 . One of the miRNAs targeted by ssc_circ_1454 was miR-16, which was dysregulated in subfertile men 36 . The network also showed 10 circRNAs that may be regulating miR-28, a miRNA (miR-28-5p) that was dysregulated in normozoospermic infertile individuals 38 . The potential miR-28 regulators included ssc_circ_1370, arisen from FAM92A, whose protein may play a role in ciliogenesis 40 , and ssc_circ_0002 from WDR7, which is associated to sperm quality in cattle 41 . Likewise, 9 circRNAs were predicted to regulate miR-26a, which has been in turn, linked to VCL, VSL and VAP motility parameters in a previous study in swine 42 . miR-26a was also targeted by ssc_circ_0361, a circRNA from ACTL6A, which is a gene with a crucial function for embryo development 43 and ssc_circ_1352 from CAGE1, an acrosomal protein with proposed roles in fertility 44 . Additional interesting interactions included ssc_circ_0345 from the hotspot gene SLC5A10 (sodium-dependent mannose and fructose transporter) ( Table 2), which regulated miR-423-5p, a miRNA that was upregulated in oligozoospermic semen 37 and let-7c, which was altered in patients with severe asthenozoozpermia 45 . Altogether, the network involved host genes and miRNAs that are connected to male fertility thereby suggesting that circRNAs may play a role in the male's reproductive ability.
The potential role of circRNAs in male reproduction is further substantiated by the associations detected between sperm motility and the abundance of some sperm circRNAs in our study and the work detailed by Chioccarelli and co-authors 18 . At least 20 of the circRNA host genes implicated in the phenotypic correlations have been previously linked to sperm biology or male fertility (Supplementary Table S4). For example, ssc_ circ_0823, which correlated with VCL (P-value = 0.009), is a circRNA hosted by CAMK4, a gene that has been implicated in sperm motility in humans 46 . ssc_circ_0780, correlated with the percentage of motile cells (P-value = 0.041), is hosted by LRGUK, a gene that is required for sperm assembly including the growth of the axonome, a key structure for the flagellar beating in sperm 47 . We successfully tested by RT-qPCR 6 of the 148 exonic circR-NAs that correlated to sperm motility. We validated their correlation with motility traits for 2 of these 6 circR-NAs, ssc_circ_1458 from LRBA and ssc_circ_1321 from PAPOLA. Both were significantly down-regulated in the ejaculates with low sperm motility (Fig. 3b). LRBA is a gene involved in coupling signal transduction and vesicle trafficking but no link with sperm function or fertility has been found thus far. ssc_circ_1321 from PAPOLA, has a human ortholog circRNA (hsa_circ_0033126) and the host gene is implicated in RNA and ATP binding. Thus, Scientific RepoRtS | (2020) 10:7985 | https://doi.org/10.1038/s41598-020-64711-z www.nature.com/scientificreports www.nature.com/scientificreports/ none of the two host genes have been previously related to sperm function or fertility. However, the association is confirmed by RT-qPCR and we cannot exclude unidentified relevant functions on sperm motility. Two other cir-cRNAs were also of high interest and tested as they were within the top 15 most abundant (Table 1) and correlated with sperm motility (Supplementary Table S4). They were ssc_circ_0839 from PAIP2, with crucial roles in spermatogenesis 29 and ssc_circ_1101 from KHDRBS3, a gene that is highly abundant in mice testes 48 . ssc_circ_0893 did not amplify when subjected to RT-qPCR and ssc_circ_1101 did not present significant differences in the RT-qPCR levels between the motility groups (Fig. 3b).
Interestingly, some circRNAs popped up as relevant in more than one of the analyses carried in this study. For example, ssc_circ_1532, from the SPATA19 hotspot gene (Table 2), which is related to sperm motility and fertility 33 , was suggested to regulate miR-99a according to the network analysis (Fig. 2). miR-99b was found to be dysregulated in the low motility sperm fraction in bull 49 and in subfertile men 36 . Another circRNA, ssc_circ_1219 from OSBPL9, a gene involved in male reproduction 50 , displayed abundance correlation with VCL (P-value: 0.04) (Supplementary Table S4) and was identified as a potential target of miR-101 (Fig. 2), a miRNA (miR-101-3p) that was altered in asthenozoospermia men 45 .
Remarkably, 4 (DENND1B, PTK2, SLC5A10 and CAMSAP1) of the 12 hotspot genes hosted a circRNA with significant abundance correlation with sperm motility. Thus, one third of the hotspot genes included circRNAs correlated with sperm motility whilst only 15.0% (148) of the 984 genes hosting the 1,598 circRNAs were correlated with motility (Supplementary Table S4). Noteworthy, the 12 circRNA hotspot genes were not hotspots in the other porcine tissues analyzed 21,23 (data not shown). Altogether, this indicates that circRNA hotspot genes may have relevant tissue-specific functions.
In conclusion, our study is the first to characterize the circRNAs present in porcine spermatozoa. We have provided a comprehensive view of the boar sperm circRNAome, which is highly sperm-specific and involves genes related to sperm biology and development and identified potential circRNAs acting as miRNA sponges. Moreover, we have detected and validated correlations between the abundance of some circRNAs and sperm motility parameters. The results described in this study may be of interest for animal breeding and for human health. On the one side, our data identified genomic regions of relevance for sperm quality that could harbor SNP markers for the early prediction of the piglets that will enter the AI studs. On the other hand, this study may spur novel research on the potential of circRNAs as noninvasive biomarkers for male fertility in human medicine. The 40 ejaculates were used for RNA-seq and 34 of them for short RNA-seq. The ages of these boars ranged between 9 and 54 months, with average and median of 15 and 11 months old, respectively. Twenty-four boars were less than 1 year old, 12 pigs were between 1 and 2 years old, 3 pigs between 2 and 3 years old (Supplementary  Table S6).

RNA isolation and library preparation.
Ejaculates were purified to remove somatic cells and immature sperm cells and purified sperm was stored at −96 °C with Trizol ® as described by Gòdia et al. 52 . Briefly, the ejaculates were suspended over a 3 ml cushion of the commercial solutions of BoviPure and BoviDiluteTM (Nidacon; Mölndal, Sweden), which includes colloidal silica particles coated with silane in an isotonic salt solution in a 15 ml RNAse-free tube. Then, the samples were subjected to centrifugation at 300 × g for 20 min at 20 °C. After centrifugation, all the upper phases were removed and the cell pellets were transferred to a new RNase-free 15 mL tube, washed with 10 mL of RNase-free PBS and centrifuged at 1,500 × g for 10 min at 20 °C. The supernatants were then removed and the pellets optically inspected using a microscope to confirm the removal of somatic cells. RNA was extracted from purified sperm cells, treated with TURBO DNA-free ™ Kit (Invitrogen) and quantified using Qubit TM RNA HS Assay kit (Invitrogen). We assessed RNA integrity with the 2100 Bioanalyzer and yield using the Agilent RNA 6000 Pico kit (Agilent Technologies). We then performed RT-qPCR assays for PRM1 and PTPRC mRNAs as well as for intergenic/genomic DNA to verify proper sperm purification.
Forty sperm RNA samples were subjected to total RNA-seq. Thirty-four of these samples were also used for small RNA-seq. For total RNA-seq libraries, ribosomal RNA (rRNA) was depleted with the Ribo-Zero Gold rRNA Removal Kit (Illumina) and libraries were constructed with the SMARTer Universal Low Input RNA library Prep kit v2 (Clontech). Resulting libraries were sequenced in a HiSeq2000/2500 system (Illumina) to generate 75 bp long paired-end reads. For small non coding RNA-seq libraries, extracted RNA without rRNA depletion was directly subjected to library preparation with the NEBNext Small RNA Library Prep Set kit (New England Biolabs) and sequenced on a HiSeq2000 (Illumina) to generate 50 bp single-end reads.

RNA-seq and bioinformatics analysis.
For the total RNA fraction, raw reads were filtered by removing adaptor sequences and low-quality reads with Trimmomatic v.0.36 53 . The identification of circRNAs was carried on these reads with the find_circ pipeline 14 with stricter filter stringency. To reduce the false positive rate in the discovery of circRNAs, we selected circular splice transcripts with at least two unique supporting reads Scientific RepoRtS | (2020) 10:7985 | https://doi.org/10.1038/s41598-020-64711-z www.nature.com/scientificreports www.nature.com/scientificreports/ in the anchor segment and with Phred quality scores of 35 or more. Moreover, only circRNAs predicted in at least 30 samples were kept. The RNA abundance of the predicted circRNAs were normalized as the number of BSJ spanning reads per million raw reads (CPM) 21 . The functional regions of circRNAs were identified based on their co-location with genomic features (e.g. exon, 3′UTR, 5′UTR, etc) from the Ensembl database (release 91) with BEDtools 54 . Our catalogue of boar sperm circRNAs was contrasted with other publically available porcine circRNA databases including heart, liver, spleen, lung, kidney, ovarium, testis, skeletal muscle, fat and fetal brains 21,23 . We also queried several human tissues (including several cell lines, brain sections placenta, muscle, fat, umbilical cord, atrium, decidua and plasma) and murine (cell lines and brain sections) available at the circBase database 22 . Genomic coordinates from the human and mouse circRNAs were liftover to Sscrofa11.1 using the UCSC liftover tool 55 .
For the small RNA-seq analysis, trimming of adaptors and low quality bases was performed with Cutadapt v1.0 56 . The mapping of sncRNAs was performed with the sRNAtoolbox v.6.17 57 with default settings and providing miRBase 58 release 21 as library dataset. Multi-adjusted read counts were then normalized by sequencing depth as CPM. We only considered the miRNAs that were detected >1 CPM in all the samples. circRNA-miRNA network visualization. To identify circRNA:miRNA interactions, we carried a Partial Correlation with Information Theory (PCIT) analysis 59 using the RNA abundance levels of the of exonic cir-cRNAs and miRNAs after stabilization with log2 transformation. Only negative RNA abundance correlations between circRNAs and miRNAs were kept. We further assessed circRNAs-miRNAs interactions using miRanda v.3.3a 26 to predict miRNA target sites in the circRNAs sequences. The potential interactions identified with both approaches were kept and visualized with Cytoscape v.3.7.0 60 .

Gene Ontology analysis and correlation with sperm motility parameters. GO analysis was carried
with PANTHER v.13.1 61 with the overrepresentation test and P-values corrected with FDR. Annotation Data Set was "GO biological process complete". Pearson correlation was used to determine associations between circRNA abundance levels and sperm motility parameters. P-values < 0.05 were considered statistically significant.
Validation of circRNAs and reverse transcription quantitative PCR (RT-qPCR). circRNAs were validated by Sanger Sequencing and quantified by RT-qPCR using divergent primers. Primers were designed using the Primer Express software (Applied Biosystems). Primer sequences are shown in Supplementary Table S7. For cDNA synthesis, 5 µl of RNA were reverse transcribed using the High Capacity cDNA Reverse Transcription kit in a final volume of 50 μL (Applied Biosystems) following the manufacturer's protocol. circRNAs were amplified and visualized in 3% high resolution agarose gel electrophoresis and confirmed by Sanger Sequencing.
The abundance level of 6 circRNAs, correlated to sperm motility parameters in the RNA-seq study, was analyzed by RT-qPCR in 36 samples, none of them included in the RNA-seq. These 36 samples belong to two groups with extreme and divergent values for sperm motility from a bank of 300 ejaculates with phenotypic records. Isolated RNA could not be treated with RNAse R due to extremely low RNA yield (see Supplementary  Table S1), sample availability and treatment optimization. Nevertheless, one of the circRNAs was designed in the BSJ (Supplementary Table S7) and amplified successfully; which is suggestive of the true circular structure of the RNA products identified by RNA-seq and validated by RT-PCR. Moreover, some reports suggest RNase R can introduce technical noise by degrading some circRNAs while some linear transcripts are resistant to digestion in their 3′ ends 62 .
Quantitative PCR reactions were performed in triplicate following the manufacturer's instructions. The final reaction (15 μL) included 7.5 μL SYBR Select Master Mix (Life Technologies -Thermo Fisher Scientific), 300 nM of each primer and 3.75 μL of cDNA 1:4 diluted on a QuantStudio 12K Flex Real-Time PCR System (Applied Biosystems). To evaluate the efficiency of the RT-qPCR assays, standard curves with 6 serial dilutions from a pool of sperm cDNA were generated. Thermal profile was set as follows: 50 °C for 2 min, 95 °C for 10 min and 40 cycles at 95 °C for 15 sec and 60 °C for 60 sec. Moreover, a melting profile (95 °C for 15 sec, 60 °C for 15 sec and a gradual increment of temperature with a ramp rate of 1% up to 95 °C) was programmed at the end of the RT-qPCR to assess the specificity of the reactions. The genes ISYNA2 and GRP137 were selected as endogenous controls following the stability values after a GeNorm pilot experiment. Their stability was determined considering a GeNorm M value <0.5. Relative expression values were calculated using the ThermoFisher Cloud software (Applied Biosystems) applying the 2 -ΔΔCt method. The same software was used to compare the biological groups. Significance was set at a P-value < 0.05.

Data availability
The datasets generated and/or analysed during the current study are available at NCBI's BioProject PRJNA520978. Sample phenotypes can be provided upon reasonable request.