Exome sequencing

DNA was extracted from spleen collected from female mice, when possible, using the Qiagen MasterPure^™ Complete DNA and RNA Purification Kit and protocol from Epicentre (Madison, WI). DNA was sheared using a Bioruptor UCD-200 with 7 rounds of sonication (7 minutes per round on high, 30s on 30s off) and genomic DNA libraries were constructed and individually barcoded using a previously described protocol designed to facilitate multiplexed exome capture (Rohland and Reich 2012). To reduce molecular interference during enrichment, we used truncated adaptors containing unique “internal” barcodes on the P5 end of genomic fragments (Rohland and Reich 2012). PCR primers were designed according to Rohland and Reich (2012).

In-solution sequence capture was performed using Nimblegen SeqCap EZ Mouse Exome probes as described in Nimblegen’s SeqCap EZ Library User’s Guide. Libraries were pooled equally to obtain 1 μg total DNA for each hybridization experiment. Libraries were then enriched using two separate capture reactions with eight libraries each, including blocking oligonucleotides specific to our custom adapters (Rohland and Reich 2012) and mouse COT-1 DNA (Invitrogen) to reduce non-specific hybridization. The capture reactions were hybridized for 68 hours at 47°C in an Eppendorf Mastercycler Pro, and then washed, eluted, and PCR-enriched. Capture enrichment success was verified using qPCR analysis of three targeted regions on pre- and post-capture library pools. Sequencing was performed using 76 bp paired-end reads on the Illumina Hi-seq 2500 platform provided by the Epigenome Center at the University of Southern California.

Illumina reads were mapped to different pseudoreference genomes dependent upon their species of origin (Sarver et al. 2017). A pseudoreference contains the backbone of the mm10 reference mouse genome (strain=C57BL/6J), but allelic states taken from representative strains of M. m. musculus, M. m. domesticus, M. m. castaneus or M. spretus are inserted into the mm10 reference genome. This approach leaves all called SNPs on a common coordinate system (mm10). All pseudoreferences were taken from Sarver et al. (2017). Preliminary analyses demonstrated that the four strains of unknown origin (BID, KAK, MPR, and TEH) fell within the M. m. musculus clade, and so were mapped to the M. m. musculus pseudoreference. The advantage of this technique is that species-specific variation can be incorporated into the reference, thus improving mapping accuracy and recovery, while preserving the coordinates of the mm10 build so that genome annotation can be used. Sequences were mapped with BWA MEM (v0.7.9a, Li and Durbin 2009), using default mapping parameters. Alignments to their respective pseudoreferences were used to identify variants using GATK HaplotypeCaller (McKenna et al. 2010), following PCR duplicate removal and indel realignment. Variants from dbSNP (version 142) were used as the training set during SNP recalibration and subjected to standard hard filtering parameters according to the GATK Best Practices recommendation (Auwera et al. 2013; DePristo et al. 2011): MQ > 56, QD > 24, FS < 12, MQRankSum < 8, ReadPosRankSum < 3, DP > 3. To assess confidence, we compared SNPs identified in the Montpellier strains that also overlapped with dbSNP and quantified how many calls agreed with known allelic states at those sites. In addition to novel SNPs, we identified 152,342 insertion/deletion mutations (Supplementary Table 4). SNPs were classified by functional categories using GRCm38.73 as annotated by EMBL-EBI Ensembl (www.ensembl.org).

Genetic relationships

We placed the 26 Montpellier exomes in the context of an additional 36 inbred strains commonly used in genetic research and compiled in dbSNP (v142) (Sherry et al. 1999; Sherry et al. 2001). Most of these are “classical inbred strains”, which share a common ancestry predominantly originating from a limited M. m. domesticus stock, with some human-mediated contribution from other species (Keane et al. 2011; Yang et al. 2007; Yang et al. 2011). Seven strains from dbSNP are considered wild-derived strains: three M. m. domesticus (LEWES/EiJ, WSB/EiJ, and ZALENDE/EiJ), one M. m. musculus (PWD/PhJ), one M. m. molossinus (MOLF/EiJ), one M. m. castaneus (CAST/EiJ), and one M. spretus (SPRET/EiJ). M. m. molossinus is a subspecies that originated through natural hybridization between M. m. musculus and M. m. castaneus (Yonekawa et al. 1986; Yonekawa et al. 1988).

