Turnover Drives High Benthic Macroinvertebrates’ Beta Diversity in a Tropical Karstic Lake District

Abstract

Beta diversity is useful to explain community assembly across landscapes with spatial variation. Its turnover and nestedness components help explain how beta diversity is structured across environmental and spatial gradients. Assessing beta diversity in freshwater ecosystems is essential to conservation, as it reveals the mechanisms that maintain regional diversity. Nonetheless, so far, no studies have examined the beta diversity patterns of benthic macroinvertebrates in tropical lakes. We aimed to examine the beta diversity patterns and components of the deep benthic macroinvertebrate communities of tropical Lakes of Montebello, Mexico, along spatial and environmental gradients. We used presence/absence data of deep benthic macroinvertebrates from 13 lakes distributed along environmental and spatial gradients. We calculated beta diversity indices and correlated them to each lake’s environmental and spatial variables. The macroinvertebrate communities of the Lakes of Montebello showed high beta diversity driven by a turnover pattern that emphasises the importance of regional-scale conservation efforts. Short distances between lakes and high environmental heterogeneity promoted species turnover, resulting in a great singularity level among lakes. We did not find significant correlations between the beta diversity components and the environmental variables, suggesting a random distribution given by the species’ high dispersal capacity in a reduced spatial extent across the lake district.

1. Introduction

Biodiversity is a major driver and stabiliser of ecosystem services [1]. Unfortunately, ecosystem services are at risk due to high rates of global biodiversity loss occurring in recent decades [2]. Conservation and restoration initiatives depend on knowing community–environmental linkages structuring biodiversity [3]. Biodiversity can be observed in three scales: alpha, beta and gamma [4]. Alpha diversity—diversity at the local scale—is related to gamma diversity—diversity at the regional scale—through species differentiation between sites known as beta diversity [4]. Beta diversity is helpful to explain community assembly across landscapes with spatial variation [5]. Two processes, turnover and nestedness [6], are comprised in beta diversity, which makes it helpful to explain the community assembly across landscapes with spatial variation [5]. The first process, turnover, is the replacement of some species by others between locations because certain species are unique to discrete places [6]. The second process, nestedness, mirrors the loss of species from discrete locations that happened when specific sites contained smaller subsets of species of other richer sites elsewhere [3]. The relative importance of beta diversity components in structuring communities is related to species dispersal abilities and local environmental characteristics [7,8]. The turnover component usually suggests environmental filtering [9] or historical constraints [10], whereas the nestedness component reflects dispersal and extinction-related phenomena that promote community disaggregation [11,12].

Spatial and local environmental predictors of beta diversity have a different range of variation. For example, some factors may vary at small spatial extents causing high species turnover in small areas, while others change only at larger scales, and patterns emerge at a large spatial extent [13,14]. Spatial predictors are mainly related to species dispersal abilities and are more complex in actively dispersing organisms like aquatic insects [7]. On the other hand, local environmental conditions are associated with habitat type and niche determination. They ultimately sort the species set that constitutes a community at a given location [15]. Some studies have found that species similarity decays with environmental or spatial distances in freshwater ecosystems [5,16,17]. However, the role of habitat type, spatial scale and performance of different types of dispersers on aquatic species distribution and beta diversity patterns is context-dependent and difficult to predict [18].

Freshwater communities can be structured by dispersal limitation, biotic interactions, environmental sorting or stochastic events [19,20,21,22]. Different studies have shown: (a) a prevalence of environmental sorting in structuring communities between water bodies [5,16,22], (b) the spatial distances and dispersal limitations strongly related to community structure [23], (c) the effect of top–down controls over the community structure [21], (d) the joint importance of spatial processes and local environmental factors as determinants of community structure [24,25,26,27] and (e) a mismatch between community composition and environmental conditions due to stochastic events [22]. The distribution patterns of benthic macroinvertebrates are affected by local factors, including lake size, depth, pH, nutrients, biotic interactions [24] and by dispersal factors, speciation and extinction processes [28]. The relative importance of the factors on turnover and the nestedness structuring process in benthic macroinvertebrates communities requires further investigation.

