Does littoral substrate affect macroinvertebrate assemblages in Mediterranean lakes?

The objective of this study was to investigate the effects of substrate type in macroinvertebrate assemblages in Mediterranean lakes. Samplings have taken place in the littoral zone of 21 lakes in Greece, between 2015 and 2018. We compared benthic macroinvertebrate assemblages among three substrate types of their littoral zones; sandy, covered with macrophytes and stony substrate. Benthic macroinvertebrate assemblages at sites with extended macrophyte cover differed only slightly in composition and abundance from the ones found in stony and sandy substrates. Coenagrionidae were indicative of sites covered with macrophytes and Oligochaeta and Erpobdellidae were representative of stony substrates. The type of substrate proved to be a statistically significant factor influencing the number of benthic macroinvertebrate taxa, the relative abundance of Oligochaeta and the relative abundance of Odonata. In the context of designing site-adapted management measures, priority could be given to the conservation and restoration of aquatic vegetation in lake littoral zones, which host rich macroinvertebrate assemblages with abundant taxa of Odonata.


Introduction
Freshwater biodiversity is worldwide declining at unprecedented rates (IPBES 2019). Factors underlying this biodiversity loss have been thoroughly studied over the years and can be categorized in six groups: hydrological alterations, habitat degradation and loss, pollution, overexploitation, invasive species and climate change (Dudgeon et al. 2006;Reid et al. 2019;Arthington 2021). Most of these factors are linked to human activities, which resulted in many restoration activities being carried out, focusing on mitigating these anthropogenic impacts (Gething et al. 2020). However, in many cases there are knowledge gaps that hamper restoration efforts. For example, although it is well known that lakes are subject to morphological and hydrological modifications, the effects of these changes have been less understood (Poikane et al. 2020a).
In Europe, we have increased our knowledge of freshwater ecosystems structure and function over the years, owing mainly to the implementation of the Water Framework Directive (WFD). In this context, benthic macroinvertebrates have been among the most widely used biological quality elements for ecological assessments purposes especially in rivers and lakes, as reported in the 2018 EC Intercalibration Decision (European Commission 2018). Recently, they have become more popular in bio-assessment methods as they have been proved sensitive not only to eutrophication and acidification but also to morphological changes and general degradation (Urbanič et al. 2012;Poikane et al. 2016Poikane et al. , 2020bMavromati et al. 2021). Ntislidou et al. (2021) characterize them as "ecosystem engineers" playing an important role to aquatic ecosystem services as they participate in various biogeochemical processes. In the Mediterranean region, one such WFD-compliant assessment system based on littoral benthic macroinvertebrates of Greek lakes has recently been developed (HeLLBI; Mavromati et al. 2021). HeLLBI was designed and applied to respond to anthropogenic impacts, and in particular eutrophication and morphological changes (artificial shoreline).
The distribution of macroinvertebrate communities in freshwater ecosystems is affected by both natural and human factors (Dou et al. 2022). The natural factors include all these environmental parameters that characterize a single site or water body, such as temperature, water depth, dissolved oxygen, pΗ, conductivity, the spatial heterogeneity of habitats (Free et al. 2009;Dou et al. 2022). Human activities affect macroinvertebrate communities either by nutrient enrichment caused by land use activities or due to shoreline modification (Ntislidou et al. 2018; Bartels et al. 2021;Mavromati et al. 2021). Differences in biological communities across various environmental parameters have been extensively studied in streams and rivers but in lakes, especially in the Mediterranean region, these studies are limited (Vinson and Hawkins 1998).
Biological diversity can be studied several ways; it can be characterized by the diversity of species within sites and quantified by measuring the number of taxa or by diversity indices (Costa and Melo 2008). Another way of measuring diversity is by describing the dissimilarities of biological assemblages between different environments (Costa and Melo 2008).
As bio-indicators, freshwater macroinvertebrate taxa respond differently to stressors and they have the ability to incorporate the effects of the stressors they are exposed to, in combination and over time. As they show spatial variation along different lake zones and depths, habitats and even lakes, finding out which natural factors affect their assemblages is essential in order to separate them from human pressures and their confounding effects (Solimini et al. 2006). As a rule of thumb, waterbodies in good ecological status could support a variety of species-indicators for undisturbed waters with high abundance values.
Except from their response to human stressors, freshwater macroinvertebrate taxa vary in response to habitat heterogeneity and substrate type. Generally, high habitat heterogeneity and coarser substrates present great invertebrate diversity and abundance (Zenker and Baier 2009). According to Graça et al. (2015), coarse substrates are more complex than fine ones, as they retain coarser particulate organic matter and abundant with microalgae with higher biomass, which can be used as food resources. Most studies accept the hypothesis that stony substrates present higher richness and abundance of benthic taxa and in particular abundance of collector-gatherers and filtercollectors (Pereira et al. 2017).
Additionally, the presence of macrophytes in the littoral zone results in higher habitat heterogeneity. Dense macrophyte stands could alter invertebrate communities in numerous ways. Not only do they serve as food resources directly and indirectly through the growth of periphyton, but they can also access and deplete sediment nutrient loads (Waters and Giovanni 2002;Tolonen et al. 2003;Zenker and Baier 2009;Salmon et al. 2022).
The overall objective of the study is to investigate the effects of substrate types in benthic macroinvertebrate assemblages in Mediterranean lakes differing in trophic status and other environmental parameters. In particular, we want to investigate the differences in benthic macroinvertebrate assemblages among the sandy substrates, substrates covered with macrophytes and stony substrates in Mediterranean natural lakes, identify the indicator taxa characterizing each substrate type and explore the role of macrophytes in shaping these macroinvertebrate assemblages.

Materials and methods
Study area and sampling procedure Benthic macroinvertebrates data were collected from 21 lakes (Fig. 1), including two transboundary lakes, Megali and Mikri Prespa, all belonging to the Greek National Water Monitoring Network (Mavromati et al. 2021).
The lakes included in this study belong to three different natural lake types, according to the mixing regime and depth gradient (Kagalou et al. 2021). As most Mediterranean lakes, they face multiple pressures including nutrient loading from point and non-point sources, water abstraction and morphological changes (Latinopoulos et al. 2016). Moreover, they seem to be affected by regional landscape characteristics, in comparison with most cold temperate and tropical lakes, by groundwater hydrology and by the Mediterranean climate (Alvarez Cobelas et al. 2005;Mavromati et al. 2018).
In total, 97 littoral sampling sites have been surveyed in the spring season during the 2015-2018 sampling campaign. Most lakes were sampled once. Six lakes were sampled for two years and one was sampled for three years. A nonparametric Kruskal-Wallis test was applied in this dataset and showed that there were no differences among the sampling years (Mavromati et al. 2021). The number and location of the sampling sites for each lake were selected according to lake size, habitat maps, and land use data of the lakes and their catchment areas.  Ismarida,5: Kastoria,6: Koroneia,7: Kourna,8: Lysimacheia,9: Megali Prespa,10: Mikri Prespa,11: Ozeros,13: Paralimni,14: Petron,15: Stymfalia,16: Trichonida,17: Vegoritida,18: Volvi,19: Voulkaria,20: Yliki,21: Zazari Samples were collected using a semi-quantitative approach, which consists of a three-minute kick/ sweep with a standard hand net (500 μm mesh size), at the littoral zone of each lake (up to 1.2 m depth of water). This particular approach was selected to cover potential effects of lakeshore modifications and it was conducted during the development of the littoral macroinvertebrate assessment method (Mavromati et al. 2021). At each sampling site, the cover with aquatic macrophytes (e.g. Phragmites australis, Potamogeton sp.) was recorded as a percentage. Furthermore, visual assessments of substrate composition were made based on the predominant substratum size using the size categories given in the Wentworth scale (Wentworth 1922). The substrate composition was further grouped into two categories: sandy substrate (< 2 mm diameter) and stony substrate (> 2 mm diameter). As a result, each site sampling site was categorised into one of the three substrate types: stony substrate (> 2 mm diameter), macrophyte and sandy substrate (< 2 mm diameter).
Sieving was carried out on site; sorting, identification and counting were carried out in the laboratory, and the samples were preserved in vials containing 70% ethanol. The littoral invertebrate fauna was identified to family level, except oligochaetes, which were identified as a subclass. Finally, the total number of individuals of each taxon was recorded and their relative abundance was calculated.

