Vivid molecular divergence over volcanic remnants: the phylogeography of Megadromus guerinii on Banks Peninsula, New Zealand

ABSTRACT

Megadromus guerinii, an endemic carabid beetle (Carabidae), is the most common carabid throughout its restricted range on Banks Peninsula, a formation of extinct volcanoes in Canterbury, New Zealand. This study characterises the small-scale phylogeographic patterns of M. guerinii across the formerly volcanically active Banks Peninsula using mitochondrial and ribosomal genes. Between the eastern and western areas of the peninsula, the mitochondrial, but not nuclear, DNA has a well-defined geographic distribution. Specifically, mitochondrial cytochrome c oxidase subunit I (CO1) identifies two distinct groups (> 6% divergence between eastern and western beetles) while ribosomal genes show no discernible pattern. Whether such a pattern represents male-biased dispersal, Wolbachia infection, a recent range expansion of a divergent lineage, or a deeper historic separation is explored. There is potential that male-biased dispersal could have occurred. Wolbachia infection was not detected. We conclude that historical processes have likely separated taxa in the eastern and western peninsula.

KEYWORDS:

• Carabidae
• Coleoptera
• gene flow
• intraspecific variation
• mitochondrial-nuclear discordance
• male-biased dispersal
• Pterostichini

Introduction

An understanding of how intraspecific phylogenetic structure is related to geography and geological processes is pivotal to increasing our understanding of biodiversity, geological history and what is required to manage species effectively (Bergstrom & Chown 1999; Avise 2000; Karl et al. 2011; Kodandaramaiah et al. 2012). Phylogeographic distributions of New Zealand’s flora and fauna are diverse and the underlying causes are often difficult to isolate (Goldberg & Trewick 2011; Trewick et al. 2011). The more recent causes of phylogeographic structure in New Zealand taxa include the uplift of mountain ranges (Buckley & Simon 2007), glaciation (McCulloch et al. 2010) and volcanism (Trewick et al. 2011).

Banks Peninsula is a volcanic landmass on the eastern South Island of New Zealand with an approximate area of 1150 km² (Figure 1). The peninsula originated as a volcanic island between 11 and 5.8 Ma. The volcanic activity had three main centres: the Lyttelton volcano (11–9.7 Ma); the Mt Herbert volcanic group (9.7–7.3 Ma); and the Diamond Harbour volcanic group (7–5.8 Ma) (Sewell 1988; Hampton & Cole 2009). Banks Peninsula was connected to the mainland intermittently as sea levels fell during periods of glaciation and alluvial fans created by the erosion of the rising Southern Alps spread eastward, extending the Canterbury Plains (Brown et al. 1988; Shulmeister et al. 1999; Carter 2005). The elevation of the peaks and ridges surrounding Akaroa has decreased substantially over time from an estimate of over 1500 m to the present average elevation of 760 m (Speight 1944; Boffa Miskell 2007). The eastern slopes of Banks Peninsula are thought to have acted as a refuge for species such as Clitarchus hookeri (White, 1846), the common stick insect, during the last glacial maximum approximately 25,000 years ago (Buckley et al. 2010; Golledge et al. 2012). Today, the peninsula hosts genetically diverse taxa and lineages (Shepherd et al. 2007; Buckley et al. 2009; Shepherd & Perrie 2011) with high levels of local endemism (Ward et al. 1999; Wiser & Buxton 2009; Perrie et al. 2010).

Figure 1. Map of Banks Peninsula with sample locations marked by black triangles.

Megadromus guerinii (Chaudoir, 1865) is a range restricted (confined to a geographic area less than 1000 km²; Townsend et al. 2007) ground beetle endemic to Banks Peninsula. The beetle is thought to be a forest specialist with few known habitat associations besides a bias towards a higher proportion of stones in the topsoil and a negative association with slope (Anderson et al. 2004). The species tends to exhibit clustered distributions; however, over 70% of the variation in its spatial patterns remains unexplained (Anderson et al. 2004). The female of the species will dig out and remain at nest sites, guarding her eggs until they hatch (Anderson et al. 2004). Since M. guerinii is the most common and widespread carabid found across Banks Peninsula, it is assumed to represent a single, interconnected metapopulation within a restricted range (Hansson 1991).