In addition to placing the Montpellier strains in the context of known variation among other inbred strains, we repeated our analyses after also combining our data with the genetic variants called by Harr et al. (2016), specifically their file named AllMouse.vcf_90_recalibrated_snps_raw_indels_reheader_PopSorted.PASS.vcf available from http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/. We only used sites covered in our Montpellier strains and dbSNP. Harr et al. (2016) performed whole genome sequencing from 27 M. m. domesticus, 22 M. m. musculus, 10 M. m. castaneus, and 8 M. spretus, to an average depth of 20.9X autosomal coverage. Although the Harr et al. (2016) data are from wild caught animals, they provide valuable context through which to view genetic relationships of the inbred strains (Montpellier + dbSNP strains).

SNPs that were found in our Montpellier panel as well as dbSNP were used to generate genealogies with the SNPHYLO program (Lee et al. 2014), which uses the maximum likelihood framework of the DNAML program in the PHYLIP package (Felsenstein 1993) to construct trees. 1000 bootstrap analyses were performed using all variants with the PHANGORN package (Schliep 2011) in R (www.r-project.org). Trees were visualized using the APE package (Paradis 2012; Paradis et al. 2004) in R. Resulting trees will be strongly affected by introgression known to occur in the history of inbreeding and therefore serve only as a rough approximation of genetic relationships.

We assessed genetic structure in our sample using STRUCTURE (Falush et al. 2007; Hubisz et al. 2009; Pritchard et al. 2000) and Principal Components Analyses (PCA). STRUCTURE was run under the admixture and correlated allele frequency model with ten independent runs of 10,000 burn-in MCMC iterations followed by 50,000 iterations for 2 to 8 clusters (k=2 to 8). Results were inspected with STRUCTURE HARVESTER (Earl 2012; Evanno et al. 2005). PCA was performed using the SNPRELATE package (Zheng et al. 2012). Unlike the above genealogical anayses, STRUCTURE and PCA require sites called in all 62 strains (i.e., no missing data across the 26 Montpellier exomes plus the 36 additional dbSNP strains). Multiallelic sites were excluded. To reduce the size of the dataset, we chose sites that were at least 100 kb apart, roughly the extent of linkage disequilibrium in wild mice (Laurie et al. 2007), resulting in 26,991 sites used in STRUCTURE and PCA analyses.

Protein coding variation

In addition to classifying SNPs into basic categories such as nonsynonymous and synonymous, we estimated codon usage bias across strains, using ENCPRIME (Novembre 2002), which tests the null hypothesis that codons within an amino acid class are used at equal frequency, after accounting for background base composition. Estimates of codon usage theoretically range from 20 (every amino acid coded by a single codon, representing maximal bias) to 61 (each amino acid coded by each of its synonymous codons at equal frequency, representing minimal bias). Background base composition was estimated from flanking non-exonic sequence, which is inevitably sequenced even when performing exome enrichment. Because base composition varies across the genome, we analyzed codon usage bias in 5 Mb windows. Lastly, we annotated early stop codons, which could represent alternatives to traditional gene knockouts if they disrupt normal gene function. To infer their functional impact, we determined the percentage of the protein truncated by the early stop codons, and whether they occurred on constitutive or facultative exons.

Gene flow between lineages and strains

