Multilocus Analysis of Variation and Speciation in the Closely Related Species Arabidopsis halleri and A. lyrata

Corresponding author: Barbara E. Stranger, Facultat de Biologia, Universitat de Barcelona, Diagonal 645, 08028 Barcelona, Spain., bstrange{at}porthos.bio.ub.es (E-mail)

Communicating editor: O. SAVOLAINEN

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Nucleotide variation in eight effectively unlinked genes was surveyed in species-wide samples of the closely related outbreeding species Arabidopsis halleri and A. lyrata ssp. petraea and in three of these genes in A. lyrata ssp. lyrata and A. thaliana. Significant genetic differentiation was observed more frequently in A. l. petraea than in A. halleri. Average estimates of nucleotide variation were highest in A. l. petraea and lowest in A. l. lyrata, reflecting differences among species in effective population size. The low level of variation in A. l. lyrata is concordant with a bottleneck effect associated with its origin. The A. halleri/A. l. petraea speciation process was studied, considering the orthologous sequences of an outgroup species (A. thaliana). The high number of ancestral mutations relative to exclusive polymorphisms detected in A. halleri and A. l. petraea, the significant results of the multilocus Fay and Wu H tests, and haplotype sharing between the species indicate introgression subsequent to speciation. Average among-population variation in A. halleri and A. l. petraea was 1.5- and 3-fold higher than that in the inbreeder A. thaliana. The detected reduction of variation in A. thaliana is less than that expected from differences in mating system alone, and therefore from selective processes related to differences in the effective recombination rate, but could be explained by differences in population structure.

MULTILOCUS analyses of nucleotide variation in closely related species not only reveal variation that has accumulated independently in each species since their divergence, but also may detect variation originally segregating in their common ancestor, which has not yet been lost. Indeed, early in the process of speciation, newly formed species share polymorphisms due to recent divergence from a common ancestor (CLARK 1997 ; WAKELEY and HEY 1997 ). As time goes on and species diverge further, polymorphisms that were once shared are lost or fixed due to selection or genetic drift, and new mutations arise independently in each species. The effects of these processes are observed late in speciation as the loss of shared variants between species and the accumulation of fixed differences between species. Between these two extremes, there is an intermediate stage where some loci share polymorphisms, while other loci exhibit fixed differences between species.

Speciation models have been developed to describe the expected behavior of genetic polymorphism within and between species following divergence (WAKELEY and HEY 1997 ; KLIMAN et al. 2000 ; MACHADO et al. 2002 ). In the isolation speciation model (HEY and KLIMAN 1993 ; HILTON et al. 1994 ), an ancestral population splits into two descendant species with no subsequent gene flow. Population and speciation parameters can be inferred using this model from the observed numbers of shared and exclusive polymorphisms and fixed differences (WAKELEY and HEY 1997 ; WANG et al. 1997 ). This type of analysis has been extensively applied to data from closely related Drosophila species (e.g., HILTON et al. 1994 ; KLIMAN et al. 2000 ) and has recently been used in tomato (BAUDRY et al. 2001 ) and in maize species (TIFFIN and GAUT 2001 ). Here, we use divergence population genetics (KLIMAN et al. 2000 ) to study speciation in closely related species of Arabidopsis.

Arabidopsis halleri and A. lyrata are closely related outcrossing species whose divergence from the inbreeder A. thaliana is estimated to have occurred 5 MYA (KOCH et al. 2000 ). A. halleri ssp. halleri (henceforth A. halleri) is present in Central Europe, whereas A. lyrata ssp. petraea (henceforth A. l. petraea) is patchily distributed in Northern and Central Europe, and A. lyrata ssp. lyrata (henceforth A. l. lyrata) is found only in North America (MITCHELL-OLDS 2001 ; MITCHELL-OLDS and CLAUSS 2002 ). A. halleri and A. l. petraea distribution areas overlap in some Central European areas (HEGI 1986 ) and have been shown to produce semifertile hybrids in the laboratory (MACNAIR et al. 1999 ). On the other hand, A. l. petraea and A. l. lyrata are two subspecies of A. lyrata that are completely interfertile but have an allopatric distribution. Here, we study nucleotide variation in eight unlinked genes in species-wide samples of the closely related species A. halleri and A. l. petraea and in three of these genes in A. l. lyrata. We aim to establish whether the species A. halleri and A. l. petraea (and also the subspecies A. l. petraea and A. l. lyrata) have remained isolated since their split.

This work also aims to quantify species-wide levels of variation in the studied outbreeding Arabidopsis species and to establish, through comparison with A. thaliana, the effect of mating system and population structure on intraspecific variation. In mutation-drift equilibrium, differing levels of variation might indicate differences in either the effective population size or the neutral mutation rate. They might also reflect differences in mating system, which influence effective population size (POLLACK 1987 ; NORDBORG and DONNELLY 1997 ) and overall rates of recombination (GOLDING and STROBECK 1980 ). The reduced effective recombination in species with a high level of inbreeding is expected to increase the effects of genetic hitchhiking accompanying selective sweeps (MAYNARD SMITH and HAIGH 1974 ; KAPLAN et al. 1989 ) and background selection (CHARLESWORTH et al. 1993 ) and may greatly influence genetic variability in inbreeding species (NORDBORG et al. 1996 ). Metapopulation dynamics also have important effects on levels of genetic variation, and variability levels may even be quite similar in selfing and outcrossing species (INGVARSSON 2002 ). Although some studies have compared some of these species using microsatellite data (e.g., VAN TREUREN et al. 1997 ; WESTMAN and KRESOVICH 1998 ; CLAUSS et al. 2002 ), to date few studies have used intraspecific DNA variation to investigate the effects of mating system in plants. In some comparisons, they have revealed the expected reduction of species-wide nucleotide polymorphism in inbreeders relative to outbreeders (LIU et al. 1998 , LIU et al. 1999 ; BAUDRY et al. 2001 ; CLAUSS and MITCHELL-OLDS 2003 ; WRIGHT et al. 2003 ), but not in others (SAVOLAINEN et al. 2000 ).

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Species sampling:
Table 1 describes the geographic origin of the individuals of each species used in this study and the genes sequenced for each individual. DNAs of most A. halleri and A. l. petraea individuals were kindly provided by Pierre Saumitou-Laprade and Helmi Kuittinen, respectively. In other cases, seeds were obtained from T. Mitchell-Olds' personal collection (some A. halleri and A. l. petraea individuals and all A. l. lyrata individuals) and from the Nottingham Arabidopsis Stock Centre (A. thaliana ecotypes Col-0, Mrk-0, Kas-1, Lip-0, Ler, and Ll-0; Nottingham, United Kingdom). Seeds were planted and grown in growth chambers under standard 16-hr light and 8-hr dark conditions. Leaves were harvested from individual plants prior to flowering, and genomic DNA was extracted as described (ROGERS and BENDICH 1985 ; AUSUBEL et al. 1997 ).

