Genetic conflict outweighs heterogametic incompatibility in the mouse hybrid zone?

BACKGROUND: The Mus musculus musculus/M. m. domesticus contact zone in Europe is characterised by sharp frequency discontinuities for sex chromosome markers at the centre of wider clines in allozyme frequencies.
RESULTS: We identify a triangular area (approximately 330 km2) where the musculus Y chromosome introgresses across this front for up to 22 km into domesticus territory. Introgression of the Y chromosome is accompanied by a perturbation of the census sex ratio: the sex ratio is significantly female biased in musculus localities and domesticus localities lacking Y chromosome introgression. In contrast, where the musculus Y is detected in domesticus localities, the sex ratio is close to parity, and significantly different from both classes of female biased localities. The geographic position of an abrupt cline in an X chromosome marker, and autosomal clines centred on the same position, seem unaffected by the musculus Y introgression.
CONCLUSION: We conclude that sex ratio distortion is playing a role in the geographic separation of speciation genes in this section of the mouse hybrid zone. We suggest that clines for genes involved in sex-ratio distortion have escaped from the centre of the mouse hybrid zone, causing a decay in the barrier to gene flow between the two house mouse taxa.

Background

Increasingly, evidence suggests that the sex chromosomes harbour more genes causing disruption of fertility and/or viability in hybrids than autosomes and hence will be under stronger selection in hybrid zones [1-8]. While the causes of this disruption have generally been assumed to follow the Dobzhansky-Muller model of accumulation of incompatibilities [8], genes involved in genetic conflicts have also been implicated [9-11]. In either case, gene flow of sex-linked markers across a hybrid zone is expected to be impeded, resulting in abrupt clines [12], the best examples of which have been reported from the contact zone between two house mouse subspecies, Mus musculus musculus and M. m. domesticus (see [13] for review).

The musculus-domesticus hybrid zone crosses the Jutland peninsula in Denmark and runs from the Baltic coast in East Holstein (northern Germany) across central Europe and the Balkans to the Black Sea (Figure 1; see also [14] and references therein). The transition of X-chromosome markers across various portions of this zone have been shown to be steep in comparison with autosomal markers (Bulgaria: [15]; southern Germany: [16,17]; Denmark: [18]; central Europe: [14]). Similarly, virtually no introgression was found for the Y chromosome in Bulgaria [19], northern Germany [20], and Denmark [21,22]. These observations have built up a remarkably consistent picture of wider autosomal clines versus more narrow sex linked clines spanning more than 2500 kilometres of the contact front between the subspecies in Europe, remarkable because it seems the dynamics of secondary contact are broadly similar over a very large geographic range relative to the dispersal of house mice [14,23]. This geographic repeatability and the fact that a model organism is involved, gives the musculus-domesticus contact zone the potential to play a central role in understanding the general features of secondary contact and the nature of barriers to gene flow between taxa, and in particular the action and nature of 'speciation' genes associated with sex chromosomes. However, a survey focused on the distribution of autosomal and sex-limited diagnostic markers across the Czech and Slovak Republics revealed an unexpectedly gradual transition of the Y chromosome compared to an X locus, and even to five of the autosomal loci analyzed [24], contradicting the results from other parts of the hybrid zone and causing us to question the generality of our understanding of the outcome of house mouse contact.

Figure 1

The . (a) The course of the zone is shown as bold line and the location of the study area as shaded rectangle. (b) The detail of the study area with sampling sites indicated.

The history of initial contact between M. m. musculus and M. m. domesticus is unclear. It has been suggested that source populations first met in the southern region of the current hybrid zone, and only more recently in central and northern Europe [25-27], with progressive contact from south to north similar to a zipper being pulled up through Europe. Secondary contact along such a front would initially produce parallel clines similar in width (i.e., concordant) and position (i.e., coincident). However, concordance would degenerate as clines at neutral loci widen compared to loci experiencing a barrier to gene flow [12]. In contrast, coincidence between cline centres is expected to be maintained by selection against hybrids. A tension zone [28], where selection against hybrids arises due to genetic rather than environmental factors, is not tied to a particular geographic area, and differential gene flow from the source regions towards a central hybrid sink will naturally push the zone toward troughs in population density and barriers in the environment. In the same way, clines at selected loci, each of which can be thought of as a barrier to gene flow, will be pushed by differential pressure of population density to coincide and their association is expected to form sharp steps at coincident cline centres [12]. In general, hybrid zone studies identify a set of loci with coincident cline centres, and an occasional outlier locus (or a few loci) with a displaced cline [29,30]. In cases of contact across a long front, such as the musculus-domesticus zone in Europe, a description of a set of loci with coincident cline centres forms the basis for further analysis. The consensus centre is used to map the path of the front, which may curve across the landscape at different geographic scales (cf Chorthippus grasshoppers: [31]; Bombina toads: [32]). This line is the geographic zero from which all distances of sampling localities are measured in order to estimate cline widths. However, if clines are not parallel the width of a non-coincident cline will be systematically overestimated and this is a potential explanation of the apparently gradual change of Y chromosome frequencies found by Munclinger et al. [24]: that the Y chromosome front is in fact following a different path from the consensus. It is not immediately clear how to estimate the positions of cline centres along a front allowing for this possibility since estimates of cline position depend on the model of cline cross-section (shape) assumed.

