Rapid buildup of sympatric species diversity in Alpine whitefish

Abstract Adaptive radiations in postglacial fish offer excellent settings to study the evolutionary mechanisms involved in the rapid buildup of sympatric species diversity from a single lineage. Here, we address this by exploring the genetic and ecological structure of the largest Alpine whitefish radiation known, that of Lakes Brienz and Thun, using microsatellite data of more than 2000 whitefish caught during extensive species‐targeted and habitat‐randomized fishing campaigns. We find six strongly genetically and ecologically differentiated species, four of which occur in both lakes, and one of which was previously unknown. These four exhibit clines of genetic differentiation that are paralleled in clines of eco‐morphological and reproductive niche differentiation, consistent with models of sympatric ecological speciation along environmental gradients. In Lake Thun, we find two additional species, a profundal specialist and a species introduced in the 1930s from another Alpine whitefish radiation. Strong genetic differentiation between this introduced species and all native species of Lake Thun suggests that reproductive isolation can evolve among allopatric whitefish species within 15,000 years and persist in secondary sympatry. Consistent with speciation theory, we find stronger correlations between genetic and ecological differentiation for sympatrically than for allopatrically evolved species.

Appendix S4. Results. Qualitative differences in depth range among sympatric species. In Lake Brienz ( Figure S7a), all C. sp. "Balchen1" were caught in the shallow littoral zone (< 11m). C. sp. "Balchen2" was absent from the littoral zone but occurred at all depths in all other lake habitats. C. sp. "Felchen" and C. albellus were common in all lake zones. In Lake Thun ( Figure S7b), all C. sp. "Balchen1" and C. sp. "Albock" were caught at rather shallow depths (max. 27m and 34m, respectively) in both littoral and pelagic zones. C. alpinus was caught at all depths of the profundal benthic zone (15-207m). C. sp. "Balchen2", C. sp. "Felchen" and C. albellus were caught at all depths, but most rarely in the littoral zone.
Appendix S5. Results. Morphological and reproductive differences in the littoral species between lakes. We found that the littoral spawning species "Balchen1" spawns earlier in Lake Brienz than in Lake Thun (Figure 1). This was already the case 125 years ago, as reported by Fatio (Fatio, 1890). One reason could be that the cold water temperatures, which are necessary for egg development of whitefish, are reached earlier in colder Lake Brienz than in Lake Thun. Such temporal reproductive isolation may have allowed these two populations to diverge even before the lakes were completely disconnected, and may have facilitated the evolution of weak differentiation in GRN. Another factor that could contribute to the between-lake divergence of this species is hatchery breeding of this and other species in Lake Thun, which occurred to a lesser extent in Lake Brienz ( Figure S15). This could cause stronger admixture (despite efforts of fisheries mangers to breed "Balchen" and "Albock" separately), fisheriesinduced selection and selection for hatchery-adapted genotypes in Lake Thun. Finally, "Balchen1" is the rarest whitefish species in both lakes (Figure 1), and its presumably small effective population size may lead to faster divergence by drift since the separation of the two lakes than for species with higher abundance and larger effective population sizes. To assess the contribution of stocked fish to the whitefish community in the lake and to estimate the degree of admixture attributable to the fertilization process in the hatchery, markrecapture studies could be conducted.
To further assess the introduction of whitefish from Lake Constance into Lake Thun, we identified private alleles among lake Constance whitefish species using microsatellite data from individuals collected during the pre-eutrophication period of this lake , and we explored the occurrence of these alleles in Lake Thun. For our question, we considered those private alleles from Lake Constance species to be informative that were found in Thun, but not in Brienz. This is because if one of the Constance private alleles is found in Thun and in Brienz, it is possible that that the allele independently evolved within this lake system, rather than being derived from stocking. Given that Lake Brienz was also at least once stocked with whitefish from Lake Constance (Heuscher, 1901), although to a much lesser extent than Lake Thun, it is possible that private alleles could have gotten into both lakes through stocking. Our test is conservative in this regard.
For each species from Lake Constance, we calculated the proportion of private alleles that fulfilled this criterion (present in Thun, absent in Brienz). We found one of two private alleles for C. wartmanni (50%), four of 16 for C. macrophthalmus (25%) and one of eight for C. gutturosus (12.5%) that fulfilled this criterion (Table S19) (No private alleles were found for C. arenicolus). These patterns are consistent with the introduction of any or all of the pelagic whitefish species from Lake Constance into Lake Thun, and thus does not identify a single most likely introduced species.
It is possible that the founder event associated with the introduction and/or admixture with native species of Lake Thun changed allele frequencies of the introduced species, which makes inference of the original species identity difficult when using few microsatellites. Furthermore, misidentification of introduced whitefish from Lake Constance is possible, as C. wartmanni and C. macrophthalmus are both pelagic fishes, and are phenotypically similar and both spawn in winter. Hence, we cannot conclusively say which species today's "Albock" from Lake Thun derives from, but both historical reports and our genetic data are most consistent with it being either one or both of these pelagic species from Lake Constance. More detailed genomic analyses would be necessary to resolve which species from Lake Constance and also from Lake Thun contribute to today's "Albock" from Lake Thun.
Appendix S7. Discussion. Taxonomic considerations Three of the four species that we found in this study in Lake Brienz were previously known: one is taxonomically described (C. albellus Fatio 1890; revisions by Kottelat, 1997) and two are historically documented (Fatio, 1890) and well-known by local fishermen (C. sp. "Felchen" and C. sp. "Balchen"). Because the fourth, previously unknown species is genetically and phenotypically similar to the known C. sp. "Felchen" and C. sp. "Balchen", we call it here C. sp. "Balchen2". To our knowledge, genetic substructure within "Balchen" has not been previously reported. But remarkably, Fatio (1890) and Steinmann (1950) already mentioned two different types of C. fatioi, of which one might correspond to one of the "Balchen" clusters from this study.
C. alpinus Fatio 1885 ("Kropfer"), is a native endemic to Lake Thun and is taxonomically described (Fatio, 1885;revisions by Kottelat, 1997), whereas C. sp. "Albock", was introduced from Lake Constance, and is not taxonomically described. Note that despite its local name "Albock" being the same as that of native C. fatioi, it is certainly not this native species. It is unclear whether the taxonomically described C. fatioi Kottelat 1997(Kottelat, 1997 corresponds to C. sp. "Balchen 2", C. sp. "Felchen" or something else (e.g. the whitefish species mentioned by Fatio (1890) that historically migrated between the two lakes, whose migration is now impossible due to the completion of the water gates in 1856.) In Lake Brienz, nine individuals (7 contemporary, 2 historical scale samples) had major genetic assignment proportions for C. alpinus, the profundal whitefish species only known from Lake Thun ( Figure S9). However, all seven contemporary genetic C. alpinus fish from Lake Brienz were sampled in shallow waters ( Figure S6) and they had very high . This is both very atypical based on what is known from this species from Lake Thun, and we consider it unlikely that these fish are actually corresponding to C. alpinus. Such a mismatch could result from extensive introgression of C. alpinus into remaining whitefish species of Lake Brienz, as has been suggested for the very similar case in Lake Constance where the profundal species (C. gutturosus) is extinct today, and the few individuals from contemporary samplings that are genetically assigned to C. gutturosus do not phenotypically correspond to what is known from this species .
The reasons for the absence of C. alpinus in Lake Brienz despite its natural presence in Lake Thun and the long shared history of the lakes can only be speculated. Higher turbidity, lower temperatures and steeper bathymetric slopes resulting in lower amounts of benthic habitat of intermediate depth in Brienz may not provide the niche requirements of the profundal benthic species, so that it never occurred there. Alternatively, C. alpinus might have been originally present in Lake Brienz, but went extinct after the separation of the two lakes.
Appendix S8. Discussion. Challenges and limitations to assessing sympatric whitefish species diversity. Both the choice of genetic markers and sampling design have the potential to bias the number of species and the genetic structure among species we can recover in our dataset.
The six species we identified in this study should be considered as a minimum estimate for the actual whitefish species diversity present in Lakes Brienz and Thun. With 10 neutral markers and no prior grouping based on phenotypes, we are limited to detect common and clearly reproductively isolated groups, and we miss rare and/or only weakly differentiated whitefish species. For example, while RADseq data show very clear separation of three whitefish species from Lakes Walen and Zürich (max. F ST =0.11, min. F ST =0.03, unpublished data, Feulner et al.), our 10 microsatellites could not resolve any of these species using the same samples (N=20 per species) with the program Structure. However, F ST s among groups were the same for both types of data (data not shown). We therefore suggest combining genomic data with ecological and morphological data to identify such rare populations and to resolve genetic structure among them.
It is well known that discontinuous sampling of a continuous distribution can produce discrete genetic clusters in structure analyses (e.g. Serre, 2004;Rosenberg et al., 2005). In our study, we can reproduce a more discontinuous distribution of genetic variation when considering only samples from targeted spawning fishing, for example, using whitefish from Lake Brienz from Bittner (2009) (Figure S9). Vice versa, when relying on samples from random, quantitative fishing alone, the unequal abundance of genetic variation makes detection of population structure using clustering programs such as Structure (Pritchard et al., 2000) very difficult. Furthermore, some genotype combinations that are not uncommon in the targeted sampling are completely absent from the random sampling (e.g. genetic intermediates between C. sp. "Balchen1" and C. sp. "Balchen2"( Figure S9)), which gives a misleading picture of the actual structure in the whitefish community. Hence, the two types of samplings are complementary, and both are needed to assess whitefish diversity and its genetic structure in these lakes.
Appendix S9. Discussion. Multiple dimensions of niche partitioning contribute to RI and coexistence in sympatry. In both lakes, we found significant differentiation in spawning depth and/or spawning time among all native species (Figure 1, Table S7, S8), whereby the former was generally stronger and predicted the degree of genetic differentiation better than the latter. This suggests a greater importance of spawning depth than of spawning time for maintaining RI among sympatric Alpine whitefish species, consistent with studies of other Alpine lakes Ingram et al., 2012;Hudson et al., 2016). Spawning depth differentiation also plays a major role for RI in other young species radiations of fish, such as haplochromine cichlid fishes in Lake Victoria (Seehausen et al., 2008) and benthic-limnetic stickleback species pairs in Canada (Hatfield and Ptolemy, 2001). The parallelism of the two dimensions of reproductive niche (Figure 1) likely strengthens total RI among species, as indicated by comparisons to Alpine whitefish radiations without spawning time segregation. In Lakes Thun and Brienz, the four species from the spawning depth gradient do not only spawn at different depth, but also at different times of the year, whereas the less strongly differentiated species from Lake Lucerne's and Neuchâtel's spawning depth gradient all spawn in winter and do not differ in spawning time anymore (Hudson et al., 2016;Vonlanthen et al., 2009). Historically, stronger temporal spawning segregation also occurred among Lake Lucerne whitefish, with C. zugensis spawning either in summer or in January and C. nobilis in summer (Steinmann, 1950). Lake eutrophication, which was stronger in Lakes Lucerne and Neuchâtel than in Thun and Brienz, caused low oxygen concentrations at great depth in summer, which could have favored individuals that spawn later in the year, which increased the chances for admixture with winter spawning species. Eutrophication could therefore have increased gene flow among species by contracting both the spatial and the temporal spawning gradient.
Our finding that differentiation in gill raker numbers explains residual genetic differentiation not explained by spawning depth differentiation ( Figure 3e) could indicate that pre-zygotic spawning segregation and post-zygotic divergent natural selection complement each other to maintain RI among sympatric whitefish species.
Appendix S10. Discussion. More interspecific gene flow in Lake Thun than in Lake Brienz. Although species are fewer in Lake Brienz than in Lake Thun, genetic distinctiveness is higher among species of Lake Brienz. We find several lines of evidence that suggest higher rates of gene flow between species in Lake Thun than in Lake Brienz. Lake Thun has more individuals with intermediate assignment likelihoods (Table S6), the four species that occur in both lakes show weaker genetic differentiation between them in Lake Thun (Table 1), and in genetic PCA, species of Lake Thun show more overlap in genotype space than species of Lake Brienz ( Figure 1). Finally, the relationship between genetic assignment and GRN between closely related species is weaker in Lake Thun than in Lake Brienz ( Figure S3). One explanation for higher rates of interspecific gene flow in Lake Thun could be that higher species numbers lead to more opportunities for hybridization. However, when considering individuals having minor genetic contributions from the two species unique to Lake Thun, C. alpinus and C. sp. "Albock" (<0.2), F ST s still remain clearly lower in Lake Thun than in Brienz (Table S20). Alternative explanations for weaker RI are stronger anthropogenic influences in Lake Thun than in Brienz (e.g. eutrophication and stocking, Vonlanthen et al., 2012; Figure S16). Disruptive selection could also just naturally be stronger in Lake Brienz. Finally, if effective population sizes of species were consistently lower in Lake Brienz than in Lake Thun, we would expect the dynamics of incomplete lineage sorting to lead to shorter coalescent times and the appearance of fewer intermediate genotypes among species of Lake Brienz.
Appendix S11. Discussion. Increased abundance of the introduced species. A comparison between genetic data from the 1950s-1970s and from 2004-2014 indicates that the importance of the introduced species changed over time for professional fisheries. Genetic data from the 1970s suggest that the introduced species was absent in catches of professional fishermen in Lake Thun ( Figure S12). Contrary, after 2000, it was common in spawning fishery catches of commercial fishermen (qualitative data only) ( Figure S14) and it made up 8.4% of individual whitefish in the catches of habitat stratified random fishing ( Figure S5, S14). Whether the introduced species indeed increased in abundance only during the last third of the 20 th century, or whether it was already common before, but not present in our samples from 1950-1970, is not clear from our data. However, if the former was the case, it could be hypothesized that its rise in abundance was mediated by the peak of mild eutrophication in Lake Thun in the 70s/80s. Appendix S12. Discussion. The buildup of local whitefish species diversity is limited by constraints to speciation The establishment and persistence of an introduced whitefish species in Lake Thun (and also in Lake Lucerne, Hudson et al., 2016) suggests that unsaturated niche space exists even in large whitefish radiations, and that carrying capacity for whitefish species richness might not be reached over postglacial times by intralacustrine speciation alone. Hence, speciation seems to be the limiting step for achieving high species richness in young adaptive radiations of whitefish. The extent of the ecological gradient likely determines the number of whitefish species that evolve sympatrically within a lake within a given time . The buildup of local species richness through secondary contact between allopatrically evolved species, on the other hand, is limited by a lake's geographical isolation, and at a larger scale by the number of isolated lakes available for allopatric speciation. Hence, if isolation among lakes is strong and ecological opportunity within lakes large, the majority of species in each lake may still derive from sympatric speciation, as seems to be the case in Alpine whitefish (this study, Hudson et al., 2016;Hudson et al., 2011) and haplochromine cichlids (Wagner et al., 2014). If connectivity among geographically isolated areas then suddenly increases after a period of isolation, be it due to human-mediated movements (whitefish this study, Hudson et al., 2016; cichlids from islands Young et al., 2009) or natural changes in lake levels (e.g. Lake Tanganyika), species richness in local communities may drastically increases. Figure S1. Schematic overview of the analysis workflow of the hierarchical structure analysis (blue) and the subsequent reference based assignment analysis (light blue). The former was used to identify genetic clusters in the full dataset, the latter to assign all individuals to the clusters identified in the hierarchical analysis. For each data subset, sample size (N) and the most likely K value are given. Hierarchical steps are labeled with circled numbers above arrows corresponding to the splits in the dataset. Solid arrows (step 1,2,3) indicate that the most likely K was >1. Dashed arrows (step 4 and 5) indicate similarly high likelihood for K=2 and K=1 and associations between genetic and ecological structure in both lakes, suggesting the presence of two biologically meaningful genetic groups within these genetic clusters (C. sp. "Felchen" and C. albellus, and C. sp. "Balchen 1" and C. sp. "Balchen 2", see Figure S2, S3). The bottom structure plot shows the result of the assignment analysis, wherein individuals are sorted by species (based on maximum assignment) and decreasing assignment likelihood within each species.

