Genomic resources for wild populations of the house mouse, Mus musculus and its close relative Mus spretus

Wild populations of the house mouse (Mus musculus) represent the raw genetic material for the classical inbred strains in biomedical research and are a major model system for evolutionary biology. We provide whole genome sequencing data of individuals representing natural populations of M. m. domesticus (24 individuals from 3 populations), M. m. helgolandicus (3 individuals), M. m. musculus (22 individuals from 3 populations) and M. spretus (8 individuals from one population). We use a single pipeline to map and call variants for these individuals and also include 10 additional individuals of M. m. castaneus for which genomic data are publically available. In addition, RNAseq data were obtained from 10 tissues of up to eight adult individuals from each of the three M. m. domesticus populations for which genomic data were collected. Data and analyses are presented via tracks viewable in the UCSC or IGV genome browsers. We also provide information on available outbred stocks and instructions on how to keep them in the laboratory.


Background & Summary
The house mouse (Mus musculus) has a long-standing history as a model system in genetics and biomedical research, with many classical inbred strains available for purchase worldwide. By comparing the genetic make-up of classical inbred strains with those of mice collected in the wild, it became clear that classical inbred strains represent complex genomic mixtures with contributions from different subspecies and species of Mus 1-3 . While some of this genomic mixture stems from captive breeding of mice from different parts of the world during the early establishment of inbred strains, admixture 4,5 and introgression of genomic material across subspecies and species [6][7][8][9] also occurs in the wild and is thus likely to contribute to the genomic complexity observed in inbred strains. Classical inbred strains were found to exhibit a much reduced amount of genetic variation compared to their wild mice ancestors 10 .
mouse subspecies and one outgroup, b) tissue-specific RNAseq data from three M. m. domesticus populations, and c) details on animal husbandry of wild house mice in a Supplementary File and d) some general analyses and browser tracks for visualization. The genomic dataset is summarized using basic descriptive statistics, such as F ST (ref. 42), π (nucleotide diversity 43 ) and Tajima's D 44 as measures of population differentiation and selection, as well as a description of copy-number variation based on sequencing read depth. The RNAseq data are summarized as normalized RNAseq read coverage for each base pair of the genome. All statistics are made available as genome browser tracks (bed and bigWig files) to allow close visual inspection of any genomic region of interest in the UCSC browser 45 or other visualization software such as IGV 46,47 .

Methods
Sampling procedure and sampling locations The location of populations used for re-sequencing in this study are depicted in Fig. 1. They include three M. m. domesticus populations from Western Europe and Iran, the island subspecies M. m. helgolandicus, three M. m. musculus populations from the Czech Republic, Kazakhstan and Afghanistan and a M. spretus population from Spain. Populations were sampled between 2003 and 2012. The DNA samples used for genome sequencing were obtained either directly from wild caught animals, or from the first or second generation of out-breeding in our animal facility, i.e., they are expected to represent full wild type variation. Some aspects of the genome sequences obtained from the three M. m. domesticus populations, as well as the island subspecies M. m. helgolandicus, have previously been described 26,38 . House mice form naturally extended family groups at a given location including breeding among relatives 48,49 . To obtain an unbiased population sample from a given region ideally requires sampling mice in a way to avoid catching related animals. Therefore, we aimed to collect only a single mouse per trapping location (or, for the purpose of setting up breeding colonies, one male and one female per location) and selected the next trapping location 500 m to 1 km apart. The whole area sampled ideally comprises a diameter of about 50 km. However, depending on the local conditions, following this sampling regime precisely was not always possible. Moreover, some samples were provided by collaborators who have only recorded the general area for trapping, but not the exact location (e.g., the mice from Kazakhstan). Other trapping locations had regional limitations. For example, the island of Heligoland is only 1.7 km 2 in size, or the military Camp Marmal, Mazar-e-Sharif (Afghanistan) is only 8.7 km 2 in size. In both cases, mice were collected in different localities on the island or military base respectively 50 . In Supplementary Table 1 we provide exact location information, as much as it is available for all animals involved in the study, as either directly having been sequenced, or as having been parent to one of the animal facility-born offspring of wild mice (see below).  To complement the genomic resources generated in our laboratory, we also re-analyzed previously published genome sequencing data from M. m. castaneus that were collected in the northwest Indian state of Himachal Pradesh 51 .
Mice were either caught in snap traps, or live traps ('Mäusewippfalle' No. 3451002, Firma Ehlert & Partner, 53859 Niederkassel, Germany) or were found dead after rodenticide-based pest management. Cervical dislocation was used to sacrifice mice caught in live traps in the field. Transportation of live mice to the animal facility, maintenance and handling were conducted in accordance with German animal welfare law (Tierschutzgesetz) and FELASA guidelines. Permits for keeping mice were obtained from the local veterinary office 'Veterinäramt Kreis Plön' (permit number: 1401-144/PLÖ-004697).

