Diversity Patterns and Assemblage Structure of Non-Biting Midges (Diptera: Chironomidae) in Urban Waterbodies

: Urban waters are often neglected in biodiversity research; nonetheless, the number of aquatic microhabitats present in a city and the surrounding urban area is impressive. Twenty-two waterbodies in the Belgrade functional urban area (FUA) were investigated for faunistic and diversity patterns and to assess the effects of environmental factors on the differentiation of Chironomidae assemblages. A total of 66 chironomid taxa within four subfamilies was identiﬁed. Water quality at the studied sites, expressed by the water pollution index (WPI), varied signiﬁcantly. K-means clustering gave four homogenous groups of chironomid assemblages, which showed clear preferences to speciﬁc habitat conditions and tolerance to anthropogenic pressures. These groups had high values of alpha and beta diversity components. The main component of beta diversity was species turnover. Waterbody type, water temperature, pH, nutrients and overall pollution were the most important factors inﬂuencing the distribution and composition of chironomid assemblages, which revealed clear preferences of each assemblage type to the category of waterbody type and tolerances to environmental pressures.


Introduction
Chironomidae larvae (non-biting midges) are a species-diverse insect group abundant in many freshwater ecosystems, which represent the main food source for many predatory invertebrates, fish and birds [1][2][3]. This group displays a great capability of adapting to a wide range of environmental conditions, typically occurring at high densities with a key ecological function in lotic freshwater communities [4]. In lotic and lentic waterbodies of temperate regions, chironomids can be one of the most abundant and diverse taxonomic groups within benthos or epiphyton. This is especially true for very productive freshwater ecosystems, where they can represent more than 60% of the macroinvertebrate community [5,6]. The quantitative and qualitative composition of chironomid assemblages can be indicative of changes in water quality. Chironomidae play a crucial ecological role in the cycling of organic matter in rivers, in the export of energy to riparian habitats, and provide a valuable model system for understanding which environmental variables drive species richness [7]. The structure of chironomid assemblages reflect changes in aquatic ecosystems that are strongly correlated with changes in water quality and habitat degradation, clearly pointing to increased saprobity or hydromorphological degradation [8,9]. Chironomids are often used in water quality bioassessment with in situ monitoring programs [9,10] and in laboratory experiments that test their resilience and reactions to heavy metals, nano-and microplastics or organic pollutants [11,12]. Together with Oligochaeta, non-biting midges are usually dominant and sometimes the only present invertebrate taxa in heavily modified waterbodies regardless of whether they are of artificial or natural origin [13][14][15]. Urbanisation and landscape conversion in sub-urban areas are a threat and the main cause of aquatic system degradation [16,17]. Some of the known stressors influencing ecological traits and system functioning of urban waterbodies can be from an extrinsic (catchment) source or created within the flow, such as untreated stormwater runoff, various point sources of pollution, septic system leakage, dams or rip-raps [18,19].
In the Belgrade FUA, there are numerous heavily modified and artificial waterbodies (reservoirs, canals and rivers) exposed to different types and intensities of anthropogenic pressures. Having various uses and origins, they differ in hydromorphological characteristics and water quality because of their position and distance from the urban area. Belgrade, with its 1.7 million inhabitants, has no system for treating municipal wastewaters [20]. Some of the running waters flow through industrial zones, agricultural and urban areas and have a role in precipitation drainage, as well as wastewater removal from several rural areas, thus being subjected to combined pressures [21,22]. Reservoirs in the Belgrade suburban area are also under high anthropogenic influence, such as wastewater discharge from surrounding settlements, fishing and other recreational activities. Canals and reservoirs are additionally under continual pressure caused by flow regulation. Consequently, they lack a natural water regimen as water is actively drained and pumped. Further, dredging of bed material and removal of aquatic and riparian vegetation contributes to the effect of a complex multi-stressor environment. Many of these reservoirs and canals lie in sanitary protection zones of the water supply system of Belgrade.
Keeping in mind the importance of such aquatic systems, and assessing the intensity of anthropogenic influence and possible degradation, we surveyed modified waterbodies in the Belgrade FUA to investigate whether assemblages of Chironomidae larvae can provide clear signs of environmental changes. We aimed to detect and describe (i) the faunistic patterns of chironomid assemblages and (ii) the patterns of diversity components in analysed chironomid assemblages and to examine (iii) the effects of environmental factors on the differentiation of chironomid assemblages.

