Response of Chironomids (Diptera, Chironomidae) to Environmental Factors at Different Spatial Scales

1. Introduction

Chironomids are one of the aquatic insect families with the highest species numbers, including up to 15000 described species, often contributing to the highest observed diversity in the sample, including species tolerant to environmental extreme conditions. [1].

It is clear that species respond to a suite of different environmental (=abiotic) factors, such as water temperature, salinity, sediment composition [2], but biotic factors (competition, predation) were also considered as relevant [3].

Interactions between different factors vary in different situations confounding the interpretation of results. Spatial and temporal contiguities are surely responsible of species abundances making the interpretation of the real causes of species distribution questionable. The situation can be still more intriguing in unstable habitats as Mediterranean waters [4] where relevant fluctuations are observable in different years, actually still more exacerbated by global climatic change [5]. Historical (biogeographic) factors influence chironomid species composition only at very large spatial scale [6], but in special habitats as headwaters it can be supposed they act at smaller spatial scale, even if a clear evidence was given only considering other insects orders as Trichoptera [7]. Despite the complexity of situations, often few factors were emphasized to account for the largest proportion of variation [8].

The scarcity of samples available, different sampling tools used, different taxonomic expertize of researchers contributed to increase the uncertainty in interpreting results.

The response of chironomid species at different spatial scales [9,10] was rarely investigated in the past, and factors operating at small spatial scales were generally considered prevalent ([11,12].

The aim of the present review is to confirm the importance of environmental factors acting at low spatial scales and to analyze the possible importance of factors acting at larger spatial scales using the most recent tools used to analyze data based on Moran’s eigenvector maps (MEM) [13], which, starting from a geographical coordinates distance matrix, are able to extract an eigenvectors matrix explaining species distribution at different spatial scales.

2. Materials and Methods

The input data matrix included 788 samples collected in 11 different habitats (Table 1) from lotic and lentic waters, the former classified according to [14], the latter according to [15,16]. The samples were collected with different tools according to their habitat type. Kick sample was generally used in kryal, krenal and rhithral, drift samples were collected in potamal, a ponar dredge was used in AL03, AL04,AL05, AL06 and ME lakes, an Ekman dredge was used in Alpine lakes (ALalp). The data are from samples collected in Italy [2,8], Austria and Germany [17], Switzerland [18,19], Montenegro [20], Greece [21], Algeria [22,23], France (Garonna River, unpublished data).

The maps of Europe and Italy were prepared using the QGIS version [24] 3.34.1 Prizern (2009) (https://www.qgis.org/it/site/).

Raw data included 19531 samples and 255 variables. Only sites including values of at least 10 variables and variables present in at least 100 samples were retained. This selection left 782 sites with 101 species and 9 environmental variables: altitude, distance from the source, conductivity, dissolved oxygen, oxygen % saturation, pH, total phosphorous, ammonium hydroxide. Four spatial extensions (regional, country, water basin, sampling site), habitat type, year and month were also included as factors in database.

Data were analyzed with a partial correspondence analysis (pCCA) after log(x+1) transformation of species matrix and standardization of environmental variables. The pCCA was carried out with species matrix as dependent variables, the 9 environmental variables as constraining variables and spatial coordinates as conditioning variables (covariates). Forward and backward selection of environmental variables was performed to select the variables accounting for the largest source of variation.

To analyze the usefulness of habitat type in driving species composition a multiple discriminant analysis was performed with species as dependent variables and habitat type as discriminant factor.

To analyze the influence of spatial factors, the sites having the same spatial coordinates, but sampled in more than one date, were pooled into a single record, leaving a matrix with 264 sites, with 54 species present in at least 50 sites. The spatial coordinates were expressed in the coordinate reference system WGS 84 / UTM zone 32.

Univariate Moran’s I spatial correlograms of each environmental variable and of each species were calculated. Multivariate Mantel’s correlograms for the full set of environmental variables and species assemblage were also calculated [25].

A spatial distribution was then analyzed with the R-packages adespatial, spdep, vegan, to calculate the Moran’s Eigenvector Maps (MEM) using the function mem [26]. The correlation coefficients between the environmental variables and Moran’s I eigenvectors matrix were then calculated. The species matrix was then submitted to a Multiscale Pattern Analysis (MSPA) [27] using the function mspa, to calculate the spatial structure of the data and identifying the main scales of spatial variation.

A Spatial Weighting Matrix (SWM) was selected using the functions mem.select, listw.candidates, listw.select included in adespatial package [28]; two SWM for which the subset of MEMs yields the highest adjusted R2 with species table were chosen.

As a last step a canonical analysis (redundancy analysis using the function pcaiv of package ade4) was used to explain the structure of species table constrained by MEM spatial variables selected starting from the best spatial weighting matrix (SWM). See Figure 1 in [27] for a detailed explanation of all these steps.

Details of MEM, MSPA and SWM calculations are can be followed in [29] at the web site: https://cran.r-project.org/web/packages/adespatial/vignettes/tutorial.html (2023-10-18).

Diversity was calculated: 1- as the mean species number present at each spatial scale (alpha diversity), 2- the total number of species present (gamma diversity), 3- ratio between gamma and alpha diversity (beta diversity).

3. Results

The full matrix of environmental and species data with sites as rows and variables as columns is in Appendix S1.

The list of the species included in data analysis is in Table 2.

A partial canonical constrained ordination (pCCA) analysis was carried out using environmental variables as constraining variables and spatial factors as conditioning variables. Results of pCCA and its inverse, that is with spatial variables as constraining variables and environmental variables as conditioning variables are in Table 3. pCCA confirmed results of previous investigations [2], with a clear separation of high altitudes sites, with low water temperature from all the other sites (Figure 1). The separation was evident along the first axis and was explained by an altitudinal – temperature gradient; the second axis emphasized a clear trophic gradient, with sites rich in total phosphorus and with a low oxygen concentration separated from sites with low TP and high oxygen concentration.

Figure 1.

The plot of sites evidenced a clear habitat separation (Figure 2), with kryal sites separated on the left and potamal sites on the right bottom. The alpine lakes at high altitude were plotted on the right of the kryal sites. The separation of lowland lakes from rhithral sites was less evident, but rhithral sites showed a crowding of sites in the center of the graph, while sites sampled in lakes surrounded them. The plot of species (Figure 3) show the species preferring the different habitats, separating species characteristic of kryal and alpine lakes (on the left of the Figure) from species characterizing potamal and rhithral. The separation of species preferring krenal or littoral of lakes is less evident here.

Figure 2.

Figure 3.

An inverse partial cca carried out using spatial coordinated as constraining variables and environmental variables as conditioning variables emphasized that the second axis separated sites and species according to latitude (Y axis), and longitude (X axis) (Figure 4).

Figure 4.

Discriminant analysis using habitat types as factors still better emphasized the separation of habitats with a very good agreement between expected and observed classification (Figure 5, Table 4, Appendix S2). It can be seen that some habitat as rhithral have more sites predicted than assigned while others as lakes AL05-AL06 have a lower number of predicted sites than initially assigned.

Figure 5.

The univariate spatial correlograms showed that spatial autocorrelation decreases rapidly with increasing distance class both with environmental variables and species; the highest autocorrelation for environmental variables was observed with conductivity (Figure 6a), for species with T. gregarius (Figure 6b). The mantel correlograms showed that Mantel’s autocorrelation vanished just after the first distance class (Figure 7a,b).

Figure 6.

Figure 7.

The most significant Moran’s eigenvectors (MEM) generated from spatial weighting matrix (SWM) were calculated and plotted in the European map. The full matrix of MEM is given in Appendix S3. The eigenvector values of the most significant MEM were divided into 5 classes of values and marked with five different colors for their representation (Figure 8). The MEM of the 1^st axis was related to latitude, the MEM of the 2^nd and 3^rd to longitude, while the interpretation of the other axes was not clear. Few significant correlations between environmental variables and MEM eigenvectors were observed (Appendix S4). The R² correlations between the species matrix and MEM spatial predictors, calculated from the best SWM (see Data analysis), were filed in decreasing order according to R² values (Table 5). In Appendix S5 the R² correlations calculated from another SWM are given [27]. Redundancy analysis (RDA) was carried out to relate the sites x species matrix with spatial MEM eigenvectors; results are in Figure 9. A separation of sites according to habitat (lake type, running water type) is still apparent, but it does not seem that spatial factors aid substantial information. The detailed results of MSPA with the scores of species, sites and MEM variables are in Appendix S6.

Figure 8.

Figure 9.

The highest beta and gamma diversity were observed in rhithral and potamal samples, the lowest in volcanic lakes, while the highest alpha diversity was observed in lakes AL06 and the lowest in alpine lakes ALalps. The highest gamma diversity (93 species) was observed in Sarca river basin, the highest beta in Adda river basin, at a lower spatial scale the highest gamma diversity (77 species) was observed in Lambro stream, the highest beta in Bormida stream (Appendix S7).

4. Discussion

The present study focused on the influence of environmental variables and spatial factors to taxonomic composition of chironomid assemblage.

Taxonomic composition was considered too much cumbersome in ecological studies by many authors in recent years, who tried to substitute or accompany taxonomic composition with functional composition descriptors [30] or with species traits analysis [31,32]. The choice influenced the interpretation of results at different spatial scales. For example it was emphasized [5] that traits species composition could be less affected by regional scale than taxonomic composition, while ecoregion and season accounted for 20.5% of the variance of functional composition, but only 10.9% of the variance of taxonomic structure [9].

There are few or no contribution that can be compared with our study. Feld & Hering [9] investigated the benthic macroinvertebrates at different spatial scales(ecoregion, catchment, reach, site) in four different countries (Sweden, the Netherlands, Germany and Poland). Mykräet al. [10] focused their study on Fennoscandia bioregion and, following a pCCA analysis, concluded that local factors are prevalent in explaining macroinvertebrate distribution, but regional factors should in any case be considered at a larger spatial scale. In the present analysis many countries were included, i.e. Italy, France, Switzerland, Austria, Germany, Greece, Montenegro, and Algeria, even if samples from Italy were largely prevalent including West, Central, East Alps, Prealps, lowland Po River, different rivers in Central and South Italy (Appendix S1).

A comparison of our results with other studies having analysis of spatial factors as the main objective [9] are not easy because of 3 important differences:

[9] analyzed only running waters samples, while here both lotic and lentic samples were included,

[9] considered the whole macrofauna, while here only the chironomid family was analyzed,

[9] included in the analysis 31 variables at 4 different spatial scales (mega, macro, meso, micro), while in the present analysis only 9 environmental variables and 2 spatial variables were included,

[9] used both functional and taxonomic variables, while here only taxonomic variables (species) were included as dependent variables.

Feld & Hering [9] concluded that interactions at different spatial scales confounded the interpretation of results, but in any case reach scales (meso scale) were the variables accounting for the largest source of variation. Despite the different experimental design in the two studies, variables at small spatial scale confirmed their prevalent effect in both cases.

It is evident from the Moran’s maps that there is a spatial separation of species according to latitude and longitude, but it is surely an indirect effect, because the sites at highest altitudes are more frequent in the west part of the sampled area (West Alps), while the separation according to latitude is clearly bound to the fact the sites at the lower latitude have higher temperatures than the northern sites and effect of temperature on chironomid distribution is well documented [33,34]. The factors more responsible of chironomids species distribution are confirmed to be a joint combination of substrate, water temperature, conductivity, current velocity, dissolved oxygen[2,8] effects, which can be efficiently described using the terms kryal, krenal, rhithral, potamal, lentic or lotic habitat, clearing summarizing preferences to different factors interacting to each other. Discriminant analysis emphasized that there is a very good agreement between the ‘a priori’ and ‘a posteriori’ site classification, even if some ubiquitous species are present together specialists. The anthropogenic stress emphasized a clear response of chironomids to oxygen shortage [2,8,18], while it is difficult to establish the effect of toxic substances [35].

The present analysis shows that spatial factors generally have only indirect effect, mediated by the environmental factors present at a low spatial scale.

A very debated point is the effect of spatial scale on diversity of communities, with gamma diversity acting at larger spatial scale than alpha diversity [36].

The sampling design of the present database is not suitable to estimate diversity, because of the differences in sampling effort in different samples. [11,12] analyzed the influence of spatial factors on macroinvertebrates diversity considering the effect of river order on diversity of chironomid assemblages and concluded that streams of intermediate order, more than large and low order ones had the highest diversity, suggesting that factors at intermediate spatial scale are more relevant that factors at large and small spatial scale, supporting the old intermediate disturbance hypothesis [37]. The present database is not suitable to estimate diversity in streams of different order; we can only observe that there are no substantial differences in alpha, beta and gamma diversity at three different spatial scales.

It is obvious that the present conclusions are a consequence of the database available. For example, it is well known that at larger spatial scales, that is considering chironomid species present in different world areas (Palaearctic, Nearctic, Oriental, Neotropic, Australian) the species composition is quite different, but within the restricted area considered, which includes countries of the Mediterranean-Alpine area, the spatial factors should have no influence, if not mediated by other factors (altitude, water temperature etc.).