M. m. domesticus breeding scheme
For the M. m. domesticus populations our aim was to generate an RNAseq dataset from the same individuals for which we generated the DNAseq data. In order to standardize mice to the same sex, age and environmental conditions, we did not use wild caught mice but bred wild caught mice for one to two generations in our animal facility. Supplementary Table 2 shows the breeding scheme for the M. m. domesticus mice in the study. In one case, we caught a pregnant female in the wild and used its male offspring born in the facility for DNA and RNA sequencing. For the Iranian mice, we included two wild caught individuals in the DNAseq study (AH15 and AH23), for which we did not generate RNAseq data. For two additional Iranian individuals in the DNA sequencing study, tissue samples were lost and thus no corresponding RNAseq data exist. To compensate for this, we included four individuals representing male offspring from male AH15 and male AH23 respectively (both individuals are part of the DNAseq dataset; see Table 1 (available online only)). Two male offspring from each of these two males are represented as biological replicates in the RNAseq dataset.

Procedures for wild mouse handling
We established wild-derived outbred populations for M. m. domesticus (France, Germany) and M. m. musculus (Kazakhstan and Czech Republic). For the first 11-14 generations live mice obtained from the wild were set up in a cyclical breeding scheme aiming at maintaining maximum genetic variability over time. They were then partly refreshed with newly collected mice from the same original area and a HAN rotational breeding scheme 52 was established. Individuals from each of the 4 outbred populations are maintained at the Max-Planck Institute in Plön, Germany, and are available from the authors upon request.
Wild mice are considerably more agile than classical inbred strain mice. Environmental enrichment is necessary and strongly reduces agitated stereotype behavior in wild mice kept under laboratory conditions. Standard mouse chow is provided ad libitum (e.g., Altromin 1,324 from ALTROMIN, 32,791 Lage, Germany). Further details on mouse handling and breeding are provided in Supplementary Material Text 1.

