Species Richness and Taxonomic Distinctness of Zooplankton in Ponds and Small Lakes from Albania and North Macedonia: The Role of Bioclimatic Factors

Resolving the contribution to biodiversity patterns of regional-scale environmental drivers is, to date, essential in the implementation of effective conservation strategies. Here, we assessed the species richness S and taxonomic distinctness Δ+ (used a proxy of phylogenetic diversity) of crustacean zooplankton assemblages from 40 ponds and small lakes located in Albania and North Macedonia and tested whether they could be predicted by waterbodies’ landscape characteristics (area, perimeter, and altitude), together with local bioclimatic conditions that were derived from Wordclim and MODIS databases. The results showed that a minimum adequate model, including the positive effects of non-arboreal vegetation cover and temperature seasonality, together with the negative influence of the mean temperature of the wettest quarter, effectively predicted assemblages’ variation in species richness. In contrast, taxonomic distinctness did not predictably respond to landscape or bioclimatic factors. Noticeably, waterbodies’ area showed a generally low prediction power for both S and Δ+. Additionally, an in-depth analysis of assemblages’ species composition indicated the occurrence of two distinct groups of waterbodies characterized by different species and different precipitation and temperature regimes. Our findings indicated that the classical species-area relationship hypothesis is inadequate in explaining the diversity of crustacean zooplankton assemblages characterizing the waterbodies under analysis. In contrast, local bioclimatic factors might affect the species richness and composition, but not their phylogenetic diversity, the latter likely to be influenced by long-term adaptation mechanisms.


