Patterns of Genetic Variation at a Chromosome 4 Locus of Drosophila melanogaster and D. simulans

MATERIALS AND METHODS

Drosophila stocks: We surveyed 38 isofemale lines of D. melanogaster, collected from three different localities. Eighteen lines were collected by M. Noor near Beltsville, Illinois, in 1991; 12 lines by C. Schloetterer near Beltsville, Illinois, in 1994; and 8 lines by M. Kreitman at Terhoon's Farm, Massachusetts, in 1989. For 32 of these lines, we extracted single fourth chromosomes using a y;bw;ci^D/ey^D marker stock. The extracted lines were made homozygous for either yellow or brown, to facilitate detection of contamination. No contamination was ever found in these homozygous lines. For the remaining 6, the extraction lines were lost, so that single flies from the original isofemale stocks were used in the PCR/DHPLC analyses described below. Thirty-three isofemale lines of D. simulans were surveyed from a single locality. These were originally collected by A. Berry, near Tempe, Arizona, in 1990. Single flies from each line were used for PCR/DHPLC. Stocks were maintained at 18°, on standard yeast-sucrose-cornmealagar medium.

Gene region: We analyzed portions of exons 2 and 4, the whole of exon 3, and the second and third introns of the ankyrin gene in both species (Figure 1). Ankyrin was localized to chromosome 4 of D. melanogaster by Dubreuil and Yu (1994), using a cDNA probe. The introns were identified during the course of PCR amplification of D. melanogaster genomic DNA, using a set of primer pairs designed from the ankyrin cDNA sequence (GenBank accession no. L35601). A primer pair (ANK+38, ANK−275; see below) near the 5′ end of the cDNA did not amplify, using a typical protocol designed to amplify 1 kb; we used this pair in a long PCR reaction (Expand Long PCR; Roche Molecular Biochemicals, Indianapolis) and obtained a product of ~6 kb. Upon sequencing (see below), this proved to contain one long (5.4-kb) intron (ankin2), a 108-bp exon (exon3), and a 363-bp intron (ankin3). The complete genomic structure of ankyrin was determined by the Drosophila genome project (gene CG1651; Adamset al. 2000) and confirms this structure, except that our annotation compared to the genome project annotation suggests an additional small intron upstream of our region. A correction has been forwarded to Flybase.

S. Assimacopoulos kindly verified the chromosomal location of the intron by in situ hybridization to polytene chromosomes, using a HindIII restriction fragment from the long PCR product. The same primer pair amplified the homologous region in D. simulans. This product was a 5-kb fragment, with structure similar to D. melanogaster; in particular, there is one long and one short intron, and exon2 is the same length in the two species.

