Managing cryptic biodiversity: Fine‐scale intralacustrine speciation along a benthic gradient in Alpine whitefish (Coregonus spp.)

Abstract Whitefish (Coregonus spp.) are an important catch for many freshwater fisheries, particularly in Switzerland. In support of this, supplemental stocking of whitefish species is carried out, despite lacking complete knowledge of the extent, distribution and origin of whitefish diversity in these lakes, potentially threatening local endemics via artificial gene flow. Here, we investigate phenotypic and genetic differentiation among coexisting whitefish species spawning along a depth gradient in a subalpine Swiss lake to better delineate intralacustrine whitefish biodiversity. We find depth‐related clines in adaptive morphology and in neutral genetic markers. This individual variation is structured in three distinct clusters with spatial overlap. Individual genetic distances correlate strongly with differences in growth rate and gill‐raker number, consistent with predictions of isolation‐by‐adaptation and ecological speciation. Genetic differentiation between species suggests reproductive isolation, despite demographic admixture on spawning grounds. Our results are consistent with clinal speciation resulting in three species coexisting in close ecological parapatry, one (C. sp. “benthic intermediate”) being previously unknown. A second unknown species spawning in close proximity was found to be of potential allochthonous origin. This study highlights the importance of taxonomically unbiased sampling strategies to both understand evolutionary mechanisms structuring biodiversity and to better inform conservation and fisheries management.

Clinal speciation may be especially important in aquatic environments. In lakes, abiotic factors such as light intensity/composition, oxygen concentration and temperature alongside dependent biotic factors such as the abundance and type of trophic resources, parasites and predators commonly change in a predictable manner: (i) along a depth axis from shallow to deep and (ii) along a benthic-limnetic axis from the sediment surface to open water areas (Seehausen & Wagner, 2014). The higher density of water relative to air also means that individuals in relatively close proximity may experience highly divergent selection regimes, potentially allowing speciation at a finely structured spatial scale, despite high rates of gene flow (Seehausen & Wagner, 2014;Seehausen et al., 2008).
Around 30 species of whitefish in the larger subalpine lakes of Switzerland have arisen following the glacial retreat ~15,000 years BP (Hudson et al., 2011;Steinmann, 1950;Vonlanthen et al., 2012). This negative impact has occurred both through the translocation of non-native whitefish species among lakes and through increased interspecific gene flow via uncontrolled crossing when carrying out supportive breeding (Eckmann, 2012;Hirsch et al., 2013;Hudson et al., 2011). Many of the Alpine lakes harbouring whitefish radiations have recovered from eutrophication ; however, uncontrolled supportive breeding is still routinely carried out. Given the fragility of evolutionarily young adaptive radiations (Seehausen, 2006), there exists strong potential for conflict between maintaining biodiversity and current fisheries management practice. This potential conflict is exacerbated by the fact that thorough, taxonomically unbiased sampling has not been carried out in many lakes, so the true extent of intralacustrine biodiversity may be underestimated. Currently in Lake Lucerne (Figure 1a), a radiation of three endemic whitefish species is known: the large sparsely rakered winter-spawning C. sp. "Bodenbalchen," the small and densely rakered summer-to-winterspawning Coregonus zugensis, and the summer-spawning Coregonus nobilis of intermediate size with a rather high number of gill-rakers (Kottelat & Freyhof, 2007). A fourth winter-spawning species of intermediate size and gill-raker number, C. sp. "Alpnacherfelchen," has also been posited, seemingly confined to Lake Alpnach, a large embayment separated by a 300-m-long channel from Lake Lucerne (Hudson et al., 2011;Svarvar & Müller, 1982). Individuals not corresponding to the phenotype of any known Lake Lucerne whitefish species have been occasionally reported, suggesting the potential for additional unknown species to reside in Lake Lucerne (Douglas & Brunner, 2002;Steinmann, 1950).
Here, we investigate the distribution and structure of adaptive phenotypic and neutral genetic variation among winter-spawning whitefish in Lake Lucerne. We sampled individuals continuously throughout the spawning season along a depth gradient from shallow littoral to deep benthic. Variation in functional morphology and genotypic data from 10 microsatellite loci was quantified and used to compare fish caught on the habitat gradient with whitefish of known species assignation, from distinct spawning aggregations across the lake. This allowed the effective delineation of Coregonus species diversity with fine-scale temporal and spatial resolution and to identify key anthropogenic impacts on whitefish diversity within Lake Lucerne. Given the known presence of C.
sp. "Bodenbalchen" and C. zugensis at the extremes of a winter-spawning depth gradient, we test whether intermediate or indeterminate individuals of whitefish, as reported by fishermen, exist on this gradient and whether these individuals are actually representative of distinct phenotypic and genetic clusters. We subsequently test whether such clusters reflect the known whitefish biodiversity of Lake Lucerne or whether they may constitute additional, previously unknown "cryptic" species. Given the occurrence of anthropogenic stocking among and within Swiss lakes, we then aim to assess the putative origins of the identified clusters, that is whether they represent introduced species or whether they evolved within the native adaptive radiation of Lake Lucerne whitefish.

