DNA Research, 2020, 0(0), 1–9 doi: 10.1093/dnares/dsaa019 Advance Access Publication Date: 15 September 2020 Research Article Research Article Whole genome sequencing identifies allelic ratio distortion in sperm involving genes related to Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 spermatogenesis in a swine model Marta Gòdia1, Joaquim Casellas2, Aurora Ruiz-Herrera3,4, Joan E. Rodrı́guez-Gil5, Anna Castelló1,2, Armand Sánchez1,2, and Alex Clop 1,6* 1 Center for Research in Agricultural Genomics (CRAG) CSIC-IRTA-UAB-UB, Cerdanyola del Vallès, Catalonia 08193, Spain, 2Department of Animal and Food Sciences, Autonomous University of Barcelona, Cerdanyola del Vallès, Catalonia 08193, Spain, 3Departament de Biologia Cel.lular, Fisiologia i Immunologia, Autonomous University of Barcelona, Cerdanyola del Vallès, Catalonia 08193, Spain, 4Genome Integrity and Instability Group, Institut de Biotecnologia i Biomedicina (IBB), Autonomous University of Barcelona, Cerdanyola del Vallès, Catalonia 08193, Spain, 5Unit of Animal Reproduction, Department of Animal Medicine and Surgery, Autonomous University of Barcelona, Cerdanyola del Vallès, Catalonia 08193, Spain, and 6Consejo Superior de Investigaciones Cientı́ficas (CSIC), Barcelona, Catalonia 08003, Spain *To whom correspondence should be addressed. Tel. þ34 93 563 66 00. Ext. 3353. Fax. þ34 93 5636601. Email: alex.clop@cragenomica.es Received 11 November 2019; Editorial decision 27 August 2020; Accepted 2 September 2020 Abstract Transmission Ratio Distortion (TRD), the uneven transmission of an allele from a parent to its offspring, can be caused by allelic differences affecting gametogenesis, fertilization or embryo- genesis. However, TRD remains vaguely studied at a genomic scale. We sequenced the diploid and haploid genomes of three boars from leukocytes and spermatozoa at 50x to shed light into the genetic basis of spermatogenesis-caused Allelic Ratio Distortion (ARD). We first developed a Binomial model to identify ARD by simultaneously analysing all three males. This led to the identification of 55 ARD SNPs, most of which were animal-specific. We then evaluated ARD indi- vidually within each pig by a Fisher’s exact test and identified two shared genes (TOP3A and UNC5B) and four shared genomic regions harbouring distinct ARD SNPs in the three boars. The shared genomic regions contained candidate genes with functions related to spermatogenesis including AK7, ARID4B, BDKRB2, GSK3B, NID1, NSMCE1, PALB2, VRK1 and ZC3H13. Using the Fisher’s test, we also identified 378 genes containing variants with protein damaging potential in at least one boar, a high proportion of which, including FAM120B, TDRD15, JAM2 or AOX4 among others, are associated to spermatogenesis. Overall, our results show that sperm is sub- jected to ARD with variants associated to a wide variety of genes involved in different stages of spermatogenesis. Key words: transmission ratio distortion, allelic ratio distortion, sperm, whole genome sequencing, swine C The Author(s) 2020. Published by Oxford University Press on behalf of Kazusa DNA Research Institute. V This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 1 2 M. Gòdia et al. 1. Introduction Raw sequencing reads were filtered to remove adaptors and low- Allelic transmission ratio distortion (TRD) can be defined as the pref- quality reads with Trimmomatic v.0.36.7 Filtered reads were aligned erential transmission of one allele from a heterozygous parent to the to the porcine reference genome (Sscrofa11.1) with the Burrows- offspring and consequently, the departure from the expected ratio of Wheeler Aligner (BWA) ‘mem’ v.0.7.128 and duplicate reads were re- 0.5:0.5 under the Mendelian law of inheritance. Despite their poten- moved using Picard v.2.18.7 (http://broadinstitute.github.io/picard/) tial implications for male fertility, both for human medicine and ani- MarkDuplicates. Variant calling was carried with GATK v.3.8.19 mal breeding, only few studies have explored TRD at a genomic level with base quality score recalibration. SNPs were discovered and fil- in mammals. Some of these studies are based on the genotypes of het- tered with standard hard filtering parameters along with a cluster fil- erozygous parents and their offspring in mouse,1 pig2 and cattle,3 ter (maximum of three variants in a cluster of 50 bp). Indels were and have led to the identification of a few hundreds of loci displaying discarded from further analysis. The resulting single nucleotide poly- TRD. As TRD studies become more powerful with large families, an- morphism (SNP) variants were then filtered for a minimum read imal models such as livestock with large pedigrees are better placed Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 depth of 20 and a maximum of two standard deviations from the av- than humans to carry this research. In swine, Casellas et al.2 scanned erage coverage. The predicted effect of the variants was assessed with the swine genome with 29,373 SNPs in 5 boars and their 352 off- SnpEff v.4.3T.10 Gene ontology analysis was carried with the spring using a Bayesian Factor tool. The authors identified 84 SNPs Cytoscape v.3.6.0 plugin CluGO v.2.5.7.11 P-values were that were heterozygous in at least one boar and displayed significant Bonferroni-corrected (q-value). TRD. As TRD can be caused by defects compromising spermatogen- esis, fertilizing ability or embryo development,4 TRD analysis could 2.1. Assessment of allelic ratio distortion in sperm become an approach complementary to genome-wide association We used two statistical approaches to analyse ARD. In the first ap- studies (GWAS) as it could help mapping genomic regions influenc- proach, we used a Binomial model adapted from Casellas et al.2 to ing reproductive performance that would, otherwise, remain unde- evaluate ARD analysing the three boars simultaneously. Taking the tected. However, the exploration of the potential impact on TRD ith SNP with two alleles (A and B) as example, the likelihood for A caused by allelic ratio distortion (ARD) in the haploid sperm due to defects in spermatogenesis has not been explored thus far. reads from WGS was defined as: A Whole Genome Sequencing (WGS) approach to study TRD in pi ðAÞ ¼ 0:5 þ ai þ bi Hij sire to offspring designs is currently near to unfeasible due to the large number of animals that would need to be sequenced individu- and p(B) ¼ 1 – p(A). Note that ai was the sire-specific ARD parame- ally. The alternative of sequencing pools of gDNA from the offspring ter for the ith SNP, Hij was the proportion of A reads in the diploid is neither a practical option because this would not allow controlling genome of the jth boar and bi was a regression coefficient aiming to for the maternal allelic contribution. This limitation does not exist accommodate technological biases from WGS technology previously when studying ARD in sperm as the sequencing of one ejaculate observed in the diploid genome. For each SNP, the Binomial model allows calculating the allelic ratio in the population of haploid sper- was solved by maximum likelihood. Statistical significance was matozoa, and thus determine the existence of this ARD. In other tested by a standard likelihood ratio test. words, each spermatozoon can be considered as a single individual It was applied to all the variants that were heterozygous in the carrying a haploid genome. blood samples for all three boars. Within each of the six sequenced The aim of this study was to identify variants under ARD in the samples, we used the number of reads carrying each allele to calcu- ejaculate of three boars from an artificial insemination stud. We have late the ratio based on the number of reads for a given allele divided sequenced the genomes of these boars from leukocytes (diploid cells) by the total number of reads in that site. ARD was calculated in and ejaculated spermatozoa (haploid cells) and used the number of sperm (haploid) after correcting its allelic ratio by the ratio in white reads carrying each allele at heterozygous sites as proxies of the alle- blood cells (diploid). The rationale behind this is that the ratio in lic frequency to estimate ARD in sperm. We hypothesize that these blood should be 0.5 and any deviation from this value should be con- SNPs displaying ARD are indicating the presence of loci influencing sidered technical and therefore may also affect sperm (bi 6¼ 0). the efficiency of spermatogenesis and that these may have an impact Moreover, all the heterozygous variants with a ratio below 0.4 or on sire to offspring TRD. above 0.6 in blood were considered to be prone to technical errors and were thus removed from the analysis. We also used the results from this model to compare this ARD with the TRD in swine2. 2. Materials and methods In the second approach, and in order to evaluate ARD indepen- gDNA from blood from three boars of the Pietrain breed from differ- dently within each animal, we first identified the heterozygous site in R ent commercial boar studs were extracted with the MaxwellV RSC each pig in blood (again within the allele ratio 0.4–0.6) and then Whole Blood DNA Kit (Promega Biotech Ibérica SL, Alcobendas, used the Fisher’s exact test to compare the allelic ratio between blood Madrid, Spain) and treated with DNAse-free RNAse (Hoffmann-La and sperm within each animal. Only variants in ARD in sperm above Roche, Basel, Switzerland). Ejaculated sperm from the same animals >0.6 or <0.4 were considered. To correct for multiple testing, a false was obtained by the hand glove method and purified as described by discovery rate (FDR) method was employed. This was applied to: Gòdia et al.5 and gDNA was extracted as in Hammoud et al.6 Blood and sperm samples were collected by specialized professionals. The i. identify coding variants in common genes affected in the three six WGS libraries were prepared with TruSeq DNA PCR-Free Kit boars. This was based on the hypothesis that ARD variants may (Illumina, Inc., San Diego, CA, USA) and sequenced to generate not be shared in the three boars but may affect common genes 150 bp paired end reads in an Illumina’s HiSeq X Ten System. The with similar functional consequences. SNPs located in coding WGS fastq files were deposited in the NCBI Sequence Read Archive regions were extracted with BEDTools intersect v.2.17.0.12 (SRA) under SRA experiment SRX7136525. Coding regions were extracted from the Ensembl (v96) porcine Identification of allelic ratio distortion in sperm 3 annotation. The variant effect on protein sequence was predicted 3.2. Analysis to detect ARD in the SNPs heterozygous with SnpEff v.4.3T.10 in the three boars ii. identify ARD regions shared (less than 1 Mbp apart) in the three Under the hypothesis that ARD variants could be common in a pop- boars which could be indicative of a common affected regulatory ulation, we applied a Binomial model in order to take statistical ad- element. ARD regions were determined by identifying these geno- vantage of analysing the three boars simultaneously. This method mic segments containing at least three ARD SNPs with consecu- allowed the identification of ARD in the polymorphic sites that were tive distances between SNPs below 1 Mbp within each pig. The heterozygous in the three boars, regardless of whether this ARD was ARD regional overlap between the three pigs was evaluated with present in one or more pigs. A total of 302,384 SNPs were heterozy- BEDTools closest and intersect v.2.17.0.12 gous in the three samples. Fifty-five SNPs displayed statistically sig- iii. identify ARD variants with moderate or high functional potential nificant ARD using this Binomial model (Fig. 1). in genes known to be related to spermatogenesis or sperm quality We then evaluated ARD independently within each animal by in each pig regardless of whether they are shared or not in these Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 comparing the allelic events of blood and sperm using the Fisher’s ex- pigs. The hypothesis here was that a large number of different act test. Most of the 55 variants identified with the Binomial model genes and biological pathways may lead to ARD and thus each presented ARD with the Fisher’s exact test in only one pig pig might have its own set of functions altered which may not be (Supplementary Table S3) and 17 did not display ARD in any animal necessarily shared in the three boars. with the Fisher’s exact test. These results suggest that most of the 55 variants are not the With the aim to assess whether our findings of ARD in sperm were ARD causal variant or that ARD is animal-specific. As a matter of the result of stochastic effects, we also employed the Fisher’s exact fact, this approach included only those variants that were heterozy- test, this time identifying first heterozygous sites in balanced allelic gous in the blood of the three pigs, thereby discarding a large propor- ratios (0.4–0.6) (considering sperm as the reference diploid genome) tion of potential candidates. One intergenic SNP presented ARD in and then evaluating the allelic ratio to identify ARD in these SNPs in all the three pigs (SNP_ARD_26; Supplementary Table S3). blood. With the exclusion of 2 SNPs that located within unplaced scaf- folds, 53 of the variants identified with the Binomial model, grouped 2.2. Variant validation by Sanger sequencing into 44 regions containing 1 or more SNPs with consecutive SNP dis- We selected 10 variants for genotype validation using PCR coupled tances below 1 Mbp (Supplementary Figure S2). Thirty-seven, 5 and with Sanger sequencing (Supplementary Table S1). We focused on 2 regions contained 1, 2 and 3 SNPs, respectively (Supplementary ARD variants mapping within genes with known function on sperm Table S3). The previous work from Casellas et al.2 identified 84 biology, spermatogenesis or meiosis. Amplification reactions ranged SNPs in TRD. Of these, seven SNPs could not be liftover into the between 1.5 or 2.5 mM dNTPs, 0.3 mM of each primer, 1.5–2.5 mM coordinates of the Sscrofa11.1 genome assembly or mapped into un- MgCl2, 30 ng of genomic DNA and 0.75 Unit of Amplitaq Gold placed scaffolds. The remaining TRD SNPs were arranged by prox- DNA Polymerase (Thermo Fisher Scientific, Barcelona, Spain). The imity in 63 regions (Supplementary Figure S2; Supplementary Table final volume of the reactions was 15 ml. The thermal profile included S3). Ten out of the 44 ARD regions, containing 12 ARD SNPs, were a denaturation step at 95 C for 10 min, followed by 35 cycles of de- less than 2 Mbp apart from a TRD segment and 1 additional ARD naturation at 95 C for 1 min, annealing at 60 C for 1 min and exten- region marked by 2 SNPs was just 2.08 Mbp from a TRD segment sion at 72 C for 1 min, plus a final extension step at 72 C for 7 min. (Table 1; Supplementary Table S3). These results suggest a possible The specific conditions for each reaction are detailed in shared biological basis and also that a proportion of the TRD may Supplementary Table S1. Amplicons with the expected size were pu- be originated during spermatogenesis. rified with the ExoSAP-IT PCR Clean-up kit (Thermo Fisher These 14 ARD SNPs were less than 100 kbp away from nine cod- Scientific, Barcelona, Spain) and sequenced with the BigDye ing genes (Table 1). Two of the ARD variants (rs1111577152 and Terminator Cycle Sequencing Kit v3.1 (Applied Biosystems, Foster rs1113494508), located 16 bp apart to each other, mapped 54 kbp City, CA, USA) and with the forward or reverse primers listed in downstream from INO80D (Table 1), a INO80 Complex Subunit Supplementary Table S1. Sequencing reactions were electrophoresed member of the chromatin-remodelling complex expressed in develop- in an ABI 3730 DNA analyzer (Applied Biosystems, Foster City, CA, ing spermatocytes, which plays a key role in DNA damage repair as USA). 3. Results and discussion 3.1. WGS, mapping and variant calling In average, 458 M PE reads were obtained per sample (Supplementary Table S2). Up to 99.5% of the reads mapped to the porcine genome (Sscrofa11.1). In average, 14.2% of reads were duplicates and were thus discarded for further analysis. A genome coverage between 46 and 55x was obtained per sample (Supplementary Table S2). The average number of SNPs per sample was 10 M and 6.3 M of these passed quality control filters. From these, in average, 2.8 M SNPs were heterozygous in the blood of each animal (Supplementary Table S2). The reference allele ratio of Figure 1. Manhattan plot of the allelic ratio distortion across the porcine chro- all the heterozygous SNPs displayed very similar distribution in both mosomes. The Binomial model identified 55 significant SNPs in allelic ratio blood and sperm in the three boars (Supplementary Figure S1). distortion. 4 M. Gòdia et al. Table 1. List of ARD regions in close proximity or overlapping to TRD segments ARD or TRD rsID Chr Position Distance between ARD Closest gene distance between and TRD regions gene and ARD SNP TRD 4 111,484,345 1.38 Mbp ARD rs327579254 4 112,863,646 None TRD 7 95,211,457 52 kbp ARD rs342810440 7 95,263,998 SIPA1L1 (90 kbp) TRD 9 23,775,901 0.82 Mbp ARD rs342042877 9 24,599,014 None TRD 11 60,336,131 1.22 Mbp ARD rs339426473 11 61,564,228 None Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 ARD novel 13 200,024,203 1.86 Mbp MORC3 (intronic) TRD 13 201,881,431 ARD rs325570178 14 56,797,388 1.14 Mbp SLC35F3 (intronic) ARD rs340156423 14 57,037,751 ARD rs337352239 14 57,200,494 KCNK1 (60 kbp) TRD 14 58,345,290 ARD rs339246273 15 691,771 0.97 Mbp NEB (intronic) TRD 15 1,657,293 ARD rs1111577152 15 109,256,963 2.08 Mbp INO80D (54 kbp) ARD rs1113494508 15 109,256,979 TRD 15 111,342,012 TRD 16 16,901,264 1.05 Mbp ARD rs325913039 16 17,955,129 GOLPH3 ortholog (51 kbp) ARD rs694882285 17 19,909,268 0.50 Mbp None TRD 17 20,406,228 TRD 18 23,809,561 1.36 Mbp TRD 18 24,134,625 ARD rs788330877 18 25,496,780 FAM3C (70 kbp) In italics, the SNPs identified in the TRD study by Casellas et al. Chr: chromosome. The rsID variant is only provided for the ARD variants identified in our study. it is essential for successful meiosis and spermatogenesis in mice.13 variants detected in TOP3A were novel whilst the variants in Other ARD variants mapped within introns of genes with no UNC5B were already annotated in dbSNP. However, as we do not reported links with spermatogenesis (Table 1). know their allelic frequency in any population, these variants could be thus rare or uncommon. The fact that only two genes harboured 3.3. ARD coding variants in common genes ARD coding variants in the three pigs highlights the complexity and As spermatogenesis includes a set of complex processes including dif- the multi-aetiological nature of ARD. ferent stages such as proliferation and differentiation of spermatogo- nia, meiosis, spermiogenesis and sperm maturation, we hypothesized that ARD could be related to a wide variety of biological pathways. 3.4. Shared ARD regions in the three boars Under this assumption, we expected ARD variants to be rare, or at We also considered the possibility that ARD variants could be rare least not common, and thus not shared between the three pigs. In and only present in one of the three pigs but affect common regula- fact, Huang et al. already suggested that TRD variants tend to be tory regions of relevance in spermatogenesis. We extracted the rare because they are wiped out from the population as one allele is regions in each pig that contained at least three SNPs with consecu- preferentially transmitted to the offspring over the other.4 We there- tive SNP distance below 1 Mbp and then selected those that over- fore sought to identify ARD variants independently in each pig using lapped or were less than 1 Mbp apart in the three pigs. We identified the Fisher’s exact test, which despite being different in the three pigs, four genomic regions in chromosomes 3, 7, 11 and 14 that contained would affect common genes or regulatory elements. a total 55 genes (Table 3; Supplementary Figure S3), several of which The three pigs presented coding variants in ARD in two genes in- play a role at different stages of spermatogenesis (Table 3). volved in spermatogenesis: TOP3A and UNC5B (Table 2). TOP3A Some of the detected genes (VRK1, GSK3B, NID1 PALB2, is a topoisomerase that plays a relevant role in meiotic recombina- ZC3H13 or NSMCE1, among others) are related to early stages of tion, as it has been found to promote the dissolution of double spermatogenesis (spermatogonia proliferation and meiosis). Defects Holliday junctions.14,15 UNC5B is an upstream effector of the in VRK1 have not only been related to spermatogonia loss and infer- Elmo1/Dock180 complex,16 which when disrupted in mice results in tility in male mice18 but also to meiosis in females.19 GSK3B contrib- aberrant seminiferous epithelium, multinucleated giant cells, utes to the induction of meiosis20 and NID1 is related to the uncleared apoptotic germ cells and decreased sperm output.17 In our distribution of meiotic crossover.21 Also related to recombination, survey, TOP3A was affected by three ARD variants. A synonymous PALB222 and ZC3H1323,24 play a role in DNA repair during ho- ARD SNP in TOP3A was shared by two boars and one of these mologous recombination, whereas NSMCE1 is relevant for meiotic boars also presented a missense ARD variant (Table 2). All the chromosome segregation.25,26 Identification of allelic ratio distortion in sperm 5 Table 2. List of ARD variants affecting a common gene in the three boars Sample Chr Start rsID Closest P-value Ratio in Ratio in snpEff Read depth Allele gene blood sperm (blood/sperm) (Ref/Alt) S2 12 60,452,676 novel TOP3A 0.03 0.60 0.35 synonymous 47/40 G/C S1 12 60,465,223 novel TOP3A 0.04 0.58 0.35 missense 48/48 A/G S1 12 60,466,709 novel TOP3A 0.05 0.43 0.65 synonymous 37/46 T/C S3 12 60,466,709 novel TOP3A 0.03 0.54 0.28 synonymous 41/43 T/C S3 14 74,186,268 rs324649834 UNC5B 0.04 0.60 0.38 synonymous 40/55 A/C S2 14 74,199,368 rs339908015 UNC5B 0.02 0.59 0.32 synonymous 41/38 T/C S1 14 74,204,519 rs337527282 UNC5B 0.04 0.60 0.38 synonymous 62/39 C/T The ratios were calculated based on the reference allele. Chr: Chromosome; S: sample. Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 Table 3. List of ARD regions in close vicinity or overlapping in the three samples Chr S1 S2 S3 Genes in the region 3 19,346,924-21,139,827 19,698,060-20,410,463 20,623,640-21,407,898 GTF3C1, NSMCE1, ERN2, & & (5) PALB2, NDUFAB1a, 22,281,850-24,847,672 24,847,578-24,911,195 EARS2, GGA2, COG7, (14) (7) ENSSSCG00000031197, USP31, IGSF6, CDR2a, PDZD9, CRYM, ZP2 7 116,030,235-117,205,191 117,052,038-117,122,913 116,439,959-118,289,816 RF00322, GLRX5, RF02192, (4) (5) (7) RF02193, TCL1B, C14orf132, BDKRB2, BDKRB1, GSK3Ba, AK7, PAPOLAa, VRK1 11 20,080,154-20,761,357 20,283,123-20,824,460 21,085,554-21,441,712 HTR2A, ESD, RUBCNL, LCP1, (5) (4) (3) ENSSSCG00000034648, CPB2, ZC3H13 14 55,926,799-57,746,832 54,287,808-55,764,609 54,474,132-57,200,642 RF00001, RF00019, HEATR1, (9) (6) (11) ERO1B, NID1, LYST, GNG4, RF00026, B3GALNT2, ARID4B, RF00425, TOMM20a, RF00397, IRF2BP2, TARBP1, RF00026, PCNX2 Columns 2, 3 and 4 indicate the ARD genomic intervals for each sample (S1, S2 and S3, respectively). The number of ARD SNPs representing these intervals in each sample is indicated between brackets. Chr: chromosome; S: sample. a Gene name from orthologous genes. Other genes that arise from this study (i.e. AK7, BDKRB2, meiotic stages of spermatogenesis, and thus, ARD variants may af- ZC3H13, ARID4B or HTR2A, among others) have been associated fect non-shared genes. For each boar, we extracted the ARD variants to spermiogenesis, the process in which spermatids mature into sper- with a predicted moderate or high damaging effect on protein se- matozoa in the epithelium of the seminiferous tubules and to the ac- quence thereby potentially altering the protein function. We identi- quisition of sperm motility in the epididymis. This is the case, for fied 408 (131, 129 and 148 for sample S1, S2 and S3, respectively) example, of AK7, which has been linked to spermatogenic failure ARD variants with moderate or high protein damaging effect, none and male infertility probably related to defects on the tail forma- of them was shared between animals and they mapped to 378 genes tion.27 BDKRB2 regulates the AQP9 water channel in the murine ep- (Supplementary Table S4). Six of these variants showed significant ididymis28 and ion transport in the vas deferens of human and pig.29 ARD after correction for multiple testing (FDR  0.05) and affected ZC3H13 is a member of the m6A methyltransferase complex, which five genes (Supplementary Table S4). Four of these (ENSSSCG is involved in the late maturation of spermatids by regulating the ex- 00000034083, ENSSSCG00000030031, ENSSSCG00000033287 pression of key genes.30 ARID4B is involved in Sertoli cell function and ENSSSCG00000039784) are novel genes with unknown func- and is linked to spermatogenic arrest at the stages of meiotic sperma- tion predicted by Ensembl Genebuild after mapping the transcripts tocytes and post-meiotic haploid spermatids,31 whereas HTR2A has and protein sequences from EMBL, GenBank, DDBJ, UniProtKB been associated to sperm count and motility, a property that is ac- and RefSeq databases to the pig genome. The remaining gene, quired by spermatozoa in the epididymis.32 FAM120B, has been shown to play a role in adipogenesis as a trans- activator of PPARG.33 Noteworthy, PPARG has been linked to 3.5. ARD in genes related to spermatogenesis within sperm function and energy metabolism in pigs34 and humans35 and each boar it has also been shown to be expressed in mouse late spermatids and Finally, we also considered the possibility that ARD may originate primary Sertoli cells previously stimulated in vitro with lipopolysac- from a large number of genes and processes throughout the post- charide.36 Thus, it seems plausible that FAM120B modulates sperm 6 M. Gòdia et al. Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 Figure 2. Comparison of the extent of ARD assessed in sperm and assessed in blood. (A) Number of heterozygous SNPs (allelic ratio 0.4–0.6) in each tissue. (B) Number of significant SNPs in ARD (allelic ratios <0.4 or >0.6; *P < 0.05) using the Fisher’s exact test. maturation at the later stages of spermatogenesis through the regula- described genes such as BRCA2,41 EME1,42 GEN1,43 HSF1,44 tion of PPARG. MEI1,45,46 RAD51B,47 RAD9B, a paralog of RAD9A, which is in- The catalogue of 378 genes was enriched for biological functions volved in DNA double-strand break repair during meiosis in mice,48 related to replication fork processing (q-value: 4.4  102), damage MSH249 and PMS2.50 DNA checkpoint (q-value: 4.8  102) and filament cytoskeleton or- The second group of genes associated to ARD variants included ganization (q-value: 2.2  102), which all are relevant processes in- genes related to the formation and maturation of sperm. That is the volved in: (i) the maintenance of genome integrity during meiosis and case of HRB, a gene that is essential for acrosome formation with de- (ii) the formation of the sperm. Of these variants, four had a pre- ficient mice showing meiosis and spermiogenesis defects leading to dicted high impact on TDRD15, JAM2, PCDHGA9 and AOX4. abnormal sperm and infertility.51,52 HIPK4, which is associated to The TDRD family is associated to piwi RNA biology which is essen- abnormal round-headed spermatozoa53 and CFAP100,54 KIF24,55 tial to keep genome stability during spermatogenesis and TDRD15 HAP154 and MARCH10,56 all linked to ciliogenesis. Other genes has been shown to be upregulated in mature versus immature horse within this group were related to sperm maturation in the epididymis testes.37 Little is known however, about the PCDHGA9 protocad- with impact on sperm motility, acrosome formation, mitochondria herin, but protocadherins have been linked to cell adhesion and in homeostasis or capacitation, such as KCNK17,57 PLA2G3,58 addition, PCDHGA9 is mainly expressed in human testes.38 JAM2 PPP3CC,59 SLC26A8,60 CCDC18961 and PINK1.62 has been directly linked to cell adhesion of Sertoli cells to form the Overall, these multiple functions described by the genes linked to blood–testis barrier and to spermatogenesis disruption.39 Finally, ARD variants are a reflex of the complexity of spermatogenesis, sug- AOX4 has been found to be upregulated in germ cells compared gesting that ARD can arise at any moment from meiosis to sperm with Sertoli cells during a synchronized first round of spermatogene- maturation. sis.40 Of note, the alternative alleles of TDRD15 and AOX4, pre- To the best of our knowledge, this pioneer study is the first to eval- dicted to cause a premature stop codon on the protein sequence of uate the potential forces of spermatogenesis that could drive TRD by these genes, were more abundant in the spermatozoa of samples S2 evaluating ARD at the sperm level using WGS. One of the advantages and S3, respectively. This phenomenon could be caused, at least, by of WGS over genotyping platforms is that it allows the interrogation two different scenarios. One possibility is that the reference allele, of practically the whole genome and has thus the potential to identify the one that is present in the reference genome assembly, is not neces- the causal variants. Moreover, WGS allows querying ARD at the sarily the most frequent or the most beneficial allele in a population. sperm level, which would be impossible with genotyping arrays. Alternatively, a detrimental allele could have been hijacked by the al- To validate our results, we first confirmed by PCR followed by lele from another SNP with a stronger influence on spermatogenesis Sanger sequencing, 10 variants identified from the different analyses in close linkage disequilibrium with the flagged SNP. The hijacking we carried and selected for mapping near or for altering the protein scenario is interesting but difficult to test as most likely, ARD alleles sequence of genes with known functions on spermatogenesis with high impact on the sequence of their host proteins would tend (Supplementary Table S1). Two amplicons did not amplify. For the to be rare and thus a large population would have to be screened to other 8 ARD variants, we confirmed all the heterozygous genotypes identify enough animals with all the existing haplotypes and robustly (Supplementary Figure S4). Moreover, one of the amplicons that did measure linkage disequilibrium. not amplify corresponded to SNP_ARD_26, the ARD variant that A careful inspection of the 378 genes associated to ARD variants appeared in the Binomial model and showed ARD in all the three yielded a large proportion of genes whose functions are relevant for pigs. Thus, we could not confirm the existence of the heterozygous spermatogenesis. The first group of genes were involved in the for- state in any of the three samples. We then carried an experiment us- mation and repair of double-strand breaks during meiosis, which re- ing the Fisher’s exact test to assess the allelic ratio in blood when sult in homologous recombination. This group included well- considering only the variants that were heterozygous in sperm with Identification of allelic ratio distortion in sperm 7 allelic ratios between 0.4 and 0.6. In other words, we queried ARD Conflict of interest in blood taking sperm as the diploid reference in which ARD should None declared. not happen. As expected, the number of heterozygous sites was very similar between blood and sperm across the three boars (Fig. 2A). However, the number of SNPs showing allelic ratio deviations (ARD Supplementary data < 0.4 or ARD > 0.6) at P  0.05 in blood versus the sperm reference, Supplementary data are available at DNARES online. was one-third of the number of SNPs in ARD detected in sperm (Fig. 2B). Moreover, while the genes harbouring coding variants in ARD in sperm were enriched for pathways that are relevant for sper- References matogenesis and meiosis (replication fork processing, damage DNA checkpoint and filament cytoskeleton organization), the genes 1. Casellas, J., Gularte, R.J., Farber, C.R., et al. 2012, Genome scans for transmission ratio distortion regions in mice, Genetics, 191, 247–59. Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 encompassing coding variants displaying ARD in blood were associ- 2. Casellas, J., Manunza, A., Mercader, A., Quintanilla, R. and Amills, M. ated to pathways not related to sperm biology, such as microtubule 2014, A flexible Bayesian model for testing for transmission ratio distor- anchoring (q-value: 0.01), maintenance of animal organ identity (q- tion, Genetics, 198, 1357–67. value: 1.6  105), vitamin transport (q-value: 0.01) and cardiac 3. Id-Lahoucine, S., Cánovas, A., Jaton, C., et al. 2019, Implementation of muscle cell contraction (q-value: 0.02). This, together with the fact Bayesian methods to identify SNP and haplotype regions with transmis- that we used stringent criteria to select the SNPs that would be sub- sion ratio distortion across the whole genome: TRDscan v.1.0, J. Dairy jected to the ARD study, suggests that our findings of ARD in sperm Sci., 102, 3175–88. are real and have a biological basis. 4. Huang, L.O., Labbe, A. and Infante-Rivard, C. 2013, Transmission ratio In conclusion, our survey using WGS at 50x depth in three boars distortion: review of concept and implications for genetic association stud- indicates the presence of ARD at the sperm level and shows that ARD ies, Hum. Genet., 132, 245–63. 5. Gòdia, M., Mayer, F.Q., Nafissi, J., et al. 2018, A technical assessment of can arise from multiple stages during spermatogenesis. Forthcoming the porcine ejaculated spermatozoa for a sperm-specific RNA-seq analysis, studies to more deeply investigate ARD in sperm should include a Syst. Biol. Reprod. Med., 64, 291–303. larger number of boars and augmented sequencing depth. This combi- 6. Hammoud, S.S., Nix, D.A., Zhang, H., Purwar, J., Carrell, D.T. and nation would drastically increase the power to identify ARD variants Cairns, B.R. 2009, Distinctive chromatin in human sperm packages genes and clarify the biological basis of spermatogenesis and its consequen- for embryo development, Nature, 460, 473–8. ces on TRD. Moreover, the variants that we identified in this study 7. Bolger, A.M., Lohse, M. and Usadel, B. 2014, Trimmomatic: a flexible should be tested in a larger sire : offspring pedigree to assess their alle- trimmer for Illumina sequence data, Bioinformatics, 30, 2114–20. lic frequency and confirm the TRD effect. Additionally, if sufficiently 8. Li, H. 2013, Aligning sequence reads, clone sequences and assembly con- tigs with BWA-MEM, arXiv, 1303.3997. frequent, these variants should be included in genetic association stud- 9. DePristo, M.A., Banks, E., Poplin, R., et al. 2011, A framework for varia- ies for sperm quality and male fertility to assess their potential implica- tion discovery and genotyping using next-generation DNA sequencing tion on the male’s reproductive ability. data, Nat. Genet., 43, 491–8. 10. Cingolani, P., Platts, A., Wang le, L., et al. 2012, A program for annotat- ing and predicting the effects of single nucleotide polymorphisms, SnpEff: Acknowledgements SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; We thank Sam Balasch (grup Gepork S.A.) for providing the blood and sperm iso-3, Fly (Austin), 6, 80–92. samples. We thank Betlem Cabrera (CRAG) for her lab support. 11. Bindea, G., Mlecnik, B., Hackl, H., et al. 2009, ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway an- notation networks, Bioinformatics, 25, 1091–3. Accession numbers 12. Quinlan, A.R. and Hall, I.M. 2010, BEDTools: a flexible suite of utilities for comparing genomic features, Bioinformatics, 26, 841–2. SRR10441782, SRR10441783, SRR10441781, SRR10441780, 13. Serber, D.W., Runge, J.S., Menon, D.U. and Magnuson, T. 2016, The SRR10441779, SRR10441778. mouse INO80 chromatin-remodeling complex is an essential meiotic fac- tor for spermatogenesis, Biol. Reprod., 94, 8. 14. Dorn, A., Rohrig, S., Papp, K., et al. 2018, The topoisomerase 3 alpha Funding zinc-finger domain T1 of Arabidopsis thaliana is required for targeting the This work was supported by the Spanish Ministry of Economy and enzyme activity to Holliday junction-like DNA repair intermediates, PloS Competitiveness (MINECO) under grant AGL2013-44978-R and grant Genet., 14, e1007674. AGL2017-86946-R and by the CERCA Programme/Generalitat de Catalunya. 15. Martin, C.A., Sarlos, K., Logan, C.V., et al. 2018, Mutations in TOP3A AGL2017-86946-R was also funded by the Spanish State Research Agency cause a bloom syndrome-like disorder, Am. J. Hum. Genet., 103, 456. (AEI) and the European Regional Development Fund (ERDF). We thank the 16. Schaker, K., Bartsch, S., Patry, C., et al. 2015, The bipartite Rac1 Agency for Management of University and Research Grants (AGAUR) of the Guanine nucleotide exchange factor engulfment and cell motility 1/dedica- Generalitat de Catalunya (Grant Numbers 2014 SGR 1528 and 2017 SGR tor of cytokinesis 180 (elmo1/dock180) protects endothelial cells from ap- 01060). We also acknowledge the support of the Spanish Ministry of optosis in blood vessel development, J. Biol. Chem., 290, 6408–18. Economy and Competitivity for the Center of Excellence Severo Ochoa 2016– 17. Elliott, M.R., Zheng, S., Park, D., et al. 2010, Unexpected requirement 2019 (Grant Number SEV-2015-0533) grant awarded to the Centre for for ELMO1 in clearance of apoptotic germ cells in vivo, Nature, 467, Research in Agricultural Genomics (CRAG). MG acknowledges a Ph.D. stu- 333–7. dentship from MINECO (Grant Number BES-2014-070560). We acknowl- 18. Wiebe, M.S., Nichols, R.J., Molitor, T.P., Lindgren, J.K. and Traktman, edge support of the publication fee by the CSIC Open Access Publication P. 2010, Mice deficient in the serine/threonine protein kinase VRK1 are Support Initiative through its Unit of Information Resources for Research infertile due to a progressive loss of spermatogonia, Biol. Reprod., 82, (URICI). 182–93. 8 M. Gòdia et al. 19. Schober, C.S., Aydiner, F., Booth, C.J., Seli, E. and Reinke, V. 2011, The 38. Paul, C. and Robaire, B. 2013, Impaired function of the blood-testis bar- kinase VRK1 is required for normal meiotic progression in mammalian rier during aging is preceded by a decline in cell adhesion proteins and oogenesis, Mech. Dev., 128, 178–90. GTPases, PloS One, 8, e84354. 20. Guo, T.B., Chan, K.C., Hakovirta, H., et al. 2003, Evidence for a role of 39. Evans, E., Hogarth, C., Mitchell, D. and Griswold, M. 2014, Riding the glycogen synthase kinase-3 beta in rodent spermatogenesis, J. Androl., 24, spermatogenic wave: profiling gene expression within neonatal germ and 332–42. sertoli cells during a synchronized initial wave of spermatogenesis in mice, 21. Jeffreys, A.J. and Neumann, R. 2005, Factors influencing recombination Biol. Reprod., 90, 108. frequency and distribution in a human meiotic crossover hotspot, Hum. 40. Sharan, S.K., Pyle, A., Coppola, V., et al. 2004, BRCA2 deficiency in mice Mol. Genet., 14, 2277–87. leads to meiotic impairment and infertility, Development, 131, 131–42. 22. Yan, Z.C., Fan, D.D., Meng, Q.J., et al. 2016, Transcription factor 41. Toledo, M., Sun, X., Brieno-Enriquez, M.A., et al. 2019, A mutation in ZFP38 is essential for meiosis prophase I in male mice, Reproduction, the endonuclease domain of mouse MLH3 reveals novel roles for MutL 152, 431–7. gamma during crossover formation in meiotic prophase I, PloS Genet., 23. Knuckles, P., Lence, T., Haussmann, I.U., et al. 2018, Zc3h13/Flacc is re- 15, e1008177. Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020 quired for adenosine methylation by bridging the mRNA-binding factor 42. Matos, J., Blanco, M.G., Maslen, S., Skehel, J.M. and West, S.C. 2011, Rbm15/Spenito to the m(6)A machinery component Wtap/Fl(2)d, Genes Regulatory control of the resolution of DNA recombination intermediates Dev., 32, 415–29. during meiosis and mitosis, Cell, 147, 158–72. 24. Lin, Z., Hsu, P.J., Xing, X.D., et al. 2017, Mettl3-/Mettl14-mediated 43. Akerfelt, M., Vihervaara, A., Laiho, A., et al. 2010, Heat shock transcrip- mRNA N-6-methyladenosine modulates murine spermatogenesis, Cell tion factor 1 localizes to sex chromatin during meiotic repression, J. Biol. Res., 27, 1216–30. Chem., 285, 34469–76. 25. Taylor, E.M., Copsey, A.C., Hudson, J.J.R., Vidot, S. and Lehmann, 44. Libby, B.J., De La Fuente, R., O’Brien, M.J., et al. 2002, The mouse mei- A.R. 2008, Identification of the proteins, including MAGEG1, that otic mutation mei1 disrupts chromosome synapsis with sexually dimor- make up the human SMC5-6 protein complex, Mol. Cell. Biol., 28, phic consequences for meiotic progression, Dev. Biol., 242, 174–87. 45. Ben Khelifa, M., Ghieh, F., Boudjenah, R., et al. 2018, A MEI1 homozy- 1197–206. gous missense mutation associated with meiotic arrest in a consanguine- 26. Hwang, G., Verver, D.E., Handel, M.A., Hamer, G. and Jordan, P.W. ous family, Hum. Reprod., 33, 1034–7. 2018, Depletion of SMC5/6 sensitizes male germ cells to DNA damage, 46. Da Ines, O., Degroote, F., Amiard, S., Goubely, C., Gallego, M.E. and Mol. Biol. Cell., 29, 3003–16. White, C.I. 2013, Effects of XRCC2 and RAD51B mutations on 27. Lorès, P., Coutton, C., El Khouri, E., et al. 2018, Homozygous missense somatic and meiotic recombination in Arabidopsis thaliana, Plant J. ., 74, mutation L673P in adenylate kinase 7 (AK7) leads to primary male infer- 959–70. tility and multiple morphological anomalies of the flagella but not to pri- 47. Lyndaker, A.M., Vasileva, A., Wolgemuth, D.J., Weiss, R.S. and mary ciliary dyskinesia, Hum. Mol. Genet., 27, 1196–211. Lieberman, H.B. 2013, Clamping down on mammalian meiosis, Cell 28. Belleannée, C., Da Silva, N., Shum, W.W., et al. 2009, Segmental Cycle, 12, 3135–45. expression of the bradykinin type 2 receptor in rat efferent ducts and 48. Kolas, N.K., Svetlanov, A., Lenzi, M.L., et al. 2005, Localization of epididymis and its role in the regulation of aquaporin 9, Biol. Reprod., MMR proteins on meiotic chromosomes in mice indicates distinct func- 80, 134–43. tions during prophase I, J. Cell Biol., 171, 447–58. 29. Pierucci-Alves, F. and Schultz, B.D. 2008, Bradykinin-stimulated cycloox- 49. Baker, S.M., Bronner, C.E., Zhang, L., et al. 1995, Male mice defective in ygenase activity stimulates vas deferens epithelial anion secretion in vitro the DNA mismatch repair gene PMS2 exhibit abnormal chromosome syn- in swine and humans, Biol. Reprod., 79, 501–9. apsis in meiosis, Cell, 82, 309–19. 30. Wu, R.C., Jiang, M., Beaudet, A.L. and Wu, M.Y. 2013, ARID4A and 50. Kang-Decker, N., Mantchev, G.T., Juneja, S.C., McNiven, M.A. and van ARID4B regulate male fertility, a functional link to the AR and RB path- Deursen, J.M. 2001, Lack of acrosome formation in Hrb-deficient mice, ways, Proc. Natl. Acad. Sci. USA, 110, 4616–21. Science, 294, 1531–3. 31. Cortés-Rodriguez, M., Royo, J.L., Reyes-Palomares, A., Lendinez, A.M., 51. Juneja, S.C. and van Deursen, J.M. 2005, A mouse model of familial oli- Ruiz-Galdon, M. and Reyes-Engel, A. 2018, Sperm count and motility are goasthenoteratozoospermia, Hum. Reprod., 20, 881–93. quantitatively affected by functional polymorphisms of HTR2A, MAOA 52. Crapster, J.A., Rack, P.G., Hellmann, Z.J., et al. 2020, HIPK4 is essential and SLC18A, Reprod. Biomed. Online, 36, 560–7. for murine spermiogenesis, Elife, 9, e50209. 32. Li, D.C., Kang, Q.H. and Wang, D.M. 2007, Constitutive coactivator of 53. Firat-Karalar, E.N., Sante, J., Elliott, S. and Stearns, T. 2014, Proteomic peroxisome proliferator-activated receptor (PPAR gamma), a novel coac- analysis of mammalian sperm cells identifies new components of the cen- tivator of PPAR gamma that promotes adipogenesis, Mol. Endocrinol., trosome, J. Cell Sci., 127, 4128–33. 21, 2320–33. 54. Kobayashi, T., Tsang, W.Y., Li, J., Lane, W. and Dynlacht, B.D. 2011, 33. Santoro, M., Guido, C., De Amicis, F., et al. 2013, Sperm metabolism in Centriolar kinesin Kif24 interacts with CP110 to remodel microtubules pigs: a role for peroxisome proliferator-activated receptor gamma (PPAR and regulate ciliogenesis, Cell, 145, 914–25. gamma), J. Exp. Biol., 216, 1085–92. 55. Iyengar, P.V., Hirota, T., Hirose, S. and Nakamura, N. 2011, 34. Aquila, S., Bonofiglio, D., Gentile, M., et al. 2006, Peroxisome Membrane-associated RING-CH 10 (MARCH10 protein) is a proliferator-activated receptor (PPAR)gamma is expressed by human microtubule-associated E3 ubiquitin ligase of the spermatid flagella, J. spermatozoa: its potential role on the sperm physiology, J. Cell. Physiol., Biol. Chem., 286, 39082–90. 209, 977–86. 56. Dube, E., Hermo, L., Chan, P.T. and Cyr, D.G. 2008, Alterations in gene 35. Wang, G., Cheng, S.T., Zhang, S.S., Zhu, Y., Xiao, Y. and Ju, L.G. 2020, expression in the caput epididymides of nonobstructive azoospermic men, LPS impairs steroidogenesis and ROS metabolism and induces PPAR tran- Biol. Reprod., 78, 342–51. scriptional activity to disturb estrogen/androgen receptor expression in 57. Sato, H., Taketomi, Y., Isogai, Y., et al. 2010, Group III secreted phos- testicular cells, Mol. Biol. Rep., 47, 1045–56. pholipase A2 regulates epididymal sperm maturation and fertility in mice, 36. Li, B., He, X., Zhao, Y., et al. 2019, Identification of piRNAs and piRNA J. Clin. Invest., 120, 1400–14. clusters in the testes of the Mongolian horse, Sci. Rep., 9, 5022. 58. Miyata, H., Satouh, Y., Mashiko, D., et al. 2015, Sperm calcineurin inhi- 37. Schmidt, T., Samaras, P., Frejno, M., et al. 2018, ProteomicsDB, Nucleic bition prevents mouse fertility with implications for male contraceptive, Acids Res., 46, D1271–D81. Science, 350, 442–5. Identification of allelic ratio distortion in sperm 9 59. Rode, B., Dirami, T., Bakouh, N., et al. 2012, The testis anion trans- 62. Horibe, A., Eid, N., Ito, Y., Otsuki, Y. and Kondo, Y. 2019, porter TAT1 (SLC26A8) physically and functionally interacts with Ethanol-induced autophagy in sertoli cells is specifically marked at the cystic fibrosis transmembrane conductance regulator channel: a androgen-dependent stages of the spermatogenic cycle: potential mecha- potential role during sperm capacitation, Hum. Mol. Genet., 21, nisms and implications, Int. J. Mol. Sci., 20, 184. 1287–98. 62. Horibe, A., Eid, N., Ito, Y., Otsuki, Y. and Kondo, Y. 2019, 60. Iso-Touru, T., Wurmser, C., Venhoranta, H., et al. 2019, A splice donor Ethanol-induced autophagy in sertoli cells is specifically marked at variant in CCDC189 is associated with asthenospermia in Nordic Red androgen-dependent stages of the spermatogenic cycle: potential mecha- dairy cattle, BMC Genomics, 20, 286. nisms and implications, Int. J. Mol. Sci., 20, 184. Downloaded from https://academic.oup.com/dnaresearch/article/27/5/dsaa019/5906030 by guest on 16 December 2020