Fine-scale parallel patterns in diversity of small benthic Arctic charr (Salvelinus alpinus) in relation to the ecology of lava/groundwater habitats

It is critical to study factors that are important for origin and maintenance of biological diversity. A comparative approach involving a large number of populations is particularly useful. We use this approach to study the relationship between ecological factors and phenotypic diversity in Icelandic Arctic charr (Salvelinus alpinus). Numerous populations of small benthic charr have evolved in lava springs in Iceland. These charr appear morphologically similar, but differ in important morphological features related to feeding. We found a clear relationship between diversity in morphology, diet, and ecological factors among populations. In particular, there were clear differences in morphology and diet between fish coming from habitats where the lava spring flowed on as a stream compared to habitats where the lava spring flowed into a pond. Our study shows that ecological factors are important for the origin and maintenance of biological diversity. The relationship between phenotype and ecological factors are observed on a fine scale, when comparing numerous populations that are phenotypically similar. This strongly suggests that for understanding, managing, and conserving biological diversity important ecological variables have to be taken into the account.


Introduction
Scientists realize that there is a strong link between ecological and microevolutionary processes in generating intraspecific diversity, both among and within populations (e.g., Caroll et al. 2007;Fussmann et al. 2007;Post and Palkovacs 2009). Ecological factors are believed to affect the phenotype of organisms through phenotypic plasticity and adaptive evolutionary processes. Commonly, similar phenotypes are found in similar environments. When these patterns are observed within related lineages they are termed parallel evolution (Schluter et al. 2004;Brakefield 2006;Arendt and Reznick 2008;Parsons et al. 2011), but they are described as convergent evolution when they are observed in unrelated lineages (Arendt and Reznick 2008;Parsons et al. 2011).
Studies of intraspecific parallel evolution have commonly focused on a relatively limited number of populations (usually <10) that show large phenotypic differences. Good examples of this are studies on marine and freshwater sticklebacks, Gasterosteus aculeatus, where armor reduction, both reduced number of plates and loss of spines, has evolved repeatedly after colonization in freshwater (Bell and Foster 1994;Schluter et al. 2004;Reimchen and Nosil 2006). Another common example of parallel evolution is the benthic/pelagic morph pairs in northern freshwater fishes (Robinson and Wilson 1994;Skúlason and Smith 1995;Robinson and Schluter 2000). Detailed studies focusing on a range of phenotypic and ecological characters are needed to understand how natural selection works on a fine scale and how parallel patterns emerge.
Northern freshwater fishes are good candidates for assessing the role of natural selection in parallel evolution. Northern freshwater systems are relatively young. They have only been colonized after the most recent glaciation about 10,000 years ago (Skúlason et al. 1999;Snorrason and Skúlason 2004). These systems have few fish species and the colonizing species are presented with a diversity of unexploited habitats and resources. Lack of interspecific competition and high intraspecific competition has created the opportunity for character release, which can result in resource polymorphism within the colonizing species (Robinson and Wilson 1994;Skúlason and Smith 1995;Skúlason 1996, Robinson andSchluter 2000;Schluter 2000;Snorrason and Skúlason 2004). It is believed that this has taken place repeatedly and independently among systems and within species such as whitefish (Prosopium spp., Coregonus spp.), sunfish (Lepomis spp.), Arctic charr (Salvelinus alpinus), and threespine stickleback (McPhail 1994;Robinson and Wilson 1994;Skúlason and Smith 1995;Smith and Skúlason 1996;Snorrason and Skúlason 2004). Phenotypically similar morphs or species are commonly found in similar habitats among different lakes. This supports the hypothesis that natural selection has influenced the evolution of these populations (Schluter 2000). Arctic charr is a good candidate species to test this hypothesis (Johnston et al. 2011;Kapralova et al. 2011;Kristjansson et al. 2011). This species has a northern circumpolar distribution and is found the farthest north of any freshwater salmonid species (McPhail and Lindsay 1970;Scott and Crossman 1976). The species demonstrates great variability in phenotypes, including anadromous, lake resident, and dwarf populations (Johnson and Burns 1984;Skúlason et al. 1992;Brunner et al. 2001;Magnan et al. 2002;Snorrason and Skúlason 2004;Sigursteinsdóttir and Kristjánsson 2005;Klemetsen 2010).
In Iceland, there is unusually high phenotypic diversity of Arctic charr Snorrason and Skúlason 2004). Iceland has a depauperate freshwater fish fauna as a result of its geographic isolation in the mid Atlantic and the short time since the end of the last glacial epoch. The diversity of freshwater habitats in Iceland is substantial and remarkable (Garðarsson 1979). Older bedrock areas in Iceland have been shaped by glaciers and are relatively impermeable to water. These areas are dominated by direct runoff rivers and lakes, which fluctuate greatly in water flow and temperature (Garðarsson 1979). Younger volcanic areas of the island are dominated by recent lava. The lava bedrock is very porous, which makes it rich in groundwater. Groundwater springs are therefore common in lava areas. These springs are usually quite constant in both water flow and temperature, and have a relatively high mineral content (Cantonati et al. 2006). The lava also provides a complex three-dimensional substrate that provides abundant hiding places for invertebrates and fishes, which often occur in high densities (Malmquist et al. 2000). This habitat complexity, in combination with low interspecific but high intraspecific competition, are the suggested causes for the unusually high level of phenotypic diversity of Arctic charr in Iceland (Snorrason and Skúlason 2004).
An interesting aspect of the phenotypic diversity of Icelandic Arctic charr is the frequent occurrence of small benthic phenotypes (<15 cm adult size; Fig. 1). Populations of small benthic Arctic charr are commonly found in groundwater springs in lava rocks of the neo-volcanic zone (Sturlaugsson et al. 1998;Sigursteinsdóttir and Kristjánsson 2005;Egilsdóttir and Kristjánsson 2008;Kristjánsson 2008;Kapralova et al. 2011;Johnston et al. 2012). These populations are found both in spring-fed streams or where springs drain into lakes and/or ponds. These fish are similar in morphology; they are small but robust in body shape, with a subterminal mouth and dark coloration. However, subtle morphological differences are apparent, for example in features related to foraging (Sigursteinsdóttir and Kristjánsson 2005). A recent genetic analysis strongly indicates that these populations have evolved independently and repeatedly (Kapralova et al. 2011). The small benthic charr in Iceland thus provide an opportunity to determine the relationships between phenotypic diversity and ecological factors on a much finer scale and on a greater number of populations than has previously been done.
Our study has three objectives: First, we estimate the phenotypic variability among small benthic Arctic charr populations. Second, we test the hypothesis that morphological patterns among small benthic Arctic charr populations in Iceland are associated with variation in ecological characteristics of these locations. Specifically, we predict differences in fish morphology between spring-fed stream (stream habitat hereafter) and lake/pond (pond habitat hereafter) habitats. Third, we examine extent to which observed morphological diversity is reflected in diet variability among populations.