PCR amplification and sequencing:
The genes studied code for enzymes in the phenylpropanoid pathway (CAD, cinnamoyl alcohol dehydrogenase; CHI, chalcone isomerase; CHS, chalcone synthase; DFR, dihydroflavonol reductase; F3H, flavanone-3-hydroxylase; and FAH-1, ferulate-5-hydroxylase), which may influence resistance to UV-radiation, and enzymes involved in glucosinolate biosynthesis (GS, glucosyltransferase; and MAM-L, methylthioalkylmalate-like synthase), and thus in plant defense against insects. The PCR amplification primers used, amplification conditions, and description of the fragment amplified for each of the eight genes studied are as described in KUITTINEN et al. 2002 . In those cases where the length of the amplified fragment was >600 bp, additional sequencing primers spaced at 400-bp intervals were designed.

In A. thaliana, sequences (one per ecotype) were obtained on both strands by direct sequencing of PCR products. In the other species, two sequences were obtained per individual. For the CHS, GS, and MAM-L loci, the PCR products were run on agarose gels, the appropriate-sized bands were excised, and the DNA was extracted from the gel slices using the QIAquick gel extraction kit (QIAGEN, Valencia, CA). The purified fragments were subsequently cloned into a pCRII vector using the TOPO-TA cloning system (Invitrogen, Carlsbad, CA). After transformation of TOP10F bacteria, plasmid DNA was extracted from individual transformed colonies using a plasmid DNA isolation kit (QIAGEN). The M13 universal and M13 reverse primers, as well as internal locus-specific primers, were used for sequencing. Both strands of 4–18 plasmids per locus per plant were sequenced using the Big Dye Terminator 2.0 kit (Applied Biosystems, Foster City, CA) and run on an ABI 3700 sequencer (Applied Biosystems). For the CAD, CHI, DFR, F3H, and FAH1 loci, sequences were obtained by direct sequencing of PCR products after their purification with QIAquick columns (QIAGEN); in all lines with two or more polymorphisms, the cloning and sequencing strategy indicated above was used.

Data analyses:
Contigs were assembled using Sequencher (Gene Codes, Ann Arbor, MI), SeqMan (DNASTAR, Madison, WI), and SeqEd version 1.0.3 program (Applied Biosystems). Multiple sequences for each gene were aligned using either a combination of methods implemented in Megalign (DNASTAR) or ClustalW (THOMPSON et al. 1994 ), with some additional manual alignment. Sequences were edited with the MacClade (version 3.06 or 4.0) program (MADDISON and MADDISON 1992 , MADDISON and MADDISON 2000 ).

