Gòdia et al. Genet Sel Evol (2020) 52:72 https://doi.org/10.1186/s12711-020-00592-0 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access A systems biology framework integrating GWAS and RNA‑seq to shed light on the molecular basis of sperm quality in swine Marta Gòdia1 , Antonio Reverter2 , Rayner González‑Prendes3 , Yuliaxis Ramayo‑Caldas4 , Anna Castelló1,5 , Joan‑Enric Rodríguez‑Gil6 , Armand Sánchez5 and Alex Clop1,7* Abstract Background: Genetic pressure in animal breeding is sparking the interest of breeders for selecting elite boars with higher sperm quality to optimize ejaculate doses and fertility rates. However, the molecular basis of sperm quality is not yet fully understood. Our aim was to identify candidate genes, pathways and DNA variants associated to sperm quality in swine by analysing 25 sperm-related phenotypes and integrating genome-wide association studies (GWAS) and RNA-seq under a systems biology framework. Results: By GWAS, we identified 12 quantitative trait loci (QTL) associated to the percentage of head and neck abnormalities, abnormal acrosomes and motile spermatozoa. Candidate genes included CHD2, KATNAL2, SLC14A2 and ABCA1. By RNA-seq, we identified a wide repertoire of mRNAs (e.g. PRM1, OAZ3, DNAJB8, TPPP2 and TNP1) and miRNAs (e.g. ssc-miR-30d, ssc-miR-34c, ssc-miR-30c-5p, ssc-miR-191, members of the let-7 family and ssc-miR-425-5p) with functions related to sperm biology. We detected 6128 significant correlations (P-value ≤ 0.05) between sperm traits and mRNA abundances. By expression (e)GWAS, we identified three trans-expression QTL involving the genes IQCJ, ACTR2 and HARS. Using the GWAS and RNA-seq data, we built a gene interaction network. We considered that the genes and interactions that were present in both the GWAS and RNA-seq networks had a higher probability of being actually involved in sperm quality and used them to build a robust gene interaction network. In addition, in the final network we included genes with RNA abundances correlated with more than four semen traits and miRNAs interact‑ ing with the genes on the network. The final network was enriched for genes involved in gamete generation and development, meiotic cell cycle, DNA repair or embryo implantation. Finally, we designed a panel of 73 SNPs based on the GWAS, eGWAS and final network data, that explains between 5% (for sperm cell concentration) and 36% (for percentage of neck abnormalities) of the phenotypic variance of the sperm traits. Conclusions: By applying a systems biology approach, we identified genes that potentially affect sperm quality and constructed a SNP panel that explains a substantial part of the phenotypic variance for semen quality in our study and that should be tested in other swine populations to evaluate its relevance for the pig breeding sector. Background Sperm carries the paternal genome and a wide reper- toire of molecules including RNAs, which are essential *Correspondence: alex.clop@cragenomica.es for fertilization and the development of a new organism. 1 Animal Genomics Group, Centre for Research in Agricultural Genomics Spermatogenesis, the process whereby germ cells pro- (CRAG) CSIC-IRTA-UAB-UB, Campus UAB, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain liferate and develop into mature spermatozoa, is con- Full list of author information is available at the end of the article trolled by multiple factors. Both DNA polymorphisms © The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat​iveco​ mmons​.org/licen​ses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creat​iveco​mmons​.org/publi​cdoma​in/ zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Gòdia et al. Genet Sel Evol (2020) 52:72 Page 2 of 21 and gene expression have been linked to sperm quality identified five and two miRNAs associated to sperm mor- and/or fertility in several mammalian species including phology and motility, respectively. Moreover, our group cattle [1] and swine [2] and review in [3]. High-quality has also identified a correlation between the abundance sperm is decisive to optimize the propagation of the best of some circRNAs [19] and piRNAs [20] with semen genetic material in livestock and the sustainability of the quality parameters in swine. pig breeding sector. For this reason, ejaculated sperm is Based on these recent studies, it is now clear that the subjected to strict quality filters in boar artificial insemi- genetic complexity of sperm quality involves several nation (AI) studs. AI farms regularly evaluate the quality molecular mechanisms and pathways that are highly of ejaculates by measuring traits such as concentration, interconnected. Complex traits are typically affected by morphology, viability and motility kinetics, as a way to a large number of genomic regions, many of which may predict their fertilizing ability [4]. Although the heritabil- explain only a small proportion of the phenotypic vari- ity of these traits is low to moderate [5–8], the molecular ance and do not reach significant levels in a GWAS or dif- processes and genetic mechanisms that control sperm ferential expression analysis. Moreover, classical GWAS quality are not yet fully understood and boar replacement or differential expression analyses carried out on a single due to insufficient sperm quality remains an economic trait [21] cannot consider the pleiotropic effects of vari- hurdle for the sector [9]. ants or the interactions between them. In recent years, Currently, there are few genetic and transcriptomic different methods such as the associated weight matrix studies that have used high-throughput techniques to (AWM) [22] and partial correlation coefficient with investigate the genetic basis of sperm quality in swine. information theory (PCIT) [23] have been developed to To date, five genome-wide association studies (GWAS) carry out analysis of gene networks from GWAS or tran- have been performed. Diniz et al. [5] identified a single scriptomics data and to identify co-associated genes for quantitative trait locus (QTL) region associated to sperm a set of correlated phenotypes [22, 24–26]. Furthermore, motility in Large White pigs. Two years later, Zhao and the integration of GWAS and RNA-seq data can be used collaborators [10] reported three multi-single nucleo- to design knowledge-based technologies such as DNA tide polymorphism (SNP) QTL regions associated with marker panels including SNPs with a high functional epididymal weight, sperm concentration and total sperm potential for their application to animal breeding [27, 28]. per ejaculate, respectively and seven singleton QTL SNPs that display a genetic association with a phenotype related to sperm motility, semen temperature, seminifer- and show functional potential (e.g. coding or regulatory ous tubule diameter and number of ejaculates in a White variants) are less likely to show spurious associations Duroc × Erhualian ­ F2 population. Marques et al. [6] than non-functional SNPs. detected 16 and six QTL in Large White and Landrace Our aim was to identify candidate genes, pathways and pigs, respectively, associated with sperm motility, num- DNA variants associated to sperm quality in pigs by inte- ber of cells per ejaculate and morphological abnormali- grating GWAS and RNA-seq results under an unprec- ties. More recently, several QTL have been identified in edented systems biology approach. Moreover, we sought a Duroc population associated to number of sperm cells, to estimate the weight of the most relevant genes and sperm motility, sperm progressive motility, total mor- DNA variants on the sperm phenotypes. phological abnormalities, coiled tail, bent tail, proximal droplets, distal droplets and distal midpiece reflex [11, Methods 12]. Sample collection and phenotype measurements The presence of RNA molecules in the boar sperm is Three hundred fresh sperm ejaculates, each from a differ- well documented [13, 14], but their relation to sperm ent Pietrain boar from commercial farms, were collected quality is very little explored. Porcine sperm RNAs are by specialized professionals between September 2014 and highly fragmented and their gene abundances are mostly January 2017. Sperm was obtained using the gloved-hand associated to prior transcriptional events linked to sper- method [29], immediately diluted (1:2) in commercial matogenesis, fertility and embryo development [13]. extender and kept at 16°C for up to 2 h until phenotype A complex suite of RNAs are present in sperm, includ- assessment. Blood samples were collected from special- ing coding (mRNA), long noncoding RNAs (e.g. circu- ists during their routine sample collection and gDNA was lar RNA—circRNA-) and short noncoding RNAs (e.g. extracted using a phenol–chloroform based method [30]. microRNA –miRNA- or Piwi interacting RNA—piRNA-) The ejaculates were purified to remove somatic cells as stored with Trizol® at − 80°C until further use. [13]. Several studies have reported a relation between described previously [14] and purified spermatozoa were RNA abundances and semen quality in mammals [15–17]. In swine, Curry et al. [18] performed quantita- Phenotypic records from fresh sperm were meas- tive RT-PCR (RT-qPCR) that targeted 10 miRNAs and ured as previously described [14] and included: sperm Gòdia et al. Genet Sel Evol (2020) 52:72 Page 3 of 21 concentration (CON), percentage of viable cells (VIAB), as a regression coefficient on the corresponding (values percentage of morphologically abnormal acrosomes − 1, 0, 1) of the SNP i ; correcting for the fixed effect of (ACRO), osmotic resistance test (ORT), percentage farm (Farmj), season and year (SeasonYear k ) and boar of morphologically abnormal sperm cells (of the head age ( Agel ); um is the infinitesimal genetic effect of indi- -HABN-, neck -NABN- and tail -TABN-) and of cells vidual m, with u ∼ N 0, Gσu2 , where G is the genomic   with cytoplasmatic droplets (proximal -PDROP- and dis- relationship matrix (GRM) calculated using the filtered tal -DDROP-). Sperm motility traits were also assessed SNPs based on the methodology described by Yang et al. using the computer-assisted semen analysis (CASA) sys- [37], and σu2 is the additive genetic variance; and eijklm is tem (Integrated Sperm Analysis System V1.0; Proiser) the residual term. and included the percentage of motile spermatozoa cells The significance of SNP associations was corrected (MT) (with average path velocity -VAP- > 10 µm/s), aver- for multiple testing with the false discovery rate (FDR) age curvilinear velocity (VCL) (µm/s), average straight- approach [38] and only significant SNPs (FDR ≤ 0.05) line velocity (VSL) (µm/s) and average VAP (µm/s). All were kept for further analysis. Significantly associated phenotypes were assessed after 5 and 90 min of incuba- SNPs with consecutive distances shorter than 5 Mbp tion of the samples at 37°C, except for sperm concen- were considered to belong to the same GWAS interval tration, ORT, sperm abnormalities and cytoplasmatic [39]. A new interval was called if the consecutive SNPs droplets, which were measured only after 5 min of incu- were more than 5 Mbp apart. SNPs that mapped to the bation at 37°C. To calculate the correlations between sex chromosomes or to unmapped scaffolds were not RNA abundance and phenotype, sperm traits were cor- considered for further analysis. Genomic heritability rected using the fixed effects of farm of origin (3 levels), was assessed with GCTA v.1.91.5 through a genomic season and year of collection (9 levels) and boar age (3 restricted maximum likelihood (GREML) approach using levels) with the "lm" function of R [31] using a linear the GRM based on the methodology from Yang et al. model. The 90 min/5 min incubation ratios were also [37]. Manhattan plots of the GWAS results displaying calculated. In total, 25 phenotypic measures per sample the genetic associations (P-value) between each SNP and were recorded. Phenotypic correlations between traits phenotype were generated with the “qqman” R package were assessed and graphically displayed with the R pack- [40]. age “corrplot” [32]. The different analyses are described below, and the RNA isolation, sequencing and gene annotation complete outline is summarized in Additional file 1: Fig- RNA isolation from 40 sperm samples was performed as ure S1. previously described [14] and included 35 samples from boars analyzed in the GWAS. The other five boars did Genome‑wide association study (GWAS) not pass the genotyping quality control and thus were using the high-density (660 K markers) Axiom™ Porcine Two hundred and eighty-eight boars were genotyped not included in the GWAS. Extracted RNA was subjected the Qubit™ RNA HS Assay kit (Invitrogen), assessment to quality control assays including quantification with Genotyping Array (Thermo Fisher Scientific). The result- ing genotype dataset was stringently filtered by exclud- of RNA integrity with the 2100 Bioanalyzer using the ing the samples with a genotype call rate lower than Agilent RNA 6000 Pico kit (Agilent Technologies), and 96%. SNP locations were converted from Sscrofa10.2 evaluation by RT-qPCR of the sperm-specific PRM1, the to Sscrofa11.1 coordinates using plink v1.9 [33]. Then, somatic PTPRC mRNA and genomic DNA to confirm we excluded SNPs that (i) had a minor allele frequency that the samples were free from somatic cell RNA and lower than 0.05, (ii) deviated from Hardy–Weinberg gDNA contaminations. equilibrium (P-value ≤ 0.001), and (iii) for which there The ribosomal RNA (rRNA) from the 40 RNA sam- were more than 5% missing genotypes. These are stand- ples was depleted with the Ribosomal RNA depletion ard parameters that are typically used in similar studies Kit (Illumina) and libraries were prepared with the [34–36]. Single-SNP association analysis was carried out SMARTer Low Input Library Prep kit v2 (Clontech) using the genome-wide complex trait analysis (GCTA) and sequenced to generate 75 bp pair-end reads on v.1.91.5 software [37] with the following model: an Illumina’s HiSeq2000/2500. Undepleted total RNA was also subjected to short noncoding RNA (sncRNA) Yijklm =µ + δSNPi + Farmj + SeasonYeark library preparation (34 of the previous 40 samples) + Agel + um + eijklm , using the NEBNext library prep kit (New England Bio- labs) and sequenced at 50 bp single-end on a Hiseq2000 where Yijkl is the phenotype modeled as a function of the (Illumina). population mean (µ), δ is the SNP allelic effect, estimated Gòdia et al. Genet Sel Evol (2020) 52:72 Page 4 of 21 Total RNA-seq reads were evaluated for quality con- were calculated using the Pearson correlation coefficient. trol with the FastQC software (https​://www.bioin​forma​ Only correlations with a P-value ≤ 0.05 were kept. tics.babra​ham.ac.uk/proje​cts/fastq​c/). Low-quality reads (phred –Q < 20 and read length < 25 bp) and sequenc- SNP calling from RNA‑seq data and linkage disequilibrium ing adaptors were trimmed with Trimmomatic v.0.36 with GWAS lead SNPs [41]. Filtered reads were mapped to the porcine genome Mapped RNA-seq reads of the 35 samples with RNA-seq (Sscrofa 11.1) using HISAT2 v.2.1.0 [42]. Duplicate reads and genotype data were subjected to SNP calling. Vari- were removed with Picard Tools v.2.18.29 (http://picar​ ant calling was performed with SAMtools mpileup and d.sourc​eforg​e.net) Markduplicates. RNA levels of the BCFtools v.1.9 [49]. Only SNP variants for which the genes annotated in the porcine genome (Ensembl v.91) alternative allele was present in at least 10 samples with were then quantified with StringTie v.1.3.4 [43]. Only the a minimum Phred quality of 25 and a minimum read genes with average RNA abundances ≥ 10 fragments per depth of 10 were kept. The effect of the SNP on protein kb of exon per million reads mapped (FPKM) were kept sequence was predicted with SnpEff v.4.3T [50] and only for further analysis with the aim to discard low abundant low, moderate and high impact variants were kept. The genes and spuriously mapped reads. new SNP genotypes were merged to the Axiom geno- The effect of external variables on gene expression was types and the linkage disequilibrium (LD) ­R2 between assessed using the following mixed effect model as in GWAS lead SNPs and RNA-seq SNPs was assessed with Reverter et al. [44]: PLINK v1.9 [33] using the default parameters, with the Yijklmn =µ + Li + Gj + GFjk exception of “-ld-window 100”, “-ld-window-kb 0” and “– ld-window-r2 0” to assess all the pair-wise LD values. + GYSjl + GAjm + GRjn + eijklmn , where Yijklmn represents the log2-transformed FPKM Expression GWAS value from library i (40 levels), gene j (4120 levels), farm Expression GWAS (eGWAS) included the 35 samples k (3 levels), year-season l (6 levels), age m (3 levels) and with RNA-seq and genotype data. The RNA abundances assay run n (4 levels). Accordingly, Yijklmn was mod- of the detected genes were taken as quantitative traits eled as a function of the mean (µ), fixed effect of library and tested for association with the genotypes that passed ( Li ) and the random effects of gene (Gj ), gene by farm quality control using a linear model. Single-SNP associa- (GF jk ), gene by year-season (GYS jl ), gene by age (GAjm) tion analysis was performed with the GCTA v.1.91.5 soft- and gene by assay run (GRjn). Random residuals in eijklmn ware [37], with the following model: were assumed to be independent and identically distrib- Yi = µ + SNPi + ei , uted. Using standard stochastic assumptions, the effects of Gj , GF jk , GYS jl , GAjm and GRjn were assumed to fol- where Yi is the log2-transformed gene abundance mod- low a normal distribution with zero mean and between- eled as a function of the population mean (µ), fixed effect gene, between-gene within-farm, between-gene within of each SNP ( SNP i ), and a residual effect (ei). year-season, between-gene within age, between-gene eGWAS significant associations (FDR ≤ 0.05) were con- within assay and within gene components of variance, sidered only if: (i) the eGWAS associated SNP was also a respectively. Restricted maximum likelihood estimates significant hit (FDR ≤ 0.05) in the GWAS for sperm qual- and solutions to model effects were obtained using VCE6 ity phenotypes and (ii) the gene’s RNA abundance corre- [45]. lated to the same phenotype as the corresponding GWAS For the sncRNA-seq data, trimming of adaptors and SNP hit. low-quality bases (phred –Q < 20 and read length < 12 bp) was performed with Cutadapt v1.0 [46]. Reads were SNP co‑association and gene co‑abundance analyses mapped to the Sus scrofa genome (Sscrofa11.1) with the We also carried out a SNP co-association analysis by sRNAtoolbox v.6.17 [47] using default settings and with building an AWM from the GWAS results [22, 51]. The the porcine miRBase [48] release 21 database. Multi- AWM was constructed from two matrices that contained adjusted read counts were normalized by library size as row-wise SNPs and column-wise phenotypes. The first counts per million (CPM). Only miRNAs with an average matrix included the P-values of the association between abundance higher than 1 CPM in all the samples were each SNP and the phenotype, and the second matrix considered. miRNA abundance was stabilized with the corresponded to the SNP z-score standardized additive log2 transformation. effect. As live cells with intact plasma membrane are The relationship between the 25 phenotypes and each essential for fertilization [52, 53], the percentage of via- of the log2-stabilized mRNA’s and miRNA’s abundances ble spermatozoa at 5 min (VIAB_5) was selected as key Gòdia et al. Genet Sel Evol (2020) 52:72 Page 5 of 21 phenotype and the associated SNPs (P-value ≤ 0.01) were genes was a better predictor of sperm quality phenotypes. included in the AWM. In the next step, the dependency For this, first we extracted 20 genes from the network. between phenotypes was estimated based on the aver- These genes were (i) correlated with at least four pheno- age number of non-key phenotypes associated (Ap) with types, (ii) did not present interactions (edges) between these SNPs (P-value ≤ 0.01) (Ap ≥ 2). Then, SNPs that them, (iii) all samples presented RNA abundance lev- were located less than 2500 bp or more than 1 Mbp from els higher than 0 FPKM, and (iv) were potentially rel- the nearest annotated gene (Ensembl v.91) were kept. The evant according to the existing literature. The RSQUARE most significant SNP from each annotated gene was kept statement of the REG procedure implemented in the to build the AWM. The standardized SNP effects across SAS software [57] was used as an exploratory model to phenotypes were computed and represented using the evaluate all possible subsets of linear regressions using hierarchical cluster analysis based on Euclidean distance unadjusted gene abundances and sperm phenotypes and with the R package “dendextend” [54]. Then, significant extract the ­R2 magnitude from each prediction. Then, we gene–gene interactions were assessed to build the SNP selected the subset of 10 genes that were most commonly network with the PCIT algorithm [23]. PCIT applies present in all the phenotype models. This subset of com- first-order partial correlation coefficients together with mon genes was then used for the STEPWISE statement an information theory approach to identify meaning- of the REG procedure implemented in the SAS software ful gene–gene associations [23]. Only significant gene [57], which performs a linear regression analysis for each co-associations (P-value ≤ 0.05) were kept in the SNP of the phenotypes to develop a model to predict the phe- network. notype based on gene RNA levels. The model is: For the RNA co-abundance analysis, significant gene– Yij = intercepti + GPEij + eij , gene interactions that were used to build the RNA net- work were also predicted with PCIT using the stabilized where Yij represents the predicted phenotype value from RNA abundances. Interactions between genes and miR- i -th phenotypes (25 levels), j-th genes (10 levels). Yij was NAs were also assessed with PCIT [23], and only signifi- modeled as a function of the intercept value for the phe- cant negative correlations (P-value ≤ 0.05) were kept. notype (intercept i ), the gene abundance by parameter estimate (GPE ij ) and a residual term (eij ). The accuracy Integration of SNPs and RNA network data and network was ascertained from the model’s goodness-of-fit and visualization based on the proportion of variance explained by the The genes and interactions that were present in both the model ­(R2). We also developed a genome-wide SNP panel GWAS and RNA-seq networks were considered to have a to identify the SNPs that could best predict the pheno- higher probability of being involved in sperm quality and typic variance of sperm-related traits. The panel included were used to build a robust gene interaction network. the lead SNPs from the GWAS and from the eGWAS The resulting network was named “shared network”. In hits, and the GWAS most significant SNP for each of the addition, the genes that were not present in the shared genes included in the network that also: (i) correlated network but that presented an abundance correlation with at least four phenotypes and (ii) were identified in with more than three phenotypes and their co-associated the shared network. The proportion of the phenotypic genes were merged with the shared network to create the variance explained by these SNPs was estimated with the so-called final network. This final network also included GREML analysis implemented in the GCTA software the interactions between miRNA and mRNA genes. using the GRM calculated with the 73 autosomal SNPs Network visualization was performed with Cytoscape based on the methodology from Yang et al. [37]. v3.6 [55] and included information on: (i) the number of phenotypes associated to a gene or miRNA, (ii) the phe- Results notype with the highest correlation for each gene, (iii) Phenotypic parameters whether the gene was annotated as a transcription factor Three hundred ejaculates were phenotyped for 25 sperm (TF) or TF co-factor, and (iv) whether the gene was pre- quality traits (Table 1). Phenotype correlations (see Addi- sent in the shared network or was only found in the final tional file 2: Figure S2) were consistent with their physi- network. TF and TF co-factors were extracted from the ological similarities. In general, SNP-based heritabilities AnimalTFDB3.0 database [56]. (Table 1) were low to moderate with motility-related traits displaying higher values. MT_90 was the most her- Development of an RNA model and SNP panel itable trait (­ h2: 0.39), whereas motility ratios, NABN and for the phenotypic prediction of sperm quality VIAB_5 showed heritability values close to 0 (Table 1). The unadjusted RNA abundance of a subset of the genes The sperm phenotypes correlated with farm, boar age in the network was used to identify which combination of and season per year (see Additional file 3: Table S1) and Gòdia et al. Genet Sel Evol (2020) 52:72 Page 6 of 21 Table 1 Descriptive statistics, genomic heritability ­(h2) and number of significant SNPs in the GWAS for sperm quality parameters (N = 300) Trait Acronym Mean (SD) h2 (SE) Number of SNPs Number of SNPs in autosomal in unplaced chromosomes scaffolds Concentration (sperm/mL) CON 141.3 (65.5) 0.13 (0.11) 0 0 Viability 5 min VIAB_5 90.1 (6.3) 1 × 10–6 (0.11) 0 0 Viability 90 min VIAB_90 77.4 (17.3) 0.14 (0.13) 0 0 Osmotic resistance test ORT 79.8 (12.5) 0.13 (0.12) 0 0 Head abnormalities HABN 2.1 (5.9) 0.16 (0.11) 41 0 Neck abnormalities NABN 3.0 (4.9) 1 × 10–6 (0.13) 18 0 Tail abnormalities TABN 2.7 (3.4) 0.09 (0.12) 0 0 Proximal droplets PDROP 3.5 (5.1) 0.12 (0.15) 1 0 Distal droplets DDROP 4.5 (4.5) 0.06 (0.11) 0 0 Motility 5 min MT_5 75.4 (18.1) 0.21 (0.15) 3 217 Motility 90 min MT_90 64.1 (22.0) 0.39 (0.14) 2 252 Average path velocity 5 min (µm/seg) VAP_5 34.0 (10.2) 0.17 (0.11) 0 0 Average path velocity 90 min (µm/seg) VAP_90 30.8 (9.5) 0.35 (0.13) 0 0 Curvilinear velocity 5 min (µm/seg) VCL_5 46.2 (12.5) 0.11 (0.10) 0 0 Curvilinear velocity 90 min (µm/seg) VCL_90 39.7 (10.2) 0.35 (0.13) 0 0 Straight line velocity 5 min (µm/seg) VSL_5 27.0 (8.3) 0.23 (0.13) 0 38 Straight line Velocity 90 min (µm/seg) VSL_90 25.9 (8.3) 0.34 (0.13) 0 0 Abnormal acrosomes 5 min ACRO_5 7.0 (5.6) 0.08 (0.11) 4 0 Abnormal acrosomes 90 min ACRO_90 16.4 (12.6) 0.06 (0.10) 0 0 Ratio motility R_MT 0.9 (0.2) 1 × 10–6 (0.11) 0 0 Ratio average path velocity R_VAP 0.9 (0.3) 1 × 10–6 (0.08) 0 0 Ratio Curvilinear velocity R_VCL 0.9 (0.3) 1 × 10–6 (0.09) 0 0 Ratio straight line velocity R_VSL 1.0 (0.3) 0.06 (0.10) 0 0 Ratio viability R_VIAB 0.9 (0.3) 0.08 (0.11) 0 0 Ratio acrosomes R_ACRO 3.4 (3.5) 0.08 (0.11) 1 0 All traits except stated are presented as a percentage Number of SNPs = GWAS number of single nucleotide polymorphisms significantly associated (FDR) with the trait The values shown are raw excepting the ratios which were previously corrected and stabilized SD standard deviation, SE standard error were thus included as fixed effects in the GWAS model displayed significant associations (FDR ≤ 0.05) for each and phenotypes were also corrected for these effects to trait is summarized in Table 2. carry out the correlation analysis. Seven sperm quality traits exhibited significant asso- ciation signals (Fig. 1a–g) and (see Additional file 4: GWAS analysis Table S2), and only one SNP was associated with more After quality control, 466,592 SNPs and 276 samples than one trait (Table 2; Fig. 1d, e) and (see Additional remained for the GWAS. In total, 324 SNPs across the file 4: Table S2). The number of SNP signals was larg- autosomal chromosomes and unplaced scaffolds dis- est for HABN and NABN with 41 and 18 associated played genetic associations (FDR ≤ 0.05) with one or SNPs, respectively (Fig. 1a, c) and (see Additional more sperm quality phenotype (Table 1) and (see Addi- file 4: Table S2). Six of the 19 QTL were represented tional file 4: Table S2). Among these 324 SNPs, 255 by one associated SNP only and were discarded from mapped to unplaced scaffolds and were not considered further analyses (Table 2; Fig. 1). The most significant for further data analysis (Additional file 4: Table S2). SNPs (rs318575212 and rs332927981) were associated Nineteen chromosomal regions tagged by 69 significant with ACRO_5 (both with FDR = 0.006 and an additive SNPs were identified on Sus scrofa (SSC) chromosomes effect = 4.11) (Table 2). 1, 3, 4, 6, 7, 9, 13 and 16. The number of SNPs that Gòdia et al. Genet Sel Evol Table 2 Summary of the results of the genome wide association analysis for sperm quality traits SSC Interval #SNP Interval Mbp Top SNP Top SNP location bp Top SNP P-value Top SNP FDR Top SNP MAF Beta Trait (2020) 52:72 1 I1 1 – rs339761632 13,501,755 4.64 × 10–8 0.02 0.06 4.84 PDROP 1 I2 8 82.90–83.49 rs81354986 82,895,619 1.69 × 10–6 0.03 0.07 5.05 HABN 1 I3 8 94.88–98.74 rs327733412 94,880,167 1.61 × 10–7 0.02 0.07 5.65 HABN 1 I4 1 – rs337166779 126,397,198 2.05 × 10–6 0.03 0.06 5.02 HABN 1 I5 11 243.86–246.44 rs343194423 246,224,386 1.72 × 10–7 0.01 0.07 3.17 NABN 1 I6 2 258.54–258.55 rs332256425 258,548,786 1.76 × 10–6 0.04 0.06 3.44 NABN 3 I1 1 – rs332055717 2,911,413 6.35 × 10–8 0.01 0.09 5.07 HABN 3 I2 3 113.75–113.84 rs328292697 113,750,595 1.09 × 10–7 0.01 0.07 3.41 NABN 4 I1 2 2.41–2.42 rs318575212 2,412,006 2.88 × 10–8 0.01 0.08 4.11 ACRO_5 rs332927981 2,415,239 6 I1 2 65.60–66.66 rs335394654 65,597,553 1.86 × 10–7 0.03 0.14 3.04 ACRO_5 7 I1 2 6.20–6.38 rs326239534 6,377,172 9.87 × 10–6 0.02 0.17 -9.15 MT_5 7 I2 2 85.73–86.88 rs336588919 86,884,279 4.13 × 10–8 0.01 0.06 3.75 NABN 9 I1 2 5.76–5.78 rs1110111787 5,776,597 1.55 × 10–7 0.02 0.07 5.43 HABN 9 I2 1 – rs342738178 28,463,580 1.53 × 10–5 0.03 0.14 − 10.42 MT_5, MT_90 9 I3 1 – rs328217450 137,959,590 4.77 × 10–8 0.02 0.18 2.36 R_ACRO 13 I1 18 25.36–28.47 rs690794887 25,535,100 3.06 × 10–7 0.02 0.14 3.78 HABN 13 I2 3 33.82–37.65 rs327865244 33,819,549 3.79 × 10–8 0.01 0.15 4.28 HABN 16 I1 1 – rs324239602 6,476,358 6.08 × 10–6 0.01 0.46 9.07 MT_90 SSC, Sus scrofa chromosome; #SNP, number of SNPs significantly associated (FDR) with the trait; Interval, region of the GWAS interval; Beta, additive effect; FDR, false discovery rate; MAF, minor allele frequency; ACRO_5, abnormal acrosomes 5 min; HABN, head abnormalities; NABN, neck abnormalities; PDROP, proximal droplets; R_ACRO, ratio acrosomes; MT_5, motility 5 min; MT_90, motility 90 min Page 7 of 21 Gòdia et al. Genet Sel Evol (2020) 52:72 Page 8 of 21 Sperm RNA isolation, RNA‑seq and bioinformatics analysis Table S5). Twenty-three SNPs were predicted to have a Isolated RNA from mature spermatozoa was free from low effect on protein sequence (synonymous variants somatic cell RNA. Total RNA-seq resulted in an aver- and 5′ UTR premature start codon), 13 SNPs showed a age of 40.7 M reads per sample and 98.2% of the reads moderate effect (missense variants) and one SNP was passed the quality control filters (see Additional file 5: predicted as a splice donor variant and thus, to have a Table S3). An average of 83% of the reads mapped to the high impact on protein sequence (see Additional file 7: porcine genome and after removal of duplicate reads and Table S5). RNA abundance filters, we identified 4120 genes (see Interval 1 (I1) on SSC13 was associated to HABN, Additional file 6: Table S4). The most abundant protein harboured 21 expressed SNPs (7 and 14 with moder- coding transcripts included PRM1, OAZ3, DNAJB8, ate and low effects, respectively). The polymorphism ANKRD35, SPATC1 and ZNRF4, among others, as well rs331304027 (a missense variant with a moderate effect as mitochondrial genes such as COX1, CYTB and ND5 on the ULK4 gene) was in LD (LD = 0.40) with the (see Additional file 6: Table S4). The variance compo- strongest GWAS SNP hit of the interval (rs690794887) nent estimated by the mixed model explained 84% (80% (Table 3). SSC13 I2, was also associated to HABN, due to the main effect of the gene) of the variation in included 11 SNPs (1 with a high, 5 with a moderate and gene abundance. Consequently, RNA abundances were 5 with a low effect on protein sequence). Of these, the not corrected for external effects. For short RNA-seq, variant with the highest LD (LD = 0.2) with the GWAS we obtained an average of 7.3 M of reads per sample. Of hit (rs327865244) was a 5′ UTR premature start codon these, 99.2% passed quality control and 81.5% mapped to gain (low effect) SNP (rs323872641) in the ABHD14A the porcine genome (see Additional file 5: Table S3). In gene (Table 3) and (see Additional file 7: Table S5). This more detail, 42% of the aligned reads corresponded to interval was the only one that presented a SNP with a sncRNAs, including piRNAs (16%), snRNA (8%), tRNA high effect (novel), a splice donor variant in the IQCF5 (9%) and miRNA (9%) (see Additional file 5: Table S3). gene, with almost no LD (LD = 0.02) with the GWAS The remaining aligned reads corresponded to mitochon- hit (see Additional file 7: Table S5). SSC7 I2 was associ- drial transfer and ribosomal RNAs (see Additional file 5: ated to NABN and encompassed two expressed SNPs Table S3). We identified 95 miRNAs out of the 306 that (both with a low effect). rs330912302 (a synonymous are annotated in swine (see Additional file 6: Table S4). SNP in the CHD2 gene) presented an LD (LD = 0.4) with The most abundant miRNAs with CPM higher than 1000 the strongest hit of the interval (rs336588919) (Table 3). were ssc-miR-30d, ssc-miR-34c, ssc-miR-30c-5p, ssc- SSC1 I3 was associated to HABN and harboured three miR-191, ssc-let-7a, ssc-let-7g, ssc-miR-28-3p and ssc- expressed SNPs (1 with a moderate and two with a low miR-425-5p (see Additional file 6: Table S4). effect) (Table 3) and (see Additional file 7: Table S5). SNP calling from RNA‑seq and linkage disequilibrium Correlation of genes’ and miRNAs’ abundances with sperm with GWAS hits quality traits Under the hypothesis that some of the GWAS hits may The correlation analysis of the 4120 genes and the tag a causal variant that alters the protein sequence and 25 phenotypes resulted in 6128 significant correla- function, and to identify additional SNPs with the poten- tions (P-value ≤ 0.05) involving 3007 genes and the 25 tial to obtain better genetic markers than those identified traits (see Additional file 8: Table S6). These genes pre- in the GWAS, we sought to identify variants in annotated sented between one and nine significant correlations genes using the RNA-seq data. As a prerequisite, these with the different semen quality traits (see Additional variants had to be in LD with the cognate GWAS hit. file 8: Table S6). Three hundred and forty-four genes After filtering, we identified 7719 expressed variants, 37 were significantly correlated with more than four traits. of which were located within the genomic intervals iden- For the miRNAs, the abundance of the 95 miRNAs tified in the GWAS (Table 2) and (see Additional file 7: and the studied phenotypes resulted in 306 significant (See figure on next page.) Fig. 1 Manhattan plots depicting the genetic associations between SNPs and the sperm quality traits that showed genome-wide significant values. Significant associations have been found with the percentage of: a Percentage of cells with head abnormalities (HABN); b percentage of cells with abnormal acrosomes after 5 min incubation at 37°C (ACRO_5); c percentage of cells with neck abnormalities (NABN); d percentage of motile spermatozoa after 5 min incubation at 37°C (MT_5); e Percentage of motile spermatozoa after 90 min incubation at 37°C (MT_90); f Percentage of cells with proximal droplets (PDROP); g Ratio of the percentage of abnormal acrosomes at 5 min versus 90 min incubation times (R_ACRO). The x-axis represents chromosome length (Mb), and the y-axis shows the negative l­og10 P-values of the genetic associations. The horizontal red line represents the significance threshold (FDR ≤ 0.05) Gòdia et al. Genet Sel Evol (2020) 52:72 Page 9 of 21 Gòdia et al. Genet Sel Evol (2020) 52:72 Page 10 of 21 Table 3 Summary of the SNPs identified from the RNA-seq datasets in genes mapping within the GWAS regions SSC Interval Top SNP # SNPs called Highest LD SNP with highest LD Genotypic # called SNP effect Gene Trait of the GWAS frequency (0/0; samples interval 0/1; 1/1) 1 I3 rs327733412 3 0.07 rs710447566 0.34; 0.54; 0.11 35 Low KATNAL2 HABN 7 I2 rs336588919 2 0.4 rs330912302 0.63; 0.12; 0.25 32 Low CHD2 NABN 13 I1 rs690794887 21 0.4 rs331304027 0.06; 0.09; 0.85 33 Moderate ULK4 HABN 13 I2 rs327865244 11 0.2 rs323872641 0.49; 0.37; 0.14 35 Low ABHD14A HABN SSC, Sus scrofa chromosome; # SNPs called, number of SNPs identified in the SNP calling analysis; LD, linkage disequilibrium; Genotypic frequency: allelic frequency for each of the genotypes; # called samples, number of samples with reads in the given SNP position; HABN, head abnormalities; NABN, neck abnormalities The columns SNP effect and gene refer to the SNP with the highest LD in the region correlations (P-value ≤ 0.05) which involved 87 miRNAs detected with the AWM explained 74.1% of the pheno- and 17 semen traits (see Additional file 9: Table S7). The typic variance of the key phenotype (VIAB_5). The SNP miRNAs presented between 1 and 9 significant correla- network predicted with PCIT [23] resulted in significant tions with the semen quality traits (see Additional file 9: correlations that involved 2648 nodes (all the genes) con- Table S7). nected by 2,984,616 edges (Fig. 2). For the RNA network analysis, the RNA levels of the Expression GWAS analysis 4120 detected genes were used to identify potential con- In order to predict whether the GWAS hits tagged a nections using PCIT [23]. The RNA network included causal variant that altered gene expression, we performed 4120 nodes (all the genes) connected by 1,173,995 edges a within-trait eGWAS with the genotypes of 464,020 (Fig. 2). PCIT also built 4539 significant interactions SNPs that passed the quality control and the normal- between 95 miRNAs and 630 genes. ized RNA abundances. Then, we focused only on the To obtain the shared network, common SNP and RNA associations between GWAS SNP hits (with FDR ≤ 0.05) network edges were extracted, thus, focusing only on the and transcripts with abundances that correlated with the shared set of interacting genes from both approaches. same phenotype. We identified 45 SNPs (FDR ≤ 0.05) that This comparison resulted in 613 nodes connected by were located in three genomic regions related to ACRO_5 16,591 edges (Fig. 2). The final network included a set and HABN (Table 4). Six SNPs had unknown positions of 344 additional genes (since they correlated with on the genome after the lift-over from Sscrofa10.2 to more than three phenotypes) and their interactions, Sscrofa11.1. The remaining eGWAS hits were on SSC4, 6 resulting in 700 genes. Moreover, the final network also and 13 (Table 4) and (see Additional file 10: Table S8). All involved 1564 edges connecting 202 genes and 94 miR- the SNPs had a trans effect on genes that were located on NAs (Fig. 2). Of the 1313 genes included in the final other chromosomes. The eQTL identified on SSC4, was network, 1135 had an abundance that correlated with at related to ACRO_5 and was associated to three genes, least one phenotype, 68 have been reported as TF and 89 NCLN, ASCC1 and AATF. The eQTL on SSC6 was also as TF co-factors (Fig. 3a). Nearly a quarter of the genes related to ACRO_5 and was associated to the IQCJ gene. (282 of the 1313) presented at least 200 edges. The genes Finally, the eQTL on SSC13 for HABN, included SNPs that presented the largest number of interactions were that were in the HARS, ACTR2, EPB41L3 and RAB1B PLCH2 (579 edges, present in the final but not in the genes. shared network and correlated with three phenotypes), CEP152 (399 edges, in the shared network and correlated Gene network analysis with four traits) and SLC41A2 (382 edges, in the shared After SNP selection, 2648 of the 466,592 SNPs were network). retained to build the AWM. Trait hierarchical cluster Gene ontology analysis of the genes included in the distributions were in agreement with the biological simi- final network presented enrichment for DNA repair (e.g. larities and phenotypic correlations (see Additional file 2: RAD51, SETX, and SOD1), meiotic cell cycle (e.g. BAG6, Figures S2 and S3). A clear separation between (i) mor- HSPA2, and RAD51), gamete generation (e.g. TSSK3, phological abnormalities and motility parameters and (ii) PRDM14, and PRKAR1A) and spermatogenesis (e.g. cell viability and ORT was observed based on the additive BAG6, CAPZA3, and HSPA2) (see Additional file 11: effects of the SNPs calculated in the association analy- Table S9). sis. Consistent with previous studies [58, 59], the SNPs Gòdia et al. Genet Sel Evol (2020) 52:72 Table 4 Summary of the results from the within-trait expression genome wide association analysis SSC Interval # SNP: Top eGWAS Top eGWAS Top eGWAS P-value Top eGWAS FDR Top Beta Trait RNA Associated gene transcripts location bp eGWAS abundance MAF correlation 4 I1 2 rs318575212 2,412,006 7.36 × 10–3 0.03 0.09 − 0.39 ACRO_5 − 0.33 NCLN rs332927981 2,415,239 2 rs318575212 2,412,006 1.83 × 10–4 0.03 0.09 − 1.8 ACRO_5 − 0.46 ASCC1 rs332927981 2,415,239 2 rs318575212 2,412,006 2.87 × 10–4 4.83 × 10–2 0.09 − 1.1 ACRO_5 − 0.4 AATF rs332927981 2,415,239 6 I1 2 rs335394654 65,597,553 5.63 × 10–5 0.02 0.11 − 1.65 ACRO_5 − 0.35 IQCJ 13 I1 31 rs328397029 25,684,259 1.84 × 10–5 2.95 × 10–3 0.09 − 1.03 HABN − 0.38 HARS, ACTR2, EPB41L3, RAB1B SSC, Sus scrofa chromosome; # SNP: transcripts, number of single nucleotide polymorphisms significantly associated to a transcript; Beta, additive effect; MAF, minor allele frequency; ACRO_5, abnormal acrosomes 5 min; HABN, head abnormalities Page 11 of 21 Gòdia et al. Genet Sel Evol (2020) 52:72 Page 12 of 21 Discussion Investigating the genomic regions and molecular pro- cesses that control sperm quality has become a focus of interest in humans and in livestock including swine, in the latter case for its relevance on the sustainability of pig breeding and production [5, 6, 10–12]. In fact, our results and those obtained by other groups [5–8], have shown that boar sperm quality has a genetic basis, which means that it can be selected for in breeding strategies. Here, we provide an exploratory analysis using multiple bioinformatics tools. Since our study was carried out on a relatively small sample size with one phenotypic ejacu- late evaluated per boar, and sperm quality traits are influ- enced to a great extent by environmental factors [60], our results should be considered as preliminary. This is the first study that explores at the genomic level the molec- Fig. 2 Number of nodes (genes) in each of the gene network ular components of sperm and semen quality using an analyses. The SNP network involved 2648 nodes connected by 2,984,616 edges (interactions). The RNA Network included 4120 integrative approach that fits GWAS and RNA-seq data. nodes connected by 1,173,995 edges. The shared network included Moreover, our study includes, for the first time, traits the 613 nodes and 16,591 edges present in both the SNP and the such as ORT, ACRO, VIAB and the dissection of mor- RNA networks. The final network included (i) the shared network, phological abnormalities of different parts of the sperm (ii) 700 additional genes corresponding to genes that correlated cell (HABN, NABN and TABN). with more than four traits and their interacting genes (iii) as well as 94 co-associated miRNAs. These miRNAs interacted with 202 nodes involving 1564 edges GWAS analysis The GWAS revealed 12 QTL that were represented by two or more significant SNPs and several positional can- Development of an RNA model and a SNP panel didate genes for HABN, NABN, ACRO_5 and MT_5 The ­R2 model predicted that the RNA levels of 20 genes (Table 2). The highest signals were observed on SSC4 for could explain between 55 to 78% of the phenotypic vari- ACRO_5 (~ 2.41–2.42 Mbp) (Table 2) and (see Additional ation across traits. The 10 selected genes that were most file 4: Table S2), ~ 69 kb upstream of the solute carrier commonly present in all the phenotype models explained family 45 member 4 (SLC45A4) gene. SLC45A4 encodes the vast majority (93 to 99%) of the phenotypic variation a proton-coupled sugar transporter that plays a role in that was predicted by the model. The final set of 10 genes the nutrition of spermatozoa during their maturation included in the linear regression model was composed in the epididymis [61] where acrosome assembly con- of MICAL3, EFHC1, TRAPPC2L, ATP9A, THADA, tinues during the post-testicular maturation phase [62]. MOBKL3, BLVRB, LARP4, CARS2, and NDUFV2. The Another solute carrier family 35 member B3 (SLC35B3) analysis resulted in significant models for 10 of the 25 was selected as a potential candidate for the MT_5 QTL phenotypes (Table 5). The most significant model was for on SSC7 (Table 2) and (see Additional file 4: Table S2). PDROP, which could predict the phenotype with an effi- SLC35B3 maps 0.6 Mbp away from this QTL. ciency of 68% (Table 5). The estimated parameters of the We detected several significant regions for HABN significant models are in Additional file 12: Table S10. (Table 2). Interestingly, HABN showed little correla- The SNP-based panel was built with 73 SNPs (18 lead tion with the other phenotypes (see Additional file 4: SNPs from GWAS hits, 2 lead SNPs from the eGWAS Table S2), but the biological rationale behind this remains hits, 53 SNPs from the shared network and correlated ≥ 4 to be elucidated. The QTL on SSC1 I2 (~ 94.9–98.8 Mbp) phenotypes) (see Additional file 13: Table S11). These included interesting candidate genes such as the kata- polymorphisms could explain between 5 and 36% of nin catalytic subunit A1 like 2 (KATNAL2). Dunleavy the phenotypic variance across the 25 traits (Table 5). A et al. [63] reported that, in mice, Katnal2 is a critical moderate proportion (> 20%) of the phenotypic variance regulator of male germ cell development by affecting could be explained for 18 of the 25 traits. The best pre- sperm head shaping, acrosome attachment and sperm dictions were for sperm abnormalities (NABN, HABN, tail growth. Other candidate genes in that region were TABN) and sperm motility related traits (e.g. MT_5, SLC14A2, encoding the urea transporter A, which has VAP_90 and VCL_90) (Table 5). been suggested to participate in sperm head forma- tion by reducing its volume though excretion of urea Gòdia et al. Genet Sel Evol (2020) 52:72 Page 13 of 21 Fig. 3 Co-association network based on the AWM and transcriptomics data. a Full network with 1313 genes and 94 miRNAs; b Subset of the network showing the transcription factor CARF and all its predicted interactions; c Subset of the network with the TRAPPC2L interactions, which included several miRNAs; d Subset of the network with the CHD2 gene interactions. The node color corresponds to the phenotype group with the highest correlation value, as follows: concentration (red), abnormal acrosomes (green), abnormalities and droplets (pink), osmotic resistance test (orange), motility (light blue) and viability (dark blue). miRNAs are depicted in yellow. Node and text sizes correspond to the number of significant phenotypes correlated with that gene or miRNA. Nodes with a black line border correspond to genes identified in the shared network. Node shape indicates classification as: triangle (TF), V (TF co-factor) and ellipse (other genes and miRNAs) [64], and the SMAD family member 2 (SMAD2), which ABCA1 has been suggested to contribute to cholesterol is involved in spermatogonial differentiation [65]. On transport and fertilization capacity [70]. SSC13 I1, we identified two candidate genes: the testis Four of our GWAS hits are located near previously and ovary-specific PAZ domain gene 1 (TOPAZ1) and reported QTL for semen quality traits. This is the the IQ motif containing F1 (IQCF1). Luangpraseuth- case for the SSC1 I6 QTL, associated to NABN, which Prosper et al. [66] demonstrated that Topaz1-knockout mapped 335 kbp downstream from a QTL associ- mice presented meiotic arrest and male infertility. As for ated to sperm abnormalities and motility in boars [6]. IQCF1, Fang et al. [67] reported that this gene localizes The QTL SSC3 I2, associated to NABN lies 350 kpb in the acrosome and that it is involved in sperm capacita- upstream from a PDROP QTL [12]. The SSC4 I1 QTL, tion in mice. Iqcf1−/− mice were significantly less fertile associated to ACRO_5, resides 655 kpb upstream from than wild type mice [67]. The QTL region on SSC13 I2 a QTL for the distal midpiece reflex [12] and the SSC7 included the candidate protein kinase C delta (PRKCD) I1 QTL, associated to MT_5 maps 123 kpb upstream gene. PRKCD is involved in spermatogenesis and embry- from a PDROP QTL [12]. These discrepancies across onic development [68] and was highlighted in a GWAS studies could arise due to different technical (e.g. sam- for semen volume in Holstein–Friesian bulls [69]. ple size, SNP arrays, QTL or phenotyping accuracy), Four QTL regions were identified for NABN (Table 2) environmental (e.g. temperature, animal husbandry or and (see Additional file 4: Table S2). The QTL on SSC1 sperm processing), or biological factors (e.g. genetic I5 included as candidate gene the transporter ATP heterogeneity). binding cassette subfamily A member 1 (ABCA1). In humans, ABCA1 localizes on the dorsal side of the SNP calling from RNA‑seq data sperm head and in the middle piece of the tail [70]. Calling genomic variants from RNA-seq data can be a complementary method to detect previously unknown or Gòdia et al. Genet Sel Evol (2020) 52:72 Page 14 of 21 Table 5 R2 and phenotypic variance for each trait (N = 35). Consequently, although we used stringent cri- from the RNA model and SNP panel teria for genotype calling, these results should be consid- Acronym RNA model SNP panel ered as merely indicative and a larger number of samples 2 should be analysed to draw more robust conclusions. R P-value Phenotypic variance For HABN, we found new genetic variants in genes explained (SE) of physiological interest (Table 3) and (see Additional file 7: Table S5). On SSC13 I1, we detected several vari- CON 0.17 0.82 0.05 (0.05) ants in the unc-51 like kinase 4 (ULK4) gene, which VIAB_5 0.43 0.06 0.27 (0.07) is in potential LD with the lead SNP of this GWAS hit VIAB_90 0.23 0.61 0.28 (0.07) (Table 3) and (see Additional file 7: Table S5). Although ORT 0.22 0.62 0.24 (0.07) ULK4 has not been directly linked to sperm traits, Liu HABN 0.16 0.84 0.29 (0.06) et al. [71] showed that this gene has an essential role NABN 0.22 0.64 0.36 (0.07) in ciliogenesis, a process that is also crucial in sperm. TABN 0.26 0.49 0.26 (0.07) The previously discussed GWAS positional and physi- PDROP 0.68 < 0.0001 0.17 (0.07) ological candidate genes CHD2 and KATNAL2, also DDROP 0.42 0.07 0.06 (0.05) presented genetic variants in putative LD with the lead MT_5 0.46 0.03 0.31 (0.07) SNPs on SSC7 I2 (low effects: rs330912302 LD = 0.4 MT_90 0.34 0.22 0.30 (0.07) and rs339719658 LD = 0.37) and SSC1 I3 (low effects: VAP_5 0.58 0.002 0.34 (0.07) rs700749617 LD = 0.01, rs710447566 LD = 0.07 or mod- VAP_90 0.55 0.005 0.34 (0.07) erate effect: rs690151450 LD = 6.9 × 10–3), respectively VCL_5 0.61 0.001 0.33 (0.07) (see Additional file 7: Table S5). Although the results on VCL_90 0.55 0.01 0.34 (0.07) SNP calling and LD evaluation should be taken with cau- VSL_5 0.36 0.16 0.31 (0.07) tion, these SNPs in CHD2 and KATNAL2 deserve further VSL_90 0.61 0.001 0.33 (0.07) investigation in larger datasets. ACRO_5 0.5 0.02 0.21 (0.06) ACRO_90 0.21 0.68 0.23 (0.07) R_MT 0.3 0.35 0.13 (0.06) The porcine sperm transcriptome R_VAP 0.18 0.79 0.18 (0.07) The transcriptome profile obtained in this study is very R_VCL 0.28 0.42 0.14 (0.07) similar to that from our previous work [13] and from R_VSL 0.21 0.68 0.21 (0.07) research in other species [72, 73]. Five of the 10 most R_VIAB 0.44 0.05 0.23 (0.07) abundant protein-coding transcripts (PRM1, OAZ3, R_ACRO 0.57 0.003 0.19 (0.07) DNAJB8, TPPP2 and TNP1) have been associated to sperm function via different mechanisms. PRM1 is a Acronym descriptions are in Table 1 protamine that replaces histones in the ultra-compacted SE standard error chromatin of sperm. In a study on bulls, the RNA levels of PRM1 were reduced in low-fertility animals [74], and ungenotyped polymorphisms in transcribed genes that in humans, PRM1/PRM2 sperm ratios differed between might carry important functional implications or may be fertile and infertile men [75]. OAZ3 plays a role in the better genetic markers for that given trait. Should these regulation of polyamine concentration during spermio- genes be involved in related phenotypes and should the genesis and has been linked to sperm function and fer- variants be: (i) in LD with the GWAS lead SNP and (ii) tility in different species such as humans [76] and mice have a predicted effect on protein sequence, these poly- [77]. DNAJB8 is a heat shock binding protein that regu- morphisms could be suggested as potential causal candi- lates the ATPase activity of HSP70, which is a crucial dates. For that purpose, we sought to identify transcribed protein for male fertility and spermatogenesis, and it variants in the QTL regions and assessed their LD with shows reduced RNA levels in infertile men [78]. TPPP2 the lead SNP hit of the QTL. Having said that, RNA-seq has been shown to affect sperm motility, probably by has some particular characteristics (namely splicing that regulating energy production, and fertility in mice [79]. makes read alignment challenging, allele specific expres- TNP1 is a spermatid specific protein that is involved in sion that could miscall a true heterozygous animal with the replacement of histones by protamines in the sperm an incorrect homozygous genotype and RNA editing chromatin [80] and defects in this gene have been shown that post-transcriptionally generates additional variation to cause male infertility [80]. not present in the DNA sequence), which make geno- Our sperm samples contained also a large and varied type calling from RNA-seq a challenging task. Moreover, population of piRNAs (see Additional file 5: Table S3) our analysis was carried on a small number of samples Gòdia et al. Genet Sel Evol (2020) 52:72 Page 15 of 21 [20] and, to a lesser extent, of miRNAs (see Additional eGWAS file 5: Table S3 and Additional file 6: Table S4). piRNAs GWAS hit SNPs may tag causal variants with regula- play an essential role in transposon silencing, are crucial tory functions on gene expression. For this reason, we for proper spermatogenesis [2], and have been anno- also performed a within-trait eGWAS by linking for tated in the sperm of multiple animals species includ- each phenotype, GWAS lead SNPs with genes that have ing humans [81], mice [82], bull [15] and boar [20]. In RNA abundance correlated with the same trait. A robust this study, we focused on the miRNA fraction since, as eGWAS would require a larger sample size. However, we reviewed by Noora Kotaja [83], their involvement in the considered that the analysis was worthwhile as it could maintenance and regulation of spermatogonial stem cell, provide indicative results, which would deserve further meiotic and post-meiotic processes and spermiogenesis investigation in larger populations. We identified three is well documented. Some of the most abundant miRNAs eQTL all with a trans-effect (Table 4) and (see Additional identified in our dataset (see Additional file 6: Table S4) file 10: Table S8). Only one of these regions included present relevant functions for spermatogenesis and genes of interest that were directly associated to sperm embryo development. As a matter of fact, six of the seven quality, i.e. the trans-eQTL on SSC13 for HABN, which miRNA with an average abundance CPM higher than was associated to several genes including actin related 1000 have been linked to sperm function or male fertility. protein 2 (ACTR2) and histidyl-TRNA synthetase (HARS) For example, miR-34c has been proposed to be essential (Table 4) and (see Additional file 10: Table S8). Heid et al. for spermatogenesis, since its absence leads to infertility [96] identified ACTR2 in the sperm head from bulls and in mouse [84], miR-30c that is upregulated in high motile suggested that it has a role in sperm capacitation and bull spermatozoa [15], and miR-191, have abundance lev- acrosome reaction. HARS has also been reported to be els that are significantly correlated with improved human involved in the attachment of histidines to their cor- embryo development [85]. Let-7 has been suggested as a responding tRNA molecules, a fundamental cellular regulator of IGF1 during the differentiation of spermato- process for the translation of mRNA into protein [97]. gonia to primary spermatocytes [86]. Recently, the level Waldron et al. [98] showed that knockout zebrafish for of miR-425-5p in the boar sperm has been linked to far- HARS presented severe defects in high proliferative cells. rowing rate and litter size [87]. Although its role in sperm remains unknown, HARS is overexpressed in sperm of low-fertility bulls [99] and Correlation between genes and miRNAs with semen traits we do not rule out a potential involvement of this gene For mRNA transcripts, the strongest correlation was in spermatogenesis. trans-eQTL hotspots (these trans- between TTC28 and HABN (− 0.71) (see Additional eQTL involving several genes) are of particular interest file 8: Table S6). TTC28 is required for the condensation since their SNPs could have important regulatory roles of spindle microtubules during mitosis and meiosis [88]. and influence gene expression, and thus are more likely Other genes of interest included ABCA3, its RNA levels to contribute to the phenotype. correlating with nine phenotypes (see Additional file 8: Table S6). This gene encodes an ABC transporter that Gene network analysis plays a role in flipin-cholesterol complexes as a mecha- In spite of the considerable number of candidate genes nism to remove cholesterol from the sperm membrane that were identified in our GWAS, many genes might [89]. Although the molecular basis induced by choles- have been missed by this traditional single-trait approach terol efflux from sperm is not well understood, it has due to the lack of an acceptable significant association been reported to be required for sperm capacitation [90]. (FDR > 0.05). After all, sperm quality is a complex poly- Another example is EFHC1 with RNA levels that corre- genic phenotype, which is also influenced by environ- lated with six phenotypes (see Additional file 8: Table S6). mental factors such as husbandry, weather, or testicular Efhc1−/− knockout mice show a reduced flagellar beating pathologies that involve an intricate network of genes frequency [91]. and molecular processes. Moreover, low allelic frequency Several miRNAs of interest including miR-23a, miR- and low LD of the GWAS SNP with the causal variant 27a and miR-122 correlated with seven, eight and eight decrease the power of the GWAS to detect genetic asso- semen quality traits, respectively (see Additional file 9: ciations. For this reason, an alternative strategy to exploit Table S7). miR-23a, is dysregulated in subfertile men [92], GWAS information is to perform an AWM analysis that abundance of miR-27a in spermatozoa is associated with extracts SNPs, which although they have a strong genetic lower progressive motility and normal morphology [93], association but lower than the significance threshold, and expression of miR-122 is associated with abnormal are also associated to a certain number of traits [22]. sperm development [94] and dysregulated in subfertile The association of one SNP to more than one trait pro- men [95]. vides additional robustness to the potential relevance of Gòdia et al. Genet Sel Evol (2020) 52:72 Page 16 of 21 that SNP to semen quality, in our case. This, followed encoded by bromodomain adjacent to zinc finger domain by a PCIT analysis to study gene–gene interactions can 1B (BAZ1B), were also identified. Its closest paralog, provide information on the relevant genes and pathways BAZ1A encodes a member of the chromatin remodeling for some phenotypes and then one can search for SNPs complex [105]. Dowdle et al. [106] showed that Baz1a−/− in these genes or that affect them. Obviously, in paral- mice were infertile because of spermatogenesis defects lel to GWAS, transcriptomics data can contribute addi- tied to changes in chromatin composition. Another TF tional valuable information in the description of these gene of interest was estrogen receptor 1 (ESR1), which genes and pathways. The integration of both sources of was present in the shared network. ESR1 has already information can also be used to improve the accuracy of been associated with pig sperm motility and cytoplas- genomic predictions [35]. We believe that the small sam- matic droplets [107]. Moreover, polymorphisms in ESR1 ple size and the limitation due to measuring only one have been suggested to influence estrogen levels which in ejaculate per boar in our study can be overcome partially turn, affect sperm motility [108]. by the AWM and PCIT approaches. For this reason, we The network comprised several new candidate genes addressed the genetics that underlie boar sperm quality for sperm quality. The trafficking protein particle com- through this integrative systems biology approach. The plex 2 like (TRAPPC2L) gene correlated with 27 miR- genetic co-association and RNA co-abundance interac- NAs including miR-30d, which was the most abundant tions revealed a number of appealing features such as miRNA in our samples (see Additional file 6: Table S4) new candidate genes, TF, TF co-factors, and miRNAs and was found to be dysregulated in oligozoospermic that belong to biological processes and relevant functions men [109] (Fig. 3c). TRAPPC2L belongs to the TRAPPC related to sperm. gene family, with a reported role in ciliogenesis [110]. The TF with the largest number of predicted interac- Interestingly, TRAPPC2L was associated in the final net- tions (129) was encoded by the calcium responsive tran- work with the spermatogenesis and centriole associated scription factor (CARF) gene, its RNA abundance being 1 (SPATC1) gene, which is localized in the neck region in turn, correlated with nine phenotypes (Fig. 3b) and of mouse and human sperm [111]. Disruption of its (see Additional file 8: Table S6). CARF acts as a tran- homolog Spatc1l in mice led to male sterility due to sepa- scriptional activator promoted by calcium influx [100]. ration of sperm heads from tails, thereby advocating for a Since calcium ions are essential in sperm function [101], role in sperm head–tail integrity [112]. The network also we cannot discard the possibility that this TF could be included DNAI2, which correlated with four phenotypes involved in pathways related to sperm maintenance and (see Additional file 8: Table S6). Mutations in DNAI2 function. Some of the CARF predicted target genes from have been associated with ciliary defects and detected in our analysis include interesting candidates such as the males with reduced fertility due to impaired sperm tail la ribonucleoprotein domain family member 4 (LARP4), function [113]. DNAI2 was also associated to boar sperm THADA armadillo repeat containing (THADA) and EF- motility in a previous GWAS [6]. CHD2 is another inter- hand domain containing 1 (EFHC1) genes. LARP4, has esting gene in the network since it was also identified as a been proposed to regulate mRNA stability and trans- candidate gene in our GWAS analysis. This gene included lation of mRNAs [102]. Blagden et al. [102] reported new DNA variants in potential LD with GWAS lead larp-knockout Drosophila mutants in which a consid- SNPs, which would be worth testing in a genetic associa- erable proportion of the spermatocytes had meiotic tion study (Fig. 3d; Table 3). CHD2 was hydroxymethyl- defects. Although the role of THADA remains uncer- ated in human sperm after exposure to bisphenol A, an tain in sperm, Moraru et al. [103] showed that in Dros- epigenetic modifier that causes spermatogenesis defects ophila, THADA modulates the calcium signalling, energy and alters sperm motility [114]. storage and thermogenesis balance. EFHC1 encodes a Of the 94 miRNAs identified in sperm and included in myoclonin1 protein, which has been detected in sperm the final network, 30 interacted with at least 20 genes. flagella in mice testis [104]. Although Efhc1-deficient Some of these 30 miRNAs correlated with sperm traits mice were fertile, mutants presented a reduced ciliary and have also been linked to sperm quality and fertility in (flagellar) beating frequency [91]. previous studies. It is worth noting that miR-16, a miRNA Other TF with a large number of interactions were that was found to be down-regulated in the semen of the SMAD family member 4 (SMAD4) gene (interacting infertile males with sperm abnormalities [115], corre- with 32 genes) and the lysine demethylase 3A (KDM3A) lated with four sperm phenotypes (see Additional file 9: gene (281 gene interactions), both potentially target- Table S7) and potentially interacted with 67 genes (e.g. ing a set of genes that are enriched for cellular macro- ATP9A, found in the shared network and included in the molecular complex assembly processes (see Additional RNA model). Similarly, miR-10b, previously associated file 11: Table S9). TF involved in DNA repair, such as that with human infertile semen samples [116], correlated Gòdia et al. Genet Sel Evol (2020) 52:72 Page 17 of 21 with a motility-related parameter (VCL) and interacted with only 73 SNPs for motility and morphological-related with 32 genes (including the previously discussed TRAP- traits, respectively (Table 5). However, we have used an PC2L that is present in the final network). integrated and informed approach based not only on the GWAS and eGWAS FDR significant associations but also Development of an RNA model and a SNP panel on a robust network built from co-associated SNPs (iden- In this study, we provide a novel and innovative approach tified at suggestive levels but across several phenotypes) to develop an RNA model to estimate the phenotypes as well as gene RNA co-abundance. Moreover, our SNPs based on gene abundances. The model, which includes 10 were chosen to minimize LD between them and thus genes, was predicted to be significant for 10 phenotypes maximize the informativity of the panel. This allowed the and performed best for PDROP and some of the motil- informed inclusion of a large number of SNPs with inde- ity-related traits in our samples (Table 5). The model for pendent marker potential and thus the development of a PDROP reported a highly significant role of the THADA more powerful panel for the prediction of semen quality gene (see Additional file 12: Table S10), which was also in pigs. present in the shared network, and its RNA levels are These results only hold in our Pietrain population with positively correlated with PDROP. THADA regulates a modest sample size and one ejaculate measured per energy metabolism via calcium signalling by binding the boar, thus the validation of the panel will require addi- sarco/ER ­Ca2+ ATPase transporter mechanism [103] tional evaluations in other populations. Nonetheless, which plays an important role in the control of sperm the integrative approach that we propose for ultimately motility acrosome reaction [117]. The CARS2 gene was building a SNP array provides compelling results of its another strong contributor in the model for PDROP and application to any type of complex trait with a genetic was also identified in the shared network (see Additional basis. This opens another avenue to improve traits that file 12: Table S10). CARS2 plays a critical role in protein are influenced by several genes that are of interest for the synthesis but no direct link to spermatogenesis or sperm animal breeding industry. function has been reported. Although SNPs have become the marker of choice for Conclusions the genetic improvement of livestock species, the devel- In summary, our results suggest that the genetic vari- opment of a SNP array for the prediction of boar sperm ants identified in the 12 QTL regions that are mapped quality remains to be done. Here, we propose a SNP to—or near—the CHD2, KATNAL2, SLC14A2, IQCF1 model with 73 SNPs including those identified through and ABCA1 genes, together with other candidate genes the GWAS, eGWAS and gene:gene interaction and based on a systems biology approach including among phenotypic correlation analysis (see Additional file 13: others, LAPR4, THADA, EFHC1, SMADA4, SPATC1 or Table S11). The model could hold promising potential for TRAPPC2L, may modulate sperm quality in pigs. This its application in animal breeding programs. This panel network also includes TF genes such as CARF, with a of 73 SNPs estimated between 5 and 36% of the pheno- large number of potential interactions with target genes typic variance across the 25 traits that were evaluated. that are likely to be key players in shaping the complex These SNPs were better predictors for the phenotypes inheritance of sperm quality traits. We have developed related to sperm abnormalities and motility (Table 5). a SNP panel based on a systems biology approach that Remarkably, when considering only the GWAS lead may be able to explain a larger amount of phenotypic SNPs, the panel explained between 4 and 26% of the phe- variance than that obtained from a stand-alone GWAS. notypic variance, and only for three traits (HABN, NABN The model included GWAS lead SNPs, top eGWAS and TABN) was the model able to predict more than 20% SNPs and SNPs from genes identified in the shared net- of the phenotypic variance. Thus, this systems biology work and could potentially explain more than 30% of approach allowed us to include an additional set of SNPs the phenotypic variance for sperm quality traits such as that increased the predictive potential of the panel. motility and morphology. Although our results are very In a previous study for sperm motility and morphologi- promising for the pig breeding sector, caution should be cal abnormalities using two porcine lines, Marques et al. taken due to the sample size of our study and the lack [6] identified several QTL that cumulatively explained of repeated measures from multiple ejaculates per boar. 10.8% of the genetic variance including 412 and 271 Future work should include the validation of the RNA SNPs for each line. Gao et al. [11] identified 20 and 16 and SNP model in a large number of pigs belonging to QTL that could explain 35.3 and 20.6% of sperm motility different breeds and populations. The implications of and morphological abnormalities traits in Duroc boars, this research are broad, ranging from applications to respectively. Our approach was able to predict 30 to 31% animal breeding strategies to modelling the biology of and 26 to 36% of the variance of the same group of traits infertility in mammals. Gòdia et al. Genet Sel Evol (2020) 52:72 Page 18 of 21 Supplementary Information Additional file 11: Table S9. Gene Ontology analysis of the genes The online version contains supplementary material available at https​://doi. included in the Final Network. GO biological process terms with significant org/10.1186/s1271​1-020-00592​-0. Bonferroni corrected P-values and their associated genes. Additional file 12: Table S10. Parameter estimates for the significant RNA Additional file 1: Figure S1. Outline of the analysis pipeline. It illustrates models. For each of the phenotypes, the model outputs the estimated val‑ the framework of the dataset, analyses and methodologies included in ues for the 10 genes obtained from the GRM regression analysis. The lower the study. the value of Pr >|t|, the higher the involvement of the gene abundance on the total phenotypic variance. Additional file 2: Figure S2. Correlation across boar sperm qual‑ ity traits. Heatmap plot of the correlations among the 25 sperm traits Additional file 13: Table S11. Description of the SNPs included in the measured on 300 boars. CON = concentration; VIAB_5 = viability 5 min; SNP panel. Chromosome, position, SNP ID and analysis from which the VIAB_90 = viability 90 min; ORT = osmotic resistance test; HABN = head SNP was extracted. abnormalities; NABN = neck abnormalities; TABN = tail abnormalities; PDROP = proximal droplets; DDROP = distal droplets; MT_5 = motility 5 min; VAP_5 = average path velocity 5 min; VCL_5 = curvilinear velocity Acknowledgements 5 min; VSL_5 = straight line velocity 5 min; MT_90 = motility 90 min; We thank Betlem Cabrera (CRAG), Dr. Fabiana Quoos Mayer (Instituto de Pes‑ VAP_90 = average path velocity 90 min; VCL_90 = curvilinear veloc‑ quisas Veterinárias Desidério Finamor) and Dr. Martina Rocco (CRAG) for their ity 90 min; VSL_90 = straight line velocity 90 min; ACRO_5 = abnormal laboratory support. We gratefully acknowledge Craig Lewis from Genus PIC acrosomes 5 min; ACRO_90 = abnormal acrosomes 90 min; R_MT = ratio and Sam Balasch from Gepork for contributing the sperm samples. motility; R_VAP = ratio average path velocity; R_VCL = ratio curvilinear velocity; R_VSL = ratio straight line velocity; R_VIAB = ratio viability; Authors’ contributions R_ACRO = ratio acrosomes. Figure S3. SNP based dendrogram for the 25 MG, AS, and AlC conceived and designed the experiments. JR-G carried the semen parameters. Dendrogram of the standardized SNP effects across phenotypic analysis. MG performed sperm purifications and RNA extractions. the 25 sperm traits. AnC designed and carried the RT-qPCR and their analyses. MG analyzed the data with support from AR, RGP and YRC. MG and AlC wrote the manuscript. Additional file 3: Table S1. Effect of external factors on sperm qual‑ All authors discussed the data and read and approved the contents of the ity traits. Effect of farm, age and season per year across the sperm manuscript. All authors read and approved the final manuscript. quality related phenotypes. * = P-value ≤ 0.05; ** = P-value ≤ 0.001; *** = P-value ≤ 0.0001; ns = not significant. Funding Additional file 4: Table S2. Details on the SNPs showing significant This work was supported by the Spanish Ministry of Economy and Competi‑ associations (FDR ≤ 0.05) in the GWAS across autosomal chromosomes tiveness (MINECO) under grant AGL2013-44978-R and grant AGL2017-86946-R and unplaced scaffolds. Chr = chromosome; BP = base pairs (location); and by the CERCA Programme/Generalitat de Catalunya. AGL2017-86946-R Beta = additive effect; MAF = minor allele frequency; FDR = false discovery was also funded by the Spanish State Research Agency (AEI) and the European rate; HABN = head abnormalities; MT_5 = percentage of motile sperma‑ Regional Development Fund (ERDF). We thank the Agency for Management tozoa at 5 min; MT_90 = percentage of motile spermatozoa at 90 min; of University and Research Grants (AGAUR) of the Generalitat de Catalunya NABN = neck abnormalities; PDROP = proximal droplets; R_ACRO = ratio (Grant Numbers 2014 SGR 1528 and 2017 SGR 1060). We also acknowledge abnormal acrosomes. the support of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo Ochoa 2016–2019 (Grant Number SEV-2015-0533) Additional file 5: Table S3. Details of the RNA-seq extraction and map‑ grant awarded to the Centre for Research in Agricultural Genomics (CRAG). ping statistics. Average and Standard Deviation (SD) for the 40 samples MG acknowledges a Ph.D. studentship from MINECO (Grant Number BES- processed, including the amount of RNA obtained and several bioinfor‑ 2014-070560) and a Short-Stay fellowship from MINECO (EEBB-I-18-12860) at matics statistics for total RNA-seq (40 samples) and short RNA-seq (34 AR’s group. samples) datasets. sncRNA = short noncoding RNA; MttRNA = mitochon‑ drial transfer RNA; piRNA = Piwi interacting RNA; snRNA = small nuclear Availability of data and materials RNA; MtrRNA = mitochondrial ribosomal RNA; tRNA = transfer RNA; The datasets generated and/or analysed during the current study are available miRNA = micro RNA; rRNA = ribosomal RNA; snoRNA = small nucleolar at NCBI’s BioProject PRJNA520978. The phenotypic and genotypic datasets RNA. used in the current study are available from the corresponding author on Additional file 6: Table S4. List of protein coding genes and miRNAs reasonable request. identified in sperm. Average and Standard Deviation (SD) for the samples processed. Protein coding and miRNA abundances are expressed in frag‑ Ethics approval and consent to participate ments per kb per million mapped reads (FPKM) and counts per million The ejaculates obtained from pigs were privately owned for non-research pur‑ (CPM), respectively. poses. The owners provided consent for the use of these samples for research. Specialized professionals at the farm obtained all the ejaculates and blood Additional file 7: Table S5. SNPs identified in the RNA-seq data mapping following standard routine monitoring procedures and relevant guidelines. within the GWAS regions. Chr = chromosome. LD = linkage disequilibrium. Genotypic frequency for each of the genotypes. # samples called = num‑ Consent for publication ber of samples with reads in the given SNP position. Not applicable. Additional file 8: Table S6. Correlations between gene abundances and phenotypes. P-values are given when (P-value ≤ 0.05). The correlation Competing interests value is indicated between brackets. ns = not significant. The authors declare that they have no competing interests. Additional file 9: Table S7. Correlations between miRNA abundances Author details and phenotypes. P-values are given when (P-value ≤ 0.05). The correlation 1 Animal Genomics Group, Centre for Research in Agricultural Genomics value is indicated between brackets. ns = not significant. (CRAG) CSIC-IRTA-UAB-UB, Campus UAB, Cerdanyola del Vallès, 08193 Barce‑ Additional file 10: Table S8. Associations identified in the within trait lona, Catalonia, Spain. 2 CSIRO Agriculture and Food, Queensland Bioscience eGWAS. Thirty-nine SNPs showed significant associations (FDR ≤ 0.05) with Precinct, 306 Carmody Rd., St. Lucia, Brisbane, QLD 4067, Australia. 3 Animal semen phenotypes in the GWAS and also displayed significant association Breeding and Genomics, Wageningen University & Research, 6708PB Wage‑ with the abundance of genes which abundance correlated with the same ningen, The Netherlands. 4 Animal Breeding and Genetics Program, Institute phenotype (P-value ≤ 0.05). Chr = chromosome. FDR = False Discovery for Research and Technology in Food and Agriculture (IRTA), Torre Mari‑ Rate; ACRO_5 = Abnormal Acrosomes 5 min; HABN = Head abnormalities. mon, 08140 Caldes de Montbui, Catalonia, Spain. 5 Unit of Animal Science, Gòdia et al. Genet Sel Evol (2020) 52:72 Page 19 of 21 Department of Animal and Food Science, Autonomous University of Barce‑ 20. Ablondi M, Gòdia M, Rodriguez-Gil JE, Sánchez A, Clop A. Characterisa‑ lona, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain. 6 Unit of Animal tion of sperm piRNAs and their correlation with semen quality traits in Reproduction, Department of Animal Medicine and Surgery, Autonomous swine. Anim Genet. 2020. https​://doi.org/10.1111/age.13022​. University of Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, 21. Fraser L, Brym P, Pareek CS, Mogielnicka-Brzozowska M, Paukszto L, Spain. 7 Consejo Superior de Investigaciones Científicas (CSIC), 08003 Barce‑ Jastrzebski JP, et al. Transcriptome analysis of boar spermatozoa with lona, Catalonia, Spain. different freezability using RNA-Seq. Theriogenology. 2020;142:400–13. 22. Fortes MR, Reverter A, Zhang Y, Collis E, Nagaraj SH, Jonsson NN, et al. Received: 2 April 2020 Accepted: 24 November 2020 Association weight matrix for the genetic dissection of puberty in beef cattle. Proc Natl Acad Sci USA. 2010;107:13642–7. 23. Reverter A, Chan EK. Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics. 2008;24:2491–7. References 24. Ramayo-Caldas Y, Ballester M, Fortes MRS, Esteve-Codina A, Castello A, 1. Boe-Hansen GB, Fortes MRS, Satake N. Morphological defects, sperm Noguera JL, et al. From SNP co-association to RNA co-expression: Novel DNA integrity, and protamination of bovine spermatozoa. Andrology. insights into gene networks for intramuscular fatty acid composition in 2018;6:627–33. porcine. BMC Genomics. 2014;15:232. 2. Gòdia M, Swanson G, Krawetz SA. A history of why fathers’ RNA matters. 25. Sanchez MP, Ramayo-Caldas Y, Wolf V, Laithier C, El Jabri M, Michenet Biol Reprod. 2018;99:147–59. A, et al. Sequence-based GWAS, network and pathway analyses reveal 3. Krausz C, Escamilla AR, Chianese C. Genetics of male infertility: from genes co-associated with milk cheese-making properties and milk research to clinic. Reproduction. 2015;150:R159–74. composition in Montbeliarde cows. Genet Sel Evol. 2019;51:34. 4. Gadea J. Sperm factors related to in vitro and in vivo porcine fertility. 26. Reverter A, Fortes MR. Breeding and Genetics Symposium: building Theriogenology. 2005;63:431–44. single nucleotide polymorphism-derived gene regulatory networks: 5. Diniz DB, Lopes MS, Broekhuijse ML, Lopes PS, Harlizius B, Guimaraes SE, Towards functional genomewide association studies. J Anim Sci. et al. A genome-wide association study reveals a novel candidate gene 2013a;91:530–6. for sperm motility in pigs. Anim Reprod Sci. 2014;151:201–7. 27. Tuggle CK, Giuffra E, White SN, Clarke L, Zhou H, Ross PJ, et al. GO- 6. Marques DBD, Bastiaansen JWM, Broekhuijse M, Lopes MS, Knol EF, FAANG meeting: a gathering on functional annotation of animal Harlizius B, et al. Weighted single-step GWAS and gene network analysis genomes. Anim Genet. 2016;47:528–33. reveal new candidate genes for semen traits in pigs. Genet Select Evol. 28. Warr A, Affara N, Aken B, Beiki H, Bickhart DM, Billis K, et al. An improved 2018;50:40. pig reference genome sequence to enable pig genetics and genomics 7. Smital J, Wolf J, De Sousa LL. Estimation of genetic parameters of semen research. Gigascience. 2020;9:giaa051. characteristics and reproductive traits in AI boars. Anim Reprod Sci. 29. King GJ, Macpherson JW. A comparison of two methods for boar 2005;86:119–30. semen collection. J Anim Sci. 1973;36:563–5. 8. Wolf J. Genetic parameters for semen traits in AI boars estimated from 30. Sambrook J, Fritsch EF, Maniatis T. Molecular cloning: A laboratory data on individual ejaculates. Reprod Domest Anim. 2009;44:338–44. manual. Cold Spring Harbor: Cold Spring Harbor Laboratory Press; 1989. 9. Robinson JA, Buhr MM. Impact of genetic selection on management of p. 9.14-9.23. boar replacement. Theriogenology. 2005;63:668–78. 31. R Developmental Core Team. R: a language and environment for statis‑ 10. Zhao X, Zhao K, Ren J, Zhang F, Jiang C, Hong Y, et al. An imputation- tical computing. Vienna: R Foundation for Statistical Computing. http:// based genome-wide association study on traits related to male www.R-proje​ct.org. reproduction in a White Duroc x Erhualian F2 population. Anim Sci J. 32. Taiyun W, Viliam S: R package "corrplot": Visualization of a correlation 2016;87:646–54. matrix; Version 0.84. 2017. https​://githu​b.com/taiyu​n/corrp​lot Accessed 11. Gao N, Chen Y, Liu X, Zhao Y, Zhu L, Liu A, et al. Weighted single-step 30 September 2020. GWAS identified candidate genes associated with semen traits in a 33. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. Duroc boar population. BMC Genomics. 2019;20:797. PLINK: a tool set for whole-genome association and population-based 12. Zhao Y, Gao N, Li X, El-Ashram S, Wang Z, Zhu L, et al. Identifying can‑ linkage analyses. Am J Hum Genet. 2007;81:559–75. didate genes associated with sperm morphology abnormalities using 34. González-Prendes R, Quintanilla R, Cánovas A, Manunza A, Figueiredo weighted single-step GWAS in a Duroc boar population. Theriogenol‑ Cardoso T, Jordana J, et al. Joint QTL mapping and gene expression ogy. 2020;141:9–15. analysis identify positional candidate genes influencing pork quality 13. Gòdia M, Estill M, Castelló A, Balasch S, Rodríguez-Gil JE, Krawetz SA, traits. Sci Rep. 2017;7:39830. et al. A RNA-Seq analysis to describe the boar sperm transcriptome and 35. Ramayo-Caldas Y, Marmol-Sanchez E, Ballester M, Sanchez JP, Gonzalez- its seasonal changes. Front Genet. 2019;10:299. Prendes R, Amills M, et al. Integrating genome-wide co-association and 14. Gòdia M, Mayer FQ, Nafissi J, Castelló A, Rodríguez-Gil JE, Sánchez gene expression to identify putative regulators and predictors of feed A, et al. A technical assessment of the porcine ejaculated sperma‑ efficiency in pigs. Genet Sel Evol. 2019;51:48. tozoa for a sperm-specific RNA-seq analysis. Syst Biol Reprod Med. 36. Pegolo S, Mach N, Ramayo-Caldas Y, Schiavon S, Bittante G, Cecchinato 2018;64:291–303. A. Integration of GWAS, pathway and network analyses reveals novel 15. Capra E, Turri F, Lazzari B, Cremonesi P, Gliozzi TM, Fojadelli I, et al. Small mechanistic insights into the synthesis of milk proteins in dairy cows. RNA sequencing of cryopreserved semen from single bull revealed Sci Rep. 2018;8:566. altered miRNAs and piRNAs expression between High- and Low-motile 37. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome- sperm populations. BMC Genomics. 2017;18:14. wide complex trait analysis. Am J Hum Genet. 2011;88:76–82. 16. Jodar M, Sendler E, Moskovtsev SI, Librach CL, Goodrich R, Swanson S, 38. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A et al. Absence of sperm RNA elements correlates with idiopathic male practical and powerful approach to multiple testing. J R Stat Soc B. infertility. Sci Transl Med. 2015;7:295re6. 1995;57:289–300. 17. Wang X, Yang C, Guo F, Zhang Y, Ju Z, Jiang Q, et al. Integrated analysis 39. Muñoz M, Bozzi R, García-Casco J, Nunez Y, Ribani A, Franci O, et al. of mRNAs and long noncoding RNAs in the semen from Holstein bulls Genomic diversity, linkage disequilibrium and selection signatures in with high and low sperm motility. Sci Rep. 2019;9:2092. European local pig breeds assessed with a high density SNP chip. Sci 18. Curry E, Safranski TJ, Pratt SL. Differential expression of porcine sperm Rep. 2019;9:13546. microRNAs and their association with sperm morphology and motility. 40. Turner SD. qqman: an R package for visualizing GWAS results using Q-Q Theriogenology. 2011;76:1532–9. and manhattan plots. bioRxiv. 2014:005165. 19. Gòdia M, Castelló A, Rocco M, Cabrera B, Rodríguez-Gil JE, Balasch S, 41. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for et al. Identification of circular RNAs in porcine sperm and evaluation of Illumina sequence data. Bioinformatics. 2014;30:2114–20. their relation to sperm motility. Sci Rep. 2020;10:7985. 42. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60. Gòdia et al. Genet Sel Evol (2020) 52:72 Page 20 of 21 43. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg 67. Fang P, Xu W, Li D, Zhao X, Dai J, Wang Z, et al. A novel acrosomal SL. StringTie enables improved reconstruction of a transcriptome from protein, IQCF1, involved in sperm capacitation and the acrosome RNA-seq reads. Nat Biotechnol. 2015;33:290–5. reaction. Andrology. 2015;3:332–44. 44. Reverter A, Barris W, McWilliam S, Byrne KA, Wang YH, Tan SH, et al. 68. Suh KS, Tatunchak TT, Crutchley JM, Edwards LE, Marin KG, Yuspa SH. Validation of alternative methods of data normalization in gene co- Genomic structure and promoter analysis of PKC-delta. Genomics. expression studies. Bioinformatics. 2005;21:1112–20. 2003;82:57–67. 45. Groeneveld E: VCE, a multivariate multimodel REML (co)variance com‑ 69. Hering DM, Olenski K, Rusc A, Kaminski S. Genome-wide association ponent estimation package. In: Proceedings of the 5th world congress study for semen volume and total number of sperm in Holstein- on genetics applied to livestock production: 7–12 August 1994; Guelph; Friesian bulls. Anim Reprod Sci. 2014;151:126–30. 1994. 70. Morales CR, Marat AL, Ni X, Yu Y, Oko R, Smith BT, et al. ATP-binding 46. Martin M. Cutadapt removes adapter sequences from high-throughput cassette transporters ABCA1, ABCA7, and ABCG1 in mouse spermato‑ sequencing reads. EMBnet J. 2011;17:10–2. zoa. Biochem Bioph Res Commun. 2008;376:472–7. 47. Rueda A, Barturen G, Lebron R, Gomez-Martin C, Alganza A, Oliver JL, 71. Liu M, Guan ZL, Shen Q, Lalor P, Fitzgerald U, O’Brien T, et al. Ulk4 is et al. sRNAtoolbox: an integrated collection of small RNA research tools. essential for ciliogenesis and CSF flow. J Neurosci. 2016;36:7589–600. Nucleic Acids Res. 2015;43:W467–73. 72. Sendler E, Johnson GD, Mao S, Goodrich RJ, Diamond MP, Hauser 48. Kozomara A, Griffiths-Jones S. miRBase: integrating microRNA annota‑ R, et al. Stability, delivery and functions of human sperm RNAs at tion and deep-sequencing data. Nucleic Acids Res. 2011;39:D152–7. fertilization. Nucleic Acids Res. 2013;41:4104–17. 49. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The 73. Selvaraju S, Parthipan S, Somashekar L, Kolte AP, Krishnan Binsila sequence alignment/map format and SAMtools. Bioinformatics. B, Arangasamy A, et al. Occurrence and functional significance 2009;25:2078–9. of the transcriptome in bovine (Bos taurus) spermatozoa. Sci Rep. 50. Cingolani P, Platts A, le Wang L, Coon M, Nguyen T, Wang L, et al. A 2017;7:42392. program for annotating and predicting the effects of single nucleotide 74. Feugang JM, Rodriguez-Osorio N, Kaya A, Wang H, Page G, Oster‑ polymorphisms, SnpEff: SNPs in the genome of Drosophila mela- meier GC, et al. Transcriptome analysis of bull spermatozoa: implica‑ nogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80–92. tions for male fertility. Reprod Biomed Online. 2010;21:312–24. 51. Reverter A, Fortes MR. Association weight matrix: a network-based 75. de Mateo S, Gazquez C, Guimera M, Balasch J, Meistrich ML, Ballesca approach towards functional genome-wide association studies. Meth‑ JL, et al. Protamine 2 precursors (Pre-P2), protamine 1 to protamine ods Mol Biol. 2013b;1019:437–47. 2 ratio (P1/P2), and assisted reproduction outcome. Fertil Steril. 52. Berger T, Anderson DL, Penedo MCT. Porcine sperm fertilizing potential 2009;91:715–22. in relationship to sperm functional capacities. Anim Reprod Sci. 76. Jodar M, Kalko S, Castillo J, Ballescà JL, Oliva R. Differential RNAs 1996;44:231–9. in the sperm cells of asthenozoospermic patients. Hum Reprod. 53. Quintero-Moreno A, Rigau T, Rodriguez-Gil JE. Regression analyses and 2012;27:1431–8. motile sperm subpopulation structure study as improving tools in boar 77. Tokuhiro K, Isotani A, Yokota S, Yano Y, Oshio S, Hirose M, et al. OAZ-t/ semen quality analysis. Theriogenology. 2004;61:673–90. OAZ3 is essential for rigid connection of sperm tails to heads in 54. Galili T. dendextend: an R package for visualizing, adjusting and com‑ mouse. PLoS Genet. 2009;5:e1000712. paring trees of hierarchical clustering. Bioinformatics. 2015;31:3718–20. 78. Montjean D, De La Grange P, Gentien D, Rapinat A, Belloc S, Cohen- 55. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Bacrie P, et al. Sperm transcriptome profiling in oligozoospermia. J Cytoscape: a software environment for integrated models of biomo‑ Assist Reprod Genet. 2012;29:3–10. lecular interaction networks. Genome Res. 2003;13:2498–504. 79. Zhu F, Yan P, Zhang J, Cui Y, Zheng M, Cheng Y, et al. Deficiency of 56. Hu H, Miao YR, Jia LH, Yu QY, Zhang Q, Guo AY. AnimalTFDB 3.0: a com‑ TPPP2, a factor linked to oligoasthenozoospermia, causes subfertility prehensive resource for annotation and prediction of animal transcrip‑ in male mice. J Cell Mol Med. 2019;23:2583–94. 57. SAS Institute Inc. 2013: SAS/ACCESS® 9.4 Interface to ADABAS: Refer‑ tion factors. Nucleic Acids Res. 2019;47:D33–8. 80. Schlicker M, Schnulle V, Schneppel L, Vorob’ev VI, Engel W. Distur‑ bances of nuclear condensation in human spermatozoa: search for ence. Cary: SAS Institute Inc; 2013. mutations in the genes for protamine 1, protamine 2 and transition 58. Ramayo-Caldas Y, Renand G, Ballester M, Saintilan R, Rocha D. Multi- protein 1. Hum Reprod. 1994;9:2313–7. breed and multi-trait co-association analysis of meat tenderness and 81. Krawetz SA, Kruger A, Lalancette C, Tagett R, Anton E, Draghici other meat quality traits in three French beef cattle breeds. Genet S, et al. A survey of small RNAs in human sperm. Hum Reprod. Select Evol. 2016;48:37. 2011;26:3401–12. 59. Snelling WM, Cushman RA, Keele JW, Maltecca C, Thomas MG, Fortes 82. Larriba E, Del Mazo J. An integrative piRNA analysis of mouse gametes MR, et al. Breeding and Genetics Symposium: networks and pathways and zygotes reveals new potential origins and gene regulatory roles. to guide genomic selection. J Anim Sci. 2013;91:537–52. Sci Rep. 2018;8:12832. 60. Lopez Rodriguez A, Van Soom A, Arsenakis I, Maes D. Boar management 83. Kotaja N. MicroRNAs and spermatogenesis. Fertil Steril. and semen handling factors affect the quality of boar extended semen. 2014;101:1552–62. Porcine Health Manag. 2017;3:15. 84. Yuan S, Tang C, Zhang Y, Wu J, Bao J, Zheng H, et al. mir-34b/c and 61. Vitavska O, Wieczorek H. Putative role of an SLC45 H(+)/sugar mir-449a/b/c are required for spermatogenesis, but not for the first cotransporter in mammalian spermatozoa. Pflug Arch Eur J Physiol. cleavage division in mice. Biol Open. 2015;4:212–23. 2017;469:1433–42. 85. Xu H, Wang X, Wang ZK, Li JH, Xu ZM, Miao MH, et al. MicroRNA expres‑ 62. Olson GE, Winfrey VP, Nagdas SK. Structural modification of the hamster sion profile analysis in sperm reveals hsa-mir-191 as an auspicious sperm acrosome during posttesticular development in the epididymis. omen of in vitro fertilization. BMC Genomics. 2020;21:165. Microsc Res Tech. 2003;61:46–55. 86. Shen G, Wu R, Liu B, Dong W, Tu Z, Yang J, et al. Upstream and 63. Dunleavy JEM, Okuda H, O’Connor AE, Merriner DJ, O’Donnell L, downstream mechanisms for the promoting effects of IGF-1 on Jamsai D, et al. Katanin-like 2 (KATNAL2) functions in multiple aspects differentiation of spermatogonia to primary spermatocytes. Life Sci. of haploid male germ cell development in the mouse. PLoS Genet. 2014;101:49–55. 2017;13:e1007078. 87. Alvarez-Rodriguez M, Martinez C, Wright D, Barranco I, Roca J, Rodri‑ 64. Li X, Chen G, Yang B. Urea transporter physiology studied in knockout guez-Martinez H. The transcriptome of pig spermatozoa, and its role in mice. Front Physiol. 2012;3:217. fertility. Int J Mol Sci. 2020;21:1572. 65. Wu FJ, Lin TY, Sung LY, Chang WF, Wu PC, Luo CW. BMP8A sustains 88. Izumiyama T, Minoshima S, Yoshida T, Shimizu N. A novel big protein spermatogenesis by activating both SMAD1/5/8 and SMAD2/3 in TPRBK possessing 25 units of TPR motif is essential for the progress of spermatogonia. Sci Signal. 2017;10:eaa1910. mitosis and cytokinesis. Gene. 2012;511:202–17. 66. Luangpraseuth-Prosper A, Lesueur E, Jouneau L, Pailhoux E, Cotinot 89. Mengerink KJ, Vacquier VD. An ATP-binding cassette transporter is C, Mandon-Pepin B. TOPAZ1, a germ cell specific factor, is essential a major glycoprotein of sea urchin sperm membranes. J Biol Chem. for male meiotic progression. Dev Biol. 2015;406:158–71. 2002;277:40729–34. Gòdia et al. Genet Sel Evol (2020) 52:72 Page 21 of 21 90. Visconti PE, Westbrook VA, Chertihin O, Demarco I, Sleight S, Diekman 106. Dowdle JA, Mehta M, Kass EM, Vuong BQ, Inagaki A, Egli D, et al. Mouse AB. Novel signaling pathways involved in sperm acquisition of fertilizing BAZ1A (ACF1) is dispensable for double-strand break repair but is capacity. J Reprod Immunol. 2002;53:133–50. essential for averting improper gene expression during spermatogen‑ 91. Suzuki T, Miyamoto H, Nakahari T, Inoue I, Suemoto T, Jiang B, et al. esis. PLoS Genet. 2013;9:e1003945. Efhc1 deficiency causes spontaneous myoclonus and increased seizure 107. Gunawan A, Kaewmala K, Uddin MJ, Cinar MU, Tesfaye D, Phatsara C, susceptibility. Hum Mol Genet. 2009;18:1099–109. et al. Association study and expression analysis of porcine ESR1 as a 92. Abu-Halima M, Ayesh BM, Hart M, Alles J, Fischer U, Hammadeh M, et al. candidate gene for boar fertility and sperm quality. Anim Reprod Sci. Differential expression of miR-23a/b-3p and its target genes in male 2011;128:11–21. patients with subfertility. Fertil Steril. 2019;112:323–35.e2. 108. Carreau S, Bourguiba S, Lambard S, Galeraud-Denis I, Genissel C, 93. Zhou JH, Zhou QZ, Yang JK, Lyu XM, Bian J, Guo WB, et al. MicroRNA- Levallet J. Reproductive system: aromatase and estrogens. Mol Cell 27a-mediated repression of cysteine-rich secretory protein 2 translation Endocrinol. 2002;193:137–43. in asthenoteratozoospermic patients. Asian J Androl. 2017;19:591–5. 109. Salas-Huetos A, Blanco J, Vidal F, Godo A, Grossmann M, Pons MC, et al. 94. Liu T, Huang Y, Liu J, Zhao Y, Jiang L, Huang Q, et al. MicroRNA-122 influ‑ Spermatozoa from patients with seminal alterations exhibit a differen‑ ences the development of sperm abnormalities from human induced tial micro-ribonucleic acid profile. Fertil Steril. 2015;104:591–601. pluripotent stem cells by regulating TNP2 expression. Stem Cells Dev. 110. Westlake CJ, Baye LM, Nachury MV, Wright KJ, Ervin KE, Phu L, et al. 2013;22:1839–50. Primary cilia membrane assembly is initiated by Rab11 and transport 95. Abu-Halima M, Hammadeh M, Schmitt J, Leidinger P, Keller A, Meese protein particle II (TRAPPII) complex-dependent trafficking of Rabin8 to E, et al. Altered microRNA expression profiles of human spermatozoa the centrosome. Proc Natl Acad Sci USA. 2011;108:2759–64. in patients with different spermatogenic impairments. Fertil Steril. 111. Goto M, O’Brien DA, Eddy EM. Speriolin is a novel human and mouse 2013;99:1249–55.e16. sperm centrosome protein. Hum Reprod. 2010;25:1884–94. 96. Heid HW, Figge U, Winter S, Kuhn C, Zimbelmann R, Franke WW. 112. Kim J, Kwon JT, Jeong J, Kim J, Hong SH, Kim J, et al. SPATC1L Novel actin-related proteins Arp-T1 and Arp-T2 as components of maintains the integrity of the sperm head-tail junction. EMBO Rep. the cytoskeletal calyx of the mammalian sperm head. Exp Cell Res. 2018;19:e45991. 2002;279:177–87. 113. Loges NT, Olbrich H, Fenske L, Mussaffi H, Horvath J, Fliegauf M, et al. 97. Ibba M, Söll D. Aminoacyl-tRNA synthesis. Annu Rev Biochem. DNAI2 mutations cause primary ciliary dyskinesia with defects in the 2000;69:617–50. outer dynein arm. Am J Hum Genet. 2008;83:547–58. 98. Waldron A, Wilcox C, Francklyn C, Ebert A. Knock-down of histidyl-tRNA 114. Zheng H, Zhou X, Li DK, Yang F, Pan H, Li T, et al. Genome-wide synthetase causes cell cycle arrest and apoptosis of neuronal progeni‑ alteration in DNA hydroxymethylation in the sperm from bisphenol tor cells in vivo. Front Cell Dev Biol. 2019;7:67. A-exposed men. PLoS One. 2017;12:e0178535. 99. Aslam MKM, Kumaresan A, Yadav S, Mohanty TK, Datta TK. Comparative 115. Liu T, Cheng W, Gao Y, Wang H, Liu Z. Microarray analysis of microRNA proteomic analysis of high- and low-fertile buffalo bull spermatozoa expression patterns in the semen of infertile men with semen abnor‑ for identification of fertility-associated proteins. Reprod Domest Anim. malities. Mol Med Rep. 2012;6:535–42. 2019;54:786–94. 116. Tian H, Li ZL, Peng D, Bai XG, Liang WB. Expression difference of miR- 100. Tao X, West AE, Chen WG, Corfas G, Greenberg ME. A calcium-respon‑ 10b and miR-135b between the fertile and infertile semen samples (p). sive transcription factor, CaRF, that regulates neuronal activity-depend‑ Forensic Sci Int-Genet Suppl Ser. 2017;6:e257–9. ent expression of BDNF. Neuron. 2002;33:383–95. 117. Harper C, Wootton L, Michelangeli F, Lefièvre L, Barratt C, Publicover S. 101. Publicover S, Harper CV, Barratt C. [Ca2+]i signalling in sperm—making Secretory pathway Ca2+-ATPase (SPCA1) Ca2+ pumps, not SERCAs, the most of what you’ve got. Nat Cell Biol. 2007;9:235–42. regulate complex [Ca2+]i signals in human spermatozoa. J Cell Sci. 102. Blagden SP, Gatt MK, Archambault V, Lada K, Ichihara K, Lilley KS, 2005;118:1673–85. et al. Drosophila Larp associates with poly(A)-binding protein and is required for male fertility and syncytial embryo development. Dev Biol. 2009;334:186–97. Publisher’s Note 103. Moraru A, Cakan-Akdogan G, Strassburger K, Males M, Mueller S, Jabs Springer Nature remains neutral with regard to jurisdictional claims in pub‑ M, et al. THADA regulates the organismal balance between energy stor‑ lished maps and institutional affiliations. age and heat production. Dev Cell. 2017;41:72–81.e6. 104. Suzuki T, Inoue I, Yamagata T, Morita N, Furuichi T, Yamakawa K. Sequen‑ tial expression of Efhc1/myoclonin1 in choroid plexus and ependymal cell cilia. Biochem Biophys Res Commun. 2008;367:226–33. 105. Racki LR, Yang JG, Naber N, Partensky PD, Acevedo A, Purcell TJ, et al. The chromatin remodeller ACF acts as a dimeric motor to space nucle‑ osomes. Nature. 2009;462:1016–21. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions