Suggested Sampling Methodology for Lake Benthic Macroinvertebrates under the Requirements of the European Water Framework Directive

: The estimation of the number of samples required for reliably monitoring lakes’ benthic macroinvertebrates is difﬁcult due to the natural variability and cost and time constraints. To de-termine a statistically robust and effective sampling design, we collected benthic macroinvertebrate samples from 15 Greek natural lakes. We compared the spatial and temporal variability of the benthic macroinvertebrate community composition to identify differences among lakes, between lake zones (sublittoral and profundal) and sampling periods. Furthermore, we examined the sampling precision and determined the number of required samples to attain maximum taxa richness. The diminution of the sampling effort was estimated and the desired precision level, considering different benthic macroinvertebrate abundances, was modelled. No temporal or spatial variation between lake zones was observed in communities’ compositions. The precision of our sampling design was adequate, and rarefaction curves revealed an adequate taxa richness (>70%). The developed model could be applied to assess the required sampling effort in lakes within the Mediterranean ecoregion with similar benthic macroinvertebrate abundances.


Introduction
The implementation of the European Water Framework Directive (WFD) [1] requires the use of different Biological Quality Elements (BQEs), including phytoplankton, macrophytes, fish, phytobenthos, and benthic macroinvertebrates, to assess the ecological quality of surface water bodies. Among these, benthic macroinvertebrates have been widely used as biological indicators for more than a century, as they meet the requirements of an effective ecological indicator due to their ability to respond to different pressures by changing their abundance, taxonomic richness, community composition, and biological traits [2]. These responses should be included and monitored in the framework of monitoring programs implemented by all European Member States (MSs). To fulfill the WFD objectives, the sampling effort to acquire biological data should be sufficient and feasible. Particularly, the number of samples should be representative and precise [3]; otherwise, it may influence the biomonitoring results. Thus, a sampling campaign needs to be designed with additional requirements for evaluating uncertainty [3]. At the same time, it should be cost and time effective [4].
Lake monitoring efforts usually focus on how benthic communities respond to different stressors in different lake zones [5]. The outcome could further use to develop ecological quality indices (i.e., [6][7][8][9][10]). Thus, monitoring assessment tools based on littoral benthic fauna is mainly related to the impacts of anthropogenic shoreline alterations (e.g., [6][7][8]), whereas sublittoral and profundal benthos is related to eutrophication (e.g., [5,9,10]). The developed indices are mainly expressed in numerical terms via metrics derived from raw abundance data [11] or encompassing, directly or indirectly, taxa richness [12]. Even if the sampling error is minimized as the number of samples increases, these indices display sample dependency [13]. The number of samples depends on having an appropriate spatial sampling scale and adequate sample size [13]. In biomonitoring programs, though, assessing the benthic macroinvertebrate sampling effort is often neglected [12].
The design of a sampling campaign should overcome the spatial and temporal variation in benthic macroinvertebrate community composition [14][15][16][17][18] that is affected by, among other factors, pollution and habitat characteristics, such as the water depth and substratum type [19,20]. Conventionally, lake benthic macroinvertebrate communities can differ among habitats and sites within or among lakes [19]. Moreover, several studies have shown that the number of sampling sites reflects the habitat heterogeneity and the cost and time constraints included in the sampling effort (e.g., [4,11,15]). An increase in the number of samples is usually associated with a low standard error and high sampling precision [21]. Studies based on only a few sites may fail to estimate the species distribution variability [13,16,22] and dissimilarities in assemblages [23]. Ideally, to achieve an accurate dataset, an appropriate number of sites should be sampled, which, in most cases, is high. In practice, this cannot be achieved as most of the routine monitoring programs are constrained in terms of time and budget regarding field sample collection and taxonomical processing of macroinvertebrates. An optimal sampling design should fit within such restrictions [24], as inadequate sampling is likely to lead to an invalid and low quality assessment and subsequent incorrect or limited management measures [12].
Several statistical approaches have been applied so far to estimate an adequate sampling effort [4]. Generally, a sufficient number of samples depends on the mean abundance, the degree of aggregation, and the desired precision [25]. The degree of aggregation, per se, further depends on the size of the sampler [26]. In the case of any sampler, the number of replicates affects the accuracy and precision of the statistics [25]. Three replicates are usually considered the minimum sample size to perform statistical analysis [27]. Taylor's law has been widely applied to optimize the sampling design and statistical interpretation of population data [28,29], using a power function between variance and mean data [26,30,31].
To the best of our knowledge, studies on the assessment of the optimal sampling effort of benthic macroinvertebrates in lakes, from the perspective of the WFD, have been narrowed in the Northern and Central Baltic ecoregions, with their majority focusing on the littoral benthic macroinvertebrates (e.g., [13,15,25]). Considering the natural Greek lakes, such research is entirely lacking for all the BQEs, apart from fish [32]. Thus, in the present study, we developed a model to estimate the optimal sampling effort required for routine monitoring of benthic macroinvertebrates in 15 natural lakes under the requirements of the WFD. We evaluated the minimum number of samples for the desired level of precision concerning a given benthic abundance. The precision and the sampling effort of benthic assemblages were examined in the sublittoral and profundal zones as well as in both zones simultaneously. Since the selected natural lakes represent the most frequent lake types present in Greece and also in most Mediterranean countries, the constructed model could (a) be further applied within the ecoregion's lakes with similar mean benthic macroinvertebrate abundances and a desirable precision level and (b) contribute to the establishment of routine monitoring programs in compliance with the objectives of the WFD.

