PRDM9 Diversity at Fine Geographical Scale Reveals Contrasting Evolutionary Patterns and Functional Constraints in Natural Populations of House Mice Covadonga Vara,1,2 Laia Capilla,1 Luca Ferretti,3 Alice Ledda,4 Rosa A. Sanchez-Guillen,1,5 Sofia I. Gabriel,6 Guillermo Albert-Lizandra,1,2 Beatriu Florit-Sabater,1,2 Judith Bello-Rodrıguez,1,2 Jacint Ventura,7 Jeremy B. Searle,8 Maria L. Mathias,6 and Aurora Ruiz-Herrera*,1,2 1 Genome Integrity and Instability Group, Institut de Biotecnologia i Biomedicina, Universitat Autonoma de Barcelona, Barcelona, Spain 2 Departament de Biologia Cel.lular, Fisiologia i Immunologia, Universitat Autonoma de Barcelona, Barcelona, Spain 3 Oxford Big Data Institute, Li Ka Shing Centre for Health Information and Discovery, Nuffield Department of Medicine, University of Oxford, Oxford, United Kingdom 4 Department for Infectious Disease Epidemiology, Imperial College London, London, United Kingdom 5 Instituto de Ecologıa AC (INECOL), Red de Biologıa Evolutiva, Xalapa, Veracruz, Mexico 6 CESAM – Centre for Environmental and Marine Studies, Departamento de Biologia Animal, Faculdade de Ci^encias da Universidade de Lisboa, Lisbon, Portugal 7 Departament de Biologia Animal, Biologia Vegetal i Ecologia, Universitat Autonoma de Barcelona, Barcelona, Spain 8 Department of Ecology and Evolutionary Biology, Cornell University, Ithaca, NY *Corresponding author: E-mail: aurora.ruizherrera@uab.cat. Associate editor: Belinda Chang Abstract One of the major challenges in evolutionary biology is the identification of the genetic basis of postzygotic reproductive isolation. Given its pivotal role in this process, here we explore the drivers that may account for the evolutionary dynamics of the PRDM9 gene between continental and island systems of chromosomal variation in house mice. Using a data set of nearly 400 wild-caught mice of Robertsonian systems, we identify the extent of PRDM9 diversity in natural house mouse populations, determine the phylogeography of PRDM9 at a local and global scale based on a new measure of pairwise genetic divergence, and analyze selective constraints. We find 57 newly described PRDM9 variants, this diversity being especially high on Madeira Island, a result that is contrary to the expectations of reduced variation for island populations. Our analysis suggest that the PRDM9 allelic variability observed in Madeira mice might be influenced by the presence of distinct chromosomal fusions resulting from a complex pattern of introgression or multiple coloni- zation events onto the island. Importantly, we detect a significant reduction in the proportion of PRDM9 heterozygotes in Article Robertsonian mice, which showed a high degree of similarity in the amino acids responsible for protein–DNA binding. Our results suggest that despite the rapid evolution of PRDM9 and the variability detected in natural populations, functional constraints could facilitate the accumulation of allelic combinations that maintain recombination hotspot symmetry. We anticipate that our study will provide the basis for examining the role of different PRDM9 genetic backgrounds in reproductive isolation in natural populations. Key words: PRDM9, Mus musculus domesticus, Robertsonian fusion, postzygotic reproductive isolation, selection, recombination. Introduction Genes that cause hybrid sterility have been described Understanding the genetic basis of speciation is a long- mostly in Drosophila and include the Odysseus-site standing quest in biology. This entails the investigation of Homeobox (OdsH) gene (Ting et al. 1998), the JYAlpha gene mechanisms responsible for postzygotic reproductive isola- (Masly et al. 2006), Hmr (Barbash et al. 2003), nucleoporin tion, which can be attributable to genic (e.g., speciation genes) Nup96 (Presgraves et al. 2003), and the Overdrive (Ovd) gene or chromosomal (e.g., inversions and fusions) factors (Coyne (Phadnis and Orr 2009). In mice, the PR domain zinc finger 9 and Orr 1998; Orr and Turelli 2001; Faria and Navarro 2010; (Prdm9) gene contributes to hybrid sterility in Mus musculus Farre et al. 2013; Capilla et al. 2014, 2016). Whether both types subspecies (Mus musculus domesticus  Mus musculus mus- of mechanisms can influence speciation processes indepen- culus) (Mihola et al. 2009). The Prdm9 gene is located within dently or in combination remains largely unexplored. the proximal centromeric regions of mouse chromosome 17 ß The Author(s) 2019. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. 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 Open Access 1686 Mol. Biol. Evol. 36(8):1686–1700 doi:10.1093/molbev/msz091 Advance Access publication April 19, 2019 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE and it is implicated in the genomic distribution of recombi- and references therein). This Rb system is characterized by the nation hotspots. It codes for a meiotic-specific histone (H3) presence of seven different metacentric chromosomes methyltransferase with a C-terminal tandem repeat zinc fin- (Rb[3.8], Rb[4.14], Rb[5.15], Rb[6.10], Rb[7.17], Rb[9.11], and ger (ZnF) domain that adds H3K4me3 marks at nucleosomes Rb[12.13]), distributed in nongeographically coincident (stag- close to double-strand breaks (DSBs) in early meiosis; it does gered) clines leading to a progressive reduction in diploid so through its recognition of a species-specific and highly numbers towards the center of the range, about 30-km polymorphic repetitive DNA motif (Mihola et al. 2009; west of the city of Barcelona (Gündüz, Lopez-Fuster, et al. Baudat et al. 2010). The high amino acid variation that char- 2001; Medarde et al. 2012). The Barcelona Rb system is also acterizes PRDM9 in mice and the repetitive nature of the noteworthy for the high levels of chromosomal polymor- DNA motif that it recognizes have suggested the existence phism, which is consistent with the absence of a metacentric of a positive molecular feedback implicated in the generation race (Medarde et al. 2012; Capilla et al. 2014; Sanchez-Guillen of recombination hotspots (Oliver et al. 2009; Jeffreys et al. et al. 2015). On the other hand, mice from the Madeira ar- 2013; Buard et al. 2014; Capilla et al. 2014; Kono et al. 2014; chipelago (including the islands of Madeira and Porto Santo) Schwartz et al. 2014; Baker et al. 2017). That is, the repair of show an extensive chromosomal radiation that includes mice DSBs in early meiosis can introduce mutations in preferred with the standard 2n ¼ 40 (i.e., the island of Porto Santo) and PRDM9 binding motifs, thus influencing the rapid evolution mice with a highly structured chromosomal race organization of novel binding domains (ZnFs) in the PRDM9 protein on the island of Madeira itself (Britton-Davidian et al. 2000). (Boulton et al. 1997; Coop and Myers 2007; Davies et al. Madeira, an island with an extreme topography, has six well- 2016). Such evolutionary turnover has important implications established chromosomal races that have been described for the establishment of reproductive postzygotic barriers within a geographical range of only 742 km2. Diploid numbers since differences in recombination landscapes in hybrids vary from 2n ¼ 22 to 2n ¼ 28, with up to nine metacentric can account for failure in chromosomal synapsis during mei- chromosomes accumulated within a maximum of 1,200 years osis (Qiao et al. 2012; Bhattacharyya et al. 2013; Smagulova (Förster et al. 2009). Importantly, most metacentric popula- et al. 2016). In this scenario, the presence of evolutionarily tions are geographically isolated and do not co-occur with distinct heterozygous combinations of PRDM9 can result in others except for some cases of marginal overlap of the chro- asymmetric DSBs (i.e., differences in genomic distribution and mosomal races. The chromosomal differences between some number of DSBs) between homologous chromosomes, which of the Madeira races are so pronounced that hybrids are at depending on the interaction between different X-linked and very low frequency, presumably due to reduced F1 fitness autosomal loci (Balcova et al. 2016) can result in subfertile and associated with meiotic impairment (Britton-Davidian et al. even sterile phenotypes (Davies et al. 2016). Since interallelic 2000). PRDM9 incompatibilities can result in hybrid sterility due to The contrasting nature of the chromosomal variation in failure in recognition of DNA-binding sites (Flachs et al. 2012), the Madeira and Barcelona Rb systems provides a unique understanding PRDM9 variability in natural populations can opportunity to investigate the mechanisms underlying reveal potential drivers behind the evolutionary dynamics of PRDM9 variation. It permits the examination of PRDM9 evo- this gene. lutionary dynamics between a continental mouse population With regards to the involvement of chromosomal rear- (with a distribution of metacentric chromosomes in a poly- rangements in evolution, house mice represent one of the morphic state; i.e., the Rb Barcelona system) and those of an most extraordinary models of chromosomal variation in island’s chromosomal system (subdivided into entirely or par- mammals. The variation that exists in mice may promote tially reproductively isolated units with metacentric chromo- speciation due to meiotic impairment and infertility associ- somes in a homozygous state; i.e., the Madeira system). Here, ated with chromosomal heterozygosity (Capilla et al. 2014; using a phylogenetic approach on a data set of 400 wild- Pavlova and Searle 2018). The standard karyotype of M. m. caught house mice from the Madeira archipelago (Madeira domesticus consists of 40 all-acrocentric chromosomes. and Porto Santo Islands) and the Barcelona system (conti- However, numerous populations in Western Europe and nental) (fig. 1 ) we 1) identify the extent of PRDM9 diversity in North Africa show high variability in diploid numbers and natural house mouse populations, 2) determine the phylo- karyotypes resulting from Robertsonian (Rb) fusions of non- geography of PRDM9 at a local and global scale using a new homologous acrocentric chromosomes and/or whole arm method for the computation of genetic distances between reciprocal translocations, giving rise to new metacentric chro- complex repeats, and 3) analyze how selection might account mosomes (Pialek et al. 2005). This includes mice from the for the PRDM9 natural diversity observed. In this context, our northeast of the Iberian Peninsula, the so-called “Barcelona Rb study provides the grounds for understanding the complexity system” (Adolph and Klein 1981; Gündüz, Lopez-Fuster, et al. of the mechanisms that can drive evolutionary dynamics of 2001; Medarde et al. 2012) and on the island of Madeira PRDM9 in natural populations. (Britton-Davidian et al. 2000), both Rb systems are character- ized by distinct patterns of chromosomal variation. On the one hand, the Barcelona Rb system extends over an area of Results 5,000 km2 within the provinces of Barcelona, Tarragona, and PRDM9 Diversity at a Local Scale Lleida (Spain) and includes individuals with diploid numbers We successfully obtained polymerase chain reaction (PCR) (2n) ranging from 2n ¼ 27 to 2n ¼ 40 (Medarde et al. 2012, products from the 395 wild-caught mice included in the 1687 Vara et al. . doi:10.1093/molbev/msz091 MBE FIG. 1. Geographical distribution and chromosomal characteristics of the Mus musculus domesticus populations analyzed. (A) Localities sampled in the Barcelona Rb system. Standard populations (blue dots) correspond to mice with 40 acrocentric chromosomes (standard karyotype) whereas Rb populations (red dots) had diploid numbers ranging from 2n ¼ 39 to 2n ¼ 28 (see supplementary table 1, Supplementary Material online, for further details of chromosomes involved in Rb fusions). (B) Localities sampled in the Madeira archipelago. Porto Santo Island is displayed in the inset; it contains standard mice with 40 acrocentric chromosomes (standard karyotype). Six metacentric races have been sampled on the Madeira Island: PEDC, PADC, PLDB, PPOD, PSAN, and PSVI (see supplementary tables 2 and 3, Supplementary Material online, for further details on the chromosomal composition of the chromosomal races). Following previous studies (Britton-Davidian et al. 2000), the chromosomal race PEDC is distributed on both on the northern (PEDC N.) and southern (PEDC S.) coasts of the Madeira Island. study, detecting considerable variation in the number of ZnF Table 1. PRDM9 Variability Found in Wild Mus musculus domesticus repeats, depending on the region sampled (table 1 and sup- including the Present Survey (Barcelona Rb system and Madeira ar- plementary tables 1–4 and figs. 1–4, Supplementary Material chipelago) and Previous Studies in Eurasia (Buard et al. 2014; Kono et al. 2014). online). In the Barcelona Rb system, alleles with 10 ZnF repeats were the most frequently detected (73%), followed Barcelona Madeira Archipelago Eurasiaa by alleles with 12 ZnF repeats (13.2%), 11 ZnF repeats (10%), Madeira Porto Santo 8 ZnF repeats (2%), and 13 ZnF repeats (1.8%). In the case of Madeira however, a higher diversity was detected when con- N 185 199 11 76 PRDM9 alleles 13 53 7 27 sidering the number of ZnF repeats. PRDM9 alleles with 11 Min ZnF 8 10 11 8 ZnF repeats were the most frequently detected (47%), fol- Max ZnF 13 16 12 17 lowed by alleles with 10 ZnF repeats (25%), 13 ZnF repeats NOTE.—N, number of wild mice included in the study; PRDM9 alleles, number of (17%), 12 ZnF repeats (6%), 15 ZnF repeats (2%), and 14 and distinct PRDM9 alleles found based on the amino acid sequence; Min ZnF, mini- 16 ZnF repeats (1% in each). mum number of ZnF repeats; and Max ZnF, maximum number of ZnF repeats. a In order to investigate the source of contrasting variability Data obtained from Kono et al. (2014) and Buard et al. (2014). patterns observed in the Madeira and Barcelona Rb systems, we sequenced the ZnF array in a total of 292 mice. This comparing standard versus Rb populations, with the 10A represented 74% of the sampled mice, a proportion that is allele the most commonly detected variant followed by 11B consistent with previous studies on wild specimens (e.g., and 12B (fig. 3A). The highest variability was found in pop- Kono et al. 2014). This resulted in the identification of 25 ulations characterized by the standard karyotype (2n ¼ 40; distinct ZnFs, four of which (#22, #26, #27, and #29) were Castellfollit del Boix, Olost, and Santa Perpètua de Mogoda) newly identified for M. m. domesticus in this study (supple- where the majority of PRDM9 variants were represented. In mentary table 5, Supplementary Material online). addition, standard populations presented private alleles (8A, In the Barcelona Rb system, we obtained PRDM9 sequen- 10B, 11B, 10C, 12F, and 12G) that were undetected in Rb ces from 132 mice (101 homozygous and 31 heterozygous) populations. In fact, only five alleles (10A, 11C, 12B, 12E, that represented 13 different PRDM9 alleles (figs. 2 and 3). and 12D) were found in Rb mice. Among these, 10A was These were classified according to both the sequence and the most frequently distributed allele, with frequencies rang- number of ZnF repeats: one allele with eight ZnF repeats ing from 23% to 100% in the localities Les Pobles, L’Ametlla de (8A), three alleles with 10 repeats (10A, 10B, and 10C), three Segarra, Sant Sadurnı d’Anoia, and Castelldefels (fig. 3A). alleles with 11 repeats (11B, 11C, and 11D), and six alleles with In the case of the Madeira archipelago, 160 mice were 12 repeats (12B, 12C, 12D, 12E, 12F, and 12G) (fig. 2). Eight of successfully sequenced (110 homozygous and 50 heterozy- the 13 different PRDM9 alleles found in the Barcelona Rb gous), 152 of which were from the Madeira Island Rb system system (8A, 10C, 11C, 11D, 12D, 12E, 12F, and 12G) were and 8 mice from Porto Santo Island (fig. 4). This represented newly identified in this study. The diversity of PRDM9 alleles 54 different alleles: six alleles with 10 ZnF repeats, 20 alleles was widely variable among localities, particularly when with 11 ZnF repeats, 7 alleles with 12 ZnF repeats, 1 allele with 1688 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE A B FIG. 2. Natural PRDM9 allelic diversity found in the Barcelona Rb system and Madeira archipelago Mus musculus domesticus populations. (A) Representation of the mouse PRDM9Dom2 protein (reference genome, C57/BL6 strain). It consists of four domains: KRAB-like, SSXRD, PR/SET, and ZnF array. The underlined ZnFs represent repeats from position ZnF3 to ZnF6 that bind to the DNA in the mouse (Baker et al. 2015; Paigen and Petkov 2018). (B) Representation of the ZnF alignments for all PRDM9 alleles found in the present study. Each ZnF is color coded based on amino acid sequence affinity found at the most variable sites (1, þ3, þ6) responsible for DNA binding. Additional information on ZnF amino acid sequences is provided in supplementary table 5, Supplementary Material online. Purple boxes encompass the more variable ZnF among individuals. Red asterisks indicate previously described PRDM9 alleles (Buard et al. 2014; Capilla et al. 2014; Kono et al. 2014). Geographical labels: (b) PRDM9 alleles found only in the Barcelona Rb system, (b/m) PRDM9 alleles found in both the Barcelona Rb system and the Madeira archipelago. PRDM9 alleles without geographical label correspond to alleles found only in the Madeira archipelago. 14 ZnF repeats, 6 alleles with 15 ZnF repeats, and 3 alleles with Moreover, some alleles were detected in all Rb races (supple- 16 ZnF repeats (figs. 2 and 4). Only five of them (10A, 11B, 11S, mentary table 2, Supplementary Material online) but absent 11V, and 12C) were previously described in Eurasian popula- in standard populations (40 all-acrocentric chromosomes) tions and the Barcelona system (Buard et al. 2014; Capilla et al. from Porto Santo (i.e., 10A and 13C). We found the greatest 2014; Kono et al. 2014). The allele 10A corresponds to the allele heterogeneity among the chromosomal races PADC alleles 48 and Ce3 described by Kono et al. (2014) and Buard (2n ¼ 24–28) and PEDC (2n ¼ 23–26), occurring in the west- et al. (2014) in M. m. domesticus and Mus musculus castaneus, ern part of the island. In contrast, populations situated in respectively. The allele 11B is homologous to the previously eastern regions (PLDB, PPOD, PSAN, and PSVI) had a more reported alleles 54, 55, and Db1 in M. m. domesticus (Buard homogeneous allele distribution (fig. 4A). Regarding standard et al. 2014; Kono et al. 2014), whereas 11S is the PRDM9 allele populations (2n ¼ 40) on the island of Porto Santo, the allele present on the t-haplotype (Kono et al. 2014). The 11V and variability was lower than other populations from eastern 12C allele corresponded to allele 3 and 16, respectively (Buard Madeira. Among the alleles reported on Porto Santo, three et al. 2014). of them (11X, 11V, and 11W) were private alleles that Among the 49 alleles newly identified in the Madeira ar- were undetected in Rb populations and two more (12C chipelago, the most commonly represented on Madeira and 11F) were present in very low frequencies in the Island was 11K, followed by 10A, 10E, 11S, and 13C (fig. 4A). Madeira Island. 1689 Vara et al. . doi:10.1093/molbev/msz091 MBE A B C D FIG. 3. Allelic frequencies, population structure and phylogeny of the ZnFs found in the Barcelona Rb system. (A) Representation of the PRDM9 allelic frequencies in each locality (see table 1 and supplementary table 1, Supplementary Material online, for further details on the chromosomal composition of each locality). Legend—N, number of specimens sequenced per locality. (B) PCA in a subset of 50 mice from 5 localities of the Barcelona Rb system (CAS, Castelldefels; OLO, Olost; BOI, Castellfollit del Boix; ANO, Sant Sadurnı d’Anoia; and MOG, Sta. Perpètua de Mogoda). (C) Plots showing the proportion of inferred ancestry for K ¼ 3 to K ¼ 5 in a subset of 50 mice from 5 localities of the Barcelona system (CAS, Castelldefels; OLO, Olost; BOI, Castellfollit del Boix; ANO, Sant Sadurnı d’Anoia; and MOG, Sta. Perpètua de Mogoda). (D) Phylogenetic tree including all mice sequenced from the Barcelona Rb system (see text for further details). PRDM9 alleles are color coded following the pattern used in panel (A). For each individual, both the karyotype (white: standard; red: Rb fusions) and the locality (see inserted color legend) are indicated. Genetic Diversity and Signal of Positive Selection on Understanding the evolutionary consequences of this ZnF the ZnF Array Domain array variability is important since polymorphisms in both Once the PRDM9 diversity was characterized at a local scale, number and sequence of the arrays might affect DNA- we analyzed genetic diversity and the signal of positive selec- binding specificity resulting in changes in the recombination tion on the ZnF array domain at a global scale, including pre- landscape of each individual (Davies et al. 2016). Thus, we viously described Eurasian populations (Buard et al. 2014; Kono analyzed whether selection played a role in shaping this var- et al. 2014). When comparing the overall genetic (nucleotide) iability by inferring dN/dS estimates averaged across all pair- variability of the ZnF array domain found in all three regions wise PRDM9 sequence alignments considering either included in the study (Barcelona, Madeira, and Eurasia), we hypervariable codons only (1, þ3, þ6) or excluding them observed that the genetic diversity in Madeira was greater than (see Material and Methods). We observed that the amino in the Barcelona Rb system and comparable to the diversity acids that recognize the specific DNA motif at positions across the whole of Eurasia (table 2). Moreover, the genetic 1, þ3, and þ6 are under strong positive selection, as shown differentiation between Madeira and the Barcelona Rb system by their extremely high dN/dS values (table 3). Position þ6 (FST ¼ 0.0295, P < 0.0001) or even between Madeira and showed less variability involving two of the residues described Eurasia was low (FST ¼ 0.0084, P < 0.0001). (Q and K) compared with positions 1 (A, E, Q, T, and V 1690 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE A B FIG. 4. Allelic frequencies and phylogeny of the ZnFs found in the Madeira archipelago. (A) Representation of the PRDM9 allelic frequencies found in each metacentric race from the Madeira Island and Porto Santo Island (see supplementary tables 2 and 3, Supplementary Material online, for further details on chromosomal composition of metacentric races). Number of specimens sequenced per locality on Madeira Island: Achadas da Cruz (N ¼ 7), Arco da Calheta (N ¼ 11), Ch~ao da Ribeira (N ¼ 16), Estreito da Calheta (N ¼ 3), Faial (N ¼ 2), Faj~a da Ovelha (N ¼ 2), Levada Grande (N ¼ 2), Lombada dos Cedros (N ¼ 1), Lombada dos Marinheiros (N ¼ 4), Lombo da Velha (N ¼ 3), Lombo das Laranjeiras (N ¼ 3), Lombo do Doutor (N ¼ 3), Lugar da Raposeira (N ¼ 2), Madalena do Mar (N ¼ 2), Maloeira (N ¼ 18), Moledos (N ¼ 2), Ponta Delgada (N ¼ 8), Ponta do Pargo (N ¼ 18), Ponta do Sol (N ¼ 3), Porto da Cruz (N ¼ 1), Porto Moniz (N ¼ 3), Prazeres (N ¼ 2), Ribeira da Janela (N ¼ 6), Ribeira da Laje (N ¼ 3), Ribeira da Vaca (N ¼ 3), Ribeira Funda (N ¼ 3), Santa (N ¼ 1), Santana (N ¼ 4), S~ao Vicente (N ¼ 3), Seixal (N ¼ 1), Sıtio da Faj~a (N ¼ 3), Socorro (N ¼ 7), Solar da Maloeira (N ¼ 1), Solar dos Prazeres (N ¼ 1), and Porto Santo Island (N ¼ 8). (B) Phylogenetic tree including all mice sequenced from the Madeira archipelago (see text for further details). PRDM9 alleles are color coded following the pattern used in panel (A). For each individual, both the diploid number (2n) and the geographical distribution are indicated (see inserted color legend). substitutions) and þ3 (D, H, N, V, and S substitutions) (fig. 2). Table 2. Genetic Diversity Estimates of the ZnF Array Found in the Interestingly, the rest of the PRDM9 sequence (the so-called Studied Areas (Madeira archipelago, Barcelona Rb system, and Eurasia). nonhypervariable sites) was similarly subjected to positive selection in both the Madeira the Barcelona Rb systems, Madeira Barcelona Eurasia whereas it was under weak negative selection across Eurasia Average genetic distance 5.406 2.307 6.123 (table 3). Moreover, we detected differences in dN/dS values Average pairwise nucleotide 0.00456 0.002856 0.00416 between areas of study, with the effects of selection being diversity excluding hypervariable pronounced in the Madeira archipelago, intermediate in the sites (21, 13, 16) Barcelona Rb system and weaker across Eurasia (table 3). To Average pairwise nucleotide diversity 0.02351 0.018626 0.02694 (nucleotide differences per base) test whether these contrasting patterns in positive selection across regions was explained by the presence of Rb fusions, we calculated dN/dS values considering mice with and without Rb fusions in the Madeira and Barcelona Rb systems. (Billings et al. 2013; Baker et al. 2015; Paigen and Petkov Remarkably, individuals with Rb fusions showed higher dN/ 2018). In the case of the house mouse, this hypothesis dS values than standard ones, both in Madeira and the implies that the ZnF repeats located at positions ZnF3 Barcelona Rb system (table 3). This was consistent for both to ZnF6 within the array would make primary contact hypervariable sites and the rest of the ZnF array domain with the DNA strand before recruiting the DSBs- (table 3). induction machinery for the formation of recombination hotspots (review by Paigen and Petkov [2018]). Since the presence of PRDM9 heterozygous allelic combinations Differences in Heterozygous PRDM9 Allele might result in subfertile phenotypes caused by asymmet- Combinations ric DSBs between homologous chromosomes, we sought Even though the PRDM9 protein can use all its ZnFs to to understand the functional constraints that might affect bind to the DNA in a sequence-specific manner, it has PRDM9 diversity in Barcelona and Madeira populations by been suggested that only a subset of ZnFs within the array analyzing the degree of conservation of positions ZnF3– contributes to the binding of the specific DNA motif ZnF6 within the array in PRDM9 heterozygous mice. 1691 Vara et al. . doi:10.1093/molbev/msz091 MBE (0.0001–0.0019) NOTE.—Ninety-five percent confidence intervals for all estimates are indicated in parentheses and italics. hs, hypervariable sites (1, þ3, þ6); nhs, nonhypervariable sites; Rb, individuals with Rb fusions; St, individuals with standard karyotype (i.e., In the case of the Madeira system, we detected a clear (0.0006–0.001) 4.79 (2.82–8.13) (0.003–0.006) 197 (83–1110) 196 (136–321) (0.15–0.19) deficit of heterozygotes with respect to Hardy–Weinberg 0.00086 0.00413 0.00086 0.169 equilibrium within Rb races, with PRDM9 heterozygotes re- All duced to about half of the expected (supplementary table 6, Supplementary Material online). These deviations from Hardy–Weinberg equilibrium were not found in standard 1.25 (0.75–1.97) (0.001–0.008) (0.002–0.005) (0.002–0.004) 53.6 (24.8–236) 293 (100–8920) 75.3 (51.5–115) (0.19–0.24) mice in Porto Santo, suggesting a relationship between the 0.00407 0.00362 0.00289 Eurasia presence of Rb fusions and the selection against PRDM9 0.218 heterozygotes on Madeira. A contrasting pattern was found in the Barcelona Rb system where both standard and Rb populations show a reduction in the number of heterozygote (1.1–05–0.0008) (0.001–0.004) individuals for PRDM9 (supplementary table 6, (0.07–0.12) 7.3 (2.2–182) St versus Rb ‘ (1–1) 0.00237 0.00032 0 (0–0) Supplementary Material online). 0.0951 As a way to explain this reduction in heterozygosity, we Table 3. Average dN, dS, and dN/dS Values between Aligned ZnF Arrays in the Studied Areas (Madeira archipelago, Barcelona Rb system, and Eurasia). observed that not all heterozygous PRDM9 combinations were equally represented among populations (fig. 5 and sup- (0.0002–0.003) plementary table 7 and fig. 4, Supplementary Material online). (0.03–0.07) ‘ (1–1) ‘ (1–1) ‘ (1–1) The most common PRDM9 allele combinations in the 0.00144 0 (0–0) 0 (0–0) 0.0516 Rb Barcelona Rb system were 10A/12B, 10B/11B, and 10A/11C, whereas in the Madeira archipelago, most frequent alleles Barcelona were 11K/13C, 10A/11K, and 10A/13C. Although in the Madeira archipelago no heterozygous combinations were (9.1–05–0.001) 142 (74.7–367) 313 (113–3180) 187 (78.5–1290) 5.16 (2.25–14.8) 7.95 (2.45–83.7) 4.34 (1.51–36.7) (0.002–0.005) (0.10–0.15) ‘ (1–1) shared between standard and Rb mice, only one allele com- 0.00298 0.00067 0 (0–0) 0.128 bination (10A/12B) was present in both standard and Rb St mice in the Barcelona Rb system (fig. 5A). In both Rb systems (Barcelona and Madeira), differences in population allelic fre- quency were significant between standard and Rb mice (2.4–05–0.0007) (0.001–0.004) (0.07–0.11) (Pearson test, P  0.05). When analyzing the conservation ‘ (1–1) 0.00226 0.00028 0 (0–0) 0.089 in amino acid sequence of the ZnF repeats located at posi- All tions ZnF3–ZnF6 along the ZnF array in the different hetero- zygous combinations, we detected higher sequence conservation in Rb mice when compared with standard (0.0006–0.003) (0.005–0.011) mice in both the Madeira and Barcelona Rb systems (fig. 5). (0.19–0.26) St versus Rb ‘ (1–1) 0.00817 0.00158 In the case of the Barcelona Rb system, an average conserva- 0 (0–0) 0.225 tion of 94.4% was found in Rb mice versus 75% in standard mice, this difference being statistically significant (Wilcoxon test, P ¼ 0.0075). These values were lower in the Madeiran (0.0002–0.0009) archipelago, where the average sequence conservation of (0.003–0.006) 9.2 (4.65–20.7) 330 (189–708) (0.15–0.20) ‘ (1–1) repeats from ZnF3–ZnF6 in Rb mice was 70.5% versus 0.00488 0.00053 0 (0–0) 0.175 44.4% in mice with no Rb fusions, though the differences in Rb conservation were statistically significant (Wilcoxon test, Madeira P ¼ 0.0252) (fig. 5 and supplementary fig. 4, Supplementary Material online). dN(nhs)/dS(nhs) 8.12 (3.99–17.9) 2.15 (0.94–7.62) (0.006–0.013) (0.001–0.008) 278 (159–566) 53.7 (27.7–179) (0.19–0.28) ‘ (1–1) 0.00929 0.00433 0 (0–0) Phylogeography of PRDM9 0.232 St In order to investigate whether the observed variability of PRDM9 was reflective of population structure, we studied absence of Rb fusions); and 1, infinite values. the phylogeography of the different alleles detected in the (0.0003–0.001) (0.004–0.007) Rb systems individually (Barcelona and Madeira) and in rela- (0.15–0.21) ‘ (1–1) 0 (0–0) 0.0053 0.0006 tion to previously described alleles in Eurasian populations 0.181 All (Buard et al. 2014; Capilla et al. 2014; Kono et al. 2014). Although previous allozyme and microsatellite analysis (Britton-Davidian et al. 2007; Förster et al. 2013) provided information regarding the origin of the chromosomal radia- dN(hs)/dS(nhs) dN(hs)/dS(hs) tion of Madeiran mice, the population structure of the Barcelona system is not entirely known at this stage. Thus, dN(nhs) dS(nhs) dN(hs) dS(hs) we first determined genetic diversity and population differ- entiation in this system based on the single nucleotide 1692 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE A B Barcelona Madeira Rb system archipelago Rb St Rb St ZnF3 ZnF4 ZnF5 ZnF6 13J/15C 13B/15E 10A 11K/12K 11B/12C 12E 10A/12G 10A/12E 10A/12B 100% conservation ZnF3-ZnF6 10A/11C 8A/12F 13C/15D Population 11K/15F frequency (%) Prediction: Hotspot symmetry in homologous chrs 11C/12D PRDM9 heterozygous 11B/12J 10A/11P 10A/11K combinations 13B/15A 10A/12D 10E/13C ZnF3-ZnF6 10A/13C Aa sequence C conservation (%) 10A/12L 10A/12C ZnF3 ZnF4 ZnF5 ZnF6 10A/11B 13C/16A 10B 12K/13C 11K/13C 11B 11B/12H 11S/13C 11K/13H 33.3% conservation ZnF3-ZnF6 10A/13G 11V/12C 11U/15B 10G/11T 10B/11B Prediction: Hotspot asymmetry in homologous chrs 10A/11S 8A/12G 11W/12C 11S/12C FIG. 5. Conservation in amino acid sequence of the ZnF repeats located in positions 3–6 (ZnF3–ZnF6) along the ZnF array in PRDM9 heterozygous combinations found in the study. (A) Bubble chart representing population frequency (size of the bubble) and the percentage of amino acid (Aa) sequence conservation of the highly variable positions (1, þ3, and þ6) of the initial repeats along the array (from ZnF3 to ZnF6) (purple intensity). Data are shown for each heterozygous combination found in Standard and Rb mice in the Barcelona Rb system and the Madeira archipelago (see supplementary table 7, Supplementary Material online, for further details). (B, C) Examples of PRDM9 alleles in heterozygous combination and the predicted distribution of recombination hotspots. Yellow starts represent the location of recombination hotspots across homologous chromosomes (red and blue lines). (B) The ZnF repeats located in positions 3–6 (ZnF3–ZnF6) in the 10A/12E combination present 100% of Aa sequence conservation. This can result in symmetric distribution of recombination hotspot in homologous chromosomes. (C) The ZnF repeats located in positions 3–6 (ZnF3–ZnF6) in the 10B/11B combination present 33.3% of Aa sequence conservation. This can result in asymmetric distribution of recombination hotspot in homologous chromosomes. polymorphism (SNP) genotyping of a subset of 50 mice using differentiation, being more pronounced in Castelldefels the Mega Mouse Universal Genotyping Array. We applied a (fig. 3B), a population with the lowest diploid number maximum likelihood method to infer the genetic ancestry of (2n ¼ 28–31) and with the presence of a high number of each individual, where the individuals are assumed to have heterozygous Rb fusions (supplementary table 1, originated from K ancestral clusters. The results of the prin- Supplementary Material online). cipal component analysis (PCA) and the plots for K ¼ 3 to Likewise, estimation of population structure revealed het- K ¼ 5 are shown in figure 3B and C. We observed clustering erogeneity of genetic clusters according to populations for the two Rb populations (Castelldefels and Sant Sadurnı (fig. 3C). The genetic structure based on the ADMIXTURE d’Anoia) when compared with the standard populations analysis (Alexander et al. 2009) revealed that mice from (Olost, Castellfollit del Boix, and Santa Perpètua de Castelldefels (one of the populations with the highest number Mogoda). In fact, the PCA revealed that Olost and of Rb fusions detected) showed high levels of ancestry Castellfollit del Boix populations, both located at the northern (0.99%), forming a cohesive cluster (fig. 3C). An exception area of the Rb system distribution, clustered together. to this pattern was a single individual (specimen 955), which Compared with these, both Rb populations (Castelldefels shared ancestry with mice from Sant Sadurnı d’Anoia. This and Sant Sadurnı d’Anoia) showed progressive genetic individual was characterized by the highest diploid number in 1693 Vara et al. . doi:10.1093/molbev/msz091 MBE Castelldefels (2n ¼ 32) and by the absence of Rb fusion 3.18, Discussion characteristic of the remaining mice from this population Distinctive Phylogeographic Patterns of ZnF (supplementary table 1, Supplementary Material online). Within the Sant Sadurnı d’Anoia population (2n ranging Variability and the Evolution of PRDM9 in the House from 35 to 39), mice were distributed into two different Mouse across Eurasia groups, with three individuals sharing ancestry with the stan- Our study includes a wide survey of nearly 400 wild-caught dard populations Olost and Castellfollit, both populations mice representing two highly distinct house mouse Rb sys- highly homogenous according to the PCA. The Santa tems: the Madeira archipelago (Madeira and Porto Santo Perpètua Mogoda population was very unusual, despite its Islands) and the Barcelona system (continental). We identi- standard condition (i.e., the absence of Rb fusions), it was fied 49 newly described alleles on the Madeira archipelago, highly differentiated from Olost and Castellfollit del Boix in and 8 new alleles in the continental (Barcelona) Rb system, both the PCA and the estimates of individual ancestry. revealing that intraspecific PRDM9 diversity in M. m. domes- We then constructed phylogenetic trees for both the ticus is far greater than previously reported. This adds sub- Madeira and Barcelona Rb systems by the Neighbor-Joining stantially to previously described allelic diversity in natural approach using a new measure of pairwise genetic divergence populations of Mus musculus subspecies (28 distinct (see Material and Methods). The key feature of our approach PRDM9 alleles in M. m. domesticus, 34 in Mus musculus is that it considers point mutations, insertions, and deletions musculus, and 37 in M. m. castaneus) (Buard et al. 2014; of whole repeat units, as well as duplications of single repeat Capilla et al. 2014; Kono et al. 2014). units as a consequence of nonhomologous recombination, Most importantly, we detected contrasting evolutionary slippage, or related biological processes. This is of relevance patterns in the continental versus island chromosomal races since single-unit duplications and indels have different of the house mouse. This was reflected by unprecedented weights in their contribution to the genetic distance. This high levels of PRDM9 diversity in Madeiran house mice, an approach allowed us to effectively estimate genetic distances insular Rb system characterized by distinct metacentric races, among the polymorphic repeats of PRDM9 detected in the when compared with the Barcelona Rb system (where Rb Madeira and the Barcelona Rb systems. The distribution of fusions are not yet fixed within populations and are present haplotypes obtained showed greater diversification on the in polymorphic state). Our analysis of genetic diversity and island of Madeira than the Barcelona system (figs. 3D and positive selection, together with the phylogeographic recon- 4B). Although the phylogenetic analysis grouped PRDM9 struction of PRDM9, suggests that the variability observed alleles into three major phylogroups in the Barcelona system could be the result of, at least, two possible (not mutually (those representing the alleles 10A, 11B, and 12B; fig. 3D), the exclusive) scenarios in the Madeira system: 1) the current allelic diversification was much higher in the Madeira Island populations of Rb mice could reflect a complex pattern of with alleles 10A, 11B, 11S, 11K, and 13C representing the most introgression or multiple colonization events into the island common haplotypes (fig. 4B). Differences between both Rb and 2) meiotic impairment on hybridization of different pop- systems were also exemplified by higher genetic diversity ulations of Rb mice allowed PRDM9 alleles to diverge among estimates of the ZnF array on Madeira (5.406) than in the metacentric populations. Barcelona Rb system (2.307), the former being comparable to Initial studies based on mtDNA lineages (Bonhomme et al. the diversity across the whole Eurasia (6.123) (table 2). When 2011) indicated a recent common ancestry of all extant house plotting the chromosomal configuration in the phylogenetic mice populations, as well as a complex history owing to trees we detected that in both Rb systems (Madeira and founder events, genetic drift and secondary admixture. This Barcelona), standard mice showed a tendency to carry private was suggestive of two phases of mouse colonization westward alleles. This was exemplified by alleles 8A, 10B, 10C, 11B, 12C, across the Mediterranean basin; an initial event during the 12G, and 12F in Olost, St. Perpètua de Mogoda, and early progression of Neolithic human expansion (starting Castellfollit del Boix (fig. 3D), and 11W and 11X in Porto 12,000 years ago) followed by a second, more recent (some Santo (fig. 4A). hundred years ago), related to maritime trade. This progres- Taking advantage of previous surveys of wild captured sion can explain the presence of Mediterranean alleles in the mice in Europe and Asia, we further compared the allele Barcelona system and subsequent in situ diversification. variation found in the Madeira and Barcelona Rb sys- Madeira, on the other hand, has a more recent history of tems with those described in Eurasia (fig. 6). When the mouse colonization originated first from Danish Vikings in geographical distribution was plotted onto the global the 9th century followed by a second incursion by 15th cen- phylogenetic tree, Madeira haplotypes showed mixed tury Portuguese settlers (Gündüz, Auffray, et al. 2001; Britton- origins. The phylogenetic reconstruction of all three Davidian et al. 2007; Förster et al. 2009, 2013). Our analyses of sampled regions (and informed by evidence that the genetic diversity reveal a complex colonization pattern in genetic differentiation between Madeira and the Madeira, much more so than in the Barcelona Rb system Barcelona Rb system or even Eurasia is small, table 2) and leading to a PRDM9 diversity comparable to the diversity suggests that Madeira was colonized by a highly across the whole of Eurasia. Additionally, the phylogenetic different population, most likely as a result of multiple reconstruction of PRDM9 in Madeira indicated that the colonization/introgression events from multiple parts of same group of alleles appeared in different clusters indepen- Eurasia. dently. Altogether, our results suggest multiple colonization/ 1694 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE Geographical 2n distribution 2n=40 2n 2n=36 2n=32 2n=24 Madeira Geographical Barcelona distribution Eurasia FIG. 6. Phylogeographic depiction of PRDM9 at a global scale. Phylogenetic reconstruction including mice sequenced from the Madeira archi- pelago, Barcelona Rb system and Eurasia (data extracted from Kono et al. [2014] and Buard et al. [2014]). For each individual, both the diploid number (2n) and the geographical distribution (Madeira, Barcelona, and Eurasia) are indicated (see inserted legends). The diploid number is not known for Eurasian samples so they are represented in gray. introgression events into Madeira from Eurasia (including the higher in Rb mice than standard mice, both in Madeira Mediterranean basin). and the Barcelona Rb systems. This was true both for Under this scenario of multiple colonization events hypervariable codons and for the rest of PRDM9 and together with the extreme topography of Madeira sequence. (i.e., chromosomal races separated by mountain bar- An association between Rb fusions and PRDM9 diversity riers), the presence of Rb fusions might have resulted was also observed in the Barcelona Rb system, although the in meiotic impairment in hybrids that subsequently con- effects were mild, probably related to a differential demo- tributed to the high rate of diversification of PRDM9. graphic history. The Barcelona system represented approxi- Our analysis of genetic diversity and positive selection mately 10% of the novel PRDM9 alleles found in our study, supports this view. Since chromosomal fusions can act as with the highest diversity observed in standard mice located postzygotic barriers and can restrict gene flow among in the northern-eastern distribution of the sampled popula- chromosomal races (Franchini et al. 2010; Gimenez et al. tions. The allelic distribution together with the pattern of 2013; Capilla et al. 2014), this can result in high genetic genetic diversity detected by the PCA and ancestry analyses differentiation between populations, thus facilitating was suggestive of a larger connectivity among populations in the appearance of new PRDM9 variants. In fact, in the the Barcelona Rb system than in the case of Madeira. This was case of Madeira, each chromosomal race is geographi- exemplified by the 10A allele, which was present in all pop- cally localized and homogeneous with respect to the ulations especially in Rb mice, representing the greatest allele constitution of Rb fusions, and chromosomal hybrids homogeneity (nearly 95% of mice from Rb populations car- are rare (Britton-Davidian et al. 2000; Nunes et al. ried 10A). Considering the factors that might be responsible 2005). On the other hand, Madeira is extremely variable for the widespread distribution of this allele in the Barcelona regarding its genetic structure, as can be appreciated system, it is possible that the PRDM9 allele frequency pattern both from the PRDM9 phylogenetic tree and from the observed in this system could result from a population bot- high genetic differentiation between chromosomal tleneck that facilitated the formation and subsequent expan- races. Moreover, the relative deficit of heterozygous for sion of Rb fusions. According to the PCA and ancestry PRDM9 (despite the strong genetic mixing) points to analysis, standard mouse populations within the Barcelona strong selection for karyotype compatibility. In this con- system, such as Olost and Castellfollit del Boix, may represent text, differences in positive selection can be explained by the origin of the genetic diversity detected out of the presence of Rb fusions since dN/dS values were which Rb populations such as Gava and Sant Sadurnı 1695 Vara et al. . doi:10.1093/molbev/msz091 MBE d’Anoia subsequently appeared as differentiated populations. variations in the amino acid sequence of ZnF repeats can Consequently, Rb mice may have diverged into different, ge- modulate PRDM9 binding affinities to specific DNA motifs netically distinct populations resulting in a homogeneous dis- (Berg et al. 2010). In fact, differences in a single ZnF repeat can tribution of PRDM9 alleles with a predominance of the 10A affect the specificity of over 70% of the meiotic DSBs, leading allele. This is consistent with suggestions by Gündüz, Lopez- to the redistribution of recombination sites in a single gener- Fuster, et al. (2001) for the origin of the Barcelona system by ation (Smagulova et al. 2016). The rapid evolution of ZnF primary intergradation. This Rb system may therefore repre- sequences would therefore lead to rapid changes in the dis- sent an example of a radiation that could eventually lead to tribution of recombination sites across the genome by creat- the formation of a metacentric race without geographic iso- ing new DNA motifs that would be recognized with stronger lation (Sans-Fuentes et al. 2010). affinity by the new PRDM9 allelic variants. This will, in turn, compensate the loss of the primary DNA motifs by gene Functional Implications of Maintaining Conserved conversion during the repairing of crossovers during the early DNA-Binding Repeats stages of meiosis (Smagulova et al. 2016). Despite the ZnF variability observed in both Rb systems, we Although little is known at this stage about the dynamics detected significant deviations from Hardy–Weinberg equi- and the mechanistic constraints that affect recombination librium (i.e., a deficit of heterozygotes) where the PRDM9 landscapes at a finer scale in metacentric wild populations, heterozygous combinations most frequently encountered our results nonetheless suggest that Rb fusions might play a maintained amino acid sequence similarity from position role in bringing new selective pressures on PRDM9. These ZnF3 to ZnF6 in the ZnF array. These results suggest that outcomes provide impetus to examine the effect of genetic observed PRDM9 variability (either as a result of phylogeo- backgrounds on PRDM9 and the development of reproduc- graphic dynamics and/or population isolation triggered by tive isolation between natural populations. meiotic impairment of hybrids), can be subjected to func- tional constraints that facilitate the accumulation of allelic combinations that maintain recombination hotspot Materials and Methods symmetry. Sampling and DNA Extraction Under this scenario, it is important to take into account Our data set comprised a total of 395 wild-caught house the particular features of the PRDM9 protein (Davies et al. mice, 185 mice from the Barcelona Rb system (Medarde 2016; Tiemann-Boege et al. 2017; Paigen and Petkov 2018), et al. 2012), 199 mice from Madeira Island, and 11 from which include 1) its DNA-binding properties, 2) the rapid Porto Santo Island (table 1, fig. 1, and supplementary tables coevolution of the ZnF domain and the DNA motif that it 1 and 2, Supplementary Material online). Mice were karyo- recognizes, 3) asymmetric DSBs due to preferential binding of typed previously (Medarde et al. 2012; Chmatal et al. 2014, PRDM9 to new target sequences, and 4) hybrid sterility trig- supplementary fig. 1, Supplementary Material online). gered by asymmetric DSBs and interallelic PRDM9 combina- Animals were housed and treated in strict accordance with tions. Previous work on house mouse inbred strains found ethical guidelines approved by the Universitat Autonoma de that male sterility is a complex process that can be influenced Barcelona (Spain) and University of Lisbon (Portugal). by the heterozygous allelic combination of PRDM9 (Dzur- Genomic DNA was extracted from tissue biopsies preserved Gejdosova et al. 2012; Gregorova et al. 2018) and their inter- in absolute ethanol or fresh tissue using proteinase K diges- action with different X-linked and autosomal loci (Balcova tion (Sambrook et al. 1989). et al. 2016). Likewise, the presence of highly divergent PRDM9 The Barcelona Rb system included mice from 15 localities alleles in a heterozygous state can be detrimental due to selected for their geographical distribution and distinctive elevated levels of asymmetric DSBs between homologous chromosomal configurations (fig. 1 and supplementary table chromosomes (Davies et al. 2016). This is especially important 1, Supplementary Material online). This covered the full ex- since allelic combinations of PRDM9 with conserved DNA- tent of the chromosomal polymorphism previously described binding repeats were maintained in metacentric populations, for the Barcelona Rb system (Medarde et al. 2012). This in- suggestive of the presence of functional constraints. cluded localities with high frequencies of almost all metacen- In the light of our findings, the rapid formation of new tric chromosomes (2n ¼ 28–39), localities in the vicinity of Rb PRDM9 alleles in the presence of Rb fusions, can act syner- areas but without chromosomal fusions (2n ¼ 40), and local- gistically with the formation of new binding motifs distrib- ities geographically located between these population groups uted across the genome in different metacentric populations. containing intermediate diploid numbers (2n ¼ 39–40) (sup- This could result in functional constraints that might facilitate plementary table 1, Supplementary Material online). Four of the accumulation of PRDM9 heterozygous combinations that the localities considered in the study had mice with standard maintain a high degree of similarity in the amino acids re- karyotypes (2n ¼ 40, Castellfollit del Boix, Santa Perpètua de sponsible for protein–DNA binding, thus resulting in recom- Mogoda, Olost, and Caldes de Montbuı) and defined here as bination hotspot symmetry in homologous chromosomes. standard populations (supplementary table 1, Supplementary Mounting evidence in humans and mice has shown that Material online). The remaining localities (N ¼ 11) included variation in both the sequence and number of ZnF repeats mice with Rb fusions, defined here as Rb populations (2n ¼ strongly influence the distribution of meiotic DSBs (Berg et al. 28–39, Badalona Les Pobles, L’Ametlla de Segarra, El Papiol, 2010; Grey et al. 2011; Brick et al. 2012; Baker et al. 2015). Small Calafell, Sant Sadurnı d’Anoia, Cubelles, La Granada, El Prat de 1696 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE Llobregat, Viladecans, and Castelldefels). Mice from Rb pop- packages were used to distinguish allelic variants and to de- ulations were characterized by having between one to seven fine haplotypes of heterozygous mice: Phase (http://stephen- Rb fusions involving 14 different chromosomes (Rb[3.8], slab.uchicago.edu/software.html#phase) and Champuru v1.0 Rb[4.14], Rb[5.15], Rb[6.10], Rb[7.17], Rb[9.11], and (http://www.mnhn.fr/jfflot/champuru/). This allowed infer- Rb[12.13]), either in heterozygous or homozygous states (sup- ences on whether a sample was homozygous or heterozygous plementary table 1, Supplementary Material online). The Rb for its amino acid composition (Capilla et al. 2014). All allele fusions were present as chromosomal polymorphisms in the sequences were translated using ExPasy: SIB bioinformatics Barcelona Rb system, with varying frequencies (supplemen- resource portal (Artimo et al. 2012). tary table 1, Supplementary Material online). Additionally, we included 118 published PRDM9 sequen- The Madeira system was represented by mice from 34 ces of wild-derived M. m. domesticus and M. m. castaneus localities that included all six chromosomal races occurring specimens drawn from various locations in Eurasia (Buard on the island of Madeira (as originally defined by Britton- et al. 2014; Capilla et al. 2014; Kono et al. 2014) for the Davidian et al. [2000]) (supplementary table 2, phylogeographic analysis of PRDM9 at a global scale. Supplementary Material online): PSAN (2n ¼ 22), PADC Overall, the phylogenetic comparative study included three (2n ¼ 24–28), PEDC (2n ¼ 23–26), PLDB (2n ¼ 24), PPOD M. m. domesticus sampling regions: the Barcelona Rb system, (2n ¼ 27–28), and PSVI (2n ¼ 26–27). A total of 20 different the Madeira archipelago (including both Madeira and Porto metacentrics in homozygous state are described in Madeira, 9 Santo Islands) and Eurasia (previously published data; Buard of which have not been detected elsewhere (Rb[2.19], Rb[4.5], et al. 2014; Kono et al. 2014) (table 1). Rb[4.16], Rb[5.18], Rb[9.18], Rb[11.12], Rb[11.19], Rb[13.17], and Rb[15.18]), highlighting the uniqueness of this system PRDM9 Allele Classification (see supplementary table 3, Supplementary Material online, Following the nomenclature used in previous studies (Buard for details on the chromosomal composition of chromosomal et al. 2014; Capilla et al. 2014; Kono et al. 2014), PRDM9 alleles races). Mice with the standard karyotype (2n ¼ 40) from the were classified using the number of ZnF repeats and the ex- neighboring island of Porto Santo (four localities) were also tent of amino acid variation in the highly variable positions included in our investigation (fig. 1 and supplementary table 1, þ3, and þ6 of each ZnF repeat (table 1 and supplemen- 2, Supplementary Material online). tary table 5, Supplementary Material online). Each ZnF repeat was identified by a number (from #3 to #34) based on its Prdm9 Sequencing amino acid sequence (supplementary table 5, Supplementary The entire ZnF array is encoded by the last exon of the Prdm9 Material online). Subsequently, the combination of different gene that extends from the second ZnF repeat towards the C- of ZnF repeats was classified as distinct PRDM9 alleles (fig. 2). terminal domain. We amplified this region in all mice by PCR Moreover, for each PRDM9 allele found in heterozygous state, using the primers ZFA-L_F (forward) and ZFA-L_R (reverse) the amino acid conservation of the ZnF repeats located in described by Kono et al. (2014) (supplementary table 4, positions ZnF3–ZnF6 of the array was calculated based on the Supplementary Material online) using the ExTaq (TaKaRa) hypervariable sites (1, þ3, and þ6; Kono et al. 2014) (sup- following the protocol in Capilla et al. (2014). Briefly, the plementary table 6, Supplementary Material online). PCR conditions were 95  C (3 min), 30 cycles of 95  C (30 s), 64  C (30 s), and followed by 72  C (90 s). The ZnF array SNP Genotyping exhibits not only length polymorphisms (variation of the SNP genotyping data included 50 mice from 5 localities of the number of ZnF repeats) but also amino acid variation (alleles Barcelona Rb system: 26 mice from three standard popula- with the same number of ZnF repeats but different amino tions (Castellfollit del Boix, N ¼ 8; Santa Perpètua de Mogoda, acid composition). As a consequence, it is necessary in each N ¼ 13; and Olost, N ¼ 5) and 24 mice from two Rb popula- case to identify the length of the allele and subsequently the tions (Sant Sadurnı d’Anoia, N ¼ 12 and Castelldefels, amino acidic composition of the sequence. To this end, PCR N ¼ 12). Genomic DNA was extracted using a standard pro- products were separated on 1% agarose gels. This allowed tocol with proteinase K (Sambrook et al. 1989). Subsequently, both the ZnF array length to be discerned (homozygous mice were genotyped using the Mega Mouse Universal and heterozygous) and optimal DNA recovery from bands Genotyping Array (Morgan et al. 2015), which consist of for subsequent sequencing (supplementary fig. 2, 77,808 evenly distributed SNP markers built on the Illumina Supplementary Material online). Infinium platform. SNPs were filtered to remove markers with To determine the sequences of ZnF repeats, PCR products missing values >5% threshold (i.e., markers that do not fit of the homozygous mice were subjected to bidirectional Hardy–Weinberg expectations, P < 105) using PLINK ver- Sanger sequencing with the same primers used for amplifica- sion 1.9 (Purcell et al. 2007). The proportion of missing data tion: primers ZFA-L_F and ZFA-L_R. Only sequences repro- and heterozygosity per locus and per sample were also cal- ducible in at least three independent amplification reactions culated to evaluate possible bias. This resulted in a final data were included in the data analyses. Bands of heterozygous set of 63,344 informative SNPs distributed across all chromo- mice were purified from the gels using the Nucleospin Gel somes, with the exception of the Y. Using this data set, the and PCR Clean-up kit (Macherey-Nagel) and subsequently ADMIXTURE software (Alexander et al. 2009) was used to sequenced. Sequences from heterozygous mice were analyzed estimate individual ancestry and admixture proportions as- using the Bioedit 7.2.5 package. Two different software suming K populations based on a maximum likelihood 1697 Vara et al. . doi:10.1093/molbev/msz091 MBE method. Analyses were run only for SNPs with a greater than in the SNP genotyping (computed from the 20 Mb around 95% genotype call. The numbers of clusters (K) evaluated here the PRDM9 nucleotide sequence) and the mean ZnF genetic ranged from 1 to 10. To further evaluate the final K value, an distance for the same individuals. Agreement between the Evanno’s DK was applied (Evanno et al. 2005). Subsequently, genetic distances was defined in terms of Pearson correlation three different K values (K ¼ 3, 4, 5) showed the lowest like- of the pairwise distances of all individuals for which both lihood values (supplementary fig. 3, Supplementary Material genotype and PRDM9 sequence data were available. The cor- online). ADMIXTURE analyses were plotted using an R frame- relation was computed numerically for a grid of weight values work. In addition, PCA was implemented using a module of with a step of 0.25. In case of ties, the combination with the PLINK 1.9. lowest ratios wi/wm and ws/wm was chosen. The final choice of weights maximizing the agreement between genetic dis- Phylogenetic Analysis of PRDM9 Alleles tances is wm ¼ 1 for nucleotide point mutations, wi ¼ 3.5 for Phylogenetic trees for the nucleotide sequences of PRDM9 unit insertions/deletions, and ws ¼ 1.75 for single-unit dupli- alleles were built via a Neighbor-Joining approach informed cations. All pairs of sequences were pairwise aligned with our by a new measure of pairwise genetic divergence between approach, using the genetic distance corresponding to these ZnF repeats. Briefly, divergence between pairs of ZnF repeats weights. was computed by first masking codons corresponding to the The PRDM9 phylogenetic trees were reconstructed from hypervariable amino acid positions (1, þ3, þ6). All repeat the molecular distance discussed above, using the Neighbor- units share a high sequence similarity and have the same Joining method implemented in the R library APE (Paradis length, it was therefore straightforward to align these units et al. 2004). Trees were rooted at midpoint using the R pack- and compute the pairwise genetic (Hamming) distances age phangorn (Schliep 2011). l(t,t0 ) between any two units t and t0 . Then, the evolution of the ZnF repeats was modeled as a Markov process involv- Genetic Diversity and Positive Selection ing nucleotide mutation, insertion and deletions of a single The PRDM9 alleles were tested for deviations from Hardy– repeat unit as well as single-unit duplication/slippage. Under Weinberg proportions in all populations with more than ten this framework, we defined the genetic distance between sequenced individuals. The significance of the deficit of het- repeats as an edit distance, that is, as the minimum cost to erozygous individuals for PRDM9 alleles with respect to change a repeat to another through a series of elementary Hardy–Weinberg predictions was assessed by an exact per- operations. The edit distance is minimized by one or possibly mutation test based on 1,000,000 permutations. several alignments between repeats r ¼ (r1; r2; r3; . . .) and r0 ¼ Pairwise nucleotide diversity was estimated as the mean (r0 1 ; r0 2 ; r0 3 ; . . .). The cost was defined as the weighted sum of number of nucleotide differences per base averaged across all all contributing elements: 1) point mutations, small indels, pairs of aligned nucleotide sequences, ignoring gaps in the and other within-unit processes (cost l[rj, r0 j0 ] for an align- pairwise sequence alignment. FST values between multiple ment between the jth and the jth units of the two repeats, population/chromosomal races were computed using the weight wm); 2) insertions and deletions of whole units (cost 1, formula described by Hudson et al. (1992). This computed weight wi); and 3) single-unit duplication/slippage (cost 1, the within-population diversity as the average of the within- weight ws; plus an extra cost l[rj, rj1] or l[r0 j0 , r0 j0 1 ] for population diversities of all population (in order to weight the the point mutations and small indels between neighboring contribution of each population equally) and permitted the duplicated units, which was weighted by wm). An R imple- removal of hypervariable codons from the sequence. mentation of the algorithms described here is provided at Selection was inferred from dN/dS estimates averaged https://github.com/lucaferretti/RepeatDistance. across all pairwise sequence alignments as follows. First, we Subsequently, each pair of ZnF repeats was aligned using a computed all the pairwise alignments using the Needleman– modified Needleman–Wunsch algorithm corresponding to Wunsch algorithm described before. For each pair, we esti- the model considered, with different costs for each type of mated dN as the number of nonsynonymous mutations per mutations (wm for point mutations, wi for insertions/dele- nonsynonymous site, considering either hypervariable codons tions, and ws for slippage). The genetic distance between two only (1, þ3, þ6) or excluding them, and ignoring gaps in sequences was defined as the weighted sum of the cost of the pairwise alignment. The fraction of nonsynonymous sites each mutation between the sequences. Our selection of costs was estimated by randomization conditional on the nucleo- was as follows: 1) Mismatch between units rju and r0 j0 : wml(rj, tide frequencies of the sequence. Then, dS was similarly esti- r0 j0 ) and 2) insertion of a whole unit rj in repeat r: min[wi, ws þ mated as the number of synonymous mutations per wml(rj, rj1), ws þ wml(r0 j0 , r0 j1 )]. The best alignment synonymous site (since there were often no synonymous according to this choice of costs corresponds precisely to mutations in hypervariable codons, we used synonymous the alignment that minimizes the sum of costs of mutations mutations in the other codons as an alternative estimate of between sequences; hence the minimum cost obtained from neutral rates). Finally, we averaged dN and dS across all se- the Needleman–Wunsch algorithm provides the distance be- quence pairs and estimated dN/dS as the ratio of these aver- tween repeats defined above. ages. This corresponds to the approximate Maximum Weight parameters were chosen among the ones that Likelihood estimate of dN/dS for Poisson-distributed muta- maximized the agreement between mean genetic distances tions. Neutral evolution would result in an expected value of in the subset of 50 mice from the Barcelona system included dN/dS around 1. The 95% confidence intervals were 1698 PRDM9 Diversity at Fine Geographical Scale . doi:10.1093/molbev/msz091 MBE computed by bootstrapping (sampling 500 random sets of Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation codons with replacement). We computed dN/dS both across of ancestry in unrelated individuals. Genome Res. 19(9):1655–1664. Artimo P, Jonnalagedda M, Arnold K, Arnold K, Baratin D, Csardi G, de all sequences and within/between specific subsets (such as Castro E, Duvaud S, Flegel V, Fortier A, et al. 2012. ExPASy: SIB sequences from individuals with standard karyotype or with bioinformatics resource portal. Nucleic Acids Res. 40:597–603. Rb fusions). All analyses of genetic diversity were imple- Baker CL, Kajita S, Walker M, Saxl RL, Raghupathy N, Choi K, Petkov PM, mented via custom R scripts. Paigen K. 2015. PRDM9 drives evolutionary erosion of hotspots in Mus musculus through haplotype-specific initiation of meiotic re- Data Accessibility combination. PLoS Genet. 11(1):e1004916. Baker Z, Schumer M, Haba Y, Bashkirova L, Holland C, Rosenthal GG, Data underlying this article are available on GeneBank (from Przeworski M. 2017. Repeated losses of PRDM9-directed recombi- MK848086 to MK848149). nation despite the conservation of PRDM9 across vertebrates. eLife 6:e24133. Supplementary Material Balcova M, Faltusova B, Gergelits V, Bhattacharyya T, Mihola O, Trachtulec Z, Knopf C, Fotopulosova V, Chvatalova I, Gregorova S, Supplementary data are available at Molecular Biology and et al. 2016. Hybrid sterility locus on chromosome X controls meiotic Evolution online. recombination rate in mouse. PLoS Genet. 12(4):e1005906. Barbash DA, Siino DF, Tarone AM, Roote J. 2003. A rapidly evolving Acknowledgments MYB-related protein causes species isolation in Drosophila. Proc Natl Acad Sci U S A. 100(9):5302–5307. C.V. and L.C. were supported by FPI predoctoral fellowships Baudat F, Buard J, Grey C, Fledel-Alon A, Ober C, Przeworski M, Coop G, from the Ministerio de Economıa y Competitividad de Massy B. 2010. Prdm9 is a major determinant of meiotic recom- (MINECO) (BES-2011-047722 and BES-2015-072924, respec- bination hotspots in human and mice. Science 21:836–840. tively). This study was partially supported by grants from Berg IL, Neumann R, Lam K-W, Sarbajna S, Odenthal-Hesse L, May CA, MINECO to J.V. (CGL2010-15243) and to A.R.-H. (CGL2010- Jeffreys AJ. 2010. PRDM9 variation strongly influences recombination hot-spot activity and meiotic instability in humans. Nat Genet. 20170, CGL2014-54317-P, CGL2017-83802-P, and BFU2015- 42(10):859–863. 71786-REDT). L.F. is supported by funding from BBSRC grant Bhattacharyya T, Gregorova S, Mihola O, Anger M, Sebestova J, Denny P. BBS/E/I/00007039. S.I.G. (fellowship SFRH/BPD/88854/2012), 2013. Mechanistic basis of infertility of mouse intersubspecific J.B.S. and M.L.M. were supported by Fundaç~ao para a Ci^encia hybrids. Proc Nat Acad Sci U S A. 110:468–477. e a Tecnologia (FCT, PTDC/BIA-EVF/116884/2010). Financial Billings T, Parvanov ED, Baker CL, Walker M, Paigen K, Petkov PM. 2013. DNA binding specificities of the long zinc-finger recombination pro- support from “Alianza 4 Universidades” to R.A.S.-G. is also tein PRDM9. Genome Biol. 14(4):R35. acknowledged. Thanks are due for the financial support of Bonhomme F, Orth A, Cucchi T, Rajabi-Maham H, Catalan J, Boursot P, CESAM (UID/AMB/50017 - POCI-01-0145-FEDER-007638), Auffray J-C, Britton-Davidian J. 2011. Genetic differentiation of the FCT/MCTES through national funds (PIDDAC), and the house mouse around the Mediterranean basin: matrilineal foot- cofunding by the FEDER, within the PT2020 Partnership prints of early and late colonization. Proc Biol Sci. 278(1708):1034–1043. Agreement and Compete 2020. The authors acknowledge Boulton A, Myers RS, Redfield RJ. 1997. The hotspot conversion paradox F. Pardo-Manuel de Villena for his assistance in SNP genotyp- and the evolution of meiotic recombination. Proc Natl Acad Sci U S ing. We also acknowledge the valuable help in field work and/ A. 94(15):8058–8063. or karyotyping and genotyping of Joaquim Tapisso, Ana Brick K, Smagulova F, Khil P, Camerini-Otero RD, Petukhova GV. 2012. Cerveira, Maria da Graça Ramalhinho, Janice Britton- Genetic recombination is directed away from functional genomic elements in mice. Nature 485(7400):642–645. Davidian, Pedro Martins da Silva, Rita Monarca, Nuria Britton-Davidian J, Catalan J, Lopez J, Ganem G, Nunes AC, Ramalhinho Medarde, Marta Pla-Bagaria, and Laura Palacios-Fernandez. MG, Auffray JC, Searle JB, Mathias ML. 2007. Patterns of genic diver- sity and structure in a species undergoing rapid chromosomal radi- Author Contributions ation: an allozyme analysis of house mice from the Madeira archipelago. Heredity 99(4):432–442. A.R.-H. conceived and devised the study. C.V., L.C., L.F., A.L., Britton-Davidian J, Catalan J, Ramalhinho MG, Ganem G, Auffray J-C, R.A.S.-G., S.I.G., M.L.M., J.B.S., and A.R.-H. contributed to the Capela R, Biscoito M, Searle JB, Mathias ML. 2000. Rapid chromo- design of the methodological approaches. C.V., L.C., R.A.S.-G., somal evolution in island mice. Nature 403(6766):158. B.F.-S., J.B.-R., and A.R.-H. contributed to PRDM9 sequencing Buard J, Rivals E, Dunoyer de Segonzac D, Garres C, Caminade P, de experiments and allele classification. C.V., G.A.-L., L.F., and Massy B, Boursot P. 2014. Diversity of Prdm9 zinc finger array in wild mice unravels new facets of the evolutionary turnover of this coding A.R.-H. contributed to SNP genotyping and phylogenetic minisatellite. PLoS One 9(1):e85021. analysis of PRDM9 alleles. L.F. and A.L. developed the ap- Capilla L, Garcia Caldes M, Ruiz-Herrera A. 2016. Mammalian meiotic proach for the computation of genetic distances between recombination: a toolbox for genome evolution. Cytogenet Genome complex repeats. S.I.G., J.V., J.B.S., M.L.M., and A.R.-H. contrib- Res. 150(1):1–16. uted to reagents and data collection. C.V., L.C., and A.R.-H. Capilla L, Medarde N, Alemany-Schmidt A, Oliver-Bonet M, Ventura J, Ruiz-Herrera A. 2014. Genetic recombination variation in wild wrote the first draft of the manuscript with contributions of Robertsonian mice: on the role of chromosomal fusions and all authors. All authors read and approved the final version of Prdm9 allelic background. Proc R Soc B 281:1–18. the manuscript. Chmatal L, Gabriel SI, Mitsainas GP, Martınez-Vargas J, Ventura J, Searle JB, Schultz RM, Lampson MA. 2014. Centromere strength provides the cell biological basis for meiotic drive and karyotype evolution in References mice. Curr Biol. 24(19):2295–2300. Adolph S, Klein J. 1981. Robertsonian variation in Mus musculus from Coop G, Myers SR. 2007. Live hot, die young: transmission distortion in Central Europe, Spain, and Scotland. J Hered. 72(3):219–221. recombination hotspots. PLoS Genet. 3:377–386. 1699 Vara et al. . doi:10.1093/molbev/msz091 MBE Coyne JA, Orr HA. 1998. The evolutionary genetics of speciation. Philos polymorphism zone in the house mouse. Heredity (Edinb). Trans R Soc Lond B Biol Sci. 353(1366):287–305. 109(2):78–89. Davies AB, Hatton E, Altemose N, Hussin JG, Pratto F, Zhang G, Hinch Mihola O, Trachtulec Z, Vlcek C, Schimenti JC, Forejt J. 2009. A mouse AG, Moralli D, Biggs D, Camerini-Otero RD, et al. 2016. Re-engineer- speciation gene encodes a meiotic histone H3 methyltransferase. ing the zinc fingers of PRDM9 reverses hybrid sterility in mice. Nature Science 323(5912):373–375. 530(7589):171–176. Morgan AP, Fu C-P, Kao C-Y, Welsh CE, Didion JP, Yadgary L, Hyacinth L, Dzur-Gejdosova M, Simecek P, Gregorova S, Bhattacharyya T, Forejt J. Ferris MT, Bell TA, Miller DR. 2015. The mouse universal genotyping 2012. Dissecting the genetic architecture of F1 hybrid sterility in array: from substrains to subspecies. G3 6:263–279. house mice. Evolution 66(11):3321–3335. Nunes AC, Britton-Davidian J, Catalan J, Ramalhinho MG, Capela R, Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters Mathias ML, Ganem G. 2005. Influence of physical environmental of individuals using the software STRUCTURE: a simulation study. characteristics and anthropogenic factors on the position and struc- Mol Ecol. 14(8):2611–2620. ture of a contact zone between two chromosomal races of the house Faria R, Navarro A. 2010. Chromosomal speciation revisited: rearranging mouse on the island of Madeira (North Atlantic, Portugal). J theory with pieces of evidence. Trends Ecol Evol (Amst). Biogeogr. 32(12):2123–2134. 25(11):660–669. Oliver PL, Goodstadt L, Bayes JJ, Birtle Z, Roach KC, Phadnis N, Beatson Farre M, Micheletti D, Ruiz-Herrera A. 2013. Recombination rates and SA, Lunter G, Malik HS, Ponting CP. 2009. Accelerated evolution of genomic shuffling in human and chimpanzee—a new twist in the the Prdm9 speciation gene across diverse metazoan taxa. PLoS chromosomal speciation theory. Mol Biol Evol. 30(4):853–864. Genet. 5(12):e1000753. Flachs P, Mihola O, Simecek P, Gregorova S, Schimenti JC, Matsui Y, Orr HA, Turelli M. 2001. The evolution of postzygotic isolation: accu- Baudat F, de Massy B, Pialek J, Forejt J, et al. 2012. Interallelic and mulating Dobzhansky–Muller incompatibilities. Evolution intergenic incompatibilities of the Prdm9 (Hst1) gene in mouse hy- 55(6):1085–1094. brid sterility. PLoS Genet. 8(11):e1003044. Paigen K, Petkov PM. 2018. PRDM9 and its role in genetic recombina- Förster DW, Günduz _I, Nunes AC, Gabriel S, Ramalhinho MG, Mathias tion. Trends Genet. 34(4):291–300. ML, Britton-Davidian J, Searle JB. 2009. Molecular insights into the Paradis E, Claude J, Strimmer K. 2004. APE: analyses of phylogenetics and colonization and chromosomal diversification of Madeiran house evolution in R language. Bioinformatics 20(2):289–290. mice. Mol Ecol. 18(21):4477–4494. Pavlova SV, Searle JB. 2018. Chromosomes and speciation in mammals. Förster DW, Mathias ML, Britton-Davidian J, Searle JB. 2013. Origin of the In: Zachos FE, Asher RJ, editors. Mammalian evolution, diversity and chromosomal radiation of Madeiran house mice: a microsatellite systematics. Berlin (Germany): De Gruyter. p. 17–38. analysis of metacentric chromosomes. Heredity (Edinb). Phadnis N, Orr HA. 2009. A single gene causes both male sterility and 110(4):380–388. segregation distortion in Drosophila hybrids. Science Franchini P, Colangelo P, Solano E, Capanna E, Verheyen E, Castiglia R. 323(5912):376–379. 2010. Reduced gene flow at pericentromeric loci in a hybrid zone Pialek J, Hauffe HC, Searle JB. 2005. Chromosomal variation in the house involving chromosomal races of the house mouse Mus musculus mouse. Biol J Linn Soc. 84(3):535–563. domesticus. Evolution 64:2020–2032. Presgraves DC, Balagopalan L, Abmayr SM, Orr HA. 2003. Adaptive Gimenez MD, White TA, Hauffe HC, Panithanarak T, Searle JB. 2013. evolution drives divergence of a hybrid inviability gene between Understanding the basis of diminished gene flow between hybrid- two species of Drosophila. Nature 423(6941):715–719. izing chromosome races of the house mouse. Evolution Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, 67(5):1446–1462. Maller J, Sklar P, de Bakker PIW, Daly MJ, et al. 2007. PLINK: a tool set Gregorova S, Gergelits V, Chvatalova I, Bhattacharyya T, Valiskova B, for whole-genome association and population-based linkage analy- Fotopulosova V, Jansa P, Wiatrowska D, Forejt J. 2018. Modulation ses. Am J Hum Genet. 81(3):559–575. of Prdm9-controlled meiotic chromosome asynapsis overrides hy- Qiao H, Chen JK, Reynolds A, Höög C, Paddy M, Hunter N. 2012. brid sterility in mice. eLife 7:e34282. Interplay between synaptonemal complex, homologous recombina- Grey C, Barthès P, Chauveau-Le Friec G, Langa F, Baudat F, de Massy B. tion, and centromeres during mammalian meiosis. PLoS Genet. 2011. Mouse PRDM9 DNA-binding specificity determines sites of 8(6):e1002790. histone H3 lysine 4 trimethylation for initiation of meiotic recom- Sambrook J, Fritsch EF, Maniatis T. 1989. Molecular cloning: a laboratory bination. PLoS Biol. 9(10):e1001176. manual. New York: Cold Spring Harbor Laboratory Press. Gündüz _I, Auffray J-C, Britton-Davidian J, Catalan J, Ganem G, Sanchez-Guillen RA, Capilla L, Reig-Viader R, Martınez-Plana M, Pardo- Ramalhinho MG, Mathias ML, Searle JB. 2001. Molecular studies Camacho C, Andres-Nieto M, Ventura J, Ruiz-Herrera A. 2015. On on the colonization of the Madeiran archipelago by house mice. the origin of Robertsonian fusions in nature: evidence of telomere Mol Ecol. 10(8):2023–2029. shortening in wild house mice. J Evol Biol. 28(1):241–249. Gündüz I, Lopez-Fuster MJ, Ventura J, Searle JB. 2001. Clinal analysis of a Sans-Fuentes MA, Garcıa-Valero J, Ventura J, Lopez-Fuster MJ. 2010. chromosomal hybrid zone in the house mouse. Genet Res. Spermatogenesis in house mouse in a Robertsonian polymorphism 77(1):41–51. zone. Reproduction 140(4):569–581. Hudson RR, Slatkin M, Maddison WP. 1992. Estimation of levels of gene Schliep KP. 2011. Phangorn: phylogenetic analysis in R. Bioinformatics flow from DNA sequence data. Genetics 132(2):583–589. 27(4):592–593. Jeffreys AJ, Cotton VE, Neumann R, Lam KG. 2013. Recombination reg- Schwartz JJ, Roach DJ, Thomas JH, Shendure J. 2014. Primate evolution of ulator PRDM9 influences the instability of its own coding sequence the recombination regulator PRDM9. Nat Commun. 5:4370. in humans. Proc Natl Acad Sci U S A. 110(2):600–605. Smagulova F, Brick K, Pu Y, Camerini-Otero RD, Petukhova GV. 2016. Kono H, Tamura M, Osada N, Suzuki H, Abe K, Moriwaki K, Ohta K, The evolutionary turnover of recombination hot spots contributes Shiroishi T. 2014. Prdm9 polymorphism unveils mouse evolutionary to speciation in mice. Genes Dev. 30(3):266–280. tracks. DNA Res. 21(3):315–326. Tiemann-Boege I, Schwarz T, Striedner Y, Heissl A. 2017. The consequen- Masly JP, Jones CD, Noor MAF, Locke J, Orr HA. 2006. Gene transposition ces of sequence erosion in the evolution of recombination hotspots. as a cause of hybrid sterility in Drosophila. Science Philos Trans R Soc B 372(1736):20160462. 313(5792):1448–1450. Ting CT, Tsaur SC, Wu ML, Wu CI. 1998. A rapidly evolving Medarde N, Lopez-Fuster MJ, Mu~ noz-Mu~noz F, Ventura J. 2012. Spatio- homeobox at the site of a hybrid sterility gene. Science temporal variation in the structure of a chromosomal 282(5393):1501–1504. 1700