Statistical analysis
The nonparametric Kruskal-Wallis test was applied to examine statistically significant differences of macroinvertebrate community indicators between different substrate types (stony substrate, macrophytes and sandy substrate), followed by Wilcoxon signed rank pairwise tests for the statistically significant indicators. Nonparametric tests were chosen, as the data did not meet the assumptions of parametric tests. Boxplots were used to show distributions of numeric values of statistically significant biological community indicators. Boxplots were prepared in R using the ggplot2 function (R Core Team 2018).
Non-metric Multidimensional Scaling (NMDS) is an ordination technique, which enables complex multivariate data to be visualised in two dimensions. This technique was employed to check for differences in benthic macroinvertebrate assemblages with samples being a priori grouped by substrate type. NMDS centroids were calculated as the centre points of all replicates for each sampling method in multidimensional space, as Gething et al. (2020) suggested, using Bray-Curtis similarity coefficients on fourth root-transformed data. According to Anderson et al. (2008), data transformation has been recommended as a way to reduce the contribution of highly abundant species in relation to less abundant ones in the calculation of Bray-Curtis measure; when the transformation is severe (e.g. fourth root), rare species will have higher contribution to the analysis (Clarke et al. 2014). This test statistic calculates a pseudo-F value, similar to the F value in ANOVA. Larger pseudo-F values indicate more pronounced group separation; however, its significance is usually of more interest than its magnitude. The ordination diagram was also enhanced with the convex polygons around each centroid (Costa and Melo 2008). Differences in the benthic fauna between the three different substrate types were statistically tested via a permutational multivariate analysis of variance (PERMANOVA) using the adonis2 function in vegan (Oksanen et al. 2012;Arbizu 2020;Gething et al. 2020).
In order to analyse the multivariate homogeneity of group dispersions (variances) of macroinvertebrates assemblages and characterise their variability, the betadisper function in vegan package was calculated based on Bray-Curtis distances (Oksanen et al. 2012). Ordination plots were prepared and statistical analyses were performed in R version 3.5.3 (R Core Team 2018).
The Similarity Percentages Analysis (SIMPER) was performed in Primer v7 software, to examine which taxa contributed most to the average similarity of sampling sites within each substrate type and the average dissimilarity between different types of substrate (Clarke and Gorley 2015). Το identify dominant taxa within each macroinvertebrate assemblage for each substrate type, the indicator value index was applied using the indval function within the ladsv package and the indicators function within the indicspecies package (De Caceres et al. 2016). The former is performed by a permutation test to assess the statistical significance of the association between species and site groups, yielding a percentage indicator value for each species (De Caceres and Legendre 2009;Legendre and Legendre 2012). An indicator value of 0.25 was accepted as being ecologically relevant (Dufrene and Legendre 1997), and all significant indicators with a fidelity value below 0.25 were removed to exclude rare taxa (De Caceres et al. 2012). All analyses were performed with the use of ladsv (Dufrene and Legendre 1997), indicspecies (De Caceres and Legendre 2009), and tidyverse (Wickam 2017) R packages in R environment version 3.5.3 (R Core Team 2018).

