Response of Chironomids to Key Environmental Factors: Perspective for Biomonitoring

Simple Summary Benthic macroinvertebrates of inland waters, including running waters and lakes, are frequently used in biomonitoring. Sometimes, environmental data associated with species lists are not available; in this situation traits or functional adaptations of species to environment can be considered as a tool to translate the list of species into a useful index to evaluate the environmental quality a body of water. Abstract Chironomids are the species-richest family among macroinvertebrates and are often used as indicators of ecological conditions in inland waters. High taxonomic expertise is needed for identification and new species are still being described even in the well-known West Palearctic region. Our Microsoft Access relational database comprises data on Chironomid species collected in rivers and lakes in Italy and some other European countries over a period of about 50 years, often associated with physical-chemical data, but in some cases, only data on Chironomids are available with no associated environmental data. The aim of the present paper was to propose the calculation of ecological traits of Chironomid species as a tool to derive information on water quality, when only data on Chironomid species composition are present, while environmental data are lacking. Traits summarizing the species’ response to environmental variables were evaluated, with emphasis on natural and man-influenced factors: current velocity, water temperature, conductivity, dissolved oxygen, and nutrients. Traits calculations were carried out in the R environment using a subset of our data, including both environmental data and Chironomid abundances. The relations between sites, Chironomid, species and traits were evaluated using correspondence analysis and other multivariate methods. The response of species showed an interaction among different factors, with the possibility of ordering species along a single environmental gradient, extending from cold running waters to warm standing waters, with few exceptions.


Introduction
The analysis of environmental factors responsible for macroinvertebrate assemblage structure has a long history. Among macroinvertebrates, Chironomids are considered a hard-to-identify group, therefore studies concerning macroinvertebrates as bioindicators have been often limited to Ephemeroptera, Plecoptera, and Trichoptera (EPT) [1]. Taxonomic problems were brought forward to justify this choice, but even though Chironomids are a difficult group, this is not a valid reason to disregard them. Chironomids include species living in almost all water bodies, sometimes present with a very large number of species, so their exclusion can lead to a serious misjudgment in water quality assessment [2].
A frequently overlooked problem in the development of a biotic index is species identification accuracy. Particularly in taxonomic hard groups, frequent mistakes in species identification were observed. It must be pointed out that different species within the same genus may show a different indicator value, therefore, an index based only on genus identification can lead to misleading conclusions with respect to an index based on the identification of species [3]. However, an intermediate level between genus and species, which we termed "morphotaxon", can be used to describe the ecological responses of Chironomidae [4].
On the other hand, the attempt to use Chironomid species as indicators of toxic substances did not make much progress, with the same tolerant/intolerant species probably being tolerant/intolerant to a set of many other different factors [17]. In contrast, studies concerning the response of Chironomid species to habitat heterogeneity or alteration were more fruitful [5][6][7].
The use of species identification in the assessment of water quality was criticized and then refined considering biological and ecological traits [18][19][20][21][22][23], suggesting that nontaxonomic aggregation of taxa as similar as possible in their species traits could aid in the interpretation of information given by a taxonomic list of species. For example, biological traits were preferable to taxonomic species lists in analyzing the response of multiple stressors in central European lowland rivers [24].
The problem is that the possibility of translating a list of species into biological and ecological traits needs basic research to prepare such a translation.
The aim of the present paper is to discuss the advantages and limitations of the use of ecological traits with respect to the taxonomic approach, testing a large database, including both lotic and lentic habitats, with multivariate data analysis. The discussion considers the situations where environmental data associated with species lists are lacking or scanty, and trait calculation is proposed as a method of overcoming the problem of missing environmental data.

Chironomid Database
A large database including data on Chironomid species collected in rivers and lakes in Italy and in some other European countries over a period of about 50 years was considered. Physical-chemical data associated with Chironomid samples were available, but only for a subset of samplings. Data were filed in a relational database in Microsoft Access© which is the property of the corresponding author. An extract of these databases is available in Table  S1. Data were stored in different Tables; the description of these Tables is here summarized.  Table S1. 1 Species: it contains a list of the variables used, including both environmental variables (morphometric, physical, chemical) and species belonging to the family Chironomidae; species were aggregated in species groups (morphotypes), each morphotype corresponding to a genus, a subgenus, a species group or single species [4]. As regard species, the table contains the species name, author, year of the original description, and taxonomic status (senior synonym, junior synonym, new combination) as additional fields.
Table S1. 2 Sites: it contains a list of the sampling localities, and other additional fields such as latitude, longitude, altitude, distance from the source (for running waters), depth (for lakes), habitat (krenal, kryal, rhithral, potamal, littoral, sublittoral, profundal, etc.). Table S1. 3 Conn: it connects each environmental variable or species with the sampling station and a numerical value; for environmental variables, this number is the measured value, for species it is an index of abundance (see below); additional fields are sampling year, month, day, sampling tool, and bibliographic source of information (when appropriate).
The samples here selected for data analysis included larvae collected with different tools, such as Surber net, kick net, hand net, grab, etc., and measures of environmental variables (e.g., water temperature, conductivity, nutrients, etc.) associated with Chironomid samples, when available. The species abundance value was the total number of specimens in the sample, after sorting under a stereomicroscope LEICA MZ12.5 magnification. However, since data on Chironomid larvae were mostly derived from quantitative or semiquantitative sampling techniques [25,26] for calculations species abundance was expressed as a total number of individuals per m 2 .
A crosstab query was then created with sites and other variables describing the sampling site as rows, and environmental variables or species as columns.
The crosstab query created produced a matrix with 9127 sampling sites, including lentic and lotic waters, sampled in different years and months, in Italy above all [4,27], but including also data from Algeria [28] and other countries in Europe [6]. A complete list of the sites examined is reported in Table S1. The same query included 160 columns, that is a row label, a sequence number, 6 factors, 11 environmental variables, and 143 Chironomid taxa. The 11 environmental variables included were: sampling year, sampling month, altitude in m a.s.l., distance from the source in km, O 2 content in mg L −1 , conductivity in µS cm −1 , pH, total phosphorous in µg l −1 , N-NO 3 in mg L −1 , N-NH 4 in µg l −1 , water temperature in • C. The 4 factors were: habitat, river basin, water body, and sampling station (Table S1).
The sampling year was included to detect potential temporal trends of species during the long time period considered (about 50 years). Month was used to indicate seasons.
The taxa included in the analysis were the morphotypes or species groups described in [3]; in the following section of the present work these taxa will be named "species" for simplicity, even if they are often taxa larger than species (i.e., genus or group of species).

Data Analysis
The crosstab query generated a matrix with n sites as rows and p species + s environmental variables as columns ( n M p+s ), which was input into an R script (Table S1). The M matrix was separated into an n L p matrix of species and an n R s matrix of environmental variables. Each environmental variable was used to calculate: 1-a correlation matrix between each species and the environmental variables p C s ; 2-a weighted mean of each environmental variable for each species, i.e., means of each environmental variable weighted according to species abundances, which can be considered the optimum for each species; 3a weighted standard deviation, which can be considered a measure of species tolerance. The weighted mean of each environmental variable for each species generated a trait matrix p U s with p species as rows and s environmental variables as columns [4,31]. The presence of missing data in the n R s matrix forced us to calculate matrices p C s and p U s matrices using only the available data. The n L p matrix, including the reduced n (=2258) sites and p (=91) species, and the p U s matrix, including the same species and s (=11) traits, were analyzed with a correspondence analysis (CA = unconstrained ordination) [31,32]. The n L p values were log(x + 1) transformed before calculation. As a second step, a canonical constrained ordination was carried out using the transpose of n L p , that is p L' n , and p U s as input matrices. As the last step, the n L p matrix was post multiplied by the p U s matrix, submitted an unconstrained ordination and compared with the previous results. The large amount of missing data in the n R s matrix hindered the ability to carry out a canonical constrained ordination between the n L p and n R s matrices.
The sites × species matrix n L p was post-multiplied by species × traits p U s matrix to obtain a site × traits matrix ( n L p U s ), i.e., a matrix with sites as rows and species traits as columns. This n L p U s matrix was also submitted to correspondence analysis.
A discriminant analysis was carried out to test the goodness of classification in different habitats using the Chironomid taxa assemblages: both the n L p and the n L p U s matrices were submitted to multiple discriminant analysis, using the habitats as a grouping factor. Finally, a cluster analysis of the site × species matrix was carried out using the complete linkage clustering method [32] to detect clusters of species.
Software: All data analyses were carried out in the R environment [33], using the packages vegan, Hmisc, mass, and scatterplot3D [34], including R scripts produced by the corresponding author, which are available on request.