The “Lagunas de Montebello” National Park (LMNP) is a tropical karstic lake district with more than 50 lakes displaying an ample range of morphometric and physicochemical characteristics, e.g., from shallow to deep, from small to large and from oligotrophic to eutrophic lakes. The heterogeneity of environmental characteristics of lakes with similar origins within a reduced area turns the Lakes of Montebello into a natural experimental site for testing community distribution hypotheses. The natural and anthropogenic variation in habitat conditions among the lakes likely determines which species can occur in each lake or the species sorting perspective of the metacommunity theory [15]. Previous studies on the LMNP evaluated the zooplankton [29,30], the littoral benthic macroinvertebrates [31], the aquatic springtails [32] and the deep benthic macroinvertebrates composition and bathymetric distribution [33,34]. Deep benthic macroinvertebrates composition showed high regional diversity and richness per lake—alpha diversity—similar to other tropical lakes [34]. However, a study is lacking that evaluates the beta diversity and distribution patterns of the deep benthic macroinvertebrates’ communities across the Lakes of Montebello.

Assessing patterns and components of beta diversity may provide important information about the assembly mechanisms and distribution of aquatic communities [35]. The global diversity of freshwater invertebrates is declining [2,36]. Knowledge about their structuring patterns is essential to develop preservation and restoration programmes, as it reveals the mechanisms that maintain beta and gamma diversity [18]. However, few studies have examined the beta diversity patterns of the benthic macroinvertebrates communities in lakes, e.g., [5,8,37,38], most of which were conducted in the littoral zone of temperate and boreal lakes. To our knowledge, this is the first study examining the beta diversity patterns of deep benthic macroinvertebrates along an environmental and spatial gradient of tropical lakes.

This study aimed to examine the beta diversity patterns of the deep benthic macroinvertebrates’ communities of the Lakes of Montebello along spatial and environmental gradients. For this purpose, we addressed the following research questions: (a) How are the deep benthic macroinvertebrates distributed at the Lakes of Montebello related to the environmental and spatial gradients? and (b) are the beta diversity components of the deep benthic macroinvertebrates correlated to the environmental and spatial variables? We proposed the following hypotheses: (1) the deep benthic macroinvertebrates dissimilarity increases with distance and environmental differences between lakes, (2) the turnover component of beta diversity mainly correlates to the water and sediment physicochemical variables and (3) the nestedness component mainly correlates to the distance between lakes and the morphometric variables. To test these hypotheses, we used presence/absence data of deep benthic macroinvertebrates from 13 lakes of the LMNP, distributed along environmental and spatial gradients, and we calculated beta diversity indices and correlated them to the environmental and spatial variables measured at each lake.

2. Materials and Methods

2.1. Study Area

The LMNP (Figure 1) is a natural protected area (1959) and Ramsar site (2003) in south-southeast Chiapas, encompassing the municipalities of La Trinitaria and Independencia. The LMNP borders with Guatemala—16°04′40″ N to 16°10′20″ N, 91°37′40″ N and 91°47′40″ N, 1200–1800 m a.s.l. [39,40]. It is a karst basin made up of sedimentary sequences of marine and transitional origin, which have undergone karstic and tectonic processes leading to the formation of numerous dolines, uvalas and poljes [40], with varied morphometric characteristics [41].

The climate in the region is type C(fm), that is, humid temperate with rains throughout the year; however, the NW end of the LMNP has an A(cm)-type climate, that is, hot humid with abundant rains during the summer. The average annual temperature is 23.6 °C, reaching a minimum average of 20.9 °C in January and a maximum average of 25.6 °C. The minimum rainfall occurs during May (40 mm), while during the rainiest month (September), rainfall reaches up to 1400 mm [39,42].

Soil types identified in the LMNP include Acrisols, Fluvisols, Greysols, Vertisols, Lithosols and Rendzinas [43]. The native vegetation of the LMNP includes pine–oak–liquidambar forests, cloud forests and riparian vegetation [44,45]. However, there has been an extensive land-use change for agricultural areas (e.g., tomatoes, corn, beans and coffee) and urban development [46].