Results
The dataset used in this study included 77 taxa (together with Oligochaete subclass). The taxonomic groups with higher numbers of taxa were Diptera (13), Gastropoda (10), Coleoptera (10) and Odonata (7). The highest number of taxa were observed in Lakes Lysimacheia and Vegoritida (23 and 21 taxa, respectively) and the lowest in Lake Koroneia (one taxon). The lake with the highest number of EPT taxa was Vegoritida (6 taxa); on the other hand, Lake Koroneia had no taxa belonging to orders Ephemeroptera, Plecoptera and Trichoptera.
According to SIMPER, the benthic fauna of our littoral sites was dominated by five taxa: Chironomidae, Oligochaeta, Corixidae, Gammaridae and Caenidae. The results of the analysis describing the macroinvertebrates assemblages of the 97 sampling sites are shown in Table 1. Chironomidae contributed most to the average similarity of all sites; especially in sites with sandy substrate they almost represented 40% of all taxa. Benthic macroinvertebrate communities which occurred at sites with extended macrophyte cover differed only slightly from the ones found in stony and sandy substrates, in composition and abundance. Caenidae were mostly found in substrates dominated by macrophytes. The average dissimilarity between stony substrates and sites covered with macrophytes was 58.08%, between stony and sandy substrates was 58.55% and between sandy substrates and sites covered with macrophytes was 60.20% ( Table 2). The taxa that contributed to the average dissimilarity between the three substrate types were more or less the same (Gammaridae, Corixidae, Oligochaeta and Caenidae) but with different contributions between the different pairs of substrate types.
Indicator species analysis revealed two indicator taxa for sites with stony substrates (Oligochaeta and Erpobdellidae) and one indicator taxon (Coenagrionidae) for sites dominated by macrophytes (Table 3). The best combinations of indicator taxa for each substrate group are shown in Table 4; the species combination of sites with stony substrate (Gammaridae and Oligochaeta) showed high sensitivity; the majority of sites covered with stones included this particular combination. Their predictive power was rather low, meaning that this taxa combination occurred in less than half of sites belonging to this group. Chironomidae and Coenagrionidae showed also high sensitivity in sites covered with macrophytes. On the other hand, sites with sandy substrate had only one representative taxon (Atyidae) and no taxa combination; half of sites belonging to this group included this particular combination (sensitivity = 0.49).
The nonparametric Kruskal-Wallis test indicated no significant differences between sampling periods (p = 0.136; Mavromati et al. 2021). The statistical analysis showed that the substrate type had a statistically significant effect on the Number of Taxa (H = 8.233, p = 0.016), the relative abundance of Oligochaeta (H = 9.646, p = 0.008) and the relative abundance of Odonata (H = 5.997, p = 0.050) ( Table 5). The pairwise comparison for taxa richness showed that the statistically different substrate types were stony substrate-macrophytes and macrophytessandy substrate (p < 0.05). The increase in taxa richness was evident with greater substrate complexity, as macrophytes supported the greatest number of taxa while sandy substrates the fewest (Fig. 2). The relative abundance of Odonata was higher in substrates associated with macrophytes; whereas, the lowest percentages were recorded in stony substrates. The same results were shown in the pairwise comparisons as statistically significant differences were observed only between stony substrate and macrophytes (p < 0.05). On the other hand, the relative abundance of Oligochaeta was higher in stony substrates, followed by macrophytes and finally by sandy substrates. The pairwise comparison showed statistically significant differences between stones-macrophytes and stones-sand. The remaining indicators of benthic macroinvertebrate assemblages, which were tested in our study showed no statistically significant differences between the three substrate types.
According to PERMANOVA test, macroinvertebrate assemblages were significantly different  No distinct clusters were observed in the ordination space, indicating only slightly distinct benthic communities (Fig. 3). The plot of the NMDS analysis revealed three polygons that defined the maximum area of each group's site scores in the two-dimensional ordination space, which seem to overlap at a certain extent. Macroinvertebrates occurring in stony substrates seem to be a subset of those occurring in the other two substrate types. The hulls are fit to the raw data as they appear in the plot, which creates angular polygons. Despite the low degrees of freedom, pseudo-F value of PERMANOVA analysis was large enough to reject the null hypothesis of no differences in the centroid locations and/or the dispersion of groups between the three substrate types. The larger the pseudo-F value, the greater the difference is supposed to be and is different from the p value. Differences in the centroids locations of the three substrates are also evident in Fig. 4 where also the level of dispersion within each substrate is shown. These three centroid locations represent the three different substrate types and according to the boxplot, it is clear that the centroid locations of sites covered with macrophytes and those of stony substrates differ substantially.
The analysis of multivariate homogeneity of group dispersions (variances) revealed that the within-group spread of macroinvertebrate communities did not differ among substrate types (F 2, 94 = 2.3302, p = 0.095). This result clarifies the nature of multivariate effects of macroinvertebrate assemblages and suggests that differences in benthic fauna were mainly due to differences in centroid locations.

