Volcanic islands represent excellent models with which to study the effect of vicariance on colonization and dispersal, particularly when the evolution of genetic diversity mirrors the sequence of geological events that led to island formation. Phylogeographic inference, however, can be particularly challenging for recent dispersal events within islands, where the antagonistic effects of land bridge formation and vicariance can affect movements of organisms with limited dispersal ability. We investigated levels of genetic divergence and recovered signatures of dispersal events for 631 Galápagos giant tortoises across the volcanoes of Sierra Negra and Cerro Azul on the island of Isabela. These volcanoes are among the most recent formations in the Galápagos (<0.7 million years), and previous studies based on genetic and morphological data could not recover a consistent pattern of lineage sorting. We integrated nested clade analysis of mitochondrial DNA control region sequences, to infer historical patterns of colonization, and a novel Bayesian multilocus genotyping method for recovering evidence of recent migration across volcanoes using eleven microsatellite loci. These genetic studies illuminate taxonomic distinctions as well as provide guidance to possible repatriation programs aimed at countering the rapid population declines of these spectacular animals.
THE evolution of genetic distinctiveness among island populations is the long-term balance of colonization history, environmental variation, and subsequent patterns of gene flow associated with geological and historical changes of geographical barriers to dispersal (Bowman et al. 1983; Berry 1984; Wagner and Funk 1995; Grant 1998; Emerson 2002; Jordan et al. 2005). This is particularly true for volcanic islands. In fact, although historical processes such as secondary contact between previously isolated species or multiple colonization events (e.g., Wright and Simovich 1985; Grant and Grant 1997) may complicate phylogeographical inference, a number of biogeographical scenarios in which species phylogenies mirror the sequence of paleogeographical events have been described (Wright 1983; Thorpe et al. 1994; Tarr and Fleischer 1995; Rassmann 1997; Caccone et al. 2002; Hormiga et al. 2003; Beheregaray et al. 2004). Inferring colonization dynamics of island biota from genetic data is more difficult for populations or lineages deriving from relatively recent dispersal events within islands, where the antagonistic effects of land bridge formation and vicariance can affect movements of organisms with limited dispersal ability (Emerson et al. 1999; Juan et al. 2000).
The Galápagos giant tortoise Geochelone nigra (or G. elephantopus; see Zug 1997) is a vivid example of taxon evolution following a single colonization event, radiation across islands as they were formed by volcanic activity, and subsequent divergence under restricted gene flow. Previous genetic studies reconstructed the evolutionary history of island populations and clarified most taxonomic issues (Caccone et al. 2002; Ciofi et al. 2002; Beheregaray et al. 2004; Russello et al. 2005). The authors identified two ancestral lineages on the islands of Española and San Cristóbal, three lineages (initially assigned to a single taxon) on the island of Santa Cruz (Beheregaray et al. 2003; Russello et al. 2005), two distinct taxa on Pinzón and Santiago, and a phylogenetically younger and more complex group distributed across the five volcanic areas on the largest island of Isabela. Although phylogenetically distinct units were recognized for central and northern volcanoes of Isabela, genetic analyses (Caccone et al. 2002; Beheregaray et al. 2004) did not recover significant lineage sorting for tortoises distributed across the southern volcanoes Sierra Negra and Cerro Azul (Figure 1). These populations most probably represent one of the most recent colonization events by Galápagos tortoises (Beheregaray et al. 2004) such that their phylogenetic patterns are more difficult to recover than are those for older populations. Moreover, demographic fluctuations linked to volcanic activity, historical exploitation of tortoises, and recurrent gene flow between populations on adjacent volcanoes might have obscured identification of past and recent events that determined current population structure and levels of genetic distinctiveness. Patterns of morphological differentiation are equally intriguing. Tortoises from Sierra Negra and Cerro Azul were originally assigned to the taxon guntheri and vicina, respectively. This classification, based on minor morphological characteristics of the carapace and on differences in body size at adulthood, is subject to debate (Pritchard 1996; Ernst et al. 2000) and probably linked to gradients of environmental conditions that vary with altitude and aspect (Colinvaux 1972; Fritts 1984).
The southern region of Isabela Island therefore offers an ideal setting in which to investigate in detail genetic patterns of colonization and dispersal linked to relatively recent geological events and to clarify unresolved taxonomic issues. Moreover, in the past century, tortoises from both Cerro Azul and Sierra Negra have come under increasing threat due to habitat degradation, human depredation, and predation by feral animals, particularly goats (Powell and Gibbs 1995; Kaiser 2001). These factors have prompted establishment of captive breeding programs (Cayot et al. 1994). Knowledge of the level and pattern of population differentiation is therefore crucial for guiding implementation of repatriation plans using tortoises from captive stocks.
In this study, we integrate analyses of microsatellite loci and mitochondrial DNA control region sequences to investigate the genetic structure of Galápagos giant tortoises from southern Isabela. Mitochondrial gene genealogies can now be analyzed using algorithms that attempt to separate the effects of past events (e.g., fragmentation and colonization) from contemporary processes, such as recurrent gene flow (Templeton 1998), and estimate timing of historical demographic changes (Slatkin and Hudson 1991; Rogers and Harpending 1992; Fu 1997). Genetic analyses of dispersal and migration based on microsatellite markers, on the other hand, have often relied on assumptions rarely met by the actual ecological and demographic history of populations (Neigel 1997; Waples 1998; Whitlock and McCauley 1998; Balloux and Lugon-Moulin 2002; Neigel 2002). However, methods based on genealogical history of alleles (e.g., Slatkin 1991; Hudson 1998), maximum-likelihood estimates of gene flow (Rannala and Hartigan 1996), and those integrating maximum-likelihood techniques and coalescent theory (Bahlo and Griffiths 2000; Beerli and Felsenstein 2001; Vitalis and Couvet 2001) have relaxed assumptions concerning either temporal or intraspecific changes of migration rates and population size. Other methods, implementing Bayesian approaches for individual assignment and immigrant identification, are not conditioned upon migration-drift equilibrium (Paetkau et al. 1995; Rannala and Mountain 1997; Cornuet et al. 1999; Gaggiotti et al. 2002) and do not require prior knowledge of the distribution of sampling locations (Pritchard et al. 2000; Dawson and Belkhir 2001; Corander et al. 2003).
We infer historical patterns of colonization, population dispersal, and demographic changes for tortoises of Sierra Negra and Cerro Azul using nested clade analysis and mismatch distribution of mitochondrial DNA sequences in relation to known geological history. Evidence of recent migration across volcanoes was recovered by microsatellite analysis, implementing a Bayesian multilocus genotyping method developed by one of the authors (Wilson and Rannala 2003). The algorithm requires fewer assumptions than other model-based Bayesian techniques developed for assigning individual genotypes to populations of origin and, in particular, relaxes the key assumption of Hardy–Weinberg equilibrium within populations. By integrating analytical processes that describe evolutionary patterns under different temporal frameworks we can compare patterns of genetic diversity between tortoise populations from different sampling sites of Sierra Negra and Cerro Azul. We also discuss suggested distinctions based on distribution of named taxa in view of possible demographic reinforcement programs and other conservation measures for the tortoises of southern Isabela.
MATERIALS AND METHODS
Sampling and genetic techniques:
Blood samples were collected from Galápagos giant tortoises on the volcanoes Sierra Negra (guntheri) and Cerro Azul (vicina) (Figure 1). At Sierra Negra, samples were obtained from La Cazuela, a site on the eastern flank of the volcano, and from Cabo Rosa and Roca Union, on the southern slope. At Cerro Azul, tortoises were sampled on the western and eastern slopes of the caldera. Samples were also obtained from the Cinco Cerros region, in the southeast, where both vicina and guntheri have been observed. Blood was drawn from the brachial vein of the front leg of tortoises, preserved at ambient temperature in a lysis buffer (0.1 m Tris–HCl, 0.1 m NaCl, 5 mm EDTA, 0.5% SDS), and subsequently used for whole DNA extraction as decribed in Ciofi et al. (2002). A total of 631 tortoises were analyzed at 11 microsatellite loci, with an average of 78.8 ± 6.9 (SE) individuals per sampling site. Part of this sample set (359 tortoises, 44.8 ± 3.5 individuals per sampling site) was used for mitochondrial DNA (mtDNA) sequences analysis. A stratified random sampling of the Cerro Azul data set reduced the number of samples analyzed for mtDNA sequences variation and minimized discrepancy in sample size across locations for nested clade analysis. The gender of sampled tortoises was estimated by looking at the shape of the plastron, tail length, and relative position of the cloaca.
Ciofi et al. (2002) described the construction of a genomic library enriched for dinucleotide repeats and documented allelic variation at 10 polymorphic loci in Galápagos tortoises. Significant deviation from Hardy–Weinberg proportions was recorded at 8 loci for which a relatively high frequency of null alleles was suggested. To avoid bias in allele frequency estimation, we used four of these loci (Gal45, Gal75, Gal136, and Gal263) with no previous evidence of possible mismatch between primers and microsatellite flanking sequences and a new set of seven microsatellites isolated from the same genomic library (Table 1). Polymerase chain reaction (PCR) conditions and thermal profiles, electrophoresis of PCR products, and assessment of allele size were as described in Ciofi et al. (2002).
A DNA fragment of 705 bp of the mitochondrial DNA control region was amplified by PCR using primers designed by Caccone et al. (1999). Sequence polymorphism among individuals was assessed by comparing single-stranded conformation polymorphisms (SSCP) among individual tortoise samples. PCR amplification, detection, and sequencing of representative SSCP gel phenotypes were carried out as described in Beheregaray et al. (2003). We assessed the reliability of our SSCP protocol by using samples from ∼200 individuals of known sequences (Caccone et al. 2002; Beheregaray et al. 2003) as controls in all gels. We sequenced multiple individuals with the same SSCP gel band and showed that samples with the same SSCP phenotype as the control had the same sequence, while individuals with different phenotypes had unique sequences (Beheregaray et al. 2004).
Analysis of mitochondrial DNA variation:
Mitochondrial sequence variation was estimated by haplotypic (gene) diversity, number of polymorphic sites, mean number of nucleotide differences between haplotypes, and nucleotide diversity (Nei 1987) using ARLEQUIN (Schneider et al. 2000). We obtained maximum-likelihood values for different models of sequence evolution using the subroutine provided by ModelTest 3.6 (Posada and Crandall 1998) into PAUP 4.0b10 (Swofford 2003). The Akaike information criterion implemented in ModelTest inferred the Tamura and Nei (TN) model (Tamura and Nei 1993) as the most likely model of DNA substitution for our data set, with γ distribution shape parameter of 0.62 and a proportion of invariable sites of 0.87. We estimated the extent of mtDNA control region differentiation from the average number of pairwise differences among sampling sites and the parameter ΦST (Excoffier et al. 1992) using the TN molecular distance measure implemented in ARLEQUIN. Significance of pairwise sequence differences and ΦST values were obtained after 10,000 halpotype permutations.
Patterns of genetic structure across volcanoes were investigated using a hierarchical analysis of molecular variance (AMOVA) implemented in ARLEQUIN. Total genetic variance was partitioned into covariance components (see Excoffier 2003) that were used to compute fixation indexes as measures of the degree of genetic differentiation. Significance of fixation indexes was tested by 10,000 permutations of individual genotypes among sampling sites or sampling sites among groups as appropriate (Schneider et al. 2000). The latter were defined by pooling sampling sites according to volcanic area and geographical distance. Analyses were performed both including and excluding the number of pairwise nucleotide differences among haplotypes (Excoffier et al. 1992).
Haplotype network and nested clade analysis:
A nested clade procedure was implemented to assess geographical associations of haplotypes and infer historical patterns of colonization and dispersal. A network of haplotypes was first constructed using statistical parsimony (Templeton et al. 1992). The method, implemented in the TCS program (Clement et al. 2000), links first haplotypes with the smaller number of differences as defined by a 95% confidence criterion and identifies the most probable ancestral haplotype according to coalescence theory (Castelloe and Templeton 1994). Haplotypes were then organized into a system of nested clades where a higher nesting level corresponds to longer evolutionary time (Templeton et al. 1992).
Geographical association between haplotypes was first assessed using the nested contingency test described in Templeton and Sing (1993) by permuting clade types within a nested category against sampling locations (considered as categorical variables). The clade distance (DC) and the nested clade distance (DN) were then calculated from geographical distances among sampling sites obtained using a 1:100,000 topographic map of southern Isabela (Instituto National Galápagos 1989). The null hypothesis of no geographical association within nested clades was tested using GEODIS 2.0 (Posada et al. 2000) by comparing observed values of DC and DN between tip clades (with one mutational connection), interior clades (with two or more connections), and interior-tip clade comparison (I-T) within a nested group of clades with a null distribution derived from 10,000 Monte Carlo permutations of clades against sampling sites. The interpretative key given by Templeton (2004) was then used to infer recurrent and historical events from patterns of statistically significant distance measures.
For each location, we tested whether the frequency distribution of the observed pairwise differences between mtDNA haplotypes was statistically different from a distribution expected under a model of demographic expansion (Rogers and Harpending 1992). The mismatch distribution is expected to follow a unimodal pattern for samples under population change and a multimodal pattern for populations at demographic equilibrium (Rogers and Harpending 1992). The sum of square deviations (SSD) between the observed and the expected distribution and the raggedness index r of the observed distribution of the mismatch classes (Harpending 1994) were computed as test statistics under the null hypothesis of population growth using the parametric bootstrap approach of Schneider and Excoffier (1999) implemented in arlequin (Schneider et al. 2000). We then estimated the number of generations passed since range expansion (t = τ/2u) using τ, a measure of time to expansion in terms of mutational differences accumulated by a random pair of individuals, and u, the number of substitutions per generation for the mtDNA control region under study (Slatkin and Hudson 1991; Rogers and Harpending 1992). Rates of mtDNA nucleotide substitution are difficult to calibrate to meaningful units because estimates of gene coalescence often predate timing of actual population splits (Arbogast et al. 2002; Graur and Martin 2004). A number of correction factors have been proposed to account for this pattern (Edwards and Beerli 2000; SigurÐardóttir et al. 2000; Ruokonen and Kvist 2002; Hoffman and Blouin 2004; Near et al. 2005). We calculated u = μL, where μ is the mutation rate per nucleotide per generation estimated assuming a rate of evolution of 3.4% per million years (Caccone et al. 2002; Beheregaray et al. 2004) and L is the length (705 bp) of the control region sequence analyzed.
Analysis of microsatellite variation:
We assessed the average number of microsatellite alleles per locus, observed heterozygosity, unbiased gene diversity (Nei 1987), and tested for heterozygote deficiency (Guo and Thompson 1992) using GENEPOP 3.3 (Raymond and Rousset 1995). Departure from Hardy–Weinberg expectations was also tested by comparing, for each sampling location, the observed FIS value to a frequency distribution of FIS indices obtained after 10,000 permutations of alleles performed with GENETIX 4.01 (Belkhir et al. 2000). A sequential Bonferroni correction adjusted critical probability values for multiple tests to minimize type I errors (Rice 1989; Sokal and Rohlf 1994). We tested for the presence of null alleles by assuming that homozygote excess due to null heterozygotes and a large number of nonamplifying DNA samples caused by null homozygotes are an indication of a high frequency of nulls at a specific locus (Beaumont et al. 2001). For those locations where evidence of heterozygote deficiency was observed we correlated FIS values of each locus to the frequency of nonamplifying samples using a Spearman rank-order correlation test.
Genetic differentiation was assessed using GENETIX by the index θ (Weir and Cockerham 1984; Belkhir et al. 2000) and by ρST, an estimator of differences in allele sizes transformed to standard deviations from the global mean of repeat unit numbers using RST 2.2 (Goodman 1997). The distribution of private (unique) alleles was assessed as a mean to describe population distinctiveness (Slatkin 1985) and levels of genetic introgression between taxa (e.g., Gottelli et al. 1994; Beaumont et al. 2001). Patterns of relative similarities among sampling sites were also resolved by computing the principal components that best represented the original relationships among sampling sites and alleles (considered as the correlated variables). Principal component analysis was performed using PCA-Gen (Goudet 1999) on multilocus genotypes. Patterns of genetic structure across volcanoes were investigated using AMOVA as for mtDNA analysis.
Evidence of recent migration events across southern Isabela was assessed using the Bayesian multilocus genotyping procedure implemented with Markov chain Monte Carlo (MCMC) methods in BayesAss (Wilson and Rannala 2003). This method does not require populations to be in either migration drift or Hardy–Weinberg equilibrium. To examine the strength of the information in the tortoise microsatellite data set, 95% confidence intervals were determined for migration rates and compared to a scenario where all proposed changes throughout the Markov chain are accepted (thereby simulating the event where any information that may exist in the data is insufficient to affect the posterior distribution of migration rates). This option is available in BayesAss 1.3 at http://www.rannala.org. The MCMC was run for a total of 3.0 × 106 iterations, with the first 106 iterations discarded as burning to allow the chain to reach stationarity. Samples were collected every 2000 iterations to infer posterior probability distributions of parameters of interest.
Thirty-six mtDNA haplotypes were identified, defined by 67 polymorphic sites. Eight haplotypes were found only once, while the most common was shared by 64 individuals in La Cazuela and eastern Cerro Azul (Table 2). Mean number of pairwise differences was 2.53 ± 0.53 with the lowest value (0.41 ± 0.22) recorded in La Cazuela and the highest one (4.07 ± 2.07) reported for Cinco Cerros. Low values of haplotype (h) and nucleotide (π) diversity and a relatively small number of polymorphic sites (p) were also recorded in La Cazuela. Conversely, haplotypes from Cinco Cerros showed relatively higher values of h, π, and p (Table 3). Most sampling locations shared either none or one haplotype with any one of the other locations. Cinco Cerros shared four haplotypes with western Cerro Azul and two haplotypes with Roca Union, but had 60% unique haplotypes.
Different levels of polymorphism were observed across microsatellite loci. Number of alleles per locus varied from 8 to 20. The level of genetic diversity in the populations from different sampling sites, on the other hand, was remarkably similar (Table 3). Average allele diversity (8.1 ± 0.2 SE) and mean gene diversity (0.70 ± 0.05) did not differ significantly across locations (ANOVA; F = 0.69, P = 0.68; F = 0.75, P = 0.63, respectively). Departure from Hardy–Weinberg expectations was detected as heterozygote deficiency for Roca Union (FIS = 0.07; P < 0.01) and all three sampling locations of Cerro Azul (FIS = 0.03–0.15; P < 0.01). Significant deviation from equilibrium was still recorded for eastern Cerro Azul and Cinco Cerros at two and three loci, respectively, after adjustment for multiple comparisons. Tortoises from Isabela, a relatively young island, represent the most recent dispersal event within the overall pattern of island colonization in the Galápagos (Beheregaray et al. 2004). They have also been subject to a plethora of stochastic and deterministic threats that have shaped and are currently affecting the demography of extant populations (see discussion). Moreover, deviation from Hardy–Weinberg expectation was observed, on average, for 50% of the loci examined, suggesting recurrent demographic changes rather than null alleles as the most likely explanation for the observed departure from equilibrium. This hypothesis was tested by correlating the number of samples that did not amplify to the FIS values for each locus of Roca Union and the sampling sites of Cerro Azul. No significant association was found for any of the samples (rs = −0.269–0.383; P = >0.05).
Mitochondrial DNA differentiation among sampling sites:
Significant genetic divergence at mtDNA control region sequences was recorded among sampling sites (P < 0.01). Average corrected (Schneider et al. 2000) percentage sequence divergences among locations varied from 0.07% between Cabo Rosa and Cinco Cerros to 1.0% between Roca Union and eastern Cerro Azul. Tortoises from La Cazuela had 0.23% average sequence divergence from eastern Cerro Azul and 0.94%–1.2% from all other sampling sites. The same pattern of genetic divergence was reflected by ΦST values, which varied from 0.11 to 0.85 and were all significant at the 1% level.
The nested AMOVA revealed low levels of mtDNA genotypic structuring across southern Isabela (Table 4). We tested for two grouping schemes. The first design had one group including sampling sites from Cerro Azul and Cinco Cerros and another group with Cabo Rosa, Roca Union, and La Cazuela assigned to Sierra Negra. The second design had Cabo Rosa pooled with Cerro Azul and Cinco Cerros given similarities between tortoises from eastern Cerro Azul and southwestern Sierra Negra identified by preliminary analyses (Ciofi et al. 2002). Variation among individuals within and among sampling sites within groups was always significantly higher (P < 0.01) than that recorded among groups for either of the two pooling schemes, indicating that genetic diversity among haplotypes was unlikely to depend on specific geographical assemblage of sampling sites. However, most mtDNA variation was recorded among sampling sites and, to a minor extent, among tortoises within sampling locations.
Colonization and patterns of dispersal in southern Isabela:
Statistical parsimony identified mtDNA haplotype 1 as ancestral, suggesting that La Cazuela, where haplotype 1 was found at a very high frequency (88%), was probably the first colonization site of southern Isabela (Table 2, Figure 2). Haplotypes 2 and 3, also sampled in La Cazuela, were just two and one mutational steps away, respectively, from the sequence with the largest root probability. A first, visual inspection of the network and the geographic distribution of haplotypes suggested a colonization event of southern Sierra Negra and of Cerro Azul (shaded and open circles in Figure 2, respectively) with no subsequent backmigration. The scattered spatial pattern of the other haplotypes (Templeton 1998) would imply instead recurrent interchange between Cerro Azul and southern Sierra Negra.
MtDNA haplotypes separated by up to 11 mutational steps were connected in a single parsimony network with 95% probability. The three internal ambiguities found in the cladogram (wiggling links in Figure 2) were solved by applying criteria derived from empirical validation of the coalescent theory (Crandall and Templeton 1993; Templeton and Sing 1993; Fetzner and Crandall 2003). Accordingly, a haplotype connection was maintained to high-frequency (older) haplotypes, with an interior position in the network (Castelloe and Templeton 1994), rather than to low-frequency ones, located in tip clades. Also, connections were preferentially given between haplotypes occurring in the same geographical area because newly formed haplotypes, in particular singletons (with a frequency of 1), have a tendency to remain close to the population of origin (Posada and Crandall 2001).
Figure 2 shows the haplotype categories arranged in a nested series of clades according to the algorithm of Templeton et al. (1987). Four nesting levels were defined, with 18 clades that could be tested for a geographical association. Clade 24 was the only observed haplotype symmetrically stranded and was grouped with the nesting category that had the smallest sample size. Stranded haplotypes that were missing intermediates were left unnested (Templeton and Sing 1993).
Nested contingency analysis recovered a strong overall association between clades and geographical distribution (Table 5). Only two clades (1.1 and 3.4) did not show significant geographical association of haplotypes. This could be due to panmixia, absence of historical demographic changes, insufficient genetic variation, or inadequate sampling (Templeton et al. 1995). The last two hypotheses are unlikely because genetic diversity indexes were relatively high compared to other island populations (Ciofi et al. 2002; Beheregaray et al. 2004). Sample size (n = 359) was adequate relative to a population size of approximately 900 tortoises estimated for southern Isabela (MacFarland et al. 1974; Pritchard 1996), and sampling design included all but two of the extant tortoise populations in southern Isabela (Cerro Paloma, in west Sierra Negra, and Gavilanes, in north Cerro Azul, both with very few individuals remaining). A more plausible explanation for the absence of significant geographical association could be low gene flow between sampling sites of clade 1.1 (eastern and western Cerro Azul) and between sites of clade 3.4 (Cinco Cerros and southwest Sierra Negra).
Distance measures and probability values for clades showing geographical and genetic variation are reported in Figure 3. Interior-tip status could not be determined a priori for the two major clades (4.1 and 4.2). Although the ancestral status of haplotype 1 from La Cazuela embedded in clade 4.1 would make the interior state of this clade more plausible, we performed two separate NCA with either of the clades as interior. The oldest inferred event, considering interior status for clade 4.1, was gene flow restricted by isolation by distance from Sierra Negra to Cerro Azul. Range expansion in the same direction was instead inferred when clade 4.2 was set as the interior clade (Table 5). Range expansion and long-distance dispersal from Sierra Negra to Cerro Azul were recovered by distance analyses for clades 4.1 and 4.2 and a similar pattern, inclusive of isolation by distance, was detected by the analysis of three-step clades. Analysis of two-step clades revealed backmigration events from Cerro Azul to southwest Sierra Negra and secondary dispersal from Sierra Negra (La Cazuela) to Cerro Azul through a northern route. Dispersal and restricted migration from Cerro Azul to southern Sierra Negra was inferred by one-step clade analysis along with a recent dispersal from southern Sierra Negra back to western Cerro Azul.
Mismatch analysis recovered a strong signal of demographic stability for one of the southern sampling sites of Sierra Negra (Cabo Rosa). The distribution of the number of pairwise differences among haplotypes of tortoises from Cabo Rosa was multimodal (data not shown) and significantly different from the unimodal pattern expected under a model of population expansion (SSD = 0.08, P < 0.01; r = 0.17, P < 0.01). A similar mismatch curve was found for western Cerro Azul, with a raggedness index larger than the one predicted for a unimodal distribution, but with an SSD value not significantly different from expectation (P = 0.09). La Cazuela, Roca Union, Cinco Cerros, and eastern Cerro Azul all showed unimodal mismatch distributions not different from the one expected under a model of demographic growth (SSD = 0.01–0.06, P > 0.1; r = 0.02–0.13, P > 0.05). The initiation of each range expansion was estimated using a generation time for Galápagos giant tortoises of 25 years (H. Snell, unpublished results from captive breeding programs). On the basis of the estimated τ values and the range of estimated mtDNA mutation rates (Caccone et al. 2002; Beheregaray et al. 2004), time since population expansion in southern Isabela was 63,000–144,000 years in Sierra Negra and 20,000–131,000 years in Cerro Azul, with the most recent expansion event estimated for tortoises in eastern Cerro Azul.
Patterns of microsatellite variation and gene flow:
Differences in microsatellite allele frequencies and allele sizes between locations revealed a pattern similar to that recorded for mtDNA. Values of θ and ρST varied between 0.01 and 0.26 (P < 0.01). A similar pattern was evident when male and female samples were analyzed separately (data not shown). Private alleles were recorded for all sampling sites, with a frequency ranging from 0.005 to 0.14. Frequencies >0.05 were observed for Sierra Negra (La Cazuela and Roca Union) for one allele at three loci. Our data set from Cerro Azul was obtained during two separate surveys, each conducted on the eastern (Las Pegas and Los Crateres) and the western (Las Pampas and Las Tablas) slopes of the caldera. We compared samples obtained during these surveys to determine whether population substructuring could have accounted for the observed departure from Hardy–Weinberg equilibrium (Walhund effect) in eastern and western Cerro Azul. No significant differentiation (P > 0.05) was observed between the samples from Las Pegas and Los Crateres (θ = 0.002; ρST = 0.0005) and those from Las Pampas and Las Tablas (θ = 0.005; ρST = 0.009).
Principal components analysis showed a strong association between geographical distance and genetic distinctiveness (Figure 4). The first three components of PCA explained >85% of the total variation, with components 1 and 2 containing 46 and 22% of the variation, respectively. A marked distinction of La Cazuela and Roca Union from the other sampling sites was evident from the first and third component, respectively, while the second component highlighted a pattern of similarity between eastern Cerro Azul, Cinco Cerros, and, to a minor extent, Cabo Rosa.
In the AMOVA (Table 4), significant genotypic diversity was detected both within and among sampling sites, while none of the groupings tested in the analysis was different from each other. This pattern corroborated the results from mitochondrial DNA; however, >93% of genotypic variation was distributed among individuals within sampling locations while differentiation among groups and among sampling sites accounted for <5% of the total molecular variance. Values of F-statistics were also much lower than Φ estimates recorded for mtDNA.
The mean posterior probabilities and 95% confidence intervals for migration rates between sampling sites are reported in Table 6. When all proposed changes were accepted in the Markov chain, simulating the effect of having no information in the data from which to estimate migration rates, we obtained a 95% confidence interval of approximately (0.67, 0.99) for nonmigration rates and of (0, 0.16) for migration rates. Confidence intervals recovered from this data set were considerably smaller than those obtained from the above scenario, suggesting that our microsatellite data set contained an appreciable amount of information from which to estimate migration rates and other parameters of interest. Most sampling locations in southern Isabela had very low proportions of migrants. Cabo Rosa, however, had a fairly large proportion of migrants from Cinco Cerros (m = 0.304). A positive value for the lower 95% confidence interval was also recorded for the migration rate from eastern Cerro Azul to Cinco Cerros, although in this case the proportion of migrants was relatively low (m = 0.083). Migration rates from Cabo Rosa to Cinco Cerros and from Cinco Cerros to eastern Cerro Azul were quite small, suggesting recent source-sink movement patterns from northwest toward Cabo Rosa.
Island formation and tortoise phylogeography:
The Galápagos archipelago is a striking example of a chronological association between orogenic events and organismal evolution, for patterns of island colonization recovered using gene genealogies have been shown to mirror sequential volcanic activities that led to the formation of the present islands (Wright 1983; Rassmann 1997; Caccone et al. 2002; Beheregaray et al. 2004). For Galápagos tortoises, dispersal appears largely to have started from the oldest, southeastern formations (Hall 1983; White et al. 1993) with little evidence of recurrent gene flow among island populations (Caccone et al. 1999; Caccone et al. 2002; Beheregaray et al. 2004; Russello et al. 2005). The island of Isabela, in particular, was subject to two colonization events. One from the island of Santiago giving rise to the population on the northernmost volcano Wolf (Caccone et al. 2002). Another from Santa Cruz resulted in the colonization of Sierra Negra (Beheregaray et al. 2004). Our results from nested clade analysis supported these findings in that they identified the most frequent haplotype of the eastern sampling site of Sierra Negra, La Cazuela (haplotype 1 in Figure 3), as ancestral in southern Isabela. In previous studies (Beheregaray et al. 2004), ancestral haplotypes were found instead at Roca Union (southeastern Sierra Negra). Differences in sample sizes, number of populations, and geographic areas analyzed probably affected network resolutions and, as a consequence, may have resulted in slightly different nested clade design in the two studies. Nevertheless, our analysis advocates for a pattern by which Sierra Negra was proposed as the source of subsequent dispersal to volcano Alcedo and Darwin, in the central and northern part of the island, and to Cerro Azul (Beheregaray et al. 2004). Tortoises from La Cazuela, as from the majority of the other sampling sites in southern Isabela, did not present a mismatch distribution of haplotypes characteristic of other populations at demographic equilibrium (Beheregaray et al. 2003). This may be explained by the relatively recent colonization of Isabela. Timing of the demographic expansion of the La Cazuela tortoise population is consistent with the estimated time of emergence of Sierra Negra. Similarly, divergence times between tortoises from La Cazuela and other sampling sites fit the inferred sequence of volcanic events in southern Isabela. According to the observed average percentage of sequence divergence and mutation rates estimated for the tortoise mtDNA control region, separation between La Cazuela and eastern Cerro Azul tortoise populations dates back 67,000 years, while 276,000–353,000 years was inferred as the split between tortoises from La Cazuela and those from all other locations. Nordlie (1973) suggested a sequence of evolutionary stages for volcanoes of Isabela ranging from submarine formation to caldera development and decline and identified Sierra Negra and Cerro Azul as the oldest and youngest formations, respectively. After volcanic emergence, lava flows and subaerial growth between the two volcanoes most probably provided stepping stones for tortoise movements (Beheregaray et al. 2004). We identified two main dispersal events, a first via southern Sierra Negra (Cabo Rosa and Roca Union) and a second, more recent one to eastern Cerro Azul (Figure 5). The latter is advocated by the occurrence of the ancestral haplotype in tortoises from eastern Cerro Azul and the relatively recent divergence time from the La Cazuela tortoise population.
The remaining events of range expansion and dispersal are probably connected to aerial volcanism shaping the superficial geology of southern Isabela in the past 0.3 million years. For instance, recurring eruptions over the southeastern apron of Sierra Negra (Reynolds et al. 1995) may explain genetic divergence between La Cazuela and Roca Union. Similarly, dispersal from Cerro Azul back to south Sierra Negra followed by two main migratory events to and from Cerro Azul (Table 5 and Figure 5) were probably dictated by the opportunity for tortoises to move across a revegetating environment between periods of volcanic eruption, particularly along the zone of coalescence of Sierra Negra with Cerro Azul (Reynolds et al. 1995). Timing of the divergence between mtDNA lineages in Cerro Azul and southern Sierra Negra varies between 294,000 years ago for Roca Union and eastern Cerro Azul east to just 20,000 years between Cinco Cerros and Cabo Rosa. Although the timing of population divergence is affected by the choice of molecular clock, these estimates appear fairly consistent with the geological history of southern Isabela. On the other hand, rates of mtDNA evolution should not substantially alter hierarchical nested analysis that, in accordance with coalescent theory, infer dispersal dynamics merely on the basis of temporal differences between interior, older clades and younger clades located at the tips of a tree. The accuracy of nested clade designs has been tested on species with prior expectations from known phylogeographical history (Templeton 2004 and references therein). Moreover, inference of historical demography has been improved to account for deterministic events such as haplotype extinction (Masta et al. 2003), stochasticity of the mutational process that may hinder detection of a particular event in time (Templeton 2004), and particular biogeographical settings that may render discrimination of similar historical events difficult (Hutchison and Templeton 1999; Knowles and Maddison 2002). In this study historical misinterpretation of nested clade analysis that can arise from inappropriate sampling design (Templeton 2001; Knowles and Maddison 2002) was minimized by an exhaustive sampling scheme including a relatively large sample size for all locations of southern Isabela known to harbor a sizeable number of giant tortoises.
Dispersal and population structuring:
For giant tortoises, dispersal between islands represents a stochastic event dependent on local oceanic currents (Pak and Zaneveld 1973; Caccone et al. 2002), so that backmigration after population establishment is unlikely to occur. Previous phylogeographic studies, in fact, identified a strong pattern of lineage sorting defining monophyletic clades for older Galápagos islands and clear genetic and demographic distinction for tortoise populations on islands of intermediate age (Caccone et al. 2002; Beheregaray et al. 2003). Significant genetic divergence was also found in the island of Isabela among tortoise populations ranging from central to northern volcanoes. Tortoises from southern Isabela, on the other hand, have so far proven difficult to distinguish on both genetic and morphological grounds (Caccone et al. 2002). The relatively young age of Cerro Azul and the pattern of dispersal events reported in this study may explain the relatively incomplete lineage sorting. Molecular analysis of variance also advocates for this hypothesis. Approximately 73 and 41% of mtDNA variation was recorded among and within sampling sites, respectively, compared to 97 and 3% of total variation, respectively, found between populations inhabiting older islands (Beheregaray et al. 2004).
Nevertheless, we revealed a significant level of genetic differentiation among sampling sites with both microsatellite and mitochondrial markers. This may be seen as characteristic of island populations where vicariant events determined by eruptive regimes of shield volcanoes can create a patchy distribution of suitable habitats, thereby promoting genetic divergence (e.g.,Vandergast et al. 2004). Genetic distinction was particularly evident for tortoises of La Cazuela and, to a minor extent, Roca Union, the two easternmost sampling sites of Sierra Negra. Similarly, relatively high differentiation was recorded between the eastern and western slopes of Cerro Azul, a result in agreement with previous, preliminary microsatellite analysis (Ciofi et al. 2002). These patterns of genetic distinction were most probably affected by magmatic activity of the past 100,000 years (Naumann and Geist 2000). Sierra Negra, for instance, has the highest eruptive rate of any volcano in the archipelago and has undergone over 90% resurfacing in the past 4500 years (Reynolds et al. 1995). In particular, the oldest superficial lava flows are found on the southern and eastern flanks of the caldera. Volcanic activities may have represented geographical barriers that allowed occasional main dispersal events as suitable intervening habitats were formed, but their recurrence may have reduced frequent gene flow. The absence of substantial past and recent tortoise movement across sites was inferred by both mtDNA analysis and a Bayesian method for estimating recent migration rates from microsatellite allele frequencies. The intensive harvests of tortoises and the conversion of natural habitats into pastures on Sierra Negra may have also affected, during the past 150 years, the isolation of La Cazuela from the southwestern range of the species, as seen from the microsatellite data.
A similar pattern may account for the levels of genetic divergence recorded across Cerro Azul. Young and intermediate basalts are <3000 years old and cover most of the lower and upper flanks of the volcano. The southeastern sector of Cerro Azul, on the other hand, which is characterized by basalts up to 5000 years old, is covered by relatively uniform vegetation. Sampling locations relative to this area are Cinco Cerros and Cabo Rosa, where the lowest levels of genetic differentiation and relatively high migration rates were recorded between sample sites. This is in agreement with previous studies where Bayesian clustering methods grouped a significant proportion of tortoise genotypes from Cabo Rosa with eastern Cerro Azul (Ciofi et al. 2002), while Roca Union had >50% of its genotypes clustered with La Cazuela. Here, a more comprehensive analysis based on different molecular markers and a larger sample set and size highlight differences among locations and limit genetic similarities to the very southern boundary between volcanic areas.
Genetic differentiation at mitochondrial and nuclear markers:
The level of population structure across southern Isabela volcanoes was higher for mtDNA control region (ΦST = 0.32) than for microsatellites (θ = 0.051; ρST = 0.069). Such differences were also clear from the pattern of molecular variance distribution among and within sampling sites. Several factors linked to differential gene flow, mutation rate, and genetic drift may account for the observed discrepancy. A first possibility is that mutation rates at microsatellites may be too high in relation to the timing of population divergence to obtain reliable estimates. This would result in size homoplasy of alleles (Estoup et al. 1995; Garza and Freimer 1996; Grimaldi and Crouau-Roy 1997; Orti et al. 1997) and a reduction of intraspecific diversity (Paetkau et al. 1997; Lu et al. 2001; Roberts et al. 2004; but see Richard and Thorpe 2001). Also, measures of differentiation between groups as measured by F-statistics may be reduced when high levels of within-population heterozygosities are recorded (Charlesworth 1998). This is often the case for highly variable loci such as microsatellites (Hedrick 1999) and, given the average values of heterozygosity per sample size for Galápagos tortoises of 0.70, may be a plausible explanation for the differences between mitochondrial and nuclear DNA estimates.
The relatively higher structuring of the mtDNA control region could also be accounted for by a stronger effect of demographic stochasticity (i.e., genetic drift) for mtDNA than microsatellite loci (Hutchison et al. 1974; Birky et al. 1989; Avise 2004). Short coalescence time of mtDNA genes (Avise et al. 1984) can be exacerbated by a reduction of mitochondrial effective size due to sex-biased founding events, multiple insemination of founding females (Beheregaray et al. 2003), or a relatively small population size, as in the case of tortoises from southern Isabela (Pritchard 1996).
Given mtDNA inheritance from females to offspring, the higher genetic structuring among sampling sites recovered by mtDNA could also be due to a higher vagility of male tortoises. A pattern of sex-biased dispersal, as recorded for other Galápagos reptiles (Rassmann et al. 1997), was not found for island populations of giant tortoises. However, although no significant difference in nuclear DNA was found between sampling sites for males and females, data from mark–recapture studies on other islands have so far recovered significant sexual variability in spatial ecology (H. Snell, unpublished data), an indication that different, gender-biased movements may well be part of the species' life history.
Patterns of genetic and morphologic distinctiveness across southern Isabela:
In a volcanic environment, different topologies of xeric and mesic habitats may succeed one another fostering the radiation and successive divergence of individuals into a variety of empty niches (faunal drift, sensu Barton 1989). Morphological variants of Galápagos dome- and saddle-back giant tortoises are a vivid example of such ecological dynamics (Pritchard 1996; Caccone et al. 2002). On Isabela, differently domed tortoises have been described in the northern and central volcanoes with fairly distinct genotypic structure (Caccone et al. 2002; Ciofi et al. 2002), and a similar pattern, based on morphological information, has been suggested across the southern range of the species. Tortoises from Sierra Negra, associated to the named taxon guntheri, have larger males with a flatter dorsal area of the carapace. The carapace is also less elevated than other dome tortoises and has small marginal and smooth, black, nonstriated scutes. The named taxon vicina, described on Cerro Azul, is characterized by a more domed and striated carapace with brown scutes and larger marginals (Pritchard 1996). This classification, however, is still subject to debate. While a multivariate analysis of morphological characters showed some degree of distinction between tortoises from Sierra Negra and Cerro Azul (Fritts 1984), the number of specimens for which clear distinctive characteristics are present may be too low to justify nomenclatural recognition (Pritchard 1996; but see Ernst et al. 2000).
In our study, we investigated whether patterns of genetic divergence provide evidence of population structuring similar to that proposed by previous morphological analysis. Apart from the relatively high degree of genetic divergence shown by tortoises from La Cazuela, we could not recover significant differences in the pattern of genetic differentiation across sampling sites in either the mitochondrial or nuclear genome that may warrant a clear taxonomic subdivision between and within volcano-associated populations. Measures of genetic divergence across sampling sites located within the boundaries of either of the two volcanic areas were not significantly different from values recorded between locations from different volcanoes.
Fritts (1983) and Pritchard (1996) further suggested that body size and details of carapace morphologies represent adaptive responses to different habitat conditions across Cerro Azul and Sierra Negra. Often, genetic similarities are related to habitat affinities rather than to main geographical barriers (Malhotra and Thorpe 1997; Lu et al. 2001; Sacks et al. 2004) and clinal variation at nuclear loci can result from altitudinal changes in fitness traits among genotypes (Storz and Dubach 2004). Tortoises from east and west Cerro Azul are indeed distributed along an altitudinal gradient with different habitat types. However, although genetic distinction was recorded between the two sides of the caldera, there were no significant differences across tortoises from different altitudes and ecotypes (e.g., Las Tablas and Las Pampas, across the western slope of Cerro Azul; Ciofi et al. 2002 and this study). Similarly, tortoises sampled at low altitudes across Sierra Negra and Cerro Azul had significantly different population genetic structure despite the fact that all sampling locations were characterized by similar habitat type, namely dry and sparse low canopy forest. Overall patterns of genetic diversity across southern Isabela appear therefore not to mirror either geographical or ecological distribution of named taxa. Similar results have been recovered for other Galápagos reptiles such as marine iguanas, for which morphological and coloration characters were thought to be the result of phenotypic plasticity (Rassmann et al. 1997).
Circumstantial evidence indicated the occurrence of both guntheri and vicina in Cinco Cerros, a southeastern Cerro Azul location characterized by a more mesic vegetation than other lowland localities (Porter 1976). However, the migration rate of tortoises from the eastern side of Cerro Azul to Cinco Cerros was only ∼8% (Table 6). Allele frequency distributions at any locus did not show a pattern that might suggest the presence of a bimodal hybrid zone, with individuals genetically similar to either of two putative parental genotypes (Jiggins and Mallet 2000). Moreover, the presence of 2 private alleles and 10 unique sequences out of 17 haplotypes did not suggest the presence of a putative hybrid zone in Cinco Cerros. On the other hand, the nearby Cabo Rosa, where tortoises were reported of a single taxon (guntheri), had a migration rate of 30% from Cinco Cerros. This occurred despite Cinco Cerros and Cabo Rosa sharing only 1 mtDNA haplotype out of 22 identified from these two sampling locations.
The observed pattern of genetic divergence and similarities did not support suggested morphological subdivision between southern locations of Cerro Azul and Sierra Negra. However, asymmetrical gene flow has been reported to occur where peripheral populations meet at range boundaries (Cicero 2004; Grant et al. 2004) and can be maintained without affecting selective forces on adaptive quantitative traits, including phenotypic variation (Sandoval 1994; Ross and Keller 1995; Hedrick 1999 ; Haavie et al. 2000; Andersson et al. 2004; Jordan et al. 2005). In addition, although a correlation between patterns of diversity at neutral loci and differences in life-history traits has been reported (Stephan and Cho 1994; Kashi et al. 1997; Li et al. 2000; Turpeinen et al. 2001), empirical investigations on neutral markers as indicators of genetic variation relevant to evolutionary potentials and population persistence are still scarce (Pearman 2001). Further studies employing variable, functional genes (Ford 2002) and a comprehensive set of phenotypic data (Hoekstra et al. 2004) may recover additional information on genotypic and morphological diversity in southern Isabela and therefore help direct conservation efforts.
Implications for population management:
In recent years, population management planning has been shifting from a discipline-bound decision-making policy to a more integrative approach, where ecological and genetic studies are combined in an effort to identify intervention priorities for conservation of threatened wildlife (Paetkau 1999; Roman et al. 1999; Crandall et al. 2000). When considering genetic diversity, the conservation of differential adaptive potential promotes the retention of population viability traits in a changing environment, while the recognition of phylogeographic distinctiveness could secure the preservation of evolutionary lineages (Moritz 2002). It has been argued that distinct evolutionary lineages, if lost, are more difficult to recover than functional diversity (Moritz 1999). Moreover, phenotypic variation can arise through the interaction between genome and environment (Lynch and Walsh 1998) as seems to be the case with tortoises from southern Isabela. It is therefore important to recognize evolutionary processes as much as ecological and phenotypical traits to conserve populations that represent both adaptive and historical diversity (Ferrier 2002; Moritz 2002). From this perspective, tortoises from different locations of Cerro Azul and Sierra Negra warrant protection efforts given their degree of genetic diversity and evolutionary history regardless of supposed degree of morphological separation.
It is well known that small, isolated populations face a higher probability of extinction than those in a metapopulation system, where migration between groups is possible (Hanski and Gilpin 1997; Frankham 1998). A decrease in gene flow exacerbates genetic drift, a process suggested by the relatively low genetic variability found at La Cazuela, and edge effects linked to human-related threats and stochastic events. Restoration should be attempted more often for populations, such as La Cazuela, that have become isolated as a result of both vicariance and/or recent anthropogenic activities (e.g., Powell and Gibbs 1995; Naumann and Geist 2000; Kaiser 2001). During the past 3 decades, Galápagos tortoises have been subject to an intensive captive breeding program (Merlen 1999), which includes a proposed demographic reinforcement of southern Isabela. It is difficult to know to what extent adaptation or coadaptation develops among different environmental gradients and constrains population viability through outbreeding (Lynch and Walsh 1998; Moritz 1999). In any case, a clear differentiation between volcanoes or altitudinal ranges based on morphological grounds is difficult to support for these tortoises. Each population of southern Isabela has a distinct genetic structure and, although some extreme circumstances may necessitate mixing of individuals from genetically different stocks (Hedrick 1995), integrative analysis of nuclear and mitochondrial DNA suggests that assignment tests should be carried out within a captive breeding framework prior to devising augmentation plans for southern Isabela.
We thank the Parque Nacional de Galápagos (PNG) for providing permission to conduct field research and the staff of the Charles Darwin Research Station and PNG for logistical support during sample collection. Michel Milinkovitch is gratefully acknowledged for providing primer sequences for microsatellite loci PCR amplification. We also thank Thomas H. Fritts and Jonathon Marshall for insightful comments on previous versions of the manuscript. The project was funded by the U.S. National Science Foundation (DEB 9322672), National Geographic Society (6800-00), and the Yale Institute for Biospherics Studies (YIBS). CITES permits were issued by the United States Fish and Wildlife Service and Instituto Ecuatoriano Forestal de Areas Naturales y Vida Silvestre. C. Ciofi and L. B. Beheregaray were supported by a Gaylord Donnelley Environmental Fellowships from YIBS at the time of this study.
Communicating editor: S. W. Schaeffer
- Received July 8, 2005.
- Accepted December 13, 2005.
- Copyright © 2006 by the Genetics Society of America