The LMNP belongs to the Río Grande de Comitán-Lagos de Montebello basin. This basin is part of Hydrological Region No. 30 Grijalva-Usumacinta. The Río Grande de Comitán basin is made up of the main tributary Río Grande. The feeding of the lake area is mainly underground, generating a system of sinkholes and other forms caused by the dissolution of limestone rocks [40].

Durán et al. [40] classified the LMNP lakes into two groups: (a) plateau lakes, on the NW, fed by surface runoff from the Río Grande de Comitán and groundwater and (b) mountain lakes, on the SE, at higher elevations in relation to the other lakes and fed mainly by groundwater through regional systems of faults and fractures.

2.2. Field Sampling

A set of 13 lakes was selected that represented the morphometric variation and trophic status of the systems that make up the LMNP lake district. These are the eutrophic Balantetic, San Lorenzo, Bosque Azul and La Encantada, and the oligotrophic Esmeralda, Ensueño, Agua Tinta, Cinco Lagos, Pojoj, Kichail, Tziscao, Patianu and Dos Lagos (Figure 1).

Five sampling dates were carried out between 2013 and 2016, covering both characteristic seasons of the tropics (i.e., cold-dry season, when the lakes were mixing, and warm-rainy season, when the lakes were stratified). Logistics prevented applying the same sample effort at each lake. Sampling took place at the deepest portion of each lake. Further details about sampling are described by Cortés-Guzmán et al. [34].

We used a water quality multiparametric sonde Hydrolab DS5 (OTT Hydromet, Loveland, CO, USA) to measure bottom water (1 m above sediments) in situ physicochemical parameters (temperature, dissolved oxygen, electrical conductivity and pH). A 5-L Uwitec GmbH (Mondsee, Austria) water sampler bottle obtained water samples for chlorophyll-a (Chl-a) concentration. Samples for Chl-a concentration were analysed in a Turner Design 10-AU fluorometer following the EPA method 445.0 [47]. Sediment samples were obtained with an Ekman-type dredge (0.0225 m²) to analyse texture (Beckman Coulter LS230 laser diffraction analyser), organic matter content (lost on ignition, LOI, at 550 °C) and carbonates (acidulation). The same dredge was used to obtain the macroinvertebrate samples, with three replicates per station. Samples were sieved in situ through 250 μm mesh size. Collembolans and amphipods were determined with the help of specialists, and the other groups were determined following specialised taxonomic keys, e.g., [48,49,50,51].

2.3. Statistical Analysis

We calculated three beta diversity indices based on the Sørensen pairwise dissimilarity developed by Baselga [6,52] using the presence/absence data pooled of all sampling dates. Sørensen dissimilarity (β_SOR) is a measure of total pairwise dissimilarity. Simpson dissimilarity (β_SIM) measures the turnover component of Sørensen dissimilarity. Nestedness-resultant dissimilarity (β_SNE) measures the dissimilarity derived from the nestedness component of the Sørensen dissimilarity. Indices were calculated using the function index.family = “sor” of the betapart package [6,52] in R [53]. To assess each species’ contribution to dissimilarity, we applied a SIMPER test in PAST. SIMPER runs pooling all lakes using Euclidean distances as a measure of dissimilarity. We performed a non-metric multidimensional scaling (non-metric MDS) analysis using the function “metaMDS” of the vegan package in R to know the species distribution among the Lakes of Montebello.

We assessed the correlation between the beta diversity components, the water and sediment physicochemical variables and the lakes’ morphometric variables through bio-env [54] and Mantel [55] tests. Water variables included temperature, dissolved oxygen concentration, pH, electrical conductivity and Chl-a concentration. Sediment variables included sand, silt, clay, carbonates and organic matter percentages. Morphometric variables included surface area, maximum depth, mean depth, maximum length, maximum width and the distance between each lake and Lake Tziscao at the extreme SE. We ran this analysis for mixing and stratification seasons separately.

First, we obtained three dissimilarity matrices (Sørensen dissimilarity, turnover and nestedness) based on the species presence/absence data using the function “beta-pair” of the betapart package [6,52] in R. Then, we obtained the best environmental distance matrices correlating each dissimilarity matrix with the physicochemical variable matrix and the morphometric variable matrix through the bio-env test. This method is based on standardised environmental variables and tests all combinations to show the strongest correlation between biological dissimilarity and environmental distance matrix. Finally, we used the Mantel test to examine the statistical significance between beta diversity components and the environmental distance matrices selected by the bio-env test. Bio-env and Mantel tests were based on Pearson correlation and 999 permutations for obtaining p-values. Bio-env tests were run with the function “bioenv” and the Mantel test was run with the function “mantel” of the Vegan package in R.

3. Results

3.1. Environmental and Morphometric Variables

Lakes of Montebello displayed an ample range of morphometric characteristics (

Table 1. Morphometric parameters of the Lakes of Montebello (Alcocer et al. [41]): L_MAX, maximum length; b_MAX, maximum width; A, surface area; Z_MAX, maximum depth; Z_MEAN, mean depth; d, straight-line distance between each lake and Lake Tziscao at the extreme SE; lakes from NW to SE.

Lake		LMAX	bMAX	A	ZMAX	ZMEAN	d
		(m)	(m)	(m2)	(m)	(m)	(km)
Balantetic	BL	814	232	135,636	3	2	12.5
San Lorenzo	SL	3091	1290	1,813,415	67	12	10.9
Bosque Azul	BA	1322	822	524,558	58	20	7.4
La Encantada	EC	385	306	82,004	89	29	6.6
Esmeralda	ES	140	113	11,241	7	4	6.4
Ensueño	EN	220	191	27,048	35	22	6.2
Agua Tinta	AT	210	196	29,661	24	15	6.2
Cinco Lagos	CL	817	595	236,574	162	43	3.5
Pojoj	PO	1057	736	436,879	198	35	3
Kichail	KI	582	445	124,865	22	10	2.7
Patianu	PA	262	184	34,004	26	11	1.9
Dos Lagos	DL	336	234	52,454	42	25	4.5
Tziscao	TZ	3201	1481	3,065,515	86	29	0

). Mean depth ranged from 2 to 43 m, while maximum depth from 3 to 198 m. Seven out of 13 lakes were deep, showing a mean depth of 20 m or higher. Lakes surface area ranged from 11,241 m² (Esmeralda) to 3,065,515 m² (Tziscao). Straight-line distances between lakes were up to 12.5 km (from Lake Balantetic at the extreme NW to Lake Tziscao at the extreme SE).

Water and sediment parameters also varied among lakes (

Table 2. Ranges of water physicochemical and sediment parameters of the Lakes of Montebello lakes: DO, dissolved oxygen concentration; K₂₅, electric conductivity; T, temperature; Chl-a, chlorophyll-a concentration; CO₃²⁻, carbonates; OM, organic matter. Lakes’ abbreviations are as in Table 1.

Lake	Water Physicochemical Parameters
	T (°C)	DO (mg L−1)	pH	K25 (µS cm−1)	Chl-a (µg L−1)
BL	20.7–22.3	4.4–18.1	7.6–8.8	474–945	42.1–148.4
SL	17.3–19.5	0.0–0.1	7.2–7.9	682–756	23.0–49.8
BA	17.3–17.6	0.0–0.0	7.2–7.3	540–578	11.1–29.1
EC	17.9–18.1	0.0–0.4	7.1–7.8	431–475	8.2–38.9
ES	19.0–23.4	5.5–6.4	7.7–8.2	341–374	0.8–1.1
EN	18.9–19.9	0.1–7.3	7.4–8.2	258–276	0.4–0.5
AT	20.9–21.4	1.9–2.2	7.3–7.4	321–345	0.1–0.8
CL	18.3–18.4	1.5–6.1	7.2–8.3	212–221	0.2–0.3
PO	18.3	0.1	7.1	275	0.9
KI	17.4	6.7	8.2	284	0.8
PA	17.3–18.0	0.0–3.8	7.1–8.0	277–284	0.5–0.9
DL	18.5–18.6	0.0–0.1	7.1–7.2	1424–1437	0.8–1.2
TZ	18.0–18.1	0.9–7.2	7.3–8.2	245–280	0.4–0.5
Lake	Sediment Parameters
	Sand (%)	Silt (%)	Clay (%)	CO32− (%)	OM (%)
BL	0	66	34	34	21
SL	0	49	51	34	23
BA	12–21	70–79	8–9	48–68	22–26
EC	20	74	6	66	29
ES	5–37	49–80	13–14	9–61	16–51
EN	13–21	57–69	18–21	22–26	20–21
AT	26–32	62–70	4–6	47–65	21–30
CL	68–90	9–31	0.5–1	74–80	42–86
PO	71	28	1	63	20
KI	57	43	1	10	54
PA	0–2	64–70	28–36	11–17	19–23
DL	41–52	46–56	2–3	40–45	35–41
TZ	5–8	77–80	12–18	40–81	21–23

). Plateau lakes (8.2–148.4 µg L⁻¹) at NW showed a higher chlorophyll-a concentration than mountain lakes (0.1–1.2 µg L⁻¹). All lakes were warm (>17 °C), neutral to slightly basic (7.1–8.8) and with moderate conductivity (212–945 µS cm⁻¹), except for Dos Lagos, which showed a high conductivity (1424–1437 µS cm⁻¹). Sediment parameters varied randomly among lakes (Table 2). Lakes showed a high concentration of carbonates (>10%) and organic matter (>15%).

3.2. Beta Diversity

The deep benthic macroinvertebrate composition of the Lakes of Montebello was previously described by Cortés et al. [34]. The Sørensen index showed a very high beta diversity (β_SOR = 0.918). The turnover (β_SIM = 0.770) was the main factor contributing to beta diversity, while nestedness contributed a minor proportion (β_SNE = 0.148). The non-metric MDS, including species collected at stratification and mixing seasons, showed a random distribution according to the distance among lakes (Figure 2).

SIMPER tests showed an average dissimilarity of 15.72% in the stratification season among all lakes. The 10 species that contributed the most to dissimilarity summed up 25.44%, and most of them were found in Tziscao and Ensueño (

Table 3. Species that contributed the most to dissimilarity among the Lakes of Montebello during the mixing and stratification seasons.

Taxa	Season	Contribution to Dissimilarity (%)	Cumulative Contribution (%)	Lakes
Procladius sp.	Stratification	3.18	3.18	AT, EN, TZ
Amphipoda genus 1		2.47	5.65	EN, TZ
Ballistura libra		2.47	8.13	BA, TZ
Harnischia sp.		2.47	10.60	DL, ES
Homochaeta sp.		2.47	13.07	BA, TZ
Chironomidae genus 1		2.47	15.55	DL, EN
Polypedilum sp.		2.47	18.02	DL, EN
Probezzia sp.		2.47	20.49	DL, EN
Chaoborus sp.		2.47	22.97	PO, TZ
Lepidocyrtus lanuginosus		2.47	25.44	BA, TZ
Sminthurides sp.	Mixing	3.40	3.40	BL, CL, DL, EN, KI, PA
Microchironomus sp.		3.31	6.71	BL, DL, EN, PA, TZ
Homochaeta sp.		3.02	9.74	CL, DL, EC, TZ
Rheotanytarsus sp.		3.02	12.76	DL, EC, EN, TZ
Chaoborus sp.		2.55	15.31	BL, EN, TZ
Procladius sp.		2.55	17.86	DL, EN, TZ
Limnodrilus sp.		2.55	20.42	BL, DL, TZ
Naididae genus 1		2.55	22.97	DL, EN, TZ
Koenikea sp.		2.55	25.52	EN, ES, KI,
Polypedilum sp.		2.55	28.07	DL, ES, SL

). The average dissimilarity in the mixing season was 16.03%. The 10 species that contributed the most to dissimilarity summed up 28.07%, and most of them were found in Dos Lagos, Tziscao, and Ensueño (Table 3). Procladius sp. (Chironomidae, Diptera), Polypedilum sp. (Chironomidae, Diptera), Chaoborus sp. (Chaoboridae, Diptera) and Homochaeta sp. (Naididae, Oligochaeta) were among the 10 species that contributed most to dissimilarity in both seasons.