As a first approach it seems judicious to focus on inferring the path of the displaced cline while making as few assumptions as possible about its shape. Here we analyze the orientation of Y chromosome change in the Czech-Bavarian section of the musculus-domesticus hybrid zone assuming only that the cross-section of its cline is monotonic and that at the scale of the field sampling, its front is approximately linear. We present for the first time strong evidence that its orientation is strikingly different from the consensus contact front. Then we show that, when it is analyzed with respect to its own orientation rather than that of the consensus front, the Y chromosome cline width estimate is significantly reduced, and of the same order of magnitude as in other portions of the musculus-domesticus hybrid zone in Europe. We also map the geographic region where the path of the Y chromosome departs from the consensus front, describing a triangular salient of musculus Y chromosome introgression pushing far into domesticus territory. Pooling all census data from the musculus side of the zone, we demonstrate a significant female bias in the census sex ratio. A similar and significant female bias is also found in those domesticus localities without introgressed Y chromosomes. In contrast, domesticus localities where the musculus Y has introgressed show a sex ratio close to parity, the proportion of males being significantly higher here than in either of the other locality classes. We note that the mitochondrial genome displays a geographic pattern of introgression similar to the Y chromosome, though weaker and more stochastic. We discuss the relationship between the Y chromosome invasion and the perturbation of the sex ratio and suggest the Y invasion is due to manipulation of the sex ratio in its favour. We point out that this manipulation need not have been "selfish": male favouring factors would increase by natural selection when countering the selfish action of female favouring factors. Finally, we consider the possibility that the spatial pattern of the mitochondrial marker, similar to that of the Y, may be a trace of a previous introgression from musculus to domesticus, associated with the spread of such a female favouring element. If this is the case, the Y chromosome introgression is not the first breakdown of the consensus hybrid zone associated with perturbation of the sex ratio, but rather the current incursion of an ongoing arms race between elements distorting the sex ratio.

Results

Allele frequencies and the orientation of change across the field area

Frequencies of musculus alleles at the Btk locus at all the sampling sites under study are depicted as pie diagrams in Figure 2a. The transition of this X-chromosome marker from the domesticus to the musculus side is rather abrupt and the position and orientation of its front is very similar to the position and orientation of a consensus over autosomal loci [14]. On the other hand, the spatial distribution of Y chromosome frequencies is very different from that of other markers (Figure 2b). Figure 3 shows log likelihood (support) profiles for monotonic change in allele frequencies from musculus-like to domesticus-like, moving through the field area following a compass bearing lying within a 90° range including due west (W) and northwest (NW). The Y chromosome is a clear outlier, its most likely change being oriented roughly 45° clockwise of the other loci. Summing log likelihood profiles over all loci produces a consensus profile (in bold black in Figure 3) under the assumption that the orientation of change in allele frequencies is the same for all the loci.

Figure 2

Allele frequencies across the field area. Frequencies of musculus (red) and domesticus (blue) alleles are shown for the X chromosome marker (a), Y chromosome (b), and mtDNA (c). The pies are proportional to sample sizes. The bold lines indicate orientation of the X and Y fronts, respectively, estimated using the PAVA algorithm; the X marker orientation is almost identical to that of the consensus over all loci (see Table 1). M and D denote musculus and domesticus regions delimited by the green line.

Figure 3

Likelihood (support) profiles for monotonic change in allele frequencies. The profiles show likelihood of monotonic change of individual markers from M. m. musculus to M. m. domesticus, moving through the field area following a compass bearing in a 90° range including due west (W) and northwest (NW). Yellow: Gpd1; pink: Abpa; red: Btk (X chromosome); brown: Mpi; orange: Idh1; purple: Es1; dark blue:Sod1; green: Np; blue: Y chromosome; bold black: consensus support over all loci. For clarity of presentation the profiles are smoothed over a window of 3°. Since the smoothing removes precise details around the maxima these may not correspond precisely to Table 1 which is based on the original (unsmoothed) profiles.

Table 1 shows the maximum likelihood estimate (MLE) and two unit support bounds for the orientation of each locus and the consensus. The latter falls almost exactly on the MLE orientation of the X marker (Btk) and within the two unit support bounds for the orientation of every other locus except Np and the Y chromosome, while the support intervals for all loci except the Y overlap with that of the consensus. Hence there is strong evidence that the Y differs from the consensus and weak evidence that Np may follow a different orientation.

Table 1