Data Analysis
The qualitative composition of chironomids was determined for each site, along with species occurrence frequencies (F = 0-1). The ASTERICS software package Version 3.1.1. [24] was used for the assessment of data and the calculation of metrics.
The modified water pollution index (WPI) [34] was used to estimate water quality classes ( Table 1). The WPI is calculated as the sum of the ratio of the measured average value and the standard threshold values for each parameter, divided by the number of used parameters. The standard threshold values for all parameters were specific for each country, given as the national legislative [23], which should minimise bias caused by ecological and geographical differences. The ASTERICS software package, version 3.1.1 [24], was used for the calculation of metrics used for the assessment of the ecological potential for these artificial waterbodies. The ecological analysis of the macroinvertebrate community structure was conducted for each site to calculate WPI. The number of taxa, ASPT (average score per taxon), BMWP (biological monitoring working party score) [35], the Saprobic index (S) [36] using bioindicator valences of each taxon according to [37], α-diversity index (H') [38] and the percentage of the subfamily Tubificinae (Oligochaeta) in macroinvertebrate communities were calculated and used in the WPI calculation. To reveal the variability patterns of analysed chironomid assemblages, we used two powerful non-hierarchical classification methods: K-means clustering [39] and Bayesian classification [40]. Hill et al. [41] emphasised that the number of misclassifications is a key parameter in assessing the analytical power of clustering methods. Contrary to numerous agglomerative and divisive methods, K-means clustering and Bayesian classification enable the allocation of misclassified assemblages to their most similar cluster. Despite their ability to minimise the number of misclassifications, both K-means clustering and Bayesian classification have two serious drawbacks. Both methods allocate assemblages into a pre-specified number of clusters. The main drawback of these methods is subjectivity in the initial selection of the number of clusters. Marinković et al. [42] proposed a simple procedure to avoid this problem. The procedure selects the number of clusters by maximising the variance ratio: where σ 2 B is the between-group variance (i.e., variance of the cluster centroids), and σ 2 W is the within-group variance (the sum of the variances within each k cluster). Maximising the variance ratio ensures that overlap of homogeneous clusters is minimised. Another drawback of both is that K-means clustering and Bayesian classification are associated with the local minima problem. Two closest points in a Euclidean space correspond to the local minimum in one-dimensional space (a linear local minimum). Three closest points in the space represent the local minimum in a two-dimensional space (a planar local minimum). Both planar and linear local minima can disintegrate compact global clusters and replace them with numerous small clusters. To avoid this undesirable effect, we specified a threshold of four points as the minimum size of initial clusters. Selection of initial clusters with at least four points specifies the spatial configuration of the clusters and eliminates the undesirable effects of linear and planar local minima.
MANOVA [43] was used to find a combination of species that maximally discriminates extracted clusters of chironomid assemblages.
Cluster analyses based on Ward's method [44] was used for grouping Chironomidae assemblage types.
For each type of chironomid assemblage, we analysed the components of alpha diversity (species richness, Shannon index and Shannon equitability). The components of beta diversity were detected using the procedures described by Baselga [45] and Podani et al. [46].
The stepwise forward selection (FS) procedure [40] was used to detect environmental variables with statistically significant effects on the chironomid assemblages. At each step of the procedure, we expanded the multiple regression model by adding an environmental variable that explains most of the residual variance (i.e., the variance of faunistic data, not explained by previously selected environmental variables). The statistical significance of the hypothesis that species assemblages are independent of the selected environmental variable was assessed using the non-parametric Monte Carlo permutation test (3000 permutations, p < 0.05).
To reduce the weighting of dominant groups, species-abundance data were transformed using the formula: where x is the number of recorded individuals. The effects of environmental variables on the faunistic differentiation of assemblages were assessed using canonical correspondence analysis (CCA) [47].
Pearson's correlation test was used to determine the correlation between WBT and WPI and their correlation with measured environmental parameters. For testing data normality, we used the Kolmogorov-Smirnov test.
Statistical analyses were performed using FLORA software [48], updated version.

Results
In total, 66 chironomid taxa in four subfamilies (Prodiamesinae, Orthocladiinae, Tanypodinae and Chironominae) were found during the study period. The most diverse and abundant was the Chironominae subfamily (21 genera, 37 taxa). Polypedilum was the most diverse genus in the Chironominae subfamily with nine species. Cricotopus was the most diverse genus from the Orthocladiinae subfamily. The highest abundance of chironomid species was detected in wadeable rivers, while in reservoirs, canals and non-wadeable rivers, the abundance of chironomid species was comparable. Cricotopus gr. sylvestris and Cricotopus bicinctus (Meigen, 1818) were found in various waterbody types, but they were most abundant in wadeable rivers. In reservoirs, the most abundant was Ablabesmyia monilis agg., while in canals, the most numerous were C. gr. sylvestris and Parachironomus gracillior (Kieffer, 1918).

Chironomid Assemblage Classification
Classification of chironomid assemblages was performed based on the faunistic similarity of the analysed assemblages. The main drawback of hierarchical classification methods is the inability to correct misclassifications. However, non-hierarchical clustering methods allow for the allocation of misclassified species assemblages to their most similar cluster. We, therefore, performed classification of the analysed assemblages using K-means clustering and Bayesian classification. These methods are the most powerful variants of non-hierarchical clustering methods [40]. The calculated variance ratio indicated that K-means clustering produced more acceptable results than Bayesian classification. The main drawback of Bayesian classification is the rigid assumption that all variables must be normally distributed, as concluded in Sekulić et al. [49]. The results of K-means clustering  Table 2. To eliminate the undesirable effects of local minima, we used four species assemblages as the minimum size of initial clusters. Since we investigated 22 chironomid assemblages, the greatest number of initial clusters with at least four assemblages was five. The greatest ratio of between-group variance to within-group variance was detected for four clusters. Therefore, we performed K-means clustering using four pre-defined clusters and named these clusters (types) with capital letters A, B, C and D. Table 2. Dependence of ordering results on the subjectively selected number of clusters. The greatest ratio of between-group variance (BGV) to within-group variance (WGV) σ 2 B /σ 2 W assures that the overlap of homogeneous clusters is minimised. In our dataset, the greatest variance ratio (bolded numbers) was obtained for four clusters.

K-Means Clustering Bayesian Classification
No. of clusters WGV σ 2 MANOVA provided a combination of species that maximally discriminated between four clusters of assemblages (high variance ratio R 2 0.452884).
Species from the Conchapelopia aggregate as well as Procladius species discriminated assemblage type A from other assemblage types. Type B stands out by the presence of Microchironomus tener (Kieffer, 1918). Polypedilum albicorne (Meigen, 1838) was the discriminating species for community type C, while Endochironomus albipennis (Meigen, 1830) and Monopelopia tenuicalcar (Kieffer, 1918) were a distinguishing feature of assemblage type D (Figure 2).
Cluster analyses using Ward's method organised the assemblages into two larger groups based on the similarity of species composition, with each encompassing two assemblage types. The first group incorporates assemblage types A and B, and the second types C and D (Figure 3).
For each group of assemblages, a set of diagnostic species was established based on their frequency of occurrence (Table S3).
Correspondence between faunistic groups (X-axis) and waterbody types (colours) is presented in a histogram (Figure 4).
Assemblage type A inhabited only wadeable rivers. Assemblage type C occurred in slow-flowing waters such as reservoirs and canals, while assemblages of type D were found mainly in canals and in some wadeable rivers. Assemblage type B was the only type detected in non-wadeable rivers. This type was also found in wadeable rivers and only in one reservoir (Figure 4). Cluster analyses using Ward's method organised the assemblages into two larger groups based on the similarity of species composition, with each encompassing two assemblage types. The first group incorporates assemblage types A and B, and the second types C and D (Figure 3).
For each group of assemblages, a set of diagnostic species was established based on their frequency of occurrence (Table S3) (Kieffer, 1918) was the diagnostic/main species for assemblage type D and was associated with Glyptotendipes pallens agg. and Xenopelopia species. Correspondence between faunistic groups (X-axis) and waterbody types (colours) is presented in a histogram (Figure 4).  . Hierarchical cluster analyses (Ward's method) of the Chironomidae assemblages. Fou assemblage types (obtained through K-means clustering) are in capital letters A, B, C, D coloured as indicated. The numbers represent study sites codes (see Figure 1).
Correspondence between faunistic groups (X-axis) and waterbody types (colours) i presented in a histogram (Figure 4). Assemblage type A inhabited only wadeable rivers. Assemblage type C occurred in slow-flowing waters such as reservoirs and canals, while assemblages of type D wer found mainly in canals and in some wadeable rivers. Assemblage type B was the only type detected in non-wadeable rivers. This type was also found in wadeable rivers and only in one reservoir (Figure 4).

Diversity Components
The highest alpha diversity (expressed as Shannon index) was recorded in assemblage types A and B, while in C and D, the values of the Shannon index were lower The same trend was observed for species richness, while the Shannon equitability index was almost identical in all assemblage types ( Figure 5).

Diversity Components
The highest alpha diversity (expressed as Shannon index) was recorded in assemblage types A and B, while in C and D, the values of the Shannon index were lower. The same trend was observed for species richness, while the Shannon equitability index was almost identical in all assemblage types ( Figure 5). High, almost identical, beta diversity values were detected in all analysed types of assemblages using Baselga's method ( Figure 6). The dominant component of beta diversity was species turnover, which was high, while nestedness was low. Analyses of beta diversity using the Podani method gave different results between different types of assemblages ( Figure 6). While overall beta diversity was similar, its components varied from type to type. Nestedness was lowest in type A and highest in type C, while in types B and D, it was similar. Species turnover was the same in types C and D, which was lower than in types A and B, with type A showing the highest value of this component of beta diversity. High, almost identical, beta diversity values were detected in all analysed types of assemblages using Baselga's method ( Figure 6). The dominant component of beta diversity was species turnover, which was high, while nestedness was low. Analyses of beta diversity using the Podani method gave different results between different types of assemblages ( Figure 6). While overall beta diversity was similar, its components varied from type to type. Nestedness was lowest in type A and highest in type C, while in types B and D, it was similar. Species turnover was the same in types C and D, which was lower than in types A and B, with type A showing the highest value of this component of beta diversity. diversity was species turnover, which was high, while nestedness was low. Analyses of beta diversity using the Podani method gave different results between different types of assemblages ( Figure 6). While overall beta diversity was similar, its components varied from type to type. Nestedness was lowest in type A and highest in type C, while in types B and D, it was similar. Species turnover was the same in types C and D, which was lower than in types A and B, with type A showing the highest value of this component of beta diversity.

Patterns in the Chironomidae-Environment Relationship
Partial forward selection analysis (Table 3) identified ten environmental variables as significant for faunistic differentiation of the analysed urban waterbodies.

Patterns in the Chironomidae-Environment Relationship
Partial forward selection analysis (Table 3) identified ten environmental variables as significant for faunistic differentiation of the analysed urban waterbodies.  CCA showed that the Chironomidae assemblage types were clearly differentiated with respect to nutrients (nitrate, nitrite and ammonium concentration), water temperature, and pH gradients. The importance of other environmental variables was lower. CCA indicated that oxygen demand (chemical and biological oxygen demand) produced effects on faunistic differentiation of the analysed Chironomidae assemblages. WBT also played a role in assemblage differentiation (Figure 7). Based on the values of WPI and water quality, the four types of assemblages were different (Figure 8). Assemblage type A was present in sites with various water qualities. Most of these sites were characterised by higher values of WPI and deterioration of water quality (classes IV, V and VI) (Table S2). Assemblage type C was present in less polluted waters (classes II and III of water quality, with lowest WPI values). Assemblage types B and D mainly occurred on sites with pure or moderately polluted water. Extremely high values of WPI (heavily impure water) were recorded at some sites inhabited by assemblage type B (Figure 8). Based on the values of WPI and water quality, the four types of assemblages were different (Figure 8). Assemblage type A was present in sites with various water qualities. Most of these sites were characterised by higher values of WPI and deterioration of water quality (classes IV, V and VI) (Table S2). Assemblage type C was present in less polluted waters (classes II and III of water quality, with lowest WPI values). Assemblage types B and D mainly occurred on sites with pure or moderately polluted water. Extremely high values of WPI (heavily impure water) were recorded at some sites inhabited by assemblage type B (Figure 8).
Since WBT and WPI are considered as derived variables, their correlations with environmental variables were analysed and are presented in Table 4. WPI and WBT were correlated. Nutrients (total phosphate, total nitrogen, nitrate, nitrite and ammonium concentrations) and water temperature showed significant correlation with these two derived variables. Distribution of four chironomid assemblage types in different classes of water quality calculated through WPI; water quality class based on WPI is in Roman numerals II-VI, from very pure to heavily impure (see Table 1); four assemblage types are in capital letters A, B, C, D coloured as indicated.
Since WBT and WPI are considered as derived variables, their correlations with environmental variables were analysed and are presented in Table 4. WPI and WBT were correlated. Nutrients (total phosphate, total nitrogen, nitrate, nitrite and ammonium concentrations) and water temperature showed significant correlation with these two derived variables.   (b) Distribution of four chironomid assemblage types in different classes of water quality calculated through WPI; water quality class based on WPI is in Roman numerals II-VI, from very pure to heavily impure (see Table 1); four assemblage types are in capital letters A, B, C, D coloured as indicated.

Discussion
Chironomids, as one of the most diverse aquatic macroinvertebrate groups, have a wide spectrum of biological and ecological preferences [10]; nevertheless, studies of chironomid assemblages mainly focus on the presence/absence of species rather than the abundance and assemblage structure [50]. The diverse fauna of Chironomidae larvae detected during this study enabled us not only to assess biodiversity within them, but also to monitor negative anthropogenic influences and degradation of studied waterbodies. The distribution of the chironomids mainly followed the a priori classification given by the non-chironomid macroinvertebrates [51]. We tried to determine whether assemblages of Chironomidae larvae in urban waters can provide clear indications of environmental changes.
Leszczyńska et al. [50] state that the knowledge of assemblage composition is necessary to assess some diversity components and other (inter-)assemblage parameters, as well as to assess assemblage dependence on environmental factors. High diversity of species characterised all four assemblage types found in Belgrade FUA while exhibiting differences in diversity components. Assemblage types are distinguished by the influence of different components on overall beta diversity. The ratio of components of beta diversity in types A and B pointed to the availability of different microhabitats in each site. A greater availability of habitats and microhabitats, such as sediment composition and the presence of macrophytes and detritus-which provides food, shelter for burrowing, mining and protection from predators-are associated with habitat heterogeneity and can harness the high diversity of chironomids [52][53][54]. Nevertheless, a huge diversity of the habitats exposed to multi stressors leads to the unclear relationship between biotic and abiotic components [55]. In the other two types (C and D), the nestedness component was slightly higher, pointing to the presence of an ecological gradient decreasing the number of species from site to site and their composition. Aquatic insects display a decrease in alpha diversity as a response to the urbanisation process, where, in sites exposed to substantial influence, only highly tolerant species prevail [16]. However, Chironomidae assemblages in Belgrade urban waters showed relatively high alpha diversity, but components of beta diversity revealed the aftereffects of urbanisation and pollution, high nestedness in canals, reservoirs, and heavily polluted rivers (type C and D) indicated that only tolerable species remained under high anthropogenic pressure.
Assemblage types A and B were similar in species composition, also preferring similar habitats such as wadeable and non-wadeable rivers. Preference for slow and stagnant waters found in canals and reservoirs was exhibited by assemblage types C and D, which were also grouped together based on species composition by cluster analyses.
Diagnostic taxa of assemblage's type A and B Rheocricotopus, Thenemanniella, Microchironomus tener and their associated species are known inhabitants of running waters and reservoirs [28,29,56], which are the types of ecosystems inhabited by the aforementioned assemblage.
Slow-flowing and stagnant waters with ample vegetation are preferred habitats of Glyptotendipes paripes, Monopelopia tenuicalcar, Parachironomus gracillior, Dicrotendipes pulsus, Kiefferulus tendipediformis, Glyptotendipes pallens agg. and Xenopelopia [29,32,33,56]. Ecological conditions in sites in which C and D assemblages types were detected corresponded with these ecological preferences. This is also in agreement with the classification of these sites as reservoirs or canals, with the addition of one wadeable river characterised by slow flow and large amounts of vegetation.
Although waterbody type is one of the most important factors determining the distribution of chironomids, the presence of pollution and other anthropogenic pressures (such as habitat degradation) can influence, to a great extent, the abundance and structure of chironomid assemblages [51]. Streams receiving waste effluents are characterised by lower chironomid species richness and the development of more dense populations of Chironomus riparius [57,58]. Leszczyńska et al. [50] found that C. riparius is the most abundant in low order streams, with low velocity and dense riparian vegetation, preferring stagnant water and soft sediments. It is also known that C. riparius may inhabit organically enriched and heavily polluted waterbodies, having efficient oxygen regulation [59]. The results presented herein show the same patterns, e.g., a high abundance of Chironomus species in heavily polluted rivers (the Topčiderska, Barička and Barajevska rivers). The Topčiderska river, which flows through the industrial zone and is the recipient of communal wastewaters, was the only site with Eukiefferiella claripennis, Parametriocnemus stylatus and Tvetenia clavescens, species found to be tolerant to habitat degradation [58,60]. Non-wadeable rivers were characterised by the dominance of the subfamily Chironominae (Dicrotendipes nervosus, Polypedilum nubeculosum and Chironomus species). Milošević et al. [9] reported a similar assemblage structure, namely one that was monotonous and driven by frequent and dominant taxa from the Chironominae subfamily, which were also reported in other studies as common dominant taxa in non-wadeable rivers [55,61].
Cricotopus bicinctus, which has already been documented as an indicator species for organic pollution [10], was abundant in wadeable rivers (the Topčiderska, Barajevska and Beljanica rivers) but also present in non-wadeable rivers.
The high values of diversity components observed in assemblage types A and B, despite the high pollution (expressed through WPI and nutrients), were likely supported by the availability of microhabitats and other favourable ecological conditions of the ecosystems they resided in (wadeable and non-wadeable rivers). Canals and reservoirs were inhabited by assemblages (C and D) that exhibited slightly lower alpha and beta diversity components. The limiting factors in these cases could be the uniformity of habitats and the limited availability of resources (food, shelter, substrate), but also the hydro-technical regime and maintenance work on these heavily modified and artificial waterbodies.

Conclusions
Urban waters in the Belgrade FUA harbour very diverse chironomid fauna. Based on their preferences for specific waterbody types and tolerance to environmental pressures, Chironomid assemblages can be grouped into several different types characterised by unique species composition.
Our study showed that chironomids could serve as useful indicators of anthropogenic pressures in various waterbody types due to the different sensitivities of the species towards the alteration of environmental conditions in their habitats.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14030187/s1, Table S1: Sampling sites and their characteristics; Table S2: Detected pressures at sampling sites and water quality classes according to the Water Pollution Index (WPI) value; Table S3: Taxa list and frequencies of species in four clusters of species assemblages.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.