Supplementary Figures
Analysis description: To find the most likely value of K, we conducted a hierarchical cluster analysis ( Figure S1; Coulon et al., 2008) using the individual-based Bayesian clustering algorithm implemented in SRUCTURE (Pritchard et al., 2000). We first determined the most likely number of clusters (K) for the full dataset, then the most likely K within each of the data subsets suggested by the initial analysis, and so forth until all subsets supported a value of K=1. To determine the most likely K at each hierarchical level, we first compared LnP(D) values from runs of different Ks, as suggested by Pritchard and Wen (2003). If this method suggested K>1 was most likely, we determined the most likely K using the delta K method of Evanno et al. (2005) (note that the Evanno et al. method cannot evaluate K=1, hence this requires to also evaluate LnP(D) values). We ran Structure for K=1 to K=15 at each hierarchical level, with 10 iterations at each K value, using 100'000 burn-in steps followed by 500'000 MCMC steps to generate the posterior sample distribution. For all runs, we used the admixture and correlated allele-frequency model. We subdivided the original dataset into the K subsets suggested by the most likely K by assigning individuals to the cluster for which they had the highest assignment likelihood in the run with the highest likelihood for this K.
To determine correspondence of genetic clusters to known species, we assessed how individuals from samplings targeted to known species were distributed amongst the clusters.
Because LnP(D) of K=2 and K=1 were very similar for the cluster from the hierarchical analysis corresponding to "Balchen" (Figure S2), we further explored the distribution of assignment likelihoods within this cluster for K=2. For each lake, we tested whether the two "Balchen" clusters (individuals were assigned based on their maximum assignment proportion obtained from Structure analysis on this entire cluster for K=2) differed in GRN and spawning depth using Wilcoxon tests, and we assessed the relationship between Structure assignment likelihoods to "Balchen" cluster 1 and GRN in spearman's rank correlation tests.
The cluster derived from the hierarchical approach containing individuals from known spawning sites of C. sp. "Felchen" and C. albellus also showed very similar values of LnP(D) for K=2 and K=1 ( Figure S2). Therefore, we explored genetic structure within this cluster in more detail. To test whether unequal sample sizes of C. sp. "Felchen" and C. albellus reduced Structure's ability to distinguish the two groups, a well-known limitation of this program (Puechmaille, 2016), we performed a Structure analysis using equal numbers of mature adult individuals from three C. albellus and three C. sp. "Felchen" spawning sites (N=84 for each species) from Lake Brienz (all individuals from these 6 spawning sites were included). We ran Structure for K=1-4 with 3 replicates per K, with 100'000 burn-in steps and 50'000 MCMC steps using the admixture and correlated allele frequency model. Additionally, we assessed the relationship between Structure assignment likelihoods to C. sp. "Felchen" (obtained from Structure analysis on the entire C.albellus-C. sp. "Felchen" cluster for K=2) and GRN in spearman rank correlation tests within each lake. To guard against the possibility that introgression from the low GRN species C. alpinus confounded a potential relationship in Lake Thun (C. alpinus is absent from Lake Brienz), we performed this analysis by including only individuals with more than 30 gill rakers. Furthermore, we tested whether the two groups (individuals were assigned based on their maximum assignment proportion) differed in GRN and spawning depth using Wilcoxon tests.
To obtain for each individual genetic assignment proportions to the clusters inferred in the hierarchical analysis, we performed Structure assignment analyses using the clusters identified by the hierarchical analysis as reference populations ( Figure S1). Genotypic data from 50 representative individuals from each of the 6 clusters identified before (named according to the known and newly recognized species they contain: C. sp. "Felchen", C. albellus, C. sp. "Balchen1", C. sp. "Balchen2", C. alpinus and C. sp. "Albock") were used to define reference populations by setting PopFlag = 1 and the corresponding population number. We therefore chose the 50 individuals with highest assignment likelihood in the corresponding clusters at each previous step in the hierarchical analysis (see text in figure for detailed criteria). For all reference populations, we only used contemporary individuals with no missing data, except for the C. sp. "Balchen 1" and C. sp. "Balchen 2" references, where we also used data from the 1950s-1970s, since these groups are rare in all samples. All remaining whitefish in the dataset (N=2088) were then assigned to the reference populations by coding them with PopFlag = 0 and population = 0, and separate analyses with subsets of individuals to assign were performed. Each of these analyses consisted of the 6*50 reference individuals plus 50 individuals to be assigned. For each of these analyses, we performed 10 replicates of K=6. We used Structure Harvester to generate input files for CLUMPP (Jakobsson and Rosenberg 2007), which we used to generate consensus percentages of assignment proportions across each of the 10 structure runs. The maximum difference in assignment likelihoods among the 10 runs was less than 0.055 for any assigned individual, and was less than 0.05 among the 420 runs for the reference individuals, indicating consistency across runs.
In Lake Brienz, nine individuals (7 contemporary, 2 historical scale samples) had major genetic assignment proportions for C. alpinus, the profundal whitefish species only known from Lake Thun ( Figure S9). However, all seven contemporary genetic C. alpinus fish from Lake Brienz were sampled in shallow waters ( Figure S6) and they had very high . This is both very atypical based on what is known from this species from Lake Thun, and we consider it unlikely that these fish are actually corresponding to C. alpinus. Such a mismatch could result from extensive introgression of C. alpinus into remaining whitefish species of Lake Brienz, as has been suggested for the very similar case in Lake Constance where the profundal species (C. gutturosus) is extinct today, and the few individuals from contemporary samplings that are genetically assigned to C. gutturosus do not phenotypically correspond to what is known from this species . Figure S2. Mean ln probability ± standard deviations and delta K of 10 runs for K=1 to K=15 for the different clusters of the hierarchical Structure analysis. Numbers correspond to hierarchical steps labeled with circled numbers in Figure S1. Figure S3. Genetic and ecological structure within the "Balchen" cluster (a,b) and within the C. sp. "Felchen"-C. albellus cluster in Lakes Thun (a,c) and Brienz (b,d). a, b) top: Genetic assignment to C. sp. "Balchen 1" (obtained from Structure analysis for K=2) is significantly related to the number of gill rakers in both lakes. Histograms of gill raker numbers for the two groups are shown along the y-axis (individuals are assigned based on their maximum assignment proportion obtained from Structure for K=2). Structure plots for K=2 are shown along the x-axis and are sorted by lake and increasing assignment to C. sp. "Balchen 1". bottom: C. sp. "Balchen 1" spawns at significantly shallower depth than C. sp. "Balchen 2" in both lakes. c,d) top: Genetic assignment to C. sp. "Felchen" (obtained from Structure analysis for K=2) is significantly related to gill raker numbers in Lake Brienz, but not in Lake Thun. Histograms of gill raker numbers for the two groups are shown along the y-axis (individuals are assigned based on the maximum assignment proportion obtained from Structure for K=2). Structure plots for K=2 are shown along the x-axis and are sorted by lake and increasing assignment to C. sp. "Felchen". bottom: In Lake Brienz, C. albellus spawns at significantly greater depth than C. sp. "Felchen", whereas in Lake Thun, C. albellus does not significantly differ in spawning depth from C. sp. "Felchen". Figure S4. Structure analysis for K=2 with reproductively mature, ripe individuals from three C. albellus spawning sites and from three C. sp. "Felchen" spawning sites from Lake Brienz. ! Figure S5. Private allele analyses for all whitefish species from Lakes Thun and Brienz. For each species, the mean number of private alleles per locus is plotted as a function of standardized sample size. a,c,e) show the frequency of private alleles for contemporary whitefish considering all individuals (species assignment based on maximum assignment LH in the Structure assignment analysis), whereas b,d,f) show the private allele frquencies only for clearly assigned individuals (structure assignment LH > 0.7). Panels a and b are for individuals from both lakes combined, c and d for Lake Brienz, e and f for Lake Thun. Figure S6. Genotypic distribution of the whitefish community from the 1950-70s in Lake Thun (left) and in Lake Brienz (right). a) Tetrahedral plot show the genotypic distribution of the whitefish community in Lake Thun (left) and Lake Brienz from the 1950-70s (right). The location of an individual whitefish is determined by its Structure assignment proportions obtained from the analysis using reference populations: The corners of the tetrahedron correspond to 100% assignment to a cluster, intermediate genotypes lie within the space framed by the corners. An individual's color corresponds to its combination of assignment proportions for the different reference clusters. b) Frequency distributions of assignment proportions from the analysis using reference populations for all pairs of genetic clusters for whitefish from Lake Thun (left and two bottom panels) and Lake Brienz (right). For each plot, only individuals that are assigned to either one of the two genetic clusters under consideration or are genetically intermediate to these and not closer to any of the other genetic clusters are used. The frequency distribution of assignment proportion of the individuals fulfilling this criterion is plotted for one of the genetic clusters under consideration (the one indicated at position 1.0). Note that the comparisons with C. alpinus are only shown for Lake Thun, since this species are absent in Lake Brienz. Figure S7. Spatial distribution of whitefish species as revealed by habitat stratified random fishing in Lake Brienz (a) and Lake Thun (b). Dots represent individual whitefish, which are colored by their maximum assignment to the 6 reference populations. (C. albellus blue; C. sp. "Felchen" red, C. sp. "Balchen 1" black; C. sp. "Balchen 2" green; C. sp. Albock" purple; C. alpinus yellow). We used the R package "beeswarm" on groups of 5m lake depth intervals for these graphs, and the function "jitter" (with jitter amount of 0.2 and 2 for x-and y-axes, respectively) to increase scatter among individuals for better visualization. Black line represents the lake bottom. Figure S8. Comparisons between differentiation (P ST ) in different ecological (gill raker numbers vs. habitat depth, a) and reproductive niches (spawning depth and spawning time, b) in each lake. Individual points are pairwise comparisons between sympatric species, results of paired t-tests are indicated on top. Figure S9. Habitat stratified random sampling (left) and targeted spawning samplings (right) are sampling different parts of the genotypic whitefish diversity of Lake Brienz. Intermediate genotypes between C. sp. "Balchen 1", C. sp. "Balchen 2" and C. sp. "Felchen" are not present in the quantitative sampling, whereas backcrosses to C. sp. "Balchen 2" of C. albellus-C. sp. "Balchen 2" hybrids are lacking in the spawning sampling. Figure S10. Results from the Structure analysis using reference populations for whitefish from Lake Brienz. The top panel shows assignment results for all samples from Lake Brienz, the middle for contemporary samples, and the bottom for historical samples. Individuals are arranged by species (determined by maximum assignment likelihood), and sorted by decreasing assignment likelihood to the corresponding species. Figure S11. Results from the Structure analysis using reference populations shown for whitefish from Lake Thun. The top panel shows assignment results for all samples from Lake Thun, the middle for contemporary samples, and the bottom for historical samples. Individuals are arranged by species (determined by maximum assignment likelihood), and sorted by decreasing assignment likelihood to the corresponding species. Figure S12. Results from the Structure analysis using reference populations shown for whitefish from Lake Brienz caught during the habitat stratified random sampling and the targeted sampling to known species during spawning time. Individuals are arranged by species (determined by maximum assignment likelihood), and sorted by decreasing assignment likelihood to the corresponding species. Figure S13. Results from the Structure analysis using reference populations shown for whitefish from Lake Thun caught during the habitat stratified random sampling and the targeted sampling to known species during spawning time. Individuals are arranged by species (determined by maximum assignment likelihood), and sorted by decreasing assignment likelihood to the corresponding species. Figure S14. Reports from the local fishery inspector and fisheries guard documenting the introduction of Lake Constance whitefish into Lake Thun in 1935 (Archiv Oberländischer Fischereiverein, Interlaken, Switzerland).  Table S1. Known whitefish species diversity in Lakes Thun (T) and Brienz (B).  Table S1 in Vonlanthen et al. (2012), except if indicated differently: † Fatio (1890) † † Local fishermen and previous reports (Steinmann, 1950;Kirchhofer, 1990)  1360 † For samples from stratified random sampling which were all genotyped in this study (N=944), total DNA was extracted from muscle or fin tissue using Chelex and proteinase K, or a QIAGEN Bio Sprint 96 extraction robot following the manufacturer's standard protocol. PCR amplification was performed using the QIAGEN Mulitplex PCR Kit in a 10 µl reaction volume of 5 µl Multiplex Master Mix, 0.5 µl Primer Mix, 3.3 µl H2O and 1.2 µl DNA extraction product according to the following protocol: initial denaturation for 15 min at 95°C, 35 cycles of 30 sec at 94°C, 90 sec at 57°C, 90 sec at 72°C; final extension for 30 min at 72°C. † † All scale samples from the 1950s--1970s were genotyped with the methods outlined in Vonlanthen et al. (2012). Table S3. Genetic differentiation (F ST ) between clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). All F ST s are highly significant (p<0.001). Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz.  Table S4. Locus by locus Fst for all pairwise comparisons of species from Lake Brienz. Values in bold indicate significance (p<0.05).  Table S5. Locus by locus Fst for all pairwise comparisons of species from Lake Thun. Values in bold indicate significance (p<0.05).  Table S6. Proportion of unclearly assigned individuals (highest assignment proportion < 0.7) for each species and over all individuals caught during quantitative fishing in Lakes Thun and Brienz.

Genetic group
Thun Brienz  Table S7. P--values of Kruskal--Wallis and post--hoc dunns test for differences in spawning depth for clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). P--values are adjusted for multiple testing using Holm's method. Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz.  Table S8. P--values of Kruskal--Wallis and post--hoc dunns test for differences in spawning time for clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). P--values are adjusted for multiple testing using Holm's method. Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz. Brienz X 2 3=82.08, p=0 Table S9. P--values of Kruskal--Wallis and post--hoc dunns test for differences in GRN for clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). P--values are adjusted for multiple testing using Holm's method. Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz.  Table S10. P--values of Kruskal--Wallis and post--hoc dunns test for differences in habitat depth in autumn in the benthic zone for clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). P--values are adjusted for multiple testing using Holm's method. Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz. Tests for differentiation in depth occupation in autumn could not be performed for C. sp. "Balchen1" and C. sp. "Balchen2", because of very small sample sizes in each lake.

Genetic group
C. sp. "Balchen1" (4) C. sp. "Balchen2" (1) C. sp. "Felchen" (12) C. albellus ( Table S11. P--values of Kruskal--Wallis and post--hoc dunns test for differences in habitat depth in autumn in the pelagic zone for clearly assigned groups (Structure assignment >0.7) within Lake Thun (below diagonal) and within Lake Brienz (above diagonal). P--values are adjusted for multiple testing using Holm's method. Sample sizes are given in the first column for Lake Thun, and in the first row for Lake Brienz. Tests for differentiation in depth occupation in autumn could not be performed for C. sp. "Balchen1" and C. sp. "Balchen2", because of very small sample sizes in each lake.

Genetic group
C. sp. "Balchen1" (2) C. sp. "Balchen2" (1) C. sp. "Felchen" (10) C. albellus ( Table S12. Mantel and partial Mantel correlations between genetic cluster membership and geographic distance, while correcting for spawning depth or/and spawning time. We only used contemporary individuals. For Lake Thun, we either considered individuals from all species or only those whose sum of assignment likelihoods for the four species also found in Lake Brienz was >0.85 Mantel or Partial mantel test between Thun all Thun 4 species Brienz cluster membership geographic distance 0.07 *** 0.04 *** 0.12 *** cluster membership geographic distance + spawning depth 0.04 *** 0.01 ns 0 ns cluster membership geographic distance + spawning time 0.06 *** 0.02 ** 0.12*** Residuals(cluster membership spawning depth)~ geographic distance + spawning time 0.03 *** ---- Table S13. Intraspecific Mantel correlations between ecological and genetic data. Pearson correlation (r) and P--values obtained from 10000 permutations are reported. Only clearly assigned individuals were used for comparisons (maximum assignment proportion >0.7). For "Balchen1" and "Balchen2", sample sizes of the fishing outside of spawning season are too low to compare habitat depth occupation.  Table S14. Between lake comparisons within species. P--values obtained from Arlequin (F ST ) or from Wilcoxon's signed rank tests are reported (adjusted for mutliple testing). Only clearly assigned individuals were used for comparisons (maximum assignment proportion >0.7). For "Balchen1" and "Balchen2", sample sizes of the fishing outside of spawning season are too low to compare habitat depth occupation.  Table S15. Multilocus pairwise F ST (8 markers) between C. sp. "Albock" from Lake Thun and the four whitefish species from Lake Constance. Data are from the pre--eutrophication period for species from Lake Constance. All F ST s are highly significant (p<0.001).   Table S17. Correlation coefficients for Manteltests between individual genetic distance and spawning depth (SD), geographic distance of spawning locations (Geo), spawning time (ST) or gill raker numbers (GRN), and partial Manteltests correcting for one of these factors. We only used contemporary individuals and for Lake Thun, only those whose sum of assignment likelihoods for the four species also found in Lake Brienz was >0.85. P--values are adjusted for multiple testing.  Table S18. Correlation coefficients for partial Manteltests between individual genetic distance and genetic cluster membership, while correcting for spawning depth, geographic distance of spawning locations, gill raker numbers or all of these factors. We only used contemporary individuals and for Lake Thun, only those whose sum of assignment likelihoods for the four species also found in Lake Brienz was >0.85.

Partial mantel test between Thun Brienz
Genetic distance cluster membership + spawning depth 0.28 *** 0.20 ***  Table S19. Number of private alleles of Lake Constance whitefish species prior to eutrophication (8 markers), and their representation in Lake Thun. Numbers in parentheses indicate the number of markers in which private alleles were found. Species specific alleles from Lake Constance species found in Lake Thun, but not in Lake Brienz 1 4 (2) 1 -- Table S20. Genetic differentiation (F ST ) between genetic groups (individuals are assigned based on their maximum assignment proportion) within Lake Thun for individuals with minor assignment to C.alpinus or "Albock" (sum of structure assignment LH to "Balchen1", "Balchen2", "Felchen", C. albellus >0.85). All F ST s are highly significant (p<0.001).