Hybrid dysfunction, a common feature of reproductive barriers between species, is often caused by negative epistasis between loci (“Dobzhansky-Muller incompatibilities”). The nature and complexity of hybrid incompatibilities remain poorly understood because identifying interacting loci that affect complex phenotypes is difficult. With subspecies in the early stages of speciation, an array of genetic tools, and detailed knowledge of reproductive biology, house mice (Mus musculus) provide a model system for dissecting hybrid incompatibilities. Male hybrids between M. musculus subspecies often show reduced fertility. Previous studies identified loci and several X chromosome-autosome interactions that contribute to sterility. To characterize the genetic basis of hybrid sterility in detail, we used a systems genetics approach, integrating mapping of gene expression traits with sterility phenotypes and QTL. We measured genome-wide testis expression in 305 male F2s from a cross between wild-derived inbred strains of M. musculus musculus and M. m. domesticus. We identified several thousand cis- and trans-acting QTL contributing to expression variation (eQTL). Many trans eQTL cluster into eleven ‘hotspots, ’ seven of which co-localize with QTL for sterility phenotypes identified in the cross. The number and clustering of trans eQTL—but not cis eQTL—were substantially lower when mapping was restricted to a ‘fertile’ subset of mice, providing evidence that trans eQTL hotspots are related to sterility. Functional annotation of transcripts with eQTL provides insights into the biological processes disrupted by sterility loci and guides prioritization of candidate genes. Using a conditional mapping approach, we identified eQTL dependent on interactions between loci, revealing a complex system of epistasis. Our results illuminate established patterns, including the role of the X chromosome in hybrid sterility. The integrated mapping approach we employed is applicable in a broad range of organisms and we advocate for widespread adoption of a network-centered approach in speciation genetics.
To understand patterns of biodiversity, it is essential to characterize the processes by which new species arise and are maintained in nature, including ecological specialization, population differentiation and reproductive isolation. Genetic dissection of reproductive isolation has proven to be an especially powerful strategy for revealing mechanisms of speciation. Many genomic regions and even specific genes that contribute to hybrid defects have been identified by genetic mapping in recombinant populations [1]–[7]. Divergence in gene regulation is expected to contribute to reproductive isolation between nascent species, and studies with F1 hybrids support this prediction [8]–[13]. Importantly, these two approaches – genetic mapping and measurement of genome-wide expression patterns in hybrids – have yet to be combined directly in the context of speciation. Hybrid sterility and hybrid inviability frequently result from negative epistasis between mutations at interacting genes [14]–[16]. This “Dobzhansky-Muller model” predicts that disruptions in gene networks should be common in hybrids. By integrating organismal phenotypes and genotypes with gene expression patterns, this prediction can be tested. Despite the identification of hybrid incompatibility genes in several species and the prevalence of the Dobzhansky-Muller model, the nature and complexity of hybrid incompatibility networks remains poorly understood. Do hybrid incompatibilities generally involve two loci or are higher order interactions common? Are incompatibilities independent or do they share some common loci? Is the genetic architecture of hybrid defects similar among taxa? Known incompatibility genes have provided the first hints about these questions, particularly in Drosophila [6], yet too few genes and taxa are represented to determine whether there are generalities underlying the speciation process. A network perspective should provide insights into the genetics of reproductive isolation that are difficult to obtain using a gene-by-gene approach. The house mouse (Mus musculus) is an excellent model for investigating speciation from a network perspective. Genomic resources are abundant for the house mouse, and reproductive processes are well characterized because the mouse is the premier model for fertility research in humans [17]. House mouse subspecies are in the early stages of speciation, showing significant but incomplete reproductive isolation. Evidence for hybrid male sterility in laboratory crosses [5], [18]–[22] and in natural zones of hybridization [23], [24] suggests it is a primary isolating barrier between these nascent species. Studies of sterility in F1 hybrids between Mus musculus domesticus and Mus musculus musculus (subsequently referred to as domesticus and musculus) revealed an important role for the X chromosome and identified several contributing autosomal loci [4], [5], [25], [26]. One of these loci is Prdm9, a histone methyltransferase [27]. Hybrids with some alleles of Prdm9 from domesticus show pachytene arrest of meiosis. The effects of sterile Prdm9 alleles appear to be due to mutations in the protein-coding sequence and there is evidence for downstream regulatory effects, but the incompatibility network involving Prdm9 has not been revealed. Genetic mapping of sterility phenotypes in F2 hybrids between M. m. domesticus and M. m. musculus recently identified an additional set of autosomal loci, which are predominantly recessive and thus contribute to sterility in second generation and subsequent hybrids. Genetic architectures of F2 sterility traits are complex, involving a moderate number of loci with a range of phenotypic effect sizes [1]. Genome-wide studies of gene expression in testis of F1 hybrids provide evidence that sterility is associated with disrupted expression [9], [10]. Like sterility phenotypes, expression patterns in hybrids depend on the origins of parental strains, and the direction of the cross. In many cases, testis expression in hybrids is intermediate between parental strains [9]–[11]. However, extensive misexpression (expression outside the range observed in parental strains) has been documented in a few crosses. Comparison of testis gene expression patterns between reciprocal F1 musculus-domesticus hybrids showed that many X-linked genes are overexpressed in sterile but not in fertile F1s [10]. To our knowledge, gene expression patterns in testes from F2 and later generation hybrids have not been described. Here, we integrate analysis of genome-wide expression in testis from F2 musculus-domesticus hybrids with results from a previous study mapping sterility phenotypes in the same individuals [1]. We show that sterility is associated with large-scale alterations in gene expression in F1s and F2s, and we identify quantitative trait loci (QTL) that cause X chromosome-wide overexpression in hybrids. We report expression quantitative trait loci (eQTL) for a large number of transcripts. We compare the locations of eQTL with sterility QTL, and identify disrupted processes during spermatogenesis based on affected networks. Using a conditional mapping approach, we pinpoint genetic interactions affecting expression. We highlight candidate pathways, processes, and interactions for several loci, which provide insight into the mechanisms underlying their contributions to sterility.
We measured levels of misexpression in F1 and F2 hybrids to identify major alterations in gene expression pattern associated with sterility in M. m. domesticus (WSB/EiJ; hereafter domesticusWSB) - M. m. musculus (PWD/PhJ; hereafter musculusPWD) hybrids. Sterility is asymmetric in these crosses: F1 males with musculusPWD mothers (hereafter MxD F1s) are almost always completely sterile whereas F1s with domesticusWSB mothers (hereafter DxM F1s) are fertile [1]. MxD F1 males showed significant differences from both parents for all reproductive traits measured. By contrast, all traits in DxM F1s (except seminiferous tubule area) were within the range observed in the parental lines. Trait measurements in MxD F1s and DxM F1s provide ‘fertile’ and ‘sterile’ examples that are useful for assessing trait distributions in F2s. Misexpression was markedly higher in testis of MxD F1s (18. 8% transcripts; Fig. 1A) than in DxM F1s (1. 6%). In both F1s, levels of misexpression were higher for X-linked transcripts than autosomal transcripts. On the X chromosome, the number of overexpressed transcripts in MxD F1s was much higher than the number of underexpressed transcripts (25. 9% over, 4. 4% under). The level of underexpression was higher on autosomes, but the difference between levels of over- and underexpression was smaller (7. 1% over, 11. 3% under). These results are consistent with previously reported differences in expression patterns between sterile and fertile F1s [10]. Misexpression in F2s varied from 0. 9–39. 0% transcripts (median 2. 1%; Fig. 1A), encompassing the levels observed in fertile and sterile F1s. There was substantial overlap between transcripts misexpressed in MxD F1s and in >5% of F2s (Fig. 1B) yet a large proportion of transcripts were misexpressed only in F1s or F2s. The relatively continuous distribution of misexpression in F2s and lack of recapitulation of the full F1 misexpression pattern indicates multiple genetic factors contribute to misexpression. Misexpression unique to F2s suggests some contributing loci act recessively. A large proportion of X-linked transcripts were negatively correlated with testis weight (lower testis weight = higher expression) – opposite of the pattern for autosomal transcripts, a majority of which was positively correlated with testis weight (Fig. 1C). This result suggests that – as in sterile F1s – the X may be broadly overexpressed in sterile F2s. To determine whether the level of misexpression was consistent throughout spermatogenesis, we compared patterns of expression in F1 and F2 hybrids among genes identified as specific/enriched to different spermatogenic cell types in previous studies [28]. Autosomal transcripts expressed in meiotic and post-meiotic cells are underexpressed in sterile MxD F1s, and transcripts specific to somatic and mitotic cells are overexpressed (Table S1). This pattern is consistent with reduced spermatogenesis, as expected based on sterility phenotypes. The X chromosome is transcriptionally silenced during meiosis (meiotic sex chromosome inactivation MSCI; [29], [30]), and thus lacks transcripts associated with meiotic cells. X-linked transcripts associated with other testis cell types showed patterns consistent with autosomal transcripts; somatic and mitotic transcripts tended to be overexpressed and the few underexpressed transcripts were predominantly postmeiotic. Misexpression patterns across spermatogenic cell types in F2 hybrids were consistent with patterns in sterile F1s. Next, we investigated the genetic basis of gene expression variation in individual transcripts. We identified 16,705–36,753 eQTL, depending on the significance criterion (Table 1). We used a permissive threshold, based on permutation of a single transcript, for downstream analyses because our goal was to identify genome-scale patterns. It is important to note that the false-positive rate among individual eQTL identified using this criterion is high, particularly for trans eQTL. The genomic positions of the eQTL and the affected transcripts are shown in Figure 2. eQTL located near the quantitative trait transcript (QTT) comprise the prominent diagonal stripe, a pattern typical of eQTL studies [32]–[34]. These proximal eQTL are likely to be cis regulatory elements [33], [35]. We refer to proximal eQTL as cis eQTL for convenience, although it is possible that they might not act solely in cis (by regulating alleles only if they are on the same DNA strand). We classified eQTL with peaks within 5 cM of the transcript (probe) position as cis eQTL and eQTL located on a different chromosome from the transcript as trans. We ignored eQTL>5 cM on the same chromosome, because this class might include long-distance cis eQTL in addition to trans eQTL. We identified cis eQTL for 60% of transcripts (14,807; Table 1) and at least one trans eQTL for 56. 7% (13,997) transcripts. The number of trans eQTL identified per transcript ranged from one (8,092; 32. 8% transcripts) to seven (3; 0. 01% transcripts). We next examined eQTL dominance and effect size. Most cis eQTL (93. 8%; Fig. S1A) were additive (mean for heterozygotes is intermediate and >2 standard errors from both homozygous means – see Methods). In contrast, a substantial proportion of trans eQTL were dominant (37. 1%), underdominant (9. 2%), or overdominant (8. 6%). Curiously, musculusPWD alleles were more likely to be dominant among cis (473/859; 55. 1%) and trans eQTL (2,850/4,580; 62. 2%). We cannot think of an experimental or biological explanation for this bias. The two categories of eQTL differed in effect size (Fig. S1B). The difference in expression level between genotype classes was larger on average for cis eQTL than for trans eQTL (t = 72. 3 (d. f. = 15931), P<2. 2×10−16). The difference in effect size is also apparent when comparing the peak LOD scores of cis (mean = 25. 05) and trans eQTL (mean = 5. 94). We tested for clustering of trans eQTL, which is commonly observed in eQTL analyses [36]–[38]. Some of these ‘trans hotspots’ are visible as vertical bands in the eQTL heatmap (Fig. 2). We identified 12 genomic regions significantly enriched for trans eQTL using a sliding window analysis (P<0. 05, permutation test; Table 2). Two adjacent hotspots on chromosome 10 were combined for simplicity in downstream analyses. The most striking pattern was observed for the X chromosome: most of the X was significantly enriched for trans eQTL and 8,286 autosomal transcripts (34. 6%) had eQTL mapped to the proximal X hotspot (0–42 cM). We discuss the massive effect of the X on gene expression in detail below, and relate this pattern to the known importance of the X in hybrid male sterility. The genomic distribution of eQTL we identified, as well as differences in dominance and effect sizes between cis and trans eQTL, are broadly consistent with patterns previously described in eQTL studies performed in a variety of (non-hybrid) organisms (e. g. humans: [37], [39]; C. elegans: [36], [40], [41]; Arabidopsis: [32], mice: [42]. This consistency indicates that misexpression and differences in expression level due to altered cell-composition associated with sterility phenotypes were not so severe that they obscured quantitative expression differences between musculusPWD and domesticusWSB. The Dobzhansky-Muller model predicts that each hybrid sterility locus will have one or more interaction partners. Mapping of genetic interactions generally requires sample sizes larger than the 305 F2s analyzed here. To increase power, we treated trans eQTL hotspots as candidate hybrid sterility loci and searched for interactions involving them. We performed conditional mapping of eQTL, using genotypes at candidate loci one at a time as covariates. Genotype covariates included the marker closest to the peak of each of the nine autosomal trans eQTL hotspots, and five markers in the X chromosome trans hotspots (Table 5). For each covariate, mapping was performed twice, including an additive effect or both an additive and interactive effect; eQTL from the full model that showed a significant increase in LOD score over the additive model were classified as significant interaction eQTL. Clustering of interaction eQTL identified by conditional mapping was even more pronounced than clustering of trans eQTL in the initial (no covariate) eQTL analysis (Fig. S3). We identified ‘interaction hotspots’ using significance thresholds from permutation for each genotype covariate. Integrating results from the conditional mapping analyses reveals a complex epistatic network showing several general patterns (Fig. 4). The large number of interactions involving the X is consistent with its substantive effect on expression pattern and sterility phenotypes. There are many interactions between loci in trans hotspots, and between trans hotspots and sterility QTL, suggesting that some incompatibilities contribute to multiple phenotypes. Overall, a large proportion of interactions are associated with sterility loci. It is important to note that many interactions may be associated with variation in gene expression unrelated to hybrid sterility. The interactions we identified include X-autosome pairs previously associated with hybrid sterility. We identified interaction hotspots in the proximal region of chromosome 17, which encompasses Prdm9, from conditional mapping using all X-linked genotype covariates; conversely, mapping conditional on Chr17@13 cM identified a hotspot on the proximal X (Fig. S4). Previous mapping of sterility phenotypes conditional on X genotypes revealed interactions between the X and six autosomal regions on four chromosomes (3,5, 7,10), contributing to five sperm morphology phenotypes [1]. We found interaction hotspots involving at least one X-linked covariate overlapping each of these autosomal regions. Each trans hotspot identified in the original analysis overlapped at least one interaction hotspot mapped with an autosomal covariate, indicating autosome-autosome interactions contribute substantially to expression variation. All of these interactions are novel. Interactions between regions with sterile alleles from the same subspecies are prevalent (Fig. 5), suggesting incompatibilities involving more than two loci are common. Conditional mapping revealed additional associations between gene expression variation and sterility. Some sterility QTL that did not overlap a trans hotspot identified in the original analysis showed evidence for interaction with one or more hotspot regions (Fig. S4). We also found interactions with sterility QTL for each of the trans hotspots that do not overlap sterility QTL. The relative contribution of loci to expression variation with detectable marginal effects versus eQTL identified only when incorporating interactions varied (Table 5). The structure of the interaction network provides additional support for the important roles of chromosomes X and 17, the major players in F1 sterility (Figs. 4; 5). By contrast, the chromosome 6 region plays a prominent role in the interaction network (Fig. 5), which was unanticipated on the basis of relatively modest enrichment of eQTL in the trans hotspot and the lack of sterility phenotype QTL on chromosome 6. We identified several novel loci that interact with multiple trans hotspots but did not have previous evidence for involvement in sterility (Fig. S4). Regions on chromosomes 7 (50–52 cM; 122. 63–125. 77 Mb), 13 (32–36 cM; 68. 47–75. 96 Mb), 14 (40–44 cM; 87. 59–97. 00 Mb) and 16 (0–4 cM; 11. 20–20. 02 Mb) had overlapping interaction hotspots identified by mapping with genotype covariates from trans hotspots on at least three chromosomes. These results indicate that some loci in the interaction network have marginal effects undetectable using single-QTL models and permutation thresholds.
The importance of evolutionary changes in transcriptional regulation for adaptation has long been recognized [e. g. 50]–[53]. Recent studies of gene expression in hybrids suggest regulatory evolution may also be an important cause of reproductive isolation between diverging populations. Misexpression has been reported in hybrids from many animal and plant taxa including Drosophila [8], [12], [54], mice [9]–[11], [55], African clawed frogs [13], [56], whitefish [57], copepods [58], maize [59], ragwort [60] and Arabidopsis [61]. Furthermore, several known hybrid incompatibility genes affect transcription of other genes, including OdsH [12] and the mouse sterility gene Prdm9 [27]. Our expression data from F1 and F2 hybrids show male sterility is associated with major alterations in genome-wide expression patterns. Clustering of trans eQTL is much less pronounced when mapping is restricted to fertile mice (Fig. S2), indicating trans hotspots in particular are associated with sterility. Each of the trans hotspots we identified overlaps a sterility QTL and/or interacts with at least one region containing a sterility QTL. One interpretation of this pattern is that divergent alleles with major effects on expression patterns are likely to cause hybrid incompatibilities. Trans regulators of gene expression must coordinate properly with cis regulators and other trans factors. The number and broad genomic distribution of regulated genes and co-factors provide many potential opportunities for incompatible interactions resulting in deleterious phenotypes in hybrids. Misexpression of a gene could result from a change in the set of positive or negative regulatory factors, or a mismatch in the spatiotemporal availability of these factors and the timing of expression. This hypothesis suggests genes in interacting regions with large cis eQTL and/or major alterations in spatiotemporal expression pattern between subspecies should be prioritized as candidates. Numerous studies of F1 hybrid sterility and evidence for reduced gene flow in hybrid zones have shown that the X chromosome plays a central role in hybrid male sterility in house mice [5], [62]–[65]. Our expression mapping results in F2s show that the X has a massive effect on testis gene expression, providing support for an important role of the X beyond the F1 generation. Most of the X chromosome is significantly enriched for QTL affecting expression of autosomal genes. The musculusPWD allele in the proximal X hotspot (10. 16 Mb–101. 19 Mb) has effects on expression suggestive of sterility (Table 2), consistent with the well-documented role of the musculus X in F1 sterility. This region harbors the largest-effect QTL identified for testis weight, sperm count, abnormal sperm head morphology, and number of offspring in X introgression experiments [25], [66]. Genes with functions related to fertility (sexual reproduction, fertilization, flagellum) were enriched among the QTT with low expression caused by the musculusPWD allele (Table 3). By contrast, the distal X hotspot shows little similarity to patterns observed in sterile F1 males. The distal hotspot overlaps several sterility QTL identified in Xmusculus introgression experiments (Supp. Table S2), but the domesticusWSB allele at hotspot eQTL is associated with the sterile expression pattern. These results reveal the presence of at least one novel locus on the X contributing to expression variation and potentially F2 sterility (Tables 2, S2). Fertility of DxM F1s, which carry the domesticusWSB X, and lack of enrichment of the distal hotspot QTT for transcripts misexpressed in F1s, indicate this locus interacts with one or more recessive musculusPWD autosomal loci. DNA-binding genes are enriched among QTT with higher expression, raising the possibility that the distal locus controls expression of regulatory genes, and the role in sterility is indirect. Variation within the trans hotspots on the X suggests each may harbor more than one sterility gene. The number of eQTL mapped, and the proportions of QTT with sterility-related characteristics, varied within the proximal and distal hotspots (Table S2). Furthermore, comparison of conditional mapping results using different markers on the X as covariates reveals differences in interaction patterns (Fig. S5). We identified a region on chromosome 17 with major effects on gene expression. Several lines of evidence implicate the known sterility gene Prdm9 as the underlying causative gene. First, the QTL for overexpression of X-linked transcripts (18. 46 Mb) and the peak in number of trans eQTL within the chromosome-17 hotspot (14. 69 Mb) are near Prdm9 (15. 68 Mb; Fig. 3B). Second, eQTL in the chromosome-17 hotspot largely show under- or overdominant effects, in contrast to trans eQTL elsewhere in the genome, which are mostly additive or dominant (Fig. 3A). This pattern is consistent with results from F1 crosses showing the most severe sterility phenotypes occur in males heterozygous at Prdm9 [67]. Finally, we find evidence for interactions between the chromosome-17 region and a musculusPWD allele on the proximal X chromosome, consistent with F1 studies [4]. If Prdm9 is the causative gene, our eQTL results provide novel insights into its role in hybrid sterility and gene regulation. In addition to the known interaction with the X chromosome, we find evidence for interaction with each autosomal locus used as a mapping covariate (Figs. S4; 5). The large number of interacting loci suggests that the DNA-binding function of Prdm9, which regulates recombination hotspots globally [73], [74], might be directly related to its role in sterility. Each Prdm9-binding site represents a potential incompatibility partner. Alternatively, disrupted regulation caused by Prdm9 might have cascading effects resulting in altered expression genome-wide. Although Prdm9 is predicted to have broad regulatory effects, previous evidence for effects on expression levels was limited to a small set of genes directly regulated by Prdm9 [27]. The combination of eQTL in the chromosome-17 hotspot (without covariates; Table 2) and eQTL dependent on interactions with eight autosomes and the X chromosome (Table 5) identifies 5,467 unique transcripts directly or indirectly affected by the region encompassing Prdm9. Chromosome 17 harbors a second, more distal sterility locus, Hstws, from musculus [18]. Hstws is necessary, in addition to the sterile Prdm9domesticus allele and the musculus X, to observe complete meiotic arrest, the most severe F1 phenotype [67]. We identified interactions between both the Prdm9 region and a distal chromosome 17 region with chromosomes 2,5, 10, and X (Fig. S4), suggesting loci on those chromosomes may be involved in the Prdm9- Hstws incompatibility. Overlap of sterility QTL with trans hotspots and/or interaction hotspots can refine estimates of the QTL position in some cases. For example, the trans hotspot on chromosome 17 is smaller than the coincident QTL for sperm count and testis weight (Fig. 3B). Moreover, the peak in number of trans eQTL is at the position closest to Prdm9. Chromosomes 5 and 10 are cases where trans eQTL and interaction eQTL patterns appear particularly useful in narrowing lists of candidate genes (Fig. S2) Functional annotation of QTT identifies affected pathways and processes associated with some hotspots, and provide clues about the mechanisms underlying sterility. Chromatin-related genes were overrepresented among QTT with lower expression associated with the sterile domesticusWSB allele at the chromosome 11 hotspot (Table 3). Mouse knockout models for two additional genes with eQTL in this region have spermatogenesis defects that might be related to chromatin; males with null alleles at the transcription factor Crem (cAMP responsive element modulator) showed defective spermiogenesis with aberrant post-meiotic gene expression [75]. Lmna (lamin A) knockouts have severely impaired spermatogenesis associated with failed chromosomal synapsis [76]. These patterns suggest prioritizing genes in the chromosome 11 hotspot with related functions. For example, 42 genes are involved in transcriptional regulation (Table 4). One of these genes (Hils1) is involved in chromatin remodeling during spermatogenesis and has evolved rapidly within rodents [77]. Males with hypomorphic Rad51c alleles are infertile due to arrest of spermatogenesis in early meiotic prophase I related to failed double-strand break repair by recombination [78]. Interactions between novel loci and better-characterized regions point to some promising candidates. For example, the chromosome-10 hotspot interacts with the proximal X and the chromosome-17 region containing Prdm9, the two loci with the most dramatic effects on expression. A gene within the chromosome-10 hotspot, Dnmt3l (DNA methyltransferase 3-like), plays a key role in epigenetic programming during spermatogenesis. Males carrying null alleles at Dnmt3l show phenotypes similar to those documented in F1s associated with the X-17 interaction, including hypogonadism, asynapsis during meiosis, abnormal formation of the sex body, and deregulation of X-linked and autosomal genes [79]–[82]. Dnmt3l does not have methyltransferase activity but shows sequence similarity to Dnmt3a and Dnmt3b, with which it interacts to promote de novo DNA methylation [83]. Misexpression of Dnmt3a was reported previously in sterile F1 hybrids [10]. Prdm9 is a histone methyltransferase; while speculative, an interaction between Dnmt3l and Prdm9 is a promising lead. Dnmt3l is essential for several epigenetic processes occurring at different stages of spermatogenesis, including paternal imprinting, transcriptional regulation, chromatin morphogenesis through meiosis, and the histone-protamine transition during spermiogenesis. Interestingly, Dnmt3l interacts with heterochromatin [80], similar to the Drosophila sterility gene OdsH [84]. Conditional mapping revealed several genomic regions involved in the interaction network that did not have previous evidence for involvement in sterility or expression (Fig. 5), indicating this mapping approach can uncover incompatibility loci without detectable marginal effects. Some of these interaction loci are very small, containing few enough genes that targeted functional evaluation would be feasible. For example interaction hotspots mapped with covariates on chromosomes 2,5, and X overlap on distal chromosome 7 (Fig. S4). This region spans 3. 1 Mb, encompassing 14 characterized RefSeq genes. We focused here on genome-wide patterns. Detailed characterization of individual loci, and analysis of gene co-expression networks including all related QTT, will yield additional information useful in pinpointing the disrupted pathways causing sterility and prioritizing candidates. The rate of accumulation of Dobzhansky-Muller incompatibilities and the evolution of reproductive barriers between incipient species depend on the genetic architecture of isolating traits. Theoretical models of DMI evolution assume that incompatibilities act independently on barrier traits [85], [86]. The complex pattern of interactions we report here violates this assumption: some sterility loci are involved in multiple incompatibilities. This aspect of the network we characterized is most consistent with branched developmental pathways [45] and gene networks models [47]. Theory that incorporates this non-independence as well as other biological characteristics of incompatibilities should continue to be pursued [45]–[47], [87]. Network characteristics are also key determinants in accurate modeling of gene-flow dynamics in zones of hybridization. Non-independence of incompatibilities due to interactions of sterility loci with multiple partners is likely to result in stronger selection and slower introgression at those loci because sterility phenotypes are expressed on a variety of genomic backgrounds. Future cline theory should incorporate epistasis with multiple partner loci. Remarkable progress in understanding the genetic basis of speciation has emerged from identification of a growing list of hybrid incompatibility genes [6] over the past 20 years. However, identification and functional characterization of hybrid incompatibility genes is feasible in only a few model organisms, and tremendous effort, time and resources are needed to identify a single gene. If this gene-by-gene approach continues as the standard in speciation genetics, it will be a long time before the number of genes and interactions identified is sufficient to reveal generalities of the speciation process. Moreover, general features of incompatibility networks, including the number and dominance of loci, types of interactions, and possibly particular developmental/regulatory pathways, are more likely to be shared among taxa than are specific incompatibility genes. The house mouse features a rich set of sophisticated genetic tools and resources, which facilitates collection of reliable genome-scale data and ultimately will enable functional characterization of candidate incompatibility genes. Although identification and characterization of reproductive barrier genes is not feasible in most species, the integrated mapping approach we employed is applicable in a broad range of organisms. For species pairs that can be crossed in the laboratory, a similar F2 intercross can be performed and sterility or inviability phenotypes can be measured. Informative marker discovery is straight-forward and relatively low cost using RADseq [88], and RNAseq or custom microarrays can be used to collect expression data from species without commercially available platforms. Functional annotation and nomination of candidate processes/pathways is possible if a genome sequence of the focal species or even a relatively distantly related taxon is available [89]. Even in species with very limited available gene annotation, the number of incompatibility loci and the nature of interactions between them can be estimated. Consequently, we suggest that network-centered approaches are powerful and have promise to substantially advance understanding of speciation.
Mice were maintained at the University of Wisconsin School of Medicine and Public Health mouse facility according to animal care protocols approved by the University of Wisconsin Animal Care and Use Committee. Reciprocal crosses of wild-derived inbred strains of M. m. domesticus (WSB/EiJ; domesticusWSB) and M. m. musculus (PWD/PhJ; musculusPWD) were performed to generate F1 hybrids. A total of 305 F2 males were generated by mating F1 siblings (294 from domesticusWSB female×musculusPWD male crosses and 11 from musculusPWD female×domesticusWSB male crosses). Male F2s were euthanized at 70 (±5) days of age. Five sterility phenotypes were quantified: testis weight, sperm count, sperm head shape, proportion of abnormal sperm, and seminiferous tubule area (see White et al 2011 for detailed methods). The left testis was flash-frozen in liquid nitrogen upon dissection and stored at −80°. Testes from musculusPWD (n = 8), domesticusWSB (n = 8), musculusPWD×domesticusWSB F1s (n = 6), and domesticusWSB×musculusPWD F1s (n = 4), were dissected using the same procedure to provide controls for expression analyses. Frozen testis samples were transferred to RNAlater-ICE buffer (Invitrogen, Grand Island, NY, USA), shipped to the Max Planck Institute in Plön and stored at −80° until processing. To identify classes of genes enriched among QTT, we used the DAVID functional annotation tool [43], [44], which integrates gene annotation information from several resources. Functionally related genes are clustered based on biological process, cellular compartment, molecular function, sequence features, protein domains, and protein interactions. To account for multiple comparisons, we used a significance threshold based on the false discovery rate (Benjamini) calculated within DAVID. We identified candidate genes in trans hotspots and among QTT that have roles in male reproduction and/or regulation of gene expression using reviews of male fertility [17] and meiosis [105] and gene ontology (GO) terms related to male reproduction, meiosis, or the regulation of gene expression: 0001059; 0001060; 0001109; 0001121; 0003006; 0006351; 0006352; 0006353; 0006354; 0006355; 0006360; 0006366; 0006383; 0006390; 0006396; 0006412; 0007127; 0007135; 0007140; 0007285; 0009008; 0009299; 0009300; 0009302; 0009304; 0010216; 0010468; 0010608; 0010628; 0010629; 0022414; 0023019; 0030724; 0030726; 0032775; 0032776; 0036206; 0040020; 0040029; 0042793; 0043046; 0043484; 0044030; 0045132; 0045835; 0045836; 0045892; 0045893; 0048133; 0048136; 0048140; 0048515; 0048610; 0050684; 0051037; 0051257; 0051604; 0070192; 0070613; 0070920; 0080188; 0090306; 0097393; 1901148; 1901311; 2000232; 2000235; 2000241; 2000242; 2000243. Many genes identified as candidates in publications were not annotated with related GO terms, highlighting the limitations of gene ontology. Moreover, genes causing sterility might not have functions obviously related to reproduction.
New species are created when barriers to reproduction form between groups of organisms that formerly interbred freely. Reduced fertility or viability of hybrid offspring is a common form of reproductive isolation. Hybrid defects are caused by negative interactions between genes that have undergone evolutionary change within each subgroup. Identifying genetic interactions causing disease or trait variation is very difficult, consequently there are few known hybrid incompatibility genes and even fewer cases where both interacting genes are known. Here, we combined mapping of gene expression levels in testis with previous results mapping male sterility traits in hybrid house mice. This new approach to finding genetic causes of reproductive barriers enabled us to identify a large number of hybrid incompatibilities, involving genomic regions with known roles in hybrid sterility and previously unknown regions. Understanding the number and type of genetic interactions is important for developing accurate models used to reconstruct speciation events. The genetics of hybrid sterility in mice may also contribute to understanding basic processes involved in male reproduction and causes of human infertility.
animal genetics
genomics
speciation
gene expression
genetics
hybridization
biology
evolutionary biology
gene networks
evolutionary processes
evolutionary genetics
Genomic Networks of Hybrid Sterility