α - and β -Diversity Patterns of Macrophytes and Freshwater Fishes are Driven by Di ﬀ erent Factors and Processes in Lakes of the Unexplored Southern Balkan Biodiversity Hotspot

: Disentangling the main drivers of species richness and community composition is a central theme in ecology. Freshwater biodiversity patterns have been poorly explored; yet, it has been shown that di ﬀ erent freshwater biota have di ﬀ erent, often contrasting responses to environmental gradients. In this study, we investigated the relative contribution of geographical and environmental (habitat-, climate- and water quality-related) factors / gradients in shaping the α - and β -diversity patterns of macrophytes and ﬁsh in sixteen natural freshwater lakes of an unexplored Balkan biodiversity hotspot, the Southern Balkan Peninsula. We employed generalized linear modeling to identify drivers of α -diversity, and generalized dissimilarity modeling to explore commonalities and dissimilarities of among-biota β -diversity. Species richness of both biota was signiﬁcantly associated with lake surface area, whereas macrophytes had an inverse response to altitude, compared to ﬁsh. Both species turnover and nestedness signiﬁcantly contributed to the total β -diversity of macrophytes. In contrast, species turnover was the most signiﬁcant contributor to the total ﬁsh β -diversity. We found that the compositional variation of macrophytes is primarily limited by dispersal and ultimately shaped by environmental drivers, resulting in spatially structured assemblages. Fish communities were primarily shaped by altitude, highlighting the role of species sorting. We conclude that among-biota diversity patterns are shaped by di ﬀ erent / contrasting factors, and, thus, e ﬀ ective / sustainable conservation strategies should encompass multiple aquatic biota.