Study Area
The Greek National Monitoring Network comprises 24 natural lakes (lake area > 0.5 km 2 ). Among them, we collected benthic macroinvertebrate samples from 15 lakes (63% of the total number) ( Figure 1); three of them are transboundary (Doirani, Megali Prespa, Mikri Prespa). We examined two types of natural lakes according to the Hellenic typology [33]: GR-SNL-shallow natural lakes (mean depth: 3-9 m, n = 8, 100% of the total number) and GR-DNL-deep natural lakes (mean depth >9 m, n = 7, 100% of the total number), which are located at altitudes ranging from 16 to 853 m.a.s.l. and differ in their limnological characteristics (Table 1). The total phosphorus concentrations varied among lakes, ranging from 7 µg/L in Lake Kourna to 670 µg/L in Lake Zazari (Table 1). Regarding land uses in their catchment area (Corine Land Cover 2018), artificial surfaces cover generally in low percentages the catchment, while agricultural areas seemed to be the major cover in many lakes (Table 1).  Table 1.

Sampling Design and Sample Processing
An adequate inventory for monitoring benthic macroinvertebrates in lakes requires biannual samplings [35]. Thus, we initially sampled in the mid-spring of 2014, before the monomictic lakes underwent the summer thermal stratification [35]. This sampling period is also recommended by the National Monitoring Program [36]. Additional samplings were conducted in autumn, when a lower water level was recorded due to intense evaporation [37][38][39]. Generally, it is suggested that winter samplings should be avoided to prevent sampling under extreme conditions [40]. In each lake, the soft sediments from the sublittoral and profundal zones were sampled using an Ekman-Birge grab (sampled area 225 cm 2 ), suitable for sampling both shallow and deep water bodies [41]. Three replicates were collected at each site to provide a measure of site variability [35]. The sampled substrate was composed of clay, silt, or sand fractions. The number of sampling sites was dependent on the maximum depth and each lake's specific characteristics (e.g., slope, substrate), following a constant sampling strategy from the sublittoral (i.e., started at depth > 2 m and extended until macrophyte growth is not possible anymore and the profundal zone starts) and profundal zones ( Table 2). Sediment samples were sieved with a 500 µm sieve. All benthic macroinvertebrates were sorted and identified at the lowest possible taxonomic level (genus or species); taxa abundances are expressed as ind/m 2 . Samples from the profundal and sublittoral zones were kept separate, and the mean abundance of the two sampling periods was calculated per lake zone.