Results
Measures of the 11 environmental variables were available for a reduced number of sites (Table 1), so the correlations, weighted means, and standard deviations of each environmental variable with each species were calculated using sites where both species and environmental data values were available (see Methods); when less than 4 records were available for the couple environmental variable-species, correlations were not calculated and mean values and standard deviations of the environmental variable calculated over all the other species were assigned to these species. Highly significant correlations (p < 0.01) between species abundance and environmental variables were observed for a reduced number of species (Tables 2 and S2).
The weighted means, considered the optimum values for each species [29,31], were used to create a species × traits matrix p U s with p species as rows and s environmental variables as columns (Table 3). Weighted standard deviation and the number of observations available are reported in Table S3.
The sites × species matrix n L p was submitted to a correspondence analysis; three major gradients (Figures 1, 2, S1 and S2, and Table S4) were evidenced, the former accounting for 7.6% of the total variance, the second 5.1%, the third 3.8% of the total variance, with eigenvalues equal to 0.71, 0.48, and 0.35, respectively. The species and sites ordered in the plane of the two axes showed the typical horseshow or arch effect [31,32]. The first gradient separated running waters from standing waters, the second separated upstream stations from downstream stations in running waters, with the following sequence ( Figures 1, 2, S1 and S2): 1-frigo-stenothermal species living in kryal were plotted in the bottom left of the graph; 2-rhithral species living in streams were plotted above the former; 3-eurithermal species, living in potamal, were plotted at the apex of the arch, extending from the top to the right part of the plot; 4-species living preferably in lentic waters, were plotted on the right part of the graph; 5-species living in springs were plotted in the central part of the area. Further separations of species from small alpine lakes, such as Paratanytarsus austriacus, Heterotrissocladius, Corynoneura, and Zavrelimyia, were plotted in the center of the area, species characterizing profundal zone of lowland large lakes, such as Micropsectra radialis and Paracladopelma, were also plotted closer to the center of the area at the right of alpine lakes species, while small prealpine and volcanic lakes were grouped on the right of the plot. This separation was still better emphasized in a 3D plot (Figures 1 and S1), where kryal, rhithral, krenal, potamal and lentic species were separated. Ablabesmyia spp.  Paracladius spp.   Insects 2022, 13, 911 10 of 22 Figure 1. plot of the species scores in the first 3 axes resulting from CA carried out from sites × species (nLp) matrix (the full set of species names is in Figure S1). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe).  Figure S1). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe).
Insects 2022, 13, 911 11 of 22 Figure 2. plot of the species scores in the first 2 axes resulting from CA carried out from sites × species (nLp) matrix and the fitted second-degree polynomial (the full set of species names is in Figure S2). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe).