Introduction
Lentic freshwaters are acknowledged to play a crucial role in regulating the global ecosystem functions e.g., carbon cycle [1] and they are among the Earth's most threatened habitats in terms of intensity of anthropogenic pressures, biodiversity loss, and non-indigenous species introduction [2][3][4]. They include an extreme variety of habitats differing in ecological characteristics and fragility [5,6]. Surface area represents one of the most apparent differentiating properties: indeed, lentic environments (304 million water bodies; 4.2 million km 2 in total area [7]) include the lake Superior (82,000 km 2 ), together with small ponds, i.e., waterbodies less than 0.05 km 2 in area [8].
Ponds and small lakes (hereafter PSL) have significant ecological functions [9,10]: among others, they provide a considerable contribution to inland water CO 2 and CH 4 emissions [11]. In addition, they are, to date, recognized as important biodiversity hotspots, especially in mountainous regions, supporting a high species richness and contributing a high degree of rare species to regional pools ( [12][13][14][15]; see also [16] for an example on planktonic Calanoida). Noticeably, PSL are threatened by a number of anthropogenic pressures, including nutrient loading, contamination, acid rain, and invasion of exotic species [17]. In addition, infilling (both natural and caused by direct habitat destruction), land drainage, decline in many of their traditional uses, and changes of function determine at a regional scale the drastic reduction in PSL number and connectivity [12].
In the last decade, several investigations have focused on the diversity of benthic invertebrates, as they are excellent bio-indicators of PSL ecological integrity [18,19]. Local factors that are related with e.g., hydroperiod, environmental harshness, water chemistry, spatial connectivity, habitat heterogeneity, and presence of predators, have been recognized to influence the biodiversity of macroinvertebrate assemblages ( [20] and literature cited). At a regional scale, attention has been primarily given to the influence of waterbodies area [21][22][23], while assuming, within the general theoretical background provided by the species-area relationship (SAR) hypothesis [24], that basins' size correlates with the number of microhabitats within the basin itself and with populations' abundance, and thence inversely correlated with the likelihood of random extinctions. However, resolving the contribution to biodiversity patterns of environmental factors acting at a regional scale is, to date, essential to the implementation of effective conservation strategies in the face of e.g. deforestation and climate change ( [25,26] and literature cited). Accordingly, several attempts have been made to model biodiversity of freshwater environments by means of regional bioclimatic factors [27,28].
In the present study, we focused on the diversity of crustacean zooplankton assemblages in 40 ponds and small lakes differing remarkably in terms of origin, extension, and altitude from a relatively wide region comprising part of Albania and North Macedonia. A recent faunal inventory focusing on ponds and lakes in the area [29] provided the starting reference information on the taxonomic characteristics of the assemblages.
A number of studies have generally indicated a positive relationship between the surface area of lacustrine environments and zooplankton diversity (e.g., [30][31][32]; but see [13]); accordingly, crustacean zooplankton has been shown to have higher species richness in small ponds as compared with lakes [16,33,34]. This notwithstanding, we hypothesized that area alone may not be an adequate predictor, and that local bioclimatic conditions may ultimately contribute in explaining diversity variations across waterbodies by influencing their physical-chemical characteristics, as observed in recent investigations on freshwater macroinvertebrates and macrophytes [27,28,35]. This could be particularly true for waterbodies in mountainous habitats, where temperature and precipitation regimes intensely reflect the chemical-physical characteristics and hydroperiod of the waterbodies themselves [36], regulating the harshness and stability of the aquatic environments and, in turn, the diversity of the biota living in them ( [22] and literature cited).
To verify the hypothesis and test whether bioclimatic factors can predict assemblages' diversity, we identified a minimum adequate model (MAM) predicting assemblages' diversity across the different waterbodies by means of a heuristic multiple regression approach and Bayesian Information Criterion model selection method while using satellite-derived bioclimatic variables as predictors. Multiple indices are, to date, available to quantify different aspects of biodiversity [37]. Here, we identified predictive MAMs estimating the diversity of crustacean zooplankton assemblages in terms of species richness and average taxonomic distinctness. Species richness is the most classical measure of biodiversity across ecosystems that has been extensively used in studies on lentic habitats (see references cited above). This index provides an incomplete understanding of biological variability, because it neglects information on the identity and taxonomic relationship among species, and it is hampered by a number of critical limitations [38,39]. Accordingly, we used the average taxonomic distinctness ∆+ [40] to compare the taxonomic relatedness of species in the crustacean assemblages of every water body. In addition, we tested the influence of bioclimatic factors on crustacean assemblages in terms of species composition. To this end, multivariate approaches that are based on a canonical analysis of principal coordinates were used to model the changes in the structure of the assemblages as affected by bioclimatic variables, and identify relationships between the latter and specific groups of zooplankton taxa.

Sampling Sites and Collection Procedures
A total of 40 sites were selected among those (53) that were surveyed between 2005 and 2017 by Belmonte and colleagues [29] in an area comprised between 39 • 55 22 -42 • Figure 1 is included, together with information of the basins origin (N = natural, A = artificial), typology (L= lake, P = pond (area < 0.05 km 2 )), elevation (in m), coordinates, number of samples collected, year of collection, species richness S, and taxonomic distinctness ∆+. Sampling procedures are described in detail elsewhere [29]. As the water bodies varied in terms of area, altitude, origin (natural, artificial), as well as in hydrology (permanent, temporary), banks morphology, and degree of aquatic vegetation cover, it was not possible to apply a standardized sampling protocol across all of the sites.

Waterbody
In the selected sites, sampling operations were carried out in spring-early summer (i.e., between April and July) always by the same operators, while using a hand-held plankton net (200 µm mesh-size, mouth diameter, 30 cm). The collection (by plankton net towing from two opposite edges of the pond) covered the whole water body when its diameter was smaller than 100 m. For larger water bodies, sample collection was carried while using a canoe. In the case of small ponds, the collection procedure was repeated three times (each sample derived from the execution of three collections). In the case of larger water bodies, a sample collection was carried out in three different stations; the three different samples were ultimately cumulated.
After collection, the samples were fixed in situ in 90-96% ethanol. In the laboratory, taxa were identified to the species level while using a compound microscope (×30-×300 magnifications) that was equipped with a camera lucida.
The quantification of the abundance of each taxon was not performed, and only presence/absence data were considered due to the huge variability of water volumes in each pond, which made impossible the comparison among the concentrations of plankton of different sites.

Landscape-Climate Variables
Together with altitude, we used area and perimeter of the water bodies together with bioclimatic factors (i.e., temperature and precipitation) to represent the landscape-climate variables. Water bodies were geo-referenced in Google Earth Pro version 7.3.2., where their surface (in km 2 ) and perimeter (in km) were measured while using the software tools. Measurements were performed by preferentially choosing images that were taken in spring or summer between 2008 and 2017, assuming that negligible variations in the water bodies size and morphology occurred during this period.
Nineteen climatic layers with a 30-second spatial resolution (0.93 × 0.93 = 0.86 km 2 at the equator; approximately 0.92 × 0.70 = 0.64 km 2 within the study area), including temperature and precipitation variables, were extracted from the WorldClim v2 data set [41] (Table A1 in online material). Besides climate, there is a growing recognition of the importance of vegetation cover in characterizing the spatial environmental heterogeneity in a given area at the meso-and topo-scales, in turn affecting climate, soil composition, hydrology and geomorphology, and, ultimately, biological processes that were related with species richness and community complexity [42,43]. Accordingly, two vegetation variables (i.e., percent tree cover and percent non-tree cover) were extracted from the Terra MODIS Vegetation Continuous Field (VCF) product (available as MOD44B v006 https://lpdaac.usgs.gov/products/mod44bv006/). Percent tree cover included all forest types and age classes, while percent non-tree cover included meadows, regeneration areas, and clear-cut areas. MODIS tiles of the study area were re-projected and re-sampled to meet the coordinate system and resolution of WorldClim layers; percent cover data were subsequently obtained for the years from 2006 to 2017, and averaged. In addition, the third VCF component of ground cover, i.e., percent bare soil (including bare soils and rocks) was extracted according to the aforementioned procedures, and was used together with tree and non-tree vegetation percent cover data to estimate the Shannon's diversity index (H) as a proxy of habitat heterogeneity.

Data Analysis
The values in the text are expressed as averages ± 1SE; for parametric statistical analysis, data were tested for conformity to assumptions of variance homogeneity (Cochran's C test) and normality (Shapiro-Wilks test) and transformed when required.
The taxonomic diversity of crustacean assemblages in each water body was estimated in terms of species richness S and average taxonomic distinctness ∆+. The index can be used as a proxy for phylogenetic diversity and it measures the mean path length through a taxonomic tree connecting every species [40]. Here, mean taxonomic distinctness values were calculated assigning equal weighting to branch lengths from a linear Linnaean classification while using eight taxonomic levels (i.e., class, subclass, order, suborder, infraorder, family, genus, and species). The taxonomic classification tree was built according with the World Register of Marine Species (WoRMS, available at https://www.marinespecies.org) and the Integrated Taxonomic Information System (ITIS, available at https://www.itis.gov).
We verified the influence on both indices of potential artefacts, due to (i) possible differences in sampling procedures and (ii) the different number of total collected samples per water body. To this end, we performed a one way permutational analysis of variance (PERMANOVA; [47]) based on Euclidean distances and 999 permutations with "sampling procedure" as a fixed factor (two levels, "hand", or "canoe") and the number of collected samples as the covariate (P and N hereafter). Both factors exerted negligible influences on the S and ∆+ estimations (S: Pseudo-F P = 4.42, P(perm) P = 0.08, Consequently, the S and ∆+ values were assumed to provide robust estimation of planktonic Crustacea diversity across the studied water bodies. PERMANOVA was further used to test the effects of the factors "origin" (two levels, "natural" and "artificial") and "typology" (two levels, "pond" and "lake") on the diversity indices. As water bodies varied greatly in elevation (Table 1), the latter was included in the analyses after log-transformation as a continuous covariate.
The georeferenced locations of the sampling sites were used to extract climatic and vegetation data from environmental layers. The final data set included 19 climatic, two vegetation (% tree cover, % non-tree cover), four geomorphological (elevation, perimeter, surface, and perimeter/surface ratio), and one habitat heterogeneity variable (Table A1). They were log-transformed and z-scaled; subsequently, their original number (25) was reduced while using an iterative variance inflation factor (VIF) analysis [48,49]. In brief, if a strong linear relationship links a variable x with at least another variable y, the correlation coefficient would be close to 1, and the VIF for x would be large. Here, diversity measures with VIF factors that were larger than 10 were excluded. Variables with VIF factors larger than 10 were discarded. The identification of a minimum adequate model (MAM hereafter; [50]) linking diversity measures with environmental variables was based on the heuristic generation of alternative regression models (see [51] for complete details on the procedure). Model selection was performed while adopting an Information Theoretic criterion [52]; the second-order Akaike Information Criterion AICc [53,54] was calculated for each combination of n explanatory variables and used to identify the best MAM among the alternative regression models that were generated by the procedure. For model comparison, AICc values were used to estimate a set of positive Akaike weights w i summing 1: With K = number of models. The model showing the highest w i was accepted as the best candidate; other candidate models were accepted if characterized by w i values within 12.5% of the highest [55][56][57]. The model building and MAM identification procedures were performed while adopting Fox and Weisberg [58] as a general reference.
Given the non-conclusive outcomes of the analyses that were performed on taxonomic distinctness (see Results section), we verified whether bioclimatic factors influenced planktonic Crustacea assemblages in terms of species composition. Species incidence data were used to calculate a Jaccard similarity matrix across the different waterbodies. In addition, a similarity matrix that was based on Euclidean distances was constructed for bioclimatic variables, and the consistency of the two matrices was tested while using the Kendall coefficient of concordance (W). Subsequently, a canonical analysis of principal coordinates (CAP) [59] was performed to model changes in assemblages composition, as affected by bioclimatic factors. The appropriate number of principal coordinates m was chosen as to minimise the P value from the permutation test based upon the trace statistic and maximizing the leave-one-out allocation success [60]. Post-hoc PERMANOVA tests were performed to confirm the results of the ordination for both bioclimatic factors and species; SIMPER analyses were performed on the Euclidean distance matrix of the former to assess the percentage contribution of each factor to the dissimilarity between the groups of waterbodies that were identified by the CAP procedure. Furthermore, the Spearman's rank correlations were estimated to identify the bioclimatic variables and the crustacean species that most effectively described the groups of waterbodies that were identified by the CAP procedure. Only the variables with a Spearman rank correlation coefficient r > 0.55 were considered.
All of the analyses were implemented in the R statistical environment v3.6.1 [61] while using a suite of packages including taxize (for taxonomic information retrieval from online databases) vegan (for diversity measures and multivariate analyses), raster, rgdal, and maptools (for environmental layers manipulation), usdm (for VIF analysis), car, leaps, and HH (for MAM identification).

General Features
The 40 water bodies analysed, varied remarkably in terms of altitude, area, and perimeter ( In the 40 water bodies, 79 Crustacea species were identified in total, being almost equally distributed between the classes Branchiopoda and Hexanauplia (41 and 38 species respectively).
Among Branchiopoda, Anomopoda outnumbered the other two orders Ctenopoda and Haplopoda (38 vs. 2 and 1 species). Daphnia, Ceriodaphnia, and Moina were the genera that were characterized by the highest number of species (nine Daphnia species, four Moina species, and three Ceriodaphnia species) together representing the majority of all the sampled Anomopoda species. Ctenopoda were represented by the congeneric Diaphanosoma brachyurum and D. lacustris while Haplopoda by the single species Leptodora kindtii. The class Hexanauplia (alias Copepoda) was dominated by species belonging to the order Cyclopoida (31) and to a minor extent Calanoida (7). The genus Cyclops (seven species) together with Acanthocyclops, Paracyclops (four species each), and Mesocyclops (three species) constituted to the majority of the species in the order; Calanoida were represented by the genera Eudiaptomus (three species), Mixodiaptomus (two species each), and by Arctodiaptomus salinus and Neodiaptomus schmackeri.   Table 2); in addition, the Cyclopoida Mesocyclops leuckarti, the Calanoida Mixodiaptomus tatricus, and the Ctenopoda Diaphanosoma brachyurum occurred in 12 sampled sites. Table 2. Summary of PERMANOVAs on species richness S and taxonomic distinctness + of crustacean zooplankton assemblages testing for the effects of waterbodies' origin and typology including elevation as a continuous covariate **: p < 0.01. As regarding high level taxa, only Anomopoda (Cladocera) were present in every site. Cyclopoida (Hexanauplia) were present in 38 sites (95% of the total), Calanoida (Hexanauplia) in 30 sites (75%), Ctenopoda (Cladocera) in four sites (10%), and Haplopoda (Cladocera) in three sites (7.5%).
The factor "origin" was the only exerting significant effects of the species richness of waterbodies (Table 2), with the six artificial water bodies being included in the study characterized by lower S values as compared with natural basins (3.5 ± 0.22 vs. 7.29 ± 0.38, respectively). Conversely, negligible effects were generally observed for the taxonomic distinctness ∆+ ( Table 2).
The 25 predictor variables (Table A1) were reduced to a set of nine characterized by negligible collinearity (Table A2). They included five climatic variables (i.e., Isothermality, Temperature Seasonality, Mean Temperature of Wettest Quarter, Annual Precipitation, and Precipitation of Coldest Quarter), % tree and % non-tree vegetation cover, habitat heterogeneity, and water body surface. Besides water body surface, the variable Mean Temperature of Wettest Quarter showed the greatest among-water bodies variability (Table A2), ranging from a minimum of −5.32 • C (SHB 2) to a maximum of 11 • C (ZVE). It was followed by % tree cover (varying between 1 and 70%, BEL and DEG, respectively) and % non-tree cover (ranging between 20 and approx. 81%, BEL-DEG and SHR 2, respectively).
The heuristic search procedure identified a Minimum Adequate Model (MAM) predicting the variation of species richness S across the different water bodies relying on the three explanatory variables % non tree cover, Temperature Seasonality, and Mean Temperature of Wettest Quarter (Figure 3; multiple r = 0.58, P = 0.002, d.f. = 3, 36). The MAM was characterized by the lowest AICc value, and by an Akaike weight w i approximately eight times larger than the second-best candidate, based on the variables Temperature Seasonality, Mean Temperature of Wettest Quarter, and Habitat heterogeneity ( Table 3). All of the predictors provided significant contributions to S variation across water bodies (minimum absolute t value = 2.34, P = 0.02 for the variable Mean Temperature of Wettest Quarter); the contributions of both % non-tree cover and Temperature Seasonality were positive (b = 0.27 ± 0.12 and 0.64 ± 0.15, respectively), while the Mean Temperature of the Wettest Quarter provided a negative contribution (b = −0.26 ± 0.14). Noticeably, none of the first ten best-performing models included water body area (Table 3).  Table  3 for predictor variables) identified adopting a multiple regression procedure (top). Dashed curves  Table 3 for predictor variables) identified adopting a multiple regression procedure (top). Dashed curves are 95% CI. For the sake of completeness the predicted taxonomic distinctness values ∆+ by the best MAM are also reported (bottom). Table 3. Summary of the heuristic multiple regression analysis followed by a parsimonious selection procedure of the Minimum Adequate Model (MAM) predicting species richness (S) and of crustacean zooplankton assemblages by means of bioclimatic variables; only the first 10 best MAMs are reported. For the sake of completeness, results of MAM analysis are reported also for taxonomic distinctness (∆+), even though the statistical power of the models was negligible (see Results). K: number of predictors included in the model; AICc: second-order Akaike Information Criterion; w i : Akaike weight. For predictor abbreviations see Table A1. In contrast with species richness, the heuristic search procedure was unable to identify a MAM with a significant predictive power for taxonomic distinctness. The single variable % tree cover, resulted the best predictor of ∆+ (Table 3); however, the correlation resulted in being non-significant (r = 0.25, P = 0.11, d.f. 1,38; see also Figure 3). Other models showed an even worst performance, independently from the number of variables involved (Table 3), indicating, in turn, that the taxonomic distinctness of the crustacean assemblages cannot be predicted by bioclimatic, landscape-scale factors.

Species
Canonical analysis of principal coordinates (CAP), followed by a confirmatory PERMANOVA test identified two main groups of waterbodies significantly different in terms of species composition (Figure 4; Pseudo-F = 7.2, P(perm) = 0.001). The variables Isothermality, Mean Temperature of Wettest Quarter, Annual Precipitation, and Precipitation of Coldest Quarter showed a correlation (r > 0.65) with the canonical axis 1 (Figure 4) and significantly differed between the two groups of waterbodies (PERMANOVA, Pseudo-F = 27.4, P(perm) = 0.001). A Simper procedure indicated that the variable Mean Temperature of Wettest Quarter contributed by 32.2% to inter-group differences, followed by Isothermality and Annual Precipitation (28.6 and 24.2%, respectively). test identified two main groups of waterbodies significantly different in terms of species composition (Figure 4; Pseudo-F = 7.2, P(perm) = 0.001). The variables Isothermality, Mean Temperature of Wettest Quarter, Annual Precipitation, and Precipitation of Coldest Quarter showed a correlation (r > 0.65) with the canonical axis 1 (Figure 4) and significantly differed between the two groups of waterbodies (PERMANOVA, Pseudo-F = 27.4, P(perm) = 0.001). A Simper procedure indicated that the variable Mean Temperature of Wettest Quarter contributed by 32.2% to inter-group differences, followed by Isothermality and Annual Precipitation (28.6 and 24.2%, respectively).  (Table 1); their color categorizes the two groups (i.e., red for group1 and blue for group2) showing significant differences in species composition (post hoc PERMANOVA, P(perm) < 0.001). Vector overlay are Spearman correlations of bioclimatic factors and species with canonical axes with r > 0. 55. In group1, the Mean Temperature of Wettest Quarter was remarkably lower than that determined for group2 (0.06 ± 0.69 vs. 8 The analysis of the relationships between species occurrences and the canonical axis 1 (see Table A2 for Spearman correlations for the complete list of species) indicated that Eucyclops serrulatus, Chydorus sphaericus, Mixodiaptomus tatricus, and Daphnia rosea in the first group were correlated mainly with Precipitation. The occurrences of Neodiaptomus schmackeri, Diapahnosoma brachyurum, Cyclops scutifer, Ergasilus sp, Bosmina longirostris, and, to a lesser extent, Mesocyclops leuckarti in the group2 were related with the variables Isothermality, Mean Temperature of Wettest Quarter, and Precipitation of Coldest Quarter.

Discussion
The earth is undergoing an accelerated rate of native ecosystem conversion and degradation and there is increased interest in measuring and modelling biodiversity while using landscape-scale, remotely-sensed predictors. Here, we made an attempt towards this direction while using crustacean zooplankton. This group of organisms, ubiquitous in lentic habitats, has been recently subjected to renewed interest as an effective bio-indicator of the environmental status of ponds and small lakes [62][63][64][65]. The heuristic procedure was used here to identify minimum adequate models predicting the diversity the crustacean zooplankton assemblages across the 40 waterbodies under analysis provided non-univocal results. On one hand, they showed that a subset of bioclimatic variables could effectively predict the variation in species richness and composition across the different waterbodies. On the other hand, they also indicated that the assemblages' taxonomic distinctness ∆+ is unrelated with landscape-scale environmental drivers.
Before discussing these results, it must be considered that the emphasis we put on landscape and bioclimatic drivers of zooplankton diversity by no means imply that other chemical, physical, and biotic characteristics of the waterbodies, such as nutrient concentration, pH, depth, predators, aquatic vegetation, etc. are to be considered of secondary importance. A number of studies have unequivocally indicated that they can directly affect zooplankton species richness ( [66,67] and references cited in the introduction; but see also further in this section). In addition, lake age [68,69], connectivity, and, in general, the spatial arrangement of the habitat have been acknowledged to play an important structuring role (e.g., [70]; see also [71] for a marine example). However, in the present study, a hierarchy of effects at different spatial and environmental scales is implicitly assumed, with landscape and bioclimatic drivers indirectly affecting the characteristics of the biota (including crustacean zooplankton) by affecting the chemical and physical conditions of the waterbodies. Indeed, the limited extension of the basins that were included in our study (Figure 1) actually implies for them a low thermal inertia, and thus the ability to rapidly respond to the external climatic conditions [72]. An indirect support to this view is also provided by the negligible predictive power of waterbodies' area for assemblages' S, ∆+, and species composition, confirming the results of investigations performed on the benthic fauna of high-altitude ponds [22].
The best MAM included as predictors the degree of non-arboreal vegetation cover of the land areas neighboring the waterbodies, temperature seasonality, and mean temperature of the wettest quarter. The positive influence of the non-arboreal vegetation cover on species richness is consistent with the results of studies that were focused on macrobenthos in lotic habitats (e.g., [73] and literature cited). This could be ascribed to a positive, indirect effect of lateral trophic enrichment on aquatic primary producers, increasing zooplankton diversity by a phytoplankton-mediated bottom up effect or, alternatively, promoting habitat heterogeneity through an increase in aquatic vegetation [66,74,75].
Noticeably, the two temperature variables had contrasting effects on species richness: the lowest number of species was predicted to occur in waterbodies that were subjected to minimum temperature variability during the year and to maximum temperatures during the wettest months, i.e., in winter. The positive influence of temperature variability on patterns of species diversity has been acknowledged for zooplankton and other freshwater invertebrates [76,77], and can generally be ascribed to an effect of habitats environmental heterogeneity on empty niches availability, and, in turn, on species richness [78]. The negative influence of the mean temperature of the wettest quarter (MeanTwet) on species richness can be explained while considering that the former scales negatively with the altitude of the water bodies (r = −0.95, P < 0.001, 38 d.f.). Thus, MeanTwet maximum values were observed for basins that were located at low altitudes, generally in highly anthropized areas. Accordingly, the lower species richness characterizing these environments might be actually determined by the interplay of a spectrum of anthropogenic perturbations, such as pollution, water caption, or introduction of predatory fish. Indeed, until 1990, the vast majority of low-altitudes waterbodies in Albania have been stocked with both native and non-indigenous fish species [79], and fish predation has been repeatedly recognized to influence the species richness of zooplankton in lentic habitats [80,81].
The "assemble first, predict later" modelling strategy that was used in the present study was successful in predicting species richness and is generally acknowledged to have several advantages, among others an enhanced capacity to synthesize complex data into a form more readily interpretable by scientists and decision-makers [82]. However, it is apparent that it presents important limitations, as testified by the failure in modeling taxonomic distinctness. ∆+ resulted in being not predictable, confirming the results of several investigations that have found weak or negligible relationships of the taxonomic distinctness of macrobenthic communities with environmental factors [27,83,84]. Species richness S and taxonomic distinctness ∆+ are not conceptually (or mechanistically) related, and they behave differently [44]. The lack of congruence between S and ∆+ and the negligible predictability of the latter is because S is likely to respond to short-term environmental changes in the waterbodies under analysis, while ∆+ is a proxy for phylogenetic diversity, reflecting a complex set of intrinsic and extrinsic traits and expressing evolutionary long-term adaptations to local environmental conditions [85]. The canonical analysis of principal coordinates allowed for us to partially overcome these limitations, applying an "assemble and predict together" strategy [82] in order to model changes in the species composition of planktonic assemblages and provide an advanced resolution of species-specific relationships with bioclimatic factors. The CAP analysis ( Figure 4) distinguished two distinct group of waterbodies, showing different climatic characteristics in terms of isothermality, mean temperature of the wettest quarter, and annual precipitation. The first (red circles in Figure 4) was mainly constituted by high-altitude ponds and lakes (elevation 1725.3 ± 123.7 m, mean ± 1SE) distributed throughout the study area characterized by Crustacea species (e.g., Eucyclops serrulatus, Mixodiaptomus tatricus) that are peculiar of pristine alpine environments [86]. The second group (blue circles in Figure 4) comprised low-altitude karst waterbodies that were located in the Dumre area (Figure 1; elevation 239.9 ± 86.2 m, mean ± 1SE), where they are subjected to several anthropogenic pressures, including agricultural and urban pollution, eutrophication, and the introduction of non-indigenous fish species [87]. Accordingly, the group is characterized by the occurrence of Neodiaptomus schmackeri, an Australasian species of Chinese origin that was recently introduced in Albanian lentic habitats through fish stocking [88] and the copepod Ergasilus sp., whose adult females are ectoparasites of fish [89].
Regarding the three isolated waterbodies in Figure 4 (i.e., JAH, PIL, and RRE), they are artificial reservoirs located at 1330, 701, and 280 m a.s.l., respectively, being heavily affected by cattle frequentation (Belmonte, personal observation). Their isolation in the CAP diagram is due to the low species richness (three species in all the waterbodies), that might be ascribed to cattle-induced eutrophication conditions [90]. However, it is worth noting that copepods vary their body size (at the community level and even for single species) inversely with the eutrophication level [91,92]. Thus, by using a plankton net with a mesh size of 200 µm, we may have underestimated the smaller component of the planktonic assemblages thus biasing the species count.
A final consideration deserves a brief mention. In this study, the spatial resolution of the bioclimatic layers (approximately 0.64 km 2 ) was lower than the area characterizing most of the ponds and lakes included in the analysis (Figure 2). In other words, the bioclimate spatial grid "matched" the dimensions of the waterbodies, the latter being completely included (and described in terms of climate and vegetation cover) within the same grid cell. Additional studies including a wider size range of waterbodies as well as bioclimatic layers resolved at different spatial resolutions are necessary to provide a more complete picture of the actual relationships linking bioclimatic factors and the diversity of lentic zooplankton at multiple regional and environmental scales. Acknowledgments: Many collaborators allowed, during 12 years, the collection of samples. Among them we are particularly grateful to Bilal Shkurtaj (Univ. of Vlore, Albania), and Giuseppe Alfonso (Univ. of Salento, Italy). The authors finally thank two anonymous reviewers whose comments greatly improved the original manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. List of the 25 bioclimatic variables used in the study. Layers coded MeanT-PrecColdQ were obtained from the Wordclim v2 dataset [41] available at http://www.worldclim.org/, and refer to average monthly climate data relative to the period 1970-2000 with a 30 arc-second spatial resolution (approx. 0.92 × 0.70 km within the study area). The vegetation layers %tr and %nontr (% tree cover and % non-tree cover) were extracted from the MODIS Vegetation Continuous Field (VCF) product (https://lpdaac.usgs.gov/products/mod44bv006/; [93]). In the table, a third VCF variable-% bare soil (%bare, in italics)-is included, as it was used together with variables %tr and %nontr only for the estimation of habitat heterogeneity (see text for further details). The R packages raster, rgdal, and maptools were used for the manipulation of bioclimatic layers [94][95][96].