Preliminary molecular work on the mitochondrial gene, cytochrome oxidase c subunit I (CO1) revealed highly divergent individuals within the town of Akaroa and Ellangowan Reserve (S He, Lincoln University, unpubl. data; K Poetz, Lincoln University, unpubl. data). Morphological examination of genitalia could not differentiate between divergent specimens (M Leggett, Lincoln University, unpubl. data), casting doubt on whether a cryptic species is present. The only morphological variant known for the species is a duplication of the anterior marginal setae on the prothoracic shield first described by Britton (1940) from three eastern individuals. High levels of intraspecific genetic variation, such as that observed within M. guerinii, do occur in the literature; for instance Funk and Omland (2003) estimate that 44% of genera contain polyphyletic species. While imperfect taxonomy or failure to recognise cryptic species is attributed as a leading cause of such phenomena (e.g. Witt & Hebert 2000; Sasakawa & Kubota 2005; McKay & Zink 2010; Gillespie et al. 2013), there are empirical cases in which intraspecific variation may be greater than variation between sister species. Numerous examples of high intraspecific genetic distances within New Zealand insects have been documented (Trewick et al. 2011); there are many possible explanations for such patterns, which could be influenced by the dynamic geological history of New Zealand since the Miocene (King 2000).

Possible causes of high intraspecific mitochondrial divergence include male-biased dispersal, with females remaining in isolated pockets and males dispersing between them thereby maintaining nuclear, but not mitochondrial, gene flow across the species (e.g. Avise & Ellis 1986; Melnick & Hoelzer 1992; Hoelzer 1997; Avise 2000; Kerth et al. 2000; Zink & Barrowclough 2008). Mitochondrial pseudogenes, characterised by sporadic indels, frame shift mutations, non-synonymous substitutions and stop codons within protein coding regions can give the illusion of exaggerated mitochondrial divergence; this is especially problematic when inferences are made from small numbers of sequences (Zhang & Hewitt 1996; Bensasson et al. 2001; Song et al. 2008). Furthermore, bacteria from the genus Wolbachia can cause drastic intraspecific divergence within mitochondrial DNA, and speciation via cytoplasmic incompatibility (Schulenburg et al. 2002; Telschow et al. 2002; Jiggins 2003; Telschow et al. 2007; Xiao et al. 2012). For example, in Ceratosolen solmsi (Mayr, 1885) (Hymenoptera), infected and uninfected individuals form two distinct clades, of which the infected group has substantially lower genetic diversity than the uninfected group (Xiao et al. 2012). Furthermore, it is easy to conceive that the past heterogeneity of habitats on Banks Peninsula could have caused populations to become periodically isolated and reconnected in such a way that strong, generalisable phylogeographic patterns were created.

This study aims to: 1. elucidate the nature of genetic variation within M. guerinii by checking for the presence of pseudogenes, Wolbachia and nuclear-mitochondrial discordance; 2. characterise the phylogeographic structure of M. guerinii across Banks Peninsula; and 3. identify candidate areas for genetically divergent taxa within the peninsula defined by historical barriers to gene flow.

Methods

Specimen sampling and molecular methods

Fifty-seven specimens of M. guerinii were collected from a range of forested localities across Banks Peninsula and identified using the taxonomic key for Pterostichini (Britton 1940). Fifty were collected in propylene glycol filled pitfall traps between 17 December 2008 and 28 February 2011, and seven were collected by hand searching between 26 and 30 April 2011. Specimens caught in pitfall traps were pinned and dried prior to the commencement of the study, whereas wild-caught specimens were kept at −20 °C until DNA was extracted. Seven sequences from a preliminary phylogenetic analysis (S He, Lincoln University, unpubl. data) were also included in this study. Genomic DNA was extracted from a crushed leg of each specimen using the Qiagen DNeasy blood and tissue kit with elutions into 100 μl of elution buffer. The mitochondrial gene CO1 was amplified using the primers HCO2198 and LCO1490 (Folmer et al. 1994). The ribosomal internal transcribed spacer 2 (ITS2) was amplified using the Megadromus-specific primers, MegaITS2 F and MegaITS2 R, developed by Fountain (2013). Two regions of 18s rDNA, V4 and V7, were amplified using the primers CV4F and CV4r, and CV7F and CV7R (Raupach et al. 2010), respectively. The D3 domain of 28s rDNA was amplified using the primers CD3F and CD3R (Raupach et al. 2010). The extraction products were screened for the presence of Wolbachia DNA by polymerase chain reaction (PCR) using the protocol of Jeyaprakash and Hoy (2000).

