In this article we describe the structure of a hybrid zone in Argyll, Scotland, between native red deer (Cervus elaphus) and introduced Japanese sika deer (Cervus nippon), on the basis of a genetic analysis using 11 microsatellite markers and mitochondrial DNA. In contrast to the findings of a previous study of the same population, we conclude that the deer fall into two distinct genetic classes, corresponding to either a sika-like or red-like phenotype. Introgression is rare at any one locus, but where the taxa overlap up to 40% of deer carry apparently introgressed alleles. While most putative hybrids are heterozygous at only one locus, there are rare multiple heterozygotes, reflecting significant linkage disequilibrium within both sika- and red-like populations. The rate of backcrossing into the sika population is estimated as H = 0.002 per generation and into red, H = 0.001 per generation. On the basis of historical evidence that red deer entered Kintyre only recently, a diffusion model evaluated by maximum likelihood shows that sika have increased at ∼9.2% yr-1 from low frequency and disperse at a rate of ∼3.7 km yr-1. Introgression into the red-like population is greater in the south, while introgression into sika varies little along the transect. For both sika- and red-like populations, the degree of introgression is 30–40% of that predicted from the rates of current hybridization inferred from linkage disequilibria; however, in neither case is this statistically significant evidence for selection against introgression.
TRANSLOCATION of an exotic species into a novel ecosystem often results in hybridization between the introduced species and related native genera. Where few barriers to gene flow exist, this frequently leads to the rapid introgression of genetic and phenotypic characters from one species into another (Abbott 1992; Rhymer and Simberloff 1996). As well as being of practical concern, hybridization with invading species gives a valuable opportunity to observe the initial stages of hybridization: in contrast, most hybrid zones are ancient (Barton and Hewitt 1985, 1989; Harrison 1993). Here, we describe a quantitative analysis of the genetic interaction between introduced sika deer (Cervus nippon) and native red deer (Cervus elephas) in Scotland and show how a low rate of hybridization can lead to substantial introgression.
Sika were first brought to the British Isles in 1860 (Powerscourt 1884). These animals were subsequently bred and distributed to parks throughout Ireland, England, and Scotland, with new introductions continuing up until the 1930s. Later, sika were either deliberately released or escaped from deer parks, quickly becoming part of the naturalized fauna (Lever 1977). Sika are congeneric with Britain’s native red deer (Cervus elaphus). They are smaller at the shoulder, having roughly half the body weight of red deer, and show a number of other differences in coat pattern, antler complexity, and breeding system (Whitehead 1964, 1972, 1993; Ratcliffe 1991; Staines 1991). Members of the genus Cervus often hybridize and produce fertile offspring (Caughley 1971; Harrington 1973, 1979, 1982; Fennesseyet al. 1991). However, during the period of sika introductions, hybridization in the wild was thought unlikely due to differences in body size. This view persisted despite reports of rapid and complete phenotypic introgression in populations of red and sika in County Wicklow, Ireland (Harrington 1973, 1982) and northwest England (Lowe and Gardiner 1975). Phenotypic hybrids were reported in many areas of Scotland from the 1950s (Whitehead 1964; Delap 1967; Harrington 1973, 1982; Lowe and Gardiner 1975). Concern over the genetic impact of sika on red deer was finally voiced by Ratcliffe (1987) while he was reviewing the status of sika in the United Kingdom and reporting on multiple cases of putative hybrids in Scotland.
In a previous genetic study, Abernethy (1994a,b) mapped allele frequencies at four nuclear loci [two protein loci, superoxide dismutase (sod-1) and 6-phosphogluconate dehydrogenase (6Pgdh); two microsatellite loci, BOVIRBP and OarFCB193] and mitochondrial DNA haplotypes across the sika phenotype range in Argyll. This survey extended ∼130 km from the south of the Kintyre peninsula, past the introduction point of sika at the Carradale estate, to the southern part of the Cowal peninsula [see Figure 1. Note that the distance scale in Abernethy (1994a,b) is inflated]. Abernethy’s survey demonstrated that alleles typical of the opposite taxa had introgressed into both species and found strong linkage disequilibrium and heterozygote deficit in the populations longest exposed to hybridization. Sika alleles at nuclear markers could be detected at high frequency further to the north than could sika mitochondrial haplotypes, suggesting that hybridization at sites far from the introduction point was initiated by dispersing sika stags.
Here, we present a new screening of the original Argyll transect, with 11 microsatellite loci plus a mitochondrial marker. This indicates that the two protein loci used in Abernethy’s study do not have completely diagnostic red and sika alleles, a conclusion that is supported by other population surveys of allozyme variation in Cervus (Gyllenstenet al. 1983; Linnell and Cross 1991; Emerson and Tate 1993). Hence, the extent of hybridization between red and sika in Argyll is lower than was inferred by Abernethy (1994a,b). Robust analysis requires that no assumption is made as to the ancestral state of the introduced and native populations. Our new data show that the deer fall into two distinct classes. These are best analyzed as separate populations that receive an influx of genes as a result of occasional hybridization with each other. If the red- and sika-like populations are considered as one, there will be strong linkage disequilibrium and heterozygote deficit. However, these reflect the relative proportions of red and sika, and do not give a useful description of hybridization and introgression (Abernethy 1994a,b). These problems are general. Most previous quantitative analyses of hybrid zones have either applied to random mating (Barton and Gale 1993) or has concentrated on associations with cytoplasmic markers (Asmussen and Arnold 1991). Moreover, analyses are often based on a “hybrid index,” which counts the number of apparently introgressed loci, on the assumption that the markers being used are fixed in alternative populations (Arnold 1997). Here, we introduce new methods that are appropriate when hybridization is rare and use these to separate the contributions of ancestral polymorphism from current hybridization.
Hybridization introduces sets of foreign alleles; these sets are broken up by recombination and may be eliminated by selection. In the long run, there is some effective rate of influx of introgressed alelles, which ultimately segregate independently and either are fixed or reach frequencies determined by their own selective disadvantage. However, in the early stages of hybridization, hybrids may be sufficiently rare that one can analyze the composition of the two interacting populations separately by treating hybridization as consisting of successive backcrosses into a base population. The base population may contain an appreciable frequency of alleles that accumulated from past hybridization or ancestral polymorphism, but these can be treated as being in linkage equilibrium. Thus, individuals can be treated either as kth generation backcrosses carrying a proportion 2-k of foreign genes or as members of the base population with genotypes in linkage equilibrium. (This view assumes no linkage: information on the distant origin of an individual could come from observing linked markers—cf. Baird 1995; Kaplanet al. 1995.)
Fortunately, one can disentangle the current rate of hybridization from its long-term consequences for allele frequencies. Linkage disequilibrium is generated by the mixing of gene pools with different allele frequencies, but associations among loosely linked markers decay rapidly at a rate equal to the recombination rate. Therefore, associations among loosely linked genes must be due to hybridization in the last few generations. Epistatic selection on rare alleles is ineffective (Turelli and Barton 1994) and so cannot contribute significantly. Thus, observing the pattern of linkage disequilibrium gives an estimate of the rate of current introgression. This argument can be seen in a different way by considering the relative frequency of successive backcrosses. In the absence of selection, these should be in the ratio 1:2:4:..., and each backcross generation should carry a proportion 1§4: 1§8: 1§16:... of introgressed genes. In later generations, these ratios will be distorted by selection, which would also cause heterogeneity in disequilibrium across loci. However, the proportions of recent backcross generations, and hence the current rate of hybridization, can be estimated from the frequency of genotypes with multiple introgressed alleles.
The effect of hybridization is to increase the frequency of introgressed alleles within each population. It is hard to separate the accumulated effects of this recent introgression from ancestral polymorphism. However, some information comes from the spatial pattern of introgression: ancient polymorphism should have diffused to uniform frequency, while recently introgressed alleles will be more frequent near regions where the hybridizing taxa overlap (cf. “demic diffusion” in man; Barbujaniet al. 1994). In our analysis, we first establish that there is significant linkage disequibrium within each population and use this to estimate the rate of current hybridization. We then use simple diffusion models for the spread of sika and red through Argyll to find whether the observed introgression is consistent with the estimated rate of current hybridization and to examine the likely genetic and ecological outcomes of hybridization.
This article describes an extension of the genetic analysis on the sika and red deer samples collected by Abernethy (1994a,b) between 1991–92. A full description of the study site and sampling is presented in Abernethy (1994b). Briefly, samples were collected from an area that covered the range of phenotypically sika-like deer in Argyll, Scotland. This extends from the sika introduction point at Carradale on the Kintyre peninsula to the southern part of the Cowal peninsula (Figure 1). Tissue samples were collected for genetic analysis from 246 deer shot during normal Forestry Commission culling operations. The culling strategy aimed for an overall reduction in deer numbers and was not targeted specifically at red or sika in different areas. Sampling should therefore reflect the relative population frequencies of red-like and sika-like deer in each area. Ten commercial forest blocks were sampled, spaced on a 140-km transect extending along the two peninsulas (Abernethy 1994b). Tissue samples were handled as previously described (Abernethy 1994b), and DNA was extracted using standard procedures (Sambrooket al. 1989).
Primer sets for the amplification of microsatellite loci identified in cattle and sheep (e.g., Bishopet al. 1994; Crawfordet al. 1995) were assayed for amplification products and polymorphism in red and sika deer. Initially 97 primer sets were tested in a panel consisting of four sika deer (two each from Kintyre and Peebles, Scotland), three red deer (one each from Galloway, the Great Glen, and Rum) and one sheep control. Microsatellite variation was assayed by PCR and polyacrylamide gel electrophoresis as described in Slate et al. (1998). MgCl2 concentration was initially kept at 2.0 mm and annealing temperatures to 50–52°. Primer sets that generated microsatellite products in deer (determined by the presence of microsatellite-associated artifacts, e.g., “stutter bands”) were identified, and those in which alleles were not shared by red and sika were selected for further investigation. Nine loci, in addition to the two selected by Abernethy (BOVIRBP and OarFCB193; Abernethy 1994b), were identified as being diagnostic because, in a larger test panel of 44 sika and 44 red of diverse geographic origins, no alleles were shared between the taxa. References for each locus and reaction conditions are listed in Table 1. All putative hybrids (either red-sika heterozygotes or red/sika homozygotes in a background of the opposing taxa) identified in the first round of screening were genotyped a second time to confirm hybrid status.
Linkage information was available for all of the microsatellite markers from either red deer (Tateet al. 1995; M. Tate, J. Slate, P. Fennessy, R. Anderson, H. Matias, M. McEwan and K. Dodds, unpublished data), ovine (Crawfordet al. 1995), or bovine (Bishopet al. 1994) linkage maps. All markers map to separate chromosomes except for BM6438 and RM95 and for TGLA387 and DRB3. In red deer BM6438 and RM95 both map to chromosome 31, but are located at opposite ends of the chromosome and show a high recombination frequency. TGLA387 and DRB3 map to chromosome 20 in sheep, but there is no evidence to suggest that they are linked in red deer or cattle (M. Tate, J. Slate, P. Fennessy, R. Anderson, H. Matias, M. McEwan and K. Dodds, unpublished data; Bishopet al. 1994).
Mitochondrial DNA haplotypes that were diagnostic between sika and red were assayed in two ways. Abernethy (1994b) originally used PCR amplification of 16s and ND1 fragments, followed by RFLP analysis. Further mitochondrial genotypes were assigned in the current study using a diagnostic 39-bp tandem repeat in the left domain of the mitochondrial control region (Cook 1993; Nagata 1995; Cooket al. 1999; Nagataet al. 1998). Red deer have a single repeat, while sika deer have multiple repeats. Diagnostic length variation in the mitochondrial control region was assayed by PCR amplification followed by agarose gel electrophoresis. PCR was performed in a reaction volume of 25 μl using 150 ng genomic DNA; 1× PARR Excellence PCR buffer (Cambio Ltd., Cambridge, UK); 2.0 mm MgCl2; 1 mm each of dCTP, dGTP, dATP, dTP (Pharmacia, Piscataway, NJ); 5 pmol of each primer H16498 and L3 (Kocheret al. 1989; Cook 1993: see Table 1); and 1 unit of Taq DNA polymerase (Applied Biosystems, Foster City, CA). Reactions were overlaid with ∼20 μl of mineral oil, and thermal cycling was carried out for 40 cycles of 95°, 1 min; 50°, 1 min; and 72°, 1 min 30 sec. Following PCR, 5 μl of loading buffer (0.25% bromophenol blue, 40% sucrose w/v in dH2O) was added directly to the reactions. PCR products were resolved by electrophoresis on 3% Nuseive GTG agarose gels (Flowgen), with 1× tris-acetate electrophoresis buffer (Sambrooket al. 1989), at ∼120 V for 2 to 4 hr. Products were visualized by ethidium bromide staining with UV illumination. Argyll red deer were fixed for a single repeat and produced a band size of 350 bp, while Argyll sika were fixed for three repeats and yielded a fragment of 430 bp. During population screening fragments were scored in relation to two marker individuals of known genotype.
ANALYSES AND RESULTS
These analyses describe the structure of the hybrid zone in terms of the contribution to red and sika populations from current hybridization and from past introgression or ancestral polymorphism. We first present allele frequencies, using maximum likelihood to allow for the presence of null alleles. We then describe an analysis based on genotype frequencies, which introduces tests of significance for linkage disequilibrium when hybridization is rare, plus maximum-likelihood estimates (MLEs) of the rate of hybridization and introgressed allele frequency in red and sika populations. Finally we describe a diffusion model that specifically assesses the influence of ecological parameters on the rates of spread and increase of deer populations and spatial variation in the rates of hybridization and introgression in red and sika. All the analyses described in this article were performed using a Mathematica 3.0 notebook (Wolfram 1996). A copy of this notebook is available on request from the authors or over the worldwide web from http://helios.bto.ed.ac.uk/evolgen/.
Estimating allele frequencies: Each of the 246 individuals could clearly be classified as red-like or sika-like, on the basis of its multilocus genotype. This division of individuals corresponded exactly to assignments based on phenotype at the time of culling. However, phenotype cannot be used as an indicator that individuals carry introgressed alleles. Table 2 shows the frequency of each allele within the 79 sika-like individuals and within the 167 red-like individuals. There is much more diversity within the red-like population than within the sika; this may be due to the limited set of alleles in the initial sika introduction (nine females and two males were introduced at Carradale in 1893) and to the diverse origins of Scottish red deer.
An allele was classified as typical of red or sika according to whether it was more frequent within red or within sika; most of the analysis was based on this pooling of alleles. This allocation was made to simplify calculations: we did not assume that the ancestral sika population contained only alleles typical of sika or the ancestral red population contained only alleles typical of red. In almost all cases, classification was straightforward. However, the origin of some alleles that are rare in both red and sika populations is unclear. For example, allele 146 at BOVIRBP has a frequency of 0.055 in sika-like individuals and 0.063 in red-like individuals. On the one hand, it is more likely that a rare allele would be found within the diverse red population than within sika. On the other, allele 146 is closest in length to allele 144, which is found in 81% of sika-like individuals; if microsatellites mutate to new alleles of similar length, this suggests an origin from the sika introduction. We therefore assign this allele as sika. Such uncertainties are rare and do not affect the assignments of those crucial individuals carrying introgressed alleles at multiple loci, which are inferred to be of recent hybrid origin (see below).
There are 246 individuals in the sample, each scored for 11 unlinked autosomal microsatellite markers plus mtDNA. All individuals are clearly sika-like or red-like; there are 54 apparently “pure” sika, 28 sika-like hybrids, 144 apparently pure red, and 25 red-like hybrids (Table 3). The mtDNA shows a similar pattern to the nuclear loci, with no indication that it introgresses to a different extent into either sika or red. Most hybrids are either heterozygous at a single nuclear locus or have discordant mitochondrial DNA. However, there are a substantial number of individuals that are introgressed at several loci. Furthermore, many individuals (especially within the red-like population) are apparently homozygous for rare introgressing alleles. Apparent red homozygotes are found within sika-like hybrids once at BOV-IRBP and once at HH064; within red-like hybrids, sika homozygotes are found once at TGLA387 and FCB048 and four times at DRB003. Homozygosity is most unlikely to be due to crossing between red and sika-like individuals, because this would give extensive heterozygosity; such intercrossing must be rare, otherwise F1-like hybrids would be seen. Homozygosity cannot be explained by matings between backcrossed hybrids, because this would only lead to homozygosity if the parents passed on introgressed alleles at the same loci. We therefore believe that the most likely explanation is that such apparent homozygotes are in fact heterozygous for null (nonamplifying or unresolveable) alleles (Callenet al. 1993; Paetkau and Strobeck 1995; Pembertonet al. 1995).
If there are null alleles at high frequency, the null homozygotes should appear as missing data. To examine this possibility, we made MLEs, for each locus in each taxon separately, of the frequency of nulls and of introgressed alleles, on the assumption that missing data are indeed null homozygotes (Table 4). (We ignore variation in introgression across sites.) Within the red population, there is a good fit if one assumes that all the missing data are null homozygotes: for none of the four loci is there a significant deviation from this assumption. Null frequencies must be high: ∼40% for DRB003, 23% for TGLA387, and 20% for FCB048. Notably, it is these loci that show the highest frequency of missing data: loci that do not show apparent homozygotes for rare alleles never show more than six missing genotypes and average 3.4 missing values out of 246 individuals scored. In sika, there are fewer missing data, and only two homozygotes for rare alleles (both in the same individual; Table 4).
These estimates are based on pooling alleles into red and sika classes (see Chakrabortyet al. 1992; Brookfield 1996). A useful check on the presence of nulls is to use the full multiallelic dataset to look for an excess of apparent homozygotes among polymorphic alleles. The MLE was found using the EM algorithm (Yasuda and Kimura 1968; Dempsteret al. 1977; see Weir (1996, pp. 67–73) for a discussion of the allele frequencies in the MNS blood group system. For further details, see the on-line Appendix to this article at http://helios.bto.ed.ac.uk/evolgen/). Null frequencies made in this way are shown in Table 2. The results are entirely consistent with the analysis based on pooled classes of alleles, strongly supporting our interpretation that homozygotes for rare alleles are in fact heterozygotes for nulls. For the remainder of this article, we simply reclassify apparent homozygotes for rare alleles as heterozygotes. The implicit assumption is that nulls derive from the red population, rather than having themselves introgressed. (Note that the high-null-allele frequency estimate in sika at BOVIRBP is based on a single null hymozygote.) This is reasonable, because at the three loci for which nulls are at high frequency in reds, no apparent homozygotes are seen within the sika-like population. Thus, neglecting introgression of nulls will cause negligible error.
Figure 2 shows the cline in frequency of sika-like individuals, together with the proportions of putative hybrids within the sika-like and red-like populations. Clines for individual loci are all similar and there is no indication that allele frequency patterns differ between loci along the transect. Although introgression is rare at any particular locus (at most 6.1% at 0 km in reds, averaged over loci), individuals carrying apparently introgressed alleles are common where the two taxa overlap (maximum 40% at 50 km). These individuals may indeed be hybrids derived from recent crosses between sika and red. Alternatively, they may be carrying alleles that were polymorphic at low frequency in the ancestral population or that entered by hybridization in the distant past. Examination of genotype frequencies allows us to distinguish these possibilities.
Genotype frequencies: Hybridization introduces sets of alleles that gradually disperse through sucessive backcrosses by segregation and recombination. Alleles that entered either population many generations back will by now have reached linkage equilibrium. We can therefore estimate the rate of hybridization over the past few generations from the linkage disequilibrium or, in other words, by estimating the excess of individuals carrying multiple introgressed alleles. One approach would be to classify each individual as being a first-, second-, or later-generation backcross, according to whether it carries 1§4, 1§8,..., of its genome introgressed (Nason and Ellstrand 1993; Boecklen and Howard 1997). However, this would be accurate only with an extremely large number of loci. Instead, we first establish that there is a significant excess of “complex” hybrids (Table 3) and then use likelihood to estimate the rate of introgression and the degree of ancestral polymorphism.
Under linkage equilibrium, the number of introgressed alleles is approximately a Poisson variable, provided these are rare; this gives a simple test for an excess of complex hybrids. First, consider just the 11 microsatellite markers. Table 5 compares observed and expected degrees of introgression into sika for the five informative sites and gives the G-statistic along with the P value. Overall, G14 = 20.51 (P = 0.0581). Table 5 also compares observed and expected degrees of introgression into red for the nine informative sites in the same way; overall, G15 = 21.31 (P = 0.213). (Note that the number of degrees of freedom here is equal to the maximum number of observed heterozygotes; this is because one is comparing classes up to this value and the higher classes, with observed numbers zero.) Because the expected numbers are so small, testing the significance of G against the asymptotic χ2-distribution is inaccurate. Table 5 also shows a randomization test, in which diploid genotypes were shuffled across individuals 1000 times. This test allows for variation in allele frequency across loci, because the null hypothesis is that the frequencies at each locus are as observed, but randomized across individuals. This test shows significant linkage disequilibrium within sika in site 3 and within red in site 2. There is also a marginally significant linkage disequilibrium within sika in site 4. The linkage disequilibrium rests, however, on only three individuals (one in site 3, sika, with five heterozygous loci; one in site 4, sika, with five introgressed alleles; and one in site 2, red, with five heterozygous loci). Examination of mitochondrial haplotype shows no indication of a strong association with nuclear markers; however, there are too few cases of introgression to test for heterogeneity of associations involving this or any other marker.
The small but significant linkage disequilibrium within red- and sika-like populations indicates a low level of introgression in the few generations preceding the sample. To estimate the rate of introgression, we must calculate the probability of sampling the observed genotypes, given that rate. In general, it is hard to describe the composition of a hybrid population: there is a total of 22n diploid genotypes, of which only 3n can be distinguished. All linkage disequilibria up to nth order must be accounted for, and in addition, there are associations between the parental genomes up to order 2n (where n is the number of loci). With a low level of hybridization the problem is simpler: most of the population is of one or other type. In the Argyll data, matters are even simpler, because we have sampled no F1 genotypes and because the markers are effectively unlinked.
To proceed, we make the crucial approximation that there is no linkage among markers or with genes that determine mate preference or fitness. We also assume symmetry across the sexes, so that the mitochondrial marker can be treated simply as another unlinked marker (albeit present in one rather than two copies). We consider all loci to be equivalent in their immediate ancestry; however, there may be different allele frequencies in the hybridizing populations due to different long-term introgression or ancestral polymorphism, leading to some variation across loci.
The neutral expectation that backcross classes are in the ratio 1:2:4..., implies an infinite accumulation of introgressed alleles: the frequency of introgressed alleles in each class decreases as 1:½:¼..., and so the total frequency in each class is the same (as it should be if introgressed alleles persist indefinitely). Backcrosses may be eliminated by selection, but even then, some foreign alleles will be incorporated in every generation (see Barton and Bengtsson 1986). Thus, introgression must be limited either by the time since hybridization or by selection on the markers themselves. We can ignore these deviations from the null model, because in practice we can only detect the effects of recent hybridization on genotype frequencies. (Ratios between backcross generations that differ by a factor of two give almost identical fits, showing that with these data, we have no power to detect selection through the rate of elimination of backcrosses.) We assume a base population in linkage equilibrium with specified allele frequencies (which may vary across loci) and superimpose backcrosses in the proportions 2-t. This series is truncated at T = 4 on the grounds that subsequent backcross generations would be undetectable with 12 loci.
The probability of observing a particular genotype can now be derived, allowing MLEs of the rate of hybridization, H (the proportion of genes that enter the population through backcrossing per generation), and of the frequencies of introgressed alleles in the base population, ui (1 ≤ i ≤ n, the number of loci). The donor population has frequency ui + Δui; we can assume to a close approximation that introgressing alleles are fixed in the other population, so that for the base population Δui = 1 - ui. The frequency of backcross classes is Hαt, with α = 2 and t = 1 labeling the F1. Counting up to t = T = 4, the total frequency of backcrosses is H((αT+1 - 1)/(α-1)) < 1. The probability of observing a haploid genotype X is then (1) where Xi = 0 or 1. Equation 1 also applies to the mitochondrial haplotype, provided that introgression is equally likely to be down the maternal or paternal line. (It would be straightforward to extend Equation 1 to models where the frequency of reciprocal crosses differed.) Note that with rare hybridization the pairwise linkage disequilibrium is equal to twice the rate of hybridization, H.
There is no indication that patterns differ across loci, and so we assume that the frequency of rare alleles in the base population is the same across all loci: ui = u. (We refer to u as the frequency of introgression, on the understanding that it might also represent ancestral polymorphism.) For sika, with parameters constant across populations, there is a significant improvement in likelihood by allowing both hybridization (H = 0.0020) and introgression (u = 0.0141). This increases log(L) by 7.55 relative to introgression alone (u = 0.0226) and by 7.96 relative to hybridization alone (H = 0.0089). Allowing both parameters to vary across populations gives log(L) = -159.46, which is an improvement of only 3.32 over assuming constant values, at a cost of 8 d.f. Results for the red-like population are similar. With parameters constant across populations, there is a significant improvement in likelihood by allowing both hybridization (H = 0.0010) and introgression (u = 0.0050). This increases log(L) by 9.20 relative to introgression alone (u = 0.0086), and by 6.67 relative to hybridization alone (H = 0.0036). Allowing both parameters to vary across populations increases log likelihood by 11.21, at a cost of 18 d.f. Thus, for both red and sika, hybridization (H) and introgression (u) are both significant; in neither case is there a suggestion of significant variation across sites. However, there are too few data to fit parameters separately for each site. Below, we fit a specific model for heterogeneity in introgression and hybridization.
Diffusion models: The frequency of hybrids derived from recent backcrosses (or equivalently, the strength of linkage disequilibrium) gives an estimate for the current influx due to hybridization. This gives a prediction for the frequency of introgressed alleles within the red- and sika-like population and its spatial pattern, which can be compared with the observed values, ui. Such a prediction requires that there be no selection against introgression and also requires assumptions as to the past distributions and abundance of red and sika, the generation time, and the way hybridization depends on red and sika densities. Here, we present a simple diffusion model for the spread of sika and red deer and consider whether the observed degree of introgression is compatible with what is known of the history of the contact and with a null model of neutral introgression. (The use of diffusion models to describe biological invasions is reviewed by Shigesada and Kawasaki 1997.)
First, consider the spread of sika, disregarding genetic variation among sika-like individuals. The wave is centered at 50 km north of the introduction point at Carradale and is ∼60 km wide; the frequency of sika never reaches 100% even at the introduction point (Figure 2). The simplest model would be to assume that the density of the whole deer population is regulated to a fixed and uniform value; that sika have a constant advantage, r, over red (expressed in terms of the difference in the rate of increase between red and sika populations); and that sika were introduced at some low density. Then, the proportion of sika, p, is governed by Fisher’s (1937) equation for the spread of an advantageous allele: (2)
The dispersal rate, σ2, is defined as the variance of distance between parent and offspring divided by the generation time. Distance, x, is measured as the distance north from the introduction point and runs from the southern end of the Kintyre peninsula at x0 = -30 km (Figure 1). This model is described by the time since release, T; the initial proportion of sika, ; the rate of increase, r; and the dispersal rate, σ2. The time since release we take to be T ∼ 77 yr before the samples were collected (around 15 sika escaped from the Carradale estate in 1914; sampling was in 1991–1992; see Introduction), and the pattern is insensitive to the initial numbers, J0. Hence, the rate of increase and the dispersal rate could be estimated from the present position and width of the cline (Figure 2). It would then be possible to predict the opportunity for hybridization from this diffusion model.
This simple model is inadequate, because red deer were absent from the Kintyre peninsula until the 1960s, when populations became established in the forestry plantations at its northern end (Whitehead 1964; Ratcliffe 1987). Therefore we must model the separate spread of red and sika. However, it is not clear how the two taxa interact. In Argyll, red and sika use the same habitat in slightly different ways; for example, sika make greater use of denser cover (Chadwicket al. 1996) and show some dietary differences (Abernethy 1994a). However, despite this partitioning, they interact through hybridization and ecological competition. A general model for the numbers of red and sika (nr, ns) is (3) Here the α-parameters are competition coefficients, defined for each class of individual (see Shigesada and Kawasaki 1997, Chap. 6). If only red deer are present, an equilibrium carrying capacity nr = 1/αrr is approached; colonization of new habitat advances at a speed σr (2rr1/2). (This is a minimal speed, which is approached asymptotically; the advance may initially be faster and occur over a broader wavefront; Fisher 1937.) The converse applies for sika alone. In another possible scenario, where interspecies competition between red and sika is weak relative to intraspecies competition (αrs ⪡ αrr, αsr ⪡ αss), then the two waves might spread past each other, with only a slight slowing and a slight reduction in numbers due to competitive interactions. At the other extreme, with strong competition sika might be able to increase from low density within the red population (αsr <αrr), in which case they would advance to displace the reds. For more symmetric levels of competition, sika might be unable to increase from low density and vice versa; nevertheless, the boundary between sika and red might advance in favor of sika.
Robust estimates of deer density in Argyll are only available for the last 3–4 yr (Forestry Commission; G. Swanson, unpublished data), and information from before this period is largely incomplete or anecdotal. Equation 3 depends on too many parameters to be fitted to these scant data. However, we are primarily concerned here with the opportunity for hybridization between sika and red and its consequences for the pattern of introgression. The contact between sika and red is constrained by what is known of the past history and so may be insensitive to the parameters in Equation 3. An explicit model of the population dynamics is useful in that it shows what information would be needed to account more fully for the spread of these two taxa and to predict their future course.
We assume that red and sika disperse at the same rate (σs = σr = σ), and that they increase at the same rate from low density (rs = rr = r). To begin with, we assume that the fitness of each type is determined by total deer density (i.e., αrr = αrs, αss = αsr); we later examine the opposite extreme, where the two types do not interact, and so spread independently. We assume that the equilibrium density of sika is higher than that of red (1/αss < 1/αrr; Chadwicket al. 1996). Sika were introduced at some low proportion J0, 77 yr before sampling. To make a rough estimate of J0, we take the Kintyre peninsula to be 10 km wide on average and assume a carrying capacity for sika of 6 km-2. The latter figure is derived from the number of sika culled in forest, taken to be 20% of the population (based on Forestry Commission records; G. Swanson and D. Hendry, unpublished data); the highest value is taken as approaching the carrying capacity and is halved to allow for roughly equal proportions of forest and open habitat. Then, a release of 15 sika corresponds to J0 = 15/(10 km) (6 km-2) = 0.25 km-2. Finally, we must estimate the time at which red deer began to advance southward into the Kintyre Peninsula. Red deer were first established south of Tarbert (30 km north of the release point) in 1964 (Ratcliffe 1987); we take this vanguard to correspond to the leading edge of the wave of advance, at 10% of carrying capacity.
With these assumptions, we have free parameters, which determine the rate of advance of sika in competition with red: the dispersal rate, σ; the initial rate of increase, r, and the relative carrying capacities for red; and sika, β = αrr/αss. We fit these parameters by maximum likelihood, based on the numbers of sika- and red-like individuals sampled across the transect (Figure 2). Equation 3 is solved numerically using a simple stepping-stone model, with a proportion m/2 of individuals exchanged in each direction; the number of demes is determined by setting m as close to 1§2 as possible. Assuming that sika have a 50% higher equilibrium density than red (β = 1.5), the best estimate is that r = 9.2% yr-1 and σ = 3.7 km yr-1/2; log[L] = -10.49. Estimates are almost identical if sika are given a different competitive advantage (β = 2, r = 8.0% yr-1, σ = 3.8 km yr-1/2, log[L] = -10.33; β = 1, r = 12% yr-1, σ = 3.9 km yr-1/2; log[L] =-10.97). These values compare well to estimates from other ecological studies (Whitehead 1964; Ratcliffe 1987; Abernethy 1994a; Chadwicket al. 1996). Sampling errors on these estimates are fairly small: for example, with β = 1.5, the 2-unit support limits for r are 8.3–10.4% and with σ, 3.21–4.64. (These limits are asymptotically equivalent to 95% confidence intervals; see Edwards 1972.) However, there is significant unexplained variation around the fitted model (log[L] = -10.49 for β = 1.5, with ∼8 d.f.); in particular, there are rather more red deer around the release point than predicted (Figure 3). This might be due to long-range dispersal or to the presence of red deer in the south before the dates indicated by Ratcliffe (1987). These statistical limits, which only reflect sampling error, are misleading: the main uncertainties stem from the many assumptions that we have made to arrive at a simple model. We have examined many variations on our model and find that these make little difference to the estimates for the initial rate of increase of sika and the dispersal rate. The estimates are insensitive to the competitive interactions between sika and red, simply because under our assumptions, the two have only recently met in Kintyre.
We can now ask whether this model for the spread of sika and red is consistent with the observed spatial pattern of genotype frequencies (Figure 5). This requires that we make some assumption about how the rate of hybridization varies with the relative abundance of sika and red. It seems simplest to suppose that the rate at which hybridization introduces foreign genes is proportional to the frequency of the opposite type; we assume that the same relationship holds both for the rate of current hybridization (denoted H) and the rate of past introgression (denoted I). Thus, if we consider introgression into the sika-like population, the rates of current hybridization and past introgression are H = H*q, I = I*q, where q is the proportion of red-like deer in the whole population and H*, I* are the ratios of H and I, respectively, to q. We express H and I in terms of H* and I* to account for the fact that in the diffusion model these parameters may vary along the transect in relation to spatial and temporal variation in the proportion of red deer. Note that H and I have different units. H is the proportion of genes that enter the population through backcrossing per generation and is measured from linkage disequilibrium; in the absence of selection, the proportion of F1’s is 2H. I is the long-term rate of increase in allele frequency per year due to hybridization, and determines the observed frequency of introgressed alleles.
First, consider introgression into sika. Let the frequency of red alleles within sika be us. The diffusion model of Equation 3 extends to (4) The first term in Equation 4 represents the diffusion of alleles within the sika-like population; the second, directional gene flow from regions where sika are more abundant (Nagylaki 1975); and the third, the influx of genes from the red-like population, at a rate assumed to be proportional to their abundance [I = I*q, where q = nr/(nr + ns)]. Initially, red and sika are assumed to be fixed for alternative alleles (ur = 1, us = 0 at t = 0).
The three coupled equations for the spread of sika (Equation 3) and for introgression from red into sika (Equation 4) and the converse equation for introgression into red are solved numerically using a stepping-stone model, in the same way as Equation 3 alone. The likelihood of the model is then calculated as the product of three components: first, the probability of the observed proportion of sika, p; second, the probability of sampling the observed sika-like genotypes, assuming a pool of alleles at linkage equilibrium at frequencies given by Equation 4, plus up to four generations of backcross hybrids generated at a rate H = H*q; and third, the probability of sampling the observed red-like genotypes. In principle, all the parameters should be estimated together, using the full likelihood. However, to a good approximation we can assume that foreign alleles are fixed in the donor population; then, parameters for introgression into the red-like population can be estimated independently of those for introgression in the opposite direction.
Under this estimation procedure, the population is divided arbitrarily into two parts: the pool of alleles at linkage equilibrium, which are assumed to have been generated by influx at a rate I; and four generations of backcrosses, which each contribute a component H to the allele frequency. In the absence of selection against introgression, Iτ = H, where τ is the generation time. Selection against introgression would be reflected in a reduction in the estimated value of I, below H/τ. [The factor by which selection reduces the net influx of neutral genes was termed the “gene flow factor” by Barton and Bengtsson (1986)].
In estimating introgression and hybridization, we tentatively adopt the model in which both sika and red increase at the same rate (r = 9.2% yr-1) and have the same dispersal rate (σ = 3.7 km yr-1/2), but sika have 50% higher carrying capacity (β = 1.5). The best estimate is that current hybridization is at a rate H* = 0.0121 gen-1 (2-unit support limits 0.0033–0.0279) and past introgression I* = 0.0012 yr-1 (0.0004–0.0019); this gives log(L) =-168.50 (Figure 4a). The estimated rate of past introgression is much smaller than the rate of current hybridization: taking the generation time as τ = 4 yr (based on Forestry Commission records; G. Swanson and D. Hendry, unpublished data), the discrepancy is I*τ/H* = 40%. However, this is not significantly different from 1 (Figure 4a), and so there is no more than a suggestion of selection against red alleles within the sika genetic background. This spatial model predicts that the frequency of introgressed alleles should decrease sharply away from the region of overlap (Figure 5a). In fact, the degree of introgression is similar throughout the sika-like population (dots in Figure 5a). Correspondingly, the spatially uniform model gave a better fit [log(L) = -162.78 with u = 0.0141, H = 0.0020]. This discrepancy could have several explanations. The “introgressed” alleles might in fact have been present at low frequency within the introduced sika (as a result of some previous hybridization with red, by is possible, but is rendered less likely by the low genetic diversity within sika (Table 2). Red might backcross into sika even when rare, so that the influx into sika would occur further to the south. Finally, some red deer may have been present in the south of Kintyre when the sika were introduced: hybridization in the early stages would lead to a spatially uniform pattern. Under all these explanations, selection against introgression would be required, because otherwise introgressed alleles would be seen at higher frequency.
Now, consider introgression into the red-like population. In this case, the frequency of alleles typical of sika does increase sharply to the south, where the two types meet (dots in Figure 5b). The best estimate is that there is no current hybridization (H* = 0) and a rate of past hybridization I* = 0.0025 yr-1; this gives log(L) = -169.84. However, this model is much less likely than a spatially uniform model; a frequency of introgressed alleles of u = 0.0050 and hybridization H = 0.0010 gen-1 gave log(L) = -151.73. The fit is poor because there is an appreciable frequency of alleles typical of sika in the northernmost samples. If we suppose that these alleles were present at some low frequency, u, within the ancestral population of red deer, the fit is much better: u = 0.0049, H* = 0.0051 (0.0009–0.0149), I* = 0.0004 (0–0.0027) gives log(L) =-152.38 (Figures 4b and 5b). While the high diversity within the red-like population (Table 2) makes ancestral polymorphism plausible, other explanations are possible. Sika might be especially likely to hybridize with red when rare, so that introgression would occur in the north of the transect. Suppose that there is no ancestral polymorphism, but rather, that the rate of backcrossing into the red-like population is almost constant when sika are above some small threshold frequency [I = I*p/(0.001 + p), H = H*p/(0.001 + p)]. The most likely estimate is then that H* = 0.0032 (0.0007–0.0072), I* = 0.0002 (0.0001–0.0004), with log(L) = -153.86. Adding the possibility of ancestral polymorphism to this model of frequency dependence gives no significant improvement [log(L) = -153.03]. With either ancestral polymorphism, or frequency-dependent hybridization, the rate of current hybridization is estimated to be much greater than the rate of past introgression. Taking the generation time to be τ = 4 yr, and assuming ancestral polymorphism, the discrepancy is I*τ/H* = 31%, suggesting some selection against introgression. However, as for introgression into sika, this is not statistically significant (Figure 4b).
Our analysis of 11 unlinked microsatellite loci and of the mitochondrial DNA shows that all 246 deer sampled in the Kintyre transect fall into two classes, which correspond to their sika- or red-like phenotype. However, many individuals carry alleles typical of the opposite population at one or more loci. This pattern might be due either to hybridization between red and sika since their recent contact or to polymorphism within the ancestral populations. Three arguments suggest that the former process best accounts for most introgressed alleles. First, the distribution of allele sizes in sika consists of one predominant allele plus several rare alleles, each of which is common in the red-like population (Table 2). This argument does not apply for hybridization in the opposite direction, where the diversity of red alleles obscures introgression from sika; also, it does not show when hybridization took place. Second, the frequency of alleles typical of sika within the red-like population increases markedly to the south, where sika are common (Figure 5b). The converse pattern is not seen within sika, though this is to be expected if introgression took place while sika were increasing from low frequency after their introduction. Finally, three individuals carry 5 apparently introgressed alleles, typical of the opposite race, out of 23 alleles sampled per individual; these are probably first- or second-generation backcrosses. There is significant linkage disequilibrium within both red- and sika-like populations, which can be accounted for by a low rate of current hybridization (H = 0.002 per generation into sika, and 0.001 per generation into red).
At first sight, the lack of F1-like genotypes in our sample is paradoxical. However, (disregarding selection and nonrandom mating) one expects twice as many first-generation backcrosses, four times as many second-generation backcrosses, and so on. Taking the average rate of hybridization as H = 0.0015, we expect to see a proportion 2H of F1’s, or 0.74 individuals in the whole sample. Thus, it is not unlikely that F1’s would be missed. Red-sika hybrids in captive herds do not show any reproductive difficulties (Harrington 1979), and our data show that F1’s in the wild must also breed successfully, because otherwise, more would be required to account for the presence of later generations of hybrids. Variation in the breeding success of recent hybrids in red- and sika-like populations could account for the higher level of introgression into sika. Body size is an important factor in dominance and reproductive success for both sexes in the cervidae (Clutton-Brocket al. 1982). Recent hybrids are intermediate in body size between red and sika (Harrington 1979; Bartos 1992) and so may show increased breeding success in sika populations compared with red populations. In addition, smaller body size in sika-like females may facilitate a reduced birth interval and therefore an enhanced reproductive rate compared to red-like deer (Chadwicket al. 1996). This could be one aspect of phenotype that is under selection; there may be a threshold body size for females, related to local environment, which determines the cut-off between different maximum reproductive rates for females. Future work combining genetic and phenotypic data should allow such questions to be addressed.
Abernethy (1994a,b) analyzed the same samples, using two protein loci, two microsatellite loci, and mtDNA, and argued for extensive hybridization between sika and red. Our interpretation is that there is only a low rate of hybridization. This discrepancy arises because the two protein loci scored by Abernethy were polymorphic within even the northernmost samples of red deer. These loci were treated as diagnostic and were included in a single “hybrid index”; thus, a high frequency of hybrids was inferred well to the north, and the clines for mtDNA and nuclear loci appeared to be staggered. Our analysis makes no assumptions as to the origin of each allele, and instead distinguishes current hybridization on the basis of spatial patterns of allele frequency and linkage disequilibria. In general, alleles should never be treated as absolutely diagnostic: even if samples at opposite ends of a transect are fixed, the most that can be said is that the actual allele frequency is likely to be less than some small value. Indeed, our analysis suggests that some sika alleles are present at low frequency within the ancestral population of red deer.
Abernethy’s (1994a,b) analysis treated the red- and sika-like deer as a single population, with strong linkage disequilibrium and heterozygote deficit. Such an approach is appropriate where extensive hybridization has generated a wide range of recombinant genotypes, and as such was applicable under the original interpretation of the data. Our approach, in which the population is divided into two parts, is valid only when hybridization is rare and when F1’s tend to backcross in one or other direction. In this case each part can be seen as consisting of a mixture of nth generation backcrosses, each containing a proportion 2-n of introgressed alleles. While individuals cannot be reliably classified to any particular backcross generation, this genetic structure allows a simple likelihood-based estimation of rates of hybridization and introgression.
To determine whether the rates of current hybridization inferred from linkage disequilibria are consistent with the observed extent and spatial pattern of introgression, we constructed a simple diffusion model for the spread of sika and red deer. Our estimates that sika increase from low density at r ∼ 9.2% yr-1, and that the rate of dispersal is σ ∼ 3.7 km in each year, are robust to assumptions about the initial introduction, the time when red deer moved south, and the nature of competition between red and sika. The rate of increase is consistent with the estimate that by the 1940’s, sika had increased to 300–400 individuals (Whitehead 1950, 1964; Ratcliffe 1987): exponential increase at a rate of 9.2% yr-1 over 30 yr, from 15 individuals initially, would give 240 individuals after 30 yr. In making these estimates, we ignored the effect of hybridization on demography. This is reasonable because hybridization is so rare and F1’s must be reasonably fertile: the spread of the sika phenotype is determined primarily by ecological rather than genetic factors.
In the current model competitive interactions between red and sika appear to have had minimal effect on the spread of the populations in Argyll. This is because the two species have met only recently. As the contact time increases, the competitive advantage of each species in different environments may become important in determining further population expansion and introgression. Outside the current sampling area there are large tracts of upland, open-hill habitat with large red deer populations, which sika have so far shown little sign of colonizing (Ratcliffe 1987). Few data are available on how sika might interact with red deer in such environments, but they are likely to differ from forestry data, where interactions are better known (Chadwicket al. 1996). If a competitive advantage to sika is an important factor, then sika may ultimately displace red deer from the environments in which they compete. If red and sika populations equilibrate, neutral loci must inevitably become completely introgressed.
Nonrandom mating has kept red and sika deer distinct during the ∼30 yr since they came into contact in Kintyre. Natural selection may also be maintaining these separate populations: our analysis suggests that introgression is 30–40% of that expected with no selection. Selection might act against hybrids (“endogenous”) or against alleles in a foreign environment (“exogenous”); in this case, the distinction is not clearcut because of the differences in habitat use between red and sika (Abernethy 1994a,b). If selection against hybrids is occurring, we speculate that it must act at the gene level because, although red and sika karyotypes differ by two Robertsonian rearrangements (2n red = 68, 2n sika = 64), there is no evidence for nondisjunction in captive hybrids (Herzog and Harrington 1991). Selection and assortative mating may maintain distinct phenotypes indefinitely, despite hybridization. However, alleles that are advantageous within the foreign genetic background will penetrate relatively rapidly, and even neutral alleles will eventually intermingle (Barton and Bengtsson 1986). Thus, even if sika do not displace red deer entirely, the introgression of adaptive traits from sika into red could have substantial consequences for the latter. The outcome depends primarily on whether single genes can confer a selective advantage within the genetic background and environment of the opposite type.
In other areas where red and sika have been in contact, such as County Wicklow in Ireland, we must assume that the barriers to gene flow posed by assortative mating and selection were not maintained. Here, complete genetic and phenotypic introgression appears to have taken place (Harrington 1973, 1979, 1982). In Wicklow, hybridization between red and sika has been taking place for ∼80 yr longer than in Argyll (Powerscourt 1884). Apart from the date of initial contact, little information is available to make a comparison between Wicklow and Argyll. With the current data, we cannot discount the possibility that introgression of red and sika will ultimately progress to the same extent in Argyll.
We are grateful to Forest Enterprise in Argyll for providing the samples used in this study. We also thank Loeske Kruuk plus the communicating editor and two anonymous referees for their helpful comments on the manuscript. This work was supported by a Natural Environment Research Council grant to N.B. and J.P. and by a University of Edinburgh postgraduate bursary to G.S.
Communicating editor: A. G. Clark
- Received July 22, 1998.
- Accepted January 19, 1999.
- Copyright © 1999 by the Genetics Society of America