Micrococcal nuclease sequencing of porcine sperm suggests enriched co-location between retained histones and genomic regions related to semen quality and early embryo development

Sample collection

Two ejaculates, each from a different healthy Pietrain boar, MN_1 (16 months of age) and MN_2 (9 months old), from two artificial insemination centres were obtained by specialists during their routine sample collection using the gloved hand method (King & Macpherson, 1973) and were immediately diluted (1:2) in commercial extender. Both ejaculates showed good semen quality as evidenced by the evaluation of several sperm parameters. The percentages of sperm cell viability, structurally altered acrosomes and morphological abnormalities were measured by staining the samples with the eosin-nigrosin technique after 5 min incubation at 37°C. These analyses showed 94.3% and 92.2% of viable cells for sample MN_1 and MN_2 respectively; 96.7% and 96.3% of cells with normal acrosomes, respectively; 2.2% and 1.1%, 3.1% and 1.1% and 2.2% and 4.6%, of head, neck and tail abnormalities, respectively. Semen samples were purified with BoviPure^TM (Nidacon, Moelndal, Sweden) to remove potentially present somatic cells and processed as described in Gòdia et al. (2018). Briefly, a variable volume of sperm according to its concentration, with a maximum of 1 billion cells and not exceeding 11 mL, were layered over 3 mL of BoviPure^TM diluted to a final ratio of 60% (v/v) with BoviDilute^TM (Nidacon, Moelndal, Sweden) and centrifuged at 300 × g for 20 min at 20 °C with slow acceleration and deceleration rates. The resulting cell pellet was washed with PBS and resuspended in 1 mL of RNase-free PBS and further pelleted by centrifugation. The pellets were stored at −80 °C until further use.

Micrococcal nuclease digestion, library construction and sequencing

Chromatin digestion was performed as previously described (Zalenskaya, Bradbury & Zalensky, 2000) with minor modifications. Forty million purified spermatozoa cells were pelleted by centrifugation (3,500 rpm for 10 min) and washed twice with 1 mL of 1x phosphate buffered saline (PBS), and 1 mM phenylmethylsulfonyl fluoride (PMSF). Pelleted cells were resuspended with 1 mL of 0.1% Lysolecithin (dissolved in 1x PBS and 1 mM PMSF), incubated on ice for 10 min and centrifuged (3,500 rpm for 10 min). The sperm head was decondensed in 1 mL of 1x PBS, 1 mM PMSF with 10mM Dithiothreitol (DTT) and incubated on ice for 12 min. After centrifugation (2,000 rpm for 10 min), nuclei were washed with 1 mL of 1x PBS and 2mM DTT and centrifuged again (3,500 rpm for 10 min). The pellet was resuspended in 100 µl of 1x PBS, 1 mM DTT and CaCl₂ was added to a final concentration of 0.6 mM. At this time, 5 U of MNase (Sigma-Aldrich) were added and the samples were incubated for 7 min at 37 °C. Digestion was stopped by the addition of ethylene glycol-bis (β-aminoethyl ether)-N,N,N′,N′-tetraacetic acid (EGTA) to a concentration of 5 mM. The resulting digested fraction corresponding to nucleosome-bound chromatin (soluble) was separated from the undigested protamine-bound fraction (pellet) by centrifugation at 20,000 rpm for 10 min. An aliquot of the digested DNA fraction (which included both the MN and SN fractions) was treated with sodium dodecyl sulfate (SDS) and proteinase K and evaluated in a 1.5% agarose 1x tris-acetate-EDTA (TAE) gel electrophoresis. Another aliquot of the nuclease digested DNA was used to extract DNA with a phenol-chloroform based protocol and was directly used for library prep after purification with the Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA). The MN and SN fractions were separated after sequencing and read mapping, using bioinformatics tools as described below. Purified DNA was subjected to quality control including quantification with the Qubit^TM DNA HS Assay kit (Invitrogen, Waltham, MA, USA) and Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA) as well as size and concentration assessment with a 2100 Bioanalyzer and the Agilent DNA 1000 Kit (Agilent Technologies, Santa Clara, CA, USA). gDNA from the two sperm samples was also extracted as done by Hammoud and colleagues (Hammoud et al., 2009). Briefly, sperm cells were thawed and pelleted at 10,000 × g for 1 min. The supernatant was removed and the pellet was resuspended in 1 mL Lysis Reagent (1% Triton X-100, 5 mM MgCl2, 320 mM Sucrose and 10 mM Tris pH 7.5). After centrifugation (20,000 × g for 5 min), cells were resuspended again with 0.5mL Lysis Reagent two times, centrifuged and the supernatant was discarded. The sperm pellet was dispersed with 305.5µl of Enzyme Master Mix (1.25 mM MgCl₂, 1.25 mM Deferoxamine Mesylate, 12.5 mM Tris pH 8.0). Then, 1.6 µl RNase A-20 mg/mL -, 1.6 µl of 1 M DTT and 6.5 µl of 5M NaCl was added to each sample. This cell suspension was kept at 37 °C for 30 min. The cells were then dispersed with gentle vertexing and subsequently, 40 µl of 10% SDS and 20 µl of Proteinase K (20 mg/µl) were added. Samples were then incubated at 50 °C for 60 min. Following this step, 1 mL of 100% Isopropanol was added and the tubes were inverted several times until a gel-like DNA precipitate appeared. Then, 0.6 mL of NaI solution was added (40 mM Tris pH 8.5, 20 mM EDTA and 7.6 M NaI), and the tubes were inverted until DNA formed a white precipitate. The DNA pellet was washed three times with 70% EtOH with a final wash with 100% EtOH. The DNA pellet was solubilized with 300 µl of H₂O. The two resulting gDNA samples were pooled to be used as input. Pooled gDNA were sheared to obtain 100 bp long fragments using a Covaris S2 instrument (Covaris Inc, Woburn, MA, USA), according to the manufacturer’s instructions. Sequencing libraries of the two MNase treated samples and the input gDNAs were prepared with the PrepX^™ DNA Library Kit (Takara, Kusatsu, Japan). The libraries were used as a template to generate 50 bp paired end (PE) reads in an Illumina HiSeq2500 system.