CO1 was amplified in 12.5 µl reactions containing 1.25 µl of DNA template, 0.25 mM of dNTPs, 0.25 mM 10× PCR buffer, 20 µM of each primer and 0.02 unit/µL of taq polymerase (itaq, iNtRON Biotechnologies). Nuclear genes were amplified in 12.5 µl PCRs containing 1.25 µl of DNA template, 0.25 mM of dNTPs, 0.25 mM 10× PCR buffer, 1 μl MgCl_2, 20 µM of each primer and 0.02 unit/µL of taq polymerase (itaq, iNtRON Biotechnologies). All PCRs were performed with negative controls.

The PCR parameters for HCO2198 and LCO1490 were an initial denaturation at 94 °C for 5 min followed by 38 cycles of 94 °C for 45 s, 48 °C for 45 s and 72 °C for 80 s, with a final elongation of 72 °C for 7 min. The cycle used for MegaITS2F and MegaITS2R was identical with the exception of the annealing temperature, which was 67 °C. CV4F and CV4r, CV7F and CV7R, CD3F and CD3R were amplified in accordance with the protocols in Raupach et al. (2010). PCR products were checked by electrophoresis in 1.5% agarose stained with SYBR Safe DNA Gel Stain (Molecular Probes) and visualised under ultraviolet light. Samples that amplified successfully were sequenced in both the 3′ and 5′ direction for each gene with ABI BigDye Terminator v3.1 and resolved in an AVANT 3100 (ABI) sequencer. As there were no indels within the nucleotide sequences, they were aligned using the MUSCLE algorithm (Edgar 2004) in Mega v5.0 (Tamura et al. 2011) and sequences were then compared to the NCBI database (Pruitt et al. 2009) using BLAST (Johnson et al. 2008); CO1 sequences were uploaded to NCBI GenBank (see accession numbers in Table S1).

Phylogenetic and phylogeographic analysis

Preliminary phylogenetic analysis was performed by the construction of neighbor-joining trees using Kimura 2 parameter distances with 10,000 bootstrap replicates in MEGA v5 (Tamura et al. 2011). For Bayesian analysis, the best-fit evolutionary models for nucleotide alignments were determined using jModeltest v0.1.1 (Posada 2008) against 24 models of evolution that were maximum likelihood optimised and subjected to correction via the corrected Akaike’s information criterion (AICc) favouring an HKY + I (Hasegawa et al. 1985) substitution model for all genes. Bayesian phylogenetic analysis was performed in BEAST v1.7.4 (Drummond et al. 2012) using the BEAGLE library (Ayres et al. 2012); three separate runs were performed using the lognormal relaxed clock and coalescent tree prior with 10,000,000 generations per Monte-Carlo Markov Chain (MCMC) (Green 1995). The convergence of each run was analysed in Tracer v1.5 (Rambaut & Drummond 2007). The substitution tree files for each MCMC run were combined in LogCombiner v1.7.4 (included in BEAST 1.7.4) with a 12.5% burnin; a maximum clade credibility tree was then created in TreeAnnotator v1.7.4 (included in BEAST 1.7.4).