Maximum-likelihood estimates (MLE) and 2LL-unit support bounds for orientation of change (in degrees from due north) of allele frequencies for individual loci; two consensus orientations were estimated, one summed over all nine loci (Consensus) and the second one excluding Np and Y chromosome (see text for details).

LL profile	Orientation estimate and 2LL-unit support			Comparison to consensus*
	MLE	Lower	Upper	ΔLL	P
Gpd1	-79.7	-94.0	-65.6	1.39	0.09
Abpa	-74.3	-90.4	-57.8	0.71	0.23
Idh1	-69.7	-88.9	-53.4	0.72	0.23
Btk (X chr.)	-67.3	-88.9	-60.4	0.03	0.81
Sod1	-64.2	-69.4	-49.9	0.67	0.25
Es1	-63.3	-69.7	-55.2	1.12	0.13
Mpi	-61.8	-88.7	-57.7	0.17	0.56
Np	-54.5	-64.2	-46.0	2.66	0.02*
Y chr.	-28.4	-32.7	-24.2	3.62	0.007**
Consensus	-67.4	-68.3	-63.8
Consensus*	-67.4	-68.3	-65.6

Treating both Np and the Y chromosome as outliers, we can estimate a reduced consensus (hereafter referred to as consensus*) likelihood profile from the remaining loci, assessing potential effects of their inclusion or exclusion. Measuring change at a locus following the consensus* MLE orientation, rather than its own MLE, causes a drop (ΔLL) in support for monotonic change and the significance (P) of these drops can be weighed by comparison of 2ΔLL to the χ2 distribution with one degree of freedom. When significance is judged at the 99% level all non-Y loci can be treated as changing in the consensus* orientation without significant reduction in their likelihood of monotonic change. Furthermore, regardless of inclusion or exclusion of both outliers, the consensus MLE orientation remains the same (cf. consensus and consensus* MLE's in Table 1), i.e., -67.4° from due north.

The coincidence of cline centres and the consensus front

Figure 4a shows the most likely monotonic change in allele frequencies at all eight non-Y loci when clines are traversed with respect to the consensus MLE orientation. While there is a range of cline shapes, this variation revolves round an interval of less than one kilometre, within which all loci pass through 50%. Thus the clines at these loci are highly coincident, and variation in cline width and symmetry forms a bowtie around a knot of coincident centres. For the purposes of illustration, the Y chromosome cline is also shown, as estimated assuming its displacement was parallel to the consensus. While this gives some notion of the scale of Y chromosome introgression (> 20 km), the width and shape of the cline will be misrepresented, because it is not being looked at in perpendicular cross-section. Measuring cline widths using the software Analyse (see Methods), the width of the Y cline along this consensus orientation is estimated at 28.1 km (lower and upper bounds: 19.5–43.5) for the sigmoid model and 19.0 (14.5–25.6) km for the asymmetric stepped model, with little justification for the complicated step hypothesis compared to the simpler sigmoid change.

Figure 4

Monotonic clines for different orientations. (a) Clines assuming the markers change in the consensus orientation -67.4° from due north. Yellow: Gpd1; pink: Abpa; red: Btk (X chromosome); brown: Mpi; orange: Idh1; purple: Es1; dark blue:Sod1; green: Np; blue: Y chromosome. (b) Monotonic clines for the X marker (red) and uncorrected Y chromosome (dotted blue) assuming the consensus orientation; solid blue line: the corrected Y chromosome cline assuming the marker changes monotonically along its own most likely orientation of change (-28.4° from due north; cf. Table 1). Note that cline centres have been arbitrarily shifted on the x-axis to allow a direct comparison of cline shapes.

The width of X and Y clines

Figure 4b highlights the X and Y clines from Figure 4a, shifting the Y cline (blue, dotted) arbitrarily to the right so that its misrepresented shape can be compared to the X chromosome (red). In addition, the corrected Y cline shape is shown (blue), measured with respect to its own most likely orientation of change (-28.4° from due north). The corrected width is 14.6 (10.5–21.5) km for the sigmoid model whereas the asymmetric stepped cline width falls within the support interval 0.1–6.4 km, with the stepped hypothesis justified in comparison to the sigmoid one (ΔLL = 8.58; df = 4 ; P = 0.002). The corrected stepped Y cline fit is significantly better than the uncorrected one (ΔLL = 11.20; df = 1; P = 2 × 10-6). The corrected Y estimates are of the same magnitude as for the Btk locus on the X chromosome: a sigmoid Btk is estimated at 9.7 (6.9–16.1), a stepped Btk at 4.7 (0.1–8.2) km, and again, the stepped hypothesis is justified in comparison to the sigmoid model (ΔLL = 7.73; df = 4; P = 0.006).

The geography of the zone, Y introgression, and sex ratios

Figure 2a shows the sampling of allele frequencies at the Btk locus, and maps the consensus front over this locus and the autosomal loci onto the field area. This front is used to divide the field area into two regions. To the east of the consensus front is the M. m. musculus side of the zone (labelled M), while to the west is the M. m. domesticus range (labelled D). Figure 2b shows the sampling of allele frequencies at the Y locus and compares the centre of change of the Y cline to that of the remaining loci. Measuring the triangular area between the two fronts allows us to quantify the geographic extent of musculus Y introgression. A line departing at right angles from the consensus front (green line) to reach the Y front (yellow line) and passing through the southwestern-most of the localities known to be introgressed, covers a distance of ~22 km, and encloses an area of ~330 km2.

Table 2 divides sampling localities into three categories: those in the M region, those in the D region where no Y chromosome introgression has been detected (DYD), and those in the D region where the musculus Y has been detected (DYM). Only localities where at least one Y chromosome has been assayed are included in the analysis.

Table 2

Frequency of the musculus Y chromosome, hybrid indices (taken as frequencies of musculus alleles averaged over eight non-Y loci) in males and females, and the census sex ratio according to category of sampling locality.

	DYD	DYM	M
N (Localities)	20	40	55
Y total	187	448	417
Y musculus	0	369	412
Freq. Y	0	0.824	0.99
Census total	509	1062	1092
Males total	226	545	491
Alleles musculus (Total)	68(2479)	1343(6682)	5124(5551)
Freq. musculus	0.03	0.20	0.92
Females total	283	517	601
Alleles musculus (Total)	54(2858)	1366(6840)	7059(7546)
Freq. musculus	0.02	0.20	0.94
Proportion of males	0.44	0.51	0.45
Pr(Bin)	0.0065	0.2037	0.0005

Pr(Bin) stands for the probability that under the binomial distribution we observe a similar or more extreme deviation from 1:1 by chance alone.

There is a clear and significant female bias in the census sex ratio in both the M and DYD localities. In contrast, the DYM localities show a sex ratio close to parity. Pairwise differences between the categories can be tested by contrasting two hypotheses: H0, two censuses are drawn from populations with the same underlying sex ratio p0; H1, two censuses are drawn from populations with different underlying sex ratios, p1,A, p1,B. In each case the underlying sex ratio is estimated as the proportion of males in the relevant set of observations. This proportion is the maximum likelihood estimate of p. The probability of a set of observations, conditioned on an underlying sex ratio, is calculated according to the binomial distribution with parameter p. The likelihood of the observations is proportional to this probability. As the hypotheses concern the same set of observations, log likelihoods for the two hypotheses can be compared using a G test with one degree of freedom. The comparisons are as follows: DYD| DYM: p0 = 771/1571; p1,A, p1,B = 226/509, 545/1062; ΔLL = 3.300; P = 0.010. DYM| M: p0 = 518/1077; p1,A, p1,B = 491/1092, 545/1062; ΔLL = 4.358; P = 0.003. DYD| M: p0 = 717/1601; p1,A, p1,B = 491/1092, 226/509; ΔLL = 0.022;P = 0.833. Assuming DYD and M have the same underlying sex ratio (MLEs 0.44, 0.45), and comparing them to DYM, we have DYM|{DYD, M}p0 = 1262/2663; p1,A, p1,B = 545/1062, 717/1601; ΔLL = 5.466; P = 0.001.

Discussion and conclusion

In this paper, we demonstrate two striking observations with regard to the Czech-Bavarian mouse contact zone: (1) a large geographic anomaly with respect to the musculus Y chromosome; and (2) a significant difference in the census sex ratio in those domesticus localities where the musculus Y chromosome is present, compared to both musculus localities and domesticus localities without Y introgression. We suggest it is unlikely that these two patterns are coincidental. We discuss each observation in turn, and then explore whether a single cause might parsimoniously explain both.

Geographic anomaly of Y introgression

We present strong evidence that the M. m. musculus Y chromosome is introgressed on the M. m. domesticus side of the Czech-Bavarian contact zone between the two house mouse subspecies, penetrating some 22 km and over an extensive region of approximately 330 km2. The introgression is shown particularly clearly when monotonic clines are fitted to the data at multiple loci (Figure 4a): all clines save the Y chromosome change abruptly across an interval of less than one kilometre, i.e., less than the estimated scale of movement of mice per generation [14,23]. It may be that monotonic cline hypotheses, which allow discrete jumps in allele frequencies between localities, are particularly suited to describing the demic nature of house mouse population density. Compared to the strong coincidence of cline centres over eight loci, the non-coincidence of the Y chromosome is striking.

Another notable characteristic of the Y chromosome transition is its width (when assessed for the correct orientation), which is comparable to other portions of the hybrid zone. This might seem surprising because the Y has escaped linkage disequilibrium with other loci at the centre of the consensus zone. Such association is expected to steepen clines [12] thus we might expect the escapee Y cline to be wider than in portions of the hybrid zone where it coincides with the consensus. However, while a narrow cline will tend to steepen wider associated clines, we do not expect the reverse to the same degree – a wide cline is likely to have little effect in steeping an associated cline that already changes at a much smaller scale. In other words, the consistent width of the Y indicates it is likely to be a driver of the steep change at other loci when in coincidence at secondary contact rather than a follower (through statistical association) of action at some other, as yet undetected, locus. This is perhaps unsurprising given evidence that the differentiation of the Y-chromosome lineages may contribute to partial reproductive isolation between Mus subspecies [33].

Two contrasting scenarios could equally explain the current geographic pattern of the Y chromosome introgression: (A) a static consensus centre and westward geographic incursion of the musculus Y chromosome, forming a bulge or salient; and (B) a dynamic consensus centre that has moved eastwards, leaving an invagination of the musculus Y chromosome in its wake. These salient/invagination scenarios are not mutually exclusive: both the consensus centre and the musculus Y front may have moved. The major evolutionary implication is, however, unaffected by these details: in the current introgression region, the musculus Y chromosome has obviously been successful relative to the domesticus Y chromosome. This is clear if we assume the musculus Y has pushed west through the consensus centre, replacing the domesticus Y on its own background, and forming a salient. However, it is also true if the consensus centre has moved east: because the coincidence of clines is expected to be maintained by a tension zone, we expect the domesticus Y to advance east along with the consensus centre. If an invagination of musculus Y chromosomes is formed during this advance, it is because they have resisted the movement of the consensus centre, again replacing domesticus Y chromosomes on their own background. Thus the evolutionary implications are the same irrespective of how the displacement between the Y front and the consensus centre came about. It is their relative positions which inform us that the musculus Y has been successful in comparison with the domesticus Y.

The musculus Y chromosome introgression is surprising for two reasons, one specific to the house mouse zone, and one more general. First, domesticus males are known to be more aggressive than musculus males [34-36], and given the male-dominated harem mating system of house mice, we might expect the domesticus Y to have an advantage over the musculus Y where the two subspecies meet. Second, and more generally, Haldane's rule [37] predicts that it is the heterogametic sex that is preferentially affected by hybrid incompatibilities. When the heterogametic sex is female, mitochondrial transmission between diverged populations should be impeded: as pointed out by Mallet [38] this may explain why mitochondrial introgression is more common in mammals and flies than in birds and bees [39]. Britton-Davidian et al. [40] demonstrated that Haldane's rule holds for hybrids between musculus and domesticus mice from Denmark: male F1's are sterile and so we might expect Y chromosome transmission between musculus and domesticus to be impeded because Y chromosomes can only be transmitted through the heterogametic sex, i.e., males. The sharp and coincident Y chromosome clines seen elsewhere in the European house mouse hybrid zone are also consistent with this argument.

The above reasons add to the degree of surprise on observing the musculus Y introgression, but they may be balanced against other lines of evidence. First, it has been shown that both musculus males and musculus females prefer individuals of the same subspecies while domesticus mice are not choosy [36,41-47], though homosubspecific preference has also been detected in domesticus individuals in some populations in Denmark and central Europe (G. Ganem and B. Bímová, unpubl. results). All other things being equal, reduced choosiness in domesticus should favour introgression of alleles from musculus into domesticus, relative to the reverse. This is not an argument for selection favouring introgression, rather for a barrier to neutral introgression (choosiness) being removed in one direction. Reduced choosiness would only confer an invasive advantage to musculus alleles if finding a mate were difficult. In our study area, this situation seems unlikely since there was little difficulty trapping mice of both sexes. Furthermore, a neutral asymmetry in prezygotic barriers seems unlikely to give rise to the observed salient of Y introgression since its extent is large relative to the movement of mice and its front is sharp: to be large and neutral, the introgression process would have to be old, and an old stochastic process would not exhibit a sharp front. Finally, it should be noted that mice toward the centre of the hybrid zone are the descendants of many generations of intercrossing making it unlikely for their mate choice behaviour to be either pure musculus-like or pure domesticus-like.

Second, although Britton-Davidian et al. [40] demonstrate Haldane's rule for musculus-domesticus hybrids, male F1's being infertile, the same study reveals that the homogametic sex may also be affected: females appear significantly underrepresented in the sex ratio at birth in crosses between Danish musculus and domesticus (in total, 231 males : 180 females; sex ratio p = 0.562; Pr(Bin) = 0.007). The significance of this bias in the birth sex ratio was not noted at the time of the study because results were assessed relative to the χ2 distribution, which only approximates to comparison with the binomial distribution if numbers are large. If only males were adversely affected in the F1 generation, Y introgression might be particularly surprising. This crossing data suggests, however, that both sexes are adversely affected, although in different ways. Third, crosses between inbred Czech musculus and inbred French domesticus, described in the same study [40], revealed that F1's fathered by Czech males were not infertile. These crosses may have involved genetic backgrounds more similar to the current Czech-Bavarian contact than the Danish crosses, and raise the possibility of a unidirectional breakdown of Haldane's rule facilitating musculus Y introgression.

Sex ratio anomaly in localities with Y introgression

We demonstrate that when localities with Y introgression are excluded the census sex ratio in both the domesticus and musculus regions of the field area is significantly, and similarly, female biased. In contrast, the sex ratio in localities with Y introgression is close to parity and significantly different from both musculus and non-introgressed domesticus (Table 2). Two potential sources of uncertainty exist in the sex ratio measurements. First, genetic studies of sex determination in laboratory mice have shown that the sex determining pathway is sensitive to genetic background [48-51]. As the musculus Y is introgressed onto a new background, we must address the possibility that the associated sex ratio perturbation is phenotypic, but not genetic, and caused by disruption of the sex determination pathway. However, we can discount this possibility because 1) the majority of individuals were typed both phenotypically and genetically, with a mismatch rate of less than 1%, consistent simply with a very low rate of human error; 2) if the musculus Y on the domesticus background did indeed cause developmental sex reversal, XY individuals being phenotypically females [48-51], then we should expect a relative excess of phenotypic females in introgressed localities, and not the observed relative female deficit. Second, the propensity of females versus males to enter traps is unknown and so we do not know if the census sex ratio is an unbiased reflection of the tertiary sex ratio. However, in a similar fashion to the comparison of relative cline positions above, the relative census sex ratios are revealing. If the trapping propensity of sexes is similar throughout the field area, it is clear that the tertiary sex ratio in introgressed localities has a higher proportion of males than other localities. It is possible that Y introgression changes the behaviour of males, making them more likely to enter traps. However, this would represent elevated risk-taking behaviour. It should be remembered that the loss of mice from the population due to once-yearly trapping by scientists is likely negligible compared to other losses. Then, elevated risk-taking behaviour should result in a deficit of males in our samples, because of the (fatal) risks they have taken throughout the year, rather than the observed relative male excess. The elevated proportion of males trapped in introgressed localities is therefore unlikely to be due to a change in male behaviour and we conclude that it is the tertiary sex ratio itself that in introgressed localities is composed of a higher proportion of males than other localities.