MNase-Seq data preprocessing and quality evaluation

Raw sequencing data was filtered to remove low quality bases, adaptor sequences and short reads (<25 bp length) with Trimmomatic v.0.36 (Bolger, Lohse & Usadel, 2014). Trimmed reads were aligned to the porcine genome (Sscrofa11.1) with Bowtie2 v.2.4.1 (Langmead & Salzberg, 2012) fitting default parameters except “—very sensitive”. Pearson correlation of the read coverages for genomic regions along the genome was calculated using the multiBamSummary file based on the two MNase samples with the tool deepTools v.3.3.2 (Ramírez et al., 2016) with options: “plotCorrelation—removeOutliers–skipZeros”. Since the correlation between the samples was very high, the reads from the two replicates (MN_1: replicate A; MN_2: replicate B) were merged and processed as a pool with the aim to increase sequencing depth and peak detection power and accuracy. Then, the MN and SN fractions were bioinformatically separated based on the genomic distance between the two paired reads. The reads from the SN fraction, with a length slightly below 100 bp, were extracted using the parameter “--maxins 110”, which specifies that the mapped paired-end reads should not exceed 110 bp. To obtain the reads from the MN fraction, with a length around 150 bp, the parameters “—minins 111” and “—maxins 1000” were applied. Subsequently, duplicate reads of these fractions and the input sample were removed with Picard-Tools MarkDuplicates v.1.56 (http://picard.sourceforge.net). Finally, for both fractions, peak calling was performed with MACS2 v.2.1.0 (Zhang et al., 2008) with “-q 0.05 -B -g hs” and compared against the input sample. Visualization of the genomic heatmaps from the MNase-Seq signals using transcription start sites (TSS) as reference points was carried with the deepTools (Ramírez et al., 2016) computeMatrix tool with parameters: “reference-point -b 500 -a 500 –skipZeros” and plotHeatmap tool using “—refPointLabel ‘TSS’”.

Peak location relative to gene annotation

Peaks were categorized by their position in relation to gene features as annotated in Ensembl (v96) using BEDtools v.2.29.2 (Quinlan & Hall, 2010) intersect. The peaks were classified as mapping to TSSs, promoters, overlapping 5′ untranslated region (UTR), 3′ UTR, exonic, intronic and intergenic as carried in Halstead et al. (2020).

In order to determine any potential positional preference or avoidance of the peaks within the gene features, 1,000 iterative permutations of the peak genomic locations with the same number of peaks following the same peak width distribution but with randomized genomic locations, were carried using BEDtools (Quinlan & Hall, 2010) shuffle. Subsequently, genomic features were assigned to the randomized peaks based on their genomic position within these features. These results were subsequently compared against our results on real data using a permutation test. This test calculates the proportion of permutations in which the number of peaks mapping to a specific feature deviates from the mean of the distribution more than the number of overlaps observed in the real data. In other words, it determines whether the peaks in the permutations are located more towards the tails of the distribution compared to the real data. The fold change (FC) was also measured. FC was calculated by dividing the number of overlaps between the real MN or SN peaks and each of the genomic features by the equivalent mean of the distribution of the number of overlaps obtained with the shuffled data. Comparisons resulting in P-value < 0.001 and FC ≥ 1.5 or FC ≤ 0.66 were classified as enriched or depleted, respectively.

Gene ontology analysis of the genes near MNase peaks

The genes mapping within or less than 5 kbp apart from the identified peaks were extracted from Ensembl v96 with BEDtools (Quinlan & Hall, 2010) closestBed and were used for Gene Ontology (GO) enrichment analysis. GO was carried out with Cytoscape v.3.8.2 plugin ClueGO v.2.5.7 (Bindea et al., 2009) using the Cytoscape’s porcine dataset and default settings. Only the significant Bonferroni corrected P-values were considered.

Motif enrichment analysis at the MN and SN underlying sequences

Motif enrichment analysis was carried out with the software HOMER v.4.10.0 (Heinz et al., 2010) findMotifsGenome.pl function with default parameters.These default parameters include: (i) random selection of the genome’s background regions to compare against the input sequence data and using the background frequency of nucleotides in the human genome, hg19; (ii) the default motif length for the search is 8-12 nucleotides long; (iii) motif search on both the forward and reverse strands of the input sequences.

Positional conservation of chromatin-associated DNA with human and cattle sperm

A comparable human sperm MNase-Seq dataset that sequenced the mono-nucleosomal fraction (GSE15690) from Hammoud et al. (2009), obtained from a pool of four donors was downloaded from the Gene Expression Omnibus database. The genomic coordinates for the human sperm MN peaks were liftover to Sscrofa11.1 using the UCSC liftover tool (Kuhn, Haussler & Kent, 2013) with default parameters, except for the parameter “-minMatch = 0.1”. To evaluate enrichment, the genomic locations of the porcine MNase peaks were randomized 1,000 times using BEDtools (Quinlan & Hall, 2010) shuffle and overlapped to the MNase peaks of human or cattle using BEDtools v.2.29.2 (Quinlan & Hall, 2010) intersect. Thereafter, the results were contrasted to the overlap of the real data with the permutation test that calculates the proportion of permutations in which the number of overlaped peaks lies further toward the tails of the distribution compared to the observed value in the real data. FC was also calculated by dividing the number of overlaps between the real MN or SN peaks and human MNase peaks by the mean of the distribution of the number of overlaps obtained from the shuffled data. To determine co-location enrichment or depletion, a threshold of P-value < 0.001 and FC ≥ 1.5 or FC ≤ 0.66, was set to determine enrichment and depletion, respectively.

The same approach was used for two cattle (GSE47843) comparable MNase-Seq datasets that sequenced the mono-nucleosomal fractions from two bulls (samples s_1_1_bovine and s_5_1_bovine) from the Holstein breed (Samans et al., 2014). Genome coordinates from bull (busTau7) genomes were first liftover to hg19, and then from hg19 to Sscr11.

Integration with other -omics data

Although mature sperm is assumed to be transcriptionally inactive, RNAs in sperm may both reflect preceding events in transcriptionally active spermatogenic precursor cells during male germ cell development and have a role in early development after fertilization (Gòdia, Swanson & Krawetz, 2018; Ostermeier et al., 2004). Thus, the existence of a relationship between the location of retained nucleosomes in sperm and the sperm transcriptome profile was hypothesized. To test this hypothesis, sperm transcriptome data from 40 healthy Pietrain boars belonging to commercial artificial insemination farms, including 40 total and 35 small RNA-seq datasets previously generated by our group (NCBI’s BioProject PRJNA520978) was used. These datasets provided a list of 4,120 protein coding RNAs (Gòdia et al., 2020b), 1,574 circular RNAs (circRNAs) mapping in pig chromosomes (Gòdia et al., 2020a) and 6,729 PIWI-interacting RNAs (piRNAs) (Ablondi et al., 2021). The averaged RNA abundances from all the samples were used for further analysis. Protein coding genes were classified according to their RNA abundances measured in Fragments per Kilobase per Million mapped reads (FPKM) as (i) absent (<1 FPKM); (ii) low abundance (≥1 to <10 FPKM); (ii) intermediate abundance (≥10 to <100 FPKM) and (iii) high abundance (≥100 FPKM). Then, the existence of positional co-location enrichment between each of these gene fractions and the MNase peaks was assessed using the Fisher Exact Test (two-tailed). Any distance below 5 kb between an RNA and a MN or SN peak was considered as co-location. Moreover, the RNA abundance of the genes that co-located with MN and SN peaks with the RNA levels of the whole set of genes annotated in the pig genome was compared by employing the Kruskal-Wallis test using RNA abundances stabilized with the log2 transformation.

The positional enrichment of circRNAs or piRNAs at MNase peaks was also studied. First, the genomic regions of the MN and SN peaks were iteratively randomized 1,000 times using BEDtools (Quinlan & Hall, 2010) shuffle. The overlap in the real data and the permutations with the circRNAs or piRNAs was determined with BEDtools v.2.29.2 (Quinlan & Hall, 2010) intersect. Then, the enrichment of the genomic coordinates of the real MNase peaks at piRNA and circRNAs was assessed by comparison with the overlap of the randomized MNase peak locations using the permutation test. The FC was also calculated by dividing the number of overlaps between the real MN or SN peaks and circRNAs or piRNAs by the mean of the distribution of the number of overlaps obtained with the shuffled data. A threshold of P-value < 0.001 and FC ≥ 1.5 (enrichment) or FC ≤ 0.66 (depletion), respectively, was set to consider co-location enrichment or depletion.

MNase profiles were also evaluated against the genomic regions showing genetic association with sperm quality traits in Pietrain pigs. These associations were identified through a Genome-Wide Association Study (GWAS) conducted by our group, which included the two samples used in this study (Gòdia et al., 2020b). This GWAS provided 18 regions showing genetic association with the percentage of sperm cells with head abnormalities (seven regions), neck abnormalities (four regions), percentage of abnormal acrosomes after 5 min incubation at 37 °C (two regions), the ratio of abnormal acrosomes after 5 min and 90 min of incubation (one region), percentage of motile spermatozoa after 5 min incubation (two regions) and after 90 min incubation (two regions, one shared with motility after 5 min incubation) and the percentage of sperm cells with proximal droplets (one region). These regions spanned in total, 16.5 Mbp (0.7%) of the porcine autosomal genome (Gòdia et al., 2020b). Again, to evaluate whether MNase peaks were enriched at GWAS hits, the genomic regions of the MN and SN peaks were iteratively randomized 1,000 times using BEDtools (Quinlan & Hall, 2010) shuffle and the overlap of both, the real observed data and the 1,000 randomized files with the GWAS hits was evaluated with BEDtools v.2.29.2 (Quinlan & Hall, 2010) intersect. This randomized overlap was compared to the overlap of the real data using the permutation test. The FC was also calculated by dividing the number of overlaps between the real MN or SN peaks and the GWAS hits by the mean of the distribution of the number of overlaps obtained with the shuffled data. A threshold of P-value < 0.001 and FC = 1.5 (enrichment) or FC = 0.66 (depletion), respectively, was set to consider co-location enrichment or depletion.

Peak location relative to repeat elements

The list of repeat elements in the pig genome from RepeatMasker was downloaded from the UCSC genome browser. The same approach as described for the other genomic features was followed to evaluate enrichment or depletion of MNase peaks at the whole set of repeat elements, long interspersed nuclear elements (LINEs) or short interspersed nuclear elements (SINEs). The genomic locations of the MN and the SN peaks were randomized 1,000 times and the overlap with the different sets of repeat elements was assessed with BEDtools v.2.29.2 (Quinlan & Hall, 2010) intersect. A permutation test was carried to evaluate enrichment (P-value < 0.001; FC ≥ 1.5) or depletion (P-value < 0.001; FC ≤ 0.66) of MN and SN peaks at repeat elements.

All scripts, including R, bash and software parameters used are available on Figshare: https://doi.org/10.6084/m9.figshare.21997523.v1.

MNase digestion, library preparation and data preprocessing

The MNase protocol yielded 95.2 and 75.6 ng of DNA for sample 1 and sample 2, respectively. The extraction of gDNA provided 1227.6 ng/ul. The MNase treatment of the pig sperm samples generated ~147 bp fragments corresponding to the MN fraction and an additional ~100 bp fragment corresponding to the SN fraction (Fig. S1).

The sequencing of the libraries yielded more than 43 million PE reads for each MNase biological replicate and 41.8 million PE reads for the input sample (Table 1). In average, 90.2% of the MNase reads and 92.7% for the input reads mapped to the porcine genome (Sscrofa11.1). The genome-wide profiles of MNase sensitivity of the two samples showed high correlation (Pearson R = 0.87) (Fig. S2) and similar mapping statistics (Table 1). Thus, the sample pooling strategy carried by Hammoud and collaborators (Hammoud et al., 2009) was replicated and the two samples were pooled with the aim to increase read depth and peak call sensitivity and accuracy. The final number of PE reads in each, the MN and SN fractions, was 62.1 and 7.1 million, respectively (Table 1).

Characterization of the MN and SN peaks

The analysis yielded 25,293 MN and 4,239 SN peaks (Tables S1 and S2). The MN peaks averaged 270 bp in width and covered 0.3% of the porcine genome. The SN peaks were 141 bp wide and covered 0.02% of the genome. Most MN (70.2%) and SN (94.5%) peaks were annotated in intronic and intergenic regions (Fig. 1). MN peaks were significantly enriched at promoters (P-value < 0.001, FC = 5.05), TSS (P-value < 0.001, FC = 11.26) (Fig. 2), 5′-UTR (P-value < 0.001, FC = 12.20), coding sequences (P-value < 0.001, FC = 4.16) and 3′-UTR (P-value < 0.001, FC = 1.56) (Table 2 and Fig. S3). This result was not replicated in the SN peaks, which showed mild depletion at 3’-UTR sites (P-value < 0.001, FC = 0.47) (Table 2 and Fig. S3). A total of 9,128 and 1,688 genes overlapped or localized less than 5 kbp apart from a MN and a SN peak, respectively (Table S3). The two MN and SN fractions co-existed in 1,145 genes (Table S3), which corresponds to 75% of the SN-associated genes.

To get an idea of the potential function of the nucleosome associated DNA, a GO analysis was conducted on the genes that overlapped with or were located within 5 kb of these peaks. For the MN peaks, the most significantly enriched terms were related to chemical perception, development (including the nervous system, multicellular organism, tissue, embryo, etc.) and cell and organ morphogenesis (Table S4). The genes located in the SN fractions were enriched for less GO terms and showed weaker significances and included nervous system process, cell projection, regulation of GTPase activity, organelle organization and Golgi vesicle transport (Table S5). To delve further into the potential functions of the nucleosome retention in sperm, the genomic sequences underlying the peak regions were subjected to sequence motif enrichment analysis. Enrichment for motifs from 55 plant and 57 animal (mostly mammalian: mouse and human) transcription factors (TFs) was identified in the MN peaks (Table S6). The animal motifs involved 15 TFs from the bHLH class (MyoD, MyoG, Myf5, NeuroG2, Olig2, Tcf21, Twist2, etc.) and other TFs related to embryo development and implantation (Atf4, CHOP, Erra, EBF1, E2F2, Foxh1, HOXA1, HOXA2, MAFK, Pax8, IRF1, IRF3, PR, Rfx1, RORg, RUNX2, Zic, Zic3). They also included TFs related to spermatogenesis or sperm function which have also been related to embryo development such as CUX1, CTCFL, PAX5 and Smad2. Finally, this list also contained Znf263 and FOXA1, two TFs that provide a direct link between the paternal gamete and the embryo (Table S6). One study identified 514 genes with paternal preferential expression during human early embryo development and proposed ZNF263 as the strongest candidate in regulating the expression of these genes (Leng et al., 2019). The pig genes that map less than 500 bp apart from a MN peak harboring a Znf263 motif were identified, and their overlap with the 514 human genes showing paternal preferential expression in the early embryo was determined. 3,264 genes mapped less than 500 bp away from a MN peak with Znf263 binding motif in pig sperm. These corresponded to 3,288 human orthologs, 92 of which (18% of 514) also showed preferential paternal expression during human early embryo development (Leng et al., 2019) (Table S7). This concordance almost reached significance, with a P-value of 0.052 as calculated using the hypergeometric test. The SN peaks presented enriched motifs for 28 animal and 46 plant TFs. Thirteen of the animal motifs enriched in the SN peaks were also enriched in the MN fraction (Table S8). The common TFs included some of the bHLH class (NeuroG2, Olig2 and Twist2), as well as TFs involved in embryo development and implantation (Atf4, CHOP, ERRA, Foxh1, HOXA1, HOXA2, IRF1 and Zic). The SN-specific enriched motifs corresponded to several TFs related to embryo development (NFX, Duxbl, PBX1, Pitx1, SF1, NFIL3, Gfi1b, HLF), meiosis (DMC1) and parent-of-origin expression driven by imprinting in embryo (Zfp57).

A comparable human sperm MNase-Seq dataset, which included pool data from four donors (Hammoud et al., 2009) and focused on MN peaks, was employed to assess the inter-species conservation of the genomic location of our MNase sites. Results in human were similar to these detected in our analysis in pig. 25,122 peaks were annotated in the human MN fraction, 20,641 of which were successfully converted to the pig genome coordinates. Twenty-six percent (5,395 peaks) of the human nucleosome associated sites overlapped with our porcine MN peaks. This is a highly significant co-location between the two species (P-value < 0.001, FC = 25.01), compared to randomization (Table 2 and Fig. S3). The two cattle samples contained 2,256 (s_1_1_bovine) and 8,446 (s_5_1_bovine) MN peaks. Of these, 1,011 (s_1_1_bovine) and 4,186 (s_5_1_bovine) were successfully liftover to pig and used to determine genomic overlap. Only 7 (s_1_1_bovine) and 16 (s_5_1_bovine) cattle MN peaks overlapped with the porcine MN peaks and no significant enrichment or depletion was observed.

To identify positional relationships between the nucleosome-associated regions and transcriptional activity, the locations of the MN and SN peaks were compared with the repertoire of RNAs that our group identified in porcine sperm (Gòdia et al., 2020b). A total of 12,125 protein coding genes were absent in sperm. On the contrary, 5,814, 3,521 and 598 genes were classified as being at low, moderate and high abundance in sperm, respectively. The location of the genes that are present in spermatozoa was, regardless of their abundance, highly enriched (low abundance: P-value = 1.8E−72, moderate abundance: Pvalue = 1.5E−91, high abundance: P-value = 1.1E−11) within or near the MN peaks when compared to the catalog of absent genes (Table 3). Similarly, the SN peaks were also enriched for the low (P-value = 5.0E−17), moderate (P-value = 3.6E−47) and highly abundant (Pvalue = 3.0E−5) genes, when compared to the set of absent genes (Table 3). In line with these results, the average RNA abundance of the genes mapping to both MN (Pvalue = 4.9E−19) and SN (P-value = 4.3E−22) peaks was significantly higher than the average abundance of all the genes annotated in the pig genome (Fig. 3).

Recently, our group annotated 6,729 sperm piRNAs in porcine sperm (Ablondi et al., 2021). The piRNA regions showed a significantly higher frequency of mapping within MN peaks (P-value < 0.001, FC = 4.47) and exhibited an overlap of 65 MN peaks (Table 2 and Fig. S3). A similar enriched co-location was also observed ((P-value < 0.001, FC = 4.70, overlap with 18 MN peaks) for the set of piRNA genomic regions involving 1,355 piRNAs (Table 2 and Fig. S3), which abundance in sperm was also found to correlate with various sperm phenotypes. These phenotypes include the percentage of motile sperm, the percentage of morphological abnormalities or the percentage of viable cells (Ablondi et al., 2021). No obvious enrichment or depletion of SN at piRNAs sites was observed (Table 2 and Fig. S3). No positive or negative co-location was observed between MN or SN peaks and the catalog of 1,574 sperm circRNAs and 148 circRNAs associated to sperm motility in swine (Gòdia et al., 2020a) (Table 2 and Fig. S3).

Our group recently published a GWAS for sperm quality traits in Pietrain boars (Gòdia et al., 2020b) which was used to evaluate the co-location of MNase sites with regions associated to these traits that could suggest a potential link between nucleosome retention and sperm quality. No significant co-location between MN or SN and the GWAS hits was found (Table 2 and Fig. S3).

Finally, the MN peaks showed depletion at LINEs (P-value < 0.001, FC = 0.48), while the SN sites presented enrichment at SINEs (P-value < 0.001, FC = 1.85) (Table 2 and Fig. S3).