Previous studies have shown that many laboratory strains have a mosaic genome that contains genetic material from multiple species (Dai et al. 2005; Ferris et al. 1982; Frazer et al. 2007; Ideraabdullah et al. 2004; Nagamine et al. 1992; Tucker et al. 1992; Yalcin et al. 2004; Yang et al. 2007; Yang et al. 2011). We tested for interspecific introgression using the ABBA-BABA test (Green et al. 2010). In this test, “A” and “B” indicate the distribution of a biallelic state across a rooted four-taxon genealogy. We tested the genealogy: (((M. m. domesticus strain #1, M. m. domesticus strain #2), M. m. musculus), M. spretus) represented by (((DOM*, LEWES/EiJ), CzechII/EiJ), SPRET/EiJ), where DOM* represents each Montpellier strain of M. m. domesticus tested individually. LEWES/EiJ and CzechII/EiJ were chosen to represent “pure” M. m. domesticus and M. m. musculus, respectively, as neither shows high levels of introgression (Yang et al. 2011, and unpublished data). We have sequenced the genome of CzechII/EiJ as part of unrelated work, only using variant calls here for the ABBA-BABA implementation. These results should be treated with caution as M. m. musculus and M. m. domesticus are closely related, can still interbreed in nature, and their genomes are not fully differentiated. An excess of either ABBA or BABA sites is best explained by introgression, in our case from M. m. musculus (represented by CzechII/EiJ) into M. m. domesticus.

We divided the genome into 489 non-overlapping 5 Mb windows, discarding any that had fewer than 20 ABBA-BABA sites with genotypes called with a minimum genotype quality (GQ) of 10. GQ is the Phred-scaled confidence in a called genotype and is strongly correlated to the number of reads that map to a particular site, the quality of the base calls from those reads, and sequencing error. For any DOM* strain that contained at least 20 windows with at least 20 ABBA-BABA sites, the null hypothesis of no introgression was evaluated using the block Jackknife procedure as implemented in the script jackKnife.R from the ANGSD package (Korneliussen et al. 2014).

Several recent studies have raised the possibility of gene flow between some wild-derived strains that may have unintentionally occurred after their establishment in the laboratory (Yang et al. 2007; Yang et al. 2009; Yang et al. 2011). To test for more recent introgression, we inspected the distribution of pairwise divergence between every possible pair of Montpellier strains having shared the same lab environment. Recent introgression should result in long stretches of the genome that are identical by descent. We used the B-SMUCE segmentation algorithm of Futschik et al. (2014) implemented in the smuceR function of the R package stepR with default settings. This multiscale segmentation algorithm looks for compositionally homogeneous segments. It estimates the best fit to a series of a minimal number of steps and provides estimates for the number of segments and their boundaries at the same time. Originally employed to find regional variation in GC (1) vs. AT (0) content, we applied this algorithm to a simple literal distance whereby two homologous SNPs yield 0 if they were identical or identically heterozygous, 1 if different and 0.5 if one of them was heterozygous. The program was run with default options to estimate the chromosomal distribution of haplotype sharing between every pair of strains.

Genetic relationships

The four M. spretus strains (Fig. 1, left panel, strains labeled in black) formed a distinct group that was distantly related to the other strains, consistent with previous phylogenetic hypotheses (Lundrigan et al. 2002; Sarver et al. 2017). M. m. castaneus and M. m. musculus formed a group (Fig. 1, left panel, strains labeled in green), that then grouped with M. m. domesticus strains (Fig. 1, left panel, strains labeled in red), also consistent with previous phylogenetic inference (Keane et al. 2011; Phifer-Rixey et al. 2012; White et al. 2009). M. m. domesticus, M. m. musculus, and M. m. castaneus are closely related lineages that diverged less than 350 thousand years ago (Boursot et al. 1996; Geraldes et al. 2008; Salcedo et al. 2007; She et al. 1990; Suzuki et al. 2004). Many genomic regions are still unsorted and lack species-specific substitutions (Geraldes et al. 2011; Salcedo et al. 2007). MOLF/EiJ falls within the other M. m. musculus strains, suggesting that even though it is a natural hybrid species between M. m. castaneus and M. m. musculus (Yonekawa et al. 1986; Yonekawa et al. 1988), a majority of its genome is derived from M. m. musculus.

Open in a separate window

Figure 1

Genealogical relationships among 26 Montpellier strains and 36 strains from dbSNP, determined by A) maximum likelihood tree, or B) STRUCTURE analysis considering k=5. Strain names on a black box indicate novel Montpellier exomes sequenced in this study. Nodes labeled with small black circles were supported with at least 95% bootstrap support.

The classical inbred strains (Fig. 1, left panel, strains labeled in blue) fell within the M. m. domesticus group, consistent with previous studies suggesting they are mostly derived from this subspecies (Yang et al. 2007; Yang et al. 2009; Yang et al. 2011). Interestingly, the classical inbred strains formed their own group within wild-derived M. m. domesticus, and even other wild-derived inbred strains from The Jackson Laboratory, like LEWES/EiJ and WSB/EiJ, fall outside of this group. The distinct grouping of the classical inbred strains reinforces what is known about their history, that they derived from a small set of founders during the early mouse fancy trade (Beck et al. 2000; Frazer et al. 2007; Moriwaki 1994; Morse 1978; Morse 2007; Silver 1995; Wade and Daly 2005; Wade et al. 2002). These general groups remained unchanged after repeating the analysis with wild-caught M. m. domesticus, M. m. castaneus, M. m. musculus, and M. spretus of Harr et al. (2016) (Supplementary Figure 1).

The four strains that were not assigned to a subspecies a priori (BID, KAK, MPR, and TEH) all fell within the M. m. musculus clade, but in a basal position within that group (Fig. 1, left panel). These strains originated from individuals collected near the center of mouse diversity in the Middle East (Pakistan and Iran) (Hardouin et al. 2015). Some of these populations may constitute independent (sub-) species that have yet to be described (Rajabi-Maham et al. 2012), or they may have captured ancestral polymorphism or secondary admixture that occurred in the field.

STRUCTURE analyses largely supported the genealogical tree, with log-likelihood values reaching a plateau at K=5 groups. The Evanno method was used to characterize the clustering, with delta K peaking at K=3 and falling to 0 at K=5. The two M. m. castaneus strains (dark green, Fig. 1, right panel) shared a large proportion of genetic variation with M. m. musculus (light green, Fig. 1, right panel). Interestingly, the allele frequencies among wild-derived M. m. domesticus strains (red, Fig. 1, right panel) differed from those of classical inbred strains (blue, Fig. 1, right panel). There is no evidence for genetic structure within classical inbred strains, which capture approximately half of the structure seen in wild derived M. m. domesticus (blue in Fig. 1, right panel). At K=3, M. spretus, M. m. castaneus + M. m. musculus, and M. m. domesticus form distinct clusters. At K=4, M. m. castaneus separates from M. m. musculus. At K=5, the classical inbreds separate from the wild-derived M. m. domesticus. K=6 does not separate out any samples from the rest.

PCA also supported the groupings seen in the genealogy (Supplementary Fig. 2). The first principal component explained 45.2% of the variation and separated the four M. spretus strains from the rest of the panel. The second principle component explained 16.7% of the genetic variation and separated the M. m. musculus and M. m. castaneus strains from the rest of the panel.

Protein coding variation

M. spretus showed low codon usage biased (corresponding to high ENCprime estimates), M. m. castaneus and M. m. musculus strains had high codon usage bias, and M. m. domesticus was intermediate (Fig. 2). Natural selection shapes codon usage bias across mammals, but shows no consistent correlation to proxies of effective population size (Kessler and Dean 2014). The effective population sizes of M. m. castaneus, M. m. domesticus, and M. m. musculus have been genetically estimated at 220K, 100K, and 60K, respectively (Geraldes et al. 2008; Geraldes et al. 2011; Phifer-Rixey et al. 2012). We might hypothesize that M. spretus has a smaller effective population size than M. musculus subspecies because it occupies a smaller geographic area and maintains lower density than M. musculus subspecies (Boursot et al. 1985; Britton and Thaler 1978; Dejager et al. 2009; Orsini et al. 1982). As predicted by a model of weak selection, M. m. castaneus has the strongest, and M. spretus the weakest, codon usage bias. However, the two intermediate taxa do not follow the predictions of weak selection. Thus, there is no consistent relationship between inferred effective population size and codon usage bias estimated from our exomes, as observed by Kessler and Dean (2014).

Open in a separate window

Figure 2

Codon usage bias among Montpellier strains.

Across all 26 exomes of wild-derived inbred strains, we identified 262 genes with at least one early stop codon segregating in at least one isoform, many of which segregate in multiple strains (Supplementary Table 5). Of these, 193 truncate less than 10% of the protein and therefore may not have a strong functional impact (Supplementary Fig. 3). However, 117 genes segregate an early stop codon that truncates more than 50% of its wild type length (Supplementary Fig. 3), and could represent a novel source of effective knockouts for future functional studies. Surprisingly, the length of wild type protein truncated by early stop codons did not differ between early stop codons that occurred in constitutively versus facultatively spliced exons (t=0.98, df=260, P=0.33) (Supplementary Fig. 4), suggesting that both types of stop codons have similar functional impacts. In fact, more stop codons affected constitutively (N=143) vs. facultatively (N=119) expressed exons.

Here, we highlight one gene, Bard1, with an early stop codon that truncates 50.6% of the wild type protein in two Montpellier strains, BIK and DDO (Supplementary Table 5). BARD1 forms a heterodimer with BRCA1, which together plays an important role in chromosomal stability and tumor suppression. A truncated BARD1 protein results in defective homologous DNA repair (Westermark et al. 2003), and knockouts for Bard1 die at an early embryonic stage (McCarthy et al. 2003). Whether or not the two Montpellier strains with truncated BARD1 have altered chromosomal stability represents one of many potentially interesting follow-up studies.

Gene flow between lineages and strains

Eight wild-derived strains of M. m. domesticus (BIK, DCP, DDO, DEB, DIK, DJO, DOT, WLA) showed a significant ABBA-BABA result (|Z-score| > 3, Supplementary Table 6), suggesting past introgression from M. m. musculus (represented by CzechII/EiJ). Four strains (22MO, BZO, DMZ, DCA) did not show a significant ABBA-BABA result. Among these 12 strains, we analyzed an average of 2,961.5 ABBA-BABA sites from an average of 79.3 windows. The remaining two strains (DGA and WMP) lacked enough data to apply an ABBA-BABA test (fewer than 20 windows with at least 20 ABBA-BABA sites).

Using the B-SMUCE segmentation algorithm (Futschik et al. 2014), two M. m. domesticus strains DJO and 22MO show “typical” segment lengths of estimated genetic distances between two independently derived strains (Supplementary Fig. 5, Supplementary Table 7). In contrast, the two M. m. domesticus strains DEB and BZO show a marked shift towards long segments with no genetic distance (Supplementary Fig. 5), a possible signature of recent introgression followed by limited recombination and inbreeding. However, the fact that all regions have pairwise genetic distances greater than zero (Supplementary Fig. 5) argue against recent introgression occurring after laboratory establishment.

We thank Charlie Nicolet and Selene Tyndale from the Epigenome Center at USC. Brent Young, Rachel Mangels, and Lorraine Provencio helped with molecular work. Matt Salomon and Rob Williams gave many helpful suggestions. Jean-Jacques Duquesne maintained the wild mouse repository in Montpellier. Funding was provided by the National Institutes of Health Grant #GM098536 (MDD), National Science Foundation Grant #1146525 (MDD), the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the National Institutes of Health Grant #HD073439 (JMG), and the University of Montana Genomics Core, supported by a grant from the M.J. Murdock Charitable Trust.