The bio-env test showed consistent results for all three beta diversity components at each season (

Table 4. Bio-env (correlation values are in brackets) and Mantel test (p-values are in brackets) results for the three beta diversity components (β_SOR Sorensen, β_SIM turnover, β_SNE nestedness), water and sediment physicochemical variables and lakes’ morphometric variables in the stratification and mixing seasons.

Season	Beta Diversity Index	Physicochemical Variables	Mantel Test	Morphometric Variables	Mantel Test
Stratification	βSOR	Chlorophyll-a, organic matter content (0.1727)	−0.26 (0.90)	Distance (0.4149)	−0.27 (0.94)
	βSIM	Chlorophyll-a, organic matter content (0.1407)	0.10 (0.46)	Distance (0.4065)	−0.19 (0.85)
	βSNE	Chlorophyll-a, organic matter content (0.1888)	0.05 (0.39)	Distance (0.4027)	0.06 (0.36)
Mixing	βSOR	Temperature, clay percentage (0.4033)	0.24 (0.07)	Distance, area (0.3861)	−0.09 (0.66)
	βSIM	Temperature, clay percentage (0.4240)	−0.34 (0.91)	Distance, area (0.3774)	0.09 (0.36)
	βSNE	Temperature, clay percentage (0.3942)	−0.03 (0.57)	Distance, area (0.3851)	0.07 (0.30)

). The beta diversity components were correlated to chlorophyll-a and organic matter content and the distance between lakes in the stratification season. Chlorophyll-a is a proxy of the lakes’ trophic status, related to the distance between lakes because there is a eutrophication gradient from NW to SE. In the mixing season, beta diversity components were correlated with the water temperature, clay percentage, distance between lakes and lakes’ surface area. However, the Mantel test showed that the physicochemical and morphometric variables did not significantly correlate to the beta diversity components (Table 4).

4. Discussion

Beta diversity of the Lakes of Montebello benthic macroinvertebrates is driven by high turnover in which species replace each other among lakes. The dominance of the turnover component in the beta diversity pattern suggests that the species filtering among lakes is driven by deterministic niche-related processes [9]. Niche determination for deep benthic macroinvertebrates is highly related to the lakes’ thermal regime (i.e., mixing and stratification) and trophic status (oligotrophic and eutrophic). Most Lakes of Montebello are deep and then thermally stratified during the warm–rainy season, when the deep bottom becomes anoxic, representing harsh conditions for the macroinvertebrates’ survival. In the cold–dry season, deep oligotrophic lakes mix, and the deep bottom reoxygenates favouring macroinvertebrates’ growth and survival [34]. Eutrophic lakes remain anoxic in the deep bottom even during the cold–dry season because of the high amount of decomposing organic matter. Seasonal changes in habitat conditions promote species turnover rather than nestedness diversity patterns because it produces mass migration or species emergence [3]. Similarly, in the coastal wetlands of The Great Laurentian Lakes (US), the benthic macroinvertebrates communities are turnover-structured because the habitats are inhospitable to most species over winter [3]. Meanwhile, spring promotes migration processes [3].

High beta diversity and a turnover-dominated pattern in the Lakes of Montebello suggest that the deep benthic macroinvertebrate community at each lake is unique and different from the other lakes. The macroinvertebrate diversity of the Lakes of Montebello is spread across the region rather than isolated in hotspots or single lakes, stressing the importance of preserving the entire landscape. Moreover, heterogenous benthic habitats yielded by the lakes’ environmental and morphometric features characterise the Lakes of Montebello. A larger environmental heterogeneity is usually a reflection of increasing spatial extent [56,57]; however, the characteristic combination of water and sediment physicochemical and morphometric variables of the Lakes of Montebello results in a great environmental heterogeneity and high beta diversity in a reduced (71.8 km) spatial extent. Regions with turnover-structured communities require conservation, restoration and protection of as many habitat locations as possible to maintain gamma diversity [3]. Habitat alteration and loss could represent species loss at the local—lake—and regional scales.