The estimated sex ratio varies among localities and it is only the strength of the effect in introgressed localities combined with massive sampling effort across more than 120 localities in total that allows a significant difference between locality categories to be revealed. Variation in census sex ratio between localities might be expected not only because of the stochastic nature of finite samples: for example, Rosenfeld et al. [52] show that female mice bias the secondary sex ratio of their offspring as a function of their diet. Food availability of commensal mice will vary naturally from farm to farm across the field area, and from year to year. Similarly, parasite prevalence and load is typically highly spatiotemporally heterogeneous, and Kaòková et al. [53] demonstrate perturbation of the secondary sex ratio during the latent phase of toxoplasmosis in mice. For such reasons it is probably unhelpful to attempt to reason in terms of a singular evolutionarily stable sex ratio in wild mouse populations. The observed sex ratio will arise as a function of local resource availability and the interaction of numerous competing strategies for investing those resources and maximising evolutionary gains. The commonality of these strategies is that each is played out to some extent through favouring numbers of individuals of one sex over the other. Where parasite genes benefit from such a strategy, the interaction is usually called host manipulation [53], where particular mouse genes benefit, genetic conflict [54,55], while in the absence of any such manipulation/conflict, strategy at the level of parental individuals is shaped only by natural selection [56]. A recent study (J. Piálek, unpubl. results) shows that significant parent-specific variation in the outcome of such strategies can be demonstrated by measuring secondary sex ratio over multiple litters. Of 258 wild caught mice pairs taken from 29 populations across the Czech-Bavarian portion of the zone and crossed in the laboratory under uniform intermediate nutrition levels, one pair from the musculus region has shown significant male bias (104 males, 76 females, sex ratio p = 0.58, Pr(Bin) = 0.022). It would be interesting to see what further variation might be revealed under high and low nutrition regimes and parasite loads.

