Morphological divergence of lake and stream Phoxinus of Northern Italy and the Danube basin based on geometric morphometric analysis

Abstract Minnows of the genus Phoxinus are promising candidates to investigate adaptive divergence, as they inhabit both still and running waters of a variety of altitudes and climatic zones in Europe. We used landmark‐based geometric morphometric methods to quantify the level of morphological variability in Phoxinus populations from streams and lakes of Northern Italy and the Danube basin. We analyzed body shape differences of populations in the dorsal, lateral, and ventral planes, using a large array of landmarks and semilandmarks. As the species identification of Phoxinus on morphological characters is ambiguous, we used two mitochondrial genes to determine the genetic background of the samples and to ensure we are comparing homogenous groups. We have found significant body shape differences between habitats: Minnow populations inhabiting streams had a deeper body and caudal peduncle and more laterally inserted pectoral fins than minnows inhabiting lakes. We have also found significant body shape differences between genetic groups: Italian minnows had deeper bodies, deeper and shorter caudal peduncles, and a shorter and wider gape than both groups from the Danube. Our results show that the morphology of Phoxinus is highly influenced by habitat and that body shape variation between habitats was within the same range as between genetic groups. These morphological differences are possibly linked to different modes of swimming and foraging in the respective habitats and are likely results of phenotypic plasticity. However, differences in shape and interlandmark distances between the groups suggest that some (though few) morphometric characters might be useful for separating Phoxinus species.

& Reznick, 2010). Divergent phenotypic adaptations to different habitats aid in an optimal utilization of local resources (Robinson & Wilson, 1994) and can facilitate phenotypic and genetic differentiation of subpopulations finally leading to speciation (Pfennig et al., 2010;Schluter, 2001;Vega-Trejo, Zúniga-Vega, & Langerhans, 2014). The role of phenotypic plasticity in evolution is still under debate (Ghalambor et al., 2015;Price, Qvarnström, & Irwin, 2003). On the one hand, it is regarded as a constraint to speciation processes (Ancel, 2000;Stearns, 1982), while on the other hand, plasticity can also be adaptive and contribute to genetic differentiation and thus speciation (Adams & Huntingford, 2004;Agrawal, 2001). The separation of intraspecific plastic responses to the environment from genetic differentiation is not only of interest for evolutionary biology, but also for conservation biology, an issue of increasing importance (Jacquemin, Martin, & Pyron, 2013).
To summarize, our aims were to (1) assess morphological differentiation between genetically homogeneous minnow populations inhabiting lakes and streams; (2) assess morphological divergence between genetically different populations inhabiting similar habitats; and (3) compare the extent of morphological differences between habitats and genetic groups.

| Sampling and preservation
We analyzed minnows from 10 European populations, from streams (S) and lakes (L) of Northern Italy (ITA) and the Danube basin (DAN; Figure 1). All captured minnows belong to the cyprinid genus Phoxinus and were divided into groups according to habitat and genetic background (see below). Characteristics of all localities and groupings are summarized in Table 1. Throughout the text, the term "lake" is used as a synonym for standing water bodies, while the term "stream" is used for flowing waters. Samples consisted of males and females. To account for possible effects of sexual dimorphism, all analyses were carried out for males and females separately.
All specimens are stored at the Fish Collection of the Museum of Natural History in Vienna (NMW) and include already existing material complemented by recently caught fish (Table 1). From the large collection of Phoxinus samples at the NMW, we selected only those that were in excellent condition, comparable to recently caught ones. The fish sampled in 2014 and 2015 were caught by electrofishing or beach seine and were anaesthetized and then killed with an overdose of MS-222 (Sigma-Aldrich Co., St. Louis, MO, USA), to minimize suffering. Subsequently, the specimens went through an ascending alcohol series and were ultimately preserved in 75% alcohol like all museum specimens. Because preservation leads to shrinkage and weight loss of the specimens, which could affect morphometry (Buchheister & Wilson, 2005;König & Borcherding, 2012;Thorstad et al., 2007), all specimens were stored in 75% alcohol for at least two months prior to scanning, to account for effects of shrinkage. Several studies have shown that the shrinkage follows an exponential function and remains virtually stable after an initial phase of approximately one month (König & Borcherding, 2012;Kristoffersen & Salvanes, 1998;Moku, Mori, & Watanabe, 2004).

| Genetic analysis
According to Kottelat and Freyhof (2007), the Italian populations should be assigned to P. lumaireul and all other examined populations to P. phoxinus. However, because of the taxonomic ambiguity within the genus Phoxinus, and uncertainty regarding morphological characteristics for species determination (see Introduction), we used mitochondrial DNA markers to ensure that we were comparing homogenous groups. Thus, three to five specimens from each population were genetically characterized for two genes, cytochrome b (cyt b) and COI. DNA was extracted from fin tissue using DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol. Polymerase chain reaction (PCR) was performed with primers GluF (5′-AACCAC CGTTGTATTCAACTACAA-3′) and ThrR (5′-ACCTCCGATCTTCGGA TTACAAGACCG-3′) for cyt b (Zardoya & Doadrio, 1999), and FishF1 (5′-TCAACCAACCACAAAGACATTGGCAC-3′) and FishR1 (5′-TAG ACTTCTGGGTGGCCAAAGAATCA-3′) for the COI (Ward, Zemlak, Innes, Last, & Hebert, 2005). Reaction volume was 25 μl, with 2.5 μl buffer, 500 μmol/L dNTPs, 1.5 μl (1.5 μmol/L) Mg2+, 0.25 μl of each primer (25pmol/μl), 0.15 μl TopTaq-polymerase (0.5 units), and 2 μl DNA (with approx. concentration of 10 ng/μl). Cycling conditions for cyt b were as follows: initial denaturation 94°C 3 min, followed by 35 cycles of initial denaturation at 94°C for 30 s, annealing at 51°C for 30 s, and extension at 72°C for 60 s. The final extension was at 72°C for 10 min. For COI, annealing was set to 55°C and extension only lasted 45 s. Purification and sequencing (in both directions) of PCR products was performed by LGC Genomics (Berlin, Germany), with primers used for PCR. The sequences were edited by eye and aligned with MEGA 5.0 (Tamura et al., 2011). The same program was used for calculation of genetic distances between the detected genetic groups.
The museum material used in this study was also genetically characterized, but because museum DNA is typically fragmented, only parts of cyt b were sequenced and used for comparison with the recent material. Fragments of cyt b, of a length between 250 to 350 bp, were amplified using in-house primers and in-house protocols, in accordance with standard procedures for museum DNA (e.g., clean room, extraction, and negative controls, etc.). More details on DNA extraction, PCR conditions, and sequencing are available on request from the authors.

| Group formation
For the genetic and morphological analysis, we compared intergroup and intragroup variations in order to justify the formation of four groups, which were then used for the geometric morphometric and statistical analyses.
On the lateral side, we also used semilandmarks to get information on curves and outlines where the homology criterion of classic landmarks cannot be met . In total, we have used 21 landmarks (+23 semilandmarks) on the lateral side, 14 on the ventral side, and seven on the dorsal side. As the dorsal and ventral sides are not as flat as the lateral side, specimens were tilted cranially, to ensure that the front part of the fish is planar and as close as possible to the scanning glass. As a result, the more caudal body parts could not be used for the analysis and landmarks behind the insertion of the dorsal fin (dorsal side) or anal fin (ventral side) were omitted. A generalized Procrustes analysis (GPA) implemented in tpsRelw (Rohlf, 2010b) was used to standardize the landmark configurations for position, orientation, and size (Mitteroecker & Gunz, 2009;Rohlf & Slice, 1990). The resulting shape differences were illustrated by thin-platespline deformation grids (Bookstein, 1989). We used between-group PCAs (bgPCA; Mitteroecker & Bookstein, 2011) to analyze and display the coordinates obtained from the GPA (Procrustes shape coordinates), which separated the groups better than a standard PCA.
In cases where PCA revealed that principal components (PCs) were influenced by preservation artefacts (i.e., no actual variation in body shape), the respective PCs were removed from the data, by projecting the shape coordinates of the specimens into the subspace perpendicular to the PC (Burnaby, 1966). These artefacts are reflected in lunate or S-shaped bendings and result from shrinkage during preservation or nonplanar storage (Kristoffersen & Salvanes, 1998;Valentin, Penin, Chanut, Sévigny, & Rohlf, 2008) and are repeatedly found in geometric morphometric studies on fishes (Cavalcanti, Monteiro, & Lopes, 1999;Leinonen, Cano, & Merilä, 2011;Sharpe et al., 2008). All further shape analyses, as well as the numbering of the reported PCs, were based on these residual data. The total variance, which is the trace of the corresponding covariance matrix, was used as a measure of shape variation and describes how morphometrically homogenous populations are. To assess statistical significances of group mean differences, we performed Monte Carlo permutation tests (10,000 permutations each), using Procrustes distance (PD) as the test statistic (Good, 2000). Level of significance was determined at α = 5%.
Statistical significances were Bonferroni-corrected, by multiplying the obtained p-values by the number of tests performed. Allometric effects were estimated using a multivariate, pooled, within-group regression of shape on centroid size (Klingenberg, 1998;Mitteroecker, Gunz, Windhager, & Schaefer, 2013). To assess the strength of the group separations, we conducted a multivariate discriminant analysis In addition, we calculated interlandmark distances (Tables S1, S2) to look for traits which may be used as distinguishing features between habitats and genetic groups. All traits were standardized for standard length. Traits on the head or caudal peduncle were also standardized for head length or caudal peduncle length, respectively. We checked whether the respective groups show overlaps of the mean value of each trait ±1 SD, and also ±0.5 SD, as an estimate of the usefulness or reliability of a trait for group delimitation. We refrained from producing p-values, as a statistically detectable significant difference of a trait among groups does not imply that the concerned trait is indeed useful as a distinguishing feature.

| Effects of size, sex, and preservation
Body size of individual fish varied among the populations, with an average length of 53.0 ± 7.8 mm (mean ± standard deviation).
According to different studies (Frost, 1943;Mills, 1987;Mills & Eloranta, 1985), the size of maturity in European minnows is around 40-50 mm; thus, the investigated fish were adults. Allometry (i.e., shape differences related to size) had a significant (p < .001) effect on shape, but explained only 1.52% of the total variation in the data (measured only lateral). Apart from slightly larger eyes of smaller individuals, no major differences were apparent. Due to the very small fraction of explained variance, allometric effects were not further considered in the study.
All body shape differences between habitats and genetic groups were equally expressed in both males and females. The exception is a slightly larger ventral region of females compared with males, which in turn leads to a deeper body. However, body shape changes between habitats or genetic groups that concern the ventral region or body depth (BD) were also expressed when only males were examined. All other traits (e.g., caudal peduncle depth [CPD], head size, length of fin bases) were unaffected by sexual dimorphism. As a consequence, males and females were pooled for all analyses. The PCA revealed that some of the PCs reflected preservation artefacts (bendings) without any actual variation in body shape. As a consequence, PCs 1 and 5 of the lateral data set, PCs 2 and 6 of the dorsal data, as well as PC 2 of the ventral data were removed, as they all depicted lunate-like or S-shaped bendings, but no shape changes that could be attributed to effects of habitat or genetic group (e.g., changes in body proportions).

| Genetic analysis and group formation
Genetic analysis revealed three genetic groups, consisting of (1) Table 2.
For the geometric morphometric analyses, we assigned the pop- to compare genetically identical groups from different environments and genetically different groups from the same environment. As DAN L and DAN S were genetically distant, we only used ITA L and ITA S for the habitat comparisons, to ensure that differences were not attributable to genetic differences.
To justify the groupings from a morphological point of view, we looked at each population individually (e.g., deformation grids, variances, PDs, PCAs) and generally found that intragroup variation was smaller than intergroup variation and that populations from the same habitat and genetic group share similar body shapes.

| Geometric morphometrics
Significant differences in shape were found between ITA L and ITA S populations for the lateral and ventral sides (p < .001), but not for the dorsal side (p = .191 before Bonferroni correction; Figure 3; Table 3).
Regarding the lateral side, ITA L had more slender bodies and caudal peduncles, shorter bases of the dorsal and anal fins, as well as larger heads and eyes. The larger head is primarily due to the increase in the eye, as both snout and postorbital head length remained unaffected.
The mouth, terminal in ITA L, was slightly subterminal in ITA S minnows ( Figure 3b). On the ventral side, ITA L had a distinctly narrower body and gape. The pectoral fins were closer together and also more vertically positioned. Furthermore, the bases of the pectoral fins were longer (Figure 3c).
Highly significant shape differences were found between ITA S and DAN S for all body planes (all p < .001; Figure 4; Table 3).
On the lateral side, DAN S minnows had more slender bodies and also more slender and longer caudal peduncles ( Figure 4b).
Ventrally, DAN S minnows had a distinctly shorter, but wider, gape and a more slender body. The bases of the pectoral fins were thus closer together. In addition, the base of pectoral fin was longer and more horizontally positioned. The distance between the pectoral and the pelvic fins was shorter; however, the distance between the anus and the anal fin was longer (Figure 4c). Minor differences were found on the dorsal side. DAN S samples exhibited a slightly longer and broader head and a shorter head to dorsal fin distance ( Figure 4a).
Highly significant shape differences were found between ITA L and DAN L for all body planes (all p < .001; Figure 5; Table 3). On the lateral side, DAN L had more slender bodies and also more slender and longer caudal peduncles. However, DAN L also exhibited bigger eyes ( Figure 5b). Ventrally, DAN L minnows had a distinctly shorter and broader gape, but a more slender body. However, the bases of the pectoral fins were further apart and of similar relative length than of ITA L. The pelvic fins of DAN L were closer together and their bases shorter. The distance between the pectoral and the pelvic fins was shorter, but the distance between the anus and the anal fin was longer ( Figure 5c). The minor differences on the dorsal side were similar to the changes between stream populations, including a slightly longer and broader gape. The distance from the nape to the dorsal fin was shorter in DAN L (Figure 5a).
A tabular summary of the main shape differences among the groups is provided in the Table S3.
The dorsal plane showed the smallest values of shape variance and also small PDs indicating more similar shapes ( Standard error estimate(s) are shown above the diagonal and were obtained by a bootstrap procedure (1,000 replicates). Analyses were conducted using the Tamura-Nei model (Tamura & Nei, 1993), calculated to be the most appropriate model for the dataset. The differences in the composition bias among sequences were considered in evolutionary comparisons (Tamura & Kumar, 2002). The analysis involved 33 (cytochrome b; cyt b) and 35 (cytochrome oxidase I; COI) nucleotide sequences with a total of 589 (cyt b) and 651 (COI) positions. Evolutionary analyses were conducted in MEGA5 (Tamura, Dudley, Nei, & Kumar, 2007).  Tables S2 and S3 for all distances).

| DISCUSSION
We analyzed body shape differences of populations of the genus Phoxinus in the dorsal, lateral, and ventral planes, using a large array of landmarks and semilandmarks. Previous studies mostly focussed on only one side (usually lateral) and used a small standard set of landmarks (Armbruster, 2012;Brinsmead & Fox, 2002;Simonović, Garner, Eastwood, Kováč, & Copp, 1999). The inclusion of sliding landmarks and all planes allowed a thorough examination of body shape changes and, to some extent, also a cross-validation of the results (e.g., if the same patterns are found on different planes). We found the lowest values for shape variance on the dorsal side, and, furthermore, the DA showed low percentages of correct classifications on this side. The shape changes were mainly attributed to ventral body regions, again supported by the DA. Thus, besides using the (popular) lateral side, we strongly suggest using also the ventral side in similar future studies. The bgPCA and PDs showed that differences in body shape between habitats and genetic groups varied, but was generally within the same range, pointing to strong environmental effects on body shape.

| Body shape differences between habitats (ITA L vs. ITA S)
The present study revealed a deeper body and caudal peduncle, as well as more laterally inserted pectoral fins of stream minnows in contrast to a more streamlined (i.e., more slender) body in lake minnows. Even though it is seemingly paradoxical that stream minnows exhibit a less streamlined body, our findings are in accordance with their habitat preferences. In streams, minnows prefer habitats characterized by slow flowing water, for example, close to the shore (Frost, and anal fins (Ehlinger & Wilson, 1988;Gerstner, 1999;Robinson & Wilson, 1995;Webb, 1984). In contrast, large open water bodies with patchily distributed food favor steady and sustained swimming and thus a more streamlined body with pectoral fins placed more medially (Ehlinger & Wilson, 1988;Langerhans & Reznick, 2010;Walker, 1997;Webb, 1984). The only sexual dimorphic character found was body depth, due to a slightly larger abdomen of the females. However, all body shape differences described above were also present when only males were examined. The general changes in body shape are thus not affected by different sex ratios, or different extents of sexual dimorphism among the populations.
Lake minnows in our study had larger heads (Figure 3), which might be due to habitat-induced changes in head structures linked to different modes of foraging (Ahnelt, Keckeis, & Mwebaza-Ndawula, 2015;Eklöv & Jonsson, 2007;Langerhans, Layman, Langerhans, & Dewitt, 2003;Svanbäck & Eklöv, 2003). However, the head is larger owing to larger eyes, which have been shown to increase the visual acuity (Land & Nilsson, 2002;Wanzenböck, Zaunreiter, Wahl, & Noakes, 1996), both in feeding  and predator avoidance (Goatley, Bellwood, & Bellwood, 2010). The lake minnows in this study were also characterized by a narrower gape. A decrease in gape width supposedly corresponds to higher foraging success in open water, while a larger gape is favorable in feeding on benthic macro-invertebrates (Svanbäck & Schluter, 2012;Walker, 1997). Studies dealing with the diet of minnows found that lake minnows feed primarily on small prey, while stream minnows prefer larger food (Frost, 1943;Michel & Oberdorff, 1995;Straskraba, Chiar, Frank, & Hruska, 1966;Tack, 1940). Consequently, we suggest that our findings reflect divergent adaptations to the respective habitat types in both lake and stream minnows, intertwined with trophic niche partitioning, as have been found for various families, such as Gasterosteidae Hendry & Taylor, 2004;Walker, 1997), Cyprinidae (Jacquemin et al., 2013), or Poeciliidae (Langerhans et al., 2007;Robinson & Wilson, 1995). Nevertheless, Collin and Fumagalli (2011), who examined minnows in Switzerland, found deeper bodies and caudal peduncles in lake populations, that is, the opposite to our findings. These authors attributed shape differences between lake and stream minnows to high predatory pressure by salmonid fishes, implying that the shape changes induced by predatory pressure can overcome those of sustained swimming in open water. Predatory pressure in open waters might also result in a deeper posterior body, which can enhance rapid acceleration (Langerhans et al., 2003(Langerhans et al., , 2007Spoljaric & Reimchen, 2007;Walker, 1997).

| Body shape differences between genetic groups (ITA S vs. DAN S; ITA L vs. DAN L)
By comparing genetically divergent groups from the same habitat, we aimed to identify morphological changes in body shape which are unaffected by the environment. We found differences in overall body shape between Italian and Danube populations, with Italian minnows having deeper bodies and deeper and shorter caudal peduncles in both habitat types. Additionally, all Italian populations had longer jaws, distinctly narrower gapes and snouts, pectoral fin bases originating in a steep angle, and a shorter distance between the anus and origin of the anal fin (Figures 4 and 5). Interestingly, Collin and Fumagalli (2015) F I G U R E 6 Between-group principal component analysis of minnow body shape. Large black dots indicate mean shapes and the ellipses the 90% equal frequency ellipses for stream (bright) and lake (dark) populations in Northern Italy (ITA; dashed blue ellipses, circles) and the Danube basin (DAN; solid orange ellipses, dots). (a) dorsal; (b) lateral; (c) ventral found similar differences in body shape between Alpine minnows (elongate body) and minnows from the Pyrenees (deep body). These identified differences might be a good orientation for morphological species assignment. In addition to differences in body shape, our study indicates that Italian minnows might be differentiated from Danube minnows using genetic methods, thus supporting the reestablishment of P. lumaireul in the Po tributary by Kottelat (2007). Because we only used mitochondrial genes and further molecular studies are needed for phylogenetic and/or taxonomic inference, here we only report trends noticed based on comparisons of genetic distances. As cited in the section (1) of the discussion, body shape differences can also be attributed to environmental factors (i.e., predation) not accounted for in this study. To ensure whether the morphological differences are due to genetic or rather ecological differences (or a mixture of both) would require a common garden experiment. Nevertheless, the changes which we detected might aid further taxonomical studies to find distinguishable characters between the Phoxinus species.

| Comparison of body shape differences between habitats (L vs. S) and genetic groups (ITA vs. DAN)
The results of the DA and the detected shape differences suggest that the extent of morphological differences between genetic groups is within the same range as between different environments (Figures 4-6; Table 4). Furthermore, the shape changes we found between ITA L and ITA S closely resemble the ones between DAN L and DAN S, although considerable genetic distance separates DAN L from DAN S.
Thus, despite the genetic distance of 9% based on cyt b and 5% based on COI (Table 2), the morphological differences detected between DAN L and DAN S appear to be a consequence of habitat. It is conceivable that, in extreme cases, habitat-induced body shape changes (e.g., head length, eye diameter, gape and body width) may mask morphological differences between species/genetic groups (Langerhans & Reznick, 2010;Langerhans et al., 2007;Lucek, Kristjánsson, Skúlason, & Seehausen, 2016;Walker, 1997). This should be accounted for in future approaches on morphological species delimitations in this genus.
Further sampling and analysis of genetically homogenous groups from mixed habitats of the Danube basin should be performed in order to draw any final conclusions.
However, some morphometric characters might be useful for separating Phoxinus species. The few traits we have found, which may be used as distinguishing features, concern either the depth of the body or caudal peduncle, or the insertion of fins. In particular, caudal peduncle depth appears to be a good trait for separating Italian populations from the Danube populations. Nevertheless, all traits generally showed some extent of overlap and may thus only be useful in combination.

ACKNOWLEDGMENTS
We are grateful to Hubert Keckeis for his support in organizing the field trips to Lake Lunz and Wien River and for providing equipment and infrastructure. We also thank Philipp Mitteröcker for his help and support with the geometric morphometric analysis. Our collections at Lake Grundlsee and Lake Lunz were permitted and supported by the Austrian Federal Forest Agency (Österreichische Bundesforste AG) and the Forstverwaltung Kupelwieser, respectively, which is appreciated very much. The collecting activities in Italy were carried out thanks to scientific permission by Servizio Caccia e Pesca della Città Metropolitana di Torino and Provincia di Cuneo, and Regional Natural Parks, Aree Protette delle Alpi Cozie e del Monviso.

CONFLICT OF INTEREST
None declared.

DATA ACCESSIBILITY
Raw landmark coordinates, used in the geometric morphometric analysis, are provided in the Supporting Information. All genetic sequences are deposited in GenBank under accession numbers KX673409-KX673485.
T A B L E 4 Percentages of correct classifications of a multivariate discriminant analysis of the relative warp scores of the first five principal components of each plane, assessed by leave-one-out cross-validation Wilk's lambda is given in the last line of each group (all discriminant functions were highly significant with p < .001). ITA, Northern Italy; DAN, Danube basin; S, stream populations; L, lake populations. See Table 1 and text for further information on group affiliations.