Teleost fish exhibit astonishing examples of adaptive evolution, such as observed in the African cichlids [1], the Neotropical Midas cichlids [2], and the limnetic and benthic stickleback morphs [3]. In general, the colonization of new environments allows rapid diversification [4] as a by-product of adaptation to divergent selection regimes [5], and can finally lead to reproductive isolation [6]. Well known examples of this process in the Northern Hemisphere include the formation of distinct ecotypes in many species pairs of postglacial freshwater fish [7]. Evidence for ecological speciation in these species pairs includes the rapid evolution of reproductive isolation (e.g., separate breeding times, paucity of morphological hybrids) and the parallel evolution of inherited morphological differences that indicate specialization for different niches [7]. Ecotype divergence in the Northern Hemisphere is apparent in numerous species pairs of freshwater fish [8-12], and also includes some anadromous (adults migrating from salt water to spawn in fresh water) and freshwater resident pairs [13,14]. In contrast, published examples of ecotype divergence in the Southern Hemisphere are sparse [11].

Evolutionary processes comparable to those observed in the Northern Hemisphere postglacial lakes have also likely occurred in oceanic island groups (e.g., Hawai'i, Falkland Islands, Tasmania and the Marquesas Islands), where recently formed lakes with a depauperate freshwater fauna were secondarily colonized by diadromous (migratory between salt and freshwater) fish species that then lost their migratory behaviour [15-17]. Compared to the Northern Hemisphere, in the Southern Hemisphere anadromy is much less common [16], and the majority of diadromous species are amphidromous – a special form of diadromy in which only larvae drift to sea and early juveniles (15–50 mm) return to freshwater [18,19]. The adaptive significance of amphidromy is the maintenance of dispersal between isolated, tectonically active island land masses, thereby maintaining gene flow among geographically distant populations [18,20-23]. Accordingly, the loss of migration in amphidromous species leads to geographic isolation and is believed to have initiated genetic and morphological diversification in many taxa [16,24,25]. Well known examples of this process are the freshwater radiations of galaxiid fishes, the diversifications of which have likely been driven by landlocking [26-30]. Consequently, extensive genetic population structuring is observed in several non-migratory species [31,32], including the New Zealand endemic Gobiomorphus breviceps [33]. Of seven New Zealand endemic Gobiomorphus species, three are obligatorily freshwater resident, three are obligatorily amphidromous [34], while only the widespread and facultatively amphidromous [35]Gobiomorphus cotidianus McDowall readily establishes non-migratory populations [34]. One of the obligatorily freshwater resident species (G. alpinus [36]) arose within the last 18,000 years [37], and is closely related to G. cotidianus [36,37]. All New Zealand representatives of the genus Gobiomorphus represent a radiation within the basal Gobioidei [38]. This island Gobiomorphus complex forms a monophyletic group [38], whose ancestor most likely arrived by means of oceanic dispersal ([39], M.I. Stevens & B.J. Hicks, unpublished cytochrome b data). The extensive genetic structuring observed in the non-migratory G. breviceps [33], as well as the non-migratory and recently evolved G. alpinus clearly suggest that the loss of the marine larval life stage facilitates diversification in the New Zealand Gobiomorphus complex.

In gobiids, the structure of the peripheral lateral line canals is an important taxonomic character [40]. In addition, the morphological patterns of these canals can be correlated with particular hydrodynamic stimuli that have direct fitness consequences for fishes (e.g., during rheotaxis, prey detection or predator avoidance) [41-44]. Canal reduction is thought to be an adaptation to distinct microhabitats with slow-flowing water conditions [45], and the congruence between genetic structure and geographic distribution of oculoscapular canal morphotypes in the tidewater goby Eucyclogobius newberry suggests that these variations can be partly inherited [46]. Gobiomorphus cotidianus with predominantly superficial lateral line neuromasts exhibit better detection of moving objects in the absence of background flow [47]. The structure of the oculoscapular lateral line canals in G. cotidianus is highly variable, including complete absence in some lake populations [34]. Furthermore, all obligatorily amphidromous New Zealand Gobiomorphus species have fully developed oculoscapular canals, while these canals are completely absent in all obligatorily non-migratory species [34], a pattern suggesting parallel evolution.

The North Island of New Zealand contains numerous recently-formed lakes [48] that have a diverse history of catastrophic events (e.g., volcanism, glaciations, and sea level changes) that have allowed subsequent colonisations by amphidromous taxa [49]. One of these lakes is Lake Tarawera (Figure 1), a large (39 km²) and deep (90 m) oligotrophic lake located in the geologically highly active central Okataina dome. Lake Tarawera was formed about 10,000 years ago [50]. It drains into the Tarawera River, which is characterized by fast-flowing upper river reaches, including a waterfall (Tarawera Falls, height 65 m, ~2 km from the lake outlet) that represents a significant upstream dispersal barrier for fish into the lake. However, dispersal out of the lake is possible, and major pulses of downstream transport of fish caused by the collapse of lava flows have likely occurred [51]. About 35 km downstream of the waterfall the river enters an extended area of flat land before reaching the Pacific Ocean. No physical barriers limit downstream dispersal of fish within the Tarawera River, whereas upstream dispersal may be limited [52]. The Kaituna River, which originates in the nearby Lake Rotoiti, drains into a coastal area before also reaching the Pacific Ocean. The river mouths of the Tarawera and Kaituna rivers are separated by 50 km of coastline (Figure 1). The Tarawera and Kaituna river systems share no freshwater connections. Like Lake Tarawera, Lake Rotoiti was also affected by volcanic activities. The Rangitaiki River originates from a separate geographic area and drains into the sea about 12 km away from the Tarawera River (Figure 1). Prior to modifications for flood protection (ca 1900) the Tarawera and Rangitaiki rivers likely shared freshwater connections under flood conditions. Following several volcanic eruptions (AD 186 and 1886) that eliminated most freshwater fauna from Lake Tarawera and the surrounding lakes [49], forage fish for trout were introduced into lakes Rotorua and Rotoiti (Figure 1) and from there into Lake Tarawera (around 1900; [53]). These fish likely included G. cotidianus specimens obtained from the Waikato River, a river system originating in the central North Island lake of Lake Taupo (616 km²; Figure 1).

Figure 1 Study area and sample locations for Gobiomorphus cotidianus. A: Lateral view on the head of a male Gobiomorphus cotidianus (© Angus McIntosh, Natural Sciences Image Library, New Zealand). B: Location of the study area in New Zealand (green square) (more ...)

The loss of the migratory life stage is a key process driving speciation in island freshwater fish species [54]. We hypothesize that the loss of the migratory life stage in facultatively amphidromous taxa facilitates intraspecific morphological, behavioural, and genetic differentiation between migratory and non-migratory ecotypes, thereby providing a mechanism for incipient speciation. To date, most studies have focused on interspecific relationships among diadromous species and their non-migratory sister taxa, while little attention has been given to incipient speciation occurring as a result of intraspecific evolutionary processes. Accordingly, the facultatively amphidromous G. cotidianus offers an excellent opportunity to study morphological, behavioural, and genetic diversification between amphidromous and freshwater resident (i.e. non-migratory) populations.

To test for morphological, reproductive, and genetic diversification between migratory and non-migratory stocks, we collected G. cotidianus from one lake and two river sample sites in the Tarawera system (Figure 1). Outgroup samples from the nearby Kaituna and Rangitaiki rivers (Figure 1) were included to examine possible differences between river systems (Table 1). Otolith microchemical analyses of ⁸⁸Sr, ¹³⁷Ba and ⁴³Ca isotopes were used to distinguish migratory and non-migratory stocks. The analysis was complemented with analyses of ¹³⁷Ba/⁴³Ca ratios, as changing levels of this isotope across the otolith indicate diadromy [55]. To test for possible morphological differences between migratory types, the otolith data were contrasted with the distribution of oculoscapular canal morphotypes (as determined from the pore openings of the canals; Figure 2). Previous work suggested the presence of distinct summer and winter spawning populations in the Tarawera River [52]. Thus, we included data on the gonadal development (gonadosomatic indices) to test for temporal reproductive isolation among migratory and non-migratory stocks. All analyses were contrasted with the genetic structure as inferred by Amplified Fragment Length Polymorphisms (AFLPs; [56]), a selectively-neutral, high resolution marker system that can generate a high number of markers distributed genome-wide [57,58]. AFLPs are capable of resolving recent evolutionary splits [59-61], such as expected between different ecotypes. This multidisciplinary approach was applied to reveal patterns of diversification occurring in the New Zealand Gobiomorphus complex as a consequence of loss of migration, and, therefore, the role of this process in driving speciation of island freshwater fish species.

Table 1 Sample details and descriptive parameters for meristic, otolith, reproduction and genetic analyses.

Figure 2 Schematic dorsal view of the oculoscapular canal section of the defined canal morphotypes. The canal structures are evident by presence of pores at their extremities. Median pores (M), lateral pores (La = anterior, Lp = posterior lateral pore) and primary (more ...)

Table 2 Morphotypes in each sample site.

(A) Non-diadromous individuals have low ⁸⁸Sr/⁴³Ca ratios in the nucleus (< 2.5; Figure 3A) with a small range (< 1.5). Their freshwater residency was also supported by a constant level of normalised ¹³⁷Ba/⁴³Ca from nucleus to edge.

Figure 3 Typical patterns of 88Sr and 137Ba counts normalised to 43Ca in otolith cross sections from the nucleus to the edge. A: non-diadromous profile (category A) in a female Gobiomorphus cotidianus (92 mm) from the Kaituna River. B: diadromous profile (category (more ...)

(B) Diadromous individuals have higher ⁸⁸Sr/⁴³Ca ratios in the nucleus (> 4.5) and a larger range between the nucleus and the edge (> 2.0; Figure 3B), illustrating a marine or estuarine larval life stage. In these profiles, the larval migration was also reflected in a characteristic signature of decreasing ⁸⁸Sr/⁴³Ca ratios and increasing ¹³⁷Ba/⁴³Ca ratios from the nucleus to the otolith edge (Figure 3B).

The ⁸⁸Sr/⁴³Ca ratios in the otolith nucleus were lower in non-diadromous fish (mean = 2.2, N = 40) than in diadromous individuals (mean = 11.4, N = 20; ANOVA F_1.58 = 139, P < 0.001). Comparisons of the individual nucleus isotope patterns can be seen in Figure 4. A few individuals could not be confidently allocated to either category based on the descriptive characters used in the scatter plot (Figure 4; solid symbols). However, all amphidromous G. cotidianus exhibited a characteristic increase in ¹³⁷Ba/⁴³Ca from nucleus to the edge, while all freshwater resident individuals showed a constant level of ¹³⁷Ba/⁴³Ca across the otolith (Figure 3). Therefore, the relative differences in individual isotope profiles were used to allocate these individuals to one of the two categories. All Lake Tarawera and upstream Tarawera specimens cluster together (Figure 4) and possess non-diadromous isotope profiles, confirming a complete freshwater life history (Table 1). Similarly, all Kaituna River samples were non-diadromous, as no otolith from this site showed evidence of a marine larval life stage (Table 1). In the downstream Tarawera site, 66% of specimens were diadromous, while from the Rangitaiki River, 73% were diadromous (Table 1).

Figure 4 Nucleus counts against range of normalised 88Sr/43Ca in G. cotidianus and two obligatorily diadromous G. gobioides (red solid dots) included in this study. Colours refer to migratory types: blue = non-diadromous, green = diadromous. Black solid symbols (more ...)

Table 3 Genetic differentiation (FST and θB) between sample sites.

Figure 5 Analyses of genetic structure among sample locations as inferred from AFLP fingerprints. A: Population dendrogram. Numbers on branches are percent bootstrap values out of 1,000,000 pseudo-replicates: Pie diagrams at branch ends illustrate mean proportions (more ...)

Figure 6 Summarized comparison of otolith, meristic, reproductive and genetic analyses. Sample locations are indicated above figure, with each characterized by four vertically arranged bar plots. Horizontal arrangements: Migratory type: light green = diadromous, (more ...)

The large number of scored AFLP fragments (732) enables us to resolve the genetic structure within and between sample sites (Figure 5), even with relatively small sample sizes from the Kaituna and Rangitaiki rivers. These analyses clearly indicate genetic similarity between the Lake Tarawera, upstream Tarawera and Kaituna River sites (Figure 5, Table 3). In all three sample sites, the genetic similarity is paralleled by a dominance of fish with reduced canals (Figure 6). Of all non-migratory individuals, the Kaituna River samples have been collected closest to the sea. Accordingly, the higher proportion of the green genotype in the non-migratory Kaituna River fish might also indicate ongoing hybridization of non-migratory lake (washed out from Lake Rotoiti) and migratory river ecotypes. This is consistent with the dendrogram results (Figure 5; Kaituna River in an intermediate position) and the genetic distances (Table 3; similar distances between Lake Tarawera-Kaituna, and between Kaituna-downstream Tarawera). This intermediate position of the Kaituna River samples is also reflected in a higher proportion, compared to the upstream Tarawera site, of fish with full canals (Type 1). In contrast, the downstream Tarawera and the Rangitaiki River sites, that are genetically close to each other but genetically distinct from the upstream Tarawera and Lake Tarawera sites (Figures 5; Table 3), are dominated by a high proportion of individuals with full canal development (Type 1) that are amphidromous (Figures 5, 6 and Table 3).

Collectively, the concordance of the canal morphotypes with the genetic structure and the reproductive timing suggests that the observed canal morphology is partly inherited rather than a completely plastic response to the environment. But, as in other gobiids [46], we expect that canal development in G. cotidianus is controlled by the interplay of both environmental and genetic factors. Hence, the intermediate pattern of canal formation observed in the upstream Tarawera and the Kaituna River populations may be affected by both a phenotypic plastic response to the ambient river environment, and unidirectional downstream gene flow of non-migratory lake stocks (from lakes Tarawera and Rotoiti respectively) followed by hybridisation with migratory river stocks.

Although we cannot conclude with certainty whether the documented ecotypes are continuing to diverge or are collapsing (such as observed in a Gasterosteus aculeatus species pair [72]), we hypothesise that they are more likely to be diverging because of the apparent temporal reproductive isolation of the two ecotypes. Accordingly, incipient speciation may be occurring, particularly in the lake-locked Lake Tarawera population.

Based on the current data, it is not possible to establish whether the observed morphological and behavioural differences evolved within Lake Tarawera, or in the Lake Taupo/Waikato River system, because it is possible that the non-migratory ecotype originates from the Waikato River and has been introduced into lakes Tarawera and Rotorua/Rotoiti via anthropogenic fish introductions [53]. It is also not yet possible to establish the age of divergence of the two ecotypes. While we cannot rule out that the non-migratory ecotype evolved elsewhere, prior to its introduction, our data clearly support morphologically, reproductively and genetically distinct ecotypes that are derived from distinct migratory and non-migratory stocks.

Our data are congruent, and provide clear and independent lines of evidence for distinct non-migratory and migratory ecotypes on a geographically small scale along a river. The morphological changes observed in the non-migratory populations closely resemble evolutionary patterns repeatedly observed during formation of freshwater resident Gobiomorphus species in New Zealand as well as in other derived gobiid species. These patterns suggest parallel evolution. The present study is, to our knowledge, the first example that clearly suggests that distinct intraspecific ecotypes of an amphidromous fish species (as determined from morphological, reproductive, behavioural, and genetic evidence) may be formed as a consequence of loss of migration and subsequent divergence. Hence, if the reproductive isolation is maintained, these processes could, in the long term, result in the formation of new island freshwater fish species. Future research could focus on several intriguing aspects, including establishing the age of divergence of the ecotypes, exploring the role of landlocking and evolutionary mechanisms (e.g. natural selection, genetic drift) in promoting phenotypic and genetic divergence between the ecotypes, identifying the drivers of parallel evolution in New Zealand Gobiomorphus, and determining the extent to which the patterns observed here are found in other G. cotidianus populations.

Transect readouts of ⁸⁸Sr, ¹³⁷Ba and ⁴³Ca concentrations were carried out in a Perkin Elmer Elan SCIEX DRCII inductively coupled mass spectrometer with a New Wave Research Nd:YAG 213 nm wave length laser at the University of Waikato's Mass Spectroscopy Suite. Laser spot size was 30 μm, with a repetition rate of 20 Hz. A transect of line spots starting from the central nucleus to the edge was performed. Spacing between spots varied between 100–200 μm depending on otolith size. Laser power was set at 50% output with a five second firing and a ten second intersite pause between spots to allow dissipation of background analytes. Between samples the ablation chamber was purged for 90 s with the argon carrier gas. Additionally, the laser was fired at 0% power to standardize against interferences from the carrier gas. Control measurements were subtracted from sample counts per second (cps) to overcome any polyatomic interference. Isotope ratios were calculated from peak-cps for each otolith. Results are presented as dimensionless units for each isotope standardised to counts of ⁴³Ca, an accepted technique in the absence of matrix-matched standards [87]. Results were expressed as line graphs, with each datum point representing the isotope ratio at that point of the otolith. Hence, each line represents the isotope profile across the otolith (Figure 3). Nucleus ⁸⁸Sr/⁴³Ca ratios against range of ⁸⁸Sr/⁴³Ca ratios of each individual profile were plotted to compare variation between individuals (Figure 4).

For each population, the percentage of polymorphic loci and the Shannon-Wiener diversity index (H_SH = -∑ (p_j ln p_j'); where p_j is the frequency of the j-th fragment) are given (Table 1). Fixation indices (F_ST) [89] based on pair-wise distance between individuals (number of shared peaks in AFLP profiles) were calculated in ARLEQUIN v3.01 [90] with significance set to P < 0.05 (tested by 50,172 permutations among groups). Additionally, to detect population structuring the F_ST analogue θ_B was inferred using the software HICKORY [91]. HICKORY calculates θ_B from dominant marker sets based on a Bayesian approach without having prior knowledge of population inbreeding [91]. Data collection was set to a burn-in of 50,000 iterations and data were collected for 250,000 runs. All HICKORY runs were duplicated to ensure repeatability and no significant differences were found between duplicates. Population dendrograms based on F_ST, Reynolds' and Nei's genetic distances were generated in AFLPsurv v1. 0n [92] and consistency of clustering was tested by 1,000,000 bootstrapped distance matrices. A majority rule consensus tree was obtained from the bootstrapped distance matrices with the program routines NEIGHBOUR and CONSENSE from the PHYLIP v3.6 software package [93]. As all distance approaches were consistent we only present the dendrogram based on Nei's genetic distance (Figure 5A). Population structure was inferred with STRUCTURE v2.1 [94] utilizing the admixture model without prior population information. STRUCTURE determines population structure from multilocus genotype data based on a Bayesian clustering approach. During the analysis, STRUCTURE first assumes a number of populations ('clusters', K), then each individual is assigned to these populations, and, finally, for each K a posterior probability (lnP(D)) is given that describes the fit of the data to the respective K (for details about the simulation procedures see [94,95]). To infer the number of clusters K, STRUCTURE was implemented with a series of clusters (K = 1–7). Here, the burn-in was set to 100,000 generations and data were collected for 1,000,000 additional steps. For each K, five independent runs were performed to ensure reproducibility. As the most likely number of clusters present in a dataset is not necessarily indicated by the highest lnP(D) (see STRUCTURE manual for additional details) we applied the method of Evanno et al. [68] to approach the most likely number of clusters K. This method looks for a maximum of the slope (ΔK) of the lnP(D) distribution among runs.

MH, MS, and BJH conceived the study and CM wrote the first draft of the manuscript. MH collected the fish. CM and AC performed the AFLP fingerprinting and CM, BJH and RT collected the isotope and meristic data. CM, BJH, KS and MH performed the data analyses. All authors were involved in writing and data interpretation, and read and approved the final manuscript.

We are grateful to Pete Lockhart and Trish McLenachan for help with lab work, and to David Penny, Bob McDowall, Harald Ahnelt and Reiner Eckmann for helpful discussions. In addition, we are grateful to two anonymous reviewers for their helpful comments and suggestions. CM is grateful to Saira Singh for help and support. This research was supported by Scion (Rotorua, New Zealand), Environment Bay of Plenty (Whakatane, New Zealand) the Allan Wilson Centre for Molecular Ecology and Evolution (Massey University, Palmerston North, New Zealand), and a DAAD (Bonn, Germany) travel grant to CM.