Molecular methods
DNA and RNA extraction. DNA was extracted from liver, spleen, or ear samples using salt extraction 53 or DNeasy kits (Qiagen, Hilden, Germany). RNA was extracted only for M. m. domesticus mice from Germany (8 individuals), France (8 individuals) and Iran (8 individuals) using mostly the same individuals, which are included in the whole genome sequencing study (see details above and Supplementary Table 2).
Mice designated for RNA extraction were housed alone after weaning and were routinely visually inspected for health and vigor. Only male mice were included in the study, sacrificed at 12 weeks. All mice were fed standard mouse chow ad libitum (Altromin 1,324 from ALTROMIN, 32,791 Lage, Germany). Mice were sacrificed by CO 2 asphyxiation. The coat was sprayed with 75% EtOH to reduce loose hair contaminating the organs during dissection of the animal. Organs were extracted in a specific order to improve comparability. Organs were shock frozen in liquid nitrogen and stored at − 80°until RNA was extracted.
The Trizol reagent (Life Technologies, Carlsbad, California, USA) was used according to manufacturers instructions to extract RNA from each organ ( Table 2 (available online only)). With the exception of the liver, for which we used only right and left medial lobe, whole organs were processed to minimize heterogeneity of the sample. The extracted RNA was quantified on a Nanodrop and analyzed for integrity on the Agilent Bioanalyzer.
DNA sequencing library preparation. All populations apart of the Afghanistan population were sequenced using the same protocol (see below) in the following batches: Batch 1 included the M. m. domesticus populations from France, Germany and Iran (locations 1, 2 and 4 in Fig. 1). Batch 2 included the 3 mice from Heligoland (location 3 in Fig. 1). Batch 3 included the M. m. musculus populations from the Czech Republic and Kazakhstan (locations 5 and 6 in Fig. 1) as well as the M. spretus population from Spain (location S in Fig. 1). Batch 4 included the M. m. musculus mice from the Afghanistan population (location 7 in Fig. 1), which were sequenced using a more recent Illumina technology (see below). For whole genome sequencing of batch 1-3 we fragmented 1 μg of DNA of each individual using the 250 bp sonication protocol (Bioruptor, Diagenode, Liège, Belgium). The fragments were end-repaired and adaptor-ligated, including incorporation of sample index barcodes. The products were then purified and amplified (10 PCR cycles) to create the final libraries. The TruSeq DNA LT Sample Prep Kit v2 was used for all steps. After validation (Agilent 2,200 TapeStation), all libraries were quantified using the Peqlab KAPA Library Quantification Kit and the Applied Biosystems 7900HT Sequence Detection System. One library was loaded on two lanes of a Hiseq2000 sequencer and sequenced with a 2 × 100 bp v3 protocol.
The DNA quality of the Afghanistan mice (batch 4) proved problematic. Therefore, to recover high molecular weight genomic DNA, we ran a 0.7% agarose gel over night and extracted the genomic DNA from the gel using the Zymoclean Large Fragment DNA Recovery Kit (Zymo Research Europe, Freiburg im Breisgau, Germany). For whole genome sequencing we used the Nextera DNA library Prep Kit (Illumina) following manufacturer's instructions and 50 ng of genomic DNA as starting material. Each sample was run on Agilent Bioanalyzer using the Agilent DNA7500 kit to verify that the fragment sizes were in the 500 bp range. To calculate the final concentration for the sequencing run, the samples were measured with the Quant-iT dsDNA BR Assay Kit on a Nanodrop 3,300 fluorometer. The samples were paired-end sequenced (76 bp) independently on a single flow cell on a NextSeq 500 using the NextSeq 500 High Output v2 150 cycles chemistry.
RNA sequencing library preparation. We used the NEBNext Ultra RNA Library Prep Kit for Illumina with the Poly(A) mRNA Magnetic Isolation module to generate the RNAseq libraries. We used a fragmentation time of 15 min (yielding~180 bp fragments) and 15 PCR cycles for library enrichment. After validation (Agilent 2,200 TapeStation), all libraries were quantified using the Peqlab KAPA Library Quantification Kit and the Applied Biosystems 7900HT Sequence Detection System. Samples were pooled such that we generated about 12 million paired reads/sample, which were loaded on one lane each of a Hiseq2000 sequencer and sequenced with a 2 × 100 bp v3 protocol.