The DnaSP version 3.98 (ROZAS and ROZAS 1999 ) and the SITES (HEY and WAKELEY 1997 ) software packages were used for most analyses of DNA polymorphism and to perform some neutrality tests. The program HKA (distributed by Jody Hey through http://lifesci.rutgers.edu/heylab) was used to perform multilocus HKA tests and to estimate population parameters. The program MEGA version 2.0 (KUMAR et al. 1994 ) was used for phylogenetic reconstruction. The program WH (as described in WANG et al. 1997 ) was used to test the isolation model of speciation.

Nucleotide variation was estimated as nucleotide diversity, (NEI 1987 ). The level of variation between locations (_B), which refers to the average number of pairwise differences between geographic locations, was used as a measure of species-wide diversity (as proposed by WAKELEY and ALIACAR 2001 ). Nucleotide diversity within locations (_w) was also estimated. Note that our _B does not refer to net variation between locations, i.e., to _T - _w (NEI 1987 ). Synonymous and nonsynonymous variation was estimated according to NEI and GOJOBORI 1986 . Insertion/deletion polymorphisms were excluded from the analyses. Hudson's S_nn statistic was used to test for population differentiation (HUDSON 2000 ).

Several tests were used to detect possible departure from the predictions of an equilibrium neutral model: tests based on either the frequency spectrum of polymorphisms or the haplotype distribution (TAJIMA 1989 ; FU and LI 1993 ; FU 1996 ; FAY and WU 2000 ) and tests based on the relationship between intraspecific and interspecific variation (HUDSON et al. 1987 ; MCDONALD and KREITMAN 1991 ). Segregating sites that included multiple hits were excluded from the calculation of the neutrality tests.

For the Fay and Wu H test (FAY and WU 2000 ), coalescent simulations that included recurrent mutation were performed, which results in a more conservative test. The coalescent simulations were done using the level of variation () rather than a fixed number of mutations. Analyses were performed with total and with silent positions, using both an unbiased and a biased transition/transversion ratio (ratio values equal to 0.5 and 2.0, respectively). Time of divergence to the outgroup was estimated taking into account the ratio of divergence to nucleotide diversity (K/). Therefore, the time to isolation (T_d) of the outgroup sequence was calculated as ((K/) - 1)/2 (HUDSON et al. 1987 ), where T_d is expressed in 4N_e generations and N_e is the effective population size of the studied population.

The nonparametric sign test (SOKAL and ROHLF 1995 ) was performed to detect positive (or negative) tendencies across loci, for a given statistic. After performing coalescent simulations for each locus, conditioned on the number of segregating sites, the median value of the test statistic under consideration (e.g., Tajima's D) was obtained for each locus. The observed number of loci with positive or negative values can be compared to the null hypothesis, where positive or negative values both have probability (p) of ¹/₂ under neutrality. The multilocus P value was then obtained by calculating the cumulative binomial under the assumption of p = 0.5.

Sampling strategy for neutrality tests:
The null hypothesis of most tests of neutrality assumes panmixia in addition to mutation-drift equilibrium. Deviations from panmixia, e.g., population subdivision, can thus lead to rejection of the null hypothesis even if neutrality holds. Since sequences in the present study were obtained by sampling within and between locations throughout each species' distribution area, we tried to minimize the effect of population subdivision in the neutrality tests. If no genetic differentiation was detected for a particular locus and species, the sample was considered to be panmictic, and all sequences were used in neutrality tests. Otherwise, a single randomly chosen sequence per location (constructed sample) was used, and all (or a subset of all) combinations were tested. This mixed sampling strategy uses all the information available whenever there is no evidence for population subdivision and minimizes the effect of subdivision on neutrality tests if subdivision is present. The use of a single sequence per deme conforms to the assumptions of a standard equilibrium model, where = 4N_eµ and N_e is rescaled as a function of population structure (WAKELEY and ALIACAR 2001 ). To determine whether migration among demes is correlated with geographic distance, a MANTEL 1967 test was used. Tests were performed using the number of nucleotide differences between samples and the corresponding geographic distance measured in kilometers. Geographic distances were obtained with Chris Michels' program (http://jan.ucc.nau.edu/cvm/latlongdist.html). The Mantel tests were performed with the zt version 1.0. program (BONNET and VAN DE PEER 2002 ).

Multilocus tests of neutrality:
Neutrality tests using multilocus data from a single species were performed following Hey (as detailed in the HKA program) for TAJIMA's (1989) D_T, FU's (1996) F_s, FU and LI's (1993) D_F and F_F, and FAY and WU's (2000) H tests. When, as in this example, loci are either not linked or loosely linked, then the history of each locus can be considered to be independent and forces affecting one particular locus should not affect other loci. Under strict neutrality and complete independence between loci, the effective population size (N_e) would be the only population parameter shared by all loci. The multilocus test statistic (m[D_T], m[F_s], m[D_F], m[F_F], and m[H]) was obtained, in each case, by summing across loci the corresponding test statistic values for each locus and dividing by the total number of loci, which gives them all the same weight. The empirical distribution of the multilocus test statistic under neutrality was obtained from coalescent-based simulations. The simulations were run independently for each locus (with the corresponding sample size) and were conditioned on the number of segregating sites. For each replicate, the multilocus test statistic (e.g., m[D_T]) was obtained by averaging the corresponding test statistic (e.g., D_T) across loci.

Inferring ancestral variation:
The multilocus sequence comparison of two closely related species allows detection of four classes of mutations (excluding sites with multiple hits, S_mhits; Fig 1A), which can be grouped into three different categories (WAKELEY and HEY 1997 ): (i) segregating sites where mutations are exclusive to one species (S_x1, S_x2), (ii) fixed differences between species (S_f), and (iii) shared polymorphisms (S_s). Polymorphisms in this latter category can either reflect ancestral polymorphisms or be the result of coincident mutations. The SITES (HEY and WAKELEY 1997 ) software package was used to obtain the observed numbers of mutations in each class.

View larger version (24K): In this window In a new window Download PPT slide   Figure 1. (A) Diagram representing two different genealogies of a sample from species 1 and 2. Dots indicate the different classes of single mutation events that occurred in the past; Sx1, polymorphism exclusive to species 1; Sx2, polymorphism exclusive to species 2; Sf, mutation fixed between species 1 and 2; Ss, polymorphisms shared between species. (B) Similar diagram but including an outgroup species in the genealogies, which allows identification of new classes of mutations: Sf1, Sf2, and Sfout represent mutations fixed in species 1, 2, and in the outgroup, respectively; Sf1x2, mutation fixed in species 1 (and therefore different from the outgroup) but polymorphic in species 2; Sf2x1, mutation fixed in species 2 (and therefore different from the outgroup) but polymorphic in species 1 (not represented in the diagram).

As depicted in Fig 1B, availability of the homologous sequences from an outgroup species allows partitioning of fixed differences into the two lineages of closely related species (S_f1, S_f2) and also detection of two new classes of mutations (S_f1x2, S_f2x1). These latter classes correspond to mutations that are fixed in one species but still segregate as polymorphisms in the other species. Unlike biallelic variation, multiple hits are grouped into a single category because of the difficulty of unambiguously establishing the direction of mutations. These two new classes of mutations and the shared polymorphisms class constitute, therefore, the third category of mutations, hereafter called ancestral mutations. The numbers of mutations in each class were obtained using the calc_mut program (available from the authors).

The hypergeometric distribution (ROZAS and AGUADE 1993 ; CLARK 1997 ) was used to obtain the probabilities that the observed numbers of ancestral mutations (S_s, S_f1x2, and S_f2x1) were due to recurrent mutation. Probabilities were obtained separately for each class of putatively ancestral mutations and for each locus, considering all positions and only silent positions. The same reasoning as in CLARK 1997 was used to obtain expected numbers of shared mutations, with the following reasoning applied to the other two classes: a unique mutation can be assigned to class S_f1x2 (or S_f2x1; Fig 1B) if it were an ancestral mutation, while it could have two possible origins (one segregating mutation in species 2 and one fixed mutation in the outgroup species, or one segregating mutation in species 2 and a fixed mutation in species 1) if it had originated by recurrent mutation. Then, the number of recurrent mutations is the number expected by chance in the categories formed by S_x2 + S_f1x2 vs. S_fout + S_f1 + S_f1x2 given the total number of nucleotide sites in the locus. The probability of having the observed number of S_f1x2 mutations is then obtained using the hypergeometric distribution. Loci with less than five S_f1x2 (and S_f2x1) mutations were considered together.

Population parameters for two extant closely related species were estimated from the observed distribution of variation within and between species (S_x1, S_x2, S_f, S_s) under the simple isolation model of speciation (here named WH) described in WANG et al. 1997 . This model postulates a split, at some time in the past, of a panmictic ancestral population into two descendant populations, with no subsequent gene flow between the descendant populations. Under the assumption that variation is neutral and that mutations occur according to the infinite-sites model, diversity in the ancestral (_A) and each descendant (₁ and ₂) species, as well as the time from speciation (T_sp, measured in N_e generations of one of the descendant species), can be estimated.

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Levels of DNA sequence polymorphism and divergence from A. thaliana:
In both A. halleri and A. l. petraea, sequences of eight independent genes (Table 2) were obtained from one to two individuals from each of three to five populations (Table 1 and Table 3). A subset of these genes was sequenced in samples from A. l. lyrata and A. thaliana (Table 1). Two sequences per individual were obtained for each locus, except in A. thaliana (see MATERIALS AND METHODS). Samples consisted of individuals collected across each species' distribution area. In general, each location was represented by one individual, although for A. halleri two individuals from the same location were sequenced in some cases (Table 1). For each species, the length of the region analyzed for each gene varied between 500 and 2000 nucleotides; these numbers are, however, lower when all species are considered (Table 2). Figures S1–S8 (at http://www.genetics.org/supplemental/) summarize intra- and interspecific variation.

For each gene and species, nucleotide diversity estimates for total and silent variation are presented in Table 3 and Table 4, while estimates of synonymous and nonsynonymous variation are given in Table 5. In A. halleri, A. l. petraea, and A. l. lyrata, estimates of variation within and among locations were obtained in addition to estimates based on all sequences (Table 3). Comparison of the averages of among-population variation between the three outcrossing species yielded similar results to those obtained using all sequences: (i) variation was higher in A. l. petraea than in A. halleri, and (ii) variation was lowest in A. l. lyrata (Table 3). Levels of polymorphism varied across loci, with DFR and GS being the overall most variable genes. Variation was lower at nonsynonymous sites than at synonymous sites both within and among species. In the three taxa, within-population variation was generally a small fraction of total variation. Indeed, significant genetic differentiation between locations was detected in four out of eight genes in A. halleri, in seven out of eight in A. l. petraea, and in all three in A. l. lyrata (S_nn, Table 3).

Table 5 shows the variation at synonymous and nonsynonymous sites, as well as the synonymous/nonsynonymous ratio for polymorphism and divergence. Average levels of divergence (number of substitutions per site) between these species and A. thaliana were 0.16 for synonymous and 0.015 for nonsynonymous, approximately the same order of magnitude as estimates reported by WRIGHT et al. 2003 . The nonsynonymous/synonymous ratios were always less than unity, as expected for purifying selection at nonsynonymous positions. Nonsynonymous/synonymous ratios were generally larger for polymorphism than for divergence.

We compared nucleotide variation in six genes (CHI, CHS, F3H, FAH1, GS, and MAM-L) between A. thaliana (KUITTINEN and AGUADE 2000 ; AGUADE 2001 ; this study) and the outbreeding species A. halleri and A. l. petraea. The average total species-wide estimates of variation in A. thaliana for the homologous regions surveyed in the outbreeding species were 0.0035 and 0.0061 for total and silent variation, respectively. These estimates were calculated among ecotypes, which are equivalent to the average estimates among locations presented here for the outcrossing species; hence these estimates are directly comparable across taxa. In A. halleri, the average estimates of nucleotide variation in the same six regions were 0.0048 and 0.0097 for total and silent variability, respectively, while they were 0.0096 and 0.0195 in A. l. petraea. We observe that the average levels of variability are 1.5–3 times lower in the inbreeding species A. thaliana than in the outcrossing species A. halleri and A. l. petraea. Thus, we clearly observe reduced variation in the inbreeder, a reduction that is concordant at least in sign with the difference in mating system. We obtained analytical results that could also explain the observed levels of variation in species with different population structure and mating system (Appendix).

Detecting ancestral variation in A. halleri and A. lyrata:
Table 6 presents the distribution of variable sites in the A. halleri/A. l. petraea and the A. l. lyrata/A. l. petraea pairs, using A. thaliana as the outgroup. In the A. halleri/A. l. petraea comparison, most genes harbored variants in the two new classes of putatively ancestral polymorphisms (S_f1x2, S_f2x1), while only three loci exhibited shared polymorphisms (S_s). The two newly proposed classes of mutations make an important contribution to the total number of ancestral mutations: 38 out of 71 mutations in the A. halleri/A. l. petraea comparison (Table 6A) and 31 out of 36 in the A. l. lyrata/A. l. petraea comparison (Table 6B). The possible difference in population size of these subspecies is also reflected in the number of segregating sites in each new category. Indeed, the total number of mutations fixed in the A. halleri branch that still segregate in A. l. petraea (24) is higher than those fixed in A. l. petraea that still segregate in A. halleri (14). This is consistent with A. halleri having a lower effective population size, since the coalescence time is shorter for A. halleri than for A. l. petraea. The same effect can be detected in the A. l. lyrata/A. l. petraea comparison (28 vs. 3; Table 6B).

The hypergeometric distribution (ROZAS and AGUADE 1993 ; CLARK 1997 ) was used to assess whether the number of putatively ancestral polymorphisms detected in the A. halleri/A. l. petraea comparison (Table 6A) could be the result of recurrent mutation. When all positions were considered, recurrent mutation was insufficient to explain the observed number of S_f1x2 mutations in all cases (Table 7A), while significant or nearly significant departures were also detected when considering only silent positions (Table 7A). Similar results were obtained for S_f2x1 mutations (Table 7A). Finally, recurrent mutation was also insufficient to explain the number of shared polymorphisms in each locus that had mutations of this class (Table 7).

In the A. l. lyrata/A. l. petraea comparison, the presence of ancestral mutations in the classes S_f1x2 and S_s (Table 6B) could not be explained by a recurrent mutation process. Indeed, significant or nearly significant departures from the null hypothesis were detected in most cases (Table 7B). The presence of only a single silent S_f2x1 mutation precluded performing the corresponding test.

Neighbor-joining trees (SAITOU and NEI 1987 ) showed alleles from a given species clustering together except in three loci: F3H, GS, and DFR (not shown). Indeed, the pattern observed at the GS locus is suggestive of recent gene flow from A. halleri into A. l. petraea (Figure S7), with one of the two main haplotype classes present in A. l. petraea being nearly identical (across 1000 bp) to sequences present in A. halleri. Although these two species shared polymorphisms and had no fixed differences between them at the DFR locus, the distribution of polymorphic sites at DFR did not reveal interspecific haplotype sharing (Figure S4). There was no clearly identifiable region in either DFR or any of the other loci with ancestral mutations that we might have attributed to introgression. However, the heterogeneous patterns observed in GS and the other loci do not preclude the possibility of other, undetected introgressed variants.

Table 8 presents the population parameter estimates from the WH isolation model of speciation (WAKELEY and HEY 1997 ). In the A. halleri/A. l. petraea comparison, estimates of population and speciation parameters were similar for most constructed combinations (results not shown). Diversity in the ancestral population (_A) is estimated to be 10⁵-fold higher than that in extant populations, which would imply a 10⁵-fold difference in effective population size between the ancestor and the descendant species (assuming mutation-drift equilibrium and homogeneous mutation rates). Therefore, the estimated time to speciation was rather low compared to estimates from nucleotide divergence (results not shown). The biological implausibility of the parameters estimated from the WH model points to some violation of the neutrality and/or complete isolation assumptions of the model in the A. halleri/A. l. petraea comparison.

The low number of segregating sites in A. l. lyrata probably affects the estimation of population and speciation parameters in the A. l. lyrata/A. l. petraea comparison (Table 8), and these estimates should therefore be interpreted with caution. Although estimates differed between the two constructed combinations, they suggest a lower effective population size for A. l. lyrata.

Neutrality tests:
The implausible parameter estimates produced by the isolation model in the A. halleri/A. l. petraea species comparison may be caused by the violation of assumptions of the model. Natural selection, nonequilibrium demographic processes, or introgression may explain these results. We used neutrality tests for individual loci, as well as several multilocus tests, to test for deviations from a neutral equilibrium model. For individual loci, most tests were nonsignificant, as summarized by the average values of the test statistics over all combinations (Table 9). The only exception was the highly significant Fay and Wu's H statistic for CHS in A. halleri. Indeed, 13 of the 15 polymorphisms in this gene were due to ancestral variants present in a single highly divergent sequence (Figure S2). Also, pairwise MK (McDonald-Kreitman) tests were performed for all loci and species pairs. A marginally significant deviation was detected for CHS in the A. halleri/A. l. petraea comparison (P < 0.049) and for CHS in the A. halleri/A. thaliana comparison (P < 0.047). However, no significant deviation from neutral expectation was detected after Bonferroni correction (results not shown).

Despite the general lack of significance observed, in some tests the deviations from equilibrium neutral patterns were in the same direction (for each particular sample combination and species) in a majority of the genes (Table 9). That was the case for the H statistic, which was negative in seven out of eight genes in A. halleri and in six out of eight in A. l. petraea, and also for the Tajima's D and Fu and Li's D and F in A. l. lyrata. These trends prompted us to perform multilocus tests.

Polymorphism in A. halleri and in A. l. petraea was compared to divergence of each species relative to A. thaliana with the multilocus HKA test. The large number of combinations of constructed samples for all loci (16² x 32² = 2.62 x 10⁵ for A. halleri and 32⁶ x 16 = 1.72 x 10¹⁰ for A. l. petraea) precluded performing all possible tests. After performing multilocus analyses for several combinations, it was clear that test results in both A. halleri and A. l. petraea were consistent for most data subsamples (results not shown). Therefore, only results for one combination per species are shown (Fig 2). In the A. halleri/A. thaliana comparison, the HKA test rejected the null hypothesis of proportionality between polymorphism and divergence (P < 0.002), while the test statistic value was close to significance in the A. l. petraea/A. thaliana comparison (P < 0.073). In both cases, the largest contributions to the overall test statistic were due to DFR, GS, and MAM-L (Fig 2). Significant values can be explained by a larger variance in the polymorphism/divergence ratio than that expected under a neutral equilibrium model. This large variance might be attributable to selection on a subsample of the loci or to violation of assumptions of the model. The multilocus HKA test was not significant in the A. l. lyrata/A. thaliana comparison.

View larger version (15K): In this window In a new window Download PPT slide   Figure 2. Summary of three multilocus HKA tests (HUDSON et al. 1987 ). Each test compares polymorphism within a single species to divergence from a single A. thaliana sequence or from multiple A. thaliana sequences where possible. Open symbols represent contribution to the overall 2 test statistic caused by polymorphism levels at each locus; solid symbols represent contributions caused by divergence. Points above the line indicate deviation in the direction of an excess of polymorphism or an excess of divergence relative to the neutral expectation. Likewise, points below the line indicate deviation in the direction of too little polymorphism or too little divergence relative to expectation. (A) A. halleri; (B) A. l. petraea; (C) A. l. lyrata.

We employed multilocus analyses incorporating across-loci averages for each statistic and nonparametric sign tests (see MATERIALS AND METHODS). As above, test results were similar for most combinations. In A. halleri, most single and multilocus tests gave nonsignificant results. Indeed, only the multilocus m[H] test for A. halleri yielded a significant result (Fig 3A), as well as the nonparametric sign test (eight H values under the median; P = 0.0039). A similar result was obtained in this species when all tests were performed for the combination that yielded the most conservative value for the m[H] statistic: only the tests using the H statistic were significant. The m[H] test was also significant when the CHS gene, which was highly significant by itself, was excluded. In A. l. petraea, only the nonparametric sign test for the H statistic was significant (seven of the eight H values under the median; P = 0.035). The m[H] value was also under the median, indicating a similar, but nonsignificant, deviation (Fig 3A). Equivalent results were obtained for Fay and Wu's H test considering recurrent variation (see MATERIALS AND METHODS), except for the sign test in the species A. l. petraea, where half of possible combinations were nonsignificant. The significant results of the tests for A. halleri and A. l. petraea indicate deviations from the neutral equilibrium model at a multilocus level rather than for single loci. Significant values of the Fay and Wu statistic have been generally explained by positive directional selection, but they can also be the result of recent gene flow (FAY and WU 2000 ). Introgression may be the best explanation for these results.

View larger version (30K): In this window In a new window Download PPT slide   Figure 3. Histograms showing the observed values and the confidence intervals for different test statistics for individual loci and across loci. Shaded boxes indicate the 95% confidence intervals for each locus and for the average of all loci. Solid dots indicate the observed values. A thick solid line indicates the median for each distribution. (A) Fay and Wu's H test statistic for the species A. halleri and A. l. petraea. (B) Tajima's D and Fu and Li's D test statistics for A. l. lyrata.

For A. lyrata ssp. lyrata, both data subsamples gave similar results, although m[D_F] and m[F_F] were significant in only one of the two possible combinations (P = 0.048 and P = 0.035). Fig 3B shows the individual and average multilocus tests for Tajima's D and Fu and Li's D statistics. The same negative pattern was found for all loci, which may indicate population expansion or other demographic processes.

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Relationship between A. halleri and A. lyrata:
Comparative analysis of nucleotide variation at eight loci in A. halleri and A. l. petraea and at three loci in A. l. petraea and A. l. lyrata can shed light on the process of speciation at different timescales. In general, species recently originated from a common ancestor should present few fixed differences, due to the long time required for mutations to accumulate and become fixed. Genetic variants shared by these species may represent polymorphisms present in their common ancestor that have been maintained either by chance or by some form of balancing selection. However, they also may result from introgression between species or from recurrent mutations.

Our comparative analysis of variation in the two subspecies of A. lyrata (A. l. petraea and A. l. lyrata) is based on only three genes, with very low levels of variation in A. l. lyrata. The observed data can generally be explained by a simple isolation model. The parameters estimated using the WH model indicated a significant reduction of variability in A. l. lyrata, suggesting a bottleneck during the colonization of the New World from Eurasia.

In the A. halleri/A. l. petraea comparison, recurrent mutation was insufficient to explain the observed number of putatively ancestral mutations, indicating that most of them are in fact the result of individual mutation events. The high number of mutations observed in the ancestral class could thus be due either to their persistence in isolated taxa since speciation or to an ancient divergence with some degree of subsequent gene flow between species. Results from the WH isolation model were biologically implausible, which suggests violation of some assumptions in this model, quite possibly neutrality or isolation of species. Also, gene trees of some loci showed that alleles from the same species were not clustering together. Specifically, the GS locus shared haplotypes between these two species, which may be explained by introgression.

By combining single and multilocus tests of neutrality we can distinguish between processes affecting individual genes and forces influencing the entire genome. The multilocus m[H] test showed a significant or nearly significant deviation from the simple neutral model in A. halleri and A. l. petraea, respectively, and significant values for the sign test in both species. Although multilocus tests indicate deviation from neutral equilibrium expectations for the complete set of loci analyzed, the negative value of the H statistic for each locus clearly shows an excess of recent mutations at high frequency at every locus. This generalized excess suggests a multilocus effect that could be the result of demographic factors, independent selective sweeps, or introgression. The sign of the test statistic, the sampling strategy used in the tests (see MATERIALS AND METHODS), and the lack of correlation between geographic and genetic distance (results of Mantel tests not shown) exclude population subdivision as the main cause of the detected deviation. It is also highly unlikely that the multilocus deviation was due to positive selection since it would require independent selective sweeps at a majority of the loci studied. Population reduction could also explain a negative tendency of the Fay and Wu H values across the genome (results not shown). If this were the case, this process should be also observed in positive Tajima's D and Fu and Li's D and F values, which is not observed (Table 9). Therefore, we also exclude that hypothesis. Finally, recurrent mutation also cannot explain the observed data, as we included that process in coalescent simulations (see RESULTS). Thus, the most parsimonious explanation of the significant results of the multilocus H test in the closely related species A. halleri and A. l. petraea is reacquisition of ancestral mutations by introgression (FAY and WU 2000 ).

Selection could also play an important role in the maintenance (or loss) of variation added by introgression. Indeed, the pattern of variation detected in A. halleri and A. l. petraea for some loci could be better explained by introgression (see above), despite that haplotype sharing was not generally observed. This observation, together with the detection in HKA tests of the same genes showing similar deviations from neutral expectations in both species (Fig 2), might be explained more easily by selection sweeping (or maintaining) the additional introgressed variation in some loci and not in others than by the stochasticity associated with genetic drift.

A hybridization event is biologically plausible, as interspecific crosses between A. halleri (pollen) and A. l. petraea have resulted in interfertile F₁ individuals (MACNAIR et al. 1999 ). Hybrids formed in the laboratory have greatly reduced vigor relative to nonhybrids of both species (M. R. MACNAIR, personal communication). Thus, although interspecific hybrids have not been reported in nature, rare hybridization events may still occur on an evolutionary timescale.

The existence of gene flow between species alters the pattern of both within- and between-species variation. Indeed, in a simple scenario with two species 1 and 2 that possess mutations only in classes S_x1, S_x2, S_f1, and S_f2 (Fig 1), asymmetrical gene flow from species 1 into species 2 would result in some mutations switching to another class. Some mutations would switch from class S_x1 to either class S_s or S_x2, some from class S_f1 to class S_f1x2, and some from class S_f2 to S_x2. There would be, therefore, a decrease of exclusive variation in species 1 (where fixed and polymorphic mutations would switch to the ancestral mutation class), and an increase in species 2 (by addition of polymorphic mutations). Therefore, one asymmetrical gene flow event between two species with similar effective population sizes (and thus levels of variation) would result in increased variability in the recipient species, a reduction in the number of exclusive polymorphisms in the donor species, and reduction of the estimated time since isolation. On the other hand, reciprocal gene flow would result in more ancestral mutations and fewer fixed mutations in both species, which would increase the apparent level of variation in the ancestral population and reduce the estimated time since isolation. This pattern is consistent with the results of the WH model for the A. halleri/A. l. petraea comparison, supporting, therefore, the hypothesis of introgression occurring after the split of these species. The inclusion of outgroup sequences in our study has dramatically increased the number of ancestral mutations detected (Table 6) and therefore the information regarding the ancestral history of the two species. Speciation models would greatly benefit both from the addition of an outgroup species and from considering migration subsequent to species divergence since this would improve the estimates of isolation time and levels of variation in the ancestral and extant species and may improve the power to detect introgression.

Levels of nucleotide variability in A. halleri and A. lyrata:
Genetic differentiation between geographic locations can be suggestive of (although not equivalent to) population subdivision (Table 3, S_nn tests). In our study, significant genetic differentiation between populations is observed more frequently in A. l. petraea than in A. halleri, indicating that the former species (with a discontinuous distribution extending throughout North and Central Europe) has a higher level of population structure than A. halleri (distributed across Central Europe).

Several factors including migration rate between demes and number of demes might influence a species' effective population size when there is population substructure and might consequently affect the species-wide level of variation (e.g., WHITLOCK and BARTON 1997 ; PANNELL and CHARLESWORTH 1999 ; PANNELL and BARRETT 2001 ; WAKELEY 2001 ; INGVARSSON 2002 ; LAPORTE and CHARLESWORTH 2002 ). Additionally, asymmetrical gene flow between closely related species can cause a differential increase of variation in the two species and, thus, of their effective population size (see above). The average level of genetic variation among populations (total and silent variation in Table 3) is higher in A. l. petraea than in A. halleri. The higher level of variation detected in A. l. petraea than in A. halleri reflects differences in effective population size between the species that could be caused both by differences in their degree of population substructure (as is shown in a simple example in the Appendix, Equation A4) and by differential levels of introgression.

The low level of variability detected in A. l. lyrata, relative to A. l. petraea, might be explained by a bottleneck accompanying the relatively recent colonization event from the Old World (MITCHELL-OLDS 2001 ). The generally negative values of most multilocus test statistics suggest that population expansion followed colonization. Low levels of variation relative to A. l. petraea have also been observed in six loci studied in two populations of A. l. lyrata (WRIGHT et al. 2003 ), in agreement with our data. On the other hand, a survey of variation in the Adh gene detected two main haplotypes in a single population of A. l. lyrata with no indication of population expansion (SAVOLAINEN et al. 2000 ). Although the differing patterns of variation among these studies could be due to differences in sampling strategy, they could also be explained if selection maintained variation in the Adh locus. Multilocus analyses of variation in large samples from multiple populations of A. l. lyrata are necessary to elucidate the demographic history of this subspecies.

Considerations regarding the levels of variation in the outbreeding species A. halleri and A. lyrata relative to the inbreeding species A. thaliana:
Comparison of within-species levels of variation in outbreeding and inbreeding species may reveal the effect of mating system on intraspecific variation. If the assumptions of a Wright-Fisher model hold for the two species being compared (except for mating system in the selfing species), a twofold reduction in the effective population size is expected under complete selfing, and consequently the level of variation should be lower in the inbreeding than in the outbreeding species (POLLACK 1987 ; NORDBORG and DONNELLY 1997 ). Indeed, INGVARSSON 2002 has shown that the ratio of within-population genetic diversity in selfing vs. outbreeding species is much lower in structured populations. In addition, the ratio of the metapopulation levels of genetic diversity in inbreeding vs. outbreeding species can also be heavily influenced by the structure of the populations. Indeed, total levels of variation could be quite similar in selfing and in outcrossing species if migration and effective population size are much lower in inbreeding species than in outcrossing species (INGVARSSON 2002 ).

We also developed a simple analytical expression of the expected ratio of species-wide variation in inbreeding vs. outbreeding species under the same hypothetical population structure in both species (Appendix). In the absence of population extinction (Appendix, Equation A3), the expected ratio is approximated by the expression (¹/₂M + 1)/(M + 1), where M is the population migration parameter 4Nm. A twofold reduction would thus be expected under high migration, with a milder reduction under low migration. On the other hand, the effective rate of recombination is greatly reduced in inbreeding species (GOLDING and STROBECK 1980 ). High levels of inbreeding are thus expected to increase the effects of genetic hitchhiking accompanying selective sweeps (MAYNARD SMITH and HAIGH 1974 ; KAPLAN et al. 1989 ) and background selection (CHARLESWORTH et al. 1993 ) and may greatly influence genetic variability in inbreeding species (NORDBORG et al. 1996 ). A higher than twofold reduction would thus be generally expected in the standing level of variation in the selfing species relative to the outbreeder.

We have not detected, however, the very drastic reduction of variation in A. thaliana that would be expected from differences in mating system, the effect of gene flow in the outbreeders' levels of variation, and the increase in hitchhiking and background selection effects resulting from effectively reduced recombination in the inbreeder. Indeed, A. thaliana has a worldwide distribution and exhibits a metapopulation structure (SHARBEL et al. 2000 ; MITCHELL-OLDS 2001 ), while both A. halleri and A. l. petraea are restricted to Europe and with possibly different degrees of population subdivision. Under neutrality and a reasonable range of parameters (INGVARSSON 2002 ; Appendix), the expected effective population size of the inbreeder could be similar to or higher than that of the outbreeder. Thus, differences in population structure can counteract the effect of mating system on variation, which could obscure the effects of selective processes such as genetic hitchhiking and background selection.

	FOOTNOTES

Sequence data from this article have been deposited with the EMBL/GenBank Data Libraries under accession nos. AJ619855, AJ619856, AJ619857, AJ619858, AJ619859, AJ619860, AJ619861, AJ619862, AJ619863, AJ619864, AJ619865, AJ619866, AJ619867, AJ619868, AJ619869, AJ619870, AJ619871, AJ619872, AJ619873, AJ619874, AJ619875, AJ619876, AJ619877, AJ619878, AJ619879, AJ619880, AJ619881, AJ619882, AJ619883, AJ619884, AJ619885, AJ619886, AJ619887, AJ619888, AJ619889, AJ619890, AJ619891, AJ619892, AJ619893, AJ619894, AJ619895, AJ619896, AJ619897, AJ619898, AJ619899, AJ619900, AJ619901, AJ619902, AJ619903, AJ619904, AJ619905, AJ619906, AJ619907, AJ619908, AJ619909, AJ619910, AJ619911, AJ619912, AJ619913, AJ619914, AJ619915, AJ619916, AJ619917, AJ619918, AJ619919, AJ619920, AJ619921, AJ619922, AJ619923, AJ619924, AJ619925, AJ619926, AJ619927, AJ619928, AJ619929, AJ619930, AJ619931, AJ619932, AJ619933, AJ619934, AJ619935, AJ619936, AJ619937, AJ619938, AJ619939 and AJ582819, AJ582820, AJ582821, AJ582822, AJ582823, AJ582824, AJ582825, AJ582826, AJ582827, AJ582828, AJ582829, AJ582830, AJ582831, AJ582832, AJ582833, AJ582834, AJ582835, AJ582836, AJ582837, AJ582838, AJ582839, AJ582840, AJ582841, AJ582842, AJ582843, AJ582844, AJ582845, AJ582846, AJ582847, AJ582848, AJ582849, AJ582850, AJ582851, AJ582852, AJ582853, AJ582854, AJ582855, AJ582856, AJ582857, AJ582858, AJ582859, AJ582860, AJ582861, AJ582862, AJ582863, AJ582864, AJ582865, AJ582866, AJ582867, AJ582868, AJ582869, AJ582870, AJ582871, AJ582872, AJ582873, AJ582874, AJ582875, AJ582876, AJ582877, AJ582878, AJ582879, AJ582880, AJ582881, AJ582882, AJ582883, AJ582884, AJ582885, AJ582886, AJ582887, AJ582888, AJ582889, AJ582890, AJ582891, AJ582892, AJ582893, AJ582894, AJ582895, AJ582896, AJ582897, AJ582898, AJ582899, AJ582900, AJ582901, AJ582902, AJ582903, AJ582904, AJ582905, AJ582906, AJ582907, AJ582908, AJ582909, AJ582910.
¹ These authors contributed equally to this article.
² Present address: Department of Genetics and Evolution, Max Planck Institute of Chemical Ecology, 07745 Jena, Germany.
³ Present address: Departament de Genètica, Facultat de Biologia, Universitat de Barcelona, 08028 Barcelona, Spain.

	ACKNOWLEDGMENTS

We thank J. Wakeley and W. Stephan for analytical suggestions, M. J. Clauss for helpful comments on the manuscript, H. Kuittinen and P. Saumitou-Laprade for providing several A. halleri and A. l. petraea DNA samples, and J. Kroymann for MAM-L amplification primers and unpublished A. thaliana MAM-L sequences. S.R.-O. and M.A. also thank G. Blasco for her excellent technical assistance and Serveis Científico-Tècnics, Universitat de Barcelona for sequencing facilities. B.S. and T.M.-O. thank D. Schnabelrauch and A. Figuth for outstanding sequencing services and M. Voigt and J. Fritsche for excellent technical assistance. This work was supported by the Bundesministerium für Bildung und Forschung, the Max-Planck Gesellschaft, European Unio