Reconciling the Y chromosome geographic and sex ratio patterns

The geographic pattern indicates the musculus Y chromosome has been successful relative to the domesticus Y on the domesticus background. The sex ratio pattern indicates that domesticus localities with musculus Y introgression have an elevated proportion of males in the tertiary sex ratio. A parsimonious explanation of both observations is that the success of the introgressing Y chromosome is because the combination of the Y and its novel genetic background increases the proportion of fertile males. As a result, where both domesticus and musculus Y chromosomes are present, the musculus Y tends to replace the domesticus alternative. An interesting question then arises: is this tendency acting in the same direction as natural selection or against the action of natural selection? The answer depends on whether, in the absence of Y introgression, the populations would be at equilibrium sex ratio regarding investment in the sexes [57], or biased in favour of females by selfish elements. In the first case, the action of the introgressing Y is selfish, i.e., in conflict with natural selection. In the second case the interests of the introgressing Y chromosome are aligned with natural selection, and its success may be speeded as it is redressing an imbalance caused by other (selfish) elements. In either case, selfish elements, and therefore genetic conflicts, are implicated. Deciding which of the alternative contexts is more likely requires information about potential differential investment in the sexes, of which little has been published. Male and female birth weights are very similar in these mice (J. Piálek, unpubl. results) and while it seems males may demand more resources during lactation, no female sex ratio bias has been reported in other areas of the mouse contact front, despite large sample sizes. These observations would suggest that differential investment in the sexes is in general not sufficient for the mouse sex ratio to differ significantly from parity. Then, the marked female trapping bias in the Czech-Bavarian zone observed outside of introgressed localities may indicate the presence of female favouring distortion elements, and the introgressing Y chromosome would be redressing a female imbalance in the sex ratio and thus be acting in the same direction as natural selection. In this case we might seek evidence of the previous spread of a selfish female biasing factor.