Data analysis
Mapping genomic reads and Single Nucleotide Polymorphism (SNP) calling. All sequencing reads (including those of the 10 previously published M. m. castaneus 51 genomes) were processed according to a single standardized pipeline, which is outlined in Fig. 2a and described in detail (including commands) in Supplementary Material Text 2. In brief, reads were mapped against the mouse mm10 genome reference sequence 54 (http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/mouse/) using bwa-mem 55 . The Picard tools software suite (http://broadinstitute.github.io/picard/) was used for sorting, marking and removing duplicates. Raw SNP and indel calls were obtained from the alignment files following precisely the GATK 56 'Best Practice' instructions on joint genotyping of all samples together. The raw .vcf files were subjected to the GATK VSQR SNP filtering step, which uses known variants as training data to predict whether a new variant is likely a true positive, or a false positive. As training data we used the file 'mgp.v5.merged.snps_all.dbSNP142.vcf' downloaded from ftp://ftp-mouse. sanger.ac.uk/current_snps/ 57 which was filtered for 'PASS' SNPs. In addition, we used very stringent hard filtering criteria on our own dataset, and included these SNPs as training sets as well (see details in Supplementary Material Text 2). Due to an absence of a reliable indel reference dataset we did not generate VSQR calls for indels. Thus, all indels called in the .vcf file should be considered 'raw'.
On our ftp website we provide a .vcf file with all SNPs and indels where we flag all SNPs within the 90% VSQR tranche as 'PASS' SNPs (this means that we accept all variants until we reach 90% of our known 'truth' set). This is a rather conservative filter on SNP quality, but users are free to perform their own filtering using their own training set and parameters for filtering on the raw .vcf file.
Transcriptome sequencing of M. m. domesticus populations. Transcriptome sequencing reads were processed according to the pipeline depicted in Fig. 2b. In short, they were trimmed according to quality values using Trimmomatic 58 , removing bases below Q20 and maintaining an average read quality above Q25. Pairs with one read below the quality thresholds were removed from the analyses. Quality-checked (QC) reads were mapped against the mouse mm10 genome reference using TopHat2 (ref. 59) and using the default settings for paired-end samples. The output alignments were sorted and indexed with samtools 60 . The sorted alignments were assigned to the version 82 of Ensemble Mouse gtf annotation 61 using featureCounts from the subreads suite 62 in paired-end (-p) exon mode (default). The gtf file contained only linear complete chromosomes (no scaffolds, no mitochondria). The unmapped pairs were counted with featureCounts from the 'ounmapped.bam>' file TopHat2 generates. Percentages are reported relative to the total number of QC reads. On average (across all tissues and all samples) 93% of the total number of QC reads could be mapped to the mm10 genome (range 86.5-96%, Table 2 (available  online only), Supplementary Table 3). This number also includes spliced reads, which Tophat2 detects. This number dropped to 57% (range 33-66%) for reads uniquely mapping to ENS 82 annotated features (i.e., exons, Supplementary Table 3).
Copy number variation (CNV) analyses. We used the sequencing read depth approach implemented in the CNVnator software 63 to predict CNV calls relative to the mouse mm10 reference assembly. We have previously experimentally confirmed that this is a reliable approach 38 . Optimal bin size for each individual was chosen such that the ratio of the average read depth signal to its standard deviation was between 4 and 5. Bin size ranged from 100-1,500 bp and was inversely proportional to genome coverage. Only linear complete chromosomes were considered. Calls intersecting annotated gaps in the reference genome were not considered. The CNV detection statistics are provided in Table 3 (available online only). Haploid copy numbers for each detected CNV, either per population or per individual, are included in the bed files for the UCSC browser tracks (available at the ftp site).
Visualization of data within and between populations. For the genomic data we set up UCSC 45 genome browser tracks for 10 kb windows of nucleotide diversity π, Tajima´s D and F ST, calculated using vcftools 64 , and CNV tracks. The 90% truth VSQR-filtered 'PASS' data were used for all vcftools calculations, allowing only bi-allelic SNPs and a maximum of 20% of missing data per SNP. The output tables were converted to bigWig format using the BigWig utilities 65 . For the CNV tracks we used bedtools 66 to intersect calls from all individuals belonging to the same population. Within each population average copy number across individuals was calculated for every given interval and transformed to log 2 values. The genomic browser tracks provided at the ftp site allow a visualization of differences between the populations and species for genomic regions of interest.
For RNAseq data, the coverage (number of reads at a given base pair) is expected to be proportional to the expression level of the gene from which the read originated. In order to compare RNAseq coverage (and thus gene expression) across individuals with slightly varying total numbers of mapped reads, we normalized the data by proportionally sub-sampling reads from the alignment file (*.bam) for the individuals with higher numbers of total mapped reads. Specifically, for each given tissue, we first To visualize the data as number of reads covering a particular location (i.e., base pair) in the genome, we converted normalized alignment files into bedgraph files, which were further compressed into bigWig format (available on the ftp site). The RNAseq based read coverage for each basepair in the genome can then be visualized in the UCSC (available as public session under 'wildmouse') or IGV browsers 46,47 . Figure 3 shows two examples of screen shots from IGV sessions. The first is a general overview across all tissues from a section of chromosome 10. The second shows only a single tissue (brain), but with read coverage information for each individual.

Data Records
The primary read files for the genome sequences are available at the European Nucleotide Archive   Table 1 (available online only). The transcriptome read files are available at ENA under project accession number PRJEB11897 (Data Citation 5). The samples and their sample designation are described in Table 2 (available online only).
The files with the mapped reads (bam), variant calling (vcf) and browser tracks (bigWig) for the genomes and transcriptomes, as well as the IGV session files for the transcriptomes are available at: http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/ where they can be accessed via ftp.

Technical Validation
We used the software angsd 67 and its -doDepth 1 command (with options -minMapQ 30 -minQ 20) to assess the quality of each DNA sequencing library with respect to good quality coverage (both mapping and base quality) of the genome. The angsd depth analysis is based on all sequenced and mapped bases, rather than only on called genotypes at variable sites. Thus, this is the most comprehensive way to assess coverage across the whole genome. The average per-base coverage for each sequenced genome is given in Table 1 (available online only) for autosomes, the X-chromosome and the Y chromosome separately. It was calculated as: P 60 i¼1 n i i genome size , where n i is number of bases sequenced at depth i and genome size is 2,395,908,738 for autosomes, 163,487,995 for the X chromosome and 88,124,698 for the Y chromosome.
The autosomal coverage is variable across individuals (both within and between sequencing batch) but should be high enough to obtain good quality SNP calls. For males, the X-chromosome coverage is expected to be half of the autosomal coverage, while for females, the coverage should be similar for X chromosome and autosomes. We can use this fact to obtain independent confirmation of the animal's sex recorded in the field. For all but one case, the genomic sex (based on X/autosome coverage ratio) matched the sex determined in the field. Individual AL42 was recorded as male in the field, but its genomic data clearly suggest it was a female. For juvenile wild mice it can sometimes be difficult to accurately determine their sex. The lower X-chromosomal coverage for males indicates that some caution should be taken when using called genotypes for this chromosome for population genetic inferences. Approaches that take the genotype likelihoods into account to estimate parameters may be better suited for the X-chromosomal data (i.e., ref. 67). The Y chromosome is extremely poorly covered (see below).
We also assessed the uniformity of coverage across the genome, aiming at identifying specific regions, where few or no reads could be mapped. We ran the angsd -doDepth 1 command for non-overlapping 100 kb windows recording the average (across individuals within subspecies) % bases covered at >10 reads/individual in 100 kb windows. As shown in Supplementary Fig. 4 using all M. m. domesticus individuals as an example, there are some regions in the genome where coverage is low or absent. This pattern was highly correlated in the other subspecies/species (data not shown), suggesting that features of the reference genome (possibly presence of repetitive elements or un-sequenced parts of the reference genome) limit mapping of reads in these regions. Regional variation in coverage is especially striking on the Y-chromosome, where coverage is limited to several short regions within the proximal 10 Mb of the chromosome. This region roughly corresponds to the male-specific short arm of the Y chromosome and its centromere.

Confirmation that mice are naturally inbred
The natural inbreeding status of wild caught mice is expected to be influenced by population history, as well as their tendency to form extended family structures with breeding among relatives 48,49 . The degree of inbreeding is expected to be higher for small populations and also for populations that recently colonized new habitats, a process which often involves a bottleneck. For the M. musculus mice included in this dataset, Iranian, Afghanistan and Indian mice are closest to the center of origin of this species and thus are expected to be least inbred, as such populations are expected to have been large and stable over time. On the other hand, M. m. helgolandicus inhabits a very small island and experienced a strong founder event during colonization 26 . For Mus spretus, we do not have a good expectation, as wild populations of this species have never been studied. It is generally assumed that M. spretus does not form extended family structures, as these mice are not human commensals but live in fields, hedges and boundaries to forests.
We used a combination of angsd 67 to calculate genotype likelihoods and ngsF 68 to calculate inbreeding coefficients for each individual based on randomly selected 1,000 autosomal 10 kb fragments (see Supplementary Material for scripts). As shown in Fig. 4, our expectations are mostly met, with ancestral populations from India and Iran showing the lowest estimated inbreeding coefficient, while Heligoland individuals are close to representing inbred lines. Surprisingly, the Mus spretus individuals from Spain are highly inbred, suggesting that they may have experienced a recent bottleneck.

Confirmation that animals cluster with their respective population
The aim of this analysis was to confirm that each sequenced individual is assigned to its respective population, based on its SNP genotypes. We used the software NgsAdmix 69 on the same randomly selected 1,000 autosomal 10 kb fragments as used in the previous analysis. As before, genotype likelihoods were generated by angsd. We ran the NgsAdmix software for K = 1 to K = 9 (Fig. 5). The likelihood of the   (with M. m. helgolandicus clustering with the German individuals, confirming the previous results based on microsatellites 26 ). For the M. m. musculus subspecies we find the individuals from the Czech Republic to split off from those from Afghanistan and Kazakhstan. The latter two populations seem to be more closely related. Generally, the genetic clustering analysis suggests that the populations are well defined and differentiated and that there are no recent immigrants from other areas among the sequenced individuals (note that the seeming admixture in K = 9 is an artifact of too high K).

Confirmation that VSQR 90 tranche filtering yields expected levels of polymorphism
We assessed the GATK VSQR 90% tranche PASS-filtered SNPs for levels of polymorphism (Watterson's θ 71 ) within each population and compared the estimates to several reference data sets, which have previously been generated using Sanger sequencing on smaller number of loci. The Indian population is especially informative, as it has been extensively sequenced 51,72 . For each chromosome and each population we determined the number of segregating sites using the software PopGenome 73 . Values were summarized over the autosomal genome and converted into Watterson's θ in % by dividing by a i 71 and the number of sites sequenced (Supplementary Table 4

Analysis of relatedness in the sample
We used the -relatedness2 option of vcftools to assess pairwise individual relatedness among all mice in the dataset, using the KING method 74 . This analysis is based on GATK called genotypes and the 90% tranche PASS-filtered SNPs. We restricted the dataset to only include autosomal SNPs, thinned to 1 SNP every 1 Mb. We also removed sites that had more than 20% missing data and only included bi-allelic markers in the analysis. As described in Table 1 Table 5). However, since we detected third degree relationship also among animals that were unequivocally caught far apart (e.g., SP39-SP68), we only consider first and second degree relatedness relevant here. No duplicate samples were detected (expected Phi = 0.5). Relatedness was only detected within populations, and was absent between them. No first or second-degree relatedness was found for the German M. m. domesticus, the Afghanistan M. m. musculus and the Indian M. m. castaneus populations. Most related animals were found in the populations from Iran and Kazakhstan. In the case of the Iranian population the increased relatedness within the sample can be explained by the fact that some breeding adults were used in multiple crosses (see Supplementary Table 2). The relatedness observed in the population from Kazakhstan is best explained by the fact that mice were collected in close proximity, rather than over a larger regional scale.
We can use the known breeding setup of the Iranian mice to confirm the inferred relatedness categories in this population. The KING method detected 2 first-degree relationships in the Iranian population. Male AH15 is indeed the father of JR5-F1C. However, JR-7F1C is the brother of the mother (i.e., uncle) of JR11 and thus a second-degree relative. The second-degree relative identified by the KING method in the Iranian sample is consistent with the breeding scheme, with JR11 and JR15 being half siblings (they have the same father). Two third degree relationships are incorrectly identified and should be second-degree relationships instead (i.e., JR5-F1C is the uncle of JR15 and AH15 is grandfather of JR15) and one third-degree relationship does not have any known breeding history confirming it. The KING method assumes Hardy-Weinberg equilibrium among SNPs with the same underlying allele frequencies. This assumption is most likely being violated in our dataset (given that wild mice are generally inbred, see above), and could explain some of the miss-assignments between categories.

Confirmation of known t-haplotype carriers and identifying t-haplotype carriers in the total dataset
The t-haplotype is complex set of 4 inversions, comprising a 20 cM (30-40 Mbp) region of the proximal third of chromosome 17 in house mice 75 . It is a selfish genetic element, which causes transmission ratio distortion, with heterozygous t-haplotype carriers predominantly (sometimes up to 99% of times) transmitting the t-haplotype carrying chromosome to their offspring. Homozygous individuals for the t-haplotype, however, die in utero. Despite their massive transmission advantage, t-haplotype carrying individuals are rare in natural populations of mice, but have been found in all recognized subspecies.
We have previously used two published primer pairs 76 to genotype individuals in our collections of mouse samples for presence/absence of the t-haplotype. Both primer pairs span t-haplotype diagnostic indels that can be analyzed on an agarose gel (the proximal locus Tcp1 spans a 175 bp indel and the distal locus Hba-4ps spans a 16 bp indel). Three individuals typed with those primers are identical to samples included in the whole genome sequencing study described here: male AL41, male CR16 and female H14. Of those, AL41 and H14 were found to be carriers (heterozygous) of the t-haplotype, while CR16 was wildtype. We use ENSEMBL Blast to determine the location of those primers in the mm10 reference sequence. All primers yielded unique hits. We then extracted all indels spanning the region between the primer locations for all samples from the .vcf-file, and searched for (combinations of) indels matching the sizes above. For the distal locus, we found a single 17 bp indel at position chr17:26,286,509 ('TACTACTATGCACTGAA'). For the proximal locus, we found one indel at position chr17:12,921,682 that was identified as 'GTTTTTTTTTTT'. Illumina sequencing is not capable to sequence through long homo-polymer stretches which is likely the reason why the identified indel is reported shorter than the expected 175 bp. However, it is the only indel >3 bp in the region and moreover, genotypes at this indel are in almost perfect linkage disequilibrium with genotypes at the distal locus (see Table 1 (available online only)), which is highly unexpected over a distance of 13.3 Mb. The two individuals that we previously found to be positive for the t-haplotype using the PCR primers are also heterozygous for the respective indels in the whole genome dataset, while CR16 was wildtype based on the whole genome sequence. Thus, indels at chr17:12,921,682 and chr17:26,286,509 were used to genotype the remaining individuals in the dataset for the presence/absence of t-haplotypes. All M. musculus individuals positive for the t-haplotype indels were heterozygous for that t-allele. Mus spretus, on the other hand, was homozygous for the proximal t-specific indel, which is consistent with a rather old origin of the t-allele (1-3 Myr 77 ), despite not having any obvious transmission ratio distortion properties in a species outside the M. musculus complex. t + /wt individuals were found in every population apart from the Iranian population and among the three individuals from Heligoland. T-haplotype carriers reached frequencies of 50% in two M. musculus populations (Afghanistan and Czech Republic). 37,5% of the French individuals carried the t-haplotype, as did 30% of Indian individuals, 25% of individuals from Kazakhstan and 12.5% of the German population. The frequency of t-haplotypes in our data is somewhat higher than reported previously (reviewed in ref. 78), however well below the theoretical expectations based on transmission ratio distortion (see ref. 79).

K-mer distributions to determine complexity of RNAseq libraries
We analyzed the k-mer (DNA sequence stretch of length k) frequency spectrum to identify potential problems with the RNAseq libraries, such as DNA contamination and overabundance (potentially due to PCR amplification) of particular sequence stretches. Ideally, libraries are complex, meaning they exhibit great diversity in unique sequences and repeated structures 80 and there should be no sign of DNA contamination among the RNA based reads generated. In total, we generated 224 RNAseq libraries from up to 10 tissues in 24 M. m. domesticus individuals. For each library, we ran the software jellyfish 81 on all forward reads with a k-mer length of k = 12. Using reverse reads yielded the same results (data not shown). For visualization in Supplementary Fig. 5 we randomly choose three individuals from each population and three tissues (brain, testis and liver). The complete set of all 224 RNAseq k-mer distributions are available on our ftp server. The k-mer distribution of the mouse mm10 DNA sequence (red in Supplementary Fig. 5) produces a characteristic line with 2 prominent humps when plotted on a log10-log10 scale. The k-mer distribution for annotated cDNAs in the mouse genome (ENSEMBL version 83 (ref. 61)) does not produce such humps and compared to the DNA mm10 sequence is more 'complex', i.e., shows proportionally more unique sequences (with low 12-mer occurrence, but high frequency, left side of plots in Supplementary Fig. 5). The k-mer distribution of the RNAseq forward reads (grey in Supplementary Fig. 5) falls in between the cDNA and genomic DNA profile. Since we generated our RNAseq libraries from poly-adenylated RNAs, comparing the RNAseq profile to annotated cDNA seems appropriate. Most notably, the RNAseq k-mer profile lacks the characteristic humps of the genomic DNA k-mer profile, suggesting that the RNA libraries are not contaminated with DNA. Moreover, the diversity of unique sequences matched the cDNA profile much better than the DNA profile. Very similar patterns have been observed for human RNAseq data, that have been deemed good quality (see Supplementary Fig. 6A in ref. 80).