Phylogenetic trees were edited in FigTree v1.4.0 (Rambaut 2012). The unrooted Bayesian tree was used to create geophylogenies with 50 m digital elevation maps sourced from LRIS (http://lris.scinfo.org.nz/) in GenGIS v2.4 (Parks et al. 2013). The degree of mitochondrial differentiation and heterogeneity between the eastern and western populations of M. guerinii was calculated with a locus-by-locus analysis of molecular variation (AMOVA) (Weir & Cockerham 1984; Excoffier et al. 1992) and fixation statistics for molecular data (Φ_st) in Arlequin 3.5.1.3 (Excoffier & Lischer 2010) using 50,000 permutations and Tamura & Nei distances (Tamura & Nei 1993) with a gamma a value of 0.011.

Secondary analysis

Nucleotide diversity (π), haplotype diversity (h) and Tajima’s D tests were calculated in DNAsp v5.10.01 (Librado & Rozas 2009) using a 100 bp sliding window analysis at 1 bp steps, ignoring gaps and missing nucleotides. Tajima’s D test was repeated in Arlequin 3.5.1.3 for verification with 10,000 simulated samples. Pairwise genetic distances were extracted from the Bayesian consensus tree and Mantel tests performed with PASSaGE v2.0.11.6 (Rosenberg & Anderson 2011). Genetic landscape interpolations were computed using Alleles in Space v1.0 (Miller 2005) and graphed using ArcGis 10.1 (ESRI 2011). Individuals were sexed by their front tarsi as males possess flattened, wider tarsal pads than females. Chi square tests were used to test for significant differences in male and female capture rates by pitfall traps and hand searching.

Results

Fifty-four individuals were sequenced for CO1. In CO1 575–665 bp were amplified; there were no indels or stop codons; therefore, we accept the CO1 sequences as genuine mtDNA and not numts. Bayesian phylogenetic analysis reached convergence within 10 × 10⁶ generations in three separate MCMC runs. The CO1 phylogeny revealed two distinct groups within M. guerinii (Figure S1). The geographic basis of the two groups was well defined between the western and far eastern peninsula with some overlap on Mt. Pearce and in Akaroa town (Figure 2). AMOVA results showed that the inter-population genetic variance is significantly greater than the intra-population variation (P < 0.00; Φ = 0.90; Table 1). The greatest source of the species diversity for CO1 was concentrated around the hills behind Akaroa (Figure 3). Nucleotide diversity (π), haplotype diversity (h) and pairwise genetic distances (HKY + I) were substantially lower in the western population than in the eastern population (Table 2). Tajima’s D statistic was not significant in the eastern population (D  = −0.62, P > 0.10); however, it was significantly negative in the western population (D = −2.39, P < 0.01). Mantel tests showed significant isolation-by-distance over all M. guerinii (P < 0.00, r = 0.24) and within the eastern group (P < 0.00, r = 0.30); however, no significant geographic correlation was found within the western population (P = 0.15, r = 0.19). Overall, eastern females had a greater correlation between geographic and genetic distance (r = 0.21, P = 0.07) than males (r = 0.12, P = 0.20).

Figure 2. Geophylogeny of Megadromus guerinii CO1 using Bayesian maximum clade credibility tree with taxa nodes mapped to sampling locations on Banks Peninsula. For simplicity, the western and eastern groups are coloured dark red and light blue, respectively.

Figure 3. Phylogeographic divergence of Megadromus guerinii CO1 interpreted as the residual phylogenetic divergence per unit distance interpolation from mid-point values across Banks Peninsula via inverse distance weighting with geographic barriers.

Table 1. Analysis of molecular variance (AMOVA) table generated in Arlequin v3.5 using K2P distances with no added gamma parameters.

Source of variation	d.f.	Sum of squares	Variance components	Percentage of variation
Inter-population	1	432.245	14.76697 Va	89.42944
Intra-population	58	100.961	1.74564 Vb	10.57056
Total	59	533.206	16.51242	100
Φst	0.89429
P-value	< 0.00000

Table 2. Diversity measures and Tajima’s D statistics for the entire population and two subpopulations of Megadromus guerinii.

Population	Nucleotide diversity π	Haplotype diversity h	Tajima’s D
Western population	0.00142 ± 0.00046	0.483 ± 0.115	−2.39149*
Eastern population	0.00833 ± 0.00082	0.944 ± 0.024	−0.61901
Entire population	0.02832 ± 0.00066	0.871 ± 0.038	1.75361

*Indicates statistically significant (P < 0.05) D values.

Sequences from V4 (n = 2), V7 (n = 4), 18s rDNA and D3 subunit of 28s (n = 4) were identical between all specimens of M. guerinii regardless of their CO1 grouping. ITS2 had negligible variation (1 bp) that was not consistent with mitochondrial groupings (n = 19). A total of 79 specimens were observed in the field and in pitfall traps, totalling 43 males and 36 females. Both chi square tests found a significant male bias in pitfall trap captures (X² = 6.13, d.f. = 1, P = 0.01). Wolbachia were not detected for any specimens.

Discussion

Two distinct and well supported mitochondrial lineages were identified in CO1 for M. guerinii. While there are examples of differences in species distribution between the eastern and western peninsula (Townsend et al. 1997), the possibility of a common, and presumably historical, barrier to invertebrate gene flow should be further investigated. The groups show evidence of geographic structuring between the eastern and western peninsula, while also possessing markedly different phylogeographic patterns—most noticeably the lack of any structure within the western group. Although a range expansion could explain the wide distribution of closely related mitochondrial genotypes in the west, the mismatch distribution fits well within the bounds of what would be expected under a constant population size. However, Tajima’s D test statistic was significantly negative, which may indicate a possible selective sweep or a population expansion in the western population.

The phylogeographic structure in the eastern population appears to follow a roughly south to north gradient with the far south east (Flea Bay) being perhaps the most isolated location. It is possible that during the Pliocene glaciation, when the snow line was approximately a kilometre lower (McDowall 2010), and combined with the greater elevation at the time (Shulmeister et al. 1999), that the higher altitude ridges may have been inhospitable for M. guerinii, acting as a barrier to dispersal. An additional possibility that may be worth exploring is ancestral hybridisation between western M. guerinii and M. antarcticus Chaudoir, 1865 causing mitochondrial introgression. Megadromus antarcticus was only found in the western pitfall traps and records of M. antarcticus from the eastern peninsula are scant at best (H Wilson, Hinewai Reserve, pers. comm.). Preliminary unpublished work based on genetic distances for mitochondrial genes suggest that western M. guerinii are closer to M. antarcticus than they are to eastern M. guerinii; these results are being investigated further.

Due to smaller effective population sizes, and less efficient error correction mechanisms, it is expected that mtDNA should diverge faster than nuclear DNA (Chesser & Baker 1996); however, this alone seems insufficient to explain both the depth of mitochondrial divergence and the relative lack of nuclear structure. The inference made from the nuclear-mitochondrial discordance is that either the CO1 divergence occurred rapidly or that, in the presence of male-biased dispersal, gene flow has occurred between the populations, keeping recombinant genes from diverging. Megadromus guerinii females exhibit parental care by digging out a nest site and guarding their eggs (Anderson et al. 2004); it is conceivable that this behaviour may bias females towards remaining in habitats that are conducive for their reproduction, thereby limiting their dispersal ability.

The male bias found in pitfall captures may be evidence of greater mobility in males; this is supported by field observations that males are more often found moving through leaf litter, whereas females are more easily found sheltering under logs (R Emberson, Lincoln University, pers. comm.). Differential activity levels and sex-based dispersal are known to shape genetic structure in other pterostichine beetles (e.g. Rivard 1965; Lagisz et al. 2010). There is some difference in the level of isolation-by-distance between males and females, as would be expected if males dispersed greater distances; however, this is not a significant difference. While the general pattern of high levels of mitochondrial-nuclear discordance are typical of male-biased dispersal (Peters et al. 2012), and there appears to be indirect support for this from pitfall capture bias, any inference made remains tentative until the hypothesis of male-biased dispersal can be empirically tested.

This study concludes that there is no cryptic species of Megadromus Motschulsky, 1865 within M. guerinii. From the combined nuclear and mitochondrial sequence information, we infer that the eastern and western mitochondrial groups may not be entirely reproductively isolated or, if they are, such isolation must have occurred recently and been associated with an elevated rate of mitochondrial evolution. Although this pattern of exaggerated CO1 divergence with nuclear discordance is similar to what may be expected under the influence of Wolbachia infection, male-biased dispersal, or even introgression from previous hybridisation events, possibly with M. antarcticus, no solid evidence is available to make a conclusive inference. The root cause for the discordant phylogeographic patterns in M. guerinii is likely prior geological isolation and perhaps some degree of male-biased dispersal.

Supplementary data

Figure S1. CO1 maximum clade credibility tree for Megadromus guerinii. The red clade are from the west and the blue clade are from the east of Banks Peninsula. Posterior probabilities greater than 0.95 are labelled. The scale bar is in HKY + I distances.

Table S1. NCBI GenBank specimen identifiers and accession numbers for CO1 sequence data.

Acknowledgements

Specimens used in this study were placed in the Lincoln University Entomological Museum or are stored in the Molecular Ecology Laboratory at Lincoln University. We thank Katharina Poetz for assistance with molecular work, John Marris (Lincoln University) for the use of the Entomological Research Museum facilities, and the Banks Peninsula landowners who allowed us to sample on their properties.

Associate Editor: Dr Cor Vink.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was funded by Lincoln University and the JS Watson Conservation Fund.