Spatial and Temporal Variability
Non-metric multidimensional scaling (MDS, PRIMER, Version 6, PRIMER-E Ltd, Plymouth, UK) [42] was applied to visualize similarities in benthic macroinvertebrate community structures between lake zones (sublittoral and profundal), sampling periods (spring and autumn), and among lakes. A one-way analysis of similarities (ANOSIM, PRIMER, Version 6, PRIMER-E Ltd, Plymouth, UK) [42] was performed to identify significant differences in the macroinvertebrate community composition between lake zones, sampling periods and among lakes using 9999 permutations. In this analysis, high R values (>0.75) indicate that groups are well separated from each other, whereas low values (R < 0.25) imply little or non-existing differences between groups [43]. All of the above analyses were based on a Bray-Curtis similarity matrix of log(x + 1) transformed abundance data.
The mean abundances of benthic macroinvertebrates for each zone and sampling period were tested for normality using a Shapiro-Wilk test (SPSS version 21, IBM Corp., Armonk, USA). As this assumption was not met, the non-parametric Kruskal-Wallis test was applied to examine the differences between lake zones and sampling periods (SPSS version 21, IBM Corp., Armonk, USA). When the above statistical tests indicated no significant differences, all further analyses were performed on data from the combination of sublittoral and profundal zones, hereafter referred to as sublittoral/profundal zone.

Precision
The level of precision was estimated using the following equation [44,45]: where se is the standard error of the abundance of benthic macroinvertebrates and mean is the mean abundance (ind/m 2 ) of benthic macroinvertebrates. A bootstrapping analysis was applied to assess the effect of the sampling effort on the precision (P) of benthic macroinvertebrates in the sublittoral/profundal zone. The mean values of the observed data were generated (1000 replacements) under an increasing sampling effort (ranging from two up to the maximum number of samples) using the R-studio package sciplot [46].