The similarity between the geographic pattern of the Y chromosome introgression and the musculus mitochondrion might seem unremarkable – once a barrier to gene flow is broken down, perhaps many elements might flow through. On reflection, however, the mitochondrial pattern is perplexing, because the mitochondrial genome and Y chromosome are never (or at most, extremely rarely) co-inherited, and so musculus mitochondrial variants quite simply could not have hitch-hiked through statistical association during the selective success of the introgressing Y. Furthermore, the interests of maternally and paternally inherited elements are in conflict with respect to the sex ratio, making it very unlikely that the mitochondrial and Y chromosome geographic patterns arose simultaneously. We are forced to conclude that one followed the other without statistical association, yet ending up covering a similar area of the field. Their conflicting interests over the sex ratio could explain this conundrum – the mitochondrion would hitch-hike with a female biasing distorter of the sex ratio, and an introgressing Y which rescues that distortion would then be favoured by natural selection, wherever the distorter had passed.

Irrespective of speculation regarding the direction of interests, selfish or otherwise, of the musculus Y chromosome, if the increased proportion of males in introgressed localities is the cause of the musculus Y success, then this is the first case where genetic conflict has been implicated as a force causing the decay of a species barrier. Taking other transects of the mouse contact hybrid zone in Europe as representative of the Czech-Bavarian contact before the Y incursion, originally there were sharp coincident changes for sex chromosome markers, and centred on these more gradual clines at autosomal loci: the actions of many barrier loci were focused in a single geographic region. Now, the Y chromosome cline has escaped the forces holding these barriers in coincidence and the total effective barrier in that region is reduced: no change in Y chromosome state occurs across the consensus hybrid zone centre in the southern half of our field area.

In the Introduction we noted that there is growing evidence that the sex chromosomes harbour more genes causing disruption of fertility and/or viability in hybrids than autosomes, and that while the causes of this disruption have generally been assumed to follow the Dobzhansky-Muller model of accumulation of incompatibilities [4], there is also increasing evidence that genes involved in genetic conflicts may be implicated [9-11,58,59]. We stated that irrespective of the causes of such disruption, gene flow of sex linked markers across a hybrid zone is expected to be impeded resulting in abrupt clines [12]. In the current study we show an abrupt Y chromosome cline has escaped the consensus centre of the mouse hybrid zone, and that Y introgression is associated with a significant perturbation in the tertiary sex ratio towards males. Either these two striking patterns are coincidental, or the results of genetic conflict are currently laid out across the landscape of the hybrid zone, having overcome the selective forces of heterogametic incompatibility that focus genetic change into a single geographic region.

Methods

Sampling

In total, 2311 house mice were trapped at 126 sites scattered across a belt 145 km long and 50 km wide, stretched from north-eastern Bavaria (Germany) to western Bohemia (Czech Republic). The sampling sites are listed in the Appendix and their position is shown in Fig 1. All trapped mice were euthanized and dissected either in a field laboratory or in the Studenec laboratory and their sex and reproductive status was inspected immediately after dissection.

Since sampled alleles (see below) may be non-independent due to relatedness and deviations from Hardy-Weinberg equilibrium, we used the effective number of alleles calculated according to [60] and [23]:

where N is the number of diploid individuals sampled at a site, FIS is the measure of deficit/excess of heterozygotes, and FST represents the measure of relatedness in each population (note that the term (1 - FST) disappeared from the original formula as a typing error in [60], a mutation propagated in [23]). This means that if all sampled individuals were completely related, the effective sample size would be 1 and, conversely, if individuals were completely unrelated, this would result in Ne ranging from N (when FIS = 1) to 2N (when FIS = 0). Thus for high FST values the influence of large samples is substantially reduced. For haploid markers FIS = 0 and so the formula (1) simplifies to:

In the case of the X chromosome, estimation of Ne is slightly more complex because the chromosome is diploid in females and haploid in males. If Nm and Nf are the number of sampled males and females, respectively, and is the measure of deficit/excess of heterozygotes in females, then the effective number of alleles is given as:

FIS was estimated for each marker and for each population separately, whereas for FST the maximum likelihood estimate (MLE) was calculated for each population as a combined estimate from genotypes at all loci (see below) using both the Beta and normal distribution options from Analyse 1.3 [61]. The consensus* orientation estimates were similar irrespective of the model used for FST estimation. Here we show only results for the Beta distribution option since this results in higher estimates of FST and, hence, is more conservative with regard to effective sample sizes.

An estimation of FIS is problematic in some extreme cases. For example, if one of the classes (pp, pq, qq) contains a zero element and the second one is close to zero, then MLE of FIS will be 1 (e.g., when pp = 50, pq = 0, qq = 1) or -1 (e.g., when pp = 0, pq = 1, qq = 50) irrespective of the "true" value. So we used the following rule of thumb: if one of the counts was zero and the second lowest count was less than 3, then we assumed we had not enough information to justify an FIS deviation from zero. Even though this choice is obviously arbitrary the resulting estimates are rather conservative and the number of outliers is reduced.

Markers

DNA was isolated from tissues using proteinase K digestion and subsequent extraction either with phenol-chloroform and ethanol precipitation [61] or using the DNeasy® 96 Tissue Kit (QIAGEN), following standard protocols or the manufacturer's instructions.

Males (N = 1050; 115 sites) were typed for the presence or absence of an 18-bp deletion, located within the last exon of the Zinc finger protein 2, Y linked, gene (Zfy2) using the method given in [62]. The deletion is fixed in M. m. musculus and absent in M. m. domesticus [23,63]. In addition, the results were verified using a Y-specific microsatellite (TTTTG repeat) at the 5'-end of the second intron of Zfy2 under the same conditions as described in [23]. Since there was no incongruence between the two Y-chromosome markers, they are henceforth referred to as a single marker.

To set the Y chromosome results within the context of the hybrid zone, we used one X-linked and seven autosomal markers. For the X-chromosome marker, we scored a B1 SINE retroelement, mapping the Btk gene (N = 2171; 126 sites): this insertion is present in M. m. domesticus and absent in M. m. musculus (see [23,64] for details). For autosomal markers, six allozyme loci (Es1, Gpd1, Idh1, Mpi, Np, Sod1; N = 1965–2137; 126 sites) were used as described in [14]. The seventh autosomal marker maps the Abpa gene coding for the alpha-subunit of the salivary Androgen binding protein (Chromosome 7; N = 2056; 126 sites) following [22].

Finally, we compared the data with the distribution of a mitochondrial DNA marker across the area under study. This marker was represented by the BamHI restriction site mapped at position 3565 of the standard mouse mtDNA sequence [65] which has been shown to be a reliable marker distinguishing the two house mouse subspecies [66,67]. 1939 specimens from 121 sites were analyzed as described in [67].

Analyses

The orientation of change in allele frequencies in two-dimensional space was estimated with the Pooled Adjacent Violators Algorithm (PAVA), originally described by Brunk [68] as a method for finding the maximum likelihood monotonic cline over a set of observations (see also [69,70]). One advantage of this method is that it is not dependent on any particular cline model: PAVA assumes no particular shape of change other than monotonic increase or decrease across a series of sampling sites. For this reason, the distance between localities does not affect the outcome of PAVA, only their ordering. An ordering is found by taking the orthogonal projection of localities onto a cross-section through the field area, and then noting the order in which the projected localities fall along this cross-section. An implicit assumption is that at the scale of observation the path of the contact front is linear. For a finite number N of sampling localities there are ≤ N(N – 1) distinct orderings so a complete likelihood profile for monotonic change across the continuous interval of orientations [0°,360°] can be calculated. No sampling or search algorithms need be applied. Likelihood profiles and their underlying monotonic cline estimates were calculated for each marker using a routine written in Mathematica [71]. In order to test whether orientations were similar across loci, we compared likelihood profiles as described in [14] and [58]. The width of a cline is normally measured as the inverse of the maximum gradient. It is not clear how best to apply this metric to PAVA monotonic clines, so widths were estimated by fitting sigmoid or stepped models using the Analyse program [60].

Authors' contributions

MM, PM, PD, BB, JP collected the material and MM, PM, and JP participated in the design of the study. MM analyzed allozyme markers and participated in the analysis of the microsatellite locus. PM, PD, BB performed laboratory experiments except those carried out by MM. JP provided unpublished data on sex ratio from breeding experiments and logistical support for the study. SJEB analyzed the data, interpreted the results and wrote the paper as it arose out of intensive discussion with MM. All authors read and approved the final version of the manuscript.

Additional file 1

A list of sampling sites. For each site, geographic coordinates (Latitude, Longitude), category (M = musculus population; DyM = domesticus population with musculus Y introgressed; DyD = domesticus population with no introgressed musculus Y), total number of individuals and the number of males trapped, and the total number of alleles and the number of domesticus alleles for each marker (mtDNA, Y chromosome, Btk, Abpa, Es1, Gpd1, Idh1, Mpi, Np, Sod1) are given.

Click here for file