Little and large: body size and genetic clines in a New Zealand gecko (Woodworthia maculata) along a coastal transect.

Clinal variation can result from primary differentiation or secondary contact and determining which of these two processes is responsible for the existence of a cline is not a trivial problem. Samples from a coastal transect of New Zealand geckos (Woodworthia maculata) identified for the first time a body size cline 7-10 km wide. The larger geckos are almost twice the mass of the small adult geckos. Clines in allele and haplotype frequency were found at two of the four genetic loci examined. Estimated width of the morphological cline was concordant with neither the narrower mtDNA cline (3-7 m) nor the wider nuclear cline (RAG-2; 34-42 km), and cline centers were not coincident. Although the body size cline is narrow compared to the entire range of the species, it is 2-3 orders of magnitude greater than estimates of dispersal distance per generation for these geckos. No evidence of assortative mating, nor of hybrid disadvantage was identified, thus there is little evidence to infer that endogenous selection is maintaining a hybrid zone. We cannot distinguish secondary contact from primary origin of this body size cline but conclude that secondary contact is likely due to the occurrence of mtDNA haplotypes from three distinct clades within the coastal transect and the presence of two frequency clines within this region.

Introduction

Clinal variation in any morphological or genetic trait can result from either primary differentiation or secondary contact, but distinguishing these two processes is notoriously difficult (Endler 1977). Secondary contact of populations that have diverged in isolation is the most common explanation for steep concordant character or frequency clines. Although hybridization (the mating and production of offspring between individuals from distinct populations) results in gene flow, it may not lead to the loss of population differentiation (Barton and Hewitt 1985; Harrison 1993; Wu 2001). At the point of contact between populations or species, a hybrid zone may form and such contacts are often stable due to divergent selection on the distinct parental populations. The width of a hybrid zone depends on the dispersal ability of individuals and the strength of selection. Two types of selection against hybrid individuals will keep a hybrid zone narrow: either endogenous selection based on genetic and/or physiological breakdown, or exogenous selection in which the environment is the agent of selection (Endler 1977; Barton and Hewitt 1985; Barton and Gale 1993; Kruuk et al. 1999; Fitzpatrick and Shaffer 2004). Although most hybrid zones result from secondary contact of populations that have diverged in isolation, it is possible for sharp spatial differentiation to form in primary contact, should (sufficiently strong) selective pressures in different environments promote divergent evolution (Endler 1977). Understanding the significance of clinal variation and distinguishing the process involved in creating and maintaining the cline requires a description of the morphological and genetic distinction among the adjacent populations and estimates of the levels of gene flow, as well as knowledge of the extent of assortative mating and potential hybrid disadvantage.

The New Zealand common gecko (Woodworthia maculata sensu lato, previously known as Hoplodactylus maculatus) is a small nocturnal gecko found throughout the country (Fig. 1). This species, as currently described is considered to comprise a complex of at least 11 cryptic species based on extensive genetic studies (Hitchmough 1997; Nielsen et al. 2011; Hitchmough, Nielsen, and Bauer, unpublished data). However, all the gecko populations investigated here will retain the name W. maculata when the group is revised. The distribution of W. maculata sensu stricto outlined in Nielsen et al. (2011) extends 900 km from northern New Zealand to the Marlborough region in South Island New Zealand. Members of the W. maculata complex are predominantly ground dwelling and are abundant where open stony/rocky habitat and/or dense low shrubs and vines provide them with protection from introduced mammalian predators (Whitaker 1982; Lettink 2007). In forest habitats, the same species can be found in hollow trees but the difficulty of locating them within forests suggests they are at lower density (RAH, personal observation). These geckos have dispersal distances in the order of 1–50 m per generation (Anastadiasis and Whitaker 1987; Lettink 2007). Within the W. maculata complex, size variation among conspecifics has been noted, with adults from coastal populations usually being smaller than those found inland. This size variation is apparently independently derived in a number of lineages (Whitaker 1982; Hitchmough 1997). For example, on the south coast of North Island New Zealand, a population of coastal W. maculata has high densities, sexual size dimorphism (larger males), and small adults compared to geckos on the coast to the northeast and conspecifics found further inland (Whitaker 1982; Hitchmough 1997).

Figure 1

Common gecko (W. maculata) at Turakirae Head, south coast of North Island, New Zealand. Photo by Andrew Blayney.

Less than 15 km separates populations of small and large geckos on the south coast of North Island New Zealand (Fig. 2). Preliminary genetic study revealed one fixed difference out of 27 allozyme loci between the coastal Cape Turakirae population and a population at Lake Pounui (20 km northeast and 5 km from the coast), plus putative frequency differences at six other loci (sample sizes were < 10; Hitchmough 1997). However, the genetic clusters did not precisely follow variation in body size; large-bodied geckos from Ocean Beach shared a General Protein (GP-4) allele with the little-bodied population at Cape Turakirae, while the large-bodied geckos at Lake Pounui were homozygous for a distinct GP-4 allele. Given the continuous nature of the coastal habitat, it was possible for this study to sample from sites at regular intervals between the little and large coastal geckos to determine whether there was a gradual or abrupt change in body size of adults and establish whether they were interbreeding. If little and large geckos are maintaining distinct populations, this could be due to nonrandom mating and/or lower fitness of hybrids relative to parentals. We sought evidence of assortative mating by testing for deviations from Hardy–Weinberg proportions and linkage disequilibrium produced by nonrandom mating. If hybrid offspring are at a selective disadvantage, we expected to detect reduction in levels of heterozygosity between juvenile and adult geckos (Bert and Arnold 1995).

Figure 2

Sampling locations for the gecko W. maculata on the south coast of North Island New Zealand, where there is a steep cline in average body size of adults.

A steep cline should not be assumed to be a hybrid zone unless there is some additional evidence, such as the cline being narrower than an environmental gradient, or that it has the same shape in a number of geographical locations, or there is more variation in individual fitness in the steepest part of the cline (Endler 1977; Barton and Gale 1993). However, hybrid zones are often recognized simply by the existence of concordant (or parallel) clines (Harrison 1993). If there exists a hybrid zone, resulting from recent secondary contact, with little barrier to introgression, concordance and coincidence of frequency clines across the zone could occur, but over time this concordance and coincidence is likely to be lost as loci introgress independently. Alternatively, a morphological cline could have arisen via divergent selection without isolation and genetic diversity could simply be the result of isolation-by-distance.

We sought to differentiate clinal variation concordant with body size variation from simple isolation-by-distance of genetic markers where gene flow is low. Thus, we set up alternative hypotheses: (a) body size variation results from differentiation during a previous period of isolation, or (b) body size variation has arisen in situ due to divergent selection. Evidence for secondary contact (hypothesis a) could come from alleles or haplotypes within the cline belonging to distinct phylogenetic clades when a wider sampling of the species is examined. Thus, we compared the phylogenetic relationships of mtDNA sequence within our transect with a wider sampling of the species (Nielsen et al. 2011). A cline resulting from secondary contact would only be detected if the contact was very recent or it was maintained by hybrid disadvantage (Barton and Gale 1993). Endogenous maintenance of clines can be identified by three main lines of evidence (Ruegg 2008). First, morphological and genetic clines are similar in width (concordant) and location (coincident), which suggests a genome-wide barrier to gene flow (Durrett et al. 2000). We genotyped geckos for three nuclear loci and a mitochondrial DNA fragment, looking for evidence of allele frequency clines and estimated cline widths and centers. Second, the morphological cline would be narrow compared to the estimated dispersal distance, indicating a barrier to gene flow is preventing cline expansion. We compared cline width to dispersal information derived from mark-recapture studies (Anastadiasis and Whitaker 1987). Third, low population density and the presence of both parent and hybrid genotypes in the zone center, suggest a “hybrid sink,” a region into which alleles flow, but are eliminated. Tests for deviations from Hardy–Weinberg proportions and linkage disequilibrium would identify the presence of parental genotypes within the cline center, and selection against hybrids (Phillips et al. 2004). Evidence contrary to these three lines suggests exogenous maintenance and one would expect isolation-by-distance, with neutral diffusion of unlinked markers (Barton and Hewitt 1985). Secondary contact with exogenous selection is difficult to distinguish from primary intergradation (hypothesis b). However, in situ divergence might be implicated if there were an absence of clinal genetic variation concordant with a morphological cline.

Material and Methods

Sampling

Five collecting sites were spaced along 15 km of coastal habitat. The coastal strip is an area of flat ground approximately 500-m wide; the ground then rises steeply at the foot of the Rimutaka Range (rising 500 m per 1 km). The main coastal vegetation is Muehlenbeckia complexa, which provides ideal gecko habitat (Whitaker 1982). However, at Ocean Beach (site 5) the small trees Kunzea ericoides (kanuka) and Coriaria arborea (tutu) dominate. Between October 2008 and June 2009, 143 geckos were caught by searching underneath stones along the coastal transect from Turakirae Head (site 1; Fig. 2) to Ocean Beach (site 5). Our animal handling and sampling protocol was approved by the Massey University Animal Ethics committee (protocol No. 08/63), and authority to study the geckos was provided by the New Zealand Department of Conservation (WE/238/RES). At each site, we aimed to randomly sample 30 geckos. Geckos were released at the capture site after measurements and tissue samples were collected. From the first 20 geckos, both genetic and morphological data were obtained; the last 10 geckos caught from each site were used only for morphology (due to permit restrictions). However, at two locations samples are smaller (site 4 [Gateway] n = 25; site 5 [Ocean Beach] n = 28). Frozen tissue samples from 14 geckos caught in 1987–88 were used only in the genetic analyses: one from site 3 (Fisherman's Rock), five from site 5 (Ocean Beach), and seven from site 6, an inland location beside Lake Pounui, about 5 km from the coastal site 5 (Ocean Beach). One gecko from Featherston (35 km inland) was added to this sample when allele sharing suggested it did not differ significantly (including or excluding this sample had no significant effect on our results or inferences). Geckos at site 6 (Lake Pounui) were found in K. ericoides and Cordyline australis (ti) trees.

For each gecko, sex, eye color, age (adult or sexually immature could be unambiguously determined using external morphological characters), tail state (intact or damaged) were recorded. The number of lamellae on the fourth toe of the hind foot were counted and recorded. Geckos were weighed using mini scales (Tanita Model 1479V). Photos were taken of each gecko to categorize body color and pattern. Six morphological measurements were taken, but given the clear body size differences, only snout-vent length was required for cline analysis (Fig. 3). Analysis of variance (ANOVA) and t-tests were performed using XLSTAT 2009. Analysis of the cline in body size used snout-vent length for adult males and females separately given sexual size dimorphism previously recorded (Whitaker 1982).

Figure 3

Body-size variation (snout-vent length) among samples of adult geckos (W. maculata) collected from five sites on the south coast of North Island, New Zealand. Collecting sites are approximately 4 km apart on a coastal transect. Females (pale gray) and males (darker gray) examined separately (mean, SE, and range indicated).

Genetics

One hundred thirty-two geckos were used for genetic analysis. Genomic DNA was extracted from tail tips (stored in 95% ethanol) using a DNeasy tissue and blood kit (QIAGEN, Hilden, Germany). DNA was diluted to 1 ng/µl for PCR reactions. Reptile primer pairs were used to amplify one mitochondrial and three nuclear loci for every gecko: mitochondrial ribosomal 16S gene (16Sd and 16Sc; Reeder 1995) and nuclear loci, Recombination activating gene-1 (RAG-1; L2408 and H2928; Vidal and Hedges 2004), Recombination activating gene-2 (RAG-2; PY1-F and PY1-R; Gamble et al. 2008), c-mos (MOS-F and MOS-R; Godinho et al. 2006). Although RAG-1 and RAG-2 are known to occur on the same chromosome region in many vertebrates (Fugmann et al. 2006), we treat them as independent loci due to lack of linkage disequilibrium in our dataset (see results). Polymerase chain reactions were performed in 20 µl volumes containing: 200 µM dNTPs, 2 mM MgCl, 1 µM primers, 0.20 U of RedHot DNA polymerase (ABgene, Epsom, UK), and 1–5 ng template DNA. Amplification cycles consisted of—denaturation at 94°C for 60 sec followed by 35 cycles of 94°C for 20 sec, 50°C for 15 sec, and 72°C for 90 sec. Cycle sequencing used Bigdye chemistry (PE) following the manufacturer's protocols, with automated reading on an ABI3730, by the Massey Genome Service (http://genome.massey.ac.nz/), using one of the amplification primers: 16Sc, L2408, PY1-F, or MOS-F.

We used sequencher 4.7 to check and align the sequences, Sequence Aligment Editor (v2.0a11) to check for stop codons and then exported the data as a NEXUS file to be used in McClade 4.0. Genotypes of geckos that were heterozygous at any of the three nuclear loci were resolved by determining which combination of alleles would result in the observed pattern of nucleotide heterozygosity in their DNA sequences. Given the low number of variable sites, this method provided unambiguous genotypes at each locus without need to clone. For example, c-mos alleles differed from each other by a maximum of only three nucleotides and were identified in homozygous individuals; thus, determining the combination of alleles present in each heterozygous gecko was trivial. Sequences were short (<330 bp), therefore, we assumed a lack of recombination within the gene fragments. Genetic diversity analyses were conducted using Arlequin 2.0 (http://lgb.unige.ch/arlequin/software/; Schneider et al. 2000) and linkage disequilibrium between mitochondria and nuclear loci tested, and population differentiation used an exact G test in GENEPOP 4.0.10 (http://genepop.curtin.edu.au/; Raymond and Rousset 1995; Rousset 2008). We compared observed genotype distributions in each site sample with Hardy–Weinberg equilibrium expectations using χ2 tests. Evidence of isolation-by-distance was sought using a mantel test of the correlation of pairwise geographic distances and pairwise F (or Φ for mtDNA) with 1000 permutations in GENEPOP 4.0.10. Each locus was examined separately, and the three nuclear loci were combined.

Heterozygosity levels between adults and juveniles were compared in order to test the prediction that selection against hybrids would result in a reduction of heterozygosity in older animals (Bert and Arnold 1995). We pooled data from all five coastal sites and compared within-site subsamples. Total number of geckos that were heterozygous at each of the three nuclear loci were used to test the prediction that hybrids would have reduced survival. We divided the geckos into four groups; homozygous at all three nuclear loci, or heterozygous at one, two, or three loci. χ2 tests compared the expected ratio of adult heterozygosity levels based on observed juvenile numbers, with the observed adult heterozygosity scores.

For each locus, the DNA sequences were compressed and the redundant taxa and characters removed leaving only the nucleotide sites that differed and one sequence for each distinct allele (or haplotype) using MacClade. Minimum spanning trees of alleles for the three nuclear loci were constructed using the principle of parsimony. The 16S haplotype network was constructed using both a medium-joining approach (Bandelt et al. 1999; http://www.fluxus-engineering.com/sharenet.htm) and Automated Nested Clade Analysis v1.0 (http://www.rubic.rdg.ac.uk/~mahesh/software.html). The resulting networks from these two approaches were identical with the exception of an additional ambiguity within one haplogroup in the medium-joining network.

Representatives of the three haplogroups identified from 16S sequence were chosen and 15 of our geckos were used to generate ND2 sequence for phylogenetic analysis. A fragment of the mitochondrial gene ND2 was amplified and sequenced using primers ND2(Hduv)Fab: CTCCCAGAAACRCTACAAGG and COIR8 (Weisrock et al. 2001), and aligned with 23 sequences from Genbank (W. maculata 16 geckos; W. chrysosiretica 4 geckos; H. duvaucelii 3 geckos). The GTR +Γ (0.236, K = 4) model of DNA evolution was chosen using the AIC criteria and FindModel http://www.hiv.lanl.gov/cgi-bin/findmodel/findmodel.cgi. We used PHYML as a plugin with Geneious 3.5 to generate a maximum likelihood evolutionary hypothesis from 38 sequences (1000 bootstraps).

Cline analysis

For mtDNA data, network analysis of 16S sequence was used to code haplotypes into one of two haplogroups (Table 1; see results), and the frequency of representatives from the two groups at each site used in cline analysis. For each of the three nuclear loci, the frequency of the most common allele was compared across the transect (allele A in each case; Table 2; Fig. 4). For morphology (snout-vent length), hybrid scores were scaled to values between 0 and 1 using the equation (population mean – min. mean)/max. mean – min. mean; Leache and Cole 2007). Cline widths and centers were estimated with Analyse 1.3 (Barton and Baird 1996). Variation in allele frequency was assumed to be predominantly due to a combination of sampling error and a smooth frequency cline (F was set to 0). Confidence intervals were based on any parameter that makes the observed data at least two log-likelihood units from the maximum likelihood (Edwards 1992). A cutoff of 1/7.4 as likely as the maximum likelihood estimates was obtained by randomly varying the parameters (width and center) using a metropolis algorithm (1000 iterations), following recommendations given with the program.

Table 1

Observed mtDNA haplotype frequencies (600 bp;) at six sites across a cline for body-size in coastal geckos Woodworthia maculata. Haplotypes are coded into two groups used for cline analysis (italics and nonitalics) based on minimum spanning network (Fig. 4)

Haplotype	Site l (n = 26)	Site 2 (n = 28)	Site 3 (n = 21)	Site 4 (n = 23)	Site 5 (n = 26)	Site 6 (n = 8)
A	12	22	3	2	3
G	3	1		3
F		5
H	4
N	2
M	2
L	2
R	1
S			1
E			4	1
C			5	7	2
B			8	7	7
D				1	4
V				1
T				1
U					1
W					1
J					3	2
I					2	1
K						2
O						1
P						1
Q						1
Gene diversity	0.76	0.36	0.78	0.81	0.87	0.93
Nucleotide diversity	0.004	0.001	0.009	0.009	0.014	0.006

Table 2

Observed allele frequencies for three nuclear loci at six sites across a cline for body-size in coastal geckos Woodworthia maculata.

Locus	Allele	Site l n = 26	Site 2 n = 28	Site 3 n = 21	Site 4 n = 23	Site 5 n = 26	Site 6 n = 8
RAG-1	A	0.62	0.55	0.64	0.44	0.44	0.50
	B	0.23	0.25	0.26	0.41	0.52	0.19
	C	0.12	0.14	0.05	0.04	0.04	–
	E	0.02	0.05	–	–	–	0.19
	F	0.02	–	–	–	–	–
	G	–	–	0.02	–	–	–
	D	–	–	0.02	0.07	–	0.06
	I	–	–	–	0.02	–	–
	H	–	–	–	–	–	0.06
RAG-2	A	0.23	0.54	0.36	0.65	0.60	0.88
	B	0.25	0.11	0.07	0.07	0.02	0.06
	D	0.14	0.16	0.43	0.09	0.15	0.06
	E	0.06	0.07	0.05	0.15	0.14	–
	C	0.33	0.13	0.07	–	0.08	–
	H	–	–	0.02	–		–
	G	–	–	–	0.04	0.02	–
C-Mos	A	0.94	0.91	0.67	0.96	0.89	0.81
	B	0.06	0.07	0.02	–	–	–
	C	–	0.02	0.31	–	0.06	–
	E	–	–	–	0.02	–	–
	I	–	–	–	0.02	–	–
	D	–	–	–	–	0.02	–
	J	–	–	–	–	0.02	–
	H	–	–	–	–	0.02	0.13
	F	–	–	–	–	–	0.06
Gene diversity		0.481	0.491	0.574	0.428	0.455	0.433

Figure 4

Minimum networks showing putative evolutionary relationships of DNA sequences for all alleles at each of four loci from geckos (W. maculata) collected on the south coast of North Island, New Zealand. Allele and haplotype frequencies within each of the six populations of the coastal transect. mtDNA (600 bp; 16S) haplotype codes use italics and nonitalics to identify two haplogroups, colours (red, yellow, and blue) as in Figure 5.

Figure 5

Evolutionary hypothesis for three New Zealand gecko species based on 1056 bp of mtDNA (ND2) using a Maximum Likelihood analysis with GRT +Γ model of DNA evolution. Map identifies where the geckos were collected, and colour coding corresponds to clades within W. maculata, and bold indicates sequences new to this study.

Results

Morphology

One hundred forty-three geckos were caught and measured from five sites on the south coast of North Island New Zealand (Figs. 1 and 2). At each location, the random sample of 25–30 animals resulted in a mixture of adult and juvenile geckos, of both sexes. Juveniles were excluded from further morphological cline analysis, leaving 101 adult geckos.

Snout-vent length provided us with enough variation to analyze independently from other morphological measures. The average snout-vent length for adult geckos was compared among populations (Fig. 3) and revealed an increase from one end of the coastal transect to the other. For female geckos, this difference is on average about 12 mm. This increase in size along the coastal transect represents an increase in average weight from approximately 4 to 7.5 g (estimates based on nonpregnant animals). Thus, the average mass of adult geckos almost doubles between site 1 (Turakirae Head) and site 5 (Ocean Beach), a distance of only 15 km. At site 3 (Fisherman's Rock), individual adult geckos are intermediate in length. Male geckos also show an increase in size from little at site 1 (Turakirae Head) to large (on average ∼8 mm; Fig. 3). The length difference is smaller than in females because adult males are larger than females at site 1 (Turakirae Head), but similar in size to females (with a nonsignificant tendency to be slightly smaller) at site 5 (Ocean Beach). Significant variation in snout-vent length is seen in six of the 10 pairwise comparisons for adult females and between four of 10 pairwise comparisons for adult males (ANOVA test, female P-value range from 0.0001 to 0.024, male P-value range from 0.008 to 0.024 for significance). In general, the nonsignificant comparisons of average snout-vent length are between adjacent sites, thus the farther geographically the populations are from each other, the larger the size difference.

Variation was seen in gecko eye color, dorsal body pattern, and number of fourth toe lamellae, however, these characters failed to distinguish sites or show clinal variation (Fitness 2010).

Mitochondrial DNA sequences

A 600-bp fragment of mitochondrial 16S was sequenced from 132 geckos. Over the six locations, 23 haplotypes were found with most locations having private haplotypes and a maximum of 0.067 divergence among haplotypes (Genbank accessions HM542429–HM542447). Nucleotide diversity ranges from a high of 0.014 (site 5, Ocean Beach) to the low of 0.001 (site 2, Barney's Whare). The most common haplotype (A) was found at every site except site 6 (Lake Pounui; Fig. 4), and over all six locations Φ = 0.78. All populations showed significant pairwise differentiation from one another for 16S (exact G test, P values ranging from 0.00 to 0.043), with the exception of site 4 (Gateway) compared to site 3 (Fisherman's Rock; P > 0.05). MtDNA differentiation followed a model of isolation-by-distance (mantel test: P = 0.008).

A network of the 16S haplotypes reveals that 10 haplotypes are just one or two mutational steps away from the common haplotype A (Fig. 4). Using the longest branch in the haplotype network, the 23 haplotypes were coded into two clusters (Fig. 3). Each haplotype is at least nine mutational steps away from haplotypes from the other haplogroup, and within six steps of haplotypes within its own group. Only two collecting sites have haplotypes from both groups (sites 4 and 5; Table 1). This coding of 16S haplotypes is used in the cline analysis (below).

A 1056-bp fragment of ND2 was sequenced from 15 geckos, representing the majority of the diversity identified from the 16S haplotypes. The phylogenetic analysis with all previous samples of W. maculata and two outgroup taxa (Nielsen et al. 2011) resolved a monophyletic W. maculata clade (bootstrap = 100%) with haplotypes from the northern extent of their range sister to the southern diversity (Fig. 5). Pairwise genetic distances within W. maculata reached 0.058 (GTR +Γ) between geckos sampled 780 km apart. Within W. maculata five clades were resolved and three of these include haplotypes found within our coastal transect. The 16S haplogroup common in the small geckos from sites 1–3 forms a clade (red) that has not previously been sampled, the 16S haplogroup common in the large geckos from sites 5 and 6 also forms a monophyletic clade (black). The third clade (yellow) within the zone contains the 16S haplotypes N, B, and S, which were common at sites 3, 4, and 5, and previously sampled from Cape Palliser, approximately 30 km east of site 5 on the south coast of North Island. ND2 haplotypes within the transect differed by as much as 0.042 (GTR +Γ).

Nuclear loci

Three nuclear loci had between seven and nine alleles (Table 2). RAG-1 (327bp) with eight (2.44%) variable sites had nine alleles; RAG-2 (268bp) with six (2.24%) variable sites had seven alleles; c-mos (294bp) with seven (2.38%) variable sites had nine alleles (Genbank accession HM542424–HM542428; HQ343302–HQ343308; only alleles found in homozygous form are submitted). There was no evidence of significant departure from the expectations of the Hardy–Weinberg equilibrium for any locus in any of the locality samples, with approximately two-thirds of the point estimates indicating an (insignificant) excess and one-third an (insignificant) deficit (11 vs. 7, respectively; Appendix A1). Site 4 (Gateway) geckos have both the highest frequency of heterozygotes (78%, RAG-1) and the lowest (9%, c-mos; Appendix A1). Pairwise linkage disequilibrium was tested for each locality sample and no evidence was found that the alleles at different loci are nonrandomly associating, with the exception of RAG-1 associated with c-mos at site 1 (Turakirae Head) and c-mos associated with RAG-1 and RAG-2 at site 3 (Fisherman's Rock). However, with Bonferroni correction for multiple tests this is not significant (P > 0.005). Linkage disequilibrium between mitochondria and nuclear loci is also nonsignificant at all sites (P > 0.005).

Appendix A1

Populations of the gecko Woodworthia maculata, from six sites on the south coast of North Island New Zealand, show Hardy–Weinberg expected frequencies of genotypes at three nuclear loci

Locus	Heterozygosity	Site 1 n = 26	Site 2 n = 28	Site 3 n = 21	Site 4 n = 23	Site 5 n = 26	Site 6 n = 8
RAG-1	number of alleles	5	4	5	6	3	5
	observed	0.500	0.5357	0.5714	0.7826	0.731	0.750
	expected	0.547	0.608	0.574	0.633	0.533	0.672
	p(0.05)/(0.005)	NS/NS	NS/NS	NS/NS	NS/NS	*/NS	NS/NS
RAG-2	number of alleles	5	5	6	5	6	3
	observed	0.615	0.648	0.619	0.609	0.654	0.250
	expected	0.756	0.660	0.675	0.538	0.601	0.226
	p(0.05)/(0.005)	*/NS	NS/NS	NS/NS	NS/NS	NS/NS	NS/NS
C-Mos	number of alleles	2	3	3	3	5	3
	observed	0.115	0.107	0.667	0.087	0.231	0.375
	expected	0.110	0.165	0.490	0.084	0.213	0.320
	p(0.05)/(0.005)	NS/NS	NS/NS	NS/NS	NS/NS	NS/NS	NS/NS

Population differentiation (exact G test) for the three nuclear loci used a pairwise distance method for all three nuclear loci pooled together and differences were greatest in nonadjacent populations. Total F over the three loci of 0.078 is significantly greater than 0 (P < 0.05). Correlation between geographic distance and genetic distance (pairwise F) was significant for only one nuclear locus (RAG-2; P = 0.038), and when combined the three loci did not provide evidence of isolation-by-distance (P = 0.168). Only one of the three nuclear loci showed clinal variation when the frequency of the common allele was compared among populations, across the zone (see below). Heterozygosity levels of adults and juveniles (pooled data) did not differ significantly (χ2 test; P > 0.1). When sites were examined separately no significant change in heterozygosity was detected.

Frequency clines were identified for 16S haplogroups and the common RAG-2 allele. Cline centers and widths differ significantly (Table 3). Male and female size clines have the same centers and widths, but these differ from both the mtDNA and RAG-2 clines. The widest cline estimate is for the nuclear marker RAG-2 (> 34 km), the narrowest for mitochondrial haplogroups (< 6.7 km). The center of the mitochondrial cline is significantly further northeast (14 km from site 1) than the other three clines that are centered about 5–8 km from site 1 (close to site 3, Fisherman's Rock). Note that the size clines do not include data for site 6 but because all measurements are taken as distance from site 1, these cline estimates are directly comparable.

Table 3

Estimates of cline-centers and cline-widths for body-size and two genetic loci, in coastal geckos, Woodworthia maculata, on the south coast of North Island New Zealand. tanH curves fitted to cline data were used to estimate centre and width, with best log likelihoods providing support values

Character	Center (km from site 1)	Confidence interval	Width (km)	Confidence interval
Adult male snout-vent length	5.32	5.1–5.6	7.77	7.4–8.8
Adult female snout-vent length	5.37	5.2–5.6	9.18	8.6–9.6
MtDNA (16S)	13.96	13.4–14.8	4.83	3.4–6.7
Nuclear (RAG-2)	7.39	6.9–8.0	38.56	34.9–42.3

Discussion

Body size variation in reptiles is frequently inferred as adaptive and ecotypes are often characterized by differences in size (Thorpe and Malhotra 1996; Leache and Cole 2007; Pincheira–Donoso et al. 2008; but see Brock et al. 2009). Within single species adult size variation can involve a doubling in snout-vent length, however, such variation is usually documented over hundreds or thousands of kilometers (Sears and Agilletta 2004). Only at contact zones between distinct species or ecotypes are body size differences of the scale recorded here observed (Leache and Cole 2007). The adult W. maculata geckos collected from the two ends of our coastal transect are clearly distinct when snout-vent lengths are compared. Adult geckos from the northeast end of the transect (Wairarapa) are longer and heavier than adult geckos from the southwest end of the transect (Wellington). Over the coastal 15 km sampled, the average adult female gecko size changes 12 mm and male geckos change approximately 8 mm. The mass for both males and females approximately doubles from site 1 (Turakirae Head) to site 5 (Ocean Beach). More revealing is the observation that individuals from site 5 are still sexually immature at the maximum size reached by individuals at site 1, and geckos in this coastal habitat commonly live more than 7 years (max > 17 years old; Whitaker 1982; Anastadiasis and Whitaker 1987). Small adult geckos are therefore not necessarily young. Although adult geckos can continue to grow after they have reached sexual maturity, age variation cannot explain the size variation of immature geckos at these sites, and therefore the average size of adult geckos among our collecting sites is apparently the result of different developmental and growth strategies in the two populations at either end of the transect.

Cline or hybrid zone?

As the size of adult geckos from sites 1 (Turakirae Head) and 5 (Ocean Beach), 15 km apart, does not overlap, size of the geckos alone can be used to distinguish these two populations. However, adult geckos collected at sites 2, 3, and 4 were frequently intermediate in size. In all populations morphological variation is high, with polymorphism of toe lamellae number, eye color, color pattern, and hue; however, none of these characters or combinations of characters can distinguish single or groups of populations. We infer that the cline in body size is due to interbreeding and looked for evidence of assortative mating and/or hybrid disadvantage. There is a high level of heterozygosity in all populations and no evidence of reduced fitness of hybrids compared to the parental. A comparison between juveniles and adults show no significant lowering of heterozygosity in adults as would be expected from hybrid disadvantage. Each locus at each site was in Hardy–Weinberg equilibrium. In addition, there is no evidence of linkage of loci as might be expected if parental genotypes are at a selective advantage or if geckos were assortatively mating. This can be compared with a bimodal hybrid zone such as that between the lizards Podarcis bocagei and P. carbonelli, where linkage disequilibrium was detected between many pairs of loci in the contact zone, but not outside it (Pinho et al. 2009). Where linkage disequilibrium is most likely in this study is at site 3 (Fisherman's rock, near the center of three clines), and our tests might lack power resulting in us rejecting a nonrandom association of alleles between c-mos and both RAG-1 and RAG-2 at this site (0.05 > P > 0.005).

Average adult gecko size is clinal and might be used to infer that a hybrid zone exists with its center near site 3 (Fisherman's Rock), where individuals are of intermediate size. Although this study has not sampled at a great distance from the cline ends, we know from previous collections that the large geckos from sites 5 and 6 are representative of the majority of inland W. maculata, and the small geckos found at site 1 are representative of coastal W. maculata forms (Hitchmough 1997). If a hybrid zone does exist it may result from primary or secondary contact.

Evidence for secondary contact could come from haplotypes within the transect belonging to distinct phylogenetic clades. Our mtDNA phylogenetic hypothesis for W. maculata reveals that the area sampled on the Wellington south coast has representatives of not two but three distinct haplotype clades. Identifying three clades is consistent with secondary contact (e.g., Gübitz et al. 2000; Thorpe and Stenson 2003). However, our sampling of W. maculata is as yet too limited for complete confidence that the three clades have distinct geographical ranges that overlap only where the size cline has been identified. Although consistent with secondary contact, the mtDNA phylogeny reveals high diversity in the Wellington region and highlights the need for more extensive sampling throughout the species’ range.

Endogenous or exogenous selection?

Endogenous selection creates a genome-wide barrier to gene flow and is expected to result in concordance and coincidence of clines. Two of the four genetic loci we examined exhibit clinal variation over the 20 km transect. However, we did not find concordance and coincidence of morphological clines with allele frequency clines. The cline in allele frequency at the nuclear locus (RAG-2) is wider than both the adult body size cline and the mtDNA-haplogroup cline. The mtDNA-haplogroup cline is narrower than the body size cline and displaced northeast toward the larger geckos at site 5 (Ocean Beach). The center of the putative cline for nuclear GP-4 allozyme locus, on the other hand, probably lies between site 5 (Ocean Beach) and site 6 (Lake Pounui; Hitchmough 1997). A shift in mtDNA cline, with respect to other markers, has been observed in many hybrid zone studies including Hispaniolan lizards (ca. 34–45 km; Gifford 2008), Swainsons Thrush (ca. 26 km; Ruegg 2008), Fence lizards (ca. 7 km; Leache and Cole 2007), and tree weta (ca. 4–11 km; Morgan–Richards and Wallis 2003). Different cline center and width for mtDNA could arise from the fact that mtDNA is haploid (with a smaller population size) and being maternal it measures only the movement of females, who may move shorter distances than males (Morgan–Richards et al. 2000; Sequeira et al. 2005). Alternatively, the mtDNA cline may have resulted from isolation-by-distance and not secondary contact. At the level of mtDNA haplotype, significant variation is seen among every locality sample.

The estimated center of the nuclear-locus cline (RAG-2) is not coincident with the body size cline for males and females and this cline is significantly wider than both the morphology and mtDNA clines. The difference in cline width may indicate weak selection on the different traits. The cline widths of both the male and female snout-vent lengths are relatively narrow compared to the nuclear DNA, which might result from selection on morphology being stronger than on this nuclear marker, or it may have an independent origin.

A cline that is narrow compared to dispersal distance indicates cline expansion is prevented by a barrier to gene flow. At Turakirae Head (site 1), an extensive mark-recapture study over 5 years found that these geckos take at least 5 years to reach maturity and juveniles disperse further than adults. However, 75% of juveniles moved less than 1 m from where they were first marked and only one juvenile was recorded moving more than 20 m (n = 248 juvenile recaptures; Whitaker 1982). A similar result has been found in a South Island study of W. maculata (sensu lato), where 90% of animals recaptured were less than 10 m away from their original capture site (Lettink 2007). Anastadiasis and Whitaker (1987) found many of Whitaker's (1982) marked individuals still resident in the same study area 7 years after completion of the first study. So dispersal in this species is limited to probably less than 50 m per generation, animals are long lived and gene flow will therefore be low. This low dispersal estimate is in keeping with the significant mtDNA differences seen among almost all samples in our coastal transect. Between the two most distant sites (sites 1 and 6; 20 km), there are no shared haplotypes. Isolation-by-distance is evident from the mtDNA data with overall Φ reaching 0.78. In contrast, although the three nuclear loci (RAG-1, RAG-2, and c-mos) have many private alleles there was allele sharing between the extreme populations. For all four loci, the genetic differentiation between site 1 (Turakirae Head) and site 6 (Lake Pounui) is significant, but of the three nuclear loci, only RAG-2 provides evidence of isolation-by-distance. We can conclude that the character and frequency clines observed are wide compared to the dispersal abilities of the individuals, providing no evidence of secondary contact, expected if there was a barrier to gene flow. However, secondary contact with exogenous selection could also result in a wide zone.

Lack of linkage disequilibrium may be partly explained by the width of the zone relative to the dispersal ability of the geckos because if the parental types do not meet, and early-generation hybrids are absent, all current populations represent the result of introgression, with time for selection against unfit genotypes (Mallet 2005). The age and stability of the cline and slowness of gene flow across it, is indicated by the presence of alleles detected at single locations, particularly in mtDNA. This may result from isolation-by-distance without selection against hybrids. Although isolation-by-distance was detected in mtDNA and RAG-2, it was not evident in RAG-1 or c-mos. Surprisingly, RAG-1 and RAG-2 show no evidence of linkage even though the two genes are adjacent on chromosomes of all vertebrates that have been studied (Fugmann et al. 2006). Selective sweeps might be expected if disruptive selection were involved in creating the body size cline (primary intergradation) and this would increase linkage. However, simulation studies have shown that primary intergradation would lead to rapid loss of clinal variation in neutral loci that are tightly linked to a locus under strong selection (Durrett et al. 2000). We did not detect linkage disequilibrium between RAG-1 and RAG-2 (although putatively adjacent in the gecko genome), as would be expected if the RAG-2 cline resulted from selection in situ. In contrast, a cline in RAG-1 (due to linkage to RAG-2) would be expected to last for thousands of generations if the RAG-2 cline were the result of secondary contact (Durrett et al. 2000).

Conclusions

We have described a body size cline in coastal geckos, the width of which is apparently not controlled by endogenous selection. There is no evidence of assortative mating or hybrid disadvantage and the cline is wide compared to the dispersal ability of the species. However, some evidence points to secondary contact being responsible for the cline formation. We found three mtDNA clades within the transect and two genetic clines. Genetic and morphological clines were neither concordant nor coincident, from which one might infer that this is an old contact zone with independent introgression and exogenous selection. Primary intergradation is difficult to distinguish from secondary contact with exogenous selection and although the allele frequency clines implicate the latter (Barton and Hewitt 1985), the presence of a cline in RAG-2 but not RAG-1 implicates primary intergradation (Durrett et al. 2000).

It is likely that body size is important in determining the fitness of a gecko in a particular habitat, and that the described transition from large to little coastal geckos in New Zealand may be evidence of selection resulting from an environmental gradient. There is no known physical or climatic gradient along the New Zealand coastal transect we sampled, so the role of the environment in limiting gene flow and applying differential selection pressure on these geckos is yet to be determined. Characterizing the habitat change over the transect is an important step in future work to provide a better understanding of selection factors that may be operating across the zone. In addition, sampling inland from the coast to detect parallel clines is important in discriminating the processes involved in cline formation (Johansson et al. 2008).

Due to the level of hybridization, sharing of nuclear alleles and widths of genetic clines, it is unlikely that this size cline is a barrier to gene flow. These populations are therefore not considered different species. No evidence was found of assortative mating (mate choice based on size). But it is possible that selection pressures might result in increased frequency; should such behavior arise, and this could lead to the separate populations further differentiating into different species (Wu 2001; Schluter 2009). We suggest that the cline in body size described here, and low dispersal of this gecko species current is indicative of selection acting on size variation.