Adequacy of the Sampling Effort
Sample-based rarefaction curves were generated to estimate the adequacy of the sampling effort by randomly subsampling the entire community 1000 times to figure out the number of taxa as a function of the sample size [47]. Additionally, the expected total taxa richness for each lake was estimated using the first-order Jackknife estimator, as it is considered appropriate for abundance data [48]. The above analyses were performed using the statistical software EstimateS 9.1.0 (Colwell, http://purl.oclc.org/estimates, accessed on 9 May 2021) [47].

Modelling
The optimal sampling design was estimated as a function of the mean abundance of benthic macroinvertebrates and the desired precision level [44]: where n is the number of samples, s 2 is the sample variance, mean is the mean abundance of benthic macroinvertebrates per lake, and P is the precision. The model was generated by combining Taylor's power law [28] with a precisionbased estimate of the essential number of samples. Taylor's power law estimates the variance (s 2 )-mean (m) relationship, based on the equation: Data that failed to meet the following criteria were excluded from the analysis: estimates of mean and variance based on less than 15 samples (N < 15), regression based on less than five sample pairs (M < 5), and a range of the log(mean) values less than one [49,50]. Parameters a and b were estimated after the log-transformation of Equation (3): The b values, which represent the slope of the regression, are referred to as an index of the aggregation [29], where b > 1 represents an aggregated population, b = 1 indicates a random distribution, and b < 1 indicates a uniform distribution. Equation (4) was untransformed and corrected for bias [51]: where MSE is the mean square error of the regression.
The number of samples required for a desired level of precision can be estimated by combining Equations (2) and (5):

Taxa Richness
In total, 93 taxa were recorded from all lakes, with a mean (±SE) of 20.20 (±2.32) taxa per lake. The highest number of taxa (36 taxa) was observed in Lake Vegoritis and the lowest (seven taxa) in Lake Kourna (Figure 2a). In all lakes, apart from Lake Mikri Prespa, the sublittoral zone contributed to the total taxa richness with one to 22 taxa (Figure 2a). The lowest mean abundance (72 ind/m 2 ) was recorded in the profundal zone of Lake Amvrakia and the highest (8274 ind/m 2 ) in the profundal zone of Lake Kastoria. A detailed taxa list is given by Ntislidou et al. [10]. The most abundant family was Chironomidae (27 species), followed by Tubificidae (24 species). Potamothrix hammoniensis (Michaelsen, 1901) and Chironomus gr. plumosus (Linnaeus, 1758) were the most frequently recorded species (87% and 73.3%, respectively).

Spatial and Temporal Variability
The MDS analysis revealed no spatial or temporal differences between lake zones concerning macroinvertebrate community composition (Figure 3a,b). The ANOSIM comparing the benthic fauna of lake zones (sublittoral and profundal) and seasons (spring and autumn) showed that R values were relatively low (R = 0.15 and R = 0.02, respectively). Among lakes, benthic macroinvertebrate community composition differed (R = 0.45, Figure 3c). Beyond these analyses, the Kruskal-Wallis test showed no statistically significant differences in benthic macroinvertebrate communities between lake zones (p = 0.852) or sampling periods (p = 0.375). Such outcomes indicate that an optimal sampling effort would combine both sublittoral and profundal zones.

Precision
Regarding the sampling effort at the sublittoral/profundal zone, the highest mean abundance of benthic macroinvertebrates was recorded in Lake Kastoria (mean abundance (±SE) 7134.11 ± 832.4 ind/m 2 ), while the lowest was recorded in Lake Ozeros (mean abundance (±SE) 135.67 ± 18.38 ind/m 2 ) (Figure 4). The sampling effort was less than the threshold of 0.2 for all cases, and the mean precision value (±SE) was 0.14 (±0.01) (Figure 4). The bootstrapping analysis showed that the level of precision was increased by increasing the sampling effort ( Figure 5). Additionally, in deep lakes (mean depth > 9 m), apart from Lake Megali Prespa, a greater sampling effort is required to obtain more precise estimates than in shallow ones ( Figure 5).

Adequacy of Sampling Effort
No asymptotes, except for the cases of the lakes Megali Prespa and Amvrakia, were evident in the rarefaction curves based on the total sampling effort ( Figure 6). The firstorder jackknife estimator showed that the observed richness represented >70% of the estimated total richness in all lakes ( Figure 6). However, 90% of the species, which is approximately three taxa less than the total observed, was obtained in 73.9% of all cases. Figure 6. Taxa rarefaction curves at the sublittoral/profundal zone. Observed (Sobs) and estimated (Sexp) taxa richness based on the first-order Jackknife estimator. The gray section represents the ±95% confidence interval.

Model Development
The model developed for optimal sampling effort, considering a given mean abundance and the desired level of precision, was based on nine pairs of variances (s 2 ) and the mean abundance of benthic macroinvertebrates. The log(s 2 )−log(mean) regression analysis was statistically significant (p < 0.0001) and explained 86% of the variance ( Table 3). The optimal sampling effort can be calculated using the following equation: where mean is the mean abundance of benthic macroinvertebrates in the lake, and P is the level of precision. The above equation shows that the number of samples increases as precision is increasing (i.e., lower values of P) (Figure 7), and it can be used to calculate the required number of samples for different abundances and precisions (Table 4). For instance, for the mean abundance of benthic macroinvertebrates of 2000 ind/m 2 , five samples are required to attain a precision level of 0.2. Table 3. Parameters of the linear regressions between the sample variance (log(s 2 )) and the mean abundance of benthic macroinvertebrates (log(mean)).  to attain the desired level of precision (P).

Discussion
The cost and time constraints are crucial in the design of monitoring programs for assessing the ecological water quality of surface waters under the requirements of the WFD. Frequently, it is necessary to lessen sampling effort and subsequently process samples in a laboratory while maintaining precision in the estimates. The investigation of optimal sampling designs is an essential part of the decision-making procedures for both routine monitoring and ecological studies. Unfortunately, similar studies have overlooked these factors in recent years, and so have tended to draw inaccurate conclusions [24]. The present study highlights a useful process for estimating precise routine monitoring sampling designs in lakes using benthic macroinvertebrates, as is required by the WFD.
The design of lake monitoring samplings using benthic macroinvertebrates generally follows a habitat-specific procedure to overcome natural heterogeneity, especially in the littoral zone [13,52,53]. However, it has been suggested that mesohabitats or proportional composite samples should be applied in the littoral zone as a cost and time efficient approach [15,52]. Sublittoral and profundal zones, which are considered "simple" habitats, mainly consist of soft-bottom substratum with low structural or background variability [54,55]. Our outcomes across a gradient of different morphological lakes indicate that profundal and sublittoral zones do not differ significantly, possibly due to their similar depths in general. In other studies, the composition of benthic macroinvertebrate communities in sublittoral and profundal zones is considered as a united subsystem as they seem to respond to the same pressures [5,53,56].
Temporal variability also prevents monitoring surveys and is dependent on the time of the year at which the study is conducted (e.g., [57,58]). The changes in the compositions of benthic macroinvertebrate communities are influenced, among other things, by their life cycles, environmental conditions, and competition [14]. The present study showed that benthic fauna of sublittoral/profundal zones in shallow and deep Greek lakes was similar in spring and autumn, indicating low temporal variability. No temporal variability in the composition of lake benthic communities is also referred [59][60][61][62]. Schreiber and Brauns [13] suggested using spring samples to provide a representative view of macroinvertebrate communities due to the negligible compositional dissimilarities between spring and autumn.
The optimal number of samples is the largest achievable; in practice, this usually is constraint by cost and time factors. Thus, it is crucial to evaluate their number that will provide an acceptable level of precision a priori. In general, the precision threshold used for benthic macroinvertebrates is 0.2 [21,26], which means 95% confidence limits of the mean values are ±0.4 or more [25]. In our study, the precision values were lower than this threshold for the sublittoral/profundal zones, revealing that the precision of sampling design was feasible. The bootstrapping analysis indicated that precision is minimized when the sampling effort is reduced and a more intense effort is required in the deep lakes. Precise estimates can also be obtained by reducing the sampling effort by up to two samples in the shallow lakes.
Rarefaction curves are a useful tool for estimating the sampling effort, and thus, they may be helpful for sampling campaigns and monitoring programs. The lack of an asymptote in the rarefaction curves indicates the presence of sporadic species in the samples. The underestimation of taxa richness in lakes, as in rivers, is very common in many cases [12,13,[63][64][65]. In freshwater ecosystems, benthic macroinvertebrates are typically a highly diverse group [12] and tend to be dominated by rare species [66,67]. However, it is considered sufficient to record >70% of the estimated species [68], as in the present study. In most cases, the curve slope was smooth, indicating that 90% of the total taxa richness was achieved at 50% of the sampling effort. Steeply increasing curves indicate that more taxa will be recorded with a higher number of samples. However, a larger number of samples is required for species inventories, but it is questionable whether rare species play a crucial role in assessing freshwater ecosystems [69].
Regarding the developed model, parameter α of Taylor's law is related to sampling conditions and depends on the size of the sampling unit [44]. In this study, the α value (0.399) was lower than the values presented in other studies, i.e., 1.515 [25] and 1.300 [26] (log(x + 1) transformed data), possibly due to the larger sampling area (225 cm 2 × three replicates) covered in the present study. Parameter b of Taylor's law is a distribution factor for populations present in specific environments, following certain sampling procedures [70]. In this study, the estimated b value indicates that benthic macroinvertebrates of Greek lakes have an aggregated distribution (b > 1). Additionally, this value (1.649) was higher than those from lakes in Northern Europe (b = 1.25) [25], but it was within the range of values observed worldwide (1.18-2.38) [26]. The b value depends on the pairs' variance-mean and the number of samples [45].
It is crucial to know the number of samples required to give an adequate level of precision [71]. In our case, 22 samples were required to collect 30 ind/m 2 , while other studies mentioned more (62 samples) [25] or less (13 samples) [26] samples for a precision level of 0.2. This difference may be due to the higher abundance of benthic macroinvertebrates recorded in this study (range: 178-8135 ind/m 2 ) compared to others (range: 205-770 ind/m 2 ) [25]. However, this difference is minimized with increasing abundance. Increasing the number of samples usually improves the level of precision, but when the abundance is low (<10 ind/m 2 ), the maximum achieved precision is low [45]. A model of optimal size could be used for monitoring programs or pilot studies to ensure that there is a tolerable level of precision to provide cost and time efficient research.

Conclusions
The WFD requires the assessment of the ecological quality of lakes using distinct BQEs, among which are benthic macroinvertebrates. However, the application of this type of assessment is hindered by a lack of information on sampling protocols [72]. Additionally, extensive sampling campaigns increase the necessary cost and time effectiveness for sample collection in the field and sample processing in the laboratory. This study proposes a simple method of estimating the number of benthic macroinvertebrate samples required at four different levels of sampling precision for the sublittoral/profundal zone. Consequently, our framework for calculating the sampling effort may improve the efficiency of routine monitoring programs in other Mediterranean lakes of the same limnological characteristics.
Its applicability to other lake types, such as reservoirs or large oligotrophic lakes from other ecoregions, should be treated with caution as richness patterns may not be the same. In such lakes, the sampling effort may be higher than in the studied lakes. Additionally, considering samplings conducted for other purposes, such as the restoration or the surveying of rare or endangered species, habitat samplings are recommended.