| Sampling
Lake Lucerne (surface area: 114 km 2 ) is a deep (maximum depth: 214 m; average depth: 104 m) subalpine lake in central Switzerland within the Reuss river system (Figure 1a). Our sampling design comprised single-time-point surveys of whitefish spawning aggregations at different geographic locations throughout the lake, combined with focused fishing of a depth gradient at multiple dates spanning the spawning time/depth range of the two winter-spawning species that were known from the main lake: C. zugensis and C. sp. "Bodenbalchen" (Table 1, Figure 1a). Depth gradient fishing was performed using benthic gill nets, each with a surface area of 250 m 2 . To cover the known body size range of the spawning fish, each net was comprised of three panels of separate mesh sizes set in series: 25, 35 and 45 mm. For each of the five sampling date iterations (19 November 2007, 26 November 2007, 05 December 2007, 11 December 2007 December 2007), nets were set at five different water depths (2, 10, 20, 30 and 40 m) covering the known spawning depths of C. zugensis and C. sp. "Bodenbalchen," as well as intermediate depths ( Figure 1b, these samples are referred to as "depth gradient" in Table 1). One additional sampling replicate using 45 mm mesh size was performed at 2 m depths to increase samples of shallow-spawning large-sized whitefish. This was carried out on 31 December 2007 in the same locality (referred to as "Supp. Large" in Table 1). For each of the main sampling events, pelagic nets with mesh sizes of 38 and 45 mm were also set at 2-7 m depth, close to the depth gradient location over at least 60 m water depth (referred to as "Pelagic" in Table 1). Sampling of other known spawning site locations within Lake Lucerne was carried out to test for isolation-by-distance (IBD) within species (Figure 1; referred to as "Large 1-3," "Int 1," "Small 1-3" in Table 1). Gill nets were set at either at 2-5 m depth using 45 mm mesh size or at 40-60 m depth for 26 mm mesh size. One sample of C. zugensis was obtained outside of spawning season (referred to as "Small n.s." in Table 1) and was only included in analyses that did not require structuring by spawning location or spawning depth. Additionally, spawning aggregates of summer-deep-spawning C. nobilis and of winter-spawning C.
sp. "Alpnacherfelchen" (referred to as "Alpnach" in Table 1) were sampled. All nets were set overnight for around 15 hrs.
All fish were weighed, had their standard length (SL) measured and had a piece of muscle tissue taken and preserved in absolute ethanol.
Ageing was carried out using annual growth rings identified on scales from above the lateral line. Gill-raker counts for each individual were determined using the first left gill arch. Only sexually ripe or almost ripe fish were used in subsequent analyses (maturation degree 5 or 6; Smolina, 1920), with the exception of the "Small n.s." sample. A detailed summary of whitefish sampled, sampling date and location for each sampling site is given in Table 1.

| Sampling data analysis
To test for an effect of sampling depth and spawning time on the structuring of ecologically important phenotypic variation among sampled depth gradient whitefish, linear regressions of standard length (SL) and gill-raker number versus water depth and sampling date of capture were performed in R 3.1.2 (R Core Team 2014). Deviations from a normal distribution in body size (SL) and gill-raker counts for all depth gradient whitefish were tested using a Shapiro-Wilk test. Where deviations from normality were observed, the fit of a mixture of two or three normal distributions was compared based on Akaike's information criterion corrected for sample size (AIC c ), using UiDcmixetHr prSgD v 0.4 (Brewer, 2003). To test whether gaps observed in fish standard length distributions from our catches could be due to body size selectivity of the different gill net mesh sizes employed, skewed normal selectivity curves were estimated for each mesh size using the R package DON and summed to estimate total selectivity across all mesh sizes (Azzalini & Capitanio, 1999;Fujimori & Tokai, 2001). To test whether the deviations from normality of standard length variation within mesh sizes result from the co-occurrence of differently aged fish of the same species or from the co-occurrence of more than one species, three mesh selectivity SL classes were defined based on the selectivity curves for each mesh size obtained above. Boundaries between selectivity SL classes were set at the overall standard length where the selectivity curves between two different mesh sizes crossed. Whitefish were assigned to these classes based on individual SL, irrespective of the actual mesh size each fish was caught in. As within-mesh standard length deviations from normality likely result from the presence of fish from two or more such selectivity SL classes within one mesh size, the age of fish from all three selectivity SL classes was compared using a Mann-Whitney U-test. If deviations from normality were due to the capture of multiple whitefish species with different growth rates within a single mesh size, no age differences would be expected.
The number of SL classes, potentially representing distinct taxon groupings differing in growth rate, was further determined using a dynamic hybrid tree cut (Langfelder, Zhang, & Horvath, 2008), following Lucek, Kristjánsson, Skúlason, and Seehausen (2016). In short, this method is based on a bottom-up algorithm that first identifies preliminary clusters within a data set, depending on a given minimal cluster size, the distance and distinctiveness of its neighbouring objects and the connectivity of branches within a cluster. In a second step, previously unassigned objects are tested for their proximity to the preliminary clusters. This method is based on tree topology without prior assumptions on the number of inferred clusters, therefore providing an unbiased estimate for the number of clusters that are present. Dynamic hybrid tree cut was first conducted on three-year-old whitefish from  Table 1). Three-year-old individuals from other winterspawning whitefish sampling locations were subsequently assigned to the closest matching cluster based on their respective SLs. A dynamic tree cut analysis was further conducted for each remaining age class (ages: 2, 4, 5 and 6), pooling all available individuals within each age group to increase sample sizes. For each analysis, the assumed minimal cluster size was set to 10% of the available individuals.  Figure 3a). Fish from other age classes were excluded, because the trimodal pattern was less evident and/or sample sizes were small ( Fig. S1). Two distinct population genetic methods were used to explicitly test whether sampled winter-spawning whitefish diversity in Lake Lucerne was best explained by the existence of two (K = 2) or three species (K = 3) and to estimate individual assignment probabilities to the identified clusters: Structure v. 2.3.4 (Hubisz, Falush, Stephens, & Pritchard, 2009) and discriminant analysis of principal components (DAPCs) implemented in U g ON et 2.01 (Jombart, 2008).

| Sources of population genetic structure
Parameters used for DetrHcetHr were the following: 100,000 burn-in length, 500,000 MCMC chain replicates, admixture model of ancestry and correlated allele frequencies. Twenty DetrHcetHr run iterations were carried out for each K. Consensus values for individual assignment probabilities across the independent iterations run at each K were calculated in CLUMPAK (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015). All available winter-spawning whitefish from the main lake were used for these analyses (Alpnach, C. nobilis excluded).
Three-year-old fish caught outside the depth gradient site at different lake-wide locations were grouped according to their standard length frequency distribution and the SL-related AMOVA, resulting in three "small-type" (Small 1-3), three "large-type" (Large 1-3) and one "intermediate-type" (Int 1) populations caught in benthic nets, and single populations of both "large-type" and "intermediate-type" fish caught in pelagic nets (see Table 1). Genetic differentiation (F ST ) between these populations and the SL class groupings of fish caught on the depth gradient was then calculated in ra qHiON using 10,000 permutations to assess significance. To test for isolation-by-distance (IBD) among three-year-old fish caught in different parts of the lake, partial mantel tests were performed with the vegan package (Oksanen et al. n/a n/a n/a n/a 38 Alpnach 46°57′52.11″ 8°19′10.49″ 02 December 2004 20 n/a n/a n/a n/a 20 Total  647  242  137  54  51  299 Given for each sampling event are the event name, geographic location, date of sampling, number of whitefish caught (sample size), the number of threeyear-old whitefish that were included in the population genetic analyses (N Year 3 ) and the sample sizes corresponding to each standard length (SL) class (N Small , N Int. , and N Large ) and the number of individuals included in genetic analyses (N Genetics ). 2013) in R, correlating individual genetic distances (Smouse & Peakall, 1999) that were calculated in g ONSUiv (Meirmans & Van Tienderen, 2004) with geographic distances among sampling locations, while controlling for genetic structure among population groupings based on either individual capture depth or individual standard length.
Significance was assessed using 1,000 permutation steps.  (Goudet, 2001). Significance levels of F IS and deviations from HWE were corrected for multiple testing using sequential Bonferroni correction (Rice, 1989). Deviations from linkage equilibrium between all pairs of loci for each species were tested for using ra qHiON. Allelic richness (A R ) and the mean number of private alleles per locus were calculated for each species in ADZE-1.0 (Szpiech, Jakobsson, & Rosenberg, 2008). Generalized

| Genetic differentiation among Lake Lucerne species
To estimate neutral marker differentiation among species groups, multilocus pairwise F ST values were calculated in ra qHiON, using 10,000 permutations. A false discovery rate (FDR) correction was subsequently applied to all pairwise comparisons (Benjamini & Hochberg, 1995).

| Tests of association between neutral genetic and adaptive phenotypic variation
Using all three-year-old individuals sampled from the depth gradient, the correlation between individual genetic differences and individual differences in standard length and gill-raker number was assessed using Mantel tests with 1,000 permutation steps, using the vegan package in R. Partial Mantel tests were further used to test for the same correlation, controlling either for the effects of geographic distance or capture depth.

| Genetic relatedness across lakes
To investigate whether historical whitefish introductions might have contributed to contemporary whitefish diversity, genetic comparisons were made between Lake Lucerne taxa and an additional 12 whitefish species representing four major lake radiations within the Alpine radiation (Lakes Brienz/Thun, Constance, Neuchatel and Walen/Zurich).
Existing microsatellite genotype data (Bittner, Excoffier, & Largiadèr, 2010;Vonlanthen et al., 2012) were combined with Lake Lucerne genotypes to create a consensus population-based neighbour-joining (NJ) tree with 1,000 bootstrap replicates using Cavalli-Sforza chord distances (D CH ) in pyaip v. 3.69 (Felsenstein, 2005). Pairwise F ST values between all species were again estimated in ra qHiON. Additionally, DAPC was used to estimate the overall probability of assignment of Lake Lucerne whitefish to allopatric whitefish populations from other Alpine lakes (Constance, Thun, Brienz, Walen, Neuchatel, Zurich).
Pairwise analyses were carried out separately between the constituent species of Lake Lucerne and those of each of the other six lakes.
For each pairwise analysis, the lake of origin of each fish was used to define prior groups for the subsequent DAPC. The lowest associated root-mean-squared error (RMSE) value obtained following cross-validation was used to identify the optimal amount of principal components (PCs) retained during DAPC. For the genotyped species from each of these lakes, generalized private allelic richness was again calculated in ADZE-1.0, to measure the number of alleles exclusive to different multispecies groupings.

| Sampling results
In total, 647 whitefish from Lake Lucerne were analysed (Table 1).

| Genetic differentiation among Lake Lucerne species
Lucerne species groupings showed significant levels of pairwise genetic differentiation from one another, following FDR multitest correction (F ST = 0.019-0.124; Table 3 (Table 3). For each species grouping, we report the standard length range at age three (SL), the sample size (N), observed (H O ) and expected heterozygosity (H E ), significance level of deviation from HWE across all loci (p-HWE), allelic richness (A R ), inbreeding coefficient (F IS ) and its significance level (p-F IS ), number of deviations from linkage equilibrium (N LD ) at a significance level of p = .05. Significance levels for HWE and F IS are Bonferroni corrected.

| Tests of association between neutral genetic and adaptive phenotypic variation
gill-rakers (r = .342, p < .001) also remained significant when controlling for the effects of geography.

| Genetic relatedness of species across lakes
Neighbour-joining trees show that whitefish species within the Alpine radiation predominately cluster with other coexisting whitefish species, forming monophyletic lake or connected-lake radiations ( Figure   S4). One exception to this pattern is that Lake Constance taxa form a nested grouping within a predominantly Lucerne clade. Here, Lucerne

Pairwise DAPC comparisons between Lake Lucerne and other
Alpine lake species flocks varied among species and among lakes in their strength of assignment to one of the two clusters specified

| DISCUSSION
By carrying out a detailed analysis of ecological spatial structuring within the endemic whitefish radiation of Lake Lucerne, we demonstrate the importance of environmental gradients for speciation and species coexistence, even at fine spatial scales. We delineate six dis-
T A B L E 3 Genetic differentiation between Lake Lucerne whitefish species based on 10 neutral microsatellite loci species found spawning along the depth gradient, we find a strong association of individual genetic variation with growth rate and gillraker counts, traits suggested to be evolving under divergent selection in whitefish (Østbye et al., 2006;Praebel et al., 2013;Rogers & Bernatchez, 2007;Vonlanthen et al., 2009). These findings are consistent with predictions of isolation-by-adaptation and ecological speciation along clines in our studied depth gradient (Funk, Nosil, & Etges, 2006;Nosil, 2012;Rundle & Nosil, 2005). Despite indications of a non-native origin to C. sp. "pelagic intermediate," the full extent of how anthropogenic impacts such as eutrophication, stocking and other fisheries management practices have influenced whitefish diversity present in Lake Lucerne remains unclear. The discovery of two previously undescribed species in this study highlights the importance of taxonomically unbiased sampling strategies to assess the whole diversity present in a system: to both understand the evolutionary mechanisms structuring contemporary biodiversity and to better inform conservation and fisheries management plans.

| The endemic whitefish diversity of Lake Lucerne
In total, we find empirical support for the presence of six genetically and phenotypically differentiated whitefish species in Lake Lucerne.
"Alpnacherfelchen" ("Alpnach"), appears to be restricted to a separate sub-basin of Lake Lucerne, Lake Alpnach (Hudson et al., 2011;Svarvar & Müller, 1982). More research is required to establish the ecological and evolutionary distinctness of the latter taxon from other Lucerne whitefish. The final two species were previously unknown or poorly characterized and were identified in this study through indepth eco-spatially informed sampling throughout the spawning season and spawning depth range that was known for C. zugensis and C.
sp. "Bodenbalchen." Along the sampled spawning depth gradient, we find three distinct whitefish species structured in a phenotypic, likely adaptive (growth rate and gill-raker counts), and neutral genotypic cline. This depthhabitat-associated cline did not consist of a continuous phenotypic and genotypic gradation as the frequency distribution of SL of three- year-old fish was trimodal. Sampling bias caused by the different gill net mesh sizes alone cannot explain the absence of fish of specific body size required to generate a unimodal standard length distribution along the spawning depth gradient. Also, the age structure of whitefish provided evidence for more than one whitefish species cooccurring along the depth gradient with small SL class whitefish older F I G U R E 4 Likelihood of assignment of individuals from sampled whitefish species to allopatric genetic clusters (1) versus a Lucerne cluster (0). Pairwise DAPC comparisons were made using Lake Lucerne whitefish and whitefish from (a) Lake Brienz, (b) Lake Constance, (c) Lake Neuchatel, (d) Lake Thun, (e) Lake Walen, (f) Lake Zurich. Box plots represent the interquartile range around the median assignment probability for each taxon to the non-Lucerne cluster in each pairwise lake comparison than the larger intermediate SL class fish, suggesting the existence of at least two whitefish species differing strongly in their growth rates.
Analysis of the structuring of neutral genetic variation along the depth gradient provided the strongest evidence for multiple coexisting whitefish species spawning in close proximity. Grouping whitefish into three distinct classes according to the trimodal standard length distribution explained the highest amount of genetic variation. Utilizing the SL classes as distinct population units, Large and Small SL class whitefish could be assigned to shallow-spawning C. sp. "Bodenbalchen" and deep-spawning C. zugensis, respectively. While phenotypic and genetic differentiation was greatest between C. zugensis and C. sp.
"Bodenbalchen" at the extremes of the depth gradient, both species were significantly differentiated in allele frequencies and gill-raker

| Evidence for ecological speciation along an environmental gradient
If the current eco-spatial structuring of whitefish diversity along the spawning depth gradient is in any way indicative of the diversity present within the ancestral population at the onset of speciation, this may have facilitated divergence. Here, any significant association of ancestral genotypes with aspects of the environment along the depth gradient could bias mating encounter rates among individuals with different trait values, facilitating the origin of linkage disequilibrium between unlinked genes coding for these traits, and hence enable the origin of divergently adapted and reproductively isolated species through clinal speciation (Gavrilets, 2004). One aspect of our study system that is not captured by classical models of clinal speciation is that the observed depth cline appears only to be limited to the spawning season. A key aspect of many clinal speciation models is that dispersal is limited along the environmental gradient, decreasing gene flow and increasing the fitness consequences of localized selection regimes (Doebeli & Dieckmann, 2003;Heinz, Mazzucco, & Dieckmann, 2009;Kawata et al., 2007). Yet, whitefish outside the spawning season disperse over much greater distances than the geographic extent of the spawning gradient and indeed this study found co-occurrence of multiple species within certain depths during the spawning season.
How has reproductive isolation evolved/been maintained along the Lucerne depth gradient, given the potentially high levels of dispersal?
One answer to this may be that dispersal is not random with respect to the phenotype of the individuals and their chosen environment through matching habitat choice (Edelaar, Siepielski, & Clobert, 2008).
With increasing depth in lacustrine environments, abiotic components of the environment such as temperature, oxygen concentration, light intensity and spectral composition change rapidly alongside concomitant changes in biotic components (Seehausen et al., 2008).
Coupled changes in these factors will create localized divergent selection regimes at different locations along lacustrine depth gradients.
Matching habitat choice would result in spatial clustering of phenotypically similar individuals along these environmental gradients and if spawning location/timing were also governed by these phenotypefitness interactions, matching habitat choice would facilitate assortative mating and the beginnings of reproductive isolation. Support for this scenario comes from parallel patterns of spawning habitat divergence in multiple independent lake flocks within the Alpine whitefish radiation: small, densely gill-rakered zooplanktivores spawning deeper than coexisting larger, sparsely gill-rakered, zoobenthos feeding ecomorphs (Steinmann, 1950;Vonlanthen et al., 2012). The ecological factors driving predictable patterns of spawning depth/time segregation among coexisting whitefish ecomorphs remain an understudied aspect of speciation and adaptive radiation in this species complex.
While speciation from an eco-spatially structured ancestral population is theoretically less constrained than from one approaching panmixis, strong disruptive selection regimes and/or strong assortative mating would likely still be required to drive population divergence, especially given the high potential for gene flow (Doebeli & Dieckmann, 2003;Gavrilets, 2004;Nosil, 2012). Indeed, several lines of evidence support an ecological speciation scenario via the action of divergent natural selection for the origin of Lake Lucerne whitefish: (i) our study confirms previous works (Douglas, Brunner, & Bernatchez, 1999;Hudson et al., 2011) in showing that similar ecomorphs have arisen repeatedly across different lake flocks in the Alpine whitefish radiation. (ii) Significant structuring of neutral genetic and adaptive phenotypic variation was observed along an ecological gradient (water depth), whereas nonsignificant isolation-by-distance was uncovered between spatially distant but phenotypically similar populations within Lake Lucerne, strong evidence for the primacy of selective rather than neutral processes in the generation of intralacustrine whitefish diversity. (iii) Species diverging along the depth gradient show significant differences in traits related to niche utilization: gill-raker number and growth rate. Gill-rakers are projections on the gill arches thought to act as cross-flow filters to improve prey handling and retention (Sanderson, Cheer, Goodrich, Graziano, & Callan, 2001). Higher densities of gill-rakers along the gill arch have been shown to increase feeding efficiency on zooplankton in Alpine whitefish (Roesch, Lundsgaard-Hansen, Vonlanthen, Taverna, & Seehausen, 2013). Alpine whitefish with lower gill-raker densities, on the other hand, have been shown to be more efficient at foraging for large benthic prey items (Lundsgaard-Hansen, Matthews, Vonlanthen, Taverna, & Seehausen, 2013). These predictable differences in feeding efficiency are suggestive of fitness trade-offs in gill-raker number between the respective niches of coexisting Alpine whitefish species, corroborated by genetic evidence that within-lake patterns of gill-raker count variation are driven by divergent natural selection regimes in the Alpine and other whitefish radiations (Hudson et al., 2013;Praebel et al., 2013;Rogers & Bernatchez, 2007;Vonlanthen et al., 2009). Whitefish growth rate is a complex physiological trait impinging on many other aspects of the overall phenotype such as body shape. Interspecific differences in growth rate among whitefish species have been shown to be heritable Rogers & Bernatchez, 2007), and differences in growth rate have been shown to affect foraging ability in different aquatic niches in Alpine whitefish, where larger bodied whitefish are more efficient at exploiting benthic food resources . Additional roles for growth rate in whitefish adaptive divergence and speciation may include escaping size-related predation windows differing among niches (Kahilainen & Lehtonen, 2002) and size assortative mate choice (Foote & Larkin, 1988). Again, genetic evidence supports the role of divergent natural selection in driving patterns of interspecific variation in growth rate in different whitefish radiations (Rogers & Bernatchez, 2007;Vonlanthen et al., 2009). While both gill-raker number and growth rate are clearly ecologically important traits involved in adaptive divergence and potentially speciation, the frequency distributions of adult body size, indicative of growth rate, better reflected individual genetic differences among sampled depth gradient whitefish than gill-raker counts.
Further work is required to establish the exact trophic differences among the three species caught along the gradient, but body size may interact with gill-raker counts by changing the overall size of the gill arch, altering gill-raker density irrespective of the actual gill-raker number: heightening species' foraging performance in their respective habitats. (iv) Despite the observed demographic admixture (as opposed to obvious genetic admixture) between the three whitefish species on their spawning grounds, low but significant genetic differentiation between them is maintained, suggesting that these whitefish species are reproductively isolated. Overall, the above observations are consistent with isolation-by-adaptation and ecological speciation along an environmental gradient (Doebeli & Dieckmann, 2003;Endler, 1977;Nosil, 2012;Schluter, 2009).

| Anthropogenic impacts on whitefish diversity: implications for management
Two interacting anthropogenic impacts may have had profound effects on shaping endemic Lake Lucerne whitefish diversity: (i) eutrophication and (ii) stocking, both of allochthonous and autochthonous whitefish. Given the relative evolutionary youth of the Alpine whitefish radiation as a whole, mating barriers between the constituent species will be largely based on behavioural and/or extrinsic pre-and postzygotic isolating mechanisms rather than intrinsic genomic incompatibilities (Eckmann, 2015;Vonlanthen et al., 2012;Woods et al., 2009). In some Swiss lakes, anthropogenic eutrophication has led to a complete collapse of whitefish diversity over the past 60 years due to an interaction between negative population growth and speciation reversal . While levels of eutrophication were not as severe in Lake Lucerne, nutrient enrichment at its peak was likely high enough to impact the available spawning depth range and therefore potentially increase overlap and gene flow among spawning whitefish species, potentially contributing to the relatively low levels of genetic divergence currently seen among the depth gradient whitefish species . One potential outcome of human-mediated gene flow is the origin of C. sp. "benthic intermediate" through homoploid hybrid speciation via secondary contact and introgression between C. zugensis and C. sp. "Bodenbalchen" (Mavárez & Linares, 2008 (Kottelat & Freyhof, 2007). Alongside phenotypic similarities: DAPC, microsatellite-based trees and shared private allele analyses support a genetic affinity between C. sp. "pelagic intermediate" and C. wartmanni individuals from Lake Constance.
Historical records also exist of whitefish fry from Lake Constance being stocked into Lake Lucerne (Steinmann, 1950). This evidence is consistent with C. wartmanni being a progenitor of C. sp. "pelagic in- This has resulted in hybridization between the invasive species and native whitefish in certain lakes (Kahilainen et al., 2011). In other lakes, competitive exclusion of native densely rakered whitefish species by the invasive has occurred, combined with massive introgression and the loss of phenotypic and genetic differentiation between the previously distinct densely rakered and large sparsely rakered whitefish species (Bhat et al., 2014). Alternatively, neutral models of coexistence (Leibold & McPeek, 2006), potentially followed by ecological character displacement or specialization within the shared niche, buttressed by the nonoverlap of spawning habitats between these species could allow the persistence of both species. While increased sampling of C.
sp. "benthic intermediate" and C. sp. "pelagic intermediate" populations across Lake Lucerne would be helpful to draw exact inferences about their origins, ecology and co-distribution within the lake, further research on these species within the context of the local radiation may provide new insight into the factors that determine the extent of adaptive radiations and the hypothesis of diversity-dependent slowdown of speciation rates (Schluter, 2000).
Given the fine-scale partitioning of the spawning depth habitat in Lake Lucerne, a more insidious threat to intralacustrine whitefish diversity might be the supplemental stocking of native whitefish. Carried out during the spawning season, ripe adult fish are caught on their spawning grounds, their gametes stripped and mixed, and the resulting fertilized eggs are raised to fingerlings before being released. This practice is carried out in many Swiss lakes, and these stocked whitefish may comprise a large proportion of a species' population (e.g., up to 83% of C. wartmanni in Lake Constance; Eckmann, 2012), despite a lack of evidence for any resultant long-term yield benefits especially in healthy lacustrine ecosystems. Alongside impairing natural recruitment, especially in nontarget whitefish species, and imposing inadvertent artificial sexual selection on target species (Eckmann, 2012), supplemental stocking may increase rates of introgressive hybridization as demographically mixed catches of individuals from different species but of broadly similar appearance (e.g., older/larger C. zugensis, younger/smaller C. sp. "benthic intermediate") are accidentally treated as a single stock for crossing. The increase in observed intralacustrine complexity highlighted in our study will require the development of more sophisticated fisheries and conservation management plans to maintain whitefish diversity for stakeholders and posterity (Lankau, Jørgensen, Harris, & Sih, 2011;Lapointe et al., 2014).

| CONCLUSIONS
Through spatially and ecologically informed sampling, our study has revealed that currently Lake Lucerne holds at least six distinct whitefish species, comprising one of the most diverse intralacustrine species flocks within the Alpine whitefish radiation and the wider C. lavaretus species complex globally. Along a depth gradient during winter-spawning time, we find a previously unknown species: C. sp.
"benthic intermediate", currently distinct in gill-raker counts, growth rate and allele frequencies from both known whitefish species that spawn in close proximity. The nondiscrete nature of phenotypic variation among whitefish species spawning along the depth gradient makes C. sp. "benthic intermediate" vulnerable to increased gene flow via fisheries management techniques such as supplemental stocking, where the true extent of coexisting whitefish diversity may be underestimated. Moreover, the tight spatial packing of species on the spawning gradient, coupled to the dependence of spawning niche segregation on the persistence of fine-scale depth-related differences in the lacustrine environment, makes these species especially vulnerable to increased gene flow via changes in the physicochemical habitat characteristics, such as those associated with changes in primary productivity. A sixth distinct species, C. sp. "pelagic intermediate," differing in growth rate and reproductive ecology from C. sp. "benthic intermediate," may be wholly or partially of non-native origin, the result of whitefish stocking between lakes and further evidence of anthropogenic impacts on Lucerne whitefish biodiversity.