For αdiversity, various historical, environmental, and geographical factors [11] have long been recognized as major drivers, with their relative importance being the key approach to understanding the mechanisms and rules of community assembly [12,13]. For β-diversity, studies have shown that environmental filtering is the strongest mechanism, responsible for high compositional dissimilarities, as species sort themselves along environmental gradients [14][15][16][17][18]. Moreover, community dissimilarity increases as the geographical distance between two communities increases (distance decay) [19], since neighboring communities are more structurally similar than those located further apart [14]. This pattern also highlights the importance of dispersal limitations, further suggesting inherent neutrality and stochasticity [20]. As total β-diversity is partitioned into two additive components, accounting for pure spatial turnover and nestedness, many hypotheses have been proposed to describe the underlying ecological processes [21], since these two components imply different mechanisms structuring species assemblages. Spatial turnover refers to the replacement of some species by others as a consequence of environmental filtering and/or spatial and historical constraints [18,[22][23][24]. Conversely, nestedness reflects the spatial pattern, in which assemblages in depauperate ecosystems are subsets of communities of successively richer ecosystems [25] and stems from species loss caused by different processes, such as selective colonization and extinction [26]. Diversity-based research has long focused on terrestrial vertebrates and vascular plants, with studies on freshwater ecosystems being limited [27]. Few attempts have been made to investigate the effects of both space and environment on taxonomic richness and composition variation of diverse freshwater biota with varying dispersal traits [28][29][30][31][32][33], despite the fact that lacustrine habitats are analogous to terrestrial ones, often conceptualized as habitat islands [34].
The Balkan Peninsula has long been recognized as a European freshwater biodiversity hotspot [35] with extraordinary high compositional dissimilarity of its aquatic biota [36]. It is characterized by recent geological formations with distinctive geographical barriers, highly fragmented hydrographic networks [37], and the existence of many lakes, remnants of the Tethys and Paratethys seas [38].
In this study, we investigated the relative contribution/importance of environmental (habitat-, climate-and water quality-related) and geographical factors/gradients on shaping the αand β-diversity patterns of two highly diverse, differing in dispersal mechanisms freshwater biota, aquatic macrophytes and fish, across sixteen natural lakes of the Southern Balkan Peninsula. Our primary hypothesis was that macrophytes and fish show similar patterns of species richness along the same environmental/geographical gradients. We further assessed total β-diversity of each biota, partitioned into species turnover and nestedness, assuming that good dispersers, such as macrophytes, will have less biogeographical variation than actively dispersing organisms (bad dispersers), such as fish. Finally, we developed Generalized Dissimilarity Models (GDMs) to evaluate the relative contribution of environmental and geographical factors in explaining variation in the three β-diversity components (total, species turnover, and nestedness), hypothesizing that passively dispersing organisms with good dispersal ability would be subject to stronger environmental control/filtering than actively dispersing organisms (bad dispersers) across the Southern Balkan lakes.

Study Area and Data Collection
Sixteen natural lakes of the southern Balkans were studied, eleven exclusively located in Greece, one in Bulgaria, and four transboundary lakes, shared between two or more countries ( Figure 1). All lakes belong to five groups of drainage basins corresponding to five distinct, major biogeographical regions of the Balkan Peninsula, although with short geographical distance between them (see [30]). The eleven exclusively Greek lakes are included in the Natura 2000 network (except for the Ziros Lake) as sites of Community importance (SCIs), Special Conservation Areas (SCAs), and Special Protected Areas (SPAs). The transboundary Prespa lakes have been defined as a Greek National Park and a wetland of international importance under the Ramsar Convention. Lake Skadar is also a Ramsar site since 1996 and a national park in Montenegro since 1983. Lakes Ohrid between Albania and North Macedonia and Srebarna in Bulgaria are both listed as World Heritage Sites. All lakes were selected based on data availability and geographical representation (see Figure 1). Water 2020, 12, x FOR PEER REVIEW 3 of 16 Thirteen environmental variables were categorized into three groups in accordance with [39]: (1) climate (n = 6); (2) water quality (n = 2); and (3) habitat heterogeneity (n = 5) ( Table 1). All habitat heterogeneity values, except for maximum depth, were calculated using the lakemorpho package in R version 3.5.2 [40] using as input the lake polygon shapefiles obtained by the CORINE land cover 2012 inventory (https://land.copernicus.eu/pan-european/corine-land-cover/clc-2012/view). Monthly air temperatures and annual rainfall and precipitation for the year 2017 were acquired from the MERRA-2 database [41]. The elevation of each lake was obtained by accessing Google Earth imagery. The remaining variables, maximum depth, electrical conductivity, and pH, were obtained from previous literature (Table S1, Supplementary material). Thirteen environmental variables were categorized into three groups in accordance with [39]: (1) climate (n = 6); (2) water quality (n = 2); and (3) habitat heterogeneity (n = 5) ( Table 1). All habitat heterogeneity values, except for maximum depth, were calculated using the lakemorpho package in R version 3.5.2 [40] using as input the lake polygon shapefiles obtained by the CORINE land cover 2012 inventory (https://land.copernicus.eu/pan-european/corine-land-cover/clc-2012/view). Monthly air temperatures and annual rainfall and precipitation for the year 2017 were acquired from the MERRA-2 database [41]. The elevation of each lake was obtained by accessing Google Earth imagery. The remaining variables, maximum depth, electrical conductivity, and pH, were obtained from previous literature (Table S1, Supplementary material).

α-Diversity
Macrophyte and fish species richness and composition were assembled through an extensive literature review spanning form 1990 to the late 2010s, including research articles, dissertation theses, databases, and gray literature [36,[42][43][44][45][46][47][48][49]. Helophytes were excluded and only macrophytes of the isoetids, elodeids, nymphaeids, lemnids, and charids growth forms were analyzed. Aquatic ferns and bryophytes were also considered in further analyses. Fish species richness and composition data were also assembled from existing literature (research articles, databases, and gray literature) [36]. Only native species were examined; species with uncertain natural presence were omitted from analysis. Overall, our database encompassed 62 macrophyte and 95 fish species.

β-Diversity
A plethora of metrics has been proposed for quantifying beta diversity (see review [5]); yet, there is still a debate among researchers on the use of beta diversity indices [50]. Following Baselga [21,26], the relative contribution of species turnover and nestedness to the total beta diversity was determined by partitioning beta diversity into two separate components of species turnover and nestedness-resultant dissimilarities using the Sørensen dissimilarity index accounting for total beta diversity (β Sør ). β Sør was thus decomposed into β Sim (spatial turnover) and β Nes (nestedness); β Sør = β Sim + β Nes .
Based on the presence/absence matrix with rows representing the different lakes and columns representing the macrophyte or fish species composition in each lake, we calculated β Sør as follows: where α is the number of species common between two lakes, b is the number of species that occur in the first lake but not in the second, and c is the number of species that occur in the second lake but not in the first. Spatial turnover was determined by calculating the Simpson dissimilarity index (β Sim ) as follows: where α, b, and c denote the same variables as defined for β Sør . The nestedness-resultant component of Sørensen dissimilarity is simply the difference between β Sør and β Sim , and it denotes the fraction of dissimilarity caused by richness differences between nested subsets [26]. The three β-diversity indices (β Sør , β Sim , and β Nes ) were calculated using package betapart in R [51].

α-Diversity Patterns-Generalized Linear Models (GLMs)
To determine the environmental variables that best explained species richness patterns in the sixteen studied lakes, we implemented Generalized Linear Models (GLMs) with a Poisson error distribution and a logarithmic link function for macrophyte and fish species richness data. GLMs are a flexible extension of linear models, which allows for the incorporation of non-normal distributions of the response variable and nonlinear combinations of the explanatory variables. Prior to the analysis, to avoid multicollinearity and redundancy of variables, the number of variables included in the model was reduced by removing highly correlated variables (r > 0.7) through significance testing for each pair using the corrplot R package [52]. The variance inflation factor (VIF) was calculated using the car R package [53], and multicollinearity was assessed based on VIF; when VIF values >10 occurred, indicating high multicollinearity, only the variable with the lowest VIF was included in the analysis [54]. Based on the above, six variables, fetch, shoreline length, average monthly air temperature, maximum monthly air temperature, minimum monthly air temperature, and rainfall, were omitted from subsequent GLM analyses.
To determine the optimal model, we started with a global model containing all the explanatory variables included after the VIF assessment (n = 7). Then, we generated sub-model sets from the global model using the dredge function implemented in the MuMIn package [55]. We used the Akaike's information criterion corrected for small sample sizes (QAICc, [56]), to select explanatory variables driving total species richness, considering both goodness-of-fit and model complexity; the final model selected was the one with the lowest AICc value. GLMs were implemented with the MASS package [57]. AICc was calculated using the AICcmodavg v1.28 R package [58].

β-Diversity Patterns-Generalised Dissimilarity Modeling (GDM)
Generalized Dissimilarity Models (GDMs) [59,60] were developed for modeling spatial variation in macrophyte and fish assemblages. GDMs are a reformulation of the Mantel approach into a GLM, and allow for predictive modeling of community dissimilarity against a set of environmental predictors while overcoming two major limitations; nonlinearities typical in ecogeographic datasets, and uneven rates of species turnover along environmental gradients [60]. GDMs were developed based on the compositional dissimilarity of native fish and macrophyte assemblages, quantified with the three indices of beta diversity (β Sør , β Sim and β Nes ) between pairs of lakes in relation to lake-by-lake distance environmental matrices and geographical distance.

Predictors of α-Diversity
Macrophyte species richness ranged from 5 in the Zazari Lake to 31 in the Skadar Lake, the largest lake in the Balkans. The best GLM highlighted lake surface, pH, and altitude as significant descriptors of macrophyte richness, with lake surface being the most influential variable, followed by pH and altitude (Table 2 and Figure S1). For fish richness, shoreline development was the most significant predictor of α-diversity, showing a linear positive effect in fish species richness ( Table 2). Lake area also had a positive significant effect whereas altitude had negative correlation with fish species richness variation ( Figure S1). Regarding altitude, we found that fish species richness decreases substantially until the elevation of 400 m, above which it shows small changes, whereas aquatic macrophyte richness increases with altitude showing highest rates of change between sea level and altitudes of about 200 m. The response of both macrophyte and fish richness to lake area was unambiguous, with very similar patterns ( Figure S1).

Patterns of Dissimilarity/GDM Models of β-Diversity
For macrophytes, geographical distance was the unique important predictor for determining the Sørensen dissimilarity component (β Sør ); the model excluded all other environmental predictors (Table 4 and Figure 2). Geographical distance contributed to a lesser degree to the Simpson dissimilarity component (β Sim ) while conductivity was the most important predictor in the GDM model (0.311), both explaining 22.019% of the total deviance of the fitted GDM model ( Table 4). The rate and magnitude of species turnover along both gradients were exponential (Figure 2). The most important environmental gradient explaining the dissimilarity due to the nestedness component for macrophytes was altitude, which increased exponentially (coefficient: 0.170), with a small contribution of geographical distance (nonlinear increase, coefficient: 0.095), jointly explaining 8.032% of the total deviance (Table 4 and Figure 2). Table 4. Summary of the Generalized Dissimilarity Model (GDM), in which significant variables and their relative importance in the GDMs are presented (relative importance is determined by summing the coefficients of the I-splines from GDM). The percentage of the total deviance explained by the fitted GDM and the null deviance are also presented.

Macrophytes
Null  Figure 2). For freshwater fishes, βSør and βSim increased nonlinearly with geographical distance and exponentially with altitude; both predictors explained a small portion of the total variance of components (βSør: 10.910% and βSim: 9.082%) ( Table 4 and Figure 3). However, altitude was the most important predictor of βSør and βSim (2.245 and 1.931, respectively). Geographical distance was identified as significant also for βSør and βSim of fishes, although to a lesser degree (0.796 and 0.158). The contribution of additional variables was negligible (Table 4).

Figure 2.
Relationships between the observed macrophyte Sorensen dissimilarity (β Sør ), the observed macrophyte Simpson dissimilarity (β Sim ), and the observed macrophyte nestedness component (β Nes ) with the regression equation from the Generalized Dissimilarity Model (GDM) (termed ecological distance) and the observed values versus the predicted dissimilarity. The fitted I-splines derived from the GDMs are shown for each predictor retained in the selected models, depicting the nonlinear rate of turnover across each environmental gradient (geographical distance, conductivity, and altitude). The maximum height reached by each curve indicates the total amount of compositional dissimilarity associated with that variable (it also indicates the relative importance of that variable in explaining beta diversity), holding all other variables constant. The shape of each function provides an indication of how the rate of compositional turnover varies along the gradient.
For freshwater fishes, β Sør and β Sim increased nonlinearly with geographical distance and exponentially with altitude; both predictors explained a small portion of the total variance of components (β Sør: 10.910% and β Sim : 9.082%) ( Table 4 and Figure 3). However, altitude was the most important predictor of β Sør and β Sim (2.245 and 1.931, respectively). Geographical distance was identified as significant also for β Sør and β Sim of fishes, although to a lesser degree (0.796 and 0.158). The contribution of additional variables was negligible (Table 4).

Figure 3.
Relationships between the observed fish Sorensen dissimilarity (βSør) and the observed fish Simpson dissimilarity (βSim) with the regression equation from the Generalized Dissimilarity Model GDM (termed ecological distance) and the observed values versus the predicted dissimilarity. The fitted I-splines derived from the GDMs (partial regression fits) are shown for each predictor retained in the selected models, depicting the nonlinear rate across each environmental gradient (geographical distance and altitude). The maximum height reached by each curve indicates the total amount of compositional dissimilarity associated with that variable (it also indicates the relative importance of that variable in explaining beta diversity), holding all other variables constant. The shape of each function provides an indication of how the rate of compositional turnover varies along the gradient. Table 4. Summary of the Generalized Dissimilarity Model (GDM), in which significant variables and their relative importance in the GDMs are presented (relative importance is determined by summing the coefficients of the I-splines from GDM). The percentage of the total deviance explained by the fitted GDM and the null deviance are also presented.

Discussion
Ecologists have long shown increased interest in disentangling αand β-diversity patterns of biota, since understanding the mechanisms responsible for the complexity and maintenance of biodiversity remains at the heart of biogeography, ecology, and conservation biology. Moreover, investigations of spatial β-diversity components are considered pivotal for protecting biodiversity and can directly assist conservation planning [61]. In particular, species turnover indicates a relatively stronger species-sorting process, whereas nestedness could enhance the spatial signal [62]. If nestedness contributes more than turnover into total βdiversity, this might suggest prioritization of a small number of ecosystems with high αdiversity. In contrast, if β-diversity is attributed to high turnover, conservation must target multiple sites though not necessarily the richest [26,63,64].
In this study, we performed GLMs and GDMs to identify if there is congruence between αand β-diversity (and its components) for two biotic groups (macrophytes and fish) in sixteen natural lakes of the unexplored biodiversity hotspot of the Southern Balkan Peninsula. Previous studies in the freshwater realm have reported contrasting results on the relationships between diversity and environment for different biota [32,33]. However, there is no previous study in the area, comparing the beta diversity patterns of multiple biota based on a common spatial scale and identical statistical analysis.

α-Diversity
In agreement with our primary hypothesis, species richness was strongly associated with lake surface area in the Balkan lentic systems. As lakes are conceived as habitat islands, a positive relationship of lake area with macrophyte richness has been well established previously [65][66][67][68] and is well founded in island biogeography reflecting habitat heterogeneity. As larger lakes will likely have more microhabitats and therefore more available niches, more species are able to find a suitable habitat [65,69]. Fish richness per lake was also significantly correlated with lake area, a result corroborating previous findings from North America and Africa [70,71], Asia, and Europe [30,36,72,73]. The hypothesis that increased habitat heterogeneity positively affects species richness is further supported by the fact that shoreline development, which indicates a complex and expanded littoral zone that hosts diverse habitats, was identified as an important predictor of fish richness [74]. However, most of the Skadar Lake fish species, which has the largest shoreline development, move into littoral and flooded zones to spawn (nursery habitats). To sum up, lakes with more habitats support more macrophyte and fish species. Our results also highlighted pH as a key environmental variable for macrophyte richness. The effects of pH and alkalinity in shaping macrophyte richness have long been pinpointed by several studies in Europe [75,76]. In this regard, Alahuhta et al. [68] also highlighted a positive response of macrophyte richness to alkalinity in lakes of the USA, while, more recently, Elo et al. [77] demonstrated a strong positive relationship between richness and a combined factor of conductivity, pH, and alkalinity. The pH level in most freshwater systems is regulated by the concentration of bicarbonate, which serves as a source of carbon for photosynthesis of submerged and floating-leaved plants that in turn influences the growth and survival of aquatic plants. Particularly, in Greece, previous works have outlined a strong association between freshwater biota diversity patterns and physical, chemical and geomorphological variables including surface area, lake depth, Secchi depth, conductivity, and temperature [30,72,78].
Macrophyte richness was positively influenced by altitude, a contrasting outcome to several previous studies [39,67,79]. For instance, Lacoul and Feedman [79] showed a clear declining trend of species richness along a wide range of altitudes (from 77 to 4750 m), being attributed to adverse climate conditions in higher altitudes that restrict the distribution of certain aquatic plants. Similarly, Alahuhta [39] underpinned the importance of higher temperatures during the growing season and winter in the lowlands compared to higher altitudes, which may favor many aquatic plants, thus explaining the negative relationship between altitude and macrophyte richness. Fewer studies have shown an opposite pattern [80] suggesting that the response of macrophyte richness to altitude varies and it may depend on the altitudinal range [68,81]. Our results are in agreement with those of Chappuis et al. [80], providing evidence of a lacking effect of physical barriers (colonization restrictions) and climatic factors, such as low temperature. Moreover, it has been reported that the predominant anthropogenic pressures in lowland lakes of the Mediterranean region (hypereutrophication, weed cutting, and water level fluctuations) could be connected to species loss and lower richness [80,81], offering a possible explanation for our results. Contrary to macrophytes, fish species richness decreased with increasing altitude where low connectivity and low temperatures set an obstacle to fish colonization, suggesting a strong climatic control on Balkan freshwater fish biodiversity, as it has been extensively reported elsewhere [71,72,82]. As altitude is negatively correlated to temperature, it was not possible to clearly disentangle the relative effect of altitude and temperature in our study. It should be noted that the altitudinal gradient also reflects the biogeographical history of fish with rich ecosystems acting as refugia.

β-Diversity
The study of β-diversity is crucial to pinpoint the factors driving biodiversity patterns, especially in highly heterogeneous environments, such as the Balkan Peninsula [36]. We found that different factors shaped the high compositional dissimilarity of communities in Balkan lentic habitats, although total β-diversity was lower for the macrophyte communities (β Sør = 0.601) than for fish (β Sør = 0.881). It was further shown that the compositional dissimilarity of the two biota is related to both geographical distance and different environmental gradients, demonstrating that dispersal history, ecological and biological characteristics and the different response of each biotic group to environmental gradients are responsible for the observed beta diversity patterns, which has been previously evidenced for gastropods [83].
For macrophytes, total beta diversity (β Sør ) was mostly influenced by the geographical distance among lakes. Although aquatic macrophytes are known for their broad distributional ranges due to their long-distance seed dispersal strategies such as anemochory, hydrochory, and zoochory [84,85], there are studies that report that some aquatic plants may disperse more or less uniformly at distances up to~200 km, beyond which habitat isolation usually becomes the primary limiting factor [86].
Other studies have shown that macrophyte communities are likely more influenced by spatial dispersal processes than environmental filtering in regions with a highly variable topography and dispersal barriers, such as complex mountainous environments [86]. These results are in agreement with our findings, which showed that total dissimilarity increased with geographical distance, even for distances smaller than 200 km, reflecting the significance of the exploration of biogeographical patterns in the Balkan Peninsula. Contrary to macrophytes, we found that total fish β-diversity in Balkan lakes is mainly constrained by the altitudinal gradient, suggesting that species with narrow ranges are mainly restricted to systems located in higher altitudes. Altitude, therefore, acts as an environmental filter, selecting fish species within the corresponding thermal niche. Thus, niche processes are dominant in shaping fish communities, whereas the weak impact of geographical distance in fish compositional dissimilarity may be due to their dispersal capabilities, which are hindered by habitat fragmentation and poor habitat connectivity [87]. Nevertheless, there is evidence that dispersal may play a significant role in determining community structure of Greek freshwater assemblages, as has been shown by previous studies (e.g., [30,72]). These discrepancies in findings may be related to the differences in data time span and spatial extent of ecosystems actually used in the analyses.
Applying an additive partitioning framework, we found that species replacement (turnover) was the major component of macrophyte and fish total β-diversity in Balkan lakes, although the effect was stronger for fish (0.371 and 0.824 for macrophytes and fish, respectively). The complementary nestedness component accounted for 23% of total β diversity for macrophytes whereas it was found extremely low for fish (β Nes < 0.1). The preponderance of the turnover component over the nestedness component has already been highlighted in previous studies [33,88,89].
Macrophyte and fish beta diversity components (β Sim and β Nes ) had differing responses to environmental gradients and geographical distance mainly due to differences in the mode of dispersal, e.g., good or bad dispersal ability. The compositional change (turnover component) in macrophyte assemblages is limited primarily by dispersal. Our results suggest that aquatic macrophyte turnover and nestedness are also governed by conductivity and altitude, respectively. Regarding conductivity, the results imply that water quality plays a role in defining the environmental conditions in which certain species are favored over others. Conductivity in particular has been highlighted as a key driver of macrophyte composition in the Mediterranean [80,90,91]. For instance, Fernández-Aláez et al. [81] showed that conductivity was the main predictor of macrophyte assemblages in a large number of Spanish ponds along an altitudinal gradient. Thus, high conductivity (>800 µS/cm) favored macrophyte assemblages dominated by few tolerant species such as Ceratophyllum demersum and Stuckenia pectinata, whereas waters with low ionic content could support communities with aquatic bryophytes and isoetids. Moreover, random species extinctions or colonization events across the altitudinal gradient may have led to medium levels of nestedness in macrophyte assemblages. Contrary, environmental filtering (the altitudinal gradient) was the main force shaping fish turnover adding evidence to the role of species sorting processes in controlling the variation in species composition [92]. In particular, lakes in high altitudes could have acted as refugia during the glacial periods or as speciation centers (Prespa Lakes' ecoregion [36]). In conclusion, endemism may have a direct effect on fish assemblages by promoting speciation (see also [93] for consistency). These results are, therefore, in line with a growing body of evidence, suggesting a positive relationship between species composition and environmental heterogeneity in fish communities [32], although other studies [30,72] have shown that dispersal limitation played a more significant role in shaping Greek freshwater fish communities than environmental heterogeneity. Last, the nestedness component of beta diversity was negligible for fish (β Nes < 0.1), as has been reported elsewhere [36].

Conservation Implications
The current study showed that α-diversity patterns across the two biological groups of macrophytes and fish are not uniform and were shaped by different drivers under the same spatial scale, suggesting differences in their evolutionary histories. We also found weak congruence in the main factors driving βdiversity of the two biological groups, probably due to the different dispersal modes. Our findings support previous studies indicating that species richness and compositional dissimilarity among freshwater biota are shaped by different/contrasting factors [94] and, thus, conservation plans and bioassessments of human disturbance need to include multiple biota. Our results further suggest that the partitioning of β-diversity components should be considered in the future for optimizing conservation strategies. We showed that species turnover was the major driver of β-diversity and, thus, conservation efforts must aim to increase the network of protected freshwater areas across the Balkan Peninsula to maximize the protection of species diversity [95]. Ultimately, this study may serve as a benchmark for future studies towards ensuring the sustainable management of the Balkan freshwater biodiversity. In this context, our findings have potentially important implications for the European research agenda regarding freshwater ecosystems.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/7/1984/s1, Table S1: Environmental data of lakes extracted from the literature and considered in the analyses. Figure S1: Partial dependence plots showing the responses of fish (left) and macrophyte richness (right) to each predictor.