Testing phylogeographic predictions on an active volcanic island: Brachyderes rugatus (Coleoptera: Curculionidae) on La Palma (Canary Islands)

Molecular Ecology (2006) 15, 449–458 doi: 10.1111/j.1365-294X.2005.02786.x Testing phylogeographic predictions on an active volcanic island: Brachyderes rugatus (Coleoptera: Curculionidae) on La Palma (Canary Islands) mBlackwell Publishing Ltd B R E N T C . E M E R S O N ,* S H A U N F O R G I E ,*‡ S A R A G O O D A C R E * and P E D R O O R O M Í † Centre for Ecology, Evolution and Conservation, School of Biological Sciences, University of East Anglia, Norwich NR4 7TJ, UK, †Departamento de Zoología, Facultad de Biología, Universidad de La Laguna, 38205 La Laguna, Tenerife, Spain Abstract Volcanic islands with well-characterized geological histories can provide ideal templates for generating and testing phylogeographic predictions. Many studies have sought to utilize these to investigate patterns of colonization and speciation within groups of closely related species across a number of islands. Here we focus attention within a single volcanic island with a well-characterized geological history to develop and test phylogeographic predictions. We develop phylogeographic predictions within the island of La Palma of the Canary Islands and test these using 69 haplotypes from 570 base pairs of mitochondrial DNA cytochrome oxidase II sequence data for 138 individuals of Brachyderes rugatus rugatus, a local endemic subspecies of curculionid beetle occurring throughout the island in the forests of Pinus canariensis. Although geological data do provide some explanatory power for the phylogeographic patterns found, our network-based analyses reveal a more complicated phylogeographic history than initial predictions generated from data on the geological history of the island. Reciprocal illumination of geological and phylogeographic history is also demonstrated with previous geological speculation gaining phylogeographic corroboration from our analyses. Keywords: colonization, geology, mitochondrial DNA, nested clade phylogeographic analysis, phylogeography Received 16 March 2005; revision received 5 August 2005; accepted 23 September 2005 Introduction Molecular phylogenetic methods have increasingly been turned to as tools for interpreting and understanding the origins of species diversity on islands (Emerson 2002). Oceanic island systems are attractive environments for studying evolution for a number of reasons: (i) they present discrete geographic entities within defined oceanic boundaries, (ii) gene flow between individual islands is reduced by oceanic barriers, (iii) their often small geographic size has made the cataloguing of flora and fauna easier than continental systems, (iv) despite their small geographic size they can contain a diversity of habitats, and (v) they are often geologically dynamic with historical and contemporary volcanic and erosional activity. With few exceptions (e.g. Correspondence: Brent Emerson, Fax: 44-01603-592250; E-mail: b.emerson@uea.ac.uk ‡Present address: HortResearch, Mt Albert Research Centre, Private Bag 92169, Auckland, New Zealand © 2006 Blackwell Publishing Ltd Juan et al. 1998; Pestano & Brown 1999; Brown et al. 2000; Holland & Hadfield 2002; Beheregaray et al. 2003) the majority of molecular phylogenetic studies have focused on the relationships among multiple species occurring on one or more islands, or one species occurring on multiple islands. However, the above features that make oceanic islands attractive for this purpose also provide an ideal template for the study of intraspecific phylogeographic patterns within individual islands. Continental areas are the typical domain for intraspecific phylogeographic studies, primarily to infer the historical movements of species, and to relate these to past climatic change events mediated by glacial cycling. Volcanically active oceanic islands, particularly geologically young ones, offer the opportunity to generate phylogeographic predictions based on documented geological events of volcanic growth and erosional decay. In this study, we use the geologically young and well-characterized island of La Palma to generate phylogeographic predictions for Brachyderes rugatus rugatus, a local endemic subspecies of curculionid beetle occurring throughout the island in the 450 B . C . E M E R S O N E T A L . Fig. 1 Map of La Palma showing sampling locations, the distribution of Pinus canariensis (darker shading) and major geological features. The geologically older terrains of the northern shield occur above the dashed line. The geologically younger terrains of Cumbre Vieja occur below the dotted line. The solid line encloses the Bejenado terrain. Large-scale erosion or landslide events are indicated by toothed lines. Inset shows La Palma in relation to the other six main islands of the Canary archipelago. See text for further details. forests of Pinus canariensis. We then test these predictions with sequence data for the mitochondrial DNA (mtDNA) cytochrome oxidae II (COII) gene. Studies of the geology of La Palma have led to a fairly complete understanding of the island’s geological history (Ancochea et al. 1994; Carracedo & Day 2002) with two welldefined edifices — the northern shield occurring centrally in the north, composed mainly of older volcanic terrains, and a southern ridge, the Cumbre Vieja constituted mainly by terrains of more recent volcanic origin (Fig. 1). Subaerial development of the northern shield began about 1.7–2.0 million years ago (Ma), and this development continued until about 0.55 Ma, with limited activity occurring until about 0.4 Ma (Fig. 1). The subaerial volcanic development of the northern shield was dominated by two volcanoes. The Garafia volcano was active from 1.7 to 1.2 Ma, and following a catastrophic landslide of its southern flank, activity was then dominated by the Taburiente volcano, which remained active from 1.2 to 0.4 Ma. From about 0.8 – 0.7 Ma, in the final stages of the development of the Taburiente volcano, the southward migration of volcanism began through the southern, or Cumbre Nueva, rift zone of the volcano (Fig. 1). At about 0.56 Ma this rift became unstable and its western flank collapsed into the sea (Fig. 1). Following this collapse the Bejenado volcano dominated activity from 0.56 to 0.49 Ma, forming what is now the southeast wall of the Caldera de Taburiente (Fig. 1), with possible minor activity continuing until about 0.2 Ma. Following the Bejenado formation erosion enlarged the Caldera de Taburiente (Fig. 1). After the end of the growth of the Bejenado volcano, the entire northern shield became volcanically quiescent, and the island may have entered a period of volcanic calm. © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 P H Y L O G E O G R A P H Y O F B R A C H Y D E R E S R U G A T U S 451 However, it is also surmised that volcanism in the south of La Palma may have built the submarine and core parts of the Cumbre Vieja that are now concealed by younger lavas of the Cumbre Vieja eruptive period (Carracedo & Day 2002). It was from 0.12 Ma to the present that intense volcanic activity in the south of La Palma led to the rapid development of the Cumbre Vieja ridge. This greatly increased the size of the island, particularly to the south, and the southern half of La Palma is dominated by these recent lavas (Fig. 1). It is a reasonable assumption that a species with limited dispersal ability, particularly a widespread species, present on La Palma during its volcanic and erosional history, would have been influenced by these events. That being the case, these events should be evidenced through genetic footprints of range expansions and range fragmentations, consistent with predictions from volcanic and erosional events across the 690-km2 area of La Palma. In this respect the flightless pine weevil, B. r. rugatus, an endemic to La Palma, is an ideal species to test phylogeographic predictions from volcanic and erosional events. La Palma is predominantly vegetated by P. canariensis, which extends from the northern shield, through Cumbre Nueva to the south of Cumbre Vieja (Fig. 1). B. r. rugatus is reliant upon P. canariensis (although it can also be found on other introduced pine species), thus the distribution of B. r. rugatus on La Palma is clearly defined. Substantial genetic variation, with a nonrandom geographic distribution, has previously been reported for B. r. rugatus on La Palma, with a colonization origin from Tenerife early in the emergence of the island’s subaerial terrains (Emerson et al. 2000). This earlier study used traditional molecular phylogenetic methods of maximum likelihood and neighbour joining to infer the colonization history of the four subspecies of B. rugatus across the archipelago. No obvious geological explanation for the geographic structuring of mitotypes was found, but the sampling strategy and analytical techniques were not designed to specifically address this issue. Here we sample 138 beetles from 18 localities across the distribution of B. r. rugatus for the construction of a haplotype network and application of nested clade phylogeographic analysis (NCPA) for 570 base pairs (bp) of sequence data for the mtDNA COII gene. Based on the volcanic and erosional history of La Palma, and the species biology of B. r. rugatus, we make several predictions regarding the geographic structuring and demographic signature of the genetic data. We predict that ancestral haplotypes should be predominantly located in the northern shield, with the newer Bejenado and Cumbre Vieja volcanic terrains possessing more derived haplotypes, located peripherally in a haplotype network. We also predict signatures of population expansion in the areas of the geologically young Cumbre Vieja volcanic terrains: haplotypes occurring in multiple sampling localities, and haplotypes with multiple derivatives differing by only one or a few mutations. © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 Materials and methods Sampling Brachyderes were collected from La Palma between 1996 and 2003. Eighteen localities were chosen, 11 in the northern shield, 2 in the Bejenado volcanic terrain, and 5 in the Cumbre Vieja volcanic terrain (Fig. 1). We aimed to collect eight individuals for each location, but were unsuccessful at Casas del Monte where only two beetles were found. All samples were stored in absolute ethanol and stored at 4 °C prior to extraction of DNA. DNA extraction, PCR amplification and DNA sequencing For each individual beetle, DNA was extracted from the head and pronotum using the QIAGEN DNeasy extraction kit (QIAGEN) following manufacturers instructions. Polymerase chain reaction amplification of a 672-bp fragment of the mtDNA COII gene was carried out in a PerkinElmer ABI 9700 thermocycler with reagents and conditions as described in Emerson et al. (2000) with the exception that a 3 mm concentration of MgCl2 was used. Reactions were then cleaned using a QIAquick PCR clean-up kit (QIAGEN) following manufacturers instructions. Sequencing reactions were performed with the forward PCR primer and an additional internal primer using the PerkinElmer BigDye terminator reaction mix with the PCR amplification primers and run on a PerkinElmer ABI 3700 automated sequencer. Data analyses DNA sequences were aligned by eye and a haplotype network constructed using the parsimony criterion as described in Templeton et al. (1992) and implemented in the program tcs version 1.18 (Clement et al. 2000). Ambiguities in the haplotype network were resolved following the two criteria suggested by Crandall & Templeton (1993): (i) within a cladogram rare haplotypes are more likely to be tip haplotypes, and common haplotypes are more likely to be interior haplotypes; and (ii) singleton haplotypes are more likely to be connected to haplotypes from the same population as opposed to haplotypes from different populations. The haplotype network was then manually converted into a nested series of clades using the rules defined in Templeton et al. (1987), and Templeton & Sing (1993) for a nested clade phylogeographic analysis (NCPA) (Templeton et al. 1995; Templeton 1998, 2004) using the program geodis version 2.2 (Posada et al. 2000) and the most recent inference key provided on the geodis webpage (http:// darwin.uvigo.es/software/geodis.html). Coalescent theory predicts that clades at the tips of a tree are younger than interior clades to which the tips are connected. Thus 452 B . C . E M E R S O N E T A L . Locality Latitude Longitude Number of haplotypes Above Fuente de Olén Montana Tagoja Fuente de Olén El Bejenado Aridane Montana de la Venta El Jable Lomo María Fuencaliente Altos de Jedey Casas del Monte Lomo Carballo Garafia Lomo Machín Pinar Garafia Below Observatory Puntagorda Tirajafe El Jesus 28.44.58 N 28.43.44 N 28.43.78 N 28.40.44 N 28.41.11 N 28.37.01 N 28.37.00 N 28.34.74 N 28.30.04 N 28.35.71 N 28.46.36 N 28.47.96 N 28.47.82 N 28.47.02 N 28.46.58 N 28.46.93 N 28.44.23 N 28.40.99 N 17.49.79 W 17.47.84 W 17.48.74 W 17.51.04 W 17.54.68 W 17.50.52 W 17.51.60 W 17.52.28 W 17.50.29 W 17.51.05 W 17.48.73 W 17.50.03 W 17.54.17 W 17.54.03 W 17.54.37 W 17.58.82 W 17.58.01 W 17.56.11 W 8 4 7 7 6 1 2 5 5 4 2 4 5 5 5 6 5 2 Table 1 Sampling localities, latitude, longiHaplotype 50, 51, 52, 53, 54, 55, 56, 57 30(4), 47(2), 48, 49 24, 25, 26(2), 27, 28, 29, 30 2(2), 3, 4, 5, 6, 7, 8 2(3), 51, 66, 67, 68, 69 2(8) 1(2), 2(8) 14(2), 15(3), 16, 17, 18 42, 43(2), 44(2), 45(2), 46 15, 34(5), 35(1), 36(1) 9, 13 9(4), 31, 32, 33(2) 9(4), 10, 11, 12, 13 37, 38, 39(3), 40(2), 41 13(2), 39(3), 58, 59, 60 13, 19(2), 20, 21(2), 22, 23 19, 39(3), 61(2), 62, 63 64(7), 65 contrasting interior clades with tip clades represents a temporal contrast of older with younger demographic events. The geographic distribution of haplotypes is also quantified for the NCPA through two distance measures (Templeton et al. 1995). The clade distance (Dc) measures the spatial spread of a clade, while the nested clade distance (Dn) measures how far a clade is from those clades with which it is nested into a higher-level clade. Through the combination of genealogical, geographic and temporal information, NCPA first tests the null hypothesis of random geographic distribution of haplotypes for all nesting levels. If this hypothesis is rejected then the fit of the data to different demographic scenarios incorporating range fragmentation, long-distance dispersal, restricted gene flow, range expansion (Templeton 2004) and secondary contact (Templeton 2001) is tested by comparing the size and significance of Dc and Dn values to predictions from coalescent theory and simulation analysis (Templeton et al. 1995). Results Sequences from 138 individual beetles yielded 69 different haplotypes, with from one to eight haplotypes recorded from a single location (Table 1). Of the 69 haplotypes, only eight were found to occur in more than one location (Fig. 2). All sequences have been deposited in the EMBL nucleotide sequence database under accession numbers (AJ389825–AJ389838 and AM072352–AM072410). Of a total of 570 sites 66 (12%) were variable, among which 41 (62%) were parsimony informative. tude, number of haplotypes in each locality, and haplotype numbers for each locality, with the number of individuals containing that haplotype in brackets if more than one. Haplotype numbers correspond to those in Fig. 4 with a connection limit of 10 mutations and eight loops, four of which were resolved by applying the criteria of Crandall & Templeton (1993). Of the four remaining loops, three could each be reduced to two possible alternatives using the criteria of Crandall & Templeton (1993). Of these, one represented a trivial arrangement and was broken arbitrarily, while both alternatives for the other two were considered for the nested clade analysis (see below). The eighth loop has three equally probable break points connecting three phylogroups that are geographically discrete, with the exception of haplotypes 15, 34, 35 and 36 in Altos de Jedey. The three possible arrangements of the phylogroups are equally probable when considering genealogical information alone, but differ when considering both the genealogy and geography of the three networks (Fig. 3). Networks B and C (Fig. 3) juxtapose phylogroups 1 and 3 by only two mutational differences implying a biotic connection between the geographically disparate areas of 1 and 3 over a short temporal period, at the exclusion of the geographically intermediate area 2. Alternatively, network A juxtaposes phylogroups 1 and 3 with phylogroup 2, a biologically realistic arrangement that does not imply longdistance movement of Brachyderes rugatus rugatus through, but not including, an already inhabited area. Figure 4 shows the haplotype network of Fig. 3A and its nesting design inferred following the rules of Templeton et al. (1987) and Templeton & Sing (1993). Nesting designs were also inferred for the alternative topologies of Fig. 3A and NCPA was performed on these arrangements as well. Nested clade phylogeographic analysis Haplotype network and nested design Network estimation with tcs resulted in a single network Table 2 shows the nested contingency analysis of geographic associations and the interpretations of the statistically © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 P H Y L O G E O G R A P H Y O F B R A C H Y D E R E S R U G A T U S 453 Fig. 2 The geographic distributions of haplotypes occurring in more than one location. White circles represent single haplotypes found in more than one location. Black circles indicate haplotypes derived from a widespread haplotype, with each black circle representing a unique haplotype. Numbers inside black circles indicate when more than one unique derived haplotype occurs in a location. Grey filled circles indicate the geographic locations of ancestral haplotypes to the widespread haplotype when these are known. Table 2 Significant clades for nested contingency analysis of geographic associations for mitochondrial DNA COII nucleotide sequence haplotypes for Brachyderes rugatus rugatus on the island of La Palma and interpretation of the geodis output using the inference key of Templeton (2004) Clade 1.1 2.1 2.3 2.14 3.2 3.6 3.7 4.1 4.2 4.3 Total cladogram Chi-squared statistic 45.72 31.54 6.00 10.00 7.20 13.00 30.00 44.14 41.46 69.16 226.51 Probability Chain of inference Inference 0.0042* 0.0532** 0.0473* 0.0219* 0.0527** 0.0010* 0.0000* 0.0000* 0.0000* 0.0000* 0.0000* 1-2-11-12 No 1-2-11-17 No 1–2 1-19-20-2 1-2-11-12 No 1–19 No 1-2-3-4 No 1-2-11-12-13 Yes 1-2-11-12 No 1-2-3-5-6-13 Yes 1-2-11-12 No contiguous range expansion inconclusive outcome inconclusive outcome inconclusive outcome contiguous range expansion allopatric fragmentation restricted gene flow with isolation by distance past fragmentation followed by range expansion contiguous range expansion past fragmentation followed by range expansion contiguous range expansion *Significant at the 0.05 level, **significant at the 0.1 level. significant clades using the updated key of Templeton (2004) for Fig. 4. At the level of the entire cladogram a history of contiguous range expansion offers the best explanation for the geographic distribution of genetic variation within B. r. rugatus. However, among the three major clades other historical processes also feature. Clade 4.1 is best explained by past fragmentation followed by range expansion, with internal clades 1.1 and 3.2 both conforming to a history of contiguous range expansion. Contiguous range expansion © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 also has the best explanatory power for clade 4.2, with none of its constituent subclades having any significant geographic structure. Overall clade 4.3 conforms to a history of past fragmentation followed by range expansion, but internally subclade 3.6 is consistent with a history of allopatric fragmentation, while the geographic distribution of subclade 3.7 is best explained by restricted gene flow with isolation by distance. Nested clade phylogeographic analysis of the possible alternatives to the Fig. 3A were also 454 B . C . E M E R S O N E T A L . Fig. 3 Three equally probable haplotype networks constructed with tcs version 1.18 (Clement et al. 2000) for 570 bp of mitochondrial DNA COII nucleotide sequence data for Brachyderes rugatus rugatus. The three arrangements, A, B and C, are derived from the resolution of a single loop. All are equally probable when only considering genealogical relationships, but network A is the single most probable network when considering the geographic placement of the three main phylogroups, boxed and numbered 1–3 (see text for details). The map indicates the geographic distribution of the phylogroups. Dashed lines ‘a’ and ‘b’ indicate the alternative connections of two additional loops, with breakpoints indicated by dotted lines. These alternatives were considered for the NCPA. Open circles indicate sampled haplotypes, and closed circles indicate unsampled or extinct haplotypes. performed. The alternative connection ‘a’ between subclades 2 and 3 of Fig. 3A yielded a similar nesting design to Fig. 4 and the same set of statistical inferences. The alternative connection ‘b’ within subclade 2 of Fig. 3B resulted in changes to the nesting design of clades 4.1 and 4.2 of Fig. 4 with subclades 1.11 and 1.25 forming part of clade 4.1. Under this arrangement clade 4.1 is now consistent with contiguous range expansion, with two subclades also reflecting this demographic history, and clade 4.2 has an inconclusive outcome. The three phylogroups of Fig. 3 are geographically well defined and nonoverlapping except for the sampling locality of Altos de Jedey where haplotypes from both phylogroups 1 and 2 can be found. One possible explanation for this would be if range expansion has occurred from Altos de Jedey in two directions, to the north (phylogroup 2) and to the south (phylogroup 1). If this were the case then the haplotypes from these two phylogroups in Altos de Jedey should be closely related and ancestral in the network. However, phylogenetically the haplotypes from these two clades are very divergent from each other. The three haplotypes from clade 4.2 (15, 34, 36) are between 18 and 22 mutations different from haplotype 35 of clade 4.1. An alternative explanation is that Altos de Jedey is a point of secondary contact between phylogroups 1 and 2. Templeton (2001) describes a strategy for quantifying secondary contact between divergent lineages within a nested clade framework. When secondary contact occurs between previously fragmented populations, haplotypes or clades with very divergent geographic centres will be represented at the same location. The strategy of Templeton (2001) involves calculating the average distance from the geographic centre of the haplotypes or clades at increasing nesting levels for the population of interest. In panmictic populations, all haplotypes and clades are expected to have the same geographic centre, so the average distance © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 P H Y L O G E O G R A P H Y O F B R A C H Y D E R E S R U G A T U S 455 Fig. 4 Nesting design for the haplotype network constructed from 570 bp of mitochondrial DNA COII nucleotide sequence data for Brachyderes rugatus rugatus. Numbered circles refer to specific haplotypes (Table 1), and smaller filled circles represent unsampled or extinct haplotypes. Grey filled circled denote haplotypes occurring in the three areas indicated on the map, collectively referred to as Olén. between the geographic centres of clades should be the same for all populations. Under a scenario of isolation by distance, lower clade levels are expected to have small positive distances between the geographic centres of the clades that approach zero at higher clade levels. However, under a scenario of secondary contact the average distance is expected to either remain high or rise with clade level, until a maximum should be reached at the clade level where the fragmentation was inferred. Altos de Jedey does indeed conform to a zone of secondary contact, exhibiting an increasing average distance between geographic centres proceeding from lower to higher nesting levels (NL): NL 0 = 1.1 km, NL 1 = 3.6 km, NL 2 = 5.7 km, NL 3 = 7.7 km, NL 4 = 9.4 km. Discussion Ancestral haplotypes and ancestral areas The phylogeographic history of Brachyderes rugatus rugatus on La Palma appears more complex than a simple extrapolation of the geological history of La Palma would suggest. Our first prediction was that ancestral haplotypes should predominantly occur in the northern shield, with derived haplotypes featuring more on the El Bejenado and Cumbre Vieja terrains. Coalescent theory predicts that © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 ancestral haplotypes will occur at high frequency, be represented in the greatest number of populations, have multiple connections to low frequency haplotypes, and be located at the interior of a network (Crandall & Templeton 1993; Posada & Crandall 2001). None of the haplotypes within our network conform to this definition. A number of haplotypes occur at a frequency higher than one, but only eight of these occur in more than one population (Fig. 2), with only four of those having multiple connections to other haplotypes, and none of these can be considered truly central within the network (Fig. 4). These four haplotypes, 9, 39, 15 and 2, can only be considered ancestral on a local scale — ancestral sequences of recent regional population expansions. The lack of an obvious ancestral sequence (or sequences) under the above criteria can be reconciled with a population history involving geographic structure checkered by regional population extinction and recolonization. Although NCPA cannot explicitly test this hypothesis, the conclusions of fragmentation, range expansion, secondary contact, and two additional features of the network in Fig. 4 support such a scenario. First, the network is composed of three geographically distinct phylogroups (Fig. 3), indicating some historical disjunction between these areas. Second, the network itself is characterized by many missing internal haplotypes (70%), and these are not randomly distributed 456 B . C . E M E R S O N E T A L . with only 37% of tip haplotypes (level 0 clades) being connected to internal unsampled sequences. When looking at level 3 clades, the percentage of missing intermediates ranges from 40% (clade 3.7) to 85% (clade 3.4). Furthermore, grouping these level 3 clades into tip clades (3.1, 3.5, 3.7, 3.8) and internal clades (3.2, 3.3, 3.4, 3.6) reveals that internal clades have more missing intermediates (72%) than tip clades (59%). This indicates that rather than occurring at higher frequency, the oldest of haplotypes are likely to be rare or extinct. Although the frequencies of ancestral haplotypes can be expected to decrease through time under a scenario of population structure with regional extinction and recolonization, their ancestry will still be identified by their internal placement within a network (Posada & Crandall 2001). Although the exact root location of the network of Fig. 4 cannot be determined, clade 3.3 approximates the centre of the network. Clade 3.3 contains 2 internal haplotypes, 53 from above Fuente Olén, and 49 from the geographically proximate Montaña Tagoja (Fig. 1). Three other tip haplotypes within this clade (29, 48, 52) are also restricted to either of these populations or the geographically intermediate Fuente Olén. Looking at the phylogenetic relationships of the remaining haplotypes from these three locations reveals that, with the single exception of haplotype 27, all form an inclusive set of connections around clade 3.3 (Fig. 3). The genealogical and geographic unity of these 17 haplotypes, combined with their interior location within the network, suggests an ancestral area in the region of the sampling locations Montaña Tagoja, Fuente Olén and above Fuente Olén (from hereon referred to as Olén). Thus, consistent with the greater antiquity of the northern shield, ancestral haplotypes are found to predominate here, but in a restricted area of the eastern flank of the shield. Ancestral areas and range expansions As part of our first hypothesis, we also predicted that the newer Bejenado and Cumbre Vieja volcanic terrains would possess more derived haplotypes, located terminally in a haplotype network, than the northern shield terrain. Additionally we hypothesized that signatures of population expansion should feature more in the newer southern terrains than the older northern terrains. This is clearly not the case, but can be reconciled with a more complex phylogeographic history of B. r. rugatus on La Palma than we predicted. The presence of three distinct phylogroups (Fig. 3) is clearly at odds with a scenario of a long-term residency in the northern shield followed by a more recent expansion into the southern terrains. Rather, the ancestral area of Olén appears to have been the source of three distinct range expansion events corresponding to the three level 4 clades of Fig. 4, one into the northern shield (clade 4.3) and two into the Cumbre Vieja terrain (clades 4.1 and 4.2). Clade 4.3 is consistent with a history of past fragmentation followed by range expansion, and several features suggest this expansion has occurred radially around the Caldera de Taburiente, clockwise from Olén, to give rise to phylogroup 3. A pattern of radial colonization through the forest crowning La Caldera has previously been reported for Brachyderes rugatus calvus on the Island of Gran Canaria (Emerson et al. 2000). Clade 4.3 connects to ancestral haplotypes from Olén, and its most interior clade, 2.11 includes one of these haplotypes that also occurs in Aridane (Fig. 2A). The three most closely related descendant haplotypes of clade 2.11 are restricted to the western localities of Puntagorda (21, 22) and Jesus (64) suggesting these areas were colonized from Olén through Aridane. Figure 2(B–E) shows haplotypes occurring in more than one population in the northern shield, as well as any sequences that have descended from these, and where it is known, the geographic location of the sequence immediately ancestral in the network. In the three instances where the ancestral location can be determined it is to the extreme west of the shield, indicating subsequent colonization toward the east across the top of the northern shield. Sequence divergence rate estimates for invertebrate mtDNA (De Salle et al. 1987; Brower 1994) approximate to 2.15% of sequence divergence per million years (Myr) (1.08% of mutational change per lineage per Myr), and a study of the colonization sequence of the four Brachyderes subspecies on the Canary Islands found their mtDNA mutation rate to be in accord with this (Emerson et al. 2000). The most recent common ancestral (MRCA) sequence of phylogroup 3 occurs in clade 1.34 (Fig. 4). This northern shield lineage is represented by 26 haplotypes with divergences from the ancestral haplotype ranging from one to eight mutations. To estimate a date for the MRCA of phylogroup 3 one should consider all extant haplotypes descending from the MRCA, as each one represents the time interval between the present and the MRCA. Among the 26 haplotypes descending from the MRCA, the average divergence from the MRCA is 0.72%, which equates to 0.66 Ma with a standard error of 0.24 Ma. Clade 4.2 is consistent with contiguous range expansion and unites phylogroup 1 (Fig. 3) with the ancestral area of Olén. Although the NCPA did not detect fragmentation, the large number of missing haplotypes connecting the haplotypes of phylogroup 1 to those from Olén (Fig. 4), combined with the presence of phylogroup 2 in the geographically intermediate area, suggests clade 4.2 was once continuously distributed but suffered extinction in the geographically intermediate area now occupied by phylogroup 2. The phylogeographic pattern and genetic divergences of clade 4.2 are consistent with the speculation of Carracedo & Day (2002) that the young Cumbre Vieja volcanic terrains of 0.12 Ma may possibly cover much older terrains. The most recent common ancestral sequence of © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 P H Y L O G E O G R A P H Y O F B R A C H Y D E R E S R U G A T U S 457 phylogroup 1 occurs in clade 1.15 (Fig. 4), and among the 12 haplotypes descending from the MRCA, the average divergence from the MRCA is 0.75%, which equates to 0.69 (± 0.2) Ma. This estimate is consistent with the existence of B. r. rugatus, and thus forests of Pinus canariensis, in the southern half of La Palma prior to the Cumbre Vieja eruptive period, supporting the conjecture of Carracedo & Day (2002) of the existence of older, but now buried, terrains in the south. Further to this, the distribution of genotypes in this southern region (Fig. 2F) suggests that the area of Lomo María is the likely point of origin of recent range expansion in this southern area. It seems likely that the recent and extensive volcanic and erosional events of the Cumbre Vieja terrain (Fig. 1) would have contributed to the isolation and restriction of a population in the southern part of the island. Clade 4.1 is consistent with a history of past fragmentation followed by range expansion, and this expansion would appear to be from Olén, through Aridane and Bejenado, and then south. Clade 4.1 connects to ancestral haplotypes from Olén, and its most interior clade, 3.2, is also exclusively composed of 10 haplotypes from this area. The descendant clade 2.2 juxtaposing clade 3.2 is composed of haplotypes occurring only in Aridane and El Bejenado, while the tip clade 2.1 includes haplotypes from these two locations as well as all haplotypes from the more southern locations of Montaña de la Venta, Altos de Jedey and El Jable. Figure 2G supports the range expansion into these southern populations to have occurred out from the area of Aridane or El Bejenado. The time of establishment of B. r. rugatus in the Bejenado terrain (Fig. 1) is consistent with this being at or near the end of the main Bejenado eruptive period from 0.56 to 0.49 Ma (Carracedo & Day 2002). The MRCA of the Aridane and El Bejenado haplotypes in clade 3.1 is the interior unsampled haplotype of clade 1.4, and among the 13 haplotypes descending from the mrca, the average divergence from the MRCA is 0.47%, which equates to 0.43 (± 0.14) Ma. It would appear that the timing of the establishment of B. r. rugatus from phylogroup 2 within the Cumbre Vieja terrain (Fig. 1) is more recent than that for phylogroup 1, with only three closely related haplotypes (1, 2 and 35) in these areas of Montaña de la Venta, Altos de Jedey and El Jable. expansions have met, and this is supported by the NCPA. From the phylogeographic pattern it is likely that a further area of secondary contact exists in the northern shield between phylogroups 2 and 3. From the clockwise radial range expansion of clade 4.3 we predict that further sampling between Olén and Casas del Monte (Fig. 1) will identify areas with ancestral haplotypes typical of the Olén area, occurring together with tip haplotypes from clade 4.3. Range expansions and secondary contact Fig. 5 Inferred patterns of range expansion within the island of La Palma for Brachyderes rugatus rugatus. Three range expansions have occurred, each originating from Olén. The broken line 1 indicates an expansion into the Cumbre Vieja terrain, followed by extinction and then subsequent range expansion from the area of Lomo María beginning approximately 0.66 –0.88 million years ago (Ma), indicated by solid arrowed lines. Line 2 denotes a range expansion into the Bejenado terrain approximately 0.33–0.56 Ma, from where a second range expansion into the Cumbre Vieja occurred approximately 0.25 Ma. Line 3 denotes range expansion into the northern shield occurring within the last 1.07 million years. Our analyses have revealed an area of secondary contact between two range expansions of phylogroups 1 and 2 in the Cumbre Vieja terrain of La Palma. An isolated population in the southern part of this region has expanded its range out from the area of Lomo María, while another range expansion has occurred from Olén, through the Bejenado terrain, and then toward the south. We have identified Altos de Jedey as an area where these two © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458 Conclusion Our analyses reveal a more complicated phylogeographic history for Brachyderes rugatus rugatus within the island of La Palma than initial predictions generated from data on the geological history of the island. Rather than a simple demographic involving an origin in the northern shield followed by a more recent expansion into the Cumbre Vieja region, a series of range expansions have occurred, one into the northern shield and two into the Cumbre Vieja (Fig. 5). All three of these range expansions show a common point of origin from the area of Olén, and is tempting to speculate that these have tracked range expansions of the host species Pinus canariensis. However, definitive evidence for this 458 B . C . E M E R S O N E T A L . requires complementary phylogeographic studies for P. canariensis itself or other species facultatively associated with P. canariensis such as Rhyncolus crassicornis (Coleoptera, Curculionidae), Buprestis bertheloti (Coleoptera, Buprestidae), Temnoscheila pini (Coleoptera, Trogossitidae), Leipaspis pinicola (Coleoptera, Trogossitidae) or Orsillus pinicanariensis (Hemiptera, Lygaeidae). Our analyses demonstrate the utility of network-based analyses of intraspecifc DNA sequence variation for testing phylogeographic predictions on active volcanic islands, with geological data providing some explanatory power for the phylogeographic patterns found. But further to this previous geological speculation of a much older age for the Cumbre Vieja region has gained phylogeographic corroboration from our analyses. Both this study and recent analyses of Galápagos giant tortoises on the island of Isabela by Beheregaray et al. (2003), demonstrate the reciprocal illumination of geological and phylogeographic history. Acknowledgements We are grateful to the Cabildo Insular de La Palma for their permission to collect and logistical support, and to Miguel de Navascúes for assistance with collecting. We also thank Antonio Camacho and Rafael García for providing some specimens. We also thank Godfrey Hewitt and three anonymous referees for helpful comments on the manuscript. This work was supported by BBSRC grant 83/G17530 (B.E.). References Ancochea E, Hernan F, Cendrero A et al. (1994) Constructive and destructive episodes in the building of a young oceanic island, La Palma, Canary Islands, and genesis of the Caldera de Taburiente. Journal of Volcanology and Geothermal Research, 60, 243–262. Beheregaray LB, Ciofi C, Geist D et al. (2003) Genes record a prehistoric eruption in the Galápagos. Science, 302, 75. Brower AVZ (1994) Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proceedings of the National Academy of Sciences, USA, 91, 6491–6495. Brown RP, Campos-Delgado R, Pestano J (2000) Mitochondrial DNA evolution and population history of the Tenerife skink, Chalcides viridanus. Molecular Ecology, 9, 1061–1069. Carracedo JC, Day S (2002) Canary Islands. Terra, Harpenden, Hertfordshire, England. Clement M, Posada D, Crandall KA (2000) tcs: a computer program to estimate gene genealogies. Molecular Ecology, 9, 1557–1659. Crandall KA, Templeton AR (1993) Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny construction. Genetics, 134, 959–969. De Salle R, Freedman T, Prager EM, Wilson AC (1987) Tempo and model of sequence evolution in mitochondrial DNA of Hawaiian Drosophila. Journal of Molecular Evolution, 26, 157–164. Emerson BC (2002) Evolution on oceanic islands: molecular phylogenetic approaches to understanding pattern and process. Molecular Ecology, 11, 951– 966. Emerson BC, Oromi P, Hewitt GM (2000) Colonisation and diversification of the species Brachyderes rugatus (Coleoptera) on the Canary Islands: evidence from mtDNA COII gene sequences. Evolution, 54, 911–923. Holland BS, Hadfield MG (2002) Islands within an island: phylogeography and conservation genetics of the endangered Hawaiian tree snail Achitinella mustelina. Molecular Ecology, 11, 365–375. Juan C, Ibrahim KM, Oromi P, Hewitt GM (1998) The phylogeography of the darkling beetle, Hegeter politus, in the eastern Canary Islands. Proceedings of the Royal Society of London. Series B, Biological Sciences, 265, 135–140. Pestano J, Brown RP (1999) Geographical structuring of mitochondrial DNA in Chalcides sexlineatus within the island of Gran Canaria. Proceedings of the Royal Society of London. Series B, Biological Sciences, 266, 805–812. Posada D, Crandall KA (2001) Intraspecific gene genealogies: trees grafting into networks. Trends in Ecology & Evolution, 16, 37 –45. Posada D, Crandall KA, Templeton AR (2000) geodis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Molecular Ecology, 9, 487–488. Templeton AR (1998) Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Molecular Ecology, 7, 381–397. Templeton AR (2001) Using phylogeographic analyses of gene trees to test species status and processes. Molecular Ecology, 10, 779–791. Templeton AR (2004) Statistical phylogeography: methods of evaluating and minimizing inference errors. Molecular Ecology, 13, 789–809. Templeton AR, Sing CF (1993) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. IV. Nested analysis with cladogram uncertainty and recombination. Genetics, 134, 659–669. Templeton AR, Beorwinkle E, Sing CF (1987) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. I. Basic theory and an analysis of alcohol dehydrogenase activity in Drosophila. Genetics, 117, 343–351. Templeton AR, Crandall KA, Sing CF (1992) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics, 132, 619–633. Templeton AR, Routman E, Phillips CA (1995) Separating population structure from population history: a cladistic analysis of the geographical distribution of mitochondrial DNA haplotypes in the tiger salamander, Ambystoma tigrinum. Genetics, 134, 659–669. Brent Emerson is a lecturer in Evolutionary Biology at the University of East Anglia with interests in the application of molecular data to interpret phylogenetic history and population dynamics, particularly within island ecosystems. Shaun Forgie is a systematic entomologist with interests in the application of molecular markers to applied questions. Sara Goodacre is a postdoctoral researcher interested in studying patterns of genetic variation within and among species that are the result of colonization, adaptation and diversification. Pedro Oromí is a Professor Titular at the University of La Laguna, his main area of research is the systematics and biogeography of Macaronesian beetles. © 2006 Blackwell Publishing Ltd, Molecular Ecology, 15, 449–458