DNA sequencing: To sequence the D. melanogaster region, a shotgun sequencing protocol developed by P. Andolfatto was used. Briefly, the 6-kb D. melanogaster product was mechanically sheared, polished, and blunt-end cloned into a derivative of vector pZero (Invitrogen, San Diego). Approximately 30 positive clones were picked for colony PCR, using universal (M13-20 and M13rev) primers. The resulting products were dyeterminator cycle sequenced using a premixed reaction (PE Biosystems/ABI) and sequenced on an automated sequencer (ABI 377). All portions of the region were sequenced at least twice. The same protocol was used on the D. simulans 5-kb product. This yielded four long contigs. To finish the sequencing, the 5-kb product was cut with a six-base blunt-cutting restriction enzyme (DraI) that cut rarely among the contigs. Subcloning and sequencing the resulting fragments provided the missing D. simulans sequence. Ultimately, we obtained 6071 bp of D. melanogaster sequence and 5057 bp of D. simulans sequence (GenBank accession nos. AY054998 and AY054997, respectively). The sequence for D. melanogaster is from line B45 (Terhoon's Farm population). The sequence for D. simulans was obtained from line s52 genomic DNA for the initial shotgun sequencing, and a cloned long PCR product from line s18 was used to fill contig gaps. The nucleotide coordinates used below are such that +1 indicates the first base of intron 2 in our notation.

PCR/DHPLC: We used DHPLC of DNA fragments (Huberet al. 1993) to identify nucleotide and length differences between homologous fragments of ankin2-3 among lines. We designed primer pairs from the intron sequence to amplify 150- to 400-bp fragments that together span the entire region (Jensen 2000). The primer pairs did not perfectly overlap, but no gap was >105 bp for either species. The pairs covered 5296 usable base pairs (i.e., excluding actual gaps between fragments and the primer sequences) in D. melanogaster and 4459 usable base pairs in D. simulans.

Template DNA for the fragment PCRs was generated as follows. Genomic DNA was extracted from a single fly for each isofemale or chromosome-extracted line. This DNA was used as template for a long PCR reaction, using a high-fidelity DNA polymerase (Expand Hi-Fidelity; Roche Molecular Biochemicals). For D. melanogaster lines, the entire ankin2-3 region was amplified using primers ANK+38 (5′-CGCTTGGTGATGTACGAGTTG-3′) and ANK−275 (5′-TGTCCACATATCCGTCCTTTG-3′). For the D. simulans lines, only the ankin2 intron was amplified, using simulans-specific primers designed from exons 2 and 3 (sankin1U57, 5′-TAATGGAATGGCTTTAGACAACAA-3′; and sankin1L169, 5′-TATGTCCGATATTTCTCCACAGTC-3′), since not every line would amplify well using the melanogaster primers. This long PCR product was used directly as template for the fragment PCRs for all D. melanogaster lines and 23 of the D. simulans lines. The long PCR product for the remaining D. simulans lines was cloned into the TOPO-4 vector (Invitrogen), according to the manufacturer's protocol. For these lines, we used cloned ankin1 from a plasmid prep as template for the fragment PCRs.

Using PCR product as template for secondary PCR reactions raises the issue of “PCR error,” or the occasional incorporation of noncomplementary bases by the thermostable DNA polymerase during PCR amplification, which may result in artifactual variants. This is of particular importance in the study of regions of very low variation like the fourth chromosome. For a high-fidelity polymerase, such as that used to produce ankin1 template here, there is little cause for concern if the PCR product from genomic DNA is used directly as template (for a detailed analysis, see Jensen 2000). But when cloned PCR product is used as template for a line, there is a distinct possibility that at least one base along its length has been changed by misincorporation. Such spurious variants should always appear as singletons and can thus be identified by checking two different clones derived from the same PCR product. In this study, three singleton variants were found among cloned D. simulans lines, of which two were rejected and one accepted.

DHPLC was used to survey fragments for variants as follows. For each fragment or primer pair, all lines were amplified using ordinary Taq polymerase (QIAGEN, Valencia, CA) in a 25-μl reaction, using 1 μl of 1/100 dilution of the template PCR described above (details are given in Jensen 2000). One line was chosen as a comparison standard and was amplified to provide 8 μl standard reaction per line. All fragment PCRs were performed with identical primer concentrations and were primer limited. Thus, final product concentrations were approximately equal for all lines. To prepare samples for DHPLC analysis, 8 μl of each PCR reaction was mixed with 8 μl of product from the standard line, heated to 98° for 3 min to denature, and cooled to 65° over 30 min in a thermal cycler. Five microliters of this denatured and reannealed product was passed over the DHPLC column in a solvent gradient and at a temperature determined from the melting characteristics of the fragment. These parameters were determined by using WAVEMaker software (Transgenomic Inc.). UV-absorbance chromatographs for each line were compared with a chromatograph of the unmixed standard reaction that had been subjected to denaturation/reannealing. Chromatographs reveal variants, because their shapes change with the presence of heteroduplex DNA in the sample. If a line contains a variant with respect to the standard line, the mixed PCR reaction will yield approximately one-half heteroduplex DNA, which under appropriate conditions will give a chromatograph that is different in shape from the unmixed standard. If the standard line and the tested line are identical in sequence, most of the DNA will be homoduplex, and the chromatograph will mimic the standard. As this study shows, the technique can unambiguously reveal single base pair substitutions, single base pair insertion-deletions (indels), and longer indels.

When variants appeared among the lines for a fragment, on the basis of a subjective observation of associated chromatograms, lines were assigned a DHPLC variant type number according to the chromatogram shape. When ambiguities arose in chromatogram interpretation, more types were assigned. Generally, slight differences between chromatograms, such as small changes in retention times of absorbance peaks, did not indicate an underlying mutation. Gross changes in chromatograms between fragments, which resulted in differences in numbers of peaks or marked width differences in single peaks, were reliable indicators of underlying sequence differences. In two instances, variant classes with gross chromatogram differences were nevertheless isosequential. It is likely that nonspecific amplification in the PCR reactions led to these results; we did not follow up these cases. Ultimately, only gross differences were scored as separate DHPLC variants. At least three lines were sequenced, if possible, within each DHPLC variant class, including the standard class, using the original fragment PCR reaction as template in a cycle-sequencing reaction. Unsequenced lines within a DHPLC variant class were assumed to contain the same sequence variants as the sequenced members of that class. Most of those fragments whose chromatographs appeared to be the same as the standard for every line were assumed to be monomorphic and were not analyzed further. This is justified by our preliminary data and the fact that no sequence variation was ever found within DHPLC variant classes in polymorphic fragments. Details are given in Jensen (2000).

It is possible that this survey will not have identified all polymorphisms in the region for these lines. We attempted to minimize the possibility of missing variation by performing DHPLC at multiple column temperatures, when the fragment was predicted to contain several melting regimes (i.e., heterogeneity in GC content); this has been shown to increase the chances of detecting point variation within high-melting-temperature tracts (Transgenomic Inc., personal communication). Representatives of DHPLC classes were sequenced to characterize the underlying sequence changes and to demonstrate the reproducibility of the chromatogram-sequence association within classes. However, it is unlikely that failure to identify a variant would depend on its population frequency, so that the broad conclusions from this study should be little affected by DHPLC inefficiency. On the other hand, there should be no spurious sequence variation introduced by the survey method, since DHPLC variant classes were conservatively assigned and checked by direct sequencing. In fact, this source of error should be reduced relative to direct sequencing alone, since DHPLC and direct sequencing provide independent ways of checking sequence identity.

Evolutionary parameter estimates: Estimates of evolutionary parameters, including θ_w [Watterson's (1975) estimator of the scaled mutation rate θ = 4N_eu], the nucleotide site diversity π estimator of θ (the mean pairwise difference per base pair; Tajima 1983), and Tajima's (1989) D statistic for measuring departure of the site frequency spectrum from neutrality, were calculated using either DNAsp (Rozas and Rozas 1997) or SITES (Hey and Wakeley 1997). Estimates of the scaled recombination rate C = 4N_ec (Hudson 1987) and F_ST (Wright 1951) for the D. melanogaster data were calculated with SITES. In addition, Dr. Jeffrey Wall kindly calculated the scaled rate of gene conversion, G = 4N_eg, on the assumption that all intragenic recombination is caused by gene conversion (Frisseet al. 2001). Here, N_e is the effective population size, u is the mutation rate per nucleotide site, and c and g are the rates of crossing over and gene conversion per nucleotide, respectively, averaged over males and females.

D. melanogaster-D. simulans sequence alignment: Relatively long regions of homology are present between the two species, but the region has been subject to many large insertion/deletion events since divergence between the species. This makes the results of typical alignment algorithms unreliable. However, we wished to avoid as much subjective alignment as possible, to get a reasonably unbiased estimate of divergence. We combined alignment by dotplot and the Needleman-Wunsch algorithm (Weir 1996, Ch. 9), as implemented in the Meg-Align program (DNAStar) as follows. We started with a sliding alignment window of 50 bp (chosen to eliminate signals from microsatellite repeats). A window alignment was rejected if the percentage base identity was significantly low in the following sense. The average divergence for noncoding regions between D. melanogaster and D. simulans is 0.061 (Moriyama and Powell 1996), so that a 50-bp region should contain three divergent sites on average. Assuming a Poisson distribution of divergent sites, the probability of eight or more divergent sites in 50 bp is <0.025.

On the basis of this procedure, the dotplot parameters were set to accept a window alignment with 86% identity or better. This resulted in the dotplot in Figure 2; the figure changes little with changes of a few percent similarity in either direction. Each aligned diagonal was passed to the Needleman-Wunsch algorithm, to assign gaps objectively. Excluding gaps, 2195 bp were aligned by this protocol. Further subalignments were performed in the analysis of the HB element insertion (see results).

Simulations: Since it has been suggested that the fourth chromosome of D. melanogaster has undergone a recent selective sweep (Berryet al. 1991; Hiltonet al. 1994), we ask whether our more extensive data set is also compatible with a selective sweep model. Many scenarios involving both strong and weak selection, which include partial sweeps, multiple sweeps, or sweeps with significant recombination, can be devised to explain this or any other data set (Gillespie 1997, 2000; Fay and Wu 2000; McVean and Charlesworth 2000). But if we consider only a simple model, we can use simulation techniques to explore the relative likelihoods of the data in a specified context and perhaps restrict the number of possible explanations of the data (see discussion).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Dotplot alignment of D. melanogaster and D. simulans ankyrin region. Arrow indicates HB-like transposable element (see text for details).

In these simulations, we consider only “catastrophic sweeps” (Perlitz and Stephan 1997) under the following implicit assumptions: (A) Selection is strong relative to mutation, such that no variant arising during the process of fixation of the selected allele is likely to be sampled; (B) selection is strong relative to recombination, such that no sampled allele is likely to have undergone a recombination event during the fixation of the selected allele. Assumption A is quite weak; assumption B is reasonable for the fourth chromosome, despite the evidence of genetic exchange in the data set, since the absence of crossing over on this chromosome suggests that this exchange results mainly from gene conversion, which causes only low levels of recombination over intervals longer than typical genes (see discussion).

Under these assumptions, a selective sweep completely eliminates variation in a given region. If it is further assumed that evolution proceeds neutrally following the sweep under an infinite-sites model, coalescent techniques could be used to simulate samples drawn from the population. Assume a given time T_s, since the last selective sweep (in units of 2N_e generations), a scaled mutation rate θ = 4N_eu, and a sample of n alleles. Coalescent events can be simulated according to the standard neutral model with constant population size (Hudson 1990), until the accumulated time is >T_s, at which point the value T_s is assigned to subsequent nodes. This effectively truncates the coalescent at the time of the sweep, creating a star-shaped genealogy. Mutations are added to each branch by sampling from a Poisson distribution with mean (θt/2), where t is the length of the branch in units of 2N_e generations.

The simulations can be used as follows to estimate the likelihood of the pair (S, K), the observed number of segregating sites, and average pairwise difference between alleles for a given pair of sweep parameters (θ, T_s). Assume that B samples are generated for each pair of sweep parameters. Write Mδ=Σj=1BIδ(S,K,j), (1) where for the jth iteration, Iδ(S,K,j)={1,if∣Kj−K∣≤δandSj=S0,otherwise.}

Here, δ is a preassigned mesh size for the continuous variable K, and S_j and K_j are the simulated number of segregating sites and mean pairwise difference between alleles for the jth replicate. Following Weiss and von Haeseler (1998), the likelihood of the sample is approximated by L(S,K,∣θ,Ts)=1BMδ(S,K). (2) In this study, B = 10⁵ and δ = 0.1, where δ was chosen to give a fairly smooth likelihood surface. Graphical rendering of the likelihoods was performed with Mathematica 3.0 (Wolfram Research, Champaign, IL).

While it is assumed for the purposes of simulation that recombination is unlikely to have occurred during the substitution of a strongly selected allele, recombination in the genealogy cannot be excluded altogether, on the basis of the evidence contained in the data for both species (see results). Recombination following a sweep is likely to make significance tests based on the joint distribution of S and K under the above model conservative, since recombination is known to reduce the variance of the distributions of S and K (Hudson 1990; Wall 1999).