Material and Methods
We collected small benthic charr from 31 populations that were widely distributed across the volcanic active zone in Iceland (Fig. 2). The fish were collected in the summer months (June-August) of the years 2004-2007. Springs are stable habitats in terms of temperature and water flow, and we do not believe that different years of collection affected the results. We collected fish from both pond habitats, where the spring flows into a pond or a lake, and stream habitats where the spring continues as a stream. We collected a minimum of 30 fish by electrofishing (usually more than 60 individuals) at each sampling location. The fish were sacrificed immediately with an overdose of phenoxylethanol and placed in plastic bags on wet ice until they were frozen (-20 • C within 5 h). The fish were thawed in the laboratory, and fork length measured to the nearest millimeter. The fish were then blotted on a paper towel and weighed (0.1 g). A high-resolution digital photograph (Nikon CoolPix 800, 3.2 megapixels, Nikon corporation, Miyagi, Japan) was taken of the left side of each fish. We digitized 22 landmarks on each digital image using the tpsdig program (The tps program package; F. James Rohlf; http://life.bio.sunysb.edu/morph.). Six of those landmarks were sliding landmarks (Fig. 3), the other 17 were fixed landmarks. Sliding landmarks are landmarks that are allowed to slide to the left or right along a curve to minimize the shape change between the procrustes average of all the specimens and each specimen. The fish were then dissected and the stomach removed.
We randomly selected 18 of the populations for study of stomach contents. In these populations (Fig. 2), we blindly selected 30 fish for analysis of stomach contents. The stomachs were opened and all food items counted and identified to the lowest possible taxonomic status. The proportion of diet categories in each stomach was calculated by dividing the number of items in that category by the total number of items in the stomach. But this minimizes any effects of size differences. To estimate diet selection, we calculated Ivlev's electivity index (Ivlev 1961) that compares proportion of diet categories in the stomach with those in the environment. The selection index was averaged over all food groups within each individual and compared among populations and between main habitat types using one-way analysis of variance (ANOVA).
At each sampling locality, we measured water temperature, conductivity (μs; ±0.1), and pH (±0.1) using a HI 98129 Hanna meter(Hanna instruments, Smithfield, Rhode Island, USA). Current velocity (m/s; ±0.1) was measured using a Marsh-McBirney Flo-Mate Model 2000 T . The proportion of the bottom surface covered by hard substrate, lava, or large (>5 cm) diameters rocks was visually estimated. The physical complexity (i.e., roughness) of the habitat was estimated using three methods. The first method measured the horizontal extent of a 5.2-m long chain (1-cm links) laid on the surface of the substrate parallel and perpendicular to the shore, the rougher the substrate the shorter the horizontal distance. The second method used a bed profiler to contour the bottom (Young 1993). The profiler (1.2 m long) had 40-cm-long pins spaced at 3-cm intervals, and held vertically above the bottom so that each pin contacted the bottom. We took a digital photograph every 50 cm to record the vertical positions of the pins along two 5-m transects (one parallel and one perpendicular to the shore) at each location. A horizontal line was added at the lowest pin height on each photograph and the distance from that point to each pin was calculated using Sigma Scan pro 5. The standard deviation (SD) of these pins was used as an indication of roughness with higher SDs indicated more roughness. The mean of the SD from all profiles c 2012 The Authors. Published by Blackwell Publishing Ltd. was used for each site. Lastly, the surface roughness of the six stones (see below) was estimated on a scale from 1 to 5, with 1 as smooth and 5 as maximum roughness (Malmquist et al. 2000) ( Table 1). All of these measurements were collected in 25 of the 31 sampling sites. The reason for this was because of logistic problems, and it was a random which stations did not have all of the measurements collected. Six randomly chosen stones were rinsed with water and cleaned with a soft brush to estimate benthic invertebrate abundance. After the first rinsing, a small amount of buffered formalin was added to the water (<0.5%) to force animals to come out of holes and crevices, and then the stones were brushed again. Where soft substrate was present, four mud samples were collected. Each sample consisted of four cylinder cores (27.5 mm diameter). Samples were sieved through a 0.250 mm sieve. The invertebrate samples from both the stones and mud were preserved in 5% buffered formalin and were transferred to 70% ethanol in the laboratory. The animals retained were counted and identified to the lowest possible taxonomic level ( Table 2).
The morphometric data were corrected for up-or down-bending of the specimens using the "unbend" module in the tpsUtil program. The module created a line between landmarks on the snout, end of caudal peduncle, and at the fork of the caudal fin. When fish are bent, this line is curved. The procedure calculates movement of each landmark so that the line becomes straight. Relative warp analysis in tps-relw was used to analyze the variation in morphology, while controlling for geometric body size. This analysis scales the landmarks from each fish to a centroid configuration (mean shape), position, and rotation. The program then defines principal warps from the centroid configuration, which are axes along which shape variation away from the centroid configuration can occur. Partial warps and two uniform components are then calculated (weight matrix) to contain a score for each fish that describes the realized amount of bending and stretching necessary for the configuration of an individual to fit the centroid configuration. The partial warps and uniform components, therefore, describe how individuals differ from the mean along a certain axis of shape variation and were used in further analyses. The program tpsSplin was used to visualize morphological changes.
Discriminant function analysis (DFA) on the weight matrix and diet matrix was used to determine if fish from different populations would be segregated based on morphology or diet. We used MANCOVA to test if population differences were statistically significant. Similar analysis (multivariate analysis of variance [MANOVA] and DFA) was used to examine differences in diet of fish from stream and pond habitats.
Average weight scores and average diet proportions were calculated for each population and the "distribution" of populations was examined using nonmetric multidimensional scaling (NMS, Kruskal 1964), employing the Sörensen distance measure and random starting configuration in PC-ORD 5 (McCune and Mefford 2006). The number of axes was not assigned a priori but was based on calculated stress and instability by the software. The program ran the analysis 10 times with the data set and 20 times with random data from the data set. To examine relationship between morphology and diet with environmental variables, Pearson nonparametric coefficients were overlaid on the vector plot displaying correlations between ordination scores and environmental variables.

Morphology
The discriminant model classifying fish to their population was significant (Wilks-λ (1170) = 0.001, P < 0.001) and classified 72% of the fish to the correct population. The best classification was 94% at Vatnsvik, Thingvallavatn, and Grímsnes, whereas the lowest success was 57% at Trússá and Húsafell area 2. The fish differed for example in trophic structure and caudal structure, as seen in the deformation grids between the average fish and fish from populations at extreme positions (Fig. 4A).
The MANCOVA also found significant differences in fish morphology between stream and pond habitats (F (40,2012) = 17.6, P < 0.01). The discriminant model (Wilks-λ (36) = 0.74, P < 0.001) was significant and correctly classified 73% of the fish. Fish that scored high on the first discriminant axis came from stream habitat. These fish had thicker bodies and caudal peduncle, shorter heads with the curve of the operculum being higher, and had a less subterminal mouth in comparison to the low-scoring fish from pond habitats (Fig. 5A).
Only two sets of habitat variables showed intercorrelation, the board correlated with the chain length (-0.86) and pH with conductivity (-0.63). The NMS on the average morphology of Arctic charr from 25 lakes gave a two-dimensional solution (stress = 16.29, instability = 0.0005, 45 iterations). The correlations between environmental variables and morphological axes were relatively low, between 0.0 and 0.53 (Table 3; Fig. 6A). Scores on NMS axis 1 had a positive correlation with density of a number of benthic invertebrate species, and density of Collembola and Coleoptera (Table 3). Scores on NMS axis 2 showed a negative correlation with temperature, conductivity, and the density of Acarina, Oligochaeta, and Ostracoda, but positive correlation with pH and percent of rock on the substrate (Table 3; Fig. 6A).

Diet
Almost all the fish (97%) had food in their stomach and their diet was diverse (Table 4). Chironomid larvae constituted the most common diet group, ranging from an average of 22% of the dietary items inÁlftavatn to 89% in Klappárós. Other common food groups were cladocerans, copepods, pupae, and flies (Table 4). As seen by the SDs, the proportion of those groups in the stomachs was often quite variable. Groundwater amphipods (Kristjánsson and Svavarsson 2007) were found in the stomachs of 13 fish. The amphipod Crymostygius thingvallensis, (Kristjánsson and Svavarsson 2004) was found in one stomach from Herðubreiðarlindir, and this represents the second finding place of this endemic, newly described species (Kristjánsson and Svavarsson 2007). Cannibalism was observed in one fish in Kaldárbotnar. The diets of fish differed among populations (F (595,7185) = 4.6, P < 0.01). The discriminant model was significant (Wilks-λ (595) = 0.007, χ 2 = 2428, P < 0.001) and correctly assigned 39% of fish to the correct population (Fig. 4B). Correct classification ranged from 80% in Miðhúsaskógur to 13% in the Laekjarbotnar population. There were clear differences among the two habitat types in distribution of scores on the first two discriminant axes. In the pond habitat, the distribution was more variable both within and among populations (Fig. 5B). To study this further, we ran a discriminant analysis within each of the habitat category. The analysis was significant for both the stream (Wilks-λ (290) = 0.06, χ 2 = 805, P < 0.001) and pond (Wilks-λ (198) = 0.02, χ 2 = 725, P < 0.001) habitats. The proportion of fish correctly classified to populations was lower in the stream habitat (42%) than the pond habitat (67%).
A DFA analysis revealed significant differences in the diet of fish in the two main habitats (F (35,480) = 11.5, P < 0.01). The discriminant model was significant (Wilks-λ (35) = 0.5, χ 2 = 302, P < 0.001) and correctly classified 82% of the fish to the correct habitat type (Fig. 5B). There was little overlap in discriminant scores among populations and Figure 4B shows similar patterns as Figure 3B. ANOVA was used to compare diet between habitat types ( Table 5). The number of categories and number of individual prey organisms in the stomachs and proportion of crustaceans (cladocerans, ostracoda, copepoda; Table 5) were all higher in the pond fish, while chironomids were in higher proportion in stream fish.
The NMS on the average proportion of diet groups of Arctic charr from 17 populations showed a one-dimensional solution (stress = 21.2, instability = 0.0003, 36 iterations). The correlation of environmental variables with the NMS axis reveals that the chain transect (0.5, Pearson correlation  A. harpae (-0.6), C. sphaericus (-0.5), Daphnia spp. (-0.5), and Alona spp. (-0.5). Lower values on this axis thus represent a more crustacean diet. Orthocladiinae (a chironomid subfamily) had a positive correlation (0.5), which suggests that higher values on the axis represent a more chironomid larvae diet (Fig. 6B).

Discussion
One of the objectives of this study was to determine whether fine-scale parallel patterns could be detected, that is, is it possible to associate morphological patterns with variation in ecological characters, both physical and biological? The results suggest that this is the case. The extensive lava fields in Iceland offer unique habitats for animals, especially in combination with freshwater (Malmquist et al. 2000). Populations of small benthic charr are common within the volcanically active zone in Iceland. Population genetics data strongly suggest that small benthic charr have evolved independently and repeatedly in different Icelandic freshwater drainages (Kapralova et al. 2011). Hence, it can be argued that their small size, cryptic color, and morphological similarities are to some extent the result of parallel evolution, where the fish have adapted to lava and spring habitats of the neo-volcanic zone. Our common garden experiments (Kristjánsson 2008) indicate that at least parts of the observed phenotypic diversity are heritable. However, the parallel evolution of small benthic charr needs further study, as the relationship between phenotype and fitness traits has not been studied. Our prediction that there are differences in the morphology of the fish from stream habitats versus pond habitats was supported. Fish from pond habitats had narrower bodies, narrower caudal peduncle, longer heads, and more subterminal mouth. These two habitat types were also the main determinants of diet of these small benthic charr, where charr from pond habitats have much more variable diet and eat more small crustaceans than those that live in stream habitats. It is likely that the observed diet and morphological differences are caused by differences in available diet (Govoni 2011) as well as small scale differences in other ecological variables, but our results showed that variables such as temperature, conductivity, and surface roughness were correlated with the observed morphological diversity. Similar relationships between lake ecology and morphological diversity have been seen in Icelandic monomorphic Arctic charr populations .
Our results were not as clear for the detection of patterns between the average morphology and diet and respective environmental characteristics in different locations. However, characteristics such as the number of invertebrate species, conductivity, temperature, and the proportion of rock on the bottom showed indications of relationship with morphological patterns. Detailed discussion on the effects of each factor would at this stage be speculative. The complex lava habitat demand more maneuverability of charr to access food. Fish from habitats with more lava were deeper bodied, had shorter and wider caudal peduncle, and ate more chironomidae than fish from habitats with less lava (Fig. 5).
Even though the small benthic charr live in a relatively cold environment that many would say is barren, the diet of these fish is diverse. Because of water currents and short water retention times, the habitat does not support any dense plankton communities and the charr feeding is thus primarily benthic. Chironomid larvae comprise the most common food, but small benthic crustaceans, such as chydorid Cladocera, copepods, and ostracods are frequently eaten. Fish also eat from the surface, for example, flies and even spiders. In many places the small benthic charr are living partly underground, where they feed on groundwater amphipods.
Small benthic charr populations clearly show parallel patterns of small adult size. Small adult size in fishes has been suggested to result from overcrowding or lack of food (Klemetsen et al. 2002), which results in slower growth and/or earlier onset of sexual maturity. At the same time favorable environmental condition may result in early maturation and small adult body size (Noakes and Balon 1982). Sizedependent predation, where larger fish suffer higher mortality than small fish, has also been suggested to cause earlier maturation and thus small adult size in fishes (Heino and Dieckmann 2008). The small benthic charr populations in Iceland mature early (two to four years; Sandlund et al. 1992;Sturlaugsson et al. 1998;Sigursteinsdóttir and Kristjánsson, 2005;Egilsdóttir and Kristjánsson 2008), which may be the decisive developmental factor causing small adult size. The fish are commonly found in very high densities, for example, where 50-60 fish can be caught within 10 m 2 in shallow waters (<1 m deep) (Kristjánsson 2008). However, this high density does not seem to result in lack of food, as almost all fish examined had food in their stomach and benthic invertebrate densities were often quite high. Low temperature in the environment of these fish is likely to result in relatively slow metabolism and slow growth. The lava habitats having a lot of small holes and fissures sets a constrain on the maximum size of these fish. It is likely that the evolution of small adult  size in these fish is the result of evolution in relation to a combination of these factors.
All populations in the present study had a significant proportion of small mature individual of both sexes, with sexu-ally mature females as small as 8 cm in some populations (B. K. Kristjánsson, personal observation). Subterminal mouth and parr marks that are observed in these populations are indications that the small benthic charr have    evolved through paedomorphosis, which has previously been suggested as an evolutionary trajectory of the benthic morphs (small and large benthic) in Thingvallavatn (Skúlason et al. 1989b;Snorrason and Skúlason 2004). The observed diversity of the small benthic charr may indicate parallel evolution where different populations show similar responses to natural selection. The observed correlation may also be due to plastic responses, which are common in Arctic charr (Skúlason et al. 1989a(Skúlason et al. , b, 1999Adams and Huntingford 2002;Klemetsen et al. 2002;Adams and Huntingford 2004;Snorrason and Skúlason 2004;Kristjánsson 2008;Parsons et al. 2011). Those two factors are, however, not independent of each other as plasticity is an evolvable trait (Pigliucci 2005;Czesak et al. 2006). It is also possible that much of the morphological differences seen among fish from different populations are caused by local genetic factors and unique evolutionary histories. For example, it was common that populations that were close to each other geographically grouped together in the results of the discriminant analysis (Fig. 3). In some of these cases the populations are found in similar habitats, but in other cases the habitat is quite different. This pattern was especially clear when looking at populations from Borgarfjordur. These populations grouped together with high scores on discriminant axis one but come from both pond and stream habitats. This might indicate that local genetic factors may be important for the evolution of the morphological diversity of these fish. These results might also indicate that these habitats were colonized after the charr evolved the small benthic morphotype. However, the results of Kapralova et al. (2011) do not support this, as they detected highly significant genetic structuring within drainages.
As the world's biodiversity faces increased threats, it is important to map diversity and explore factors related to its origin and maintenance. It is now increasingly recognized that ecological and evolutionary processes often act on similar time scales and that microevolutionary responses can be rapid (Hendry and Kinnison 1999;Hairston et al. 2005). The current study has increased our understanding of how key ecological factors may have driven the evolution of small benthic charr in Iceland in a relatively short time. At the same time, the study emphasizes that geological and ecological factors are generating biological diversity. Conservation of biological diversity has most often focused on the conservation of species or other evolutionary significant units. The current findings show that ecological factors are also important on a small scale and how they work to generate diversity within a species. This clearly demonstrates that conservation of biological diversity has to aim for conservation of both the ecological and evolutionary processes (Noakes 2008), both equally important for the origin and maintenance of biological diversity. Therefore, for conservation efforts to be successful we have to conserve habitats and the processes of evolution, not only just focus on the conservation of species.