Discussion
The results presented in this paper allowed us to assess the role of the littoral substrates in benthic macroinvertebrate fauna. We found that the studied substrate types supported to a certain degree distinct macroinvertebrate assemblages resulting from different levels of habitat complexity. However, the presence of ubiquitous taxa across all sites was evident in our results and was highlighted in two different analyses, the SIMPER analysis and the nonparametric Kruskal-Wallis test (Gething et al. 2020), providing mixed results concerning community composition and diversity.
According to SIMPER results, Chironomidae showed the highest contribution to the similarity of all three groups of substrates; its percentage to the overall community structure is rather high, especially in sandy substrates. It is well documented that Chironomidae include several species with broad environmental preferences including different substrate types (Verdonschot 2006;Lencioni et al. 2018;Dorić et al. 2020). Čerba et al. (2022) concluded in their study that substrate types significantly affected Chironomidae community composition and abundance. A fauna list of Chironomidae larvae in mainland Greece highlights exactly this: most species were indicative of distinct climatic, geological and hydrochemical features (Płóciennik and Karaouzas 2013). On the other hand, Rossaro et al. (2014) refer to Chironomidae taxa as opportunistic, being found in different habitats rather than restricted in a single habitat only. They use the term "preference" to describe their ideal habitat type instead of "exclusivism" (Rossaro et al. 2014). The fauna of sandy substrate of lakes is mainly composed of Chironomidae taxa accompanied with a few other taxa, which explains their high contribution to the similarity results of SIMPER analysis. Ntislidou et al. (2021) studied the macroinvertebrate fauna of the profundal and sublittoral zones of three Greek eutrophic lakes and found that the assemblages of the muddy substrates of these zones were dominated by Chironomidae and Oligochaeta taxa.
Caenidae discriminated the group of sites covered with macrophytes from the other two groups and was found to be the taxon associated with the group dissimilarities. This result agrees neither with the study of Pilotto et al. (2015) which associates Caenidae species with stony substrates and the presence of zebra mussels which they use as food resources, nor with the results of McGoff and Irvine (2009) who found a negative correlation between Caenis luctuosa, the macrophyte Percentage Volume Inhabited (PVI) and the extent of macrophytes lakewards.
The percentage of dissimilarity was higher than 55% among all three groups of substrate types; nonetheless the taxa that contributed mostly to the average dissimilarity were the same with different percentages of contribution at each group. The IndVal analysis (De Caceres and Legendre 2009) showed more clear results and revealed indicator species for stony substrates and for substrates covered with macrophytes. No indicator species were found for sandy substrates, probably because of their unique macroinvertebrate assemblage in comparison with the other groups of sites and the relatively low number of taxa To overcome this obstacle the extended version of indicator species analysis was used to give us indicator species also for sandy substrates. The predictive values were quite similar in all three groups of sites. Sites with stony substrate highlighted Gammaridae and Oligochaeta as a combination of taxa with high sensitivity, which is an estimation of the frequency of the families at these sites. Studies indicate that species belonging to Gammaridae show a clear preference in stony substrates, especially when they are colonized by Dreissena sp. (Stewart et al. 1998;Hesselschwerdt et al. 2008). High sensitivity of the taxa combination of Chironomidae and Coenagrionidae was also found in the groups of sites covered with macrophytes. Pilotto et al. (2015) suggested that species of Coenagrionidae such as Ischnura elegans prefer sites with submerged macrophytes and natural littoral habitats.
Differences in community composition and diversity were evident across different substrate types; Gething et al. (2020) characterize substrate type as surrogates for increasing habitat complexity. In our results, sites covered with macrophytes supported greater number of taxa followed by sites with stony substrates and finally the ones covered with sand.
Our results showed that there was an increase in taxa richness in coarser substrates which agrees with the initial hypothesis of Graça et al. (2015). They analysed the macroinvertebrate community in streams and characterized fine substrates poorer in terms of abundance of macroinvertebrates and taxa richness. Surface complexity seems to be positively correlated with diversity and abundance of benthic fauna especially in studies concerning streams (Taniguchi and Tokeshi 2004;Barnes et al. 2013). On the other hand, our results showed that the relative abundance of Oligochaeta was greater in stony substrates in comparison with the other two types, which is a surprising result according to the literature (Rieradevall et al. 1999;Graca et al. 2015;Buendia et al. 2011). In particular, Graça et al. (2015) argue that finer sediments are more appropriate for taxa like Chironomidae and Oligochaeta, which prefer spending time within substrate particles. The benthic fauna of sandy substrates in our dataset was mainly composed of Chironomidae, resulting in low abundances of their remaining taxa, in comparison with the other two groups of sites. Sychra et al. (2010) suggested in their study that some species of Oligochaeta, such as Naididae with phytal preferences are quite abundant in reedbeds near the shore. Oligochaeta and Chironomidae are considered to be generally tolerant organisms and according to Mavromati et al. (2021) they are the dominant taxa in sites with high proportion of artificial shoreline in poor and bad ecological quality status in Mediterranean lakes.
Odonata larvae followed the opposite pattern and exhibited greater abundances in sites covered with macrophytes, followed by sandy and stony substrates. Waters and Giovanni (2002) associate Lestidae (Odonata) with vascular macrophytes present in depositional habitats. The same pattern was evident in Graça et al. (2015) where it is argued that Odonata and Trichoptera inhabit coarser, well sorted substrates. The links between aquatic vegetation and macroinvertebrate abundance and diversity are evident in the majority of studies but the results are still confounding with the type of sediment that aquatic macrophytes need for growth (Waters and Giovanni 2002;Zenker and Baier 2009). It is crucial to check in littoral zones which one is actually affecting macroinvertebrate assemblages; macrophytes or the sediment that they rely on? Another study focused only in Chironomidae community composition, concluded Fig. 4 Boxplots of the level of dispersion within the three substrate types which was calculated as the mean distance of each sample point to the centroid of the respective substrate type. The centre line in the box displays the median and the margins of the box specify the 25th and 75th percentile. Whiskers extend to the smallest (lower whisker) or the largest (upper whisker) value within the range of 1.5 × interquartile range that sites covered with macrophytes provided greater availability of food resources, microhabitats to inhabit and shelter from predators (Čerba et al. 2022). McGoff and Irvine (2009) associated positively littoral macroinvertebrate abundance with both macrophyte PVI (Percentage Volume Inhabited) and the extent of macrophytes lakewards in their study between Lake Habitat Quality Assessment and macroinvertebrate community structure.
Overall, the substrate type seems to be a primary factor affecting macroinvertebrate assemblages (along with other natural factors such as lake size and depth) (Timm and Möls 2012). We recorded this spatially variability of benthic fauna but our results were rather weak. Statistically speaking, the distinct macroinvertebrate assemblages were a result of differences in location of their centroids and not homogeneity of dispersions but still we could consider other factors affecting our dataset such as the dissimilarity measures (Bray-Curtis, Jaccard etc.) which could have altered the interpretation of the results (Anderson et al. 2006). As Barnes et al. (2013) pointed out, species richness is mainly affected by complexity and not heterogeneity but still they advocate to take into consideration not only habitat structure but the processes involved in these relationships. The effect of substrate types alone on species richness in lakes is less studied and the fact that we did not find any strong relationships could be attributed to the lack of more detailed categories (macrophyte cover, substrate size). Macrophytes seem to play an important role in regulating the balance of all trophic relationships occurring in the littoral zone of lakes. Tolonen et al. (2003) suggested that in an oligo-mesotrophic lake with a well-established macrophyte zone, the composition and size of zoobenthos differed significantly along the gradient of vegetation density horizontally from shore to open water. On the other hand, they suggested that shallow eutrophic lakes might be affected by other natural factors such as water level fluctuation and nutrient loading. Sediment deposition can significantly alter the composition of the bottom substrate, enriching it with organic matter. Jurca et al. (2021) pointed out that there was an absence of species with specific mesohabitat preferences in eutrophic lakes, as there is a clear linkage of macrophyte presence and species richness with eutrophication. All these factors should be taken into consideration when sampling for benthic macroinvertebrates in littoral zones of lakes in order to be able to draw conclusions without the confounding effects.
It is crucial to acknowledge that some macroinvertebrate taxa are depending on certain substrate types; this can be useful when developing and implementing site-adapted conservation and restoration measures (Brauns et al. 2007). In this context, priority is advised to be placed to the conservation and/or restoration of the natural littoral habitats that host rich macroinvertebrate assemblages. The need to conserve and were appropriate restore aquatic vegetation is obvious as our results suggested that sites dominated with macrophytes exhibit the greatest number of taxa and the most diverse group of Odonata. Controlling the extent and frequency of vegetation cutting and dredging could lead to habitat heterogeneity resulting on the dispersal of the macroinvertebrate community (Gething et al. 2020).
Acknowledgements The present study was conducted in the frame of the Greek National Water Monitoring Network, according to the JMD 140384/2011, implemented by The Goulandris Natural History Museum, Greek Biotope/Wetland Centre (EKBY). The network is supervised by the General Directorate for Waters of the Ministry of Environment and Energy. The data used in this research come from Act MIS 5001204 financed by the European Union Cohesion Fund (Partnership Agreement 2014-2020), and from Acts MIS 371010, 371138, 371140, 371145 financed by the European Regional Development Fund (National Strategic Reference Framework 2007-2013. K. Argiriou, M. Bozatzidou, G. Nakas, and V. Navrozidou contributed to sorting and taxa identification. H. Hadjicharalambous contributed to data analysis. EKBY's personnel conducted samplings and contributed to data analysis.

Data availability
The dataset generated during and analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.