ORIGINAL RESEARCH published: 19 July 2021 doi: 10.3389/fvets.2021.668158 Characterization of the Impact of Density Gradient Centrifugation on the Profile of the Pig Sperm Transcriptome by RNA-Seq Yu Lian 1† , Marta Gòdia 1†‡ , Anna Castello 1,2 , Joan Enric Rodriguez-Gil 3 , Sam Balasch 4 , Armand Sanchez 1,2 and Alex Clop 1,5* 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus UAB, Barcelona, Spain, 2 Unit of Edited by: Animal Science, Department of Animal and Food Science, Autonomous University of Barcelona, Barcelona, Spain, 3 Unit of Isabel Barranco, Animal Reproduction, Department of Animal Medicine and Surgery, Autonomous University of Barcelona, Barcelona, Spain, University of Bologna, Italy 4 Grup Gepork S.A., Barcelona, Spain, 5 Consejo Superior de Investigaciones Científicas, Barcelona, Spain Reviewed by: Jean Feugang, RNA-Seq data from human semen suggests that the study of the sperm transcriptome Mississippi State University, United States requires the previous elimination from the ejaculates of somatic cells carrying a larger Margarida Fardilha, load of RNA. Semen purification is also carried to study the sperm transcriptome in University of Aveiro, Portugal other species including swine and it is often done by density gradient centrifugation to *Correspondence: Alex Clop obtain viable spermatozoa from fresh ejaculates or artificial insemination doses, thereby alex.clop@cragenomica.es limiting the throughput and remoteness of the samples that can be processed in one † These authors have contributed study. The aim of this work was to evaluate the impact of purification with density equally to this work and share first gradient centrifugation by BoviPureTM on porcine sperm. Four boar ejaculates were authorship purified with BoviPureTM and their transcriptome sequenced by RNA-Seq was compared ‡ Present address: with the RNA-Seq profiles of their paired non-purified sample. Seven thousand five Marta Gòdia, Animal Breeding and Genomics, hundred and nineteen protein coding genes were identified. Correlation, cluster, and Wageningen University and Research, principal component analysis indicated high—although not complete—similarity between Wageningen, Netherlands the purified and the paired non-purified ejaculates. 372 genes displayed differentially Specialty section: abundant RNA levels between treatments. Most of these genes had lower abundances This article was submitted to after purification and were mostly related to translation, transcription and metabolic Animal Reproduction–Theriogenology, a section of the journal processes. We detected a significant change in the proportion of genes of epididymal Frontiers in Veterinary Science origin within the differentially abundant genes (1.3%) when compared with the catalog of Received: 15 February 2021 unaltered genes (0.2%). In contrast, the proportion of testis-specific genes was higher in Accepted: 22 June 2021 the group of unaltered genes (4%) when compared to the list of differentially abundant Published: 19 July 2021 genes (0%). No proportion differences were identified for prostate, white blood, lymph Citation: Lian Y, Gòdia M, Castello A, node, tonsil, duodenum, skeletal muscle, liver, and mammary gland. Altogether, these Rodriguez-Gil JE, Balasch S, results suggest that the purification impacts on the RNA levels of a small number of Sanchez A and Clop A (2021) Characterization of the Impact of genes which are most likely caused by the removal of epididymal epithelial cells but also Density Gradient Centrifugation on the premature germinal cells, immature or abnormal spermatozoa or seminal exosomes with Profile of the Pig Sperm Transcriptome a distinct load of RNAs. by RNA-Seq. Front. Vet. Sci. 8:668158. Keywords: sperm RNA, RNA-Seq, sperm purification, differentially abundant gene, somatic cell, germline cell, doi: 10.3389/fvets.2021.668158 exosome Frontiers in Veterinary Science | www.frontiersin.org 1 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq INTRODUCTION fresh commercial extender and stored at 16◦ C. The ejaculates were privately owned for non-research purposes and the Despite being a matter of debate for many years, the presence and owners provided consent for their use for research. Specialized role of sperm RNA is beginning to be elucidated. The analysis of professionals at the farm collected the ejaculates following the spermatozoon transcriptome can shed light on the previous, standard routine procedures and guidelines. present or future roles of the genes involved in spermatogenesis, early embryo development, and transgenerational inheritance Sample Processing and RNA Extraction (1). The sperm’s transcriptomic landscape has been profiled For each ejaculate, we collected two aliquots. One aliquot was in a large number of animal species including human, mice, purified (P) with the BoviPureTM (Nidacon; Mölndal, Sweden) cattle, pig, and horse (1), and several research groups focus their colloid centrifugation method. The other aliquot was not purified efforts to identify molecular markers associated to semen quality (NP). For the P samples, a maximum of 11 mL of ejaculate and and fertility in human (2), cattle (3), horse (4) and swine (5– 1 billion cells were placed over 3 mL of BoviPureTM in 15 mL 10). Transcriptomic evaluations of the ejaculated spermatozoa RNase-free tubes. The tubes were then centrifuged at 300 × g for are often preceded by a purification step to remove somatic 20 min at 20◦ C. After centrifugation, all the upper phases were and prokaryotic cells often by using a commercial colloidal removed. The cell pellet were transferred to a new 15 mL RNase- silica suspension (e.g., PureSperm R , BoviPureTM , Nidacon, free tube with 10 mL RNase-free PBS and centrifuged at 1,500 Sweden) in order to prepare density gradients and purify the × g for 10 min at 20◦ C. The supernatant was removed and the sperm cells (11–14). This step is considered to be necessary pellets were stored at −80◦ C in 1 mL Trizol R for further use since spermatozoa carry small amounts of RNA that are also for RNA extraction. A detailed description of the purification highly fragmented when compared with the somatic cells that protocol with BoviPureTM is provided by Gòdia et al. (5). For could be present in the ejaculate. The sperm transcriptome the NP aliquots, 1.5 mL of ejaculate was centrifuged at 14,000 g profile could be thus overshadowed by the profiles from these for 4 min at 20◦ C and the supernatant was removed. The pellet somatic cells. However, the purification step with gradient was then eluted in 1 mL of Trizol R and stored at −80◦ C until cell centrifugation requires viable spermatozoa only available RNA extraction. in fresh ejaculates or in frozen straws prepared for artificial Total RNA was extracted from the eight samples as described insemination. The use of fresh ejaculates limits the number of in Gòdia et al. (5). Briefly, between 47 and 300 million cells samples that can be simultaneously processed in an experiment were pre-lysed using a 5 mL sterile syringe with a 25 G needle and the remoteness of their geographical location. In terms for 5 and 2 min of vigorous vortex. After adding 200 µL of of logistics, studies requiring large number of samples would chloroform, the samples were incubated for 3 min at room benefit if they could extract RNA from frozen material. In temperature. Then, the samples were centrifuged at 12,000 × g pigs, ejaculates from artificial insemination studs often contain for 15 min and the supernatants were transferred to new RNase- small amounts of epithelial cells but the presence of somatic free tubes, in which 500 µL of isopropanol was added. These and non-spermatozoa germline cells is rare and associated to samples were then centrifuged at 12,000 × g for 10 min. The inflammatory processes of the reproductive tract (Michael Kleve- supernatants were removed carefully and the pellets were washed Feld, Personal Communication). Our group has always purified with 500 µL of 75% (v/v) ethanol solution. The samples were the porcine ejaculates with BoviPureTM when characterizing the centrifuged at 13,000 × g for 5 min and dried out at room boar sperm transcriptome (5, 6, 8, 15) to be on the safe side since temperature for 10 min. The dried samples were resuspended in the potential impact of an even tiny presence of somatic cells on 30 µL of ultrapure water. All the centrifugations were performed the semen transcriptome profile is unknown. The objective of this at 4◦ C. All the RNA samples were subjected to DNase treatment pilot study was to assess the changes on the porcine ejaculate’s with the Turbo DNA-freeTM kit (Thermo Fisher Scientific; transcriptome after BoviPureTM purification and discuss the CA, USA) following the manufacturer’s instructions. The RNAs potential underlying causes. To achieve this goal, we compared were quantified with QubitTM RNA HS Assay kit (Invitrogen; the transcriptomes of purified vs. paired non-purified ejaculates CA, USA). from four male pigs. In order to have an initial evaluation of the purity of the extracted RNA defined as RNA originating exclusively from sperm cells and devoid of DNA, for each sample, we used MATERIALS AND METHODS three qPCR assays that assess the abundance of the sperm Sample Collection specific PRM1 RNA, the somatic-cell specific PTPRC RNA and Four fresh ejaculates showing good semen quality parameters the absence of genomic DNA (gDNA) using SYBR R Select (Supplementary Table 1) were obtained each from a different Master Mix (Thermo Fisher Scientific; CA, USA) and done in Pietrain boar from a commercial farm using the gloved hand triplicate as previously described (5). Total RNA abundance method (16), diluted (1:2) immediately into the prepared was quantified with Qubit RNA HS Assay kit (Invitrogen; CA, USA). Abbreviations: Cq, quantification cycle; DAG, differentially abundant gene; FPKM, fragments per kilobase of transcript per million mapped reads; GO, Gene RNA-Seq Library Preparation Ontology; NP, non-purified sample; P, purified sample; PCA, Principal component Sequencing libraries were prepared with the SMARTer Stranded analysis; TIN, transcript integrity number; tRNA, transfer RNA. Total RNA-Seq Kit v2—Pico Input Mammalian (Takara Bio) Frontiers in Veterinary Science | www.frontiersin.org 2 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq using Dual Indexing strategy, starting with 10 ng of RNA Evaluation of the Potential Tissue of Origin quantified with the Qubit NRA HS Assay Kit (Thermo Fisher of the Differentially Abundant Genes Scientific) and following the manufacturer’s instructions. The Between Purified and Non-purified final libraries were quality controlled on a Bioanalyzer (Agilent) using a High Sensitivity DNA Kit. The libraries were sequenced Samples in an Illumina’s HiSeq4000 system to generate 76 bp long We carried an analysis to determine the potential tissue or cell paired-end reads following the manufacturer’s protocol for dual type of origin of differentially abundant genes (DAG) between indexing. Image analysis, base calling and quality scoring of NP and P. We used a two-tailed Fisher’s Exact Test to compare the run were processed using the manufacturer’s software Real the proportion of genes with tissue-specific expression between Time Analysis (v 2.7.7) and followed by generation of FASTQ the DAG and the non-DAG groups. We used an RNA-Seq dataset sequence files. from an experiment carried by the Roslin Institute on 27 porcine tissues as part of the FAANG project (https://www.faang.org) with the NCBI’s BioProject Accession Number PRJEB19386. The RNA-Seq Mapping and Analysis gene expression values (FPKM) of this dataset are available at Total RNA-Seq reads were evaluated for quality control the European Bioinformatics Institute ArrayExpress expression by FastQC software (v, 0.11.5) (https://www.bioinformatics. atlas (https://www.ebi.ac.uk/gxa/home) archive under file E- babraham.ac.uk). Low-quality reads (Phred –Q < 20 and read MTAB-5895. The catalog included two male and two female length < 25 bp) and sequencing adaptors were trimmed with Duroc juveniles. All pigs were between four and six months Trimmomatic v.0.33 (17). Trimmed reads were aligned to the old. The file provided tissue expression results from males pig reference genome (Sscrofa 11.1) using HISAT2 v.2.1.0 (18) and females separately. This dataset includes epididymis, penis, with the default parameters. Duplicate reads were removed with tonsil, amygdala, mesenteric lymph node, lung alveolus, liver, PicardTools v11.0.5 MarkDuplicates (http://picard.sourceforge. heart, adipose tissue, omentum, several central nervous system net). RNA levels of the genes annotated in the porcine genome areas, skeletal muscle, spleen, several sections of the digestive (Ensembl v.101) were quantified as Fragments Per Kilobase of tract, cortex of kidney, pituitary gland, and uterus. Tonsil was transcript per Million mapped reads (FPKM), with StringTie only analyzed in males. Obviously, the reproductive tissues were v.1.3.4 (19). The genes with FPKM ≤ 1 were excluded from only represented by the one corresponding sex. We averaged further analysis. Differential gene abundance between P and the FPKM expression of the tissues analyzed in both males and NP was carried out by DESeq2 (20) and only these genes females. We interrogated the genes with specific expression in showing absolute log2 fold-change ≥ 1.5 and False Discovery epididymis as a male reproductive tissue and tonsil and lymph Rate (FDR) ≤ 0.05 between P and NP were considered to be node as a tissue representing leukocytes which can be also be differentially abundant. present in sperm. We also queried duodenum and skeletal muscle To determine whether the RNA-Seq quality metrics of the as two control tissues that are not expected to contribute cells P and NP samples were significantly different, we used R to to sperm. carry the Wilcoxon signed-rank test. The pairwise relationships We considered that a gene has a tissue-specific expression between NP and the paired P transcriptomes were assessed by when its RNA FPKM in the target tissue was >20 and linear regression and with the Pearson correlation coefficient the FPKM in all the other tissues in the dataset was <5. considering the purification treatment and the sample identifier. Duodenum, tonsils and skeletal muscle did not show any The hierarchical cluster dendrogram and the corresponding tissue-specific gene. Thus, for these tissues we allowed up to Selective Inference p-values were carried with the R package three additional tissues with FPKM > 5. The genes identified “pvclust” (21), and the Principal Component Analysis (PCA) following these conditions were considered as nearly tissue- and the correlation plot with the R package “ggplot2” (22). RNA specific. Then, we compared the abundance of these tissue- transcript integrity (TIN) was calculated with RseQC v.2.6.4 specific genes in the list of DAGs and the list of non-DAGs from (23) using the Ensembl v.101 pig annotation. TIN indicates the our study. proportion of a gene that is covered by reads ranging from We repeated this analysis using the GTEx Portal (dbGaP 0 (no coverage) to 100 (fully covered transcript). The TINs accession number phs000424.v8.ps), which contains RNA-Seq across paired samples was compared with the Wilcoxon signed- data from 54 human tissues including the prostate and testis rank test. male reproductive tissues as well as several representative of To assess whether the differences in gene abundance identified adipose tissue, artery, central nervous system, digestive tract, between P and NP were really caused by the purification step and female reproductive organs, whole blood, EBV transformed were not stochastic, we randomly shuffled the 8 samples into two lymphocytes, heart, kidney, skeletal muscle, skin, salivary gland, groups of four samples each, 10 times and compared each time, spleen, and thyroid. We used Ensembl BioMart (version 104) the transcriptome profiles of the two groups with DESeq2 (20). to identify the human orthologs of all the genes identified in Gene Ontology (GO) enrichment analysis was performed with our sperm (DAG and non-DAG) experiment. We queried the Cytoscape v.3.8.2 (24) plugin BiNGO v.3.0.4 using EBI porcine GTEx dataset using these human orthologous genes. We targeted Gene Ontology Annotation Database (release: 2021-02-01) with prostate and testis as male reproductive tissues, EBV-transformed default settings. Only the significant corrected p-values with FDR lymphocytes, and whole blood representing lymphocytes, which correction were considered. can be also found in sperm and liver, nerve tibial and mammary Frontiers in Veterinary Science | www.frontiersin.org 3 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq gland as control tissues which should not contain cell types that whole blood and PPARA and TRMU have ubiquitous expression. could be also found in sperm. On the contrary, FAM118A displays highest expression in testis, where spermatogenesis takes part, and prostate, which, RESULTS AND DISCUSSION at least in humans (25) and horse (26), contributes exosomes and can also provide epithelial cells to the human ejaculates RNA Extraction and qPCR Results (27). TTC38’s and FAM118A’s functions are not well known. RNA extractions showed very similar yields between the PPARA is a transcription factor involved in energy metabolism P and NP groups (Supplementary Table 2). Moreover, all and mitochondrial and peroxisomal function (28). TRMU is a the minus reverse transcription of the control samples mitochondrial tRNA modifying gene (29) thus related to protein showed no amplification of PRM1 and PTPRC and the synthesis. Sample B4 showed the largest disparity between NP qPCR on the intergenic region was undetectable in all P and and P of all the samples in all these analyses. NP samples, thereby indicating no DNA contamination Contrarily to what happens in somatic cells, most sperm RNAs (Supplementary Table 2). The eight samples presented are fragmented. Thus, we also evaluated RNA fragmentation quantification cycles (Cq) ranging between 15.7 and 17.6 by measuring TIN, which ranges between 0 and 100 indicating for the PRM1 sperm-specific gene (Supplementary Table 2). from full fragmentation to full integrity of the transcript, In contrast, the Cq for PTPRC blood-specific gene for the P and found similar (p = 0.37) TIN in the NP [mean ± samples ranged between 35.6 and 37.2, and between 32.6 and standard deviation (SD) = 19.44 ± 16.04] and the P (mean 37.6 in NP (Supplementary Table 2). The 1CqPTPRC−PRM1 , ± SD = 17.89 ± 15.45) samples thereby indicating no major calculated as the Cq for PTPRC minus the Cq for PRM1 differences between both groups (Supplementary Table 6). This ranged between 18.1 and 21.0 in P, and 17 to 20.4 in NP, which is probably due to the fact that the DAGs total RNAs represent means that, in average and as expected, the PTPRC assay had a small proportion of the whole sperm transcriptome abundance lower signal when compared to the PRM1 assay in P than in (Supplementary Table 4). NP (Supplementary Table 2). The NP aliquot of sample B2 showed more PTPRC signal (1CqPTPRC−PRM1 ) than any other of the eight samples thereby suggesting a degree of somatic cell Differential RNA Abundance Analysis presence in this sample. Between P and NP We also carried differential abundance analysis to identify Sequencing Metrics these genes with different RNA levels in P when compared We sequenced an average of 67.6 million reads per sample. to NP, a scenario that could be potentially caused by the Of these, an average of 7.6 million reads were uniquely removal of somatic cells, premature germline cells, immature, mapped to the porcine genome and used for downstream and morphologically abnormal spermatozoa or seminal analysis (Supplementary Table 3). There were no significant plasma exosomes from the ejaculates. We used the 7,519 differences between the sequencing metrics of P vs. NP samples genes with FPKM ≥ 1 for the differential analysis using (Supplementary Table 3). a paired sample design. Three hundred and seventy-two of the genes with FPKM ≥ 1 and 81 of the genes with Description and General Comparison of FPKM ≥ 10 showed significant differential RNA abundance the P and NP Transcriptomes between P and NP (Supplementary Table 4) and were thus Using HISAT2, we identified an average of 7,519 protein coding considered DAGs. To check whether the list of DAGs was genes with FPKM ≥ 1 and 1,870 genes with FPKM ≥ 10 stochastic or really owed to technical differences produced (Supplementary Table 4). The genes displaying most abundant by the purification step, we randomly shuffled 10 times RNA levels included PRM1, HSPB9, OAZ3, TSSK6, and TPPP2, the eight samples into two four-sample groups and carried among others and were mostly related to sperm biology differential gene abundance analysis between both groups. (Supplementary Table 4). These analyses showed between none and five DAGs in each Linear regression and correlation analysis based on the RNA comparison (Supplementary Table 7), thereby suggesting that levels of the 7,519 protein coding genes showed high correlations the purification step has an impact on the RNA abundance (between 0.93 in B4 and 0.99 in B3) between P and NP (Figure 1). of some genes. Most of the 372 DAGs showed a decreased Hierarchical clustering analysis of the 8 samples showed that RNA abundance in the P samples (Supplementary Table 4) the samples do not tend to group by the purification treatment which indicates that some cell types, organelles or vesicles were but by their sample of origin (Figure 2A). Nevertheless, PCA removed during the purification step. Only TNNI3, CTXN1, seemed to separate NP and P into two different groups SLC27A5, HSPA12B, and ENSSSCG00000032730 presented based on their principal component 1 (Figure 2B). Principal increased RNA levels in P in at least three of the four samples component 1 biggest contributors were TTC38 followed by (Supplementary Table 4). According to the GTEx portal (dbGaP PPARA, TRMU, and FAM118A, which together explained 57.2% accession number phs000424.v8.ps), TNNI3, a component of the of this component (Supplementary Table 5). We searched the striated muscle filaments, is exclusively expressed in the heart GTEx portal (dbGaP accession number phs000424.v8.ps) to although one study linked it to prostate cancer (30). CTXN1 is determine the expression of these genes in human tissues. most highly expressed in the brain but it is also present at lower TTC38 is mostly expressed in liver, small intestine, colon and levels in several other tissues including testis, and its function is Frontiers in Veterinary Science | www.frontiersin.org 4 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq FIGURE 1 | Linear regression and Pearson correlation of the transcriptome profiles between the purified and non-purified aliquots for the four ejaculates. The Pearson correlation and its p-value are annotated at the top left of each plot. related to brain biology (31). SLC27A5 is exclusively expressed clear presence in the sperm and an increase in P, we cannot in liver and is associated to lipid metabolism and bile synthesis exclude the possibility that this was a spurious result. The 372 (32). HSPA12B is highly abundant in spleen, lung, breast, and DAGs were enriched for GO terms related to translation and adipose tissues and has been related to heart injury (33) and lung transcription and metabolic processes (Supplementary Table 8) cancer (34). Thus, no particular link with the cells expected to and in fact, the list contained a large number of ribosomal be found in the ejaculates could be drew for these genes. TNNI3, proteins (e.g., RPS3A, RPS27A, and RPL7). The DAG group SLC27A5, HSPA12B, and ENSSSCG00000032730 presented was clearly enriched over the non-DAG list for RPL and RPS low average abundance across the eight samples (FPKM < 5) ribosomal protein genes (Two-Tailed Fisher’s Exact Test p < and even if significant, this data, should be considered with 1.9E-53), which showed decreased levels in P when compared caution as it could be spurious. CTXN1 displayed a much to NP. The majority of the DAGs showed an average abundance larger average abundance (FPKM = 85) and it is more likely below 10 FPKM (Supplementary Table 4). Among the list of to be real. Its increased levels in P indicate that a particular DAGs with average abundance above 10 FPKM, the top 10 that entity in the ejaculate—perhaps a particular sub-population of showed most significant differences included four ribosomal sperm cells—carrying this gene are more likely to be present proteins and six genes with distinct functions (Table 1). These after purification. However, as this is the only gene showing 10 genes are expressed in a wide catalog of tissues (GTEx portal; Frontiers in Veterinary Science | www.frontiersin.org 5 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq FIGURE 2 | Grouping of the sperm transcriptomes of the eight samples. (A) Cluster analysis and the Selective Inference p-values. (B) Principal Component Analysis biplot. dbGaP accession number phs000424.v8.ps), but preferentially component of the 40S ribosomal subunit (44). The interaction in EBV-transformed lymphocytes and most of them also in the between RACK1 and the AChE-R acetylcholinesterase splice female reproductive organs (Table 1). PTMA is related to cell variant in the spermatocytes of transgenic mice overexpressing proliferation (35), immunity (36), and chromatin remodeling AChE-R and displaying reduced sperm differentiation and (37) and it is also expressed in rat spermatocytes and spermatids sperm count showed that the interactions between AChE-R with but its functions in these cells are unknown (38). NANS is RACK1 and enolase-α may be related to these sperm related involved in the synthesis of sialic acid in multiple cell types phenotypes (45). and cellular processes (39) and within the catalog of human To get a glimpse of the tissue or cell type of origin of the tissues, it is most highly expressed in the prostate but it is DAGs group, we compared the proportion of the 372 DAGs also present in other tissues including sperm. Sialic acid is with preferential expression in the porcine epididymis, a male a key component of the sperm glycocalyx assembled during reproductive organ known to contribute RNAs to the ejaculate sperm development, maturation and upon contact with the (46), the tonsils and lymph node, that contain leukocytes, which seminal fluid and is related to sperm. TMSB4X is linked to actin could be also present in sperm. We also queried duodenum and polymerization and cytoskeleton organization (40) and is most skeletal muscle as control tissues that should not contain cell highly expressed in EBV transformed lymphocytes. EEF1A1 types also present in sperm. Three hundred and sixty-six DAG delivers tRNA to the ribosome for translation elongation (41). and 6,992 non-DAG genes were also present in the E-MTAB- One study carried RNA-Seq on chromatoid bodies isolated from 5895 porcine gene expression catalog. Fifty-seven genes displayed testicles from GRTH knock-in mice with a sterile phenotype specific expression in the epididymis, 5 and 13 of which were lacking elongated spermatids and spermatozoa (42). Chromatoid DAG and non-DAG, respectively (Table 2). Thus, 1.3% of the bodies are typically present in the cytoplasm of spermatocytes DAGs and 0.2% of the non-DAGs showed epididymis-specific and round spermatids but are absent in elongating spermatids expression. The two-tailed Fisher’s Exact Test showed a clear and spermatozoa and seem to be crucial for spermatogenesis. difference (p = 0.001) between the abundance of epididymis In this study, the authors identified an increased presence of specific genes in the DAG list when compared to the non-DAG genes related to transcript regulation, including EEF1A1 in the group. Skeletal muscle and duodenum showed 48 and 5 tissue- mutant when compared to the wild type mice. RACK1 acts as a specific genes, respectively. Five and one genes showed tissue- scaffold assisting and modulating protein-protein interactions specific expression in the tonsils and lymph node, respectively, to recruit, assemble or regulate signaling molecules (43) and is a but none of these was present in our porcine sperm samples. Frontiers in Veterinary Science | www.frontiersin.org 6 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq TABLE 1 | Ten most differentially abundant genes with FPKM ≥ 10 between the purified and paired non-purified samples. Ensembl identifier Gene name FDR Average B1 B2 B3 B4 Average Human tissues with highest FPKM P/NP P/NP P/NP P/NP P/NP expression according to the GTEx Portal ENSSSCG00000041806 – 5.74E-14 40 0.06 0.11 0.29 0.17 0.16 ENSSSCG00000037274 Prothymosin Alpha 2.31E-12 32 0.08 0.31 0.14 0.07 0.15 Quite ubiquotous but mostly (PTMA) expressed in EBV-transformed lymphocytes. Low expression in testis. ENSSSCG00000005373 N- 2.31E-12 16 0.03 0.23 0.30 0.10 0.16 Highest expression in the prostate Acetylneuraminate and colon. Synthase (NANS) ENSSSCG00000009019 Ribosomal Protein 2.92E-12 13 0.10 0.26 0.17 0.06 0.15 Mostly expressed in the female S3A (RPS3A) reproductive organs and EBV-transformed lymphocytes. ENSSSCG00000032111 Ribosomal Protein 1.86E-11 15 0.09 0.47 0.14 0.14 0.21 Mostly expressed in the female L7 (RPL7) reproductive organs and EBV-transformed lymphocytes. ENSSSCG00000012119 Thymosin Beta 4 3.66E-11 16 0.09 0.16 0.12 0.20 0.14 Mostly expressed in EBV-transformed X-Linked lymphocytes, spleen, lung, and whole (TMSB4X) blood. ENSSSCG00000004489 Eukaryotic 4.45E-11 93 0.07 0.48 0.29 0.10 0.24 Mostly expressed in cultured Translation fibroblasts, EBV-transformed Elongation Factor lymphocytes, ovary, and uterus. 1 Alpha 1 (EEF1A1) ENSSSCG00000033019 Ribosomal Protein 5.8E-11 16 0.05 0.20 0.17 0.10 0.13 Most highly expressed in L12 (RPL12) EBV-transformed lymphocytes, ovary, cultured fibroblasts, and cervix. ENSSSCG00000029724 Receptor For 6.6E-11 19 0.09 0.26 0.28 0.17 0.20 Most highly expressed in ovary, Activated C Kinase EBV-transformed lymphocytes, 1 (RACK1) cultured fibroblasts, and cervix. ENSSSCG00000031088 Ribosomal Protein 7.0E-11 81 0.10 0.58 0.33 0.28 0.32 Most highly expressed in ovary, S7 (RPS7) EBV-transformed lymphocytes and cultured fibroblasts. B1 P/NP, Ratio of the FPKM in the purified and the FPKM in the non-purified in B1; B2 P/NP, Ratio of the FPKM in purified and the FPKM in non-purified in B2; B3 P/NP, Ratio of the FPKM in purified and the FPKM in non-purified in B3; B4 P/NP, Ratio of the FPKM in purified and the FPKM in non-purified in B4; Average P/NP, Average of the ratio of the FPKM in purified and the FPKM in non-purified in the four samples. For these three tissues we allowed up to three additional tissues The genes HBB, CDH1, PTPRC, KRT1, KRT10, CXCL8, with FPKM > 5. None of these four tissues showed a statistical and KLK3 previously used to determine the presence of difference between DAGs and non-DAGs (Table 2). leukocytes, epithelial cells, and prostate in human (47) and pig We carried the same analysis on the dataset of gene expression (5) showed no detectable levels of RNA in P and NP. Only in 54 human tissues from the GTEx Portal (dbGaP accession HBB presented moderate abundance (FPKM = 36) in the NP number phs000424.v8.ps) querying testis and prostate as male sample from B4 (Supplementary Table 4). This data indicates reproductive organs, EBV transformed lymphocytes and whole that most NP samples are free of detectable levels of RNAs blood as representatives of leukocytes and liver, nerve tibial and from somatic cell origin but that B4 may have contained a mammary gland as negative controls. Three hundred and twenty- sufficient proportion of leukocytes to provide detectable levels seven and 6,304 DAGs and non-DAGs had a human ortholog of HBB. This is in line with the results on the cluster and present in the GTEx catalog (Table 2). Only testis, liver, and nerve PCA analyses which showed that the NP sample of B4 clustered tibial showed tissue-specific genes (Table 2). For prostate, EBV and mapped apart from all the other samples and the fact transformed lymphocytes, whole blood and mammary gland, we that this sample showed the highest abundance for 266 of the selected those genes which showed FPKM > 20 in the target 372 DAGs. tissue and FPKM > 5 in not more than three additional tissues. In light of these results, we wanted to explore whether The list of testis-specific genes were the only ones that showed a the contribution of each sample to the list of DAGs was statistically significant difference between DAGs and non-DAGs homogeneous or on the contrary, there was high variability (p = 5.4E-07) but these genes were enriched in the non-DAG between samples. We carried four differential abundance group (Table 2). analyses between P and the paired NP samples but this Frontiers in Veterinary Science | www.frontiersin.org 7 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq TABLE 2 | Number of tissue-specific or nearly tissue-specific genes in the DAG and non-DAG groups and p-value of the comparison of the proportion of these genes in both groups. Tissue Number of Number of Number of Number of Number in Number in P-value tissue-specific tissue-specific nearly nearly DAG non-DAG genes genes present in tissue-specific tissue-specific pig sperm genes genes present in pig sperm Epididymis 57 18 N/A N/A 5 13 0.001 Skeletal muscle 48 6 N/A N/A 0 6 1 Duodenum 5 2 N/A N/A 0 2 1 Lymph node 1 0 14 1 0 1 1 Tonsils 5 0 50 10 2 8 0.08 Testis 908 295 N/A N/A 0 295 5.40E-07 Prostate 5 0 16 3 1 2 0.14 EBV transformed 13 0 149 40 0 40 0.26 lymphocytes Whole blood 7 0 105 6 0 6 1 Mammary 0 0 7 0 0 0 1 Liver 85 3 N/A N/A 0 3 1 Nerve tibial 2 1 N/A N/A 0 1 1 Number in DAG, number of tissue-specific or nearly tissue-specific genes in the DAG group. Number in non-DAG, number of tissue-specific or nearly tissue-specific genes in the non-DAG group. P-value, P-value of the two-tailed Fisher’s Exact Test comparing the proportion of tissue specific or nearly tissue-specific genes in the DAG vs. the non-DAG groups. N/A, non applicable. time, removing each time, a different sample and thus forming the ribosomal sub-units. If true, this would indicate that carrying the comparison on only three ejaculates. The results the DAGs are contributed mostly by non-sperm entities. demonstrated that the removal of sample B4 caused the most This is also supported by the enriched abundance of human dramatic reduction of the number of DAGs when compared testis specific genes in the non-DAG list. We hypothesize to the original differential abundance analysis carried on the that these non-DAG genes are specific to the later stages of four samples (Supplementary Table 9). Despite this outcome, spermatogenesis and they are present in the mature spermatozoa, the correlation between the NP and P transcriptomes of which are not removed by the BoviPureTM treatment. As a matter B4 was still high (R2 = 0.93) which shows that even in of fact, the list of 295 sperm-specific genes present in the non- this sample, the BoviPureTM purification did not have a DAG include important genes related to sperm biology and dramatic impact of the overall transcriptome profile of this function such as PRM1, PRM2, CATSPER1, CATSPERD, DAZL, sample. Remarkably, the removal of B2 did not show any and several members of the SPATA and SPACA families among relevant change in the number of DAGs which, contrarily others (Supplementary Table 10). to what was indicated by the 1CqPTPRC−PRM1 , suggests As reviewed by Fedder in 1996, research in human semen no particular contamination of somatic cell RNAs. The has identified the presence of leukocytes, squamous epithelial semen quality phenotypes of the fresh ejaculate of B4 cells, Sertoli cells, spermatogonia, spermatocytes, and spermatids before purification showed lower ejaculate volume but all (49). In addition, the presence of epithelial cells from the the other parameters (sperm cell concentration, sperm cell epididymis (50) and the prostate (27) have been also confirmed. viability, percentage of spermatozoa with head, tail, or neck In swine, sperm from artificial insemination studs rarely contain abnormalities, percentage of sperm cells with proximal or leukocytes but low levels of epithelial cells are frequently distal droplets and percentage of motile spermatozoa), did not found (Michael Kleve-Feld, Personal Communication). We deviate from the values observed in the other three samples detected 57 pig epididymis-specific genes, 18 of which were (Supplementary Table 1). Therefore, we could not link the present in our boar sperm samples. The DAG group was difference on the transcriptome of this sample to the semen enriched for such genes, which indicates the removal of quality phenotypes. epididysomal cells or exosomes after the purification. On the A somehow surprising finding was the clear enrichment of other side, we did not find statistically significant differences RPL and RPS ribosomal protein genes in the DAG group over in the proportion of human prostate-specific genes between the non-DAG list (Two-Tailed Fisher’s Exact Test p < 1.9E- the DAG and the non-DAG groups. This could be indicating 53), which showed decreased abundance in P when compared that the prostate did not contribute a detectable level of to NP. As intact ribosomal RNAs (rRNAs) are depleted in epithelial cells or exosomes to the pig sperm but it could also sperm to warrant translational silencing (48), we hypothesize be caused by the small number of human prostate’s nearly that spermatozoa may also have reduced RNA levels of proteins specific genes, only 16, that we identified in the study, with Frontiers in Veterinary Science | www.frontiersin.org 8 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq resulted in only one and two genes present in the DAG and non-purified samples. In any case, the results identified in and the non-DAG lists, respectively. Notwithstanding, the swine cannot be translated to humans because the pig semen absence of detectable levels of prostate cells in sperm is also contains a much lower proportion of somatic cells than the supported by the null presence of KLK3 in all the NP and men’s ejaculate. P samples. The presence of epididymal material in sperm could be due to either epithelial cells or exosomes or a combination of both. Seminal exosomes have been found in the boar sperm (51). RNA- CONCLUSION Seq data from human (52) and porcine (46) seminal exosomes Three of the four ejaculates had undetectable levels of somatic showed that the RNA cargo of these vesicles is mostly made cell RNA markers before purification, thereby showing an up of short RNAs but also include protein-coding RNAs. In agreement with the low levels of non-sperm cells typically pigs, RNA-Seq of the small RNA fraction of seminal exosomes found in samples from boar artificial insemination studs. We from ejaculates displaying normal semen quality parameters found strong indications that the purification with BoviPureTM showed that 25.3% of the sequencing reads corresponded to has a mild but noticeable impact on the RNA abundance messenger RNAs (mRNAs). As the objective of this study was of some genes and that this was originated by the removal on microRNAs, the authors did not provide additional detail of non-sperm cell entities with the most likely contribution on the messenger RNAs present in the boar exosomes (46). of epididymal epithelial cells. The evaluation of the 10 most Whether the BoviPureTM purification removes the seminal differentially abundant genes also indicated the removal of vesicles is unknown. From the one side, BoviPureTM separates premature germline cells—even immature or morphologically the seminal fluid from the mature spermatozoa which should abnormal spermatozoa—as well as leukocytes with a distinct eliminate the exosomes present in the fluid. But on the other RNA cargo. In light of these results, our group feels comfortable side, seminal exosomes also bind to the sperm’s cell membrane with the use of non-purified samples for the interrogation and their removal after gradient-based protocols is uncertain. of the relationship between spermatozoon RNAs and semen We don’t know whether these exosomes unbind from the quality traits. However, this should be decided by each research sperm during the purification or whether the sperm cells bound group in accordance to the experimental design and objective of to exosomes have a distinct density that could remove them their study. after purification. This also happens for the NP samples as they were subjected to a 14,000 g centrifugation for 4 minutes. Exosomes sedimentation often uses ultra-centrifugations of at least 120,000 g for 70 min. However, whether the sperm bound DATA AVAILABILITY STATEMENT exosomes will sediment alongside with the sperm is unsure. In consequence, we hypothesize that the most likely largest The data presented in the study are deposited in the NCBI’s contributor to the DAG list are epididysomal epithelial cells. A short read archive repository, accession numbers SRR14117413, somatic cell has been estimated to contain in average between SRR14117412, SRR14117411, SRR14117410, SRR14117409, 10 and 30 picograms of RNA, in contrast with an average SRR14117408, SRR14117407, and SRR14117406. of 15 femtograms for a human spermatozoon. In addition, our own estimation suggests that the porcine sperm contains in average 1.6 femtograms per cell (5). Hence, the presence of an even tiny proportion of these epithelial cells in the ETHICS STATEMENT sperm might be able to contribute detectable levels of RNA in Ethical review and approval was not required for the animal study the sample. because this study was done on ejaculates that were privately In order to further identify the origin of the RNAs impacted owned by a commercial farm for non-research purposes. The by BoviPureTM , future studies should be carried. To determine owners provided consent for their use for research. Specialized the contribution of abnormal or immature spermatozoa, the professionals at the farm collected the ejaculates following comparison should include an additional control consisting of standard routine procedures and guidelines. an aliquot treated with somatic cell lysis buffer or sonication. The comparison between the three conditions should be still taken with caution as the SDS present in the lysis buffer could solubilise the sperm’s membrane and remove its midpiece, AUTHOR CONTRIBUTIONS which contains a large number of mitochondria thereby altering the resulting transcriptome (13). The contribution AS and ACl conceived and designed the experiments. SB of the seminal exosomes could be further elucidated by collected the samples. JR-G carried the phenotypic analysis. ACa purifying these vesicles by either ultra-centrifugation, or filtration carried the qPCR analyses. MG performed sperm purifications by gradient based protocols and comparing the resulting and RNA extractions. YL and MG made the bioinformatics and transcriptome with that of the NP and P samples. This statistics analysis. YL, MG, and ACl analyzed the data. ACl, YL, experiment should be accompanied by the exploration under JR-G, and MG wrote the manuscript. All authors discussed the the electron microscope of extracellular vesicles in the purified data and read and approved the contents of the manuscript. Frontiers in Veterinary Science | www.frontiersin.org 9 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq FUNDING Non-purified sample of B2; P_B3: Purified sample of B3; NP_B3: Non-purified sample of B3; P_B4: Purified sample of B4; NP_B4: Non-purified sample of B4. This work was supported by the Spanish Ministry of Economy Supplementary Table 4 | Gene abundance across the eight samples for the and Competitiveness (MINECO) under grant AGL2013- 7,519 protein coding genes with average FPKM ≥ 1. DAG, Differentially abundant 44978-R and grant AGL2017-86946-R and by the CERCA gene; FDR, False Discovery Rate; P_B1, Purified sample of B1; NP_B1, Programme/Generalitat de Catalunya. AGL2017-86946-R was Non-purified sample of B1; P_B2, Purified sample of B2; NP_B2, Non-purified sample of B2; P_B3, Purified sample of B3; NP_B3, Non-purified sample of B3; also funded by the Spanish State Research Agency (AEI) and P_B4, Purified sample of B4; NP_B4, Non-purified sample of B4. the European Regional Development Fund (ERDF). We thank Supplementary Table 5 | Major gene contributors to principal component 1 of the Agency for Management of University and Research Grants the Principal Component Analysis. (AGAUR) of the Generalitat de Catalunya (Grant Numbers 2014 SGR 1528 and 2017 SGR 1060). We also acknowledge financial Supplementary Table 6 | Transcript Integrity Number values across the eight support from the Spanish Ministry of Science and Innovation, samples. P_B1, Purified sample of B1; NP_B1, Non-purified sample of B1; P_B2, Purified sample of B2; NP_B2, Non-purified sample of B2; P_B3, Purified sample through the Severo Ochoa Programme for Centers of Excellence of B3; NP_B3, Non-purified sample of B3; P_B4, Purified sample of B4; NP_B4, in R&D SEV-2015-0533 and CEX2019-000902-S. YL thanks Non-purified sample of B4. the China Scholarship Council for her PhD scholarship. We acknowledge support of the publication fee by the CSIC Open Supplementary Table 7 | Results of the 10 shuffles to organize the eight samples into 2 four-sample random groups. P_B1, Purified sample of B1; NP_B1, Access Publication Support Initiative through its Unit of Non-purified sample of B1; P_B2, Purified sample of B2; NP_B2, Non-purified Information Resources for Research (URICI). sample of B2; P_B3, Purified sample of B3; NP_B3, Non-purified sample of B3; P_B4, Purified sample of B4; NP_B4, Non-purified sample of B4. Random shuffles are identified as Random1, Random2, Random3, Random4, Random5, SUPPLEMENTARY MATERIAL Random6, Random7, Random8, Random9, Random10. A and B denote the randomly generated groups in each shuffle. The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets. Supplementary Table 8 | Results of the Gene Ontology analysis of the 372 2021.668158/full#supplementary-material differentially abundant genes between the purified and the non-purified samples. Supplementary Table 1 | Phenotypic values for several semen quality Supplementary Table 9 | List of differentially abundant genes when including the parameters of the four ejaculates used in this study. four ejaculates and when removing the ejaculates B1, B2, B3, or B4. The second Supplementary Table 2 | RNA extraction and qPCR data. P_B1, Purified sample row shows the number of DAGs identified in that analyses. Down, gene showing of B1; NP_B1, Non-purified sample of B1; P_B2, Purified sample of B2; NP_B2, significantly decreased RNAs levels in P when compared to NP. Up, gene showing Non-purified sample of B2; P_B3, Purified sample of B3; NP_B3, Non-purified significantly increased RNAs levels in P when compared to NP. sample of B3; P_B4, Purified sample of B4; NP_B4, Non-purified sample of B4; Supplementary Table 10 | List of tissue-specific or nearly tissue-specific genes P_average, Average value of the purified samples; NP_average, Average value of for pig epididymis, lymph node, tonsil, duodenum and skeletal muscle and for the non-purified samples. human prostate, testis, EBV transformed lymphocytes, whole blood, liver, and Supplementary Table 3 | Sequencing metrics. P_B1: Purified sample of B1; nerve tibial that are also present in the pig sperm. The last column shows whether NP_B1: Non-purified sample of B1; P_B2: Purified sample of B2; NP_B2: that particular gene belonged to the DAG or the non-DAG group. REFERENCES 8. Ablondi M, Gòdia M, Rodriguez-Gil JE, Sánchez A, Clop A. Characterisation of sperm piRNAs and their correlation with semen quality traits in swine. 1. Gòdia M, Swanson G, Krawetz SA. A history of why fathers’ RNA matters. Biol Anim Genet. (2021) 52:114–20. doi: 10.1111/age.13022 Reprod. (2018) 99:147–59. doi: 10.1093/biolre/ioy007 9. Fraser L, Brym P, Pareek CS, Mogielnicka-Brzozowska M, Paukszto, 2. Jodar M, Sendler E, Moskovtsev SI, Librach CL, Goodrich R, Swanson S, et al. Jastrzebski JP, et al. Transcriptome analysis of boar spermatozoa with Absence of sperm RNA elements correlates with idiopathic male infertility. Sci different freezability using RNA-Seq. Theriogenology. (2020) 142:400– Transl Med. (2015) 7:295re6. doi: 10.1126/scitranslmed.aab1287 13. doi: 10.1016/j.theriogenology.2019.11.001 3. Keles E, Malama E, Bozukova S, Siuda M, Wyck S, Witschi U, et al. 10. Alvarez-rodriguez M, Martinez C, Wright D, Barranco I, Roca J, Rodriguez- The micro-RNA content of unsorted cryopreserved bovine sperm and its martinez H. The transcriptome of pig spermatozoa, and its role in fertility. Int relation to the fertility of sperm after sex-sorting. BMC Genomics. (2021) J Mol Sci. (2020) 21:1572. doi: 10.3390/ijms21051572 22:30. doi: 10.1186/s12864-020-07280-9 11. Das PJ, McCarthy F, Vishnoi M, Paria N, Gresham C, Li G, et al. Stallion 4. Suliman Y, Becker F, Wimmers K. Implication of transcriptome profiling sperm transcriptome comprises functionally coherent coding and regulatory of spermatozoa for stallion fertility. Reprod Fertil Dev. (2018) 30:1087– RNAs as revealed by microarray analysis and RNA-seq. PLoS One. (2013) 98. doi: 10.1071/RD17188 8:e56535. doi: 10.1371/journal.pone.0056535 5. Gòdia M, Mayer FQ, Nafissi J, Castelló A, Rodríguez-Gil JE, Sánchez A, 12. Gòdia M, Reverter A, González-Prendes R, Ramayo-Caldas Y, Castelló A, et al. A technical assessment of the porcine ejaculated spermatozoa for Rodríguez-Gil JE, et al. A systems biology framework integrating GWAS and a sperm-specific RNA-seq analysis. Syst Biol Reprod Med. (2018) 64:291– RNA-seq to shed light on the molecular basis of sperm quality in swine. Genet 303. doi: 10.1080/19396368.2018.1464610 Sel Evol. (2020) 52:1–21. doi: 10.1186/s12711-020-00592-0 6. Gòdia M, Estill M, Castelló A, Balasch S, Rodríguez-Gil JE, Krawetz SA, et al. 13. Mao S, Goodrich RJ, Hauser R, Schrader SM, Chen Z, Krawetz SA. Evaluation A RNA-seq analysis to describe the boar sperm transcriptome and its seasonal of the effectiveness of semen storage and sperm purification methods for changes. Front Genet. (2019) 10:299. doi: 10.3389/fgene.2019.00299 spermatozoa transcript profiling. Syst Biol Reprod Med. (2013) 59:287– 7. Gòdia M, Castelló A, Rocco M, Cabrera B, Rodríguez-Gil JE, Balasch 95. doi: 10.3109/19396368.2013.817626 S, et al. Identification of circular RNAs in porcine sperm and 14. Nätt D, Kugelberg U, Casas E, Nedstrand E, Zalavary S, Henriksson P, evaluation of their relation to sperm motility. Sci Rep. (2020) 10:7985. et al. Human sperm displays rapid responses to diet. PLoS Biol. (2019) doi: 10.1038/s41598-020-64711-z 17:e3000559. doi: 10.1371/journal.pbio.3000559 Frontiers in Veterinary Science | www.frontiersin.org 10 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq 15. Godia M, Castello A, Rocco M, Cabrera B, Rodriguez-Gil J, Sanchez A, et al. 34. Ma H, Lu T, Zhang X, Li C, Xiong J, Huang L, et Identification of circular RNAs in porcine sperm and their relation to sperm al. HSPA12B: a novel facilitator of lung tumor growth. motility. Sci. Rep. (2019). doi: 10.1101/608026 Oncotarget. (2015) 6:9924–36. doi: 10.18632/oncotarget. 16. King GJ, W MJ. A comparison of two methods for boar semen collection. J 3533 Anim Sci. (1973) 36:563–5. doi: 10.2527/jas1973.363563x 35. Letsas KP, Frangou-Lazaridis M. Surfing on prothymosin alpha proliferation 17. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible and anti-apoptotic properties. Neoplasma. (2006) 53:92–6. trimmer for Illumina sequence data. Bioinformatics. (2014) 36. Ioannou K, Derhovanessian E, Tsakiri E, Samara P, Kalbacher H, 30:2114–20. doi: 10.1093/bioinformatics/btu170 Voelter W, et al. Prothymosin α and a prothymosin α-derived 18. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced peptide enhance TH1-type immune responses against defined HER- aligner with low memory requirements. Nat Methods. (2015) 2/neu epitopes. BMC Immunol. (2013) 14:43. doi: 10.1186/1471-2172 12:357–60. doi: 10.1038/nmeth.3317 -14-43 19. Pertea M, Pertea GM, Antonescu CM, Chang T, Mendell JT, Salzberg SL. 37. Gómez-Márquez J. Function of prothymosin α in chromatin decondensation StringTie enables improved reconstruction of a transcriptome from RNA-seq and expression of thymosin β-4 linked to angiogenesis and synaptic reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122 plasticity. Ann N Y Acad Sci. (2007) 1112:201–9. doi: 10.1196/annals.14 20. Love MI, Huber W, Anders S. Moderated estimation of fold change 15.020 and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 38. Ferrara D, Izzo G, Pariante P, Donizetti A, D’Istria M, 15:550. doi: 10.1186/s13059-014-0550-8 Aniello F, et al. Expression of prothymosin alpha in meiotic 21. Suzuki R, Shimodaira H. Pvclust: An R package for assessing the and post-meiotic germ cells during the first wave of rat uncertainty in hierarchical clustering. Bioinformatics. (2006) 22:1540– spermatogenesis. J Cell Physiol. (2010) 224:362–8. doi: 10.1002/jcp. 2. doi: 10.1093/bioinformatics/btl117 22131 22. Wickham H. ggplot2: Elegant Graphics for Data Analysis. (2016) Available 39. Schnaar RL, Gerardy-Schahn R, Hildebrandt H. Sialic acids online at: https://ggplot2.tidyverse.org (accessed April 23, 2021). in the brain: gangliosides and polysialic acid in nervous 23. Wang L, Nie J, Sicotte H, Li Y, Eckel-Passow JE, Dasari S, et al. system development, stability, disease, and regeneration. Measure transcript integrity using RNA-seq data. BMC Bioinformatics. (2016) Physiol Rev. (2014) 94:461–518. doi: 10.1152/physrev.0003 17:58. doi: 10.1186/s12859-016-0922-z 3.2013 24. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. 40. Cassimeris L, Safer D, Nachmias VT, Zigmond SH. Thymosin beta 4 Cytoscape: a software environment for integrated models. Genome Res. sequesters the majority of G-actin in resting human polymorphonuclear (1971) 13:426. leukocytes. J Cell Biol. (1992) 119:1261–70. doi: 10.1083/jcb.119.5 25. Aalberts M, van Dissel-Emiliani FMF, van Adrichem NPH, van .1261 Wijnen M, Wauben MHM, Stout TAE, et al. Identification of distinct 41. Knight JRP, Gavin Garland G, Pöyry T, Mead E, Vlahov N, populations of prostasomes that differentially express prostate stem Sfakianos A, et al. Control of translation elongation in health and cell antigen, annexin A1, and GLIPR2 in humans. Biol Reprod. (2012) disease. Dis Model Mech. (2020) 13:dmm043208. doi: 10.1242/dmm. 86:1–8. doi: 10.1095/biolreprod.111.095760 043208 26. Aalberts M, Sostaric E, Wubbolts R, Wauben MWM, Nolte-’T 42. Anbazhagan R, Kavarthapu R, Coon SL, Dufau ML. Role Hoen ENM, Gadella BM, et al. Spermatozoa recruit prostasomes of phosphorylated gonadotropin-regulated testicular RNA in response to capacitation induction. Biochim Biophys Acta helicase (GRTH/DDX25) in the regulation of germ cell specific Proteins Proteomics. (2013) 1834:2326–35. doi: 10.1016/j.bbapap.2012. mRNAs in chromatoid bodies during spermatogenesis. Front 08.008 Cell Dev Biol. (2020) 8:580019. doi: 10.3389/fcell.2020.58 27. Andrade-Rocha FT. Assessment of exfoliated prostate cells in 0019 semen: Relationship with the secretory function of the prostate. 43. Schechtman D, Mochly-Rosen D. Adaptor proteins in protein Am J Clin Pathol. (2007) 128:788–93. doi: 10.1309/9ELYPH61DBD kinase C-mediated signal transduction. Oncogene. (2001) 20:6339– 8AXYC 47. doi: 10.1038/sj.onc.1204778 28. Vamecq J, Latruffe N. Medical significance of peroxisome proliferator- 44. Rabl J, Leibundgut M, Aitade SF, Haag A, Ban N. Crystal structure of the activated receptors. Lancet. (1999) 354:141–8. doi: 10.1016/S0140-6736(98)10 eukaryotic 60S ribosomal subunit in complex with initiation factor 1. Science. 364-1 (2011) 331:730–6. doi: 10.1126/science.1198308 29. Sasarman F, Antonicka H, Horvath R, Shoubridge EA. The 45. Mor I, Sklan EH, Podoly E, Pick M, Kirschner M, Yogev 2-thiouridylase function of the human MTU1 (TRMU) L, et al. Acetylcholinesterase-R increases germ cell apoptosis enzyme is dispensable for mitochondrial translation. but enhances sperm motility. J Cell Mol Med. (2008) 12:479– Hum Mol Genet. (2011) 20:4634–43. doi: 10.1093/hmg/ 95. doi: 10.1111/j.1582-4934.2008.00231.x ddr397 46. Xu Z, Xie Y, Zhou C, Hu Q, Gu T, Yang J, et al. Expression pattern of seminal 30. Gupta A, Shukla N, Nehra M, Gupta S, Malik B, Mishra plasma extracellular vesicle small RNAs in boar semen. Front Vet Sci. (2020) AK, et al. A pilot study on the whole exome sequencing 7:585276. doi: 10.3389/fvets.2020.585276 of prostate cancer in the indian phenotype reveals distinct 47. Jodar M, Sendler E, Moskovtsev SI, Librach CL, Goodrich R, Swanson polymorphisms. Front Genet. (2020) 11:874. doi: 10.3389/fgene.2020. S, et al. Response to Comment on “Absence of sperm RNA elements 00874 correlates with idiopathic male infertility.” Sci Transl Med. (2016) 31. Levchuk LA, Ivanova SA, Semke VY. Effects of neuroprotector 8:353tr1. doi: 10.1126/scitranslmed.aaf4550 cortexin on the dynamics of neuroendocrine system parameters 48. Kierszenbaum AL, Tres LL, Hill C, Carolina N. Structural and transcriptional in patients with organic emotionally labile (asthenic) disorders. features of the mouse spermatid genome. J Cell Biol. (1975) 65:258– Bull Exp Biol Med. (2013) 155:75–7. doi: 10.1007/s10517-013-2 70. doi: 10.1083/jcb.65.2.258 083-7 49. Fedder J. Nonsperm cells in human semen: with special reference to seminal 32. Sharma V, Hiller M. Loss of enzymes in the bile acid synthesis leukocytes and their possible influence on fertility. Arch Androl. (1996) pathway explains differences in bile composition among 36:41–65. doi: 10.3109/01485019608987883 mammals. Genome Biol Evol. (2018) 10:3211–7. doi: 10.1093/gbe/ 50. Weiser D, Mietens A, Stadler B, Ježek D, Schuler G, Middendorff R. evy243 Contractions transport exfoliated epithelial cells through the neonatal 33. Zhao Y, Liu C, Liu J, Kong Q, Mao Y, Cheng H, et al. HSPA12B promotes epididymis. Reproduction. (2020) 160:109–16. doi: 10.1530/REP-19-0617 functional recovery after ischaemic stroke through an eNOS-dependent 51. Du J, Shen J, Wang Y, Pan C, Pang W, Diao H, et al. Boar seminal plasma mechanism. J Cell Mol Med. (2018) 22:2252–62. doi: 10.1111/jcmm. exosomes maintain sperm function by infiltrating into the sperm membrane. 13507 Oncotarget. (2016) 7:58832–47. doi: 10.18632/oncotarget.11315 Frontiers in Veterinary Science | www.frontiersin.org 11 July 2021 | Volume 8 | Article 668158 Lian et al. Pig’s Sperm Purification and RNA-Seq 52. Vojtech L, Woo S, Hughes S, Levy C, Ballweber L, Sauteraud The remaining authors declare that the research was conducted in the absence of RP, et al. Exosomes in human semen carry a distinctive any commercial or financial relationships that could be construed as a potential repertoire of small non-coding RNAs with potential regulatory conflict of interest. functions. Nucleic Acids Res. (2014) 42:7290–304. doi: 10.1093/nar/ gku347 Copyright © 2021 Lian, Gòdia, Castello, Rodriguez-Gil, Balasch, Sanchez and Clop. This is an open-access article distributed under the terms of the Creative Commons Conflict of Interest: SB was employed by the company Group Gepork S.A. Attribution License (CC BY). The use, distribution or reproduction in other forums The authors declare that this study received funding from MINECO, AGAUR, is permitted, provided the original author(s) and the copyright owner(s) are credited Generalitat de Catalunya and MICINN. The funder was not involved in the study and that the original publication in this journal is cited, in accordance with accepted design, collection, analysis, interpretation of data, the writing of this article or the academic practice. No use, distribution or reproduction is permitted which does not decision to submit it for publication. comply with these terms. Frontiers in Veterinary Science | www.frontiersin.org 12 July 2021 | Volume 8 | Article 668158