Species distribution in the Lakes of Montebello did not show a distance structuring pattern. Similarly, Soininen et al. [58] did not find significant spatial structuring in aquatic communities across 20 boreal lakes within a Finnish drainage basin. However, they found significant spatial structuring across five drainage basins at a total spatial extent of about 700 km. Within drainage basins, communities dominated by strong active, such as dipterans, or passive overland dispersers, such as oligochaetes in the Lakes of Montebello, are not spatially structured because organisms have no significant dispersal limitations [18,59]. Maximum straight-line distances between the Lakes of Montebello are 11 km and 71.8 km from the extreme NW to SE, while dispersal distances in freshwater macroinvertebrates are up to 1100 km [60,61,62,63]; therefore, distance is not a limitation for dispersal, explaining a random distribution or mass effect of the macroinvertebrates among the Lakes of Montebello.

Most species contributing to dissimilarity at Lakes of Montebello were dipterans—chironomids and chaoborids—with high dispersal capacity, as their adult stage is a flying insect. Increased species dispersal ability has been suggested to promote species turnover structuring [52,64,65,66]. Chironomids are dominant in the benthic macroinvertebrate communities of Lakes of Montebello [34]. The ability of generalist and mobile families, such as chironomids, to disperse and occupy many habitats results in turnover-structured communities due to the high replacement of species between locations [67]. Active overland dispersers (e.g., flying insects) tend to distribute over the landscape [68], particularly in relatively small drainage basins, such as the Lakes of Montebello sub-basin (≈580 km²). Consequently, strong dispersers usually show no significant spatial structuring [8,68].

The beta diversity components did not significantly correlate to the morphometric and physicochemical variables of the Lakes of Montebello. Although trophic status in the Lakes of Montebello shows a clear gradient from eutrophic lakes at NW to oligotrophic lakes at SE, most environmental features did not show clear gradients (e.g., substratum particle size, organic matter content) but chlorophyll-a concentration. If species are distributed across lakes through environmental filtering, high environmental variability is likely to produce high beta diversity [69] but randomly distributed if gradients are lacking.

The lack of correlation between diversity patterns and environmental variables in the macroinvertebrate communities of the Lakes of Montebello might be associated with the complexity of the distribution patterns and the missing of some important (e.g., biological) explanatory variables. Several studies have found low predictability of benthic macroinvertebrates communities based on environmental variables [22,27,68,70,71,72,73]. The low predictability of environmental variables may stem from two possibilities. First, the environmental variables measured may not account for variation in benthic macroinvertebrates composition [22]. We included water and sediment variables typically measured by researchers to explain the benthic macroinvertebrate distribution (e.g., the sediment grain size, organic matter content, dissolved oxygen concentration, temperature, trophic status, [74,75,76]). However, dispersal processes and stochastic events also affect species distributions across lakes that decrease the match between community composition and local environmental conditions [22,70]. Thus, sampling of biological and environmental variables fails to reveal strong distribution relationships [68,77,78]. Secondly, we did not include biotic interactions that might be important for structuring aquatic communities, mainly top–down control processes [70]. For instance, fish predation has been shown to modify the structure of macroinvertebrate communities by reducing or removing specific taxa in lakes [79,80]. The importance of trophic interactions might overcome environmental filtering, as Declerck et al. [21] demonstrated with a phytoplankton community strongly structured by the food web structure, independently of an environmental gradient.

To summarise, the deep benthic macroinvertebrate communities in the Lakes of Montebello showed high beta diversity driven by a turnover pattern that emphasises the importance of regional-scale conservation efforts. Short distances between lakes and high environmental heterogeneity promoted species turnover, resulting in a great singularity level among lakes. We did not find significant correlations between the beta diversity components and environmental variables; therefore, we cannot exclude the possibility of randomly distributed communities dominated by species’ high dispersal capacities and stochastic events.