Figure 2.
Plot of the species scores in the first 2 axes resulting from CA carried out from sites × species ( n L p ) matrix and the fitted second-degree polynomial (the full set of species names is in Figure  S2). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe).
A polynomial of the second degree was fitted to species scores of the two first axes (Figure 2), resulting in a multiple R-squared of 0.6845, and an adjusted R-squared of 0.6773 (F-statistic = 95.47 with 2 and 88 degrees of freedom, p-value = 2.2 * 10 −16 , residual standard error = 0.7344 with 88 degrees of freedom). The species more distant from the parabolic curve are visible in Figure 2 and are also evident in Figure S2, where all species names are plotted. Species from small Alpine lakes and from profundal zones of large lakes are the ones most deviating from the parabolic curve.
The environmental variables were included as passive variables in the map and were converted into factors with 6 different levels. When missing data were present, a level, plotted as void circles, grouped these data. The factors included were habitat, station (Figure 3), altitude, distance from the source (Figure 4), temperature, conductivity ( Figure 5), oxygen, total phosphorous ( Figure 6), nitrate, and ammonium nitrogen (Figure 7).  Table  S1) and the abbreviated name of the sampling station (see Table S1).   Table S1) and the abbreviated name of the sampling station (see Table S1).   Table  S1) and the abbreviated name of the sampling station (see Table S1).  The species × traits matrix p U s was also subjected to a correspondence analysis ( Figure 8, Table S5). The first 2 axes accounted for 70% and 21% of the total variance. Eigenvalues were 0.14 and 0.04, respectively. The first gradient separated species according to an upstream-downstream gradient, with the extreme scores assigned to altitude, distance from the source and conductivity. The second gradient separated species according to a trophic gradient, with the extreme scores assigned to oxygen, and N-NH 4       The sites × species matrix was transposed ( p L' n ) and a canonical constrained ordination (CCA) was carried out relating this matrix with the species × traits matrix p U s ( p L' n~p U s ) (Figures 9 and S3, Table S6). The first and second axis accounted for 7% and 5% of the total variance, with eigenvalues of 0.69 and 0.46, respectively. The scores of each species calculated according to the left (sites) and right (traits) set were joined by a line in the figure. The species showing preferences for the cold sites at high altitude were plotted at the bottom right of the graphs, the ones present in sites with high oxygen content at the bottom left, tolerant species such as Chironomus riparius (= thummi), Cricotopus (Cricotopus) trifascia and Virgatanytarsus present in high N-NO 3 , TP, N-NH 4, and low oxygen content waters were plotted at the top part of the graph, Rheopelopia, Uresipedilum, Tanypus from sites with high temperature and conductivity were mapped at the top left part. An arch/horseshoe effect was also visible here, with species preferring lentic waters plotted on the left, kryal and cold spring species on the bottom right, and species characterizing potamon at the top right.  indicate their values increase with water quality (blue) or with water pollution (red), according to EQR colors [35].
The sites × species matrix was transposed (pL'n) and a canonical constrained ordination (CCA) was carried out relating this matrix with the species × traits matrix pUs (pL'n~pUs) (Figures 9 and S3, Table S6). The first and second axis accounted for 7% and 5% of the total variance, with eigenvalues of 0.69 and 0.46, respectively. The scores of each species calculated according to the left (sites) and right (traits) set were joined by a line in the figure. The species showing preferences for the cold sites at high altitude were plotted  indicate their values increase with water quality (blue) or with water pollution (red), according to EQR colors [35].
The sites × species matrix was transposed (pL'n) and a canonical constrained ordination (CCA) was carried out relating this matrix with the species × traits matrix pUs (pL'n~pUs) (Figures 9 and S3, Table S6). The first and second axis accounted for 7% and 5% of the total variance, with eigenvalues of 0.69 and 0.46, respectively. The scores of each species calculated according to the left (sites) and right (traits) set were joined by a line in the figure. The species showing preferences for the cold sites at high altitude were plotted A comparison between factor loadings of species in canonical constrained and unconstrained ordination showed a good agreement in the species ordination, except for a few species such as Diamesa dampfi, Micropsectra notescens, Paratrissocladius, Paracricotopus, Psectrocladius sordidellus, Heterotrissocladius, which showed different scores in the CA first axis (calculated from n L p matrix) and in the CCA first axis (calculated from p L' n~p U s matrices) ( Figure S4, Table S6) and, as a consequence, were plotted at some distance from the regression line.
at the bottom right of the graphs, the ones present in sites with high oxygen content at the bottom left, tolerant species such as Chironomus riparius (= thummi), Cricotopus (Cricotopus) trifascia and Virgatanytarsus present in high N-NO3, TP, N-NH4, and low oxygen content waters were plotted at the top part of the graph, Rheopelopia, Uresipedilum, Tanypus from sites with high temperature and conductivity were mapped at the top left part. An arch/horseshoe effect was also visible here, with species preferring lentic waters plotted on the left, kryal and cold spring species on the bottom right, and species characterizing potamon at the top right. Figure 9. plot of the species scores (left) of the trait scores (right) in the first 2 axes resulting from CCA analysis carried out from pL'n~pUs matrices; the scores of the same species obtained with the first and second matrix are joined with a line (see Figure S4 for the full set of species names). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe). The different colors of environmental variables indicate their values increase with water quality (blue) or with water pollution (red), according to EQR colors [35].
A comparison between factor loadings of species in canonical constrained and unconstrained ordination showed a good agreement in the species ordination, except for a few species such as Diamesa dampfi, Micropsectra notescens, Paratrissocladius, Paracricotopus, Psectrocladius sordidellus, Heterotrissocladius, which showed different scores in the CA first axis (calculated from nLp matrix) and in the CCA first axis (calculated from pL'n ~pUs matrices) ( Figure S4, Table S6) and, as a consequence, were plotted at some distance from the regression line.
The sites × species matrix nLp was post-multiplied by species × traits pUs matrix to obtain a site × traits matrix (nL pUs). This nLpUs was also subjected to correspondence analysis ( Figure 10, Table S7). In this case, sites were rows and traits were columns. The first two axes accounted for 72% and 24% of the total variance, with eigenvalues of 0.05 and 0.02, respectively. The first axis reproduced an upstream-downstream and a water temperature gradient, and the second axis reproduced a water quality gradient ( Figure  10). This analysis does not allow the mapping of species, because the species (columns of the first matrix and rows of the second) do not appear in the product matrix. Figure 9. Plot of the species scores (left) of the trait scores (right) in the first 2 axes resulting from CCA analysis carried out from p L' n~p U s matrices; the scores of the same species obtained with the first and second matrix are joined with a line (see Figure S4 for the full set of species names). Different subfamilies are plotted with different colors: brown: Tanypodinae; blue: Diamesinae and Prodiamesinae; green: Orthocladiinae; orange: Tanytarsini (tribe); red: Chironomini (tribe). The different colors of environmental variables indicate their values increase with water quality (blue) or with water pollution (red), according to EQR colors [35].
The sites × species matrix n L p was post-multiplied by species × traits p U s matrix to obtain a site × traits matrix ( n L p U s ). This n L p U s was also subjected to correspondence analysis ( Figure 10, Table S7). In this case, sites were rows and traits were columns. The first two axes accounted for 72% and 24% of the total variance, with eigenvalues of 0.05 and 0.02, respectively. The first axis reproduced an upstream-downstream and a water temperature gradient, and the second axis reproduced a water quality gradient ( Figure 10). This analysis does not allow the mapping of species, because the species (columns of the first matrix and rows of the second) do not appear in the product matrix.
A discriminant analysis was carried out to test the goodness of classification of sites in different habitats when Chironomid taxa assemblages are used to discriminate among habitats (Tables 4 and S8). Both n L p and n L p U s matrices were subjected to multiple discriminant analysis, using habitat as a grouping factor. Percent of correct classifications was 46% for the n L p matrix and 47% for the n L p U s , emphasizing that the addition of the trait matrix does not improve the classification significantly. In any case, the result is that Chironomid assemblages are good discriminators between the different habitats.  Table S8.   Table S1). The different colors of environmental variables in the right figure indicate their increase with water quality (blue) or with water pollution (red), according to EQR colors [35].

ALA
A discriminant analysis was carried out to test the goodness of classification of sites in different habitats when Chironomid taxa assemblages are used to discriminate among habitats ( Table 4, Table S8). Both nLp and nLpUs matrices were subjected to multiple discriminant analysis, using habitat as a grouping factor. Percent of correct classifications was 46% for the nLp matrix and 47% for the nL pUs, emphasizing that the addition of the trait matrix does not improve the classification significantly. In any case, the result is that Chironomid assemblages are good discriminators between the different habitats.
A cluster analysis of species confirmed that the separation of species clusters is in agreement with different habitats (Figure 11). Table 4. Results of discriminant analysis: hits and misses in samples classification according to taxonomic and traits analysis. ALA: alpine lakes, B: brackish waters, K: kryal, LL: large lakes, LS: small lakes, ME: Mediterranean lakes, P: potamal, R: rhithral, S: krenal, V: volcanic lakes. Detailed results of Discriminant Analysis are provided in Table S8.   Table S1). The different colors of environmental variables in the right figure indicate their increase with water quality (blue) or with water pollution (red), according to EQR colors [35].

ALA
A cluster analysis of species confirmed that the separation of species clusters is in agreement with different habitats (Figure 11).

Discussion
Chironomid species distribution in the environment was confirmed to be related to ecological conditions. The distribution of Chironomids linked to biogeographic factors was never observed within the western Palearctic area, except for the species linked to glacial areas [36], so biogeographic factors are not considered in the present discussion.
Chironomids have been frequently used as indicators of past climatic change [37], while it is impossible to establish the occurrence of alien species [38], even if it is expected. Figure 11. Cluster analysis of species from sites × species ( n L p ) matrix.

Discussion
Chironomid species distribution in the environment was confirmed to be related to ecological conditions. The distribution of Chironomids linked to biogeographic factors was never observed within the western Palearctic area, except for the species linked to glacial areas [36], so biogeographic factors are not considered in the present discussion.
Chironomids have been frequently used as indicators of past climatic change [37], while it is impossible to establish the occurrence of alien species [38], even if it is expected. Some species such as Polypedilum nubifer are probably invaders [39], but it is impossible to state if and when they reached the West Palearctic region. It is well known that Chironomid distribution is related to ecological factors, such as water temperature [40,41], so an extension or reduction of the home range of a species is expected in relation to global warming [42]. Considering that different factors (e.g., winds) are good carriers of Chironomid species, and even long geographical barriers can be overcome [43], it is clear that species distribution is a quite dynamic process and at present alien species among Chironomids are expected, but their identity cannot be established.
However, with the ecological niche being known, it is possible to transpose the information given by each species into information about its habitat. From a mathematical point of view, the ecological niche can be expressed as a vector whose elements are the optimum values of the species for each factor, expressed as a weighted mean, while the measure of niche extension can be expressed as a weighted standard deviation [31]. The vectors can be aggregated to create a trait matrix p U s with p species as rows and s traits as columns. This p U s matrix was firstly proposed for calculating aquatic beetle traits and a fuzzy coding analysis was suggested to allow the inclusion of diverse kinds of biological information [44]. Species abundances can be expressed as a matrix n L p with n samples as rows and p species as columns. Matrix multiplication of the matrix n L p by the p U s matrix generates a n L p U s product matrix, with n sites as rows and s traits as columns; this approach was proposed for vegetation studies [45], and it was used for invertebrates living in running waters [46] and extended to Chironomids [22,23,47]: n M s = n L p * p U s This approach allows us to transpose the information given by a species list into ecological traits, allowing the construction of an index of environmental quality. Attempts to create the p U s matrix for Chironomids and other benthic invertebrates were a matter of many efforts [4,[48][49][50], but the results were obviously dependent on the database used for calculations. In the present paper, we tried to develop a new trait matrix considering the largest database available from collections of larval samples from both lotic and lentic habitats. This can produce results in disagreement with the ones produced considering a database constructed using data from lotic or lentic habitats alone. Indeed, traits of chironomids were often assigned without well-founded supporting information. For example, this was underlined in estimating the recovery of lakes after measures of restoration from acidification [51]. Significant differences were observed between traits developed for North American and European species [52] and between Scandinavian and Mediterranean species [53]. Lack of information may lead to apparently contradictory results. For example, hemoglobin content, tube building ability, feeding habit, voltinism, and body size of Chironomid larvae suggested that hemoglobin-rich species, with tube building capacity and short generation time, may be dominant in disturbed sites, while the reverse should be expected in less disturbed sites. However, this approach gave some unexpected results, such as the presence of: 1-hemoglobin-rich species in less disturbed sites; 2-species with long generation time in disturbed sites [47], and/or 3-small body-sized species in less disturbed habitats [24]. These apparently conflicting results were explained by positing that oxygen deficit was not the only factor determining disturbed conditions. It was posited that not all hemoglobin-rich species are tolerant to low oxygen levels, and also, within the same site, that different species at different depths may show variation in hemoglobin concentrations [14]. For example, species belonging to Polypedilum may be responsible for this conflicting result, because this hemoglobin-rich genus is often present in undisturbed sites, possibly due to the presence of small oxygen-poor microhabitats included in large oxygenrich habitats. Chironomini genera (Chironomus, Glyptotendipes, Polypedilum, Paratendipes, Microtendipes, etc.) are all hemoglobin-rich [14], but they show very different responses to pollution. The same is true for body size: the large Chironomus and Propsilocerus often prevail in disturbed sites, while it is expected that the small body-sized traits prevail in disturbed sites [54].
Another attractive approach is the so-called 4th corner solution problem [55,56], where the sites × species matrix n L p, the species × traits matrix p U s , and the sites × environmental variables matrix n R q are combined to produce a q D s = q R' n L p U s matrix, which allows a comparison between an expected and an observed community [56]. In the present case, the n R q matrix presents a lot of missing data, so this analysis was not performed. Caution is suggested in using this matrix approach to evaluate the ecological status because incomplete information available about the ecology of single taxa can lead to misleading results or false representations. This approach could be useful in the future when more accurate information will be available about different Chironomid species.
In the present study, as in many others [4,46,56], it is evident that Chironomid species respond to a limited number of factors, so they can be ordered according to a few gradients. We preferred to start the analysis ordering taxa with an unconstrained ordination method [31], because environmental data supporting the description of sampled sites were incomplete. Moreover, it is well known that the presence-absence of a species is not bound to the instantaneous point water condition, rather it is the result of an integration of factors over a relatively long time period, information that cannot be given by physical-chemical analysis.
Despite these limitations, the ordination of sites, based only on Chironomid species assemblages available in the present database, emphasized few major gradients responsible of the observed distributions: 1-a gradient separating lotic from lentic habitats, with species living in fast-running waters separated from species living in standing waters; 2a gradient emphasizing an upstream-downstream gradient in running waters, separating: (a) intolerant species living at high altitudes, low water temperatures, high oxygen concentrations, low conductivity, from (b) tolerant species living downstream, at higher temperatures, lower oxygen concentrations, higher conductivity, and salinity; 3-a trophic gradient separating species living in oligotrophic nutrient-poor waters from species living in organic-rich or eutrophic waters. Each of these gradients does not necessarily coincide with the principal axes resulting from canonical ordination. In the present case, the first axis separated lotic from lentic habitats, the second axis was explained as an oxygentemperature gradient, and the ordering of sites resulted in the classic arch or horseshoe effect [31,32]. This effect observed in the correspondence analysis [31] is generated by species data having unimodal distribution along a single gradient [32]; in the present case, it was a gradient from high altitude, cold, oxygen-rich, fast-flowing running waters observed in glacial streams, toward lowland, warmer, oxygen-poor, slow-flowing waters observed in lowland rivers, and continuing in still slow flowing, but cooling down and oxygen enriching waters, as observed in large lakes with increasing depth. Conductivity and nutrients were often included in this principal gradient, in several possible interactions. In relation to this principal gradient, each species can adjust with its own peculiarities, moving more or less far from this gradient. For example, species living in small-sized cold waters lakes at high altitudes (Zavrelimyia, Heterotrissocladius, Corynoneura, P. austriacus) and species living at high depth in large lakes (M. radialis, Paracladopelma) were displaced toward the center of the plot (Figure 2 and Figure S2).
Species cannot be clustered in well-defined groups, because only a few species are restricted to well-defined habitats, while most species are opportunistic. For example, few species belonging to Diamesa are restricted to kryal (Diamesa laticauda), but most (Diamesa tonsa, Diamesa zernyi) colonize different types of cold waters; some Orthocladiinae genera (Eukiefferiella, Rheocricotopus, Euorthocladius, Orthocladius, Cricotopus) characterize rhithral streams with moderate or fast current, but can be collected also in slow flowing waters; many Tanytarsini are typical of oligotrophic lakes, but are also common in springs and streams; many Chironomini genera (e.g., Dicrotendipes, Chironomus) characterize eutrophic lakes, but many of them live also in potamal river stretches and in littoral zones of lakes associated to vegetation (Endochironomus, Glyptotendipes) or to sand banks (Cryptochironomus, Harnischia) [57].
In conclusion, the key factors separating Chironomid species are confirmed to be substrate type, current velocity, water temperature, dissolved oxygen, conductivity, and nutrients, but these factors are differently related in various situations and anthropogenic stress and can contribute to creating other more complex interactions [15].
The advantage of having a matrix of ecological traits available ( p U s ) is the possibility of using only taxonomic assemblage structure information to evaluate the ecological status of an ecosystem, without the support of environmental data. This is a necessity when sampling campaigns include only the monitoring of macrobenthos; in this case, if a trait matrix is available, taxonomic information can be transposed into water quality assessment.

Conclusions
It is often stated that functional traits analysis is better than taxonomic composition analysis [23]. Indeed, this statement stresses the obvious, because the use of functional traits requires having a trait matrix available, and the development of a trait matrix implies the existence of sound taxonomic knowledge, needed to create the trait matrix. It is more appropriate to state that when a trait matrix is available, less thorough taxonomic knowledge is needed to evaluate the ecological status of a water body. In other words, a species groups list, instead of a more thorough species list, can be sufficient to analyze the system. The traits matrix approach has the advantage that a taxonomic species list can provide information comparable with the one given by a physical-chemical analysis when a trait matrix is available. If both a traits matrix p U s and an environmental variables matrix n R q are available, a further step can be done, calculating an expected ecological status and comparing it with an observed one [58].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13100911/s1. Figure S1 plot of the species scores in the first 3 axes resulting from CA carried out from sites × species matrix. As Figure 1, but with all species plotted. Figure S2 plot of the species scores in the first 2 axes resulting from CA carried out from sites × species matrix and the fitted second degree polynomial. As Figure 2, but with all species plotted. Figure S3 plot of the species scores in the first 2 axes resulting from CCA analysis carried out from p L' n and p U s matrices; the scores of the same species obtained with the first and second matrix are joined with a line. As Figure 9 (left), but with all species plotted. Figure S4 plot of the species scores in the first axis resulting from CA analysis of nLp (abscissa) against species scores in the first axis resulting from CCA analysis of p L' n~p U s (ordinate) with the fitted regression line, Table S1 input  data matrix Table S2 correlations between species and environmental variables, p values and number of samples.  Tables S5-S7. Table S5 correspondence analysis (CA) results of species × traits matrix p U s . Table S6 canonical constrained ordination (CCA) results between species × sites p L' n matrix and species × traits p U s matrix = p L' n~p U s . Table S7 correspondence analysis (CA) results of sites × traits n L p U s , matrix = sites × species multiplied by species × traits). Table S8 discriminant analysis results using habitat as grouping factor of the first two CA scores, calculated from from n L p and n L p U s matrices.