|Home | About | Journals | Submit | Contact Us | Français|
Several lines of evidence suggest that the X chromosome plays a large role in intrinsic postzygotic isolation. The role of the Z chromosome in speciation is much less understood. To explore the role of the Z chromosome in reproductive isolation, we studied nucleotide variation in two closely related bird species, the Thrush Nightingale (Luscinia luscinia) and the Common Nightingale (L. megarhynchos). These species are isolated by incomplete prezygotic isolation and female hybrid sterility. We sequenced introns of four Z-linked and eight autosomal loci and analyzed patterns of polymorphism and divergence using a divergence-with-gene flow framework. Our results suggest that the nightingale species diverged approximately 1.8 Mya. We found strong evidence of gene flow after divergence in both directions, although more introgression occurred from L. megarhynchos into L. luscinia. Gene flow was significantly higher on the autosomes than on the Z chromosome. Our results support the idea that the Z chromosome plays an important role in intrinsic postzygotic isolation in birds, although it may also contribute to the evolution of prezygotic isolation through sexual selection. This highlights the similarities in the genetic basis of reproductive isolation between organisms with heterogametic males and organisms with heterogametic females during the early stages of speciation.
Understanding the genetic basis of reproductive isolation may provide important insights into the mechanisms responsible for the origin of new species. Previous studies have suggested that the X chromosome contributes disproportionately to the origin of intrinsic postzygotic isolation, including both hybrid sterility and hybrid inviability. This has been demonstrated in laboratory crosses (e.g., True et al. 1996; Turelli and Begun 1997; Tao et al. 2003; Oka et al. 2004; Storchová et al. 2004; Masly and Presgraves 2007; Good et al. 2008) as well as in studies of hybrid zones (Payseur et al. 2004; Macholán et al. 2007; Teeter et al. 2008). The disproportionately large role of the X chromosome in intrinsic postzygotic isolation may have several distinct causes including exposure of recessive mutations in hemizygous individuals, or a higher density of genes contributing to postzygotic isolation on the X (Masly and Presgraves 2007). A higher density of X-linked genes contributing to postzygotic isolation may, in turn, be caused by one of several mechanisms including faster X evolution, the accumulation of sex-chromosome meiotic drive mechanisms, or disruption of epigenetic regulation of the X chromosome in the male germline (reviewed in Coyne and Orr 2004; Presgraves 2008).
Despite the extensive literature on the role of the X chromosome in reproductive isolation, relatively few studies have focused on the role of the Z chromosome in speciation (e.g., Jiggins et al. 2001; Sætre et al. 2003, 2007; Carling and Brumfield 2008, 2009). However, several lines of evidence suggest that the Z chromosome might be even more important in reproductive isolation than the X chromosome. The first is that faster evolution of the Z compared to the autosomes is clear in birds (Mank et al. 2007; Ellegren 2009), whereas evidence for “faster X evolution” is equivocal (e.g., Betancourt et al. 2002; Thornton and Long 2002; Torgerson and Singh 2003; Counterman et al. 2004; Connallon 2007; Baines et al. 2008). A higher rate of adaptive evolution on the Z chromosome relative to autosomes could underlie either intrinsic or extrinsic reproductive isolation. Second, the Z chromosome may promote speciation through sexual selection (Ritchie 2007). Several studies suggest that sexual selection is facilitated in species with female heterogamety, because sons inherit the Z chromosome directly from their father (Hastings 1994; Reeve and Pfennig 2003; Kirkpatrick and Hall 2004; Albert and Otto 2005). An important prediction of these studies is that male traits and female preferences will often be Z-linked, and indeed, there is some evidence for such Z-linkage in birds (Sætre et al. 2003; Sæther et al. 2007) and Lepidoptera (Iyengar et al. 2002). Finally, genetic association between genes involved in premating and postzygotic isolation, which can arise if both are on the Z chromosome, is expected to facilitate the evolution of reinforcement (Servedio and Sætre 2003; Hall and Kirkpatrick 2006). In Ficedula flycatchers, genes underlying female preference and low hybrid fitness are both Z-linked (Sætre et al. 2003). Thus, there are a variety of reasons for expecting that the Z chromosome might be involved in reproductive isolation, including prezygotic isolation and either intrinsic or extrinsic postzygotic isolation.
Two closely related bird species, the Thrush Nightingale (Luscinia luscinia) and the Common Nightingale (L. megarhynchos), are well suited for studying the role of the Z chromosome in speciation. These species are believed to have diverged in geographical isolation during the Pleistocene and have subsequently come into secondary contact in central and Eastern Europe where they form a narrow hybrid zone (Sorjonen 1986) (Fig. 1). Despite their overall morphological and ecological similarity, the two species can be distinguished by wing feather characteristics and slightly different songs (Cramp 1988). However, males with a mixed song are frequently observed in areas in which both species co-occur (Becker 2007). Although there is strong assortative mating between the species, mixed pairs as well as interspecific hybrids are occasionally observed in the region of sympatry (Becker 2007; Kverek et al. 2008). Morphological as well as genetic studies indicate that approximately 5% of birds in sympatry are hybrids (Becker 2007; R. Storchová et al. unpubl. ms.). Interestingly, interspecific pairs show slightly lower hatching performance compared to conspecific pairs, indicating that hybrids might be partially inviable (Becker 2007). Crosses between the two species in captivity have shown that, in accordance with Haldane’s rule, F1 hybrid females are sterile, whereas F1 hybrid males are fertile (Stadie 1991). Overall, these observations indicate that prezygotic isolation exists but is not complete, and that intrinsic postzygotic isolation also contributes to the reproductive barrier between the Common Nightingale and the Thrush Nightingale.
To investigate the role of the Z chromosome in reproductive isolation between nightingales, we examined patterns of polymorphisms and divergence at multiple autosomal and Z-linked loci in allopatric populations of both species. Two other pairs of recently diverged species of birds have been studied for patterns of differentiation at Z-linked and autosomal loci. Carling and Brumfield (2008) studied gene flow across a hybrid zone between Lazuli and Indigo Buntings in North America using four autosomal and two Z-linked markers, and showed that the two Z-linked markers had an average cline width of 309 km whereas the four autosomal markers had an average cline width of 466 km. Sætre and colleagues have studied Pied and Collared Flycatchers in Europe in both sympatric and allopatric populations and found that the Z chromosome shows greater differentiation than the autosomes, likely as a consequence of reduced introgression of the Z chromosome (Sætre et al. 2003) as well as positive selection at some Z-linked loci (Borge et al. 2005). Here, we expand on these studies by using DNA sequence data and an explicit divergence-with-gene-flow model to assess the contribution of gene flow and selection to observed patterns of differentiation at autosomal and Z-linked loci in nightingales.
Luscinia megarhychos and L. luscinia are sister species (Jønsson and Fjeldså, 2006) and they are not known to hybridize with other species (see the database of bird hybrids:www.bird-hybrids.com). Samples from 25 males of L. megarhynchos and 25 males of L. luscinia were analyzed in this study. Birds were caught during the breeding season (from late April until mid June) in 2006 and 2007. Two different populations from each species were sampled. According to their distance from the hybrid zone, we refer to these populations as close allopatry (Czech Republic and SW Poland for L. megarhynchos, NE Poland for L. luscinia) and distal allopatry (Spain for L. megarhynchos, Finland for L. luscinia) (Fig. 1). Within each population, birds were captured at different localities several kilometers apart. The exact positions of the localities are provided in Table A1. Additionally, one male Bluethroat (Luscinia svecica) was used as an outgroup. Luscinia svecica is the sister taxon to the clade containing L. megarhynchos and L. luscinia (Jønsson and Fjeldså, 2006). Genetic divergence at cytochrome b between L. svecica and either L. megarhynchos or L. luscinia is approximately 14% (Wink et al. 2002), which corresponds to a divergence time of 7 million years assuming approximately 2% sequence divergence per million years (Wilson et al. 1985). From each individual, a blood sample was collected by brachial vein puncture and stored in pure ethanol for further extraction of genomic DNA.
Genomic DNA was purified by DNeasy Tissue Kit (Qiagen) and used for PCR amplification and sequencing of eight autosomal and four Z-linked loci. Primer sequences were obtained from published studies (Borge et al. 2005; Backström et al. 2006, 2008). These primers were designed in conserved exonic regions of the chicken genome to amplify intronic sequences. PCR conditions as well as the length of PCR products are given in Table S1. In a few cases, nightingale-specific primers were designed from initial sequence to reduce nonspecific PCR amplification. Primer sequences are provided in Table S2. All PCR products were sequenced in both directions. The chromosomal location of each locus (i.e., whether autosomal or Z-linked) was determined on the basis of the conserved synteny between the chicken genome and the genome of passerines (Backström et al. 2006; Itoh et al. 2006; Nanda et al. 2008). We were also interested in knowing whether any of these loci were closely linked to each other. In the absence of a physical or genetic map for nightingales, we determined the chromosomal position of each locus in two bird species with known genome sequences, Gallus gallus and Taeniopygia guttata, using the Ensembl database (Table S3).We reasoned that if linkage relationships were similar in chickens and in zebra finches, they might be conserved in nightingales as well given the generally conserved karyotypes of birds (Price 2007). Most of the autosomal genes lie on different chromosomes both in the G. gallus and T. guttata genomes. Of the four Z-linked genes, PPWD1 and ADAMTS6 are close physically in both species (and thus probably also in nightingales), and therefore may not have independent evolutionary histories. However, the other two loci on the Z chromosome are separated from each other and from PPWD1 and ADAMTS6 by > 5 Mb in both species. The recombination rate on the Z chromosome is 1.43 cM/Mb in T. guttata (Pigozzi 2008) and 2.7 cM/Mb in G. gallus (Wahlberg et al. 2007). Thus, these loci are likely separated by at least 7 cM in Luscinia and are unlikely to be in linkage disequilibrium.
Sequences were edited using phred/phrap/consed/polyphred (Nickerson et al. 1997; Ewing and Green 1998; Ewing et al. 1998; Gordon et al. 1998) together with auxiliary shell scripts and Perl programs kindly provided by August Woerner. Alignments were generated by ClustalW (Thompson et al. 1994) as implemented in the program BioEdit (Hall 1999). All alignments were visually checked and manually adjusted. Exonic sequences as well as indel polymorphisms were excluded from the analyses. Individuals and/or positions with missing data were eliminated from the dataset. Sequences have been deposited in GenBank under accession numbers FJ695620-FJ696173.
Diploid sequences at each locus were separated into two haplotypes computationally using the program PHASE, version 2.1 (Stephens et al. 2001; Stephens and Donnelly 2003) with the following parameters: number of iterations = 10,000, thinning interval = 1, burnin = 1000 (see PHASE 2.1 documentation for details). We used the default recombination model, which is the general model for varying recombination rate from Li and Stephens (2003). For each dataset, we applied the algorithm five times with different random seeds, and we checked for consistency of results across independent runs. We obtained identical haplotypes across all runs for 11 of 12 loci. The only locus that showed slightly different haplotypes in different runs was 25149. We note that haplotype reconstruction can be affected by several demographic and evolutionary factors, including population subdivision, recent population expansion, or selection (Andrés et al. 2007). However, we find no evidence for population size changes or selection in these data (see Results).
Basic population genetic analyses of polymorphism, divergence, recombination, and tests of neutrality based on the allele frequency spectrum were performed with the programs DnaSP (Rozas et al. 2003) and SITES (Hey and Wakeley 1997). Genealogical relationships among haplotypes for each locus were reconstructed with the Neighbor-Joining method (Saitou and Nei 1987) using the program MEGA4 (Tamura et al. 2007). Hudson-Kreitman Aguadé (HKA) tests (Hudson et al. 1987) were performed using the HKA program (http://lifesci.rutgers.edu/~heylab/). The KST test for detecting geographic subdivision (Hudson et al. 1992) was performed using DnaSP.
We looked for evidence of gene flow by calculating the amount of linkage disequilibrium (LD) among specific classes of segregating sites as described in Machado et al. (2002) using the software SITES. D′ (Lewontin 1964) was used as a measure of LD. The test of Machado et al. (2002) is based on comparison of LD among pairs of shared polymorphisms (DSS), LD among pairs of exclusive polymorphisms (DXX), and LD among pairs of sites where one member is a shared and the other an exclusive polymorphism (DSX). If gene flow occurs after divergence, the difference, DSS − DXX or DSS − DSX, which is called x, is expected to be positive because polymorphisms that are shared due to common ancestry are expected to be older, on average, than polymorphisms that are shared due to postdivergence gene flow. Less LD is expected between older polymorphisms because recombination will have had more time to break down associations. Statistical significance of the observed x values was obtained by simulating datasets using the program WH (Wang et al. 1997; http://lifesci.rutgers.edu/~heylab/). This program simulates a model of isolation without gene flow. Because our dataset exceeded the maximum input for WH (which is 32 chromosomes per locus), we divided the dataset into three subsamples with the appropriate sample size and performed three independent simulations. Parameter estimates for all three subsamples are given in Table S4. The results of the LD test are shown only for x = DSS − DXX; however, the results based on x = DSS − DSX were similar.
The data were fit to the isolation-with-migration model (IM) using the programs IM and IMa (Nielsen and Wakeley 2001; Hey and Nielsen 2004, 2007). These programs estimate six parameters by Markov chain Monte Carlo simulations of genealogies. These parameters include the population-split time, the effective population size for the ancestral population and for the current populations, and the migration rates in both directions. Because IM and IMa assume no recombination within loci, we determined the longest region without observed recombination for each locus using the program IMgc (Woerner et al. 2007). This program removes either sites or haplotypes to produce the most information-rich contiguous DNA sequence segment that passes the four-gamete test. Nonrecombinant regions represented 79% of the length of each locus, on average, and were used as an input for the IM and IMa programs. We ran the programs separately for (1) all 12 loci, (2) the eight autosomal loci and (3) the four Z-linked loci. For each dataset, the programs were run three or four times with identical starting conditions, with the exception of the random number seed, to assess convergence. To facilitate mixing of the Markov chains we used Metropolis coupling with 15–20 chains and a geometric heating scheme. All runs began with a burn-in period of 100,000 steps and were allowed to continue for 1–10 millions steps. We were able to achieve adequate mixing of the Markov chains as indicated by trend line plots and effective sample size (ESS) values (for most runs, all ESS values were higher then 500 and no ESS value was lower than 50). Independent runs converged to the same result (e.g., maximum-likelihood estimates and marginal posterior probability distributions for the demographic parameters were essentially the same for all three runs; see Fig. S1). All estimated parameters of the IM model are scaled to the mutation rate. To convert these parameters to biologically more meaningful quantities (i.e., Ne, effective population size in number of individuals; m, migration rate per year; 2Nm, population migration rate; t, divergence time in years) we calculated the neutral mutation rate for each locus (μ) using divergence to the outgroup. The neutral mutation rate per year was calculated from the formula D = 2μt1 where D is the estimated Dxy (Nei 1987) between L. megarhynchos/luscinia and L. svecica for each locus and t1 is the divergence time. We assumed that the divergence time between L. megarhynchos or L. luscinia and L. svecica is 7 million years as estimated on the basis of cytochrome b sequence divergence (Wink et al. 2002). The geometric mean of the locus-specific mutation rates was then used for the parameter conversion. In analyses involving only Z-linked loci or only autosomal loci, the geometric mean of locus-specific mutation rates was calculated using only Z-linked or autosomal loci, respectively. To estimate mutation rates per generation, we need to know the generation time. Nightingales, like most passerines, first reproduce one year after hatching and thus a generation time of one year was assumed. This may be an underestimate because some individuals may reproduce over multiple years. Uncertainty in generation time will affect estimates of Ne, but will not affect estimates of 2Nm or t.
We sequenced eight autosomal and four Z-linked loci in close allopatric and distal allopatric populations of both species. For each locus, sequence from 38 to 50 chromosomes in L. megarhynchos and 34 to 50 chromosomes in L. luscinia was obtained (Table 1). The same loci were also sequenced in one male L. svecica, representing an outgroup. No significant genetic differentiation was revealed between close and distal allopatric populations within species for the vast majority of loci (KST ranged from − 0.1964 to 0.0335, P > 0.05). The only significant differentiation was detected for locus 25149 in L. megarhynchos (KST = 0.0454, P = 0.007) and for locus 25374 in L. luscinia (KST = 0.0580, P = 0.022). However, after applying a Bonferroni correction for multiple testing, population subdivision was not significant for any locus. Thus, we pooled data from close and distal allopatric populations in all subsequent analyses. The absence of strong geographic subdivision may be related to high dispersal rates due to migration.
Levels of nucleotide variation were quite high for both species, although slightly higher for L. luscinia than for L. megarhynchos (on autosomes, average π = 0.850% for L. luscinia and 0.796% for L. megarhynchos; average θ = 0.994% for L. luscinia and 0.679% for L. megarhynchos) (Table 1). This suggests a large effective population size for both species, consistent with their high abundances and large geographic ranges. Roughly similar values of nucleotide diversity have been observed, for example, in Drosophila (Moriyama and Powell 1996) and rabbits (Carneiro et al. 2009).
Nucleotide variation was lower on the Z chromosome than on the autosomes, even after we corrected for the different effective population size of the Z chromosome, which is three-fourth that of the autosomes assuming an equal sex ratio (after correction, the Z/A ratio for π was 0.33 in L. megarhynchos and 0.54 in L. luscinia; the Z/A ratio for θ was 0.56 in L. megarhynchos and 0.65 in L. luscinia). To test whether differences in mutation rate could account for reduced Z chromosome variation within species, we divided θ by divergence (Dxy) to L. svecica (which can by used as a surrogate for the mutation rate). The ratio θ/Dxy was still slightly lower for Z-linked than for autosomal loci (Z/A ratio was 0.58 in L. megarhynchos, 0.72 in L. luscinia).We then applied a 12-locus HKA test as a more formal approach to asking if the Z chromosome shows reduced diversity relative to the autosomes. This test takes into account differences in mutation rate (by using an outgroup) as well as differences in the effective population size of the Z and autosomes. The test was applied separately to each species and to both species together, in each case using a single sequence from L. svecica as an outgroup. None of these three tests rejected a neutral model (P > 0.10; 10,000 permutations). We also performed two-locus HKA tests with the pooled data from the autosomes as one locus and the pooled data from the Z chromosome as a second locus to test for differences in nucleotide variation between the autosomes and the Z chromosome, without testing for variation among individual loci. This was done for each species separately as well as both species together. None of these three tests were significant (P > 0.10; 10,000 permutations) (Table S5). Together, these results suggest that there is only a modest (nonsignificant) reduction of variation on the Z chromosome compared to the autosomes.
We also looked for evidence of selection by comparing the distribution of allele frequencies with the expectations under a neutral equilibrium model using Tajima’s D (TD) (Table 1). None of the 24 values of TD within species were significantly different from the neutral expectation of 0. There was a very slight tendency toward negative values of TD (17 of these 24 values were negative), and this may reflect mild population expansions.
Levels of differentiation between the two species of nightingales are summarized in Table 2. We observed a striking contrast in the level of differentiation between the Z chromosome and the autosomes. The mean FST as well as net divergence (Da) were consistently and significantly higher for Z-linked loci (FST = 0.783, Da = 0.976%) than for autosomal loci (FST = 0.298, Da = 0.327%) (Mann–Whitney U-test, P < 0.01). A similar contrast was found in the number of fixed differences and shared polymorphisms. In total, the Z-linked loci showed 14 fixed differences and only one shared polymorphism. On the other hand, the autosomal loci showed 0 fixed differences and 52 shared polymorphisms (Fisher exact test, P < 0.0001).
Another way of examining differentiation is to look at gene genealogies. Figure 2 shows neighbor-joining trees for each locus. Despite generally low bootstrap values, the genealogies suggest that lineage sorting has been completed on the Z chromosome, with the exception of the PPWD1 locus. A closer inspection of the sequences of PPWD1 gene revealed that one haplotype of L. luscinia seems to be introgressed into L. megarhynchos. The length of this introgressed haplotype is 344bp and includes six polymorphic sites. Because the length of this introgression is a little bit shorter than the length of the whole sequence, the neighbor-joining tree shows this sequence somewhere between L. luscinia and L. megarhynchos. In contrast, the autosomal loci did not show clear phylogenetic differentiation between species. The lack of concordance between genealogical relationships and species identity for autosomal loci might indicate that lineage sorting is not yet complete or that interspecific gene flow has occurred.
To reveal whether interspecific gene flow occurred after the divergence of the two species, we performed two tests. The first test is based on patterns of LD within genes (Machado et al. 2002). Under a scenario of postdivergence gene flow, observed values of the x statistic are expected to be more positive than simulated values (see Materials and Methods). The observed values of x were consistently higher than the simulated values in both species, and in four of five genes, x was significant in at least one run (Table 3). These results suggest that gene flow has occurred after the divergence of both species.
To further explore patterns of interspecific gene flow and to estimate divergence time as well as effective population sizes for both species, we fitted an IM model to our data using the IMa program. The program was run separately for (1) all 12 loci, (2) the eight autosomal loci and (3) the four Z-linked loci. Maximum-likelihood estimates (MLE) and 90% highest posterior density (HPD) intervals of all model parameters are given in Table 4, and posterior distributions are shown in Figure 3. For the whole dataset, the estimated Ne was 448,584 for L. megarhynchos, 853,268 for L. luscinia, and 220,888 for the ancestral population. Similar estimates were obtained for autosomal and Z-linked loci, although estimates of Ne for Z-linked genes were slightly lower. Ne of the ancestral population could not be estimated accurately for the Z-linked loci due to the absence of shared polymorphisms.
We estimated the time since divergence between the two species to be 1.8 Mya for the whole dataset. Interestingly, the estimated divergence time was two times higher for Z-linked loci (2.6 Mya) than for autosomal loci (1.3 Mya). This discrepancy might be partly explained by an underestimation of Ne of the ancestral population for Z-linked loci due to the absence of shared polymorphisms. However, because estimates of Ne of the ancestral population from autosomal loci are still quite low, the interesting possibility remains that reproductive isolation and thus species divergence occurred earlier on the Z chromosome than on the autosomes.
Both migration rate parameters revealed evidence of gene flow between the species. Gene flow was generally higher from L. megarhynchos to L. luscinia (2N2m2 = 0.326), than in the opposite direction (2N1m1 = 0.119). The estimates of gene flow for autosomal loci were quite high in both directions. In contrast, no gene flow was detected on the Z chromosome from L. megarhynchos to L. luscinia and limited gene flow was detected in the opposite direction (corresponding to the introgression of only one haplotype at PPWD1; Fig. 2).
To test whether models with gene flow are a significantly better fit to the data than models without gene flow, we performed log-likelihood ratio tests of nested models (Hey and Nielsen 2007) as implemented in the IMa program, where migration parameters, m1 and/or m2, were set to zero (Table 5). For the whole dataset, a model with gene flow was a significantly better fit to the data than all models without gene flow. For autosomal loci, a model with gene flow was significantly better than a model without gene flow in both directions. It was also significantly better than a model in which gene flow occurs only from L. luscinia into L. megarhynchos, however, it was not significantly better than a model with gene flow only in the opposite direction. In contrast, for Z-linked loci, a model with gene flow was not a significantly better fit to the data than any model without gene flow. Finally, we tested whether a model with asymmetrical introgression is better than models where gene flow is the same in both directions. This test, however, did not yield significant results in any of the comparisons (Table 5).
We also estimated migration rates separately for each locus using the IM program (Table 6). Not all of these estimates were recovered with high confidence, either because the tails of the posterior probability densities did not reach zero or because the 90% HPD intervals were quite large (the posterior probability distributions for locus-specific migration rates are shown in Fig. S2). Nonetheless, the estimates were remarkably consistent with classic estimators of gene flow (FST-based estimates of Nm), suggesting that introgression is the major cause of the contrasting patterns of differentiation between the Z chromosome and the autosomes. Average migration rates for Z-linked loci were significantly lower than for autosomal loci (Mann-Whitney U test, P = 0.021857). This difference was caused mainly by introgression from L. megarhynchos into L. luscinia (Table 6).
To study the role of the Z chromosome in speciation, we analyzed patterns of polymorphism and divergence at eight autosomal and four Z-linked loci in allopatric populations of two nightingale species. Four main results stem from our study. (1) We estimated that the two species diverged approximately 1.8 Mya and have large effective population sizes (L. megarhynchos Ne = 448,584, L. luscinia Ne = 853,268). (2) We observed strikingly different patterns of polymorphism and divergence between the Z chromosome and the autosomes. Z-linked loci showed slightly, but not significantly, lower levels of variation within species compared to autosomes. On the other hand, the Z chromosome was characterized by significantly higher levels of differentiation between species than were the autosomes. This was seen in levels of FST, net nucleotide divergence, and in the ratio of fixed differences to shared polymorphisms, all of which were higher for Z-linked loci than for autosomal loci. (3) We found strong evidence of gene flow following divergence of the two species. This gene flow occurred in both directions, although more introgression was from L. megarhynchos into L. luscinia. (4) Gene flow was significantly higher on the autosomes than on the Z chromosome. Below, we discuss our results in light of recent theories of speciation, with respect to the large effect of the sex chromosomes in reproductive isolation and adaptive evolution.
The Common Nightingale and the Thrush Nightingale are thought to have diverged in allopatry during the Pleistocene climatic oscillations (Sorjonen 1986) which started approximately 2.4–2.6 Mya (Hewitt 2000; Provan and Bennett 2008). Our estimated divergence time of 1.8 Mya is consistent with this idea. However, it is important to bear in mind that estimates of divergence time depend on estimates of mutation rate, which should be considered approximate owing to the uncertainty of divergence time to the outgroup.
Interestingly, previous estimates based on mitochondrial cytochrome b sequences were considerably higher (3.3 Mya; Wink et al. 2002). The high estimated divergence times from mitochondrial sequences in nightingales and other species were used to question the importance of Pleistocene climatic oscillations in the origin of many sister species of birds (Klicka and Zink 1997). However, single-locus estimates of age can be unreliable, and if mitochondrial divergence times have generally been overestimated, as seems to have been the case in nightingales, we may have to re-evaluate the importance of ice ages in bird speciation (Weir and Schluter 2004).
Interestingly, our estimate of divergence time was much higher for the Z chromosome (2.6 Mya) than for the autosomes (1.3 Mya). This might indicate that the two species did not diverge in strict allopatry, but exchanged genes during the speciation process either in sympatry or in parapatry. Under such a scenario, divergence times might vary across the genome, because interspecific gene flow is expected to be restricted earlier in genomic regions harboring loci causing hybrid incompatibilities than in genomic regions without such loci (Wu 2001). Indeed, we found strong evidence of post divergence gene flow between the nightingale species. This gene flow was significantly higher on the autosomes compared to the Z chromosome, which is consistent with lower estimates of divergence time for the autosomes. Glacial fluctuations were a major factor in shaping bird distributions during the Quaternary. It is possible that nightingales underwent repeated cycles of allopatry followed by periods of gene flow during interglacial events.
Although we do not know the exact position of the four analyzed Z-linked genes in nightingales, PPWD1 and ADAMTS6 are physically close to each other in both the chicken and in the zebra finch, and thus they may also be physically close in nightingales. The other two Z-linked loci are farther apart in both chicken and zebra finch, and thus may lie in different parts of the Z chromosome in nightingales (Table S3). Loci that are closely linked may have correlated histories, whereas loci that are far apart can provide independent information about the levels of divergence and gene flow. Genetic linkage can be caused not only by close physical proximity of genes on the chromosome, but also by chromosomal rearrangements, such as inversions. We cannot exclude the possibility that some chromosomal rearrangements differ between the nightingale species and might result in linkage of multiple Z-linked loci across a large part of the Z chromosome. However, this seems unlikely given the low frequency of karyotypic changes in birds (Price 2007) and the young age of these nightingale species.
The sex chromosomes are known to play a prominent role in reproductive isolation in species with heterogametic males (e.g., Coyne and Orr 1989; Masly and Presgraves 2007). Here, we show that in species with heterogametic females, the Z chromosome is also likely to be important in reproductive isolation. We found that Z-linked loci systematically showed reduced interspecific introgression compared to autosomal loci. Similar results have been demonstrated in two other bird systems. Genetic analysis of a sympatric population of two naturally hybridizing Ficedula flycatcher species revealed that interspecific gene flow occurs at autosomal, but not at Z-linked loci (Sætre et al. 2003). Recently, reduced introgression of the Z chromosome has also been observed across the hybrid zone between Indigo and Lazuli Bunting (Carling and Brumfield 2008, 2009). Involvement of the Z chromosome in reproductive isolation has also been found between different species of Heliconius butterflies (Jiggins et al. 2001; Naisbit et al. 2002). Together, these results indicate that the importance of the Z chromosome in speciation might represent a general phenomenon.
The Z chromosome could play a large role in postzygotic as well as prezygotic isolation (Qvarnström and Bailey 2008). The Z chromosome could underlie postzygotic barriers, such as hybrid sterility and inviability, due to the exposure of recessive mutations in the hemizygous individuals or due to a higher density of genes contributing to reproductive isolation on the Z. The higher density of isolation genes on Z may be caused by faster Z evolution, accumulation of meiotic drive elements on the Z or disruption of meiotic Z chromosome inactivation that has been recently reported in birds (Schoenmakers et al. 2009). The meiotic drive hypothesis for the large-Z effect may be unlikely because sex ratio distorters are less likely to evolve on the Z chromosome compared to the X chromosome because male-biased populations are especially prone to extinction (Tao et al. 2007). All these mechanisms lead to preferential occurrence of hybrid sterility or inviability in the heterogametic sex (Coyne and Orr 2004). Nightingales (Stadie 1991), Ficedula flycatchers (Gelter et al. 1992) as well as Heliconius butterflies (Jiggins et al. 2001) all show female-limited hybrid sterility, supporting the idea that the Z chromosome is involved in the evolution of intrinsic postzygotic isolation.
On the other hand, the Z chromosome could promote the evolution of prezygotic isolation through sexual selection and/or reinforcement. Sexual selection may facilitate speciation because it can cause rapid divergence of male-mating signals and female preferences (Ritchie 2007). According to the theory of reinforcement, natural selection against unfit hybrids could enhance prezygotic isolation in areas in which two species coexist (Noor 1999). Both theoretical and empirical studies suggest that sexual selection as well as reinforcement are facilitated in species with female heterogamety if genes for male traits and female preferences are Z-linked (Hastings 1994; Reeve and Pfennig 2003; Servedio and Sætre 2003; Kirkpatrick and Hall 2004; Albert and Otto 2005; Hall and Kirkpatrick 2006). Prezygotic isolation plays an important role in nightingale speciation as shown by strong assortative mating between species and the low frequency of interspecific hybrids. Both the Common Nightingale and the Thrush Nightingale have dull brown plumage without any differences between the sexes. Sexual selection based on plumage color is thus unlikely to contribute to speciation. However, nightingales are known for their outstanding song, which differs slightly between the species, and might be under sexual selection (Cramp 1988). If reinforcement is responsible for the evolution of prezygotic isolation, we would expect more pronounced interspecific differences in male-mating signals in sympatry than in allopatry. This has been shown for example in Ficedula flycatchers, where sympatric males show increased differences in plumage compared to allopatric males (Sætre et al. 1997). Empirical observations in nightingales, however, show just the opposite pattern; songs of both nightingale species tend to converge in sympatry (Sorjonen 1986; Becker 2007). This indicates that, contrary to the situation in Ficedula flycatchers, reinforcement is unlikely to contribute to the evolution of prezygotic isolation in nightingales. It remains unclear, however, whether divergence of male song and female preferences in allopatry might contribute to reduced introgression of the Z chromosome.
The greater differentiation observed on the Z chromosome relative to the autosomes could, in principle, be explained by reduced introgression on the Z (e.g., Sætre et al. 2003) or positive selection that increases interspecific divergence by eliminating shared variants (e.g., Borge et al. 2005).We sought to disentangle these two effects by using tests for positive selection and by using several models that explicitly test for postdivergence gene flow. Neither the HKA test (Hudson et al. 1987) nor tests based on the allele frequency spectrum (Tajima 1989) provided evidence for positive selection on the Z. In contrast, substantial gene flow on the autosomes but not on the Z was revealed by divergence-with-gene-flow models (Nielsen and Wakeley 2001; Hey and Nielsen 2004, 2007). Thus, these data provide no evidence for recent strong selection on the Z, although we cannot rule out weak or older selection. On the other hand, the absence of gene flow on the Z in the face of substantial gene flow on the autosomes is consistent with the presence of Z-linked genes contributing to reproductive isolation. The fact that these species are separated by hybrid female sterility and show Z-linked genes contributing to isolation is analogous to X-linked male sterility in species with heterogametic males. This highlights the similarities in the genetic basis of reproductive isolation between organisms with heterogametic males and organisms with heterogametic females during the early stages of speciation.
We are grateful to M. Ahola, M. Antczak, D. B. Campas, A. Cabrera, P. A. Del Baño, P. Dołata, K. L. Evans, P. Kverek, T. Polo Aparisi, and J. Vokurková for help with the sample collection. V. Pavel kindly provided the sample from the Bluethroat. We thank members of the Nachman lab, especially A. Geraldes and M. Carneiro, for many useful discussions. J. Good, A. Geraldes, M. Sans, D. Storch, J. Lindell, and two anonymous reviewers provided critical comments on early versions of the manuscript. The research was supported by junior grant GAAV (KJB501110701) and post-doctoral grant GACR (206/08/P160) to RS, by the Ministry of Education, Youth and Sport of the Czech Republic (MSM0021620828), by a Fulbright fellowship to RS, and by NSF and NIH grants to MWN.
Sequence data from this article have been deposited with the GenBank Data Libraries under accession nos. FJ695620-FJ696173.
The following supporting information is available for this article:
Figure S1: Marginal posterior probability of model parameters (scaled by the neutral mutation rate) for all loci (All), autosomal loci (Chr A), and Z-linked loci (Chr Z).
Figure S2: Marginal posterior probability distributions for the locus-specific migration rate estimates (scaled by the neutral mutation rate).
Table S1: Primers and PCR conditions.
Table S2: Nightingale-specific primers.
Table S3: Ensembl chromosomal position of best blast hit for eight autosomal and four Z-linked genes.
Table S4: Isolation model fitting.
Table S5: HKA test of positive selection.
Supporting Information may be found in the online version of this article.
(This link